Discontinuous Galerkin (DG) Schemes

Updated 11 June 2026

Discontinuous Galerkin schemes are finite element methods that use discontinuous piecewise polynomials to approximate PDEs, ensuring high-order accuracy and local mass conservation.
They utilize numerical fluxes and penalty terms to balance stability and efficiency, effectively handling hyperbolic, elliptic, and mixed-type problems.
DG methods offer robust parallel scalability and adaptability in complex geometries, making them applicable in fluid dynamics, astrophysics, and engineering simulations.

Discontinuous Galerkin (DG) schemes are a broad class of finite element methods for the discretization of partial differential equations (PDEs). Characterized by piecewise polynomial approximations that are allowed to be discontinuous across element boundaries, DG methods combine advantageous features of high-order finite volume and finite element frameworks. These methods are widely used for hyperbolic and parabolic PDEs, elliptic problems, and mixed systems, with particular relevance for variable coefficients, complex geometry, and adaptivity requirements. DG schemes unify mathematical rigor, flexibility, and robust parallel scalability, making them foundational for modern computational physics, engineering, and applied mathematics.

1. Mathematical Formulation and Core Principles

DG schemes operate by partitioning the computational domain into non-overlapping elements, over each of which the solution is approximated by polynomials up to degree $p$ . No continuity is enforced at inter-element boundaries, so trial and test functions are discontinuous across faces. For a general conservation law or first-order system,

$\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{u}) = \mathbf{S}(\mathbf{u}),$

the semi-discrete DG formulation on each element $K$ reads: $\int_K \partial_t \mathbf{u}_h \, v_h\,dx - \int_K \mathbf{F}(\mathbf{u}_h) \cdot \nabla v_h\,dx + \int_{\partial K} \widehat{\mathbf{F}} \cdot n\, v_h\,ds = \int_K \mathbf{S}(\mathbf{u}_h) v_h\,dx \quad \forall v_h \in \mathcal{V}_h^p,$ where $\widehat{\mathbf{F}}$ is a numerical flux (e.g., upwind, Lax–Friedrichs, HLL), $n$ is the local unit normal, and $\mathcal{V}_h^p$ denotes the space of discontinuous, degree- $p$ polynomials per element (Persson, 2012).

Higher-order elliptic problems are handled by rewriting the PDE in first-order form and incorporating additional variables (e.g., local gradients $q=\nabla u$ ), then coupling the equations using local DG variants such as the LDG method. Penalty terms or numerical fluxes enforce weak continuity and control inter-element jumps.

The key structural properties of DG methods include:

Local mass conservation due to the element-wise weak formulation.
High-order accuracy achievable via arbitrarily high polynomial degrees.
Matrix sparsity—each DOF couples only to a compact element-neighborhood.
Ease of $hp$ -adaptivity and geometric flexibility, due to the element-local nature of the basis functions.

2. Fluxes, Penalty Terms, and Inter-element Coupling

The design of numerical fluxes is central to the stability and consistency of DG schemes:

Flux Type	Description	Typical Usage
Upwind / Riemann	Respects physics, stability, monotonicity	Hyperbolic problems
Central (no dissipation)	Energy-conserving, can be oscillatory	Special (e.g., Hamiltonian) flows
Interior penalty	Enforces weak continuity for elliptic/parabolic	Elliptic, viscous problems

For second-order or higher PDEs, numerical fluxes are constructed to handle both the solution and its gradients or auxiliary variables on interfaces. Modern SIPDG (symmetric interior penalty DG) formulations use geometry-robust penalties such as $\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{u}) = \mathbf{S}(\mathbf{u}),$ 0, with $\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{u}) = \mathbf{S}(\mathbf{u}),$ 1, to guarantee coercivity without requiring mesh shape-regularity (Kashiwabara et al., 2020).

LDG-type (local DG) schemes for diffusion employ alternating “switch” functions and up/downwinding across global lines, yielding minimal-dissipation or over-stabilized variants (Persson, 2012). Fluxes for nonlinear, nonconservative systems (e.g., GRMHD) follow DLM path-conservative theory, introducing jump integrals over admissible paths in state space (Fambri et al., 2018).

3. High-order Accuracy and Error Analysis

DG schemes achieve design accuracy via:

High-degree local polynomials: $\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{u}) = \mathbf{S}(\mathbf{u}),$ 2th order for degree $\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{u}) = \mathbf{S}(\mathbf{u}),$ 3.
Exact or consistent quadrature: Volume and surface integrals are evaluated using suitably high-order Gauss rules or quadrature-free approaches (e.g., VEM-DG via projection) (Boscheri et al., 2023).
Proper penalization: Adequate penalty parameters ensure optimal convergence even on anisotropic meshes (Kashiwabara et al., 2020).

For smooth solutions, theoretical and computational studies confirm optimal error rates: $\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{u}) = \mathbf{S}(\mathbf{u}),$ 4 For over-stabilized LDG (nonzero $\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{u}) = \mathbf{S}(\mathbf{u}),$ 5), the gradient variable exhibits superconvergent behavior, enabling hybrid postprocessing (Persson, 2012).

Adaptive strategies such as $\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{u}) = \mathbf{S}(\mathbf{u}),$ 6-AMR, VEM projections, multiscale spectral GFEM, and a posteriori limiting (subcell FV, MOOD) further enhance convergence in the presence of localized features or singularities (Alber et al., 24 Oct 2025, Boscheri et al., 2016).

4. Numerical Methods, Time Integration, and Hybridizations

A broad array of time-integration and hybridization techniques interface with spatial DG discretizations:

Explicit Runge-Kutta (RKDG): Standard for explicit hyperbolic problems, though CFL-limited.
Implicit/IMEX RK and Multiderivative Integrators: Allow larger timesteps for diffusive, stiff, or multiscale problems; maintain high order while reducing stage count (Jaust et al., 2015, Liu et al., 2021).
ADER-DG: Single-step, arbitrarily high order methods unified in space and time, supporting communication-avoiding, vectorized, single-stage updates ideal for exascale/high-performance computing (Dumbser et al., 2018, Fambri et al., 2018).
Semi-Lagrangian DG: Unconditionally stable methods for transport and diffusive equations via weak-form characteristics, explicit or semi-implicit updates, and high-order weak Taylor expansions (Bokanowski et al., 2015, Tavelli et al., 2024).
Stage-dependent RKDG (sdRKDG): Employs different polynomial degrees at various RK stages, enlarging stability regions and reducing computational cost without loss of accuracy (Chen et al., 2024).

Hybridized DG (HDG) and VEM-DG approaches further reduce the global coupling by introducing face or virtual unknowns, enhancing both efficiency and robustness on general meshes (Boscheri et al., 2023, Jaust et al., 2015).

5. Adaptivity, Stability, and Limiting Strategies

Advanced DG schemes incorporate adaptive mechanisms to address robustness, efficiency, and accuracy in underresolved flow or shock-dominated regimes:

Volume term adaptivity (v-adaptivity): Selectively switches, per element and RK stage, between standard weak-form (fast) and entropy-conservative/dissipative (robust but expensive) volume discretizations. Entropy indicators govern the selection, ensuring entropy-admissibility and stability without sacrificing performance (Doehring et al., 25 Mar 2026).
A posteriori subcell limiting: Employs subgrid finite-volume reconstructions when DG predictions fail physical or numerical admissibility criteria (e.g., positivity, maximum principle violation). This approach confines limiting to troubled cells, preserving high-order accuracy elsewhere (Boscheri et al., 2016, Fambri et al., 2018).
Slope and flux-corrected limiters: Embedded DG schemes inject partially continuous fields into a discontinuous space, apply bounded DG scheme and project back, ensuring $\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{u}) = \mathbf{S}(\mathbf{u}),$ 7 stability and local boundedness. Two-level limiting ensures no creation of new extrema during projection (Cotter et al., 2015).

DG schemes are unconditionally energy stable when formulated with appropriate time integrators and spatial discretization (e.g., penalty-free EQ-DG-RK for fourth-order flows, IEQ-DG for Cahn-Hilliard) (Liu et al., 2021, Liu et al., 2019). Entropy-stable DG construction based on summation-by-parts and entropy-conserving/dissipative fluxes ensure discrete entropy inequalities for nonlinear hyperbolic–mixed systems, with local or global time-stepping (Bhoriya et al., 2023).

6. Efficiency, Sparsity, and Scalability

DG schemes are especially suited for hardware-efficient, parallel computations due to:

Sparsity of global Jacobians and stiffness matrices: Line-based DG (Line-DG) achieves minimal per-DOF coupling ( $\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{u}) = \mathbf{S}(\mathbf{u}),$ 8 in 3D) versus standard nodal or DGSEM methods, leading to smaller memory footprints and faster solvers for implicit and steady-state computations (Persson, 2012).
Locality of updates: Both explicit/implicit and predictor–corrector one-step methods restrict updates to element-local or face-neighbor communication (ADER, ALE, VEM, MOOD), enabling strong scaling to exascale (Dumbser et al., 2018, Boscheri et al., 2016, Boscheri et al., 2023).
High arithmetic intensity and vectorization: Fine-grained cache optimization, on-the-fly data layout switching (AoS–SoA), and batch small-matrix operations accelerate computations, key for deep memory architectures (Dumbser et al., 2018).
Efficient parallel solvers: Schwarz multigrid preconditioners, matrix-free conjugate gradient, and hybrid face-based unknowns (HDG/VEM) enable robust linear/nonlinear solution even in extreme anisotropy or grid distortion (Vu, 2024).

DG schemes generalize seamlessly to unstructured, curved, and moving meshes (ALE-DG), remain stable and convergent under extreme stretching, and are appropriate for problems such as binary black hole initial data with stretching factors up to $\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{u}) = \mathbf{S}(\mathbf{u}),$ 9 (Vu, 2024).

7. Applications and Extensions

DG methods have established themselves as state-of-the-art for a wide spectrum of applications:

Compressible and incompressible flow, Euler/Navier-Stokes, MHD: High-order accuracy is critical for capturing vortex dynamics, turbulence, and acoustic phenomena without excessive diffusion.
General-relativistic hydrodynamics, two-fluid plasma, fusion modeling: Path-conservative ADER-DG and entropy-stable DG are leveraged for complex, coupled nonlinear systems (Fambri et al., 2018, Bhoriya et al., 2023).
Elliptic and mixed-type PDEs: Primal, hybridized, virtual and generalized spectral versions of DG address challenging elliptic systems even on highly stretched or anisotropic grids, with proven nearly-exponential decay in the MS-GFEM context (Alber et al., 24 Oct 2025).
Phase field models, gradient flows: Energy- and entropy-stable nonlinearly implicit DG formulations guarantee thermodynamic admissibility on nontrivial spatiotemporal domains (Liu et al., 2019, Liu et al., 2021).
Weather/climate, geophysical flows: Embedded DG and compatible DG (div–curl–grad) frameworks permit bounded, locally conservative transport of dependent variables in vertically continuous, horizontally discontinuous finite element spaces (Cotter et al., 2015, Abgrall et al., 31 Mar 2025).

Ongoing research continues to refine DG schemes in terms of adaptivity, entropy-stability, virtual element generalizations, exascale scalability, and efficient hybridizations, reinforcing DG's centrality in scientific computing.