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

Abstract

We develop arbitrarily high-order, stationarity-preserving stabilized finite element methods for multidimensional nonlinear hyperbolic balance laws on Cartesian grids. We aim at approximating all the steady states of the problem at hand, including non-trivial genuinely multidimensional equilibria, with a level of accuracy higher than the nominal one of the underlying scheme. We formalize more precisely the meaning of stationarity preservation, providing some technical conditions for its realizability. We then recast the multidimensional global-flux quadrature of Barsukow et al. (Num. Meth. PDEs, 2025) as a local preprocessing of the physical fluxes that maps continuous polynomial vector fields to a local space with Raviart--Thomas-type structure. Both the Galerkin and SUPG formulations are recast in this setting. The resulting methods extend the stationarity-preserving finite-volume approach of Barsukow et al. (J. Comput. Phys., 2026) to high-order continuous finite elements and Barsukow et al. (Num. Meth. PDEs, 2025) to nonlinear balance laws. We analyze key properties of the proposed schemes, including local conservation and nodal superconvergence of the discrete steady kernel, and we discuss their relation to low-Mach-compliant discretizations. We apply the framework to the compressible Euler equations with gravity. A simple source-term reformulation yields machine-precision preservation of isothermal hydrostatic equilibria. Extensive numerical benchmarks, including moving equilibrium, near-equilibrium, and instability-dominated regimes, demonstrate clear improvements in robustness and accuracy over standard SUPG and reference finite-volume methods.

Figures (12)

• Figure 1: Euler equations: $L^2$ error on momentum and density for the moving isentropic vortex. Comparison between the SUPG-Std and SUPG-GFQ methods. Left: $\rho v$. Right: $\rho$.
• Figure 2: Isentropic steady vortex with Mach = $0.7$ for the Euler equations at $t_{\mathrm{end}} = 50\,\mathrm{s}$. Top row: results on a $15 \times 15$ mesh; bottom row: results on a $35 \times 35$ mesh, both with $\mathbb{Q}^1$ elements. Left: standard SUPG. Middle: SUPG–GFQ. Right: exact solution.
• Figure 3: Convergence of the relative error on $\rho v$ for $\mathbb Q^1$ SUPG-Std and SUPG-GFQ and for FV-HLLC (left) and $\mathbb Q^2$ and $\mathbb Q^3$ (right) at time 50 for the steady isentropic vortex test case with Mach numbers $M = 10^{-2}$ (top), $M = 10^{-4}$ (center) and $M = 10^{-6}$ (bottom).
• Figure 4: Low-Mach-number stationary isentropic vortex. Isocontours of the velocity magnitude at the long-time limit ($t_{\mathrm{end}} = 2000\,\mathrm{s}$), obtained with the second-order HLLC scheme (left), SUPG-Std with $Q_1$ elements (middle), and SUPG-GFQ (right). Top: $\mathrm{Ma} = 10^{-2}$. Middle: $\mathrm{Ma} = 10^{-4}$. Bottom: $\mathrm{Ma} = 10^{-6}$. Computations performed on a $25 \times 25$ mesh.
• Figure 5: Kelvin–Helmholtz instability. Density field at final time computed with $\mathbb{Q}^1$ elements on successively refined meshes. Left column: second-order HLLC scheme. Middle column: standard SUPG method. Right column: SUPG-GFQ method.

Theorems & Definitions (18)

• Definition 1: Stationarity preserving (SP) discretization
• Remark 3: Dual residuals implication for stationarity: necessary condition for SP
• Remark 4: SP and boundary conditions
• Lemma 5: Local and global second derivative operators
• proof
• Proposition 6: Second and first derivative kernel equivalence conditions
• proof
• Proposition 7: Flux potentials and LobattoIIIA integration
• Remark 8
• Proposition 9: Objectivity