Bound Preserving Lax-Wendroff Flux Reconstruction Method for Special Relativistic Hydrodynamics

Abstract: Lax-Wendroff flux reconstruction (LWFR) schemes have high order of accuracy in both space and time despite having a single internal time step. Here, we design a Jacobian-free LWFR type scheme to solve the special relativistic hydrodynamics equations on Cartesian grids. We then blend the scheme with a first-order finite volume scheme to control the oscillations near discontinuities. We also use a scaling limiter to preserve the physical admissibility of the solution after ensuring the scheme is admissible in means. A particular focus is given to designing a discontinuity indicator model to detect the local non-smoothness in the solution of the highly non-linear relativistic hydrodynamics equations. Finally, we present numerical results for a wide range of test cases to show the robustness and efficiency of the proposed scheme.

• The paper introduces a Jacobian-free LWFR scheme that maintains high-order accuracy while preserving density positivity and subluminal velocities.
• It employs a blending limiter and a novel discontinuity indicator to effectively control oscillations near shocks and discontinuities.
• The method is validated with various 1D and 2D tests, achieving up to fifth-order convergence and robust performance in extreme relativistic regimes.

Bound-Preserving Lax-Wendroff Flux Reconstruction for Special Relativistic Hydrodynamics

Introduction and Motivation

The "Bound Preserving Lax-Wendroff Flux Reconstruction Method for Special Relativistic Hydrodynamics" (2409.15805) addresses a critical challenge in the numerical integration of special relativistic hydrodynamics (RHD): maintaining high-order accuracy, robustness, and physical admissibility of the solution. In astrophysics and high-energy physics, RHD is ubiquitous due to the relativistic speeds and strong gravitational effects involved in many phenomena such as jets, black holes, and gamma-ray bursts. Existing high-order schemes, including WENO, PPM, DG, and their variants, often fail to guarantee positivity of density and pressure and enforcement of the velocity bound ($|v|<1$), especially near discontinuities or in low-pressure regimes, which leads to nonphysical states and eventual breakdown of the computation.

The presented work introduces a high-order, Jacobian-free Lax-Wendroff Flux Reconstruction (LWFR) framework for the RHD equations. The scheme achieves spatial and temporal accuracy in a single step, blends with a first-order admissibility-preserving finite volume scheme for oscillation and unphysical state control, and employs a novel scaling limiter and robust discontinuity indicator.

Relativistic Hydrodynamics and Admissible State Enforcement

The governing equations are the conservation laws of mass, momentum, and energy for an ideal fluid, with fluxes strongly coupled through the Lorentz factor. Hyperbolicity (guaranteed real eigenvalues) and well-posedness require density $\rho>0$, pressure $p>0$, and $|v|<1$. Directly enforcing these conditions is nontrivial due to the lack of explicit expressions in conservative variables. The authors reformulate the admissibility region using a convex set $U' = \{u: D>0,\, E-\sqrt{D^2+|m|^2}>0\}$ where $D$ is relativistic density and $E$ is total energy, so that admissibility can be efficiently checked and enforced after each update.

Lax-Wendroff Flux Reconstruction Framework

The LWFR method constructs a polynomial representation of the solution in each element using Gauss-Legendre solution points and applies a single-stage high-order temporal update via Taylor series expansion. Spatial derivatives are efficiently computed using approximate finite difference formulae rather than directly evaluating the Jacobian, resulting in computational savings especially relevant for high-order schemes.

The time-averaged fluxes are reconstructed within each element accompanied by correction functions (Radau polynomials), ensuring continuity and compatibility with nodal DG approaches. Numerical fluxes at inter-element boundaries employ a Rusanov-type dissipation mechanism.

Figure 1: Division of the element $\Omega_{pq}$ into $4\times 4$ sub-elements, with internal solution points for flux reconstruction.

Subcell Blending for Oscillation Control

High-order methods are susceptible to Gibbs phenomena near discontinuities, which can lead to nonphysical overshoots and negative values. Traditional TVD/TVB or polynomial limiters either degrade the solution to first-order or fail to reliably distinguish true discontinuities from smooth extrema. The approach taken here employs a blending limiter: the solution is updated as a convex combination of the high-order LWFR update and a robust first-order finite volume evolution, acting on a subcell partition of each element.

The blending coefficient is governed by local solution regularity, ideally retaining high-order accuracy in smooth regions and suppressing oscillations near discontinuities.

Novel Discontinuity Indicator

Detecting non-smoothness in RHD is especially challenging due to the compound nonlinearity of the equations and the influence of the Lorentz factor. The paper introduces a smoothness indicator based on the Legendre modal energy content of the product $K = \rho p \Gamma$, incorporating density, pressure, and Lorentz effects. By extending the support to boundary points (not only interior nodes) and exploiting modal decay, the indicator robustly distinguishes discontinuities—mitigating both “false positives” at smooth extrema and failures to flag sharp but small-amplitude features.

Figure 2: Density field in the Kelvin–Helmholtz instability at $t=3$, using the new discontinuity indicator, showing fine-scale structures and oscillation suppression.

Bound Preservation and Scaling Limiter

To guarantee positivity of density, pressure, and subluminal velocities after each update, the scheme enforces admissibility in both means (element-averaged state) and at each solution point. The authors adapt the scaling limiter of Zhang and Shu, but extend it for RHD by operating in the convex “admissible” region in conservative variables—a property ensured by the aforementioned reformulation. The blended and flux-corrected scheme enforces bounds at the subcell level, with further convex limiting as required.

Figure 3: $\ln \rho$ (log-density) field for the two-dimensional Riemann problem at $t=0.4$ using the fourth-order scheme, with 25 contour lines in $[-4.30,2.09]$, highlighting robust shock-capturing without spurious oscillations.

Numerical Results and Empirical Demonstration

A comprehensive battery of one- and two-dimensional tests is presented, including smooth flows for convergence studies, Riemann problems with strong shocks, blast and jet problems with ultrarelativistic speeds, Kelvin–Helmholtz instability, and interactions with density bubbles.

Key empirical results:

• Optimal order of accuracy is observed for smooth solutions in both 1D and 2D, including up to fifth-order ($O(\Delta x^5)$) convergence (see Tables 1 and 2 in the text).
• Strong shock problems: The blended LWFR scheme accurately resolves shock profiles, contacts, and rarefactions, with sharp resolution and no negative density/pressure, outperforming classic TVB-limited schemes in preserving smooth extrema and controlling post-shock oscillations.
• Jet and blast interaction cases: The high Lorentz-factor regimes are handled without loss of admissibility or appearance of nonphysical state variables.
• Kelvin–Helmholtz instability: The method resolves secondary structures without excessive numerical diffusion and suppresses unphysical “ringing,” even at high resolutions.

Comparisons with previous blending and limiting strategies (including those in [hennemann2021provably]) highlight the superior selectivity and robustness of the proposed indicator and subcell blending.

Figure 4: Density field zoomed in $[0.50,0.53]$ in the 1D blast wave configuration, showing resolution of narrow structures at high order.

Practical and Theoretical Implications

The work establishes that Jacobian-free high-order schemes with carefully designed blending and admissibility enforcement are effective for challenging RHD systems, especially when coupled with problem-adapted smoothness indicators. The method provides a practical pathway for high-fidelity yet physically robust simulation of relativistic flows on modern parallel architectures, both CPU and GPU, due to the single-step temporal update and polynomial quadrature-free flux calculation.

On the theoretical side, the introduction of concave-admissible regions directly in conservative variables and the modular design of the bound-preserving limiter generalize to other hyperbolic systems with physical constraints.

Figure 5: $\ln \rho$ field for a relativistic jet test, demonstrating robust jet structure capture even at very high Lorentz factors.

Future Directions

Potential extensions include:

• Application to general relativistic hydrodynamics (GRHD) with dynamic spacetime metrics.
• Adaptive mesh refinement exploiting the modular subcell blending structure for efficient resource allocation.
• Extension to relativistic magnetohydrodynamics (RMHD) with rigorous divergence control of magnetic field.
• Incorporation of more sophisticated, possibly ML-based, data-driven discontinuity indicators.
• Deployment on exascale platforms leveraging the low-communication, local-update character of the method.

Conclusion

This work demonstrates a comprehensive framework for high-order, physically admissible, and robust simulation of relativistic hydrodynamics via a bound-preserving Lax-Wendroff flux reconstruction scheme. The key contributions are the combination of a Jacobian-free update, a new smoothness-sensitive subcell blending strategy, a tailored discontinuity indicator, and mathematically rigorous admissibility enforcement. The scheme achieves high resolution and resilience against spurious oscillations across a wide set of demanding RHD applications, showing both practical value for computational astrophysics and significant methodological advancement.