Discontinuous Galerkin Methods

• Discontinuous Galerkin (DG) methods are finite element techniques characterized by the use of discontinuous, element-local polynomial approximations for solving partial differential equations.
• They employ numerical fluxes across element interfaces to ensure stability and accurately impose boundary conditions for various types of PDEs.
• DG methods enable high-order convergence, hp-adaptivity, and efficient parallel computing, making them ideal for complex simulations like hyperbolic, elliptic, and parabolic problems.

A discontinuous Galerkin (DG) method is a finite element approach for the discretization of partial differential equations (PDEs), characterized by the use of discontinuous, element-local trial and test spaces. This local discontinuity enables highly flexible, high-order, and locally conservative approximations for both linear and nonlinear PDEs across a wide spectrum of application areas, including hyperbolic, parabolic, and elliptic equations. The DG framework supports arbitrary unstructured meshes, nonconforming element types, hp-adaptivity, and is well-suited for parallel computing.

1. Formulation and Core Principles

In DG discretizations, the computational domain is partitioned into non-overlapping elements. Within each element, approximate solutions are represented by polynomials (typically of degree ≤k). These local approximations are not enforced to be continuous across element boundaries, in contrast to standard continuous Galerkin finite element methods. Instead, information between elements is communicated through numerical fluxes defined on the inter-element interfaces.

Consider a model hyperbolic system

$$\partial_t u + \nabla \cdot f(u) = 0, \quad x \in \Omega, \ t>0,$$

with $\Omega$ partitioned into elements $\{K\}$. The semi-discrete DG method seeks $u_h \in \prod_K P^k(K)$ such that for all test functions $\varphi_h$ of the same space,

$$\int_K \partial_t u_h \, \varphi_h \, dx - \int_K f(u_h) \cdot \nabla \varphi_h \, dx + \int_{\partial K} \hat{f}(u_h^{-}, u_h^{+}) \cdot n_K \varphi_h \, ds = 0,$$

where $\hat{f}$ is the numerical flux, and $u_h^{-}$, $u_h^{+}$ denote values from inside $K$ and from the adjacent element, respectively. The trace space, numerical flux construction, and interface jump conditions are central algorithmic components (Du et al., 2023, Fernandez et al., 2018).

2. Numerical Fluxes and Stability

The design of numerical fluxes at inter-element boundaries is critical for stability, accuracy, and the ability to impose boundary or interface conditions. For hyperbolic problems, canonical choices include upwind, Lax–Friedrichs, or Roe-type fluxes, which are consistent and dissipative. For instance, the upwind flux provides strict energy dissipation, whereas the central flux is energy-conserving. The flexibility of the DG framework allows the introduction of parameterized fluxes that control conservation and dissipation properties (Du et al., 2023, Duru et al., 2018). For physical problems such as elastic or electromagnetic waves, the flux may be constructed to exactly impose physical interface or boundary conditions and preserve specific conserved quantities.

The energy stability of the scheme is typically analyzed via a discrete energy functional,

$E^h(t) = \frac{1}{2} \sum_K \int_K |u_h|^2\, dx,$

with the evolution governed by interface fluxes. Under suitable flux choices and parameter constraints, one can rigorously guarantee non-increasing discrete energy, which is essential for long-time simulations and robustness (Du et al., 2023, Duru et al., 2018).

3. Error Analysis and Convergence

DG methods achieve arbitrarily high-order convergence, governed by the polynomial degree $k$ and smoothness of the exact solution. For linear and semilinear problems, the $L^2$-error for the standard upwind flux typically satisfies,

$\|u(t) - u_h(t)\|_{L^2} \leq C h^{k+1},$

given appropriate initial projection and under a global Lipschitz condition for any nonlinearities (Du et al., 2023, Duru et al., 2018). In practice, error analysis utilizes specialized projections (e.g., Gauss–Radau), energy arguments, and pseudo-time Gronwall inequalities. Superconvergence phenomena can appear for certain fluxes or polynomial degrees, particularly central fluxes in the even-degree setting. Error constants and convergence can depend on the geometric regularity of the mesh, with recent advances ensuring optimal rates even on highly anisotropic or degenerate meshes by judiciously modifying penalty terms (Kashiwabara et al., 2020).

4. Temporal Discretization and Compact Schemes

Time discretization in DG methods is typically achieved via explicit Runge–Kutta (RKDG), ADER, or Lax–Wendroff schemes, matching the spatial order. A notable advance is the compact Runge–Kutta DG (cRKDG) method, which replaces standard multi-stage coupling with a hybrid operator: inner stages use a local (element-wise) operator, and only the final update uses full DG coupling, thus maintaining a compact stencil regardless of RK stage number (Chen et al., 2023). This compactness reduces inter-element communication and is advantageous for large-scale parallel computation. Slope/WENO limiters, when required for shock capturing, can often be applied only at the final stage without compromising robustness.

5. Extensions: Hybridization, Local Approaches, and GPU Acceleration

Hybridized DG (HDG) and embedded DG (EDG/IEDG) methods reduce the number of globally coupled degrees of freedom by introducing face-based trace unknowns, leading to highly efficient static condensation and scalable solvers. HDG methods attain optimal or superconvergent rates for both primary variables and their gradients and are readily generalized to a broad class of PDEs, including wave propagation, Navier–Stokes, and Maxwell's equations (Fernandez et al., 2018). Specialized variants further reduce stencil sizes—such as line-based DG, which assembles the multidimensional problem from 1D DG operations along coordinate lines, yielding extremely sparse Jacobians essential for large implicit time steps or matrix-based solvers (Persson, 2012).

DG methods are also particularly well-suited to high-performance GPU acceleration. The elementwise independence enables efficient thread mapping, coalesced memory access, and hardware-specific autotuning, achieving high arithmetic intensity and bandwidth utilization in exascale computing environments (Chan et al., 2015, Klöckner et al., 2012, Fuhry et al., 2016).

6. Applications and Recent Advances

DG methods are widely used for nonlinear hyperbolic conservation laws, wave propagation, kinetic and Vlasov–Poisson systems, Maxwell and elastodynamic waves, convection-diffusion-reaction problems, and magnetohydrodynamics. Recent technical advances include:

• Energy-stable and Hamiltonian-conservative schemes for highly nonlinear or bi-Hamiltonian models, such as the Ostrovsky–Vakhnenko equation (Zhang et al., 2019).
• Semi-Lagrangian and local DG (LDG) strategies for convection-diffusion problems, enabling unconditionally stable, high-order, and locally conservative schemes even with large time steps (Ding et al., 2019, Boyana et al., 2024).
• Arbitrary-Lagrangian-Eulerian (ALE) DG with a posteriori subcell limiting for moving and deforming domains, vital for multiphysics applications with large deformations (Boscheri et al., 2016).
• Fourier continuation DG employing periodic spectral bases within the DG framework, yielding spectrally accurate and dispersion-minimizing approximations for wave and oscillatory problems (Appelo et al., 2021).

7. Summary Table: Core DG Scheme Ingredients

Ingredient	Function	Example References
Polynomial Space	Local, degree-$k$ discontinuous on each element	(Du et al., 2023, Persson, 2012)
Numerical Flux	Interface coupling for stability and boundary conditions	(Du et al., 2023, Duru et al., 2018)
Time Discretization	Explicit RKDG, ADER, or cRKDG for high-order in time	(Chen et al., 2023, Boscheri et al., 2016)
Limiting/Stabilization	Slope/WENO limiters, penalty terms, subcell FV limiting	(Boscheri et al., 2016, Ding et al., 2019)
Hybridization	Face-based trace unknowns, static condensation	(Fernandez et al., 2018)
Sparse Assembly	Line-DG/local 1D operators for minimal connectivity	(Persson, 2012)
GPU Implementation	Local kernels, batched memory layouts, autotuning	(Chan et al., 2015, Klöckner et al., 2012)

This comprehensive modularity, together with rigorous stability and error control properties, positions DG methods as a cornerstone technology for modern high-order simulation of PDEs in academic and industrial research (Du et al., 2023, Duru et al., 2018, Chen et al., 2023, Fernandez et al., 2018, Chan et al., 2015, Kashiwabara et al., 2020, Persson, 2012, Boyana et al., 2024, Boscheri et al., 2016, McLean, 2020, Ding et al., 2019, Fuhry et al., 2016, Feng et al., 2013, Klöckner et al., 2012, Mocz et al., 2013, Zhang et al., 2019, Appelo et al., 2021).