Efficient optimization-based invariant-domain-preserving limiters in solving gas dynamics equations

Basic Information

• Paper ID: 2510.21080
• Title: Efficient optimization-based invariant-domain-preserving limiters in solving gas dynamics equations
• Authors: Chen Liu (University of Arkansas), Dionysis Milesis (Boston University), Chi-Wang Shu (Brown University), Xiangxiong Zhang (Purdue University)
• Classification: math.NA (Numerical Analysis), cs.NA (Computational Science)
• Submission Date: October 24, 2025
• Paper Link: https://arxiv.org/abs/2510.21080v1

Abstract

This paper proposes efficient splitting methods based on optimization for invariant-domain-preserving limiters applied to high-order numerical schemes for gas dynamics equations. Core techniques include: (1) concise and efficient explicit formulas for invariant domain projection; (2) appropriate application of classical Douglas-Rachford splitting (DRS) and its extension Davis-Yin splitting (DYS). The method is applicable to various numerical schemes and can construct high-order accurate, globally conservative solvers for compressible flow equations that preserve the invariant domain. The authors validate their approach on high-order discontinuous Galerkin (DG) schemes, demonstrating the robustness and performance of ℓ¹ and ℓ² norm minimization limiters through stringent benchmark tests.

Research Background and Motivation

Core Problem

The compressible Euler and Navier-Stokes equations are fundamental models in gas dynamics with widespread applications in aerospace and astrophysics. Numerical solutions must guarantee positivity of density and internal energy (positivity), which is not only a physical requirement but also crucial for achieving nonlinear stability, particularly for extreme applications involving low density and low pressure (such as high-speed shock waves and detonation waves).

Problem Significance

For ideal gases, conservative variables are density ρ, momentum m, and total energy E, with internal energy satisfying ρe = E - ||m||²/(2ρ) and pressure p = (γ-1)ρe. The invariant domain that physical solutions must satisfy is: $G = \{U = [\rho, m^T, E]^T : \rho > 0, \rho e(U) = E - \frac{||m||^2}{2\rho} > 0\}$

Since ρe(U) is concave with respect to U, by Jensen's inequality, the set G is convex.

Limitations of Existing Methods

1. Restrictions of explicit methods: Most positivity-preserving methods (such as Zhang-Shu limiters) use fully explicit time discretization. For compressible NS equations, the time step is restricted to Δt = O(ReΔx²), applicable only for large Reynolds numbers.
2. Difficulty in high-order extension: Semi-implicit and fully implicit positivity-preserving schemes can use larger time steps (Δt = O(Δx)), but extending to arbitrary high-order accuracy is very difficult.
3. Inadequacy of existing optimization methods: Current optimization methods primarily handle bound-preserving problems for scalar variables and have not sufficiently addressed invariant domain constraints for vector variables.

Research Motivation

This paper proposes an optimization-based approach that solves constrained minimization problems to find minimal corrections to given numerical solutions while maintaining global conservation and invariant domain constraints. The key challenge is how to efficiently solve such constrained minimization problems (required at each time step).

Core Contributions

1. Explicit projection formula: First derivation of efficient explicit formulas for projecting vector variables onto the gas dynamics invariant domain Gε (by solving cubic equation roots), which is the foundation for implementing efficient splitting methods.
2. DYS method for ℓ² limiters: Proposes using Davis-Yin three-operator splitting (DYS) to efficiently solve ℓ² norm optimization problems without parameter tuning, typically converging to machine precision in a few iterations.
3. Nested splitting method for ℓ¹ limiters: Designs a method combining DRS with nested DYS to solve ℓ¹ norm optimization problems, using Douglas-Rachford splitting in the outer layer and DYS for numerically computing the proximal operator in the inner layer.
4. Theoretical accuracy guarantee: Proves a theorem (Theorem 1) that the ℓ² limiter improves DG solution accuracy: in the L² norm sense, the limited solution is closer to the exact solution than the original solution.
5. Broad applicability verification: Validates the method on high-order DG schemes with non-SSP Runge-Kutta time formats, demonstrating applicability to various time advancement schemes.

Method Details

Task Definition

Input: Cell average values Ūh of DG solution Uh (possibly violating invariant domain Gε)
Output: Corrected cell average values rh ∈ Gε
Constraints:

• Global conservation: ∫Ω rh = ∫Ω Ūh
• Invariant domain preservation: rh|Ki ∈ Gε, for all cells Ki

The numerical invariant domain is defined as (ε > 0 is a small constant): $G_\varepsilon = \{U = [\rho, m^T, E]^T : \rho \geq \varepsilon, \rho e(U) = E - \frac{||m||^2}{2\rho} \geq \varepsilon\}$

Optimization Problem Formulation

ℓ² model (Equation 8): $\min_{U_h} ||U_h - \bar{U}_h||_{L^2}^2 + \iota_{\Lambda_1}(U_h) + \iota_{\Lambda_2}(U_h)$

ℓ¹ model (Equation 9): $\min_{U_h} ||U_h - \bar{U}_h||_{L^1} + \iota_{\Lambda_1}(U_h) + \iota_{\Lambda_2}(U_h)$

Where:

• Λ₁ = {Uh : ∫Ω Uh = ∫Ω Ūh} (conservation constraint)
• Λ₂ = {Uh : Uh|Ki ∈ Gε, ∀i} (invariant domain constraint)
• ιΛ is the indicator function (0 inside the set, +∞ outside)

Core Technical Innovations

1. Explicit Formula for Invariant Domain Projection

This is the most critical technical contribution of the paper. Through Karush-Kuhn-Tucker (KKT) conditions, the projection problem is transformed into:

One-dimensional case (Appendix B): Given u, v, wᵀ ∉ Gε, find the projection ρ, m, Eᵀ. Based on the signs of Lagrange multipliers λ and μ, four cases are distinguished:

• Case 1 (μ=0, λ>0): ρ, m, Eᵀ = ε, v, wᵀ
• Case 2 (μ=0, λ=0): The point itself lies in Gε
• Case 3 (μ>0, λ>0): Solve cubic equation m³ + (4ε² - 2εw)m - 2ε²v = 0
• Case 4 (μ>0, λ=0): Solve quadratic equation for ρ, then compute m and E

Two-dimensional case (Appendix C): Similar analysis but handling two momentum components, ultimately reducing to solving cubic and quadratic equations.

Key observation: All real roots can be obtained using Cardano's formula (Appendix D) through real arithmetic, avoiding complex arithmetic and simplifying implementation.

2. Davis-Yin Splitting for ℓ² Problems

For the ℓ² model, function decomposition (Equation 36):

• f(X) = ιΛ₁(X) (conservation constraint)
• g(X) = ιΛ₂(X) (invariant domain constraint)
• h(X) = (1/2α)||X - U||²F (ℓ² term)

DYS iteration scheme (step size γ = α = 1/L):

X^{k+1/2} = prox_γf(z^k)           // conservation constraint projection
X^{k+1} = prox_γg(2X^{k+1/2} - z^k - γ∇h(X^{k+1/2}))  // invariant domain projection
z^{k+1} = z^k + X^{k+1} - X^{k+1/2}

Advantages:

• No parameter tuning required (fixed step size γ = α)
• Fast convergence (typically <20 iterations)
• Only two projections per iteration

3. Nested DRS-DYS for ℓ¹ Problems

For the ℓ¹ model, function decomposition (Equation 45):

• f(X) = ||X - U||L¹
• h(X) = ιΛ₁(X) + ιΛ₂(X)

Outer layer DRS:

w^{k+1} = λ prox_γf(2x^k - w^k) + w^k - λx^k
x^{k+1} = prox_γh(w^{k+1})

Inner layer DYS: prox_γh is computed numerically by solving the ℓ² subproblem (Equation 49): $\min_Z \frac{1}{2\gamma}||Z - X||²_F + \iota_{\Lambda_1}(Z) + \iota_{\Lambda_2}(Z)$

Proximal operator formula:

• prox_γf corresponds to the shrinkage operator: $[\text{prox}_{\gamma f}(x)]_i = u_i + S_\gamma(x_i - u_i)$ where $S_\gamma(a) = \text{sgn}(a)\max\{|a| - \gamma, 0\}$

Algorithm Flow

Complete algorithm (applied after each RK stage):

1. Compute cell average values Ūh of DG solution
2. Check if any cell average values violate Gε
3. If violations exist, solve optimization problem (8 or 9) to obtain rh
4. Correct DG polynomial: Uh ← (Uh - Ūh) + rh
5. Apply Zhang-Shu limiter to ensure quadrature points also lie in Gε

Stopping criterion: ||z^{k+1} - z^k||₂ₕ < ε (ε = 10⁻¹³)

Experimental Setup

Datasets and Test Problems

One-dimensional comprehensive tests

1. Traveling wave problem (Example 5.1): Linear advection equation with triangular and square waves, verifying bound preservation
2. Perturbed Lax shock tube (Example 5.2): Perturbed Euler equation exact solution, testing invariant domain preservation

Two-dimensional benchmark tests

1. Convergence study (Example 5.3): Method of manufactured solutions, verifying spatial convergence order
   • Domain: Ω = 0,1², final time T = 0.1
   • Meshes: Δx = 1/25, 1/50, 1/100
   • Basis functions: P² and P³
2. Sedov blast wave (Example 5.4): Strong shock wave from point source explosion
   • Domain: Ω = 0, 1.1², T = 1
   • Mesh: 160×160 uniform grid, P² basis functions
   • CFL = 0.2 (larger than required for weak positivity)
3. Mach 2000 astrophysical jet (Example 5.5): Extremely high Mach number flow
   • Domain: Ω = 0,1×-0.5,0.5, T = 0.001
   • Mesh: 640×640, Q³ DG spectral element method
   • Non-SSP classical fourth-order RK4 time format

Evaluation Metrics

1. Iteration count: Number of DYS/DRS iterations for convergence
2. Projection count: Total number of invariant domain projection operator evaluations (for fair cost comparison)
3. Convergence rate: Spatial discretization error convergence order (L² and L¹ norms)
4. Accuracy error: ||ρⁿₕ - ρ(tⁿ)||{L²ₕ} and ||ρⁿₕ - ρ(tⁿ)||{L¹ₕ}

Implementation Details

• Parameter settings:
  • ε = 10⁻¹³ (most tests), ε = 10⁻⁸ (astrophysical jet)
  • Convergence tolerance: ε = 10⁻¹³ (or 10⁻⁸)
  • DRS step size: γ = 10⁻⁷ (2D), selected through tuning
  • DYS step size: γ = 1/L = α (fixed, no tuning required)
• Comparison methods:
  • ClipAndAssuredSum (direct method, ℓ¹ scalar case only)
  • Zhang-Shu limiter (applied to quadrature points)

Experimental Results

Main Results

1. One-dimensional traveling wave test (Figure 1-2)

• ℓ² limiter (DYS): Converges within 60 iterations (all time steps)
• ℓ¹ limiter (DRS): Converges within 200 iterations
• Projection count: ℓ¹ method significantly higher due to inner DYS
• Solution quality: Three methods (ℓ², ℓ¹-DRS, ℓ¹-direct) produce consistent results within precision ε
• Convergence: Both methods exhibit asymptotic linear convergence

2. Perturbed Lax shock tube (Figure 3)

• ℓ² limiter: Converges within 20 iterations
• ℓ¹ limiter: 200 DRS iterations, but significantly higher total projection count
• Key finding: ℓ¹ minimizer is not sparse (unlike other applications)

3. Convergence study (Tables 1-2)

P² basis functions (ℓ² limiter):

P³ basis functions (ℓ² limiter): Achieves 7th-order convergence rate (initial mesh refinement)

Conclusion: Optimization limiters preserve optimal convergence order without degrading high-order accuracy.

4. Sedov blast wave (Figure 4)

• Trigger frequency: Cell average limiter triggered only at 1 time step (step 74669) throughout time evolution
• ℓ² limiter: Approximately 10 DYS iterations required at that time step
• ℓ¹ limiter: Approximately 30 DRS iterations required, higher total projection count
• Physical result: Shock position correctly captured, reasonable density field (50 contour lines, range 0.001-6)

5. Mach 2000 astrophysical jet (Figure 5)

Key findings:

• Trigger frequency difference: ℓ¹ limiter triggered significantly fewer times during time evolution than ℓ² limiter
• Single trigger cost: ℓ¹ limiter higher per-trigger cost (more projections required)
• Total cost: Due to fewer triggers, ℓ¹ total cost comparable to ℓ²
• Applicability verification: Non-SSP RK4 method successfully applied, demonstrating broad applicability

Ablation Studies

While not explicitly labeled as "ablation studies," the paper analyzes component contributions through the following comparisons:

1. ℓ² vs ℓ¹ norm choice:
   • ℓ²: Faster computation, accuracy improvement theoretically guaranteed (Theorem 1)
   • ℓ¹: Fewer triggers in some problems (e.g., astrophysical jet)
2. Splitting method choice:
   • DYS (ℓ²): Parameter-free, fast convergence
   • DRS nested DYS (ℓ¹): Flexible but higher computational cost
3. Limiter domain choice (Equation 56):
   • Apply limiter only to shock-affected regions
   • Improves efficiency and robustness

Experimental Findings

1. ℓ¹ non-sparsity (Remark 1): Unlike many other applications, ℓ¹ minimizers in this problem are not sparser than ℓ², with similar numbers of modified cells.
2. Accuracy improvement theorem verification: Theorem 1's theoretical prediction (ℓ² limiter improves accuracy) is verified in numerical experiments.
3. Parameter sensitivity: DRS requires tuning (γ), while DYS with fixed step size γ = 1/L performs well.
4. Computational efficiency:
   • One-dimensional scalar: Existing direct methods (ClipAndAssuredSum) are optimal
   • Vector invariant domain: DYS is the first efficient practical method

Related Work

Main Research Directions

1. Traditional positivity-preserving methods

• Zhang-Shu limiter 45,46: Simple scaling limiter, requires SSP time format, accuracy proven
• FCT method 20,43: Convex combination of high-low order fluxes, flexible but difficult accuracy proof
• Convex limiting and sub-cell methods 22,32,37,25: Recent developments, broad applicability

2. Optimization methods (scalar variables)

• Guba et al. 17: Optimization limiters for spectral element methods
• Bochev et al. 2,3: Constrained optimization for data transfer
• Bradley et al. 5: Communication-efficient tracer transport
• Liu et al. 29,27,28: Bound-preserving for Cahn-Hilliard-NS systems, compressible NS

Limitations: All existing optimization methods handle only scalar variables, not directly addressing invariant domain constraints for vector variables.

3. Implicit and semi-implicit methods

• Guermond et al. 18: Second-order invariant domain preserving approximation
• Liu-Zhang 31: Positivity-preserving semi-implicit DG

Advantages of This Paper

1. First to address vector invariant domains: Directly constrains Gε rather than decomposing into scalar constraints
2. Explicit projection formula: First derivation of efficient projection for gas dynamics invariant domain
3. First application of three-operator splitting: DYS applied to this class of problems for the first time
4. Broad applicability: Applicable to non-SSP time formats

Conclusions and Discussion

Main Conclusions

1. Method effectiveness: The proposed optimization-based invariant-domain-preserving limiters are effective in high-order DG schemes, applicable to various time advancement formats (including non-SSP).
2. ℓ² limiter advantages:
   • Lower computational cost (DYS converges quickly)
   • Theoretical accuracy guarantee (Theorem 1)
   • First choice in most cases
3. ℓ¹ limiter value:
   • Fewer triggers in some problems (e.g., astrophysical jets)
   • Total cost comparable to ℓ²
   • Superior in specific application scenarios
4. Global conservation: Although only maintaining global conservation (not local), numerical tests show no erroneous shock positions.

Limitations

1. Lack of local conservation: Optimization limiters maintain only global conservation, potentially limiting certain theoretical analyses.
2. Parameter tuning (ℓ¹): DRS method requires tuning, lacking optimal parameter formulas for vector cases (formulas exist for scalar cases 29).
3. Computational cost: Compared to traditional limiters (e.g., Zhang-Shu), optimization methods have higher computational cost but offer greater flexibility.
4. Three-dimensional extension: While the method directly extends to 3D, projection formula derivation is more complex, not detailed in the paper.
5. Sparsity: ℓ¹ limiter does not produce sparse solutions in this problem, unlike many other optimization applications.

Future Directions

1. Parameter optimization theory: Establish optimal parameter selection theory for DRS with vector invariant domains (similar to scalar case 29).
2. Local conservation improvement: Explore optimization limiter variants preserving local conservation.
3. Other PDE systems: Extend to other systems with invariant domain constraints (e.g., MHD equations).
4. Acceleration techniques: Research preconditioning and acceleration techniques to further improve convergence speed.
5. Adaptive methods: Automatically select ℓ¹ or ℓ² limiters based on problem characteristics.

In-Depth Evaluation

Strengths

1. Significant theoretical contributions

• High originality: First systematic study of optimization limiters for vector variable invariant domains
• Rigorous formula derivation: Complete derivation of projection formulas through KKT conditions (Appendices B, C), mathematically rigorous
• Accuracy theorem: Theorem 1 provides theoretical guarantee for ℓ² limiter accuracy improvement, beyond simply not degrading accuracy (Equation 5 vs 6)

2. Clever method design

• Appropriate splitting strategy: DYS for ℓ², nested DRS-DYS for ℓ¹, fully leveraging each method's advantages
• Efficient projection formula: Reduces to cubic equation root-finding, avoids complex arithmetic, simple implementation
• Parameter-free: DYS requires no tuning, strong practical utility

3. Comprehensive experimental design

• Multi-scale verification: From one-dimensional comprehensive tests to two-dimensional extreme flows (Mach 2000)
• Rigorous convergence verification: Manufactured solution tests confirm optimal convergence order
• Application-oriented: Sedov blast wave and astrophysical jet are recognized stringent tests in the field

4. Clear and standard writing

• Logical structure: Background → Method → Experiments, clear logic
• Complete details: Projection algorithms, parameter settings, implementation details all thorough
• Rich appendices: Technical derivations in appendices maintain readability of main text

Weaknesses

1. Limited theoretical analysis

• ℓ¹ convergence: No convergence rate analysis for ℓ¹ limiter, only numerical observations
• Parameter selection: DRS parameter γ depends on empirical tuning, lacking theoretical guidance
• Insufficient sparsity explanation: ℓ¹ non-sparsity phenomenon lacks deep theoretical analysis (only brief Remark 1)

2. Experimental setup defects

• Missing 3D: All tests are one- or two-dimensional, insufficient verification for 3D practical applications
• Lack of statistical significance: No multiple-run statistics reported (e.g., standard deviations)
• Insufficient cost quantification: Only iteration and projection counts reported, lacking actual CPU time comparisons
• Missing comparisons: No direct comparison with recently popular FCT-type methods

3. Method limitations

• Computational cost: Significantly higher per-step cost compared to Zhang-Shu and other traditional methods
• Scalability unknown: Large-scale parallel efficiency, GPU acceleration not addressed
• Robustness boundaries: Performance near extreme cases (e.g., near vacuum) insufficiently tested

4. Technical detail issues

• Limiter domain strategy (Equation 56): Ad-hoc design lacking theoretical support
• ε selection: Unclear criteria for choosing numerical invariant domain parameter ε
• Inner DYS stopping criterion: How inner precision affects outer convergence in nested methods not discussed

Impact

Contribution to the field

• Pioneering work: Establishes foundational framework for optimization limiters with vector invariant domain constraints
• Methodological significance: Demonstrates potential of convex optimization splitting methods in numerical PDEs
• Practical value: Provides positivity-preserving tool applicable to non-SSP formats

Potential influence

• Citation value: Likely to become important reference for gas dynamics numerical methods
• Extensibility: Method framework generalizable to other conservation law systems
• Software implementation: Clear algorithms in appendices facilitate community implementation and reproduction

Limitations

• Computational cost barrier: May limit application to scenarios requiring extreme accuracy
• Parameter tuning burden: ℓ¹ method parameter selection may hinder broad adoption
• 3D extension workload: Practical applications require additional 3D implementation and verification

Applicable Scenarios

Ideal application scenarios

1. High accuracy requirements: Scientific computing prioritizing accuracy over computational efficiency
2. Non-SSP time formats: Problems requiring implicit or non-SSP high-order RK methods
3. Extreme flows: Low-density, high-Mach situations where traditional methods fail
4. Research codes: Prototype algorithm development and method verification

Unsuitable scenarios

1. Industrial large-scale computing: Computational cost potentially prohibitive
2. Real-time simulation: Optimization iteration overhead unacceptable
3. Low-accuracy applications: Traditional methods (e.g., Zhang-Shu) more economical

Improvement suggestions

1. Hybrid strategy: Use traditional methods most time steps, enable optimization limiters only when necessary
2. Preconditioning: Research problem-specific preconditioning techniques to accelerate convergence
3. Machine learning: Use ML to predict optimal parameters or initial guesses

Key References

1. 46 Zhang & Shu (2010): On positivity-preserving high order DG schemes for compressible Euler equations - Classical Zhang-Shu limiter
2. 44 Zhang (2017): On positivity-preserving high order DG schemes for compressible NS equations - Accuracy proof and weak positivity theory
3. 9 Davis & Yin (2017): A three-operator splitting scheme - Theoretical foundation of DYS method
4. 26 Lions & Mercier (1979): Splitting algorithms for the sum of two nonlinear operators - Convex optimization extension of DRS
5. 5 Bradley et al. (2019): Communication-efficient property preservation - Theoretical results for ℓ¹ scalar case
6. 29 Liu et al. (2024): Optimization-based bound-preserving limiter for Cahn-Hilliard-NS - Optimal DRS parameters for scalar case
7. 18 Guermond et al. (2021): Second-order invariant domain preserving approximation - Theoretical framework for invariant domain preservation

Overall Assessment: This is a high-quality numerical analysis paper making pioneering contributions to optimization limiters with vector invariant domain constraints. The method design is clever, theoretical derivations rigorous, and experimental verification comprehensive. Main innovations lie in explicit projection formulas and efficient splitting method combinations. Limitations include computational cost, parameter tuning, and insufficient 3D verification. Valuable for high-accuracy scientific computing and method research, though industrial large-scale applications may require further optimization. The paper is clearly written with complete technical details and good reproducibility. Recommended for researchers working on high-order numerical methods and compressible flow computation.