Structure-preserving nodal DG method for Euler equations with gravity II: general equilibrium states

This paper presents a novel entropy-stable nodal discontinuous Galerkin scheme for the Euler equations with gravity that achieves well-balancing for general hydrostatic and moving equilibrium states through a linear entropy correction to the source term, while maintaining compatibility with positivity-preserving limiters and demonstrating robustness in numerical experiments.

The Gist

Imagine you are trying to simulate the weather on a planet, or the swirling gas around a black hole. You are using a super-computer to solve a complex set of rules (the Euler equations) that describe how air and gas move, heat up, and push against each other.

But there's a catch: Gravity.

Gravity is a constant, invisible hand pulling everything down. In the real world, the atmosphere doesn't just collapse; it finds a perfect balance where the downward pull of gravity is exactly matched by the upward push of air pressure. This is called equilibrium.

The Problem: The "Jittery" Simulation

For decades, computer scientists have struggled to simulate this. When they tried to model these balanced states, their computers would get "jittery."

Think of it like trying to balance a broomstick on your finger. If your hand moves even a tiny bit (due to computer rounding errors), the broomstick falls. In a simulation, the computer thinks the air is slightly out of balance, so it starts creating fake wind, fake heat, or fake pressure waves that don't exist in reality. Eventually, the numbers get so crazy (negative pressure, infinite speed) that the simulation crashes (blows up).

Previous methods could handle a "still" atmosphere (like a calm lake), but they failed miserably when the air was moving (like a jet stream or a rotating star).

The Solution: The "Triple-Threat" Algorithm

The authors of this paper (Liu, Guo, Jiang, and Zhang) have built a new, super-smart computer program (a Nodal Discontinuous Galerkin method) that acts like a master juggler. It keeps three balls in the air simultaneously, without dropping any:

1. The "Well-Balanced" Ball (Stability):

   • The Metaphor: Imagine a perfectly still pond. If you drop a pebble, the ripples should be tiny and real. If you don't drop a pebble, the water should remain perfectly still.
   • The Innovation: Their method is so precise that if the atmosphere is in a perfect balance (even if it's moving!), the computer sees it as "perfectly still" and doesn't invent fake ripples. It preserves the equilibrium state down to the tiniest decimal point. This works for both calm air and fast-moving winds.

2. The "Entropy" Ball (Thermodynamics):

   • The Metaphor: Think of entropy as "messiness" or "heat." In the real world, you can't un-mix a cup of coffee and milk, and you can't create energy out of nothing. A simulation must respect these laws, or it creates "ghost energy" that makes the numbers explode.
   • The Innovation: They added a special "correction term" (like a thermostat) that constantly checks the math. If the simulation tries to create fake energy, this term snatches it away, ensuring the laws of physics are never broken.

3. The "Positivity" Ball (Survival):

   • The Metaphor: You cannot have negative air density (you can't have "less than nothing" air) or negative pressure. If a simulation calculates a negative number, it's like a car driving backward into a wall—it crashes immediately.
   • The Innovation: They built a "safety net" (a limiter). If the math starts to drift toward a negative number, the safety net gently pulls the numbers back to a safe, positive range before the simulation crashes.

How It Works (The Secret Sauce)

The authors realized that previous methods tried to fix these problems one by one, but fixing one often broke another.

• Old Way: "Let's fix the balance!" -> Oops, now the energy is wrong.
• Old Way: "Let's fix the energy!" -> Oops, now the air density is negative.

Their New Way: They rewrote the math so that the "gravity" part of the equation is treated exactly like the "movement" part. They used a clever trick where they subtracted the known "perfect balance" from the equation and only calculated the difference.

Think of it like this: Instead of trying to calculate the weight of a giant elephant (the whole atmosphere) and a tiny ant (a small storm) separately, they calculate the weight of the elephant perfectly, and then only calculate the tiny wiggles of the ant. This makes the math incredibly stable.

Why Does This Matter?

This isn't just about better math homework. This is crucial for:

• Astrophysics: Simulating how stars form, how gas swirls around black holes, and how galaxies evolve over billions of years.
• Atmospheric Science: Predicting weather patterns and climate change with higher accuracy, especially when looking at small changes in a massive, moving atmosphere.

The Bottom Line

The authors have created a "Swiss Army Knife" for fluid simulations. It's the first tool that can handle moving atmospheres, keep the laws of physics intact, and ensure the numbers never crash, all at the same time. It's a massive leap forward for anyone trying to model the universe on a computer.

Technical Summary

Here is a detailed technical summary of the paper "Structure-preserving nodal DG method for Euler equations with gravity II: general equilibrium states" by Liu et al.

1. Problem Statement

The paper addresses the numerical simulation of the Euler equations with gravitational source terms, a system of nonlinear balance laws critical in astrophysics and atmospheric science. The primary challenge lies in designing numerical schemes that simultaneously preserve three fundamental physical properties at the discrete level:

1. Well-Balancedness (WB): The ability to exactly preserve equilibrium states (both hydrostatic and moving) without generating spurious numerical waves.
2. Entropy Stability (ES): Ensuring the discrete solution satisfies the entropy inequality, preventing non-physical oscillations and ensuring stability.
3. Positivity-Preserving (PP): Guaranteeing that physical variables (density $\rho$ and pressure $p$) remain positive throughout the simulation to avoid breakdown.

Specific Gap: While previous work (including the authors' prior study [26]) achieved these properties for hydrostatic (zero-velocity) equilibria, extending these methods to general moving equilibrium states (where velocity $u \neq 0$) is significantly more difficult. In moving equilibria, the orthogonality between entropy variables and source terms (which holds in static cases) is lost, making it hard to maintain entropy stability while preserving the equilibrium.

2. Methodology

The authors propose a structure-preserving nodal Discontinuous Galerkin (DG) scheme for both 1D and 2D cases. The method is built upon the following core components:

A. Reformulation of the Source Term

To achieve well-balancedness for a known general equilibrium state $U_e$, the governing equation is rewritten by subtracting the equilibrium balance:
$$U_t + \nabla \cdot F(U) = S(U, x) + \nabla \cdot F(U_e) - S(U_e, x)$$
This allows the discretization of the source term to exactly cancel the flux divergence when the solution is at equilibrium ($U = U_e$).

B. Entropy Conservative Fluxes

The scheme utilizes entropy conservative numerical fluxes ($F^S$) for the interior flux terms. This ensures that the flux contribution to entropy production is zero in the absence of dissipation.

C. Novel Entropy Correction Term

This is the paper's key innovation for handling moving equilibria.

• The Issue: In moving equilibria, the source term $S(U_e, x)$ is not orthogonal to the entropy variable $V$, potentially causing entropy production even when the solution is at equilibrium.
• The Solution: The authors introduce a local, high-order entropy correction term ($S_{corr}$) into the semi-discrete formulation:
  $$S_{corr} = \sigma (V - \bar{V})$$
  where $\sigma$ is a scalar coefficient calculated to balance the entropy production from the source term discretization. This term is designed to be conservative (zero cell average) and does not affect the order of accuracy or mass conservation.

D. Positivity-Preserving Limiter

The scheme incorporates a scaling-based limiter (based on Zhang et al.) applied at each stage of the time integration (SSP-RK). This limiter scales the nodal values towards the cell average to ensure $\rho > 0$ and $p > 0$ without destroying the high-order accuracy or the well-balanced property.

3. Key Contributions

1. Simultaneous Achievement of WB, ES, and PP: This is the first DG scheme in the literature to simultaneously achieve well-balancedness for general moving equilibrium states, entropy stability, and positivity preservation.
2. Generalization to Moving Equilibria: Unlike previous methods restricted to hydrostatic states, this method handles moving equilibria (e.g., rotating disks, supersonic flows) by carefully modifying the source term discretization and adding the entropy correction.
3. Rigorous Theoretical Analysis: The authors provide proofs for:
   • Mass Conservation: The scheme conserves mass exactly.
   • High-Order Accuracy: The scheme maintains optimal convergence rates ($k+1$) for smooth solutions.
   • Well-Balancedness: The scheme preserves arbitrary known equilibrium states (isothermal, isentropic, moving) up to machine precision.
   • Entropy Stability: The semi-discrete scheme is entropy conservative within an element and entropy stable globally.
   • Positivity: Under specific CFL conditions, the fully discrete scheme preserves positivity.
4. Compatibility: The entropy correction term is proven to be compatible with the PP limiter, ensuring that the enforcement of physical bounds does not violate the entropy or well-balanced properties.

4. Numerical Results

Extensive numerical experiments in 1D and 2D validate the theoretical claims:

• Well-Balancedness Tests:
  • 1D: Tested on hydrostatic, subsonic, and supersonic moving equilibria. The proposed scheme ("WBESPP") maintained errors at the machine precision level ($\sim 10^{-15}$), whereas non-WB schemes failed to preserve the equilibrium.
  • 2D: Tested on a Keplerian disk (rotating equilibrium) with both smooth and discontinuous density profiles. The scheme preserved the equilibrium perfectly, while non-WB schemes introduced large errors.
• Accuracy Tests:
  • Smooth flow tests in 1D and 2D confirmed optimal $(k+1)$-th order convergence rates for density.
• Robustness and Stability:
  • Sod Shock Tube & Double Rarefaction: The scheme successfully captured shocks and rarefaction waves. Non-ES schemes blew up immediately, and non-PP schemes failed in low-density regimes.
  • Kelvin-Helmholtz Instability: In a 2D astrophysical simulation, the scheme accurately resolved small-scale vortices while maintaining the underlying moving equilibrium.
• Time Step Analysis: The study showed that the time step restriction imposed by the source term (for positivity) is significantly weaker than the CFL restriction imposed by the flux term, making the scheme efficient.

5. Significance

This work represents a major advancement in computational fluid dynamics for astrophysical and atmospheric modeling.

• Unified Framework: It unifies the treatment of multiple physical constraints (balance, entropy, positivity) without sacrificing accuracy.
• Long-Time Simulations: The ability to preserve moving equilibria without spurious noise is crucial for long-time simulations where small errors can accumulate and destroy the physical structure of the flow (e.g., star formation, atmospheric dynamics).
• Algorithmic Innovation: The introduction of the local entropy correction term provides a new pathway for designing structure-preserving schemes for complex balance laws where source terms and fluxes interact non-trivially.

The authors conclude that this framework is highly suitable for complex, multi-dimensional astrophysical simulations and suggest future extensions to unstructured meshes and other equations of state.