Robust Discontinuous Galerkin Methods Maintaining Physical Constraints for General Relativistic Hydrodynamics

Numerically simulating general relativistic hydrodynamics (GRHD) presents significant challenges, including handling curved spacetime, achieving non-oscillatory shock-capturing and high-order accuracy, and maintaining essential physical constraints (such as positive density and pressure, and subluminal fluid velocity) under strong nonlinear coupling. This paper develops high-order accurate, physical-constraint-preserving, oscillation-eliminating discontinuous Galerkin (PCP-OEDG) schemes with the Harten–Lax–van Leer flux for GRHD. To suppress spurious oscillations near discontinuities, we incorporate an oscillation-eliminating (OE) procedure after each Runge–Kutta stage. This OE procedure, based on the exact solution operator of a novel linear damping equation, is computationally efficient and avoids the need for complicated characteristic decomposition. It ensures effective oscillation suppression while preserving the high-order accuracy and conservation properties of the DG method. To further enhance the stability and robustness of the DG method, we develop fully physical-constraint-preserving (PCP) schemes. First, we utilize the W-form of GRHD equations, which reformulates the (3+1) Arnowitt–Deser–Misner formalism via the Cholesky decomposition of the spatial metric. This addresses the challenge of the non-equivalence of admissible state sets at different points in curved spacetime, enabling the construction of provably PCP schemes via convexity techniques. Second, we rigorously prove the PCP property of cell averages using highly technical estimates and the Geometric Quasi-Linearization (GQL) approach (Wu and Shu, 2023) [57], which equivalently casts complex nonlinear constraints into linear ones by introducing auxiliary variables. Our proof shows that, with the enforcement of a simple PCP limiter, the updated cell averages of the OEDG solutions remain physically admissible throughout the simulation. Finally, we present provably convergent PCP iterative algorithms for the robust recovery of primitive variables, ensuring that these variables, approximately recovered from the evolved variables, satisfy the physical constraints throughout the iterative process. The resulting PCP-OEDG method is validated through extensive numerical experiments, including classical test problems in flat spacetime, axisymmetric ultra-relativistic jet flows, and accretion onto rotating black holes in the Kerr metric. These results demonstrate our method's robustness, accuracy, and ability to handle extreme GRHD scenarios involving strong shocks, high Lorentz factors, and strong gravitational fields.