Finite-Element and Discontinuous Galerkin Methods

Finite-element methods begin from a weak form and approximate the solution in elementwise polynomial spaces. Continuous Galerkin methods tie neighboring traces together; discontinuous Galerkin methods keep separate traces and couple them with a numerical flux. The comparison is not “smooth versus discontinuous solutions.” It concerns continuity, conservation, communication, stabilization, geometric flexibility, and how error is controlled on a given mesh.

Required background. Ordinary Differential Equations supplies the semidiscrete systems and time integration used after spatial approximation.

Helpful background. Formulations of General Relativity supplies the continuum evolution systems, while Partial Differential Equations and Green Functions provides weak formulations, characteristics, and boundary terms.

Weak form and interface fluxes

For a conservation law $\partial_tu+\nabla\!\cdot F(u)=S$, multiply by a test function $v$ and integrate over an element $K$:

$$\int_Kv\,\partial_tu\,dV -\int_K\nabla v\!\cdot F(u)\,dV +\int_{\partial K}v\,F^*\!\cdot n\,dA =\int_KvS\,dV.$$

In a discontinuous space the interior and exterior traces differ, so $F^*$ selects a consistent, stable coupling. The same interface flux, evaluated once with opposite normals, makes neighboring contributions cancel and yields elementwise conservation. Teukolsky derives the strong and weak DG forms and their flux relation in Teukolsky 2016, § 2.2, Eqs. 2.8–2.13.

To see what is actually evolved, expand $u_h|_K=\sum_jU_j(t)\phi_j$ and choose $v=\phi_i$. The time term becomes $\sum_jM_{ij}\dot U_j$, with $M_{ij}=\int_K\phi_i\phi_j\,dV$. Volume and face quadrature build the right-hand side. Each element solves its own mass-matrix system; communication enters through the face flux. Taking the test function $v=1$ removes the volume derivative and recovers the finite-volume balance for the element average.

Continuous finite elements identify shared degrees of freedom and usually have fewer unknowns for the same polynomial degree and mesh. They fit elliptic problems naturally and yield globally coupled matrices. DG has more degrees of freedom and interface work, but a block-local volume operator, explicit conservation, natural upwinding, and flexible nonconforming $h$ or $p$ interfaces.

Comparison axes for relativistic PDEs

For a smooth wave, both continuous and discontinuous high-order methods can converge rapidly. DG’s numerical flux adds controlled dissipation; continuous Galerkin often needs another stabilization for advection-dominated problems. For shocks, neither an unlimited high-order polynomial nor continuity is sufficient: DG requires limiting, filtering, artificial viscosity, or a subcell method, while continuous schemes require their own shock stabilization.

On unstructured or curved meshes, both methods approximate geometry through element maps. DG localizes mesh communication and supports elementwise polynomial order, which suits parallel adaptivity. Continuous spaces reduce interface unknowns but make hanging-node constraints and global assembly part of the implementation. Compare accuracy at equal degrees of freedom, memory, and wall time; comparing only equal element counts favors neither method fairly.

Hesthaven and Warburton construct nodal bases, mappings, fluxes, and time integration in Nodal Discontinuous Galerkin Methods, Chapters 2–7, beginning at pp. 19, 43, 75, 115, 169, and 243. Their algorithms assume matched choices of interpolation nodes, quadrature, mass matrices, and surface operators; mixing those choices can destroy the discrete energy argument.

Curved-wave benchmark

On the computational torus $0\leq\xi,\eta\lt 2\pi$, use the smooth map

$$x=\xi+0.1\sin\xi\sin\eta, \qquad y=\eta+0.1\sin\xi\sin\eta.$$

Its Jacobian determinant is $J=1+0.1\sin(\xi+\eta)\geq0.9$, so the map is nonsingular. The deformation is periodic and the map has degree one, giving a global torus coordinate map as well as local invertibility. Solve physical-space advection

$$\partial_tu+\partial_xu+\partial_yu=0, \qquad u(t,x,y)=\sin(x+y-2t),$$

through one temporal period, $t=\pi$. Use a nodal DG discretization on tensor-product Legendre–Gauss–Lobatto nodes, an upwind face flux, polynomial degrees $p=2,3,4$, and uniform computational meshes of $4^2$, $8^2$, and $16^2$ elements. State the mapped conservative formulation, quadrature, and discrete metric construction. Integrate with DOP853 using relative and absolute tolerances $10^{-12}$ and maximum step $0.05h_{\min}/(2p+1)$, shortening the last step to land at $t=\pi$. Report physical $L^2$ and $L^\infty$ error, phase and amplitude error, conservation residual, energy history, and work.

The physical $(1,1)$ Fourier coefficient is

$$\widehat u_{11}(t)=\frac{1}{(2\pi)^2} \int_0^{2\pi}\!\int_0^{2\pi} u_h(t,x(\xi,\eta),y(\xi,\eta)) e^{-i(x(\xi,\eta)+y(\xi,\eta))}J\,d\xi\,d\eta.$$

Its exact value is $e^{-2it}/(2i)$. Compare this quantity using independent, sufficiently accurate quadrature. The pulled-back solution $\sin[\xi+\eta+0.2\sin\xi\sin\eta-2t]$ is not a single computational Fourier mode; projecting onto $e^{-i(\xi+\eta)}$ would measure a different observable. Physical energy and $L^2$ norms likewise include $J$.

For smooth data, refine $p$ at fixed mesh and $h$ at fixed $p$. The observed convergence should match the implemented quadrature and mapping order; halve the maximum timestep once at fixed mesh to verify that time error is subdominant. Then initialize $u=1$ on the same curved meshes, with the identical map, metrics, flux, and time integrator, and evolve to $t=\pi$. It must remain constant to roundoff: this free-stream test probes discrete metric identities and face geometry rather than wave accuracy.

On each internal face, export the two flux contributions and verify cancellation. Teukolsky shows how mapped formulations and mass-lumping choices affect equivalence of strong and weak forms in Teukolsky 2016, § 4, Eqs. 4.5 and 4.17, and § 6, Eq. 6.3. Free-stream preservation depends on the chosen form: discrete metric identities are needed in the transformed conservative form, while the alternative form analyzed there preserves a constant differently. A free stream that drifts when the mesh is curved calls for inspection of the actual formulation, metric factors, normals, and quadrature even if a Cartesian test passes.

Aliasing, boundaries, and hp refinement

Nonlinear fluxes are not represented exactly by the solution polynomial. Underintegration aliases unresolved high modes into resolved ones and can drive instability. Overintegration, split forms, filtering, or modal limiting can control it, but each changes cost and sometimes conservation. Record quadrature order and every stabilization coefficient.

At physical boundaries, use the characteristic data prescribed by the continuum problem to build the exterior state or penalty. At $hp$ interfaces, ensure both sides share a conservative mortar or projection. Refinement indicators may use residuals or spectral tails, but validate them on a known localized feature before trusting adaptive placement.

Reproducible environment and failure criteria

Pin the implementation revision, mesh and generator, geometry order, basis and nodes, quadrature, mass-matrix treatment, numerical flux, limiter or filter, boundary state, time integrator, precision, and parallel partition. A clean rerun exports the mesh, curved metrics, face fluxes, error and conservation norms, refinement tables, and free-stream residual.

The validity regime assumes a well-posed continuum problem and solution regularity compatible with the claimed order. Failure modes include violation of metric identities, mismatched face normals, aliasing instability, underintegrated nonlinear terms, loss of positivity, and temporal order reduction. Reject convergence seen only on Cartesian meshes or in a field norm while conservation or constraints grow. A limiter that activates on the smooth benchmark also requires explanation because it may set the measured order.

Check your understanding

A DG code converges at the designed order on a Cartesian mesh but evolves a constant state on a curved mesh into a nonconstant one. What should be checked first?

Solution

Check the discrete geometric identities: mapped Jacobians, contravariant metric factors, face normals, quadrature, and the equality of opposing face fluxes. A constant state has no physical gradient, so drift identifies geometric or interface inconsistency rather than inadequate resolution. Fix the free-stream test before interpreting curved-wave convergence.

References

• Hesthaven, J. S., and Warburton, T. (2008), Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, Springer. DOI.
• Teukolsky, S. A. (2016), “Formulation of Discontinuous Galerkin Methods for Relativistic Astrophysics,” Journal of Computational Physics 312, 333–356. DOI. Open PDF.