PaperlayerView exact revision on arXiv

Source-linked AI summary

Entropy-Stable and Physical-Constraint-Preserving DGSEM for Symmetry-Reduced General-Relativistic Hydrodynamics on Stationary Spacetimes

Guosheng Fu, Jian-Guo Liu

arXiv:2608.29229v1math.NAgr-qcphysics.comp-ph

TL;DR

The paper addresses stable and physically admissible high-order DG discretization of GRHD on curved stationary spacetimes. It separates SRHD fluid algebra from geometry, couples compatible entropy and PCP mechanisms, and reports high-order accuracy and robust behavior across smooth, shocked, jet, and accretion tests.

  • Problem

    High-order GRHD discretizations must control nonlinear stability and preserve positive density, pressure, and subluminal velocity while handling geometry-dependent fluxes and sources.

  • Method

    The method uses a local orthonormal transformation, compatible entropy-conservative flux and geometric-source discretization, metric-independent PCP limiting, and a common geometry-only causal speed.

  • Results

    The fully discrete DGSEM recovers high-order accuracy for smooth SRHD and Michel flows and robustly resolves multidimensional shocks, an axisymmetric jet, and Schwarzschild and Kerr accretion flows.

  • Takeaways & Limitations

    The results support the framework as a practical high-order method for relativistic flows on curved stationary backgrounds.

  • Takeaways & Limitations

    The analysis is restricted to prescribed stationary spacetimes, conforming affine tensor-product elements, two active coordinates, and periodic semidiscrete entropy analysis.

Abstract

from arXiv · show

We develop an entropy-stable and physical-constraint-preserving discontinuous Galerkin spectral element method for symmetry-reduced general-relativistic hydrodynamics on prescribed stationary spacetimes. Using a local orthonormal transformation, the fluid variables are expressed in a form for which the relativistic hydrodynamic algebra and the admissible set are independent of the spatial metric, while the spacetime geometry enters through stationary coefficients. This separation allows entropy-conservative special-relativistic fluxes to be combined with a compatible discretization of the geometric source terms. On affine tensor-product meshes, the resulting DGSEM is conservative and satisfies a semidiscrete entropy inequality, while the transformed variables provide a convex framework for physical-constraint preservation. For practical stabilization, we use a geometry-only causal speed that is sufficient for both classical local Lax--Friedrichs entropy dissipation and the physical-constraint-preserving Lax--Friedrichs splitting. The fully discrete method combines this stabilization with SSP Runge--Kutta time stepping, oscillation elimination, and conservative local-orthonormal-state scaling. Numerical experiments cover smooth and strongly shocked special-relativistic flows, an axisymmetric jet, stationary Michel accretion, Schwarzschild Bondi--Hoyle flow, and four Kerr accretion cases. The results demonstrate the designed high-order accuracy in smooth regimes and robust performance for demanding relativistic flows on curved stationary backgrounds.

1. Introduction

The paper develops an entropy-stable and physical-constraint-preserving DGSEM for symmetry-reduced GRHD by separating metric-independent fluid algebra from stationary geometry. Compatible flux–source discretization, transformed-state limiting, and common causal-speed stabilization support conservation, entropy control, admissibility, and robust numerical tests.

  • Motivation: The method targets nonlinear stability, physical admissibility, and spatially varying geometry in high-order DG discretizations of GRHD.GRHD models relativistic fluids on curved spacetimes, while DG methods provide high-order approximation, conservation, and parallel scalability.
  • Framework: A local orthonormal transformation gives the fluid variables standard SRHD algebra while stationary spacetime geometry enters through coefficients and the reduced conservative measure.The suppressed spatial direction contributes a metric factor through the weighted reduced state.
  • Entropy stability: Matching the entropy-conservative flux and geometric-source discretizations with SBP operators and nodal geometry produces the required discrete flux–source contraction.On conforming affine tensor-product elements, the DGSEM is conservative and satisfies a semidiscrete entropy inequality with an entropy-dissipative interface flux.
  • Physical-constraint preservation: The transformed variables define a metric-independent admissible cone for a conservative weighted PCP limiter that preserves cell averages while enforcing density, cone, and relative-margin constraints.The reduction weight may vary spatially but is assumed positive at nodes.
  • Stabilization: A geometry-only causal speed supplies both the Lax–Friedrichs splitting required by PCP analysis and sufficient dissipation for the discrete entropy condition.This provides a common stabilization scale for entropy stability and physical-constraint preservation.
  • Validation and scope: The fully discrete scheme combines SSPRK(3, 3), oscillation elimination, and conservative local-state scaling, and performs accurately and robustly across smooth, shocked, jet, and black-hole accretion tests.The analysis is restricted to prescribed stationary spacetimes, periodic semidiscrete entropy analysis, and conforming affine tensor-product elements in two active coordinates.

2. GRHD model and symmetry reduction

The model rewrites GRHD on stationary spacetimes in local orthonormal variables, separating SRHD fluid algebra from geometry while retaining the suppressed-direction measure and source. The resulting two-dimensional framework covers planar, axisymmetric, Schwarzschild, and equatorial Kerr reductions.

  • The Valencia GRHD balance law evolves relativistic fluids on prescribed stationary curved spacetimes with geometry-dependent fluxes and sources.
  • From covariant GRHD to the W-form: The local orthonormal transformation makes the fluid flux take standard SRHD form while geometry enters through stationary coefficients.
  • Two-dimensional W-form reductions: Suppressing one spatial direction retains its metric scale in the conservative measure and generates an additional weighted geometric source.
  • Model geometries: The framework treats planar, cylindrical, and Schwarzschild cases as true symmetry reductions, while Kerr is an infinitesimally thin equatorial restriction.
  • Model geometries: Planar SRHD has unit reduction weight and vanishing sources, whereas axisymmetric cylindrical SRHD uses w = r and a nonzero suppressed-direction source.
  • Scope and vanishing metric scales: The reduced four-component model assumes zero suppressed velocity; nonzero suppressed velocity requires an enlarged state outside the present scope.

3. Entropy-stable DGSEM on affine tensor-product meshes

The DGSEM pairs entropy-conservative SRHD volume fluxes with a geometry-compatible source discretization and LLF surface fluxes on affine tensor-product elements. This construction aligns the discrete source contraction with the entropy analysis.

  • The construction is formulated on conforming affine tensor-product elements with positive reduction weight at solution nodes.
  • Entropy-conservative volume flux and compatible source: The entropy variables and potentials are evaluated from local orthonormal states, enabling SRHD entropy-conservative two-point fluxes in the reduced system.
  • Entropy-conservative volume flux and compatible source: The numerical source uses complete weighted geometry coefficients so its discrete entropy contraction matches the continuous source identity.
  • LLF surface flux and semidiscrete scheme: The surface discretization uses the classical local Lax–Friedrichs flux with a geometry-only causal speed.
  • LLF surface flux and semidiscrete scheme: Combining the LLF interface flux, entropy-conservative volume flux, and compatible source defines the single-element semidiscretization.

The scheme (94) satisfies the conservative balance

For admissible nodal states on conforming affine tensor-product meshes, the scheme combines causal LLF stabilization with compatible sources to obtain semidiscrete entropy stability and weighted rest-mass conservation under periodic boundaries.

  • The causal LLF speed satisfies the entropy-stability condition for admissible states sharing the same prescribed face geometry.
  • The same causal bound also provides the Lax–Friedrichs splitting required by the physical-constraint-preserving analysis.
  • Under periodic boundaries, the affine-mesh DGSEM satisfies a semidiscrete entropy inequality when nodal states are admissible and geometry is single-valued at shared faces.
  • The scheme conserves global weighted rest mass through cancellation of shared numerical fluxes.

4. Physical-constraint-preserving analysis

The PCP analysis uses a metric-independent admissible cone and proves sufficient forward-Euler preservation conditions for cell averages. A conservative three-stage scaling limiter enforces density, cone, and relative-margin constraints while preserving conservative averages.

  • Metric-independent admissible cone: The admissible set is an open convex cone, and its closure is also convex.Positive scaling and convexity support the cone-based preservation arguments.
  • Metric-independent admissible cone: The transformed local, intrinsic, and evolved states share a metric-independent admissible cone.For the Gamma-law gas, cone membership corresponds to positive density and pressure with subluminal velocity.
  • Forward-Euler PCP condition: The PCP Lax–Friedrichs splitting follows from the coordinate-direction causal bound and SRHD splitting properties.The proof uses rotational invariance, subluminal characteristic speeds, and cone scaling.
  • Forward-Euler PCP condition: A sufficient forward-Euler condition preserves admissible cell averages through a convex decomposition of volume, face, and source contributions.The condition requires admissible nodal states, positive reduction weights, compatible face geometry, and a positive residual coefficient.
  • Conservative weight-compatible scaling: The analysis applies positive numerical floors and timestep halving without altering admissible conservative cell averages.Theorem 4.3 guarantees a positive admissible timestep threshold, while scaling can collapse a polynomial to its admissible cell-average anchor when needed.
  • Conservative weight-compatible scaling: The three-stage limiter preserves the conservative cell average while imposing density, cone, and relative cone-margin constraints.The final stage retains at least a fixed fraction of the anchor state's normalized cone margin, improving robustness near the admissible-set boundary.
  • Conservative weight-compatible scaling: The conservative scaling framework extends to axis cells where the reduction weight has a simple normal zero.This includes the polar-axis cells used in the numerical experiments.

5. Fully discrete scheme

The fully discrete scheme combines geometry-aware residual evaluation with entropy-stable LLF interface stabilization, directional oscillation elimination, and conservative PCP scaling. Cross-line normalization avoids domain-wide amplitude contamination, while the complete stage map preserves conservative cell averages.

  • Scheme components: The implementation assembles primitive recovery, LLF fluxes, stage-wise oscillation elimination, PCP scaling, SSPRK(3, 3), and geometry-based timestep control.These components complete the fully discrete DGSEM used in the numerical experiments.
  • Residual evaluation and interface stabilization: The volume discretization uses an entropy-conservative flux, while interfaces use a classical LLF flux with a geometry-only dissipation coefficient.The same coefficient supplies dissipation for both the entropy inequality and PCP Lax–Friedrichs splitting.
  • Conservative local-state oscillation elimination: Oscillation elimination acts on the regular local orthonormal conservative state rather than the densitized state.This prevents the stationary geometric factor from being interpreted as fluid variation.
  • Directional oscillation elimination: The directional OE indicator uses a geometry-based causal rate and an effective DG resolution length for order-dependent damping.Directional jump amplitudes are normalized separately by coordinate direction.
  • Cross-line normalization: Domain-wide normalization can make jumps artificially small when near-hole conservative excursions dominate the global scale.This is especially problematic for black-hole accretion with large variation between near-horizon and outer-flow regions.
  • Cross-line normalization: Cross-line scaling assigns direction-matched normalization along structured coordinate lines, localizing the indicator without restricting it to a single element.The construction is used for all reported OE calculations and is invariant under componentwise affine rescaling.
  • Conservative OE map: The OE map retains hierarchical increments with damping controlled by sOE, using 0.01 ≤sOE ≤0.05 for the Q2 experiments.The study supports sOE = 0.02 as a reasonable default for these calculations.
  • Conservative OE map: The damped OE increments have zero weighted mean, so conservative OE preserves the stored conservative cell average exactly.The subsequent PCP scaling acts on the same regular local orthonormal state and preserves the conservative element average.

6. Numerical experiments

The experiments test the fully discrete DGSEM from smooth accuracy problems through shocks, axisymmetric jets, stationary Michel flow, and black-hole accretion. Results retain high-order accuracy in smooth regimes while preserving admissibility and resolving demanding relativistic structures.

  • Experiment design: The test suite spans smooth SRHD accuracy, strongly nonsmooth relativistic flows, axisymmetric discretizations, stationary Michel flow, and Schwarzschild and Kerr accretion.The configurations use Q2 approximations, conforming affine meshes, PCP, classical LLF fluxes, geometry-based dissipation, and CFL = 0.8.
  • 6.1. Two-dimensional smooth-wave accuracy: Third-order convergence is retained for every tested OE strength in the smooth traveling-wave problem.Finest-mesh observed orders range from 3.04 to 3.34 in L1, 3.04 to 3.35 in L2, and 3.04 to 3.31 in L∞.
  • 6.2. Two-dimensional Riemann problems: The three Riemann problems preserve admissibility while resolving vortex-sheet roll-up, curved shocks, and interacting contacts, including ultra-relativistic states.All cases use Q2 elements on 400 × 400 meshes and evolve to t = 0.4 with sOE = 0.02.
  • 6.3. Shock–bubble interaction: Shock–bubble computations resolve transmitted and reflected waves, compressed interfaces, vortical structures, complex wakes, and smaller-scale waves for both light and heavy bubbles.Across OE strengths, the principal wake geometry and paired vortical structures remain consistent, while stronger damping broadens interfaces and smooths spiral tips.
  • 6.4. Axisymmetric relativistic jet: The axisymmetric jet develops a bow shock, expanding cocoon, and shear-driven vortical structures that remain sharply resolved through t = 100.The jethead propagation, cocoon shape, and internal morphology agree well with prior simulations and subsequent high-order calculations.
  • 6.5. Steady Michel accretion and 6.7. Equatorial Kerr–Schild Bondi–Hoyle accretion: OE improves stationary Michel accuracy while retaining high-order convergence, and Kerr accretion produces well-defined tail shocks whose large-scale cone changes little with spin.For Michel flow, all tested OE strengths reduce errors relative to the undamped scheme; for Kerr cases, rotational influence is concentrated near the black hole.

7. Conclusions

The work develops a compatible entropy-stable and physical-constraint-preserving DGSEM for symmetry-reduced GRHD, with numerical tests supporting high-order accuracy and robust performance on stationary curved backgrounds. Its present scope is limited to prescribed stationary geometries and specified mesh and dimensional settings.

  • The local orthonormal transformation separates standard SRHD fluid algebra from stationary geometry while retaining the physical reduction weight in the conservative measure.
  • Using shared nodal geometry for two-point fluxes and geometric sources yields the required flux–source cancellation and a semidiscrete entropy result.
  • SSPRK(3, 3) time stepping with OE and PCP scaling supports high-order accuracy for smooth flows and robust resolution of shocks, jets, and Schwarzschild and Kerr accretion.
  • Curvilinear meshes, nonconforming interfaces, evolving spacetimes, and three-dimensional discretizations remain identified extensions.

Appendix A. Explicit source terms for the black-hole geometries

The appendix specifies source construction for Schwarzschild and equatorial Kerr systems by combining active Valencia sources with suppressed-direction contributions and local-frame derivatives. The resulting substitutions completely determine the geometric sources for both models.

  • Schwarzschild source terms: The Schwarzschild source combines the intrinsic shared source form with the explicitly supplied suppressed-direction contribution S⊥.
  • Schwarzschild source terms: For Schwarzschild, metric derivatives determine the active energy source, while the active θ-momentum source vanishes because the active spacetime block is independent of θ.
  • Schwarzschild source terms: The Schwarzschild source is completed by evaluating the orthonormal factor and substituting the listed metric and frame terms into the shared source expression.
  • Equatorial Kerr source terms: The equatorial Kerr source likewise combines the intrinsic active Valencia source with the explicitly given suppressed-direction contribution S⊥.
  • Equatorial Kerr source terms: For Kerr, the azimuthal momentum source vanishes because the metric is independent of eϕ, and the active energy source follows from the metric derivatives.
  • Equatorial Kerr source terms: The Kerr source is completed by evaluating the radial derivative of the local frame and substituting the listed terms into the shared source expression.

Appendix B.2. Directional spatial EC flux

The directional spatial entropy-conservative flux is assembled from symmetric state averages and local directional quantities. Its construction is dimension-independent under orthonormal-index contractions but requires admissible endpoint primitives and stable logarithmic-mean evaluation.

  • The directional construction uses averaged velocity, Lorentz-factor-related quantities, and local directional data to assemble the flux components.
  • The assembled vector is symmetric and consistent and satisfies the stated entropy-conservative relation.
  • Only contractions over the orthonormal index are required, so the formula has the same form in two and three dimensions.
  • The formulas require admissible endpoint primitives and stable evaluation of every logarithmic mean.

Appendix C. Flux-adapted geometry variables

The appendix introduces flux-adapted geometry variables and their differentials for the active ADM geometry and separate physical reduction weight. The resulting coordinate map is smooth and one-to-one when j2 > 0.

  • The active ADM geometry is paired with a separate prescribed scalar physical reduction weight w.
  • The appendix defines a flux-adapted geometry state and forms complete flux coefficients in the forward direction.
  • The map between the original and flux-adapted active geometry coordinates is smooth and one-to-one for j2 > 0.
  • The physical reduction weight modifies the geometry coefficients and their differentials through product-rule terms involving dIw.
  • Directional increments of the intrinsic geometry state yield algebraic differentials of the local geometry map.

Appendix D. Causal bound for classical-LLF entropy stability

The appendix establishes that the geometry-only causal speed yields classical-LLF entropy stability by combining the convex entropy-variable domain with an entropy-defect bound.

  • Causal-speed bound: The geometry-only causal speed bounds the local orthonormal SRHD characteristic speeds and supplies the estimate required along every entropy-variable path.The characteristic analysis uses rotational invariance of SRHD and the Gamma-law gas assumption.
  • Entropy-defect estimate: Strict entropy convexity makes the relative-entropy defect positive for distinct states, reducing the LLF proof to bounding the associated entropy production.The appendix formulates the classical LLF entropy defect using the entropy variables and entropy potential.
  • Convex entropy-variable domain: For 1 < Γ ≤ 2, admissible entropy variables form the convex set V_D ∈ R and −V_E > |V_m|.This convexity keeps straight entropy-variable paths between admissible states inside the admissible domain.
  • Entropy-variable path: The entropy-variable path connects two admissible states through V(θ) = V_L + θ[V], with U(θ) recovered from the inverse entropy-variable map.The path remains admissible because the entropy-variable domain is convex.
  • Conclusion: Combining the pathwise entropy-defect lemma with the characteristic-speed estimate proves the desired classical-LLF entropy inequality.The result establishes Proposition 3.3 for the fixed positive-weight face geometry.

Appendix E. Axis cells with vanishing reduction weight

The appendix gives a conservative treatment of axis cells where the reduction weight vanishes, preserving regular intrinsic states, stored axis values, and PCP scaling constraints.

  • Axis regularity: At symmetry axes, the reduction weight vanishes while the weighted state remains regular and the underlying local orthonormal state can have a finite nonzero limit.The analysis otherwise assumes positive reduction weight at solution nodes.
  • Regular axis trace: The stored densitized state is kept at the exact axis value W⋆ = 0, while a derivative-ratio reconstruction supplies the regular local state when needed.The reconstruction uses the normal SBP derivative along the same tensor-product nodal line and is auxiliary rather than a replacement for W⋆.
  • Flux treatment: Because the weighted flux coefficients vanish on the axis, the weighted axis flux is zero and no axis Riemann problem is formed.The reconstructed regular state is used in adjacent volume and source evaluations and during stage stabilization.
  • Conservative axis repair: A conservative axis repair restores W⋆ = 0 after explicit stages by redistributing the axis defect among positive-weight nodes on the same normal line.The repair preserves the linewise GLL quadrature moment and therefore the conservative element average.
  • PCP scaling: Trace-aware PCP scaling contracts positive-weight nodal states and regular axis traces while preserving the stored axis constraint and conservative average.The weighted Gram matrix is positive definite below the highest mode, whereas the axis-face nodal mode lies in the weighted null space and motivates the special convention.
  • Theorem: Under the stated simple-zero and admissibility assumptions, the three-stage trace-aware scaling preserves W_K and places positive-weight states and regular axis traces in the certification set.If certification conditions already hold, all three scaling factors equal one.
Loading 2608.29229v1…