Invariant domain preserving limiting of time explicit and time implicit discretizations for systems of conservation laws

Technical Summary: Invariant Domain Preserving Limiting of Time Explicit and Time Implicit Discretizations

Problem Statement
The paper addresses the challenge of designing high-order (HO) numerical approximations for nonlinear hyperbolic systems of conservation laws that preserve invariant domains (ID). While HO methods accurately capture smooth solutions, they often fail to satisfy physical constraints at the discrete level, such as the positivity of density and internal energy or the satisfaction of entropy inequalities. Conversely, low-order (LO) schemes (e.g., Godunov, Lax-Friedrichs) are known to preserve these invariant domains but suffer from excessive numerical diffusion. The objective is to construct a limiting technique that restricts the HO solution to remain within the convex set of admissible states $D$ (e.g., positive density, positive internal energy, entropy bounds) without sacrificing conservation or the high-order accuracy in smooth regions.

Methodology
The authors propose a general limiting framework applicable to any conservative spatial discretization (Finite Volume, Discontinuous Galerkin) and a wide range of time integration schemes, including explicit and implicit Runge-Kutta (RK), multistep methods, and Discontinuous Galerkin (DG) time integration.

1. Core Concept: The method limits the HO solution around a LO solution that is guaranteed to be ID-preserving. It generalizes the Flux-Corrected Transport (FCT) framework by limiting antidiffusive fluxes (the difference between HO and LO residuals) but defines the limiting coefficients to express the final solution as a convex combination of ID-preserving quantities, similar to the convex limiting framework.
2. Limiting Coefficients: Unlike standard FCT which imposes bounds on solution components, this approach defines coefficients $l_{\kappa,\kappa'}^{i,j}$ such that the limited solution $u_h^{n+1}$ is a convex combination of the LO solution and the HO update. The coefficients are determined by finding the largest scaling factors $k$ such that the intermediate states remain within the invariant domain $D$.
3. Iterative Procedure: The limiter is embedded in an iterative algorithm. Starting from the LO solution, the method iteratively refines the limited solution by re-applying the limiter to the difference between the HO solution and the current limited solution. This process continues until a tolerance is met or a maximum number of iterations is reached, aiming to recover the HO solution as closely as possible while maintaining IDP properties.
4. Acceleration: A heuristic acceleration factor $\beta \geq 1$ is introduced to scale the limiting coefficients, accelerating the convergence of the iterative process toward the HO solution.
5. Time Integration Adaptation: The framework adapts to various time schemes by defining a single-step LO scheme (Forward-Euler for explicit, Backward-Euler for implicit) to compute the necessary residuals and antidiffusive fluxes, avoiding the need to recompute LO solutions at every stage of multi-stage or multi-step methods.
6. Spatial Discretization: Specific derivations for antidiffusive fluxes are provided for Finite Volume (FV) schemes, general DG schemes, and the Discontinuous Galerkin Spectral Element Method (DGSEM). For DGSEM, the method combines cell-averaged limiting with a linear scaling limiter to ensure nodal values also satisfy the invariant domain constraints.

Key Contributions

• Generalization to Systems and Implicit Schemes: The method extends convex limiting concepts to systems of conservation laws and, crucially, applies to time implicit discretizations, a domain where previous convex limiting approaches were restricted to explicit schemes.
• Flux-Based Limiting without Bar States: The approach defines limiting coefficients to form convex combinations of LO and HO solutions without relying on "bar states" (intermediate states used in some convex limiting formulations). This allows the use of any consistent and conservative numerical flux in the HO scheme.
• Iterative Refinement: The introduction of an iterative procedure allows the limited solution to converge toward the HO solution, recovering accuracy in smooth regions while strictly preserving invariant domains.
• Acceleration Heuristic: A simple scaling factor $\beta$ is proposed to significantly reduce the number of iterations required for convergence.
• Broad Applicability: The derivation is presented for general conservative discretizations, with specific implementations detailed for FV, DG, and DGSEM, coupled with explicit/implicit RK, multistep, and DG time integrators.

Numerical Results
The paper validates the method through numerical experiments on 1D and 2D problems involving scalar hyperbolic equations and the compressible Euler equations:

• Riemann Problems (Strong Shocks): For the compressible Euler equations, the HO scheme without limiting fails due to negative pressure/density. The proposed limiter successfully resolves strong shocks sharply while maintaining positivity and entropy bounds for both explicit (ERK3) and implicit (DIRK33) time stepping, even at large CFL numbers (up to 10).
• Kelvin-Helmholtz Instability: In 2D simulations with large time steps, the LO scheme damps small scales excessively. The limited scheme recovers the high resolution of the HO scheme, preserving fine-scale structures.
• Shock Diffraction and Step Flow: The limiter demonstrates robustness in complex 2D geometries, maintaining physical admissibility where the unlimited HO scheme fails. It achieves resolution comparable to explicit schemes with small time steps but with the stability of implicit methods.
• Scalar Equations: Tests on linear transport, Burgers' equation, and the KPP problem (non-convex flux) show that the limiter effectively removes spurious oscillations and enforces the maximum principle and entropy inequalities.
• Convergence: The iterative limiter with the acceleration factor $\beta=2$ is shown to converge to the HO solution significantly faster than the standard approach ($\beta=1$).

Significance and Claims
The authors claim that this work provides a fairly general limiter capable of handling any hyperbolic conservation law on general meshes with a wide range of time integration schemes. The primary significance lies in bridging the gap between the robustness of low-order implicit schemes and the accuracy of high-order methods. By ensuring that the numerical scheme preserves invariant domains at the discrete level, the method enables the use of large time steps in implicit simulations without the risk of producing non-physical solutions (e.g., negative density). The paper emphasizes that the limiter maintains conservation and can be applied iteratively to refine the solution, offering a practical tool for robust, high-order simulations of complex fluid dynamics problems. The work concludes by noting that future extensions will consider linearized backward-Euler time stepping for steady-state problems.