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

This paper introduces efficient splitting methods, specifically Douglas-Rachford and Davis-Yin, to implement optimization-based limiters that enforce invariant-domain preservation in high-order numerical schemes for gas dynamics, demonstrating their robustness and performance through applications to discontinuous Galerkin methods on compressible flow benchmarks.

The Gist

Imagine you are a chef trying to bake the perfect cake (solving complex gas dynamics equations) in a high-pressure kitchen. You have a very precise recipe (mathematical equations) that tells you exactly how much flour, sugar, and eggs (density, momentum, and energy) to use.

However, your oven is a bit glitchy. Sometimes, when you mix the ingredients quickly (using high-speed computer simulations), the math gets confused and suggests you have negative sugar or negative eggs. In the real world, you can't have negative mass or energy. If your simulation tries to use "negative eggs," the whole cake collapses, the numbers explode, and the computer crashes. This is called losing the "invariant domain."

For years, chefs (scientists) have had to be very careful, cooking slowly (using tiny time steps) to avoid this, or using a blunt tool to just chop off any "bad" ingredients, which often ruined the taste (accuracy) of the cake.

This paper introduces a new, smart, efficient kitchen assistant that fixes the ingredients after the mixing is done but before the cake goes into the oven. Here is how it works, broken down into simple concepts:

1. The Problem: The "Bad Ingredient" Crisis

When simulating things like supersonic jets or explosions, the computer calculates the average amount of "stuff" in each little box of space. Sometimes, due to the speed of the calculation, a box might end up with a negative density (which is impossible).

• The Old Way: Stop the simulation, shrink the time step, and try again. This is like stopping the oven every 5 seconds to check the cake. It's safe, but it takes forever.
• The New Way: Let the computer cook fast, but when it makes a mistake (a negative number), this new assistant instantly fixes it without ruining the flavor.

2. The Solution: The "Optimization Limiter"

The authors created a method called an Optimization-Based Limiter. Think of this as a "Smart Corrector."

When the computer produces a "bad" box of ingredients (e.g., negative density), the Smart Corrector asks:

"What is the smallest, most polite change I can make to these numbers to make them legal (positive) again, while making sure I don't accidentally steal or create extra ingredients from other boxes?"

It solves a mini-math puzzle to find the perfect fix. It's like a chef who, upon seeing a bowl with -2 eggs, gently adds just enough positive eggs to make it 0, without taking eggs away from the flour bowl next to it.

3. The Secret Sauce: "Splitting Methods"

The hard part is solving that mini-math puzzle quickly. If you try to solve it all at once, it takes too long. The authors used a clever trick called Splitting Methods (specifically Douglas-Rachford and Davis-Yin splitting).

The Analogy: The Two-Step Dance
Imagine you are trying to park a car in a tight spot (finding the correct numbers).

1. Step A: You pull forward as close as you can to the front bumper (satisfying one rule, like keeping total energy constant).
2. Step B: You back up just enough to fit between the lines (satisfying the second rule, like keeping density positive).
3. Repeat: You do this dance forward and backward a few times. Each step gets you closer to the perfect parking spot.

The paper shows that by breaking the complex math problem into these simple, alternating steps, the computer can find the perfect solution in a fraction of a second. They even derived a special formula (like a shortcut) to handle the "parking" when the rules get really tricky (involving 3D space and complex gas physics).

4. Two Types of "Smart Correctors"

The paper tests two versions of this assistant:

• The "Smooth" Corrector ($L_2$ norm): This one tries to spread the correction out evenly. If one box is bad, it might slightly tweak a few neighbors to fix it. It's very fast and accurate for most cakes.
• The "Targeted" Corrector ($L_1$ norm): This one tries to fix the problem by changing as few boxes as possible. It's like saying, "I'll only touch the one bowl that has the negative eggs, and leave everything else exactly as is."
  • Surprise: For most cakes, the "Smooth" one is faster. But for a specific type of chaotic, high-speed jet (the "Astrophysical Jet"), the "Targeted" one is actually better because it triggers the fix less often, keeping the simulation running smoother overall.

5. Why This Matters

• Speed: You can run simulations much faster because you don't have to slow down the computer to avoid errors.
• Robustness: You can simulate extreme events (like supernovas or hypersonic missiles) that used to crash computers because the numbers went negative.
• Accuracy: The "Smart Corrector" is so gentle that it doesn't ruin the taste of the cake. The final result is still mathematically precise.

Summary

This paper gives us a mathematical safety net. It allows scientists to run high-speed, high-accuracy simulations of gas and explosions without worrying about the numbers going "off the rails." It uses a clever "dance" of math steps to instantly fix impossible numbers, ensuring the simulation stays physically realistic while running at top speed.

Technical Summary

Here is a detailed technical summary of the paper "Efficient optimization-based invariant-domain-preserving limiters in solving gas dynamics equations."

1. Problem Statement

High-order numerical schemes (such as Discontinuous Galerkin, DG) for solving gas dynamics equations (Compressible Euler and Navier-Stokes) often suffer from a critical issue: they can produce non-physical solutions where density ($\rho$) or internal energy ($e$) become negative or zero. This leads to a breakdown of the simulation, particularly in demanding scenarios involving low density, high-speed shocks, and blast waves.

To ensure nonlinear stability and physical meaningfulness, numerical solutions must remain within an invariant domain $G_\varepsilon$, defined as the set of states where density and internal energy are strictly positive (bounded below by a small $\varepsilon > 0$).

While existing methods (like the Zhang-Shu scaling limiter or Flux Corrected Transport) exist, they often rely on specific time discretizations (e.g., Strong Stability Preserving - SSP) or struggle with arbitrary high-order extensions for complex systems. The paper addresses the need for a robust, optimization-based limiter that:

1. Enforces the invariant domain constraint for vector variables (density, momentum, energy) simultaneously, not just scalar bounds.
2. Preserves global conservation.
3. Maintains high-order accuracy.
4. Works with arbitrary time-stepping schemes (including non-SSP methods).

2. Methodology

The core approach formulates the limiter as a constrained optimization problem. Given a numerical solution $U_h$ with cell averages that may violate the invariant domain, the method seeks a modified solution $X_h$ that minimizes the distance to $U_h$ while satisfying conservation and invariant domain constraints.

2.1. Optimization Formulation

The problem is defined as:
$$\min_{X_h} \|X_h - U_h\| \quad \text{subject to} \quad \int_\Omega X_h = \int_\Omega U_h \quad \text{and} \quad X_h|_{K_i} \in G_\varepsilon$$
The authors investigate two norms:

• $L_2$-norm minimization: Minimizes the Euclidean distance.
• $L_1$-norm minimization: Minimizes the sum of absolute differences (potentially leading to sparser modifications, though the paper notes this is not always the case for vector variables).

2.2. Operator Splitting Algorithms

To solve these constrained problems efficiently at every time step, the authors employ operator splitting methods:

• Douglas-Rachford Splitting (DRS): Used for $L_1$ minimization and as an outer loop for nested problems.
• Davis-Yin Splitting (DYS): A three-operator splitting method used for $L_2$ minimization. DYS is highlighted for its parameter-free nature (using step size $\gamma = 1/L$) and rapid convergence.

2.3. Key Technical Innovation: Explicit Projection

The most significant computational challenge is the projection onto the invariant domain set $G_\varepsilon$. Unlike simple scalar bounds (clipping), projecting a vector $(\rho, \mathbf{m}, E)$ onto the convex set defined by $\rho \ge \varepsilon$ and $E - \|\mathbf{m}\|^2/(2\rho) \ge \varepsilon$ is non-trivial.

• The authors derive an explicit analytical formula for this projection.
• The derivation uses Karush-Kuhn-Tucker (KKT) conditions, reducing the projection to solving a cubic equation (for the momentum component) and a quadratic equation (for density).
• These formulas are provided for 1D, 2D, and 3D cases in the appendices, allowing for efficient implementation without iterative sub-solvers for the projection itself.

2.4. Nested Approach for $L_1$ Norm

For the $L_1$ minimization, the proximal operator for the intersection of constraints is not analytically simple. The authors propose a nested splitting strategy:

• Outer Loop: DRS handles the $L_1$ norm and global conservation.
• Inner Loop: DYS is used to numerically compute the proximal operator required by the outer loop (which effectively solves the $L_2$ projection subproblem).

3. Key Contributions

1. Explicit Projection Formula: The derivation of an efficient, explicit formula for projecting a vector state onto the gas dynamics invariant domain $G_\varepsilon$. This is the first such derivation for the full vector system, enabling the use of splitting methods.
2. Splitting Method Framework: The successful application of DRS and DYS to vector-valued conservation laws. Specifically, the use of DYS for $L_2$ minimization and a DRS-DYS nested scheme for $L_1$ minimization.
3. General Applicability: The method is designed to work with any time-marching scheme (including non-SSP Runge-Kutta methods) and any high-order spatial discretization (DG, Finite Volume), provided global conservation is maintained.
4. Accuracy Preservation: Theoretical proof (Theorem 1) showing that the $L_2$-based limiter improves the accuracy of the DG solution (in the $L_2$ norm) compared to the unmodified solution, provided the exact solution is within the domain.

4. Numerical Results

The authors validated the method on a series of benchmarks:

• 1D Synthetic Tests:
  • Traveling Waves: Demonstrated that both $L_1$ and $L_2$ limiters successfully enforce bounds. The $L_2$ solver (DYS) converged faster (fewer iterations) than the $L_1$ solver (DRS).
  • Lax Shock Tube: Showed robustness in handling perturbations near shock waves.
• 2D Benchmarks:
  • Convergence Study: Verified that the limiters preserve the designed order of accuracy ($P_2$ and $P_3$) for smooth manufactured solutions.
  • Sedov Blast Wave: A strong shock test with low density. The method preserved the invariant domain and global conservation. The $L_2$ limiter was significantly cheaper computationally than the $L_1$ limiter.
  • Mach 2000 Astrophysical Jet: A highly demanding test with extreme Mach numbers using a non-SSP RK4 scheme.
    • Result: The method successfully prevented simulation blow-up.
    • Observation: While the $L_1$ solver was computationally more expensive per trigger, it was triggered less frequently during the time evolution compared to the $L_2$ limiter. This suggests $L_1$ might be preferable for specific high-speed flow problems despite the higher per-step cost.

5. Significance and Conclusion

This paper provides a unified and robust framework for enforcing physical constraints in high-order gas dynamics simulations.

• Robustness: It enables the use of high-order schemes in extreme conditions (low density, high Mach) where traditional limiters might fail or require overly restrictive time steps.
• Flexibility: By decoupling the limiter from the time-stepping scheme (working as a post-processor), it allows researchers to use advanced, non-SSP time integrators without sacrificing stability.
• Efficiency: The derivation of the explicit projection formula makes the computationally expensive optimization step feasible for large-scale 3D simulations.
• Trade-off: The study clarifies the trade-off between $L_1$ and $L_2$ limiters. While $L_2$ (via DYS) is generally faster and improves accuracy, $L_1$ (via nested DRS) may reduce the frequency of limiter activation in highly transient, high-speed flows, potentially offering better overall performance for specific astrophysical applications.

In summary, the work bridges the gap between convex optimization theory and computational fluid dynamics, offering a practical, mathematically rigorous tool for next-generation high-order gas dynamics solvers.