Computing Radially-Symmetric Solutions of the Ultra-Relativistic Euler Equations with Entropy-Stable Discontinuous Galerkin Methods

This paper derives an entropy-stable flux for the ultra-relativistic Euler equations by computing the main field and potentials, and validates the resulting Discontinuous Galerkin method through 2D and 3D simulations of radially symmetric problems involving shock waves and pressure blow-up.

The Gist

Imagine a universe where the usual rules of "heavy" matter don't apply. Instead, we are dealing with a gas so hot and energetic that its heat (thermal energy) is the main character, completely overshadowing the mass of the particles. This is the world of ultra-relativistic fluids. Think of it like a crowd of people running so fast that their speed and the energy of their movement matter far more than their actual body weight.

This paper is about building a better, safer, and more accurate "calculator" (a computer simulation) to predict how this super-hot gas behaves, especially when it moves in a circle or sphere (like an explosion or a bubble).

Here is a breakdown of what the authors did, using simple analogies:

1. The Problem: Predicting the Unpredictable

When these super-fast gases move, they can do wild things. They can suddenly form shock waves (like a sonic boom from a jet, but in a fluid) or experience a pressure blow-up (where the pressure gets infinitely high in a tiny spot, like a balloon popping but with extreme force).

Previous computer programs could simulate this, but they were like a shaky camera: sometimes they got the big picture right, but they might miss the tiny, dangerous details or even crash when things got too chaotic. The authors wanted to build a camera that never shakes and never crashes, even when the gas is doing its wildest dance.

2. The Solution: The "Entropy-Stable" Rulebook

The authors created a new set of rules for their computer program called an "entropy-stable" method.

• The Analogy: Imagine you are trying to keep a messy room tidy. "Entropy" is a measure of how messy the room is. In physics, nature generally wants to get messier over time (like a room getting messy if you don't clean it).
• The Innovation: The authors designed a specific "flux" (a way of calculating how the gas moves from one spot to another) that respects this rule of messiness. They proved mathematically that their new rulebook ensures the simulation never gets "too clean" or "too messy" in a way that breaks physics. It keeps the simulation stable, preventing the computer from crashing when the gas gets violent.

They derived this new rulebook from scratch, creating a "two-point flux" (a way to calculate the flow between two neighboring points) that acts like a perfectly balanced scale.

3. The Tools: A High-Definition Camera (DG Methods)

To run these simulations, they used a technique called Discontinuous Galerkin (DG) methods.

• The Analogy: Imagine trying to draw a picture of a stormy ocean. A low-resolution map might just show a blue blob. A high-resolution map breaks the ocean into millions of tiny tiles.
• How it works: Their method breaks the space into tiny 3D blocks (like LEGO bricks). Inside each block, they use complex math to describe the gas. They also use a trick called "flux differencing," which is like checking the math between every single pair of neighboring bricks to ensure the energy balance is perfect.

4. The Safety Net: Shock Capturing

Even with a perfect rulebook, some things happen so fast (like a shock wave hitting a wall) that the computer needs a safety net.

• The Analogy: Think of a high-speed race car. On a smooth track, it uses a high-performance engine (the complex math). But if it hits a bump, it switches to a sturdy, slower suspension (a simpler, robust method) so it doesn't flip over.
• The Implementation: Their program automatically detects when the gas is getting too chaotic (a "shock") and temporarily switches to a simpler, sturdier calculation method just for that tiny area, then switches back to the high-performance math once the danger passes.

5. The Test Drive: 2D vs. 3D

The authors tested their new calculator on five different scenarios, comparing their new 3D simulation against a trusted 1D "radial" solver (a specialized tool that only looks at the center of the explosion).

• The Scenarios: They simulated things like:
  • A shock wave moving through gas.
  • A bubble expanding into a vacuum.
  • A bubble collapsing (imploding).
  • Waves moving in a sine pattern.
• The Results:
  • In 2D (Flat): The new calculator matched the trusted tool perfectly. It captured the shock waves and pressure spikes exactly as expected.
  • In 3D (Real World): This is the big achievement. They are the first to show these results in full 3D. However, they noted a limitation: 3D is incredibly expensive to compute. While the 2D simulation could see a pressure spike of nearly 300, the 3D simulation (running on a standard computer) only saw a spike of about 289.
  • The Takeaway: The 3D results were still excellent and matched the 2D trends, but the extreme "peaks" of pressure were slightly smoothed out because the computer needed to use a slightly coarser grid to finish the job in a reasonable amount of time.

Summary

The authors built a super-stable, high-definition computer simulator for ultra-hot, super-fast gases. They created a new mathematical "rulebook" that prevents the simulation from breaking when things get violent. They proved it works perfectly in 2D and successfully ran it for the first time in full 3D, showing that while 3D is harder to compute, their method captures the essential physics of shock waves and pressure explosions accurately.

They also made sure to share all their code and data, so anyone else can try to reproduce their results, ensuring the science is transparent and verifiable.

Technical Summary

Technical Summary: Computing Radially-Symmetric Solutions of the Ultra-Relativistic Euler Equations with Entropy-Stable Discontinuous Galerkin Methods

Problem Statement
The paper addresses the numerical simulation of the ultra-relativistic Euler equations, which describe ideal gases where thermal energy dominates. These equations are formulated within the framework of hyperbolic conservation laws in flat Minkowski space-time. The specific challenge lies in developing robust, high-order numerical methods capable of handling the complex wave structures inherent in these systems, particularly shock waves and pressure blow-ups, in both two and three spatial dimensions. While previous works (e.g., Kunik et al.) have established one-dimensional radially symmetric solvers and benchmark problems, there is a need for genuine multi-dimensional solvers that preserve the structural properties of the equations, specifically entropy stability, to ensure physical correctness and numerical robustness.

Methodology
The authors develop a high-order Discontinuous Galerkin (DG) method based on the flux differencing framework combined with Summation-By-Parts (SBP) operators. The core of the methodology involves several key components:

1. Entropy-Stable Flux Derivation: The authors derive a new entropy-conservative two-point flux for the ultra-relativistic Euler equations. Unlike earlier entropy-conservative fluxes based on integral averages, this derivation utilizes differential averages and specific variable transformations (involving $u_i p^{-1/4}$ and $p$) to satisfy the entropy condition $\langle \llbracket \mathbf{h} \rrbracket, \tilde{\mathbf{F}}_k \rangle - \llbracket \psi_k \rrbracket = 0$. The resulting flux is symmetric, consistent, and entropy-conservative.
2. Numerical Scheme Construction: The entropy-conservative flux is employed in the volume terms of a nodal DG spectral element method (DGSEM) using Gauss-Lobatto-Legendre nodes. To ensure robustness in the presence of discontinuities (shocks), the high-order flux differencing scheme is blended with a first-order finite volume method using a local Lax-Friedrichs/Rusanov flux. This blending is controlled by a shock indicator applied to the variable $\tilde{p}\sqrt{1+|\mathbf{u}|^2}$.
3. Implementation: The simulations are performed using the open-source Julia framework Trixi.jl. The implementation includes adaptive mesh refinement (AMR) based on $p4est$ and a positivity-preserving limiter for pressure. Time integration is handled by a third-order, four-stage explicit strong stability preserving Runge-Kutta method with adaptive step sizing.
4. Validation Strategy: The multi-dimensional results are compared against a specific one-dimensional radially symmetric solver (RadSymS) developed in previous literature, which provides reference solutions for radially symmetric problems. Additionally, for self-similar test cases, solutions are compared against reference solutions obtained by solving the corresponding ordinary differential equations (ODEs).

Key Contributions

• Derivation of Entropy-Conservative Flux: The primary theoretical contribution is the first derivation of an entropy-conservative numerical flux specifically for the ultra-relativistic Euler equations. This includes the calculation of the main field (entropy variables), entropy flux, and entropy flux potential.
• High-Order Entropy-Stable DG Scheme: The authors present a high-order entropy-stable DG scheme for these equations, extending the Tadmor framework to the ultra-relativistic regime.
• Genuine Multi-Dimensional Benchmarks: The paper provides the first genuine three-dimensional simulation results for the benchmark problems proposed by Kunik et al. (2024). These benchmarks include radially symmetric problems with shock formation and pressure blow-up.
• Reproducibility: All Julia source code and data required to reproduce the numerical results are made available in an accompanying repository.

Numerical Results
The paper presents results for five benchmark examples in 2D and 3D:

1. Shock and Stationary Part: A problem involving a shock wave and a stationary region. The multi-dimensional solver (Trixi.jl) shows excellent agreement with the 1D solver (RadSymS) and the ODE reference solution. A convergence test confirms the high-order accuracy (EOC $\approx$ 4) of the scheme for smooth solutions.
2. Self-Similar Expansion: A smooth rarefaction wave solution. Again, the 2D and 3D results align closely with the 1D reference and ODE solutions.
3. Expansion of a Spherical Bubble: A problem where high-pressure gas expands into a low-pressure region, leading to shock formation and a subsequent pressure blow-up upon reflection at the origin. The 2D results show excellent agreement between the DG method and the 1D solver. In 3D, while the qualitative behavior is captured, the peak pressure values are under-resolved compared to the 1D reference due to computational cost constraints limiting the refinement level.
4. Collapse of a Spherical Bubble: Similar to Example 3 but with an initial low-pressure core collapsing. The results mirror the findings of Example 3, with high agreement in 2D and resolution limitations affecting the peak pressure capture in 3D.
5. Initially Sine-Shaped Radial Velocity: A complex wave structure problem. The results confirm the scheme's ability to handle complex interactions, with the same trend of high fidelity in 2D and resolution-dependent discrepancies in 3D regarding peak pressures.

In all cases, the total entropy is shown to increase over time (as expected for entropy-stable schemes with dissipation), and entropy is conserved up to machine precision in smooth, non-shock scenarios.

Significance and Claims
The authors claim that this work represents the first time genuine three-dimensional results for these specific ultra-relativistic Euler benchmarks have been obtained and presented. The study demonstrates that the derived entropy-stable flux and the associated DG framework are effective for simulating ultra-relativistic flows in multiple dimensions. The paper emphasizes that while the 2D results match reference solutions with high precision, the 3D results are currently limited by computational complexity and workstation resources, preventing the resolution of extremely localized pressure peaks seen in the 1D references. The authors suggest that these benchmarks serve as challenging test cases for other multi-dimensional solvers and note that extending the work to different space-times is a potential future direction. The work is funded by the Deutsche Forschungsgemeinschaft (DFG) under specific projects related to hyperbolic balance laws.