Discretizing Relativistic PDEs

Relativistic PDEs mix wave propagation, constraints, conservation laws, shocks, geometric source terms, and artificial boundaries. A discretization must respect the structure relevant to its problem: an accurate interior stencil is not enough if its boundary closure is unstable, and a conservative shock method is not enough if primitive states become unphysical. This chapter develops the pieces as one coupled numerical design.

Helpful background. Formulations of General Relativity supplies continuum hyperbolicity, gauges, and constraints; Ordinary Differential Equations supplies time integrators and semidiscrete stability; and Partial Differential Equations and Green Functions supplies weak solutions, characteristics, energy estimates, and boundary data.

From continuum structure to a discrete claim

Start by naming the continuum property the algorithm must preserve or approximate: derivative accuracy, a discrete energy estimate, cell-integrated conservation, weak consistency, constraint control, or outgoing boundary behavior. Then choose spatial and temporal operators that realize it together. Finally design a benchmark that can fail that property independently of visual solution quality.

To make the calculation reproducible, specify the exact PDE and domain, initial and boundary data, formulation and variables, mesh and mapping, spatial operator, flux or penalty, dissipation, time integrator, precision, Courant prescription, refinement sequence, field and constraint norms, conservation or energy history, and a rerun command or complete executable sequence. The stability framework for initial-boundary difference methods is developed in Gustafsson, Kreiss, and Oliger 2013, Chapters 4–5 and 8–13.

The six-page route

1. Finite-Difference Operators moves from Taylor consistency to summation by parts, penalty boundaries, and grid-scale noise.
2. Finite-Volume Conservation Laws derives flux cancellation and tests reconstruction, Riemann solvers, and physical admissibility on a relativistic shock tube.
3. Finite-Element and Discontinuous Galerkin Methods compares continuous and discontinuous weak forms, curved geometry, interface fluxes, and $hp$ refinement; the nodal DG construction is developed in Hesthaven and Warburton 2008, Chapters 2–7, pp. 19–269.
4. Method of Lines and Time Integration matches the semidiscrete spectrum to explicit, implicit, or IMEX stability regions and separates spatial from temporal order.
5. Stability, Constraints, and CFL Limits distinguishes bounded field evolution from preservation of an auxiliary constraint and scans the Courant threshold with random noise.
6. Numerical Boundary Conditions prescribes incoming characteristic sectors, measures reflection, and tests constraint-preserving radiative closures.

The route is cumulative. Spatial operators create a semidiscrete spectrum; time integration turns that spectrum into an amplification problem; constraints add another monitored sector; and boundaries complete the evolution operator. Testing any one in isolation is necessary but not sufficient.

A compact benchmark matrix

Periodic sine advection isolates smooth interior phase and amplitude error. A finite-interval pulse tests energy and reflections. The relativistic shock tube tests conservation, discontinuities, positivity, and primitive recovery. A curved-mesh free stream tests metric identities. A wave system with $\mathcal C=\Phi-\partial_xu$ tests whether constraints remain controlled. Each should be run at three or more resolutions with a fixed or explicitly varied Courant prescription.

Finite-volume wave propagation, nonlinear Riemann problems, and high-resolution reconstruction are developed in LeVeque 2002, Chapters 4–6, 10–12, and 17. Smooth and discontinuous benchmarks must be reported separately because the designed order in smooth regions need not survive a shock in a global norm.

Negative tests make the matrix discriminating. Centered space with Forward Euler amplifies modes with nonzero derivative eigenvalue and is unstable at fixed Courant number. Its Nyquist mode is a separate blind spot. A wrong SBP penalty sign can grow boundary energy. Unlimited reconstruction can oscillate near a discontinuity. A curved mesh tests the geometric identities required by the chosen DG form. A Courant factor above threshold should produce resolution-amplified noise even when a discrete auxiliary constraint stays exactly preserved. The boundary page therefore tests both zero incoming data and a prescribed nonzero reflection coefficient; a permanently zero, decoupled incoming field alone does not exercise the return path.

How to interpret convergence

Hold the Courant factor fixed when refining a fully discrete explicit scheme, then perform separate timestep refinement at fixed grid. Report both smooth-region and global norms for shocks. Track budgets and constraints alongside field error; a field can converge while a conservation defect or auxiliary constraint grows.

The asymptotic regime may begin only after a wave is resolved, but a stability test should not be shortened to avoid visible growth. Predefine final time and acceptance thresholds. Vary dissipation rather than treating it as an invisible constant, and show that it damps unresolved modes without setting the measured physical answer.

Chapter review

A scheme shows fourth-order field convergence for a smooth wave, but its conserved integral drifts at the same nonzero rate on every grid. Which parts of the method should be inspected, and what evidence is missing?

Solution

Inspect whether the PDE was discretized in flux form, whether neighboring cells or elements share exactly opposing interface fluxes, how physical boundaries contribute, and whether sources or refinement interfaces are balanced. The field norm alone does not test conservation. Add a flux budget, boundary-flux integral, and refinement study of the conservation residual; for adaptive meshes, verify refluxing or conservative projection.

References

• Gustafsson, B., Kreiss, H.-O., and Oliger, J. (2013), Time-Dependent Problems and Difference Methods, 2nd edition, Wiley. DOI.
• Hesthaven, J. S., and Warburton, T. (2008), Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, Springer. DOI.
• LeVeque, R. J. (2002), Finite Volume Methods for Hyperbolic Problems, Cambridge University Press. DOI.