← Back to Mesh-based solvers Module · scimba_jax

DG solver

In the Discontinuous Galerkin (DG) method, each cell has its own degrees of freedom, which are not shared with neighbouring cells. Continuity between cells is enforced weakly, by a numerical flux computed at every face.

1. 2D Laplacian problem, with one block per cell

\(-\Delta u = f\)

Compared to the FEM solver, the mesh and the basis family are the same, but the Variables class is different: VariablesDG stores a full set of degrees of freedom per cell, and nodes are not shared. The cells are instead coupled by a flux, added on every interior and boundary face. Here, this flux is SIPGFlux (symmetric interior penalty). 

import jax.numpy as jnp
from scimba_jax.linear_approximation.basis.analytic_bases import local_lagrange_basis
from scimba_jax.linear_approximation.basis.general_bases import AnalyticBasis
from scimba_jax.linear_approximation.galerkin.dg.elliptic_dg_scheme import (
    EllipticDGscheme,
)
from scimba_jax.linear_approximation.galerkin.dg.flux import SIPGFlux
from scimba_jax.linear_approximation.meshes.mesh import Mesh
from scimba_jax.linear_approximation.quad.gauss_quad import UnitSquareTensorized
from scimba_jax.linear_approximation.variables.variables_dg import VariablesDG
from scimba_jax.mapping.mapping import InvertibleFunction, Mapping
from scimba_jax.physical_models.abstract_physical_weak_model import (
    AbstractPhysicalWeakModel,
)
from scimba_jax.physical_models.classical_weakform.laplacian_weak_form import (
    LaplacianWeakForm,
)
order = 2
mesh = Mesh(
    dim=2,
    n_cells=(24, 24),
    ref_quad=UnitSquareTensorized(dim=2, order=2 * order + 1),
    mapping=Mapping(mappings=[InvertibleFunction(lambda x: x, lambda y: y)]),
)
basis = AnalyticBasis(
    nb_basis=(order + 1) ** 2,
    out_dim=1,
    mesh=mesh,
    local_basis=lambda y, i, m: local_lagrange_basis(y, i, m, order=order, out_dim=1),
    basis_type="scalar",
)
# One block of DOFs per cell, not shared between cells: the jumps across
# faces are controlled by the flux term.
variables = VariablesDG(basis=basis, nb_variables=1)
def source(x):  # same manufactured solution as the FEM page
    return 2 * jnp.pi**2 * jnp.sin(jnp.pi * x[0]) * jnp.sin(jnp.pi * x[1])
pde = LaplacianWeakForm(dim=2, f=source)
model = AbstractPhysicalWeakModel.from_weak_form(pde, dirichlet=lambda x: jnp.zeros(1))
# Continuity is enforced by a penalty, whose strength scales like 1/h.
flux = SIPGFlux(sigma=10.0 * order * (order + 1))
scheme = EllipticDGscheme(model, variables, flux, use_scan_quad=False)
scheme = EllipticDGscheme.solve(scheme, matrix_free=True, tol=1e-10)
# Same problem, order and mesh as the FEM page: 1.7e-05 relative L2 error.

2. Numerical fluxes

\(\widehat{F} = \{\nabla u\}\cdot n - \dfrac{\sigma}{h}\,[\![u]\!]\)

All DG schemes use fluxes in the same way. Given the values on both sides of a face (trace, gradient and outward normal), the flux returns the quantity exchanged between the two sides. SIPGFlux is the default flux for diffusion terms. An upwind flux follows the same interface for advection terms, and so does a Riemann solver for hyperbolic systems. 

class AbstractFlux(ScimbaPytree):
    """Interface of a DG numerical flux."""
    def __call__(self, varL, varR):
        """Contributions of both sides on an interior face."""
        raise NotImplementedError
    def boundary_call(
        self, var_interior, n_interior, dirichlet_val
    ):
        """Contribution of the interior side on a boundary face."""
        raise NotImplementedError
# SIPGFlux is consistent, symmetric (the operator stays SPD)
# and penalised by
# a local sigma/h, which is required on non-uniform meshes.
# An upwind flux uses the same interface, for advection terms.

3. Time-dependent problems

\(\partial_t u - \Delta u = 0\)

TimeDiscreteDGscheme solves an elliptic DG problem at every stage of a Butcher tableau, like TimeDiscreteFEscheme on the FEM page. When choosing a tableau, note that the interior penalty flux makes the problem stiffer than the operator alone. Therefore, on the same mesh, an explicit tableau requires a noticeably smaller time step than for FEM. An implicit tableau does not have this restriction. 

from scimba_jax.linear_approximation.basis.analytic_bases import local_taylor_basis
from scimba_jax.linear_approximation.galerkin.dg.time_discrete_dg_scheme import (
    TimeDiscreteDGscheme,
)
from scimba_jax.time_discrete.butcher_tableau import build_implicit_euler_tableau
# Same heat equation as on the FEM page. The penalty makes the explicit Euler
# time step ~5x smaller than for FEM, so implicit Euler is used here.
n_cells, order = 16, 1
mesh_1d = Mesh(
    dim=1,
    n_cells=(n_cells,),
    ref_quad=UnitSquareTensorized(dim=1, order=4),
    mapping=Mapping(mappings=[InvertibleFunction(lambda x: x, lambda y: y)]),
)
basis_1d = AnalyticBasis(
    nb_basis=order + 1,
    out_dim=1,
    mesh=mesh_1d,
    local_basis=lambda y, i, m: local_taylor_basis(y, i, m, order=order, out_dim=1),
    basis_type="scalar",
)
scheme = TimeDiscreteDGscheme(
    spatial_weak_form_factory=LaplacianWeakForm(dim=1, f=lambda x: jnp.zeros(())),
    variables=VariablesDG(basis=basis_1d, nb_variables=1),
    flux=SIPGFlux(sigma=(order + 1) * (order + 2), h=1.0 / n_cells),
    butcher_tableau=build_implicit_euler_tableau(),
    dt=0.1 / 8,
    dirichlet=lambda x: jnp.zeros(1),
)
dofsl_init = scheme.initialize(lambda x: jnp.sin(jnp.pi * x))
dofsl_final, history = scheme.solve(dofsl_init, t0=0.0, nt=8)
# Same 8 steps, order-1 Taylor basis: 5.5e-02 relative L2 error, as for FEM
# with the same order and mesh.

The block structure of DG also affects the choice of an efficient solver. The Multigrid page presents a block-Jacobi smoother, which uses the dense block of each cell instead of a single diagonal entry, and the Nonlinear solvers page covers nonlinear fluxes.