Entropy-Stable Relaxation Solvers

• Entropy-Stable Relaxation-Based Solvers are advanced numerical schemes that use adaptive relaxation parameters to enforce discrete entropy control in nonlinear hyperbolic systems.
• They combine relaxation approximations with high-order discretization methods like DGSEM and SBP to achieve robust stability, conservation, and accuracy.
• These solvers are applied to complex flows such as the compressible Euler and Navier-Stokes equations, ensuring positivity preservation and shock stability.

Entropy-stable relaxation-based solvers are a class of high-resolution numerical schemes for systems of conservation and balance laws, which achieve discrete entropy conservation or dissipation by coupling relaxation approximations with carefully designed time and space discretizations. Central to this approach is the enforcement of entropy stability at the fully-discrete level, often via relaxation parameters that are solution-adaptive and computed to ensure that prescribed convex entropy functionals are non-increasing (or exactly conserved) in time. These solvers deliver high-order accuracy, nonlinear robustness, and preservation of positivity and primary invariants for complex nonlinear systems, including the compressible Euler and Navier-Stokes equations, hyperbolic systems with source terms, and problems on unstructured meshes in multiple spatial dimensions (Renac, 2020, Bellotti et al., 2024, Ranocha et al., 2019, Fernandez et al., 2019, Yan et al., 2021, Ranocha et al., 2020, Izgin et al., 2 Apr 2026, Miroshnikov et al., 2014, Tzavaras, 2014, Ranocha et al., 2023, Doehring et al., 7 Jul 2025, Kang et al., 2021).

1. Theoretical Foundations: Entropy, Relaxation, and Stability

Entropy-stable relaxation-based methods are built upon the existence of a convex entropy-entropy flux pair $(\eta(u), q(u))$ for the underlying hyperbolic system of conservation or balance laws. This structure requires that

$$\nabla \eta(u) \cdot \nabla f^j(u) = \nabla G^j(u), \quad j=1,\dots,d,$$

where $f^j(u)$ are the fluxes and $G^j(u)$ the corresponding entropy fluxes (Renac, 2020, Bellotti et al., 2024). The entropy variables, $w = \partial_u \eta(u)$, provide a nonlinear symmetrizer for the system.

Relaxation approximations recast the original system as a larger augmented system with additional variables and fast relaxation rates that drive the solution toward equilibrium. Entropy-stable relaxation systems are constructed so that the entropy function for the augmented variables is globally convex (or polyconvex), and entropy dissipation is enforced via source (relaxation) terms (Miroshnikov et al., 2014, Tzavaras, 2014).

In discrete time and space, enforcing an entropy inequality at the numerical level typically requires introducing a single or multiple relaxation parameters per time step (or per element/subdomain) that scale the high-order update to ensure $\eta(u^{n+1}) \leq \eta(u^n)$. For dissipative systems, this ensures unconditional discrete entropy stability, whereas for conservative systems, the method yields exact discrete entropy conservation (Ranocha et al., 2019, Ranocha et al., 2023, Ranocha et al., 2020).

2. Relaxation-Based Fluxes and Finite Volume Schemes

A central ingredient in entropy-stable relaxation-based spatial discretizations is the use of relaxation-based approximate Riemann solvers. For the multicomponent compressible Euler system, such as in Renac (Renac, 2020), a pressure-relaxation (energy-relaxation) system is constructed, admitting a convex extended entropy and allowing for the definition of entropy-conservative two-point and entropy-stable three-point numerical fluxes.

• Three-point finite volume update: For the relaxation system,

$$W_j^{n+1} = W_j^n - \frac{\Delta t}{\Delta x}\left[H(W_j, W_{j + 1}) - H(W_{j-1}, W_j)\right],$$

where $H$ is constructed from the exact solution to the relaxation Riemann problem and satisfies a discrete entropy inequality (Renac, 2020).

• Discrete entropy inequality: The scheme satisfies

$$\rho \zeta(W_j^{n+1}) - \rho\zeta(W_j^n) + \frac{\Delta t}{\Delta x}\left[Z(W_j, W_{j + 1}) - Z(W_{j-1}, W_j)\right] \leq 0,$$

ensuring entropy dissipation at the volume level.

The subcharacteristic condition and positivity of cell averages are enforced at the flux (interface) level and through local CFL-type step-size constraints (Renac, 2020).

3. High-Order Entropy-Stable DG and SBP Frameworks

High-order entropy-stable solvers marry relaxation-based fluxes with discontinuous Galerkin spectral element methods (DGSEM), summation-by-parts (SBP) finite difference methods, or discontinuous Galerkin difference (DGD) operators (Renac, 2020, Fernandez et al., 2019, Yan et al., 2021). Central to these approaches are:

• SBP operators on Gauss-Lobatto grids: These provide discrete mimetic integration by parts and enable high-order, provably conservative and entropy-stable spatial discretizations.
• Interface coupling and entropy-dissipative SATs: Across element boundaries (including $h/p$-nonconforming interfaces), entropy-stable interface coupling is achieved using decoupled interpolations or interface flux corrections. Dissipative simultaneous-approximation-term (SAT) penalties ensure entropy-stability and conservation across arbitrary mesh topologies (Fernandez et al., 2019, Yan et al., 2021).
• Volume and interface flux structure:
  • Inside elements: entropy-conservative two-point fluxes (satisfying Tadmor's condition) are used for volume integrals.
  • At interfaces: relaxation-based three-point entropy-stable fluxes are deployed, yielding exact entropy cancellation in the semi-discrete analysis (Renac, 2020).

These architectures generalize to multidimensional, unstructured, and nonconforming grids, with formal proofs of semi-discrete and fully-discrete entropy inequalities (Fernandez et al., 2019, Yan et al., 2021).

4. Relaxation-Based Entropy-Stable Time Integration

Time integration uses entropy-stable relaxation-parameter techniques, either as modifications of standard explicit and implicit Runge-Kutta (RK) methods or as multiderivative time integrators (Ranocha et al., 2019, Ranocha et al., 2023, Ranocha et al., 2020, Kang et al., 2021). At each step:

• Standard RK update: $$\nabla \eta(u) \cdot \nabla f^j(u) = \nabla G^j(u), \quad j=1,\dots,d,$$0.
• Relaxed update: $$\nabla \eta(u) \cdot \nabla f^j(u) = \nabla G^j(u), \quad j=1,\dots,d,$$1, where $$\nabla \eta(u) \cdot \nabla f^j(u) = \nabla G^j(u), \quad j=1,\dots,d,$$2 is computed by solving the one-dimensional nonlinear equation:

$$\nabla \eta(u) \cdot \nabla f^j(u) = \nabla G^j(u), \quad j=1,\dots,d,$$3

with $$\nabla \eta(u) \cdot \nabla f^j(u) = \nabla G^j(u), \quad j=1,\dots,d,$$4 and $$\nabla \eta(u) \cdot \nabla f^j(u) = \nabla G^j(u), \quad j=1,\dots,d,$$5 (Ranocha et al., 2019).

• Multiderivative and positivity-preserving extensions: The relaxation approach extends to MDRK methods, IMEX and multirate integrators, and positivity-preserving Patankar-type schemes, where relaxation parameters can be coupled to ensure strict positivity as well as discrete entropy conservation (Ranocha et al., 2023, Izgin et al., 2 Apr 2026, Kang et al., 2021).
• Entropy guarantee: If $$\nabla \eta(u) \cdot \nabla f^j(u) = \nabla G^j(u), \quad j=1,\dots,d,$$6, the method is strictly entropy dissipative; for conservative problems, setting the right-hand side to zero yields machine-precision conservation.

Relaxation parameters are always solved at low (scalar or small system) cost per time-step, and strong stability preservation (SSP) and primary conservation properties of the base integrators are retained (Ranocha et al., 2019, Doehring et al., 7 Jul 2025).

5. Applications and Generalizations

The entropy-stable relaxation-based blueprint is applied to a range of complex hyperbolic systems, including:

• Compressible multicomponent Euler and Navier-Stokes equations: With strict positivity and entropy dissipation for partial densities and energy, robustly resolving material interfaces and shocks (Renac, 2020, Fernandez et al., 2019, Ranocha et al., 2020).
• General hyperbolic balance laws and source terms: Via convex extension of entropy structures and careful handling of stiff relaxation and source terms, including applications to polyconvex elastodynamics and reactive flows (Miroshnikov et al., 2014, Tzavaras, 2014).
• High-order lattice Boltzmann methods: Fourth-order entropy-stable composition via palindromic operator splitting, with adaptive stage-wise relaxation parameter selection (Bellotti et al., 2024).
• h/p-nonconforming and multirate/partitioned time integration: Allowing for efficient, scalable, and entropy-stable simulations on highly adaptive and heterogeneous meshes (Fernandez et al., 2019, Doehring et al., 7 Jul 2025).
• Positivity-preserving relaxation for ODEs/PDEs: Coupling classical Patankar and nonstandard splitting schemes with entropy-stable relaxation envelopes for unconditionally positive, conservative, and entropy-dissipative updates (Izgin et al., 2 Apr 2026).

In all cases, practical schemes achieve optimal high-order convergence in smooth regimes and stability in the presence of discontinuities, as verified in canonical and realistic multicomponent compressible flows, viscous shocks, turbulence, and complex geometries (Renac, 2020, Fernandez et al., 2019, Ranocha et al., 2020, Yan et al., 2021, Doehring et al., 7 Jul 2025).

6. Main Theorems and Theoretical Guarantees

The primary analytical and algorithmic guarantees established in entropy-stable relaxation-based solvers include:

Theorem	Guarantee	Reference
Existence/uniqueness of entropy relaxation	There is a unique $$\nabla \eta(u) \cdot \nabla f^j(u) = \nabla G^j(u), \quad j=1,\dots,d,$$7 enforcing $$\nabla \eta(u) \cdot \nabla f^j(u) = \nabla G^j(u), \quad j=1,\dots,d,$$8 for small enough $$\nabla \eta(u) \cdot \nabla f^j(u) = \nabla G^j(u), \quad j=1,\dots,d,$$9	(Ranocha et al., 2019, Ranocha et al., 2023, Izgin et al., 2 Apr 2026)
Discrete entropy-stability/inequality	Fully-discrete update is entropy-dissipative at each step (global/local)	(Renac, 2020, Ranocha et al., 2020, Yan et al., 2021, Fernandez et al., 2019)
Conservation of positivity	Positivity of key variables (density, energy, species) is preserved under CFL/step restrictions	(Renac, 2020, Izgin et al., 2 Apr 2026)
High-order convergence	Relaxed methods preserve the (base) high order of accuracy in smooth regimes	(Ranocha et al., 2023, Fernandez et al., 2019, Ranocha et al., 2020, Doehring et al., 7 Jul 2025)
Local conservation	Primary conserved quantities (mass, momentum, energy) are maintained exactly	(Ranocha et al., 2020, Ranocha et al., 2019, Fernandez et al., 2019, Yan et al., 2021, Doehring et al., 7 Jul 2025)
Robustness to shocks/discontinuities	Entropy-stable relaxation suppresses unphysical oscillations and ensures stability across discontinuities	(Renac, 2020, Bellotti et al., 2024, Ranocha et al., 2020, Doehring et al., 7 Jul 2025)

These results are underpinned by the convexity of the entropy, proper subcharacteristic and positivity conditions, and the asymptotic stability properties derived in the theory of relative entropy and relaxation systems (Miroshnikov et al., 2014, Tzavaras, 2014, Beretta, 2013).

7. Implementation and Practical Considerations

Implementing an entropy-stable relaxation-based solver involves:

• Replacing standard volume and interface fluxes with entropy-conservative and relaxation-based entropy-stable fluxes, respectively.
• Introducing and solving low-dimensional (often scalar) nonlinear equations for relaxation parameters at each time-step/subdomain/element, using standard root-finding algorithms. These computations are embarrassingly parallel across the mesh (Ranocha et al., 2020, Fernandez et al., 2019, Izgin et al., 2 Apr 2026).
• Employing existing high-order explicit (SSP) Runge-Kutta, multistage, MDRK, IMEX, or multirate integrators, modified only in the final update via the relaxation envelope (Ranocha et al., 2019, Ranocha et al., 2023, Doehring et al., 7 Jul 2025, Kang et al., 2021).
• For positivity, using dense-output style updates or Patankar/M-matrix corrections (Izgin et al., 2 Apr 2026).
• Ensuring that conservative and dissipative properties are rigorously enforced at every intermediate step.

Algorithmic overhead is minor compared to flux and residual evaluations, and stability/robustness is significantly enhanced in all reported benchmarks (Ranocha et al., 2020, Doehring et al., 7 Jul 2025, Izgin et al., 2 Apr 2026, Renac, 2020). The methods are compatible with low-storage and parallel implementations, including GPU-resident codes (Bellotti et al., 2024).

---

These ingredients establish entropy-stable relaxation-based solvers as a comprehensive, theoretically grounded, and practically robust framework for high-order, entropy-controlled simulations of nonlinear hyperbolic systems. Their widespread applicability across physical models, mesh types, and time-integration strategies makes them a central tool for advanced computational fluid dynamics and applied mathematics (Renac, 2020, Fernandez et al., 2019, Ranocha et al., 2019, Ranocha et al., 2020, Yan et al., 2021, Izgin et al., 2 Apr 2026, Miroshnikov et al., 2014, Ranocha et al., 2023, Doehring et al., 7 Jul 2025, Kang et al., 2021).