Fact-checked by Grok 2 weeks ago

Discontinuous Galerkin method

The Discontinuous Galerkin (DG) method is a class of numerical techniques for solving partial differential equations (PDEs), particularly laws, by approximating solutions with piecewise polynomial functions that are allowed to be discontinuous across element interfaces in a . This approach combines elements of finite volume and finite element methods, formulating weak solutions element-wise while incorporating numerical fluxes at interfaces to enforce and . Originating from efforts to handle transport equations, DG methods enable high-order accuracy and flexibility on unstructured grids without requiring continuous basis functions globally. The method was first introduced in 1973 by Reed and Hill for solving the neutron transport equation, a linear problem, using triangular meshes and piecewise constant approximations. Initial theoretical analysis followed in the 1970s by LeSaint and Raviart, who established error estimates for one-dimensional cases. Significant advancements occurred in the through the work of Cockburn and , who developed the Runge-Kutta DG (RKDG) framework by pairing DG spatial discretization with explicit Runge-Kutta time-stepping schemes, extending it to nonlinear time-dependent systems like the Euler equations. Variants such as the Local DG (LDG) method, introduced by Cockburn and in 1998, further adapted the approach for convection-diffusion problems by introducing auxiliary variables for higher-order derivatives. Key features of DG methods include their use of numerical fluxes—such as upwind or central fluxes—at element boundaries to couple discontinuous solutions, ensuring properties like local conservation and non-negativity for scalar equations. They achieve optimal L² error estimates of order k+1 for of degree k on smooth solutions, with provable via energy inequalities even for nonlinear conservation laws. Slope limiters, like the minmod limiter, are often applied to prevent oscillations near discontinuities while preserving high-order accuracy in smooth regions. These methods support -adaptivity, allowing variable degrees and refinements, and are highly parallelizable due to their local computations. DG methods find broad applications in for simulating compressible flows on complex geometries, electromagnetic wave propagation via , and seismic modeling on tetrahedral meshes. They excel in handling shocks and steep gradients in hyperbolic problems, such as and , and have been extended to elliptic and parabolic PDEs through interior penalty formulations. Ongoing developments focus on entropy-stable variants and asynchronous implementations for on GPUs.

Introduction

Definition and motivation

The discontinuous Galerkin (DG) method is a class of numerical methods for solving partial differential equations (PDEs) that utilizes discontinuous piecewise approximations, permitting the solution to exhibit discontinuities across element interfaces while employing high-order polynomials within individual elements. This approach integrates features of both finite volume and finite element methods, enforcing local conservation on each mesh element through carefully designed treatments. The motivation for adopting DG methods arises from their exceptional flexibility in managing unstructured meshes, arbitrary triangulations, and adaptive refinement strategies, including the accommodation of hanging nodes and varying polynomial degrees per element (hp-adaptivity). These properties make DG particularly well-suited for PDEs featuring shocks, discontinuities, or strong gradients, such as those in convection-dominated flows, where traditional continuous methods falter due to enforced global continuity constraints. DG methods offer several key advantages, including high-order accuracy, inherent local conservation, and enhanced parallelizability stemming from the locality of computations, which minimizes inter-processor communication. For problems, stability is bolstered by built-in upwinding mechanisms via numerical fluxes at interfaces. Compared to methods, DG provides superior handling of complex geometries without structured grid restrictions; relative to finite volume approaches, it achieves higher-order precision on irregular meshes without relying on intricate techniques; and in contrast to continuous finite element methods, it avoids global coupling, enabling more robust solutions for discontinuous phenomena.

Historical background

The Discontinuous Galerkin (DG) method originated in the early as a numerical technique for solving the equation, a arising in . In 1973, W. H. Reed and T. R. introduced the foundational approach, employing piecewise polynomial basis functions that are discontinuous across element interfaces on triangular meshes, coupled with upwind numerical fluxes to stabilize the and capture transport physics accurately. This innovation allowed for flexible, high-order approximations without enforcing global continuity, marking the first explicit use of discontinuous finite elements in a Galerkin framework for hyperbolic problems. Initial theoretical analysis followed in the 1970s by P. LeSaint and P. A. Raviart, who established error estimates for one-dimensional cases. Following this initial proposal, DG methods saw limited adoption during the and , partly due to computational constraints and a focus on continuous Galerkin and methods. A significant revival occurred in the late and 1990s, driven by Bernardo Cockburn and Chi-Wang Shu, who adapted and analyzed DG for time-dependent hyperbolic conservation laws. Their seminal series began with a 1989 paper developing total variation bounded (TVB) Runge-Kutta local projection DG schemes for scalar conservation laws, emphasizing explicit time integration and limiters for shock-capturing. Subsequent works in extended these Runge-Kutta DG (RKDG) methods to systems of conservation laws, establishing error estimates and demonstrating high-order accuracy on unstructured grids, which propelled DG into broader applications. The 2000s marked expansions of DG beyond problems, particularly to elliptic and parabolic equations. For diffusion-dominated phenomena, F. Bassi and S. Rebay proposed in 2000 a high-order DG for the compressible Navier-Stokes equations, introducing a compact of viscous fluxes via and interface corrections to ensure stability and accuracy. Concurrently, interior penalty DG (IPDG) methods gained prominence for second-order elliptic problems, refining earlier penalty-based ideas from the into symmetric and nonsymmetric variants analyzed for optimal convergence. By the late , hybridizable DG (HDG) variants emerged, pioneered by Cockburn, N. C. Nguyen, and J. Peraire in 2009 for convection-diffusion systems, enabling efficient postprocessing to continuous solutions and reducing through hybridization. These advancements drew inspiration from finite volume methods' flux handling and spectral elements' high-order bases, fostering DG's versatility across PDE types.

Fundamental Concepts

Discontinuous function spaces

In the discontinuous Galerkin (DG) method, the trial and test functions are sought within a finite-dimensional subspace consisting of piecewise polynomials that are discontinuous across element interfaces. Specifically, let \mathcal{T}_h denote a simplicial or tensor-product mesh partitioning the computational domain \Omega into elements K, and let \mathcal{P}^p(K) be the space of polynomials of degree at most p on each element K. The discontinuous function space is then defined as the broken Sobolev space V_h = \{ v_h \in L^2(\Omega) : v_h|_K \in \mathcal{P}^p(K) \ \forall K \in \mathcal{T}_h \}, where no inter-element continuity is imposed. This construction allows functions in V_h to exhibit arbitrary jumps at interfaces, enabling the method to handle solutions with discontinuities or sharp gradients without enforcing global regularity. The discontinuities in V_h are quantified by the jump operator across interior interfaces e \in \mathcal{E}_h, defined for a v_h \in V_h as [[v_h]] = v_h^+ - v_h^- on e, where v_h^+ and v_h^- denote the limiting values from the adjacent elements sharing e. Similarly, trace operators extract the values of v_h on element boundaries \partial K, facilitating the evaluation of numerical fluxes at interfaces. To initialize approximations, the L^2-orthogonal projection \Pi: L^2(\Omega) \to V_h is commonly employed, satisfying \int_K (\Pi u - u) v_h \, dx = 0 for all v_h \in \mathcal{P}^p(K) and each K \in \mathcal{T}_h. These operators ensure that initial data or source terms can be consistently projected onto the discontinuous space while preserving key properties like mass conservation in the discrete setting. For and , a natural on V_h incorporates both and contributions, known as the DG : \|v_h\|_{DG}^2 = \sum_{K \in \mathcal{T}_h} \|v_h\|_{L^2(K)}^2 + \sum_{e \in \mathcal{E}_h} h_e^{-1} \| [[v_h]] \|_{L^2(e)}^2, where h_e is the (or area) of e. This controls the jumps explicitly, providing a measure of the solution's regularity across elements and enabling estimates in the bilinear forms used by DG methods. terms may be added for incomplete meshes, but the interior structure emphasizes the role of discontinuities in maintaining method .

Numerical fluxes and interface conditions

In discontinuous Galerkin (DG) methods, numerical fluxes play a crucial role in coupling the solutions across element interfaces where the approximate solution u_h is discontinuous. The numerical flux \hat{f}(u_h^-, u_h^+, \mathbf{n}) is defined as an approximation to the physical flux f(u) \cdot \mathbf{n} at the interface e between adjacent elements, depending on the traces u_h^- from the interior of one element and u_h^+ from the neighboring element, as well as the outward unit normal \mathbf{n}. These traces arise from the discontinuous piecewise polynomial function spaces used in DG formulations. For hyperbolic problems, several common numerical fluxes are employed to ensure stability and accuracy. The upwind flux selects the flux based on the wind direction: for a scalar equation with flux f(u) and characteristic speed \lambda, it is given by \hat{f}(u^-, u^+, \mathbf{n}) = f(u^-) \cdot \mathbf{n} if \lambda > 0, and f(u^+) \cdot \mathbf{n} otherwise, mimicking the direction of information propagation. The central flux averages the states, \hat{f}(u^-, u^+, \mathbf{n}) = \frac{1}{2} [f(u^-) + f(u^+)] \cdot \mathbf{n}, providing less dissipation but potentially requiring stabilization. The Lax-Friedrichs flux adds artificial viscosity, \hat{f}(u^-, u^+, \mathbf{n}) = \frac{1}{2} [f(u^-) + f(u^+)] \cdot \mathbf{n} - \frac{\alpha}{2} (u^+ - u^-), where \alpha is a scaling factor like the maximum wave speed. These fluxes were adapted from finite volume methods into the DG framework in seminal works on high-order schemes for conservation laws. Interface conditions in DG methods are enforced weakly through integration by parts in the variational formulation over each element K. This yields boundary terms on the interfaces e \subset \partial K, resulting in conditions of the form \int_e (\hat{f}(u_h^-, u_h^+, \mathbf{n}) - f(u_h) \cdot \mathbf{n}) v_h \, ds = 0 for test functions v_h in the discontinuous space, ensuring the numerical flux replaces the exact flux to handle discontinuities without requiring continuity constraints. Key properties of numerical fluxes include and . Consistency requires that \hat{f}(u, u, \mathbf{n}) = f(u) \cdot \mathbf{n} when the traces agree, allowing the scheme to recover the exact solution for smooth problems. is achieved by defining a single-valued flux at each , such that the flux outgoing from one equals the incoming flux to the adjacent , preserving the global balance of the . These properties, essential for the method's reliability, were established in the foundational DG developments for systems.

Formulation for Hyperbolic Problems

Scalar conservation laws

The discontinuous Galerkin (DG) method provides a framework for numerically solving scalar laws, serving as a foundational example for hyperbolic problems. Consider the one-dimensional model \partial_t u + \partial_x f(u) = 0, \quad x \in \Omega, \ t > 0, equipped with initial condition u(x,0) = u_0(x) and suitable boundary conditions on the domain \Omega. This models phenomena such as inviscid fluid flow or , where u represents a and f(u) is the function. To derive the DG formulation, the \Omega into a of non-overlapping K_i = [x_{i-1/2}, x_{i+1/2}] for i = 1, \dots, N, with mesh size h_i = x_{i+1/2} - x_{i-1/2}. The approximate solution u_h(x,t) belongs to the discontinuous finite element consisting of polynomials of at most p on each K_i, allowing discontinuities across element interfaces. The semi-discrete scheme seeks u_h such that, for each element K_i and all test functions v_h \in \mathcal{P}^p(K_i), \int_{K_i} \partial_t u_h \, v_h \, dx - \int_{K_i} f(u_h) \partial_x v_h \, dx + \hat{f}_{i+1/2} \, v_h(x_{i+1/2}^-) - \hat{f}_{i-1/2} \, v_h(x_{i-1/2}^+) = 0, where \hat{f} denotes a single-valued numerical flux approximating f(u) at interfaces, ensuring conservation and stability. This weak form integrates the conservation law by parts on each element, replacing the exact flux at boundaries with the numerical flux \hat{f} to couple adjacent elements. On each element K_i, expand u_h(x,t) = \sum_{j=1}^{p+1} u_j^i(t) \, \phi_j(x), where \{\phi_j\}_{j=1}^{p+1} form a basis for \mathcal{P}^p(K_i), such as shifted . Substituting this expansion and the test functions v_h = \phi_k for k = 1, \dots, p+1 yields a system of ordinary differential equations (ODEs) for the coefficients \mathbf{u}^i(t) = (u_1^i(t), \dots, u_{p+1}^i(t))^T: M^i \, \dot{\mathbf{u}}^i(t) = \mathbf{r}^i(\mathbf{u}(t)), where M^i_{kj} = \int_{K_i} \phi_k \, \phi_j \, dx is the local (diagonal in the orthogonal basis case), and the residual vector \mathbf{r}^i incorporates the volume integrals -\int_{K_i} f(u_h) \partial_x \phi_k \, dx and the interface contributions from \hat{f}. This ODE system must be solved for each , with \hat{f} depending on values from neighboring elements. A prototypical example is linear , where f(u) = a u with constant speed a > 0. The upwind numerical flux is then \hat{f}_{i+1/2} = a \, u_h(x_{i+1/2}^-), selecting the solution from the left (upwind) side to enforce . This choice ensures the scheme is monotone and L^2-stable under a suitable Courant-Friedrichs-Lewy (CFL) condition.

Systems of conservation laws

The discontinuous Galerkin (DG) method extends naturally to systems of nonlinear conservation laws in multiple spatial dimensions, which model phenomena such as compressible fluid flows. The governing equations take the form \partial_t \mathbf{u} + \nabla \cdot \mathbf{F}(\mathbf{u}) = 0, where \mathbf{u} \in \mathbb{R}^m is the vector of conserved variables, \mathbf{F}: \mathbb{R}^m \to \mathbb{R}^{d \times m} is the tensor with d the spatial dimension, and the equation holds over a \Omega \subset \mathbb{R}^d with appropriate initial and boundary conditions. This formulation generalizes scalar conservation laws by allowing vector-valued states and nonlinear dependencies, capturing wave propagation in systems like gas dynamics. To derive the DG scheme, multiply the conservation law by a smooth vector test function \mathbf{v}_h and integrate by parts over each element K of a tessellation \mathcal{T}_h of \Omega: \int_K \partial_t \mathbf{u}_h \cdot \mathbf{v}_h \, d\mathbf{x} - \int_K \mathbf{F}(\mathbf{u}_h) : \nabla \mathbf{v}_h \, d\mathbf{x} + \int_{\partial K} \hat{\mathbf{F}} \cdot \mathbf{v}_h \, ds = 0, where \mathbf{u}_h is a piecewise polynomial approximation in a discontinuous finite element space, \hat{\mathbf{F}} is a single-valued numerical flux consistent with \mathbf{F}, and the colon denotes the double contraction. The scheme is advanced in time using explicit Runge-Kutta methods, ensuring total variation diminishing stability under a CFL condition. This element-local formulation allows high-order accuracy while handling discontinuities through the choice of \hat{\mathbf{F}}. For systems, the numerical flux \hat{\mathbf{F}} must resolve multiple wave families, typically via approximate Riemann solvers. Common choices include the Roe solver, which linearizes the flux at an upwind-biased average state to compute the interface flux exactly for the linearized problem, and the HLLC solver, which approximates the Riemann fan with two intermediate waves to better capture contact discontinuities. These solvers, originally developed for finite volume methods, integrate seamlessly into DG by evaluating at element interfaces using left and right limit states. A representative application is the two-dimensional Euler equations for inviscid , \partial_t \begin{pmatrix} \rho \\ \rho \mathbf{v} \\ E \end{pmatrix} + \nabla \cdot \begin{pmatrix} \rho \mathbf{v} \\ \rho \mathbf{v} \otimes \mathbf{v} + p \mathbf{I} \\ (E + p) \mathbf{v} \end{pmatrix} = 0, with \rho > 0, \mathbf{v}, total energy E, and p from an . To prevent oscillations near shocks, limiters project the solution onto characteristic variables using the local eigenstructure of the flux , applying scalar TVD or TVB limiters componentwise before reconstruction. This characteristic decomposition enhances robustness for nonlinear systems while preserving high-order accuracy in smooth regions.

Formulation for Elliptic and Parabolic Problems

Interior penalty methods

The interior penalty discontinuous Galerkin (IPDG) methods form a class of formulations designed for second-order elliptic problems, enforcing weakly across interfaces through penalty terms and flux averaging. These methods utilize discontinuous piecewise spaces and yield symmetric positive definite systems when the bilinear form is symmetric, facilitating efficient solution strategies. Originating from early work on mixed methods and evolving into a unified , IPDG approaches balance accuracy and for diffusion-dominated problems. Consider the model second-order elliptic problem -\nabla \cdot (a \nabla u) + c u = f in a bounded domain \Omega \subset \mathbb{R}^d with a(x) a symmetric positive definite matrix, c \geq 0, and suitable boundary conditions, such as Dirichlet u = g on \partial \Omega. Let \mathcal{T}_h be a simplicial mesh of \Omega with elements K, edges e, and mesh size h = \max_K \operatorname{diam}(K). The IPDG method seeks u_h \in V_h = \{ v \in L^2(\Omega) : v|_K \in \mathbb{P}^p(K), \forall K \in \mathcal{T}_h \}, where \mathbb{P}^p(K) denotes polynomials of degree at most p, satisfying the variational formulation: find u_h \in V_h such that B(u_h, v_h) = \ell(v_h) for all v_h \in V_h, where \ell(v_h) incorporates the right-hand side and boundary data. The for the IPDG method is given by \begin{split} B(u_h, v_h) &= \sum_{K \in \mathcal{T}_h} \int_K a \nabla u_h \cdot \nabla v_h \, dx \\ &\quad - \sum_{e \in \mathcal{E}_h^i \cup \mathcal{E}_h^b} \int_e \{ a \nabla u_h \cdot \mathbf{n} \} [[v_h]] \, ds - \sum_{e \in \mathcal{E}_h^i \cup \mathcal{E}_h^b} \int_e \{ a \nabla v_h \cdot \mathbf{n} \} [[u_h]] \, ds \\ &\quad + \sum_{e \in \mathcal{E}_h^i \cup \mathcal{E}_h^b} \frac{\sigma}{h_e} \int_e [[u_h]] [[v_h]] \, ds, \end{split} where \mathcal{E}_h^i and \mathcal{E}_h^b denote interior and edges, respectively; \{ \cdot \} is the ; [[\cdot]] is the jump ; \mathbf{n} is a unit normal; h_e is the edge ; and \sigma > 0 is the penalty parameter. This formulation adapts numerical fluxes for diffusive terms by averaging gradients across interfaces while penalizing jumps to control discontinuities. Several variants of the IPDG method arise from modifications to the interface terms. The symmetric interior penalty Galerkin (SIPG) method employs the full symmetric form above with positive \sigma, ensuring B(\cdot, \cdot) is symmetric and coercive under suitable \sigma. The incomplete interior penalty Galerkin (IIPG) omits the second interface term \int_e \{ a \nabla v_h \cdot \mathbf{n} \} [[u_h]] \, ds, reducing computational cost but maintaining stability for certain problems. The nonsymmetric interior penalty Galerkin (NIPG) reverses the sign of the second interface term, leading to a nonsymmetric that may offer advantages in conditioning but lacks certain consistency properties. The penalty parameter \sigma is crucial for stability and convergence; it must be chosen sufficiently large, typically scaling with the polynomial degree p raised to a power depending on the dimension d (e.g., O(p^2) in 2D and O(p^3) in 3D), and proportional to \max_{\Omega} \|a\|_\infty, based on inverse trace inequalities bounding jumps on element boundaries. Smaller \sigma risks instability, while excessively large values increase stiffness without improving accuracy. Computable lower bounds for \sigma can be obtained via eigenvalue problems on reference elements, ensuring robustness across mesh refinements and polynomial degrees. The symmetric IPDG variants (SIPG and IIPG) exhibit both consistency and adjoint consistency, essential for optimal L^2-error estimates of order O(h^{p+1}). Consistency follows from the Galerkin orthogonality property, where the exact solution satisfies the variational form weakly due to and consistent numerical fluxes. Adjoint consistency, preserved by the symmetric interface terms, ensures the error satisfies a problem , enabling higher-order without suboptimal factors. In contrast, the NIPG variant is consistent but not adjoint consistent, potentially leading to reduced L^2-rates.

Local discontinuous Galerkin methods

The local discontinuous Galerkin (LDG) method addresses elliptic and parabolic problems by reformulating them as mixed first-order systems, enabling the use of alternating numerical fluxes to handle both diffusive and convective terms effectively. This approach is particularly suited for convection-diffusion equations of the form -\epsilon \Delta u + \mathbf{b} \cdot \nabla u = f in a \Omega, where \epsilon > 0 is the and \mathbf{b} is the . To apply the LDG framework, the equation is rewritten as a system: \mathbf{q} = -\epsilon \nabla u and \nabla \cdot (\mathbf{q} + \mathbf{b} u) = f, introducing the auxiliary flux variable \mathbf{q} to approximate the separately. In the LDG scheme, independent approximations u_h and \mathbf{q}_h are sought in a discontinuous finite element space V_h, typically consisting of piecewise polynomials of degree k on a triangulation of \Omega. The method employs numerical fluxes at element interfaces to ensure weak continuity and stability. Alternating fluxes are used: the flux for u is taken from the "interior" side, \hat{u} = u_h^-, while the normal component of the flux for \mathbf{q} incorporates the opposite side with a penalty term, \hat{\mathbf{q}} \cdot \mathbf{n} = \mathbf{q}_h^+ \cdot \mathbf{n} + \frac{\epsilon}{h} [[u_h]], where h is the local mesh size and [[\cdot]] denotes the jump across the interface (with directions defined relative to a fixed orientation or unified globally). These choices promote upwinding for convection and central differencing for diffusion, adjusted by the penalty to control oscillations. The weak formulation proceeds elementwise. For the u-equation on each element K, integrate by parts to obtain \int_K (\mathbf{q}_h - \mathbf{b} u_h) \cdot \nabla v_h \, dx + \int_{\partial K} \hat{\mathbf{h}} \cdot \mathbf{n} \, v_h \, ds = \int_K f v_h \, dx for all test functions v_h \in V_h, where \hat{\mathbf{h}} is the combined numerical flux incorporating convective and diffusive contributions (with appropriate upwind for \mathbf{b} \cdot \mathbf{n} > 0). For the \mathbf{q}-equation, the form is \int_K \mathbf{q}_h \cdot \mathbf{r}_h \, dx + \epsilon \int_K u_h \nabla \cdot \mathbf{r}_h \, dx - \epsilon \int_{\partial K} \hat{u} \, \mathbf{r}_h \cdot \mathbf{n} \, ds = 0 for test functions \mathbf{r}_h \in [V_h]^d. Boundary conditions are imposed weakly through appropriate flux modifications on \partial \Omega. This mixed formulation yields L^2-stability for linear problems and optimal-order convergence in L^2 norm using degree-k polynomials. LDG methods excel in convection-dominated regimes (\epsilon \ll \|\mathbf{b}\|) by naturally incorporating upwind stabilization via the alternating fluxes, reducing non-physical oscillations compared to methods while maintaining high-order accuracy and local conservation properties. Their element-local structure also facilitates efficient parallel implementation and adaptation to unstructured meshes.

Advanced Variants

Direct discontinuous Galerkin method

The direct discontinuous Galerkin (DDG) method provides a framework for discretizing second-order partial differential equations within the discontinuous Galerkin paradigm, particularly suited for elliptic and parabolic problems, by directly handling the second-order terms without introducing auxiliary variables for fluxes. Originally developed for problems, the method extends naturally to more general second-order operators through a that incorporates carefully defined numerical traces at element s to ensure consistency, conservation, and stability. This approach contrasts with methods like the local discontinuous Galerkin (LDG) by avoiding mixed formulations and instead relying on interface corrections to approximate derivatives. The method applies to the general second-order elliptic equation -\nabla \cdot (A \nabla u) + \mathbf{b} \cdot \nabla u + c u = f on a bounded domain \Omega \subset \mathbb{R}^d, where A is a symmetric positive definite diffusion tensor, \mathbf{b} is the convection vector, c \geq 0 is the reaction coefficient, and f is the source term, subject to suitable boundary conditions (e.g., Dirichlet or Neumann). The discrete solution u_h is sought in a discontinuous piecewise polynomial space V_h of degree k on a simplicial mesh. The semi-discrete DDG formulation is obtained by integrating by parts the diffusion and convection terms within each element K, yielding volume integrals over K and boundary integrals over faces e \in \partial K. Specifically, for test functions v_h \in V_h, the scheme reads \sum_K \int_K \left( A \nabla u_h \cdot \nabla v_h + \mathbf{b} \cdot \nabla u_h \, v_h + c u_h v_h - f v_h \right) \, dx + \sum_e \int_e \left( \widehat{A \nabla u_h} \cdot \mathbf{n} \, [[v_h]] \right) \, ds = 0, where \mathbf{n} is the unit normal, [[\cdot]] denotes the jump across e, \{\cdot\} the average, and boundary terms are adjusted accordingly (e.g., \widehat{A \nabla u_h} \cdot \mathbf{n} = - (A \nabla u_h \cdot \mathbf{n}) for outflow boundaries). The convection term employs standard numerical fluxes, such as upwind or central approximations consistent with \mathbf{b} \cdot \mathbf{n}, while the reaction term remains local. Stability is enforced through parameters in the numerical fluxes chosen sufficiently large (e.g., \beta_0 > C for some constant C > 0) to control jumps. Central to the DDG method is the numerical diffusive flux \widehat{A \nabla u_h} = A \tilde{\nabla} u_h, where \tilde{\nabla} u_h is a consistent numerical gradient defined elementwise as \tilde{\nabla} u_h|_K = \nabla u_h + \sum_{e \in \partial K} \mathbf{l}_e \frac{[[u_h]]}{h_e} (with \mathbf{l}_e lifting operators), but at interfaces, it takes the form \tilde{\nabla} u_h = \{\nabla u_h\} + C_{11} \frac{[[u_h]] \mathbf{n}}{h_e} + C_{12} [[ \nabla u_h \cdot \mathbf{n} ]] \mathbf{n}, with stabilization parameters C_{11}, C_{12} (e.g., C_{11} = O(1), C_{12} = O(1/h_e) for higher-order corrections if k \geq 2). These correction terms ensure consistency by recovering the exact flux for smooth solutions and provide stability analogous to penalty methods, though without explicit symmetrization. For the convection part, Prager-Synge-type fluxes—drawing from hypercircle principles for bounding errors via complementary energy—can be employed to enhance robustness, particularly in convection-dominated regimes, by defining monotone or centered traces that minimize dissipation while maintaining L^2-stability. Numerical experiments demonstrate optimal convergence rates of O(h^{k+1}) in the L^2-norm for smooth solutions. Unlike interior penalty discontinuous Galerkin (IPDG) methods, which enforce weakly through symmetric bilinear forms with explicit penalties on both solution jumps [[u_h]] and gradient inconsistencies \{\nabla u_h\} \cdot \mathbf{n}, the DDG approach directly discretizes the second-order operator via the numerical , resulting in a potentially non-symmetric form that simplifies implementation for nonlinear or variable-coefficient problems. This direct treatment reduces computational overhead compared to mixed methods and has been shown to yield superconvergence at certain points (e.g., downwind points in 1D) under specific choices. The method's efficacy is illustrated in applications like convection-diffusion equations, where it preserves maximum principles with appropriate limiters and achieves high-order accuracy on unstructured meshes.

Hybridizable discontinuous Galerkin

The hybridizable discontinuous Galerkin (HDG) introduces a hybrid variable \hat{u}_h defined on the interfaces between elements to enhance computational efficiency over standard discontinuous Galerkin approaches. This variable represents the of the solution on element boundaries, allowing the formulation to decouple the problem into local solves within each element K, which are then coupled globally only through the traces. By enforcing of the hybrid variable across interfaces, the method reduces the global system size significantly while maintaining the flexibility of discontinuous approximations inside elements. This hybridization builds on precursors like the local discontinuous Galerkin but achieves better scalability for higher-order polynomials. For diffusion problems governed by -\nabla \cdot (\kappa \nabla u) = f, the HDG approximates both the scalar u_h and the \mathbf{q}_h = -\kappa \nabla u_h using discontinuous piecewise polynomials of degree p within each element K. The local problems enforce the constitutive relation \mathbf{q}_h + \kappa \nabla u_h = 0 in K and the balance equation \nabla \cdot \mathbf{q}_h = f in K, with the condition u_h - \hat{u}_h = O(h^{p+1}) imposed weakly on the \partial K via a stabilization term. Globally, of the normal \mathbf{q}_h \cdot \mathbf{n} across interfaces is ensured through the hybrid variable, leading to a conservative scheme that couples only the traces. In the HDG scheme, solving the local element problems yields an explicit expression for the interior solution: u_h = \Pi_K \left( \hat{u}_h - \tau^{-1} \mathbf{q}_h \cdot \mathbf{n} \right), where \Pi_K is the L^2- onto the space of degree p+1 in K, and \tau > 0 is a stabilization . This allows static condensation of the local , resulting in a reduced global system solved solely for the hybrid traces \hat{u}_h on the mesh skeleton. The resulting is symmetric positive definite for appropriate \tau and sparse, facilitating efficient iterative solvers. A key feature of HDG methods is their superconvergence : after solving for \hat{u}_h, a simple local post-processing step recovers an enhanced approximation \tilde{u}_h of degree p+1 that converges at O(h^{p+2}) in the L^2-norm, one order higher than the standard O(h^{p+1}) for u_h. This post-processing solves a local Neumann problem using the computed flux and traces, providing higher accuracy without additional global cost. The primary advantages of HDG include a drastic reduction in global unknowns—by approximately one compared to standard DG for polynomial degree p, as only face dofs are globally coupled—making it particularly suitable for implicit time-stepping in time-dependent problems and high-order simulations. This efficiency stems from the hybridization, which minimizes communication in implementations while preserving optimal rates.

Theoretical Aspects

Stability analysis

The stability of discontinuous Galerkin (DG) methods for problems is often established through estimates in the L^2 norm, demonstrating that the numerical solution remains bounded by its initial data. For the semi-discrete DG formulation applied to linear equations with periodic or inflow conditions, the scheme satisfies an of the form \frac{d}{dt} \|u_h\|^2 + \sum_e \int_e (f(u_h^+) - f(u_h^-)) \cdot \mathbf{n} \, ds \leq 0, where u_h is the DG approximation, f is the function, and the sum is over all interfaces e. This estimate arises from integrating the DG weak form against u_h itself and using within each element, with the interface terms providing dissipation when upwind numerical fluxes are employed, such as the Godunov flux for scalar cases or Roe's flux for systems. The upwind choice ensures that the jump terms (f(u_h^+) - f(u_h^-)) \cdot \mathbf{n} are non-positive, leading to \frac{d}{dt} \|u_h\|^2 \leq 0 and thus \|u_h(\cdot, t)\|_{L^2} \leq \|u_h(\cdot, 0)\|_{L^2} for all t > 0. This framework, originally introduced for and extended to general systems, underpins the L^2- of Runge-Kutta DG (RKDG) methods. For elliptic problems, stability in interior penalty DG (IPDG) methods is analyzed through the coercivity of the bilinear form defining the discrete problem. The symmetric IPDG formulation for the Poisson equation -\Delta u = g yields a bilinear form B(v_h, w_h) that includes volume terms \int_\Omega \nabla_h v_h \cdot \nabla_h w_h \, dx, antisymmetric interface terms -\int_\Gamma (\llbracket v_h \rrbracket \cdot \{\nabla_h w_h\} + \{\nabla_h v_h\} \cdot \llbracket w_h \rrbracket) \, ds, and a symmetric penalty term \sum_e \frac{\sigma_e}{h_e} \int_e \llbracket v_h \rrbracket \cdot \llbracket w_h \rrbracket \, ds, where \nabla_h and \{\cdot\}, \llbracket \cdot \rrbracket denote broken gradient, average, and jump operators, respectively, \sigma_e > 0 is the penalty parameter, and h_e is the element size. Coercivity holds when \sigma_e is sufficiently large, satisfying B(v_h, v_h) \geq \alpha \|v_h\|_{DG}^2 for all v_h in the DG space, where the DG norm is \|v_h\|_{DG}^2 = \|\nabla_h v_h\|_{L^2}^2 + \sum_e \frac{\sigma_e}{h_e} \|\llbracket v_h \rrbracket\|_{L^2(\partial K_e)}^2 and \alpha > 0 depends on the minimal \sigma_e (typically \alpha \sim \min(\sigma_e, 1)). This bound ensures well-posedness and stability independent of mesh size, as established in the unified analysis of DG methods for second-order elliptic problems. Entropy stability for nonlinear scalar conservation laws extends these ideas by adapting Kružkov entropy inequalities to the DG framework, ensuring that the numerical solution respects physical dissipation. For the DG applied to u_t + \nabla \cdot f(u) = 0, an pair (\eta(u), q(u)) with \eta'' > 0 and \nabla q = f' \nabla \eta satisfies a discrete \frac{d}{dt} \int_\Omega \eta(u_h) \, dx + \sum_e \int_e (\hat{q}^+ - \hat{q}^-) \cdot \mathbf{n} \, ds \leq 0, where \hat{q} is an entropy-consistent numerical flux (e.g., satisfying \hat{q}(v,v) = q(v) and upwind-biased for ). For the Kružkov \eta_k(u) = |u - k| and corresponding q_k, this yields cell-wise bounds that propagate to global total variation diminishing (TVD) properties or -boundedness, preventing non-physical oscillations in shock-capturing simulations. This adaptation, building on finite volume stability, has been rigorously proven for high-order DG schemes using summation-by-parts operators to mimic exactly in the discrete sense. Seminal developments include -stable flux formulations that ensure the for convex entropies, applicable to both scalar and Euler systems. For linear hyperbolic problems, von Neumann analysis provides insight into by examining eigenvalue bounds in Fourier space, revealing and properties of DG schemes. Assuming a uniform mesh and solutions u_h(x,t) = \hat{u}(t) e^{i \kappa x}, the semi-discrete DG yields a eigenvalue problem whose must lie in the left half-plane for under explicit time-stepping (CFL \Delta t \leq C h / \lambda_{\max}, where \lambda_{\max} is the ). For polynomial degree p, the eigenvalues satisfy |\lambda| \leq 1 in the upwind case, with increasing for higher modes, though central fluxes may require stabilization. This analysis confirms L^2- for RKDG methods and guides parameter choices for minimal in smooth regions. A notable limitation of high-order DG methods is their failure to satisfy a discrete without additional limiters, particularly for convection-dominated problems. Unlike low-order schemes, high-order polynomials within elements can produce overshoots or undershoots near discontinuities, violating bounds like \min u_0 \leq u_h \leq \max u_0 even with upwind fluxes, as the local approximation lacks monotonicity. This issue arises because the DG solution is not necessarily sign-preserving or bound-preserving in the L^\infty sense for p \geq 1, necessitating slope/flux limiters to enforce the principle while maintaining high-order accuracy in smooth flows.

Error estimates and convergence

Error estimates for discontinuous Galerkin (DG) methods play a crucial role in assessing their accuracy and guiding refinement strategies. For problems, such as scalar laws, the standard a priori yields a suboptimal rate in the L^2 norm. Specifically, the satisfies \|u - u_h\|_{L^2} \leq C h^{p + 1/2} \|u\|_{H^{p+1}(\Omega)}, where u is the exact solution, u_h is the DG approximation using polynomials of degree p, h is the maximum element diameter, and C is a constant independent of h but depending on the regularity of u and the geometry. This estimate arises from the of Runge-Kutta DG schemes and holds under conditions derived from estimates. To achieve optimal , post-processing techniques can be applied to the DG solution, recovering an bound of O(h^{p+1}) in the L^2 norm. These methods involve local projections onto higher-degree polynomials, leveraging superconvergence properties at certain points within elements. For elliptic problems, interior penalty DG (IPDG) methods provide optimal-order a priori error estimates in the DG norm, which incorporates both the broken and interface jumps. The error bound is \|u - u_h\|_{\rm DG} \leq C h^p \|u\|_{H^{p+1}(\Omega)}, where the DG norm is defined as \|v\|_{\rm DG}^2 = \|\nabla_h v\|_{L^2(\Omega)}^2 + \sum_e \frac{\sigma}{h_e} \|\llbracket v \rrbracket\|_{L^2(e)}^2, with e denoting s, \nabla_h the broken , \llbracket \cdot \rrbracket the jump operator, and \sigma > 0 the penalty . The constant C remains bounded independently of p provided the penalty parameter \sigma is fixed and sufficiently large relative to the coefficient. This optimality holds across various IPDG variants, including symmetric, incomplete, and nonsymmetric forms, and extends to convection-dominated regimes with appropriate stabilization. A posteriori error estimates enable adaptive refinement by providing computable indicators of local error contributions. Residual-based estimators are commonly employed, decomposing the global error into element-wise terms such as \eta_K = h_K \|r(u_h)\|_{L^2(K)} + \left( h_K^{1/2} \|[[ \mathbf{q}_h + \mathbf{n} \cdot \nabla u_h ]]\|_{L^2(\partial K)} + h_K^{-1/2} \|[[u_h]]\|_{L^2(\partial K)} + \cdots \right), where r(u_h) is the interior , \mathbf{q}_h is the numerical for the , and additional terms account for contributions. These estimators are reliable and efficient, satisfying c \|\eta\| \leq \|u - u_h\|_{\rm DG} \leq C \|\eta\| for positive constants c, C, and support goal-oriented adaptations for both elliptic and problems. In hybridizable DG (HDG) methods, convergence properties exhibit distinct behaviors between traces and volumetric variables. The traces on element interfaces achieve optimal O(h^{p+1}) convergence in the L^2 norm, while the solution within elements demonstrates superconvergence, often reaching O(h^{p+3/2}) or higher under structured meshes and specific flux choices. This separation facilitates efficient solvers via hybridization and enables post-processing to elevate the overall accuracy to O(h^{p+2}). Recent advancements since 2010 have established unified frameworks for a priori estimates in convection-diffusion problems, accommodating both diffusion-dominated and limits through parameter-robust analyses. These frameworks extend to hp-versions, yielding quasi-optimal rates O(N^{-(p+1)/d}) in d dimensions, where N is the , and incorporate adaptivity for singularly perturbed cases. Such results build on stability analyses to ensure robustness across degrees and mesh refinements.

Applications

Fluid dynamics

The discontinuous Galerkin (DG) method has been extensively applied to the compressible Navier-Stokes equations, particularly for simulating flows involving shocks and , where high-order accuracy is essential for resolving multi-scale phenomena. In these applications, DG formulations incorporate slope such as the total variation bounded (TVB) minmod to prevent oscillations near discontinuities while maintaining high-order accuracy in regions. For instance, the TVB , parameterized to balance resolution and stability, has been shown to effectively capture shock-turbulence interactions in two-dimensional simulations of conservation laws, extending naturally to the Navier-Stokes system. This approach allows DG to handle the hyperbolic-parabolic nature of the equations without excessive numerical dissipation, enabling robust solutions for and supersonic regimes. High-order DG methods excel in aeroacoustics and detonation wave simulations, where precise capture of wave propagation and nonlinear structures is critical. In , DG schemes on unstructured meshes have been used to model noise generation from airfoil trailing edges, achieving low-dissipation propagation of over long distances. For waves, Runge-Kutta DG (RKDG) methods resolve unstable cellular structures in two- and three-dimensional reactive flows, demonstrating superior performance in detailing interactions compared to lower-order methods. A representative example is the simulation of turbulent flow over a three-dimensional NACA 0012 at 6 million, where high-order DG on high-density unstructured meshes captures separation and wake vortices with minimal artificial dissipation. Key advantages of DG in include its adaptability to unstructured meshes, which facilitates simulations of complex industrial geometries like turbine blades or fuselages, and its suitability for implicit large-eddy simulation () without explicit subgrid-scale models. The method's inherent upwinding and local conservation properties act as a natural filter for under-resolved , providing accurate predictions of energy spectra in transitional flows. For example, hybridized DG variants have demonstrated -consistent implicit LES for turbulent boundary layers, achieving grid-independent results on meshes with O(10^6) elements. However, challenges arise in ensuring stability for under-resolved turbulent flows, where traditional DG schemes may violate physical inequalities, leading to instabilities. Recent developments since 2015, such as -stable DG formulations using summation-by-parts operators and -conserving fluxes, address this by guaranteeing discrete dissipation, improving robustness in high-Mach-number simulations. Implementations of DG for (CFD) are supported by open-source libraries like deal.II and FEniCS, which provide flexible frameworks for high-order discretizations. The deal.II library, through extensions like exaDG, enables efficient parallel simulations of the compressible Navier-Stokes equations on unstructured grids, including shock-capturing via limiters. Similarly, FEniCS facilitates DG solvers for both compressible and incompressible flows, with built-in support for hybridized variants and adaptive mesh refinement in turbulent simulations.

Wave propagation and acoustics

The discontinuous Galerkin (DG) method has been extensively applied to the linear \partial_{tt} u - c^2 \Delta u = 0, where c denotes the wave speed, due to its ability to handle high-order approximations and complex geometries while maintaining for hyperbolic problems. Early formulations employed space-time DG approaches, which discretize both spatial and temporal domains simultaneously using discontinuous bases, enabling exact enforcement of laws and estimates across element interfaces. More recent space-time DG methods extend this to viscoelastic wave propagation, incorporating models while preserving high-order accuracy in both dimensions. For time marching, DG schemes often couple spatial discretization with implicit or explicit time-stepping schemes, such as the Newmark method, which provides second-order accuracy and unconditional for linear systems when combined with appropriate treatments. In acoustics, high-order DG methods are particularly effective for aeroacoustics simulations, where they resolve fine-scale wave features in turbulent flows using unstructured meshes and polynomial degrees up to order 10 or higher. A key challenge in frequency-domain acoustics is the pollution effect in the Helmholtz equation, which causes phase errors to accumulate with increasing wave number k, degrading solution accuracy beyond the pre-asymptotic regime. To mitigate this, hp-adaptive DG strategies dynamically adjust polynomial order p and mesh size h based on a posteriori error indicators, achieving quasi-optimal convergence even for high-frequency problems by refining elements in regions of rapid phase variation. These adaptations ensure that the relative error remains bounded independently of k, as demonstrated in two-dimensional Helmholtz benchmarks with point sources. Elliptic formulations for the frequency-domain Helmholtz equation can be referenced in such contexts to handle time-harmonic acoustics. For electromagnetics, DG methods excel in time-domain simulations of , formulating the system as a hyperbolic set and using upwind numerical fluxes to capture wave discontinuities and ensure hyperbolicity. This approach supports explicit time integration on unstructured grids, making it suitable for transient electromagnetic wave propagation in heterogeneous media. To simulate open-domain problems, perfectly matched layers (PMLs) are integrated via stretched-coordinate formulations within the DG framework, absorbing outgoing waves with minimal by modifying the metric in the PML region while preserving the upwind flux structure. Numerical validations show reflection coefficients below 10^{-4} for plane waves at grazing incidence, confirming the method's efficacy for simulations. DG methods also address nonlinear wave equations, such as the i \partial_t u + \Delta u + |u|^2 u = 0, through local DG (LDG) formulations that rewrite higher-order derivatives as mixed systems and employ alternating fluxes for . These schemes exhibit dispersion-preserving properties, accurately capturing propagation and maintaining long-time numerical phase errors below 10^{-6} over thousands of periods due to the high-order polynomial basis and conservative interface treatments. Similarly, for the Korteweg-de Vries (KdV) equation \partial_t u + \partial_x^3 u + u \partial_x u = 0, direct DG methods with optimized numerical fluxes preserve the dispersive nature of solitary waves, achieving optimal L^2-error convergence of order k+1 for polynomial k while conserving mass and energy in the semidiscrete sense. Recent advances in the 2020s have incorporated DG frameworks for (UQ) in wave propagation, treating parameters like wave speed or medium properties as random fields via stochastic Galerkin projections. These methods expand solutions in a tensor-product basis of spatial polynomials and orthogonal polynomials, enabling efficient propagation of input uncertainties to output quantities like wave amplitudes, with variance reductions up to two orders of magnitude compared to sampling. For instance, multi-order DG hybrids have been applied to in random , providing and variance estimates with computational costs scaling as O(\epsilon^{-2}) for tolerance \epsilon, outperforming traditional stochastic in high-dimensional spaces. Such approaches are particularly valuable for seismic or wave modeling under material variability.

Numerical Implementation

Choice of basis functions

In the discontinuous Galerkin (DG) method, the choice of basis functions within each element is crucial for achieving a balance between numerical accuracy, computational efficiency, and stability. The approximation space on each element is typically spanned by polynomials of degree at most p, denoted as \mathcal{P}^p(K) for element K, allowing for discontinuous functions across element interfaces. Polynomial bases are broadly categorized into modal and nodal types. Modal bases represent the solution as a linear combination of global modes, such as orthogonal polynomials, which promote sparsity in the mass matrix due to their orthogonality properties. In contrast, nodal bases use Lagrange interpolation polynomials defined at specific points (nodes) within the element, facilitating straightforward interpolation of pointwise values but often resulting in denser matrices. Modal approaches are preferred for their efficiency in spectral-like methods, while nodal bases simplify handling of boundary conditions and data interpolation. Orthogonal polynomial bases, particularly , are widely adopted in DG methods, especially on reference intervals or tensor-product elements like . The Legendre basis functions \phi_j(\xi) = \sqrt{2j+1} P_j(\xi), where P_j are the standard Legendre polynomials on [-1,1], ensure L^2-orthogonality, leading to a M = \diag(w_i) after appropriate scaling. This enables fast mass lumping, simplifying the inversion required in explicit time-stepping schemes and reducing computational cost. For elements, tensor-product bases formed by products of one-dimensional Legendre polynomials extend this efficiency, allowing separable computations via sum factorization. On simplicial elements such as , non-standard bases like the Dubiner basis address challenges in high-order approximations. The Dubiner basis consists of hierarchical, orthogonal polynomials constructed by warping one-dimensional onto the reference , defined as \phi_{kl}(\xi,\eta) = (1-\eta)^k P_k^{(0,0)}\left(\frac{2\xi}{1-\eta} - 1\right) P_l^{(2k+1,0)}(2\eta - 1), where \xi + \eta \leq 1. This basis maintains orthogonality in the L^2 norm on the , avoiding the severe ill-conditioning that plagues bases (e.g., \xi^i \eta^j) for p > 7, where basis functions become nearly linearly dependent. By preserving conditioning up to higher degrees, the Dubiner basis supports efficient high-order DG simulations on unstructured simplicial meshes. In hp-adaptive DG methods, the polynomial degree p varies per element, combined with local mesh refinement (h-adaptivity), to exploit solution regularity. For smooth solutions, this hp-strategy yields exponential convergence rates in terms of the number of degrees of freedom, far superior to uniform p-refinement, as the error decreases as O(e^{-c N^{1/d}}) where N is the total degrees of freedom and d the dimension. Adaptivity algorithms use a posteriori error estimators to dynamically adjust p and h, optimizing resource allocation near singularities while maintaining high order elsewhere. Key considerations in basis selection include the conditioning of the mass matrix and the choice of quadrature rules. Orthogonal bases like Legendre minimize issues compared to nodal Lagrange bases on Gauss-Lobatto points, which can exhibit growth as O(p^2) or worse in high dimensions. For exact integration of bilinear forms in DG (which involve products up to degree $2p), Gauss quadrature with (p+1) points per direction suffices on intervals, ensuring no integration error for the volume terms while maintaining efficiency.

Time discretization and integration

For time-dependent problems, the discontinuous Galerkin (DG) method often employs explicit time-stepping schemes to evolve the semi-discrete in space. A prominent approach is the Runge-Kutta DG (RKDG) method, which combines a bounded (TVB) with a third-order strong stability-preserving Runge-Kutta scheme consisting of three forward Euler stages. This explicit scheme was initially developed for scalar hyperbolic conservation laws and extended to systems, achieving high-order accuracy while maintaining stability under a Courant-Friedrichs-Lewy (CFL) condition of order O(1/p^2), where p is the polynomial degree. Implicit methods are preferred for problems involving stiff terms, such as diffusion-dominated flows, where explicit schemes suffer from severe time-step restrictions. Backward Euler or (BDF) schemes are commonly used for temporal discretization, treating diffusive fluxes implicitly to allow larger time steps. For nonlinear systems arising from these implicit formulations, Newton-Krylov solvers are applied to resolve the coupled algebraic equations efficiently. Space-time DG methods extend the discontinuity to both spatial and temporal variables, enabling fully discontinuous approximations for enhanced flexibility. These formulations facilitate adaptive time-stepping and handle non-local time dependencies, such as in wave equations with memory effects, by formulating the problem over space-time slabs. High-order time integration in DG often relies on Gauss-Legendre points within each time interval to achieve accuracy. To update the solution, the resulting from the DG weak form must be inverted; techniques like mass lumping—via rules such as Gauss-Lobatto—or exact diagonalization of block-diagonal matrices are employed to facilitate efficient explicit or semi-implicit updates. Key challenges in these time discretizations include enforcing the Courant condition in explicit RKDG schemes, which limits the time step to ensure proportional to the square of the inverse . For implicit methods, effective preconditioning is essential to accelerate iterative solvers; in hybridizable DG (HDG) variants, multigrid preconditioners operating on variables reduce condition numbers and enable scalable performance for large-scale problems.