🔢 mathematics

Arbitrary order stationarity preserving stabilized finite elements for multidimensional nonlinear hyperbolic problems. Application to the Euler equations with gravity

This paper develops arbitrarily high-order, stationarity-preserving stabilized finite element methods for multidimensional nonlinear hyperbolic balance laws, which reformulate global-flux quadrature as a local preprocessing step to achieve machine-precision accuracy for steady states and demonstrate superior robustness and accuracy over standard schemes when applied to the Euler equations with gravity.

Original authors: Moussa Ziggaf, Davide Torlo, Mario Ricchiuto

Published 2026-03-25

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Moussa Ziggaf, Davide Torlo, Mario Ricchiuto

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to simulate the weather, the flow of blood in an artery, or the movement of air around an airplane wing using a computer. These are all examples of fluid dynamics. To do this, scientists use complex math equations (the Euler equations) that describe how things like density, speed, and pressure change over time.

However, computers aren't perfect. They chop the world up into a grid of tiny squares (like a pixelated image) to do the math. When you try to simulate a situation where the fluid is perfectly still or in a delicate balance (like a calm atmosphere with gravity), standard computer methods often fail. They introduce tiny, fake errors—like a camera sensor adding "noise" to a dark photo. Over time, these tiny errors grow, making a calm breeze look like a violent storm, or causing a floating balloon to suddenly sink or fly away on its own.

This paper introduces a new, smarter way to do these calculations. Here is the breakdown using simple analogies:

1. The Problem: The "Noisy Microphone"

Think of a standard computer simulation as a microphone recording a whisper in a noisy room. Even if the room is perfectly quiet, the microphone picks up static.

The Scenario: Imagine a pot of water sitting perfectly still on a stove. Gravity pulls it down, but the pressure pushes it up. They balance perfectly.
The Failure: Standard methods (like the "SUPG" method mentioned in the paper) are like that noisy microphone. They can't tell the difference between a true balance and a fake one caused by their own internal "static." They might accidentally create a tiny wave in the water that shouldn't be there. If you run the simulation for a long time, that tiny wave grows into a tsunami, ruining the result.

2. The Solution: The "Stationarity-Preserving" Method

The authors developed a new method called Stationarity-Preserving (SP) Finite Elements.

The Analogy: Imagine you are trying to balance a stack of books. A standard method tries to balance them by guessing the weight of each book and adjusting. If your guess is slightly off, the stack wobbles.
The New Method: This new method is like having a magical scale that knows exactly what a balanced stack looks like. It doesn't just guess; it forces the math to respect the rules of balance. If the real world says "the water is still," the computer math says, "Okay, the water is definitely still," and it refuses to add any fake noise.

3. How It Works: The "Global Flux" Trick

The secret sauce of this new method is something called Global Flux Quadrature.

The Old Way: Imagine you are measuring the flow of traffic in a city. You stand at one intersection, count cars, then move to the next, count again, and add them up. If you miss one car at the first intersection, your total is wrong.
The New Way: Instead of counting cars at individual intersections, the new method looks at the entire route at once. It calculates the "potential" of the flow (like measuring the total volume of water in a river from source to sea) before it even starts counting.
The Metaphor: Think of it like a Ravioli Maker. Standard methods try to stuff the dough and fill the pasta separately, which often leads to tears or leaks. This new method pre-shapes the dough (the math) specifically to hold the filling (the physics) perfectly, ensuring no leaks (errors) happen, even when the filling is complex.

4. Why It Matters: The "Low Mach" Challenge

The paper specifically tests this on low-speed flows (called "Low Mach" flows), like wind moving gently or air in a room.

The Challenge: In slow flows, the "noise" from standard methods is often louder than the actual wind. It's like trying to hear a whisper while someone is shouting.
The Result: The new method is so quiet and precise that it can hear the whisper. The authors tested it on:
- Moving Vortices: Swirling winds that should keep their shape. The new method kept them perfect; the old ones made them dissolve.
- Hydrostatic Equilibrium: Air sitting still under gravity. The new method kept it perfectly still for hours; the old ones made it wiggle and break apart.
- Instabilities: Things like the "Kelvin-Helmholtz" instability (think of the ripples you see when oil floats on water). The new method captured the beautiful, swirling patterns clearly, while the old methods blurred them into a mess.

5. The "Gravity" Bonus

The authors also added a special trick to handle gravity perfectly.

The Analogy: Imagine a scale that is perfectly balanced with a heavy rock on one side. If you add a tiny feather, the scale should tip slightly. Standard methods often tip the scale too much because they miscalculate the weight of the rock.
The Fix: This new method has a "calibration mode" for gravity. It knows exactly how to balance the weight of the air against the pull of gravity. If the air is in a perfect "hydrostatic" state (like the atmosphere on a calm day), the computer keeps it exactly that way, down to the last decimal point (machine precision).

Summary

This paper presents a super-accurate, high-order calculator for fluid dynamics.

Old Method: Good for big, fast, chaotic storms, but terrible at simulating calm, balanced states. It adds fake noise that ruins long-term predictions.
New Method: It is "stationarity-preserving," meaning it respects the natural balance of the universe. It can simulate calm air, gentle winds, and delicate instabilities with incredible precision, without adding fake noise.

It's like upgrading from a standard camera that adds grain to every photo, to a high-end scientific camera that can capture a perfectly still pond without a single ripple of digital noise. This allows scientists to study complex weather patterns, atmospheric physics, and engineering designs with a level of trust and accuracy that was previously impossible.

1. Problem Statement

The paper addresses the numerical simulation of multidimensional nonlinear hyperbolic balance laws, specifically the compressible Euler equations with gravitational source terms. A critical challenge in this domain is the accurate preservation of stationary states (equilibria), such as hydrostatic atmospheres, isentropic vortices, and moving equilibria.

Standard numerical methods, including high-order Finite Element (FE) and Finite Volume (FV) schemes, often fail to preserve these equilibria exactly due to:

Numerical Dissipation: Standard stabilization techniques (like Streamline Upwind/Petrov-Galerkin, SUPG) introduce dissipation that destroys steady states unless the solution is constant.
Incompatibility with Source Terms: The discrete approximation of the flux divergence ( $\nabla \cdot F$ ) and the source term ( $S$ ) often do not cancel each other out exactly at the discrete level, leading to spurious waves and loss of accuracy, particularly in low Mach number regimes.
Lack of Multidimensional Capability: Many "well-balanced" schemes are restricted to 1D or specific 2D configurations (e.g., hydrostatic only) and cannot handle genuinely multidimensional moving equilibria (e.g., rotating vortices).

2. Methodology

The authors propose a high-order, stationarity-preserving (SP) stabilized finite element method based on a Global Flux Quadrature (GFQ) approach.

A. Theoretical Framework: Stationarity Preservation

The paper formalizes the definition of Stationarity Preservation (SP). A scheme is SP if it admits a rich set of discrete steady states that converge to the continuous equilibria as the mesh size $h \to 0$ .

Condition 2 (Dual Residuals): The authors establish a necessary condition for SP: there must exist a set of local discrete residuals (defined on a dual mesh or sub-cells) such that if these residuals vanish, the global evolution operator vanishes.
Failure of Standard SUPG: It is proven that standard SUPG methods fail this condition because their numerical dissipation (based on Laplacian-like structures) does not vanish for non-constant steady states, and the discrete kernel of the divergence operator is too restrictive on standard continuous Lagrange elements.

B. The Global Flux Quadrature (GFQ) Approach

To overcome these limitations, the authors adapt the GFQ concept (previously used for linear acoustics and shallow water) to nonlinear systems:

Flux Pre-processing: Instead of directly discretizing $\nabla \cdot F$ , the method introduces flux potentials ( $F_1, F_2$ ) such that $F = (\partial_{x_2} F_1, \partial_{x_1} F_2)$ .
Integration Strategy: The potentials are computed by integrating the physical fluxes along coordinate lines using LobattoIIIA collocation (high-order quadrature). This maps the continuous polynomial vector fields into a local space with a Raviart-Thomas (RT) structure, ensuring the divergence-conforming property.
Source Term Treatment: The source term $S$ is similarly integrated to define a potential $S_h$ , allowing the balance law to be written as $\partial_{x_1}\partial_{x_2}(F_1 + F_2 + S) = 0$ .
Stabilization: The SUPG stabilization term is reformulated using these potentials. The stabilization now acts on the integrated residuals rather than the raw fluxes, ensuring that the dissipation vanishes exactly when the discrete stationarity conditions are met.
Time Integration: The method uses a Defect Correction (DeC) strategy to achieve explicit, high-order time stepping without inverting the non-diagonal mass matrix introduced by the stabilization.

C. Well-Balanced Correction for Gravity

For the Euler equations with gravity, a specific modification is introduced to preserve isothermal hydrostatic equilibria exactly (up to machine precision). This involves reformulating the source term using reference states ( $\bar{\rho}, \bar{p}$ ) and integrating the potential gradient in a way that matches the discrete pressure gradient.

3. Key Contributions

Arbitrary Order Extension: The first extension of stationarity-preserving global flux quadrature to arbitrary high-order continuous finite elements for nonlinear hyperbolic balance laws.
Formal Definition & Conditions: A rigorous definition of SP and the derivation of necessary conditions (Dual Residuals) to verify if a scheme is SP.
Raviart-Thomas Connection: Demonstrating that the GFQ pre-processing effectively maps standard $Q_K$ elements to a space with Raviart-Thomas properties, enabling the preservation of solenoidal/divergence-free structures without changing the basis functions.
Super-convergence: Proof that the method achieves super-convergence (order $K+2$ or higher) at steady states for polynomial degrees $K \ge 2$ , a property inherited from the LobattoIIIA collocation method.
Low-Mach Compliance: The method is shown to be inherently compatible with low Mach number flows, avoiding the numerical stiffness and dissipation issues common in standard upwind schemes.

4. Numerical Results

The authors validate the method through extensive benchmarks on the Euler equations with gravity:

Moving Isentropic Vortex: The GFQ method shows lower error magnitudes (1.3x to 2x improvement) compared to standard SUPG, even for non-stationary moving solutions.
Steady Isentropic Vortex:
- Standard SUPG fails to preserve the vortex structure, introducing spurious oscillations and dissipation.
- SUPG-GFQ preserves the vortex perfectly.
- Super-convergence: While standard SUPG converges at order $K+1$ , SUPG-GFQ achieves order $K+2$ (or higher) for $K \ge 2$ .
Low Mach Number Regimes ( $Ma = 10^{-2}$ to $10^{-6}$ ):
- SUPG-GFQ achieves machine-precision accuracy for density and momentum on coarse meshes where standard methods fail.
- Long-time simulations ( $t=2000$ s) show SUPG-GFQ preserves the initial structure, whereas standard methods and FV-HLLC suffer from excessive numerical diffusion.
Hydrodynamic Instabilities (Kelvin-Helmholtz & Rayleigh-Taylor):
- The method resolves fine-scale vortical structures with minimal mesh imprinting.
- Standard methods show "glitches" and mesh-aligned artifacts; SUPG-GFQ captures the physics accurately even on coarse grids.
Hydrostatic Equilibria:
- The "Well-Balanced" (WB) variant preserves isothermal hydrostatic equilibria to machine precision ( $10^{-15}$ ).
- Small perturbations on these equilibria are captured without generating spurious waves, outperforming both standard SUPG and non-WB GFQ.
Shock-Vortex Interaction: The method remains stable and accurate for mild shocks ( $M_s=1.1$ ), capturing the interaction with high-order smoothness and reduced mesh imprinting.

5. Significance

This work represents a significant advancement in computational fluid dynamics (CFD) for several reasons:

Unified High-Order Accuracy: It bridges the gap between high-order accuracy and the preservation of complex physical equilibria, a problem that has historically forced a trade-off between the two.
Robustness in Low Mach Regimes: It provides a robust framework for atmospheric and astrophysical flows where low Mach numbers and gravity are dominant, eliminating the need for ad-hoc preconditioning or low-Mach specific fixes.
Generalizability: The methodology is not limited to the Euler equations; the framework of "flux pre-processing" and "dual residuals" can be applied to other balance laws (e.g., magnetohydrodynamics, shallow water, multiphase flows).
Efficiency: By using continuous finite elements and DeC time stepping, the method avoids the complexity of discontinuous Galerkin (DG) or finite volume reconstruction while achieving comparable or superior accuracy for steady and near-steady problems.

In conclusion, the proposed SUPG-GFQ method offers a superior alternative to standard stabilized finite elements for problems involving equilibria, low Mach numbers, and complex multidimensional dynamics, achieving machine-precision preservation of steady states and super-convergent accuracy.