A dual-pairing summation-by-parts finite difference framework for nonlinear conservation laws

1. Problem Statement

The numerical simulation of nonlinear hyperbolic conservation laws (e.g., Euler equations, shallow water equations) presents a fundamental challenge: balancing high-order accuracy with robustness.

• Low-order methods (e.g., finite volume with Riemann solvers) are robust but suffer from excessive numerical dissipation and dispersion, failing to resolve fine-scale features like turbulence and vortices.
• High-order methods (e.g., spectral element, high-order finite difference) offer high accuracy and low dissipation but often lack robustness. They are prone to instability, spurious oscillations, and simulation crashes when encountering shocks, discontinuities, or under-resolved nonlinear modes.
• Entropy Stability: While linear stability is often proven, it is insufficient for nonlinear problems. Entropy stability (ensuring the discrete solution satisfies a discrete entropy inequality) is required to guarantee convergence to the correct weak solution and physical robustness.
• The Gap: Existing high-order entropy-stable methods (often based on Summation-by-Parts or SBP operators) typically rely on central difference stencils. These induce dissipation only at element interfaces via numerical fluxes. They often lack local linear energy stability within the element volume, making them susceptible to instabilities in under-resolved regions (e.g., turbulence) without ad-hoc filtering or limiting.

2. Methodology

The authors propose a novel framework combining Dual-Pairing (DP) operators, Upwind Summation-by-Parts (SBP) finite difference (FD) and discontinuous Galerkin (DG) methods, and skew-symmetric flux splitting.

A. Dual-Pairing (DP) SBP Operators

Unlike traditional SBP operators that use a single central derivative operator $D$, the DP framework employs a pair of forward ($D_+$) and backward ($D_-$) derivative operators.

• SBP Property: Together, they satisfy a discrete integration-by-parts property: $\langle D_+ f, g \rangle_H + \langle f, D_- g \rangle_H = \text{boundary terms}$.
• Upwind Property: The difference operator $(D_+ - D_-)$ is designed to be negative semi-definite. This introduces built-in volume dissipation (upwinding) directly inside the elements, rather than relying solely on interface fluxes.
• Flexibility: This dual structure provides extra degrees of freedom to design operators that are both high-order accurate and inherently dissipative for unresolved modes.

B. Skew-Symmetric Reformulation

To achieve entropy stability without relying on the chain rule (which is difficult to replicate discretely), the conservation laws are rewritten in a skew-symmetric form.

• The flux divergence $\nabla \cdot f(u)$ is split into a convex combination of the flux form and the advective form.
• This allows the energy method (entropy analysis) to rely solely on integration-by-parts, ensuring that the discrete scheme mimics the continuous entropy conservation properties.

C. The Semi-Discrete Scheme

The proposed semi-discrete equation for a system of conservation laws is:
$$\frac{d}{dt}u + F(u, f(u), D_x) = \underbrace{\frac{1}{2}H^{-1}(\alpha \otimes (B_I + B_n))g}_{\text{Surface Terms (Interface Flux)}} + \underbrace{\frac{1}{2}\Gamma (I \otimes (D_+ - D_-))g}_{\text{Volume Terms (Upwind Dissipation)}}$$

• $F(u, f(u), D_x)$: The skew-symmetric discrete flux operator.
• Surface Terms: Handle interface coupling and boundary conditions using upwind numerical fluxes.
• Volume Terms: The term involving $(D_+ - D_-)$ acts as an intrinsic high-order filter. It dissipates entropy locally where the solution is non-smooth or under-resolved, acting as a "shock detector" and stabilizer without explicit limiting.

3. Key Contributions

1. Novel Framework: Introduction of a systematic framework combining DP SBP FD/DG operators with skew-symmetric and upwind flux splitting. This is the first reported method of this specific combination for nonlinear systems.
2. Intrinsic Volume Dissipation: Unlike traditional SBP/DGSEM methods that dissipate only at interfaces, this method induces dissipation within the elements via the upwind DP operators. This provides local linear energy stability, crucial for robustness in turbulent flows.
3. Provable Entropy Stability: The authors prove that the semi-discrete scheme is entropy-stable (entropy non-increasing) and conservative for arbitrary nonlinear hyperbolic conservation laws, provided the DP operators satisfy specific assumptions.
4. Automatic Filtering: The volume upwind term acts as an intrinsic, high-order accurate filter. It automatically detects and resolves regions of poor resolution (discontinuities, shocks) without requiring artificial viscosity, limiters, or ad-hoc stabilization strategies.
5. Unified FD/DG Approach: The framework is applicable to both Finite Difference (on uniform grids) and Discontinuous Galerkin (on Legendre-Gauss-Lobatto nodes) formulations.

4. Numerical Results

The framework was validated against standard benchmarks using the Julia programming language, comparing the proposed DP DG/FD (with volume upwinding, $\Gamma > 0$) against:

• Linearly stable schemes (Lax-Friedrichs splitting).
• Standard entropy-stable SBP/DGSEM (without volume upwinding, $\Gamma = 0$).

Key Findings:

• Inviscid Burgers' Equation:
  • Standard entropy-stable schemes ($\Gamma=0$) and linearly stable schemes crashed or produced severe oscillations when shocks formed.
  • The proposed DP scheme ($\Gamma>0$) remained stable, suppressed oscillations, and maintained conservation.
• Compressible Euler Equations (Sod Shock Tube):
  • All other schemes crashed before the final simulation time for various mesh resolutions and polynomial orders.
  • The DP scheme completed all simulations successfully, demonstrating superior robustness.
• Kelvin-Helmholtz Instability (2D Turbulence):
  • Tested with high Atwood numbers (strong density contrast).
  • Standard schemes failed to control enstrophy growth and crashed due to numerical noise.
  • The DP scheme successfully resolved the turbulent mixing layers, controlled enstrophy, and ran to the final time for all configurations, including extreme density ratios ($A=0.99$).
• Accuracy:
  • Convergence rates for smooth solutions (e.g., Isentropic Vortex) matched or slightly exceeded theoretical expectations.
  • The method preserved "lake at rest" steady states to machine precision (well-balanced property).

5. Significance

This work addresses a critical bottleneck in computational fluid dynamics: the trade-off between accuracy and stability in high-order methods.

• Robustness without Limiters: By embedding dissipation directly into the derivative operators, the method achieves robustness comparable to low-order schemes while retaining high-order accuracy, eliminating the need for complex, non-linear limiters that degrade accuracy near smooth regions.
• Theoretical Guarantee: The provable entropy stability ensures that the numerical solutions converge to the physically correct weak solutions, a requirement often missing in high-order methods.
• Scalability: The framework is designed for modern hardware architectures and is applicable to complex geometries (via multi-block extensions) and various physical systems (Euler, Shallow Water, Magnetohydrodynamics).

In summary, the paper presents a mathematically rigorous and numerically superior approach for simulating nonlinear conservation laws, offering a "best of both worlds" solution that is simultaneously high-order accurate, entropy-stable, and robust against shocks and turbulence.