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

TL;DR

This paper develops an entropy-stable discontinuous Galerkin method for solving ultra-relativistic Euler equations, enabling accurate simulation of shock waves and pressure blow-up in radially symmetric relativistic gases.

Contribution

It introduces an entropy-stable flux and corresponding variables for the ultra-relativistic Euler equations, advancing numerical stability and accuracy in relativistic fluid simulations.

Findings

Successful 2D and 3D simulations of shock formation

Stable solutions with pressure blow-up behavior

Validation against benchmark problems

Abstract

The ultra--relativistic Euler equations describe gases in the relativistic case when the thermal energy dominates. These equations for an ideal gas are given in terms of the pressure, the spatial part of the dimensionless four-velocity, and the particle density. Kunik et al.\ (2024, https://doi.org/10.1016/j.jcp.2024.113330) proposed genuine multi--dimensional benchmark problems for the ultra--relativistic Euler equations. In particular, they compared full two-dimensional discontinuous Galerkin simulations for radially symmetric problems with solutions computed using a specific one-dimensional scheme. Of particular interest in the solutions are the formation of shock waves and a pressure blow-up. In the present work we derive an entropy-stable flux for the ultra--relativistic Euler equations. Therefore, we derive the main field (or entropy variables) and the corresponding potentials. We then present the entropy-stable flux and conclude with simulation results for different test cases both in 2D and in 3D.

Tables

Table 1. Table 1: Discretization parameters (domain and maximum refinement level L L ) for Trixi.jl

	Ex. 1		Ex. 2		Ex. 3		Ex. 4		Ex. 5
	$d=2$	$d=3$	$d=2$	$d=3$	$d=2$	$d=3$	$d=2$	$d=3$	$d=2$	$d=3$
	${[-2,2]}^d$		${[-2,2]}^d$		${[-6,6]}^d$		${[-6,6]}^d$		${[-5,5]}^d$
$L$	11	10	14	13	13	9	13	10	14	9

Table 2. Table 2: Convergence results for the 2D version of Trixi.jl. The discrete L 2 L^{2} errors are computed at the final time for the conserved variables.

Level	$L^2$ error $w_1$	EOC $w_1$	$L^2$ error $w_2$	EOC $w_2$	$L^2$ error $w_3$	EOC $w_3$
2	$7.52\text{×}{10}^{-03}$		$7.52\text{×}{10}^{-03}$		$1.10\text{×}{10}^{-02}$
3	$4.87\text{×}{10}^{-04}$	3.95	$4.87\text{×}{10}^{-04}$	3.95	$8.51\text{×}{10}^{-04}$	3.69
4	$3.29\text{×}{10}^{-05}$	3.89	$3.29\text{×}{10}^{-05}$	3.89	$6.67\text{×}{10}^{-05}$	3.67
5	$1.93\text{×}{10}^{-06}$	4.09	$1.93\text{×}{10}^{-06}$	4.09	$3.43\text{×}{10}^{-06}$	4.28
6	$1.22\text{×}{10}^{-07}$	3.99	$1.22\text{×}{10}^{-07}$	3.99	$2.30\text{×}{10}^{-07}$	3.90

Table 3. Table 3: Convergence results for the 3D version of Trixi.jl. The discrete L 2 L^{2} errors are computed at the final time for the conserved variables.

Level	$L^2$ error $w_1$	EOC $w_1$	$L^2$ error $w_2$	EOC $w_2$	$L^2$ error $w_3$	EOC $w_3$	$L^2$ error $w_4$	EOC $w_4$
2	$6.98\text{×}{10}^{-03}$		$6.98\text{×}{10}^{-03}$		$6.98\text{×}{10}^{-03}$		$1.18\text{×}{10}^{-02}$
3	$4.78\text{×}{10}^{-04}$	3.87	$4.78\text{×}{10}^{-04}$	3.87	$4.78\text{×}{10}^{-04}$	3.87	$1.00\text{×}{10}^{-03}$	3.56
4	$3.29\text{×}{10}^{-05}$	3.86	$3.29\text{×}{10}^{-05}$	3.86	$3.29\text{×}{10}^{-05}$	3.86	$7.48\text{×}{10}^{-05}$	3.74
5	$1.78\text{×}{10}^{-06}$	4.21	$1.78\text{×}{10}^{-06}$	4.21	$1.78\text{×}{10}^{-06}$	4.21	$4.08\text{×}{10}^{-06}$	4.20
6	$1.14\text{×}{10}^{-07}$	3.96	$1.14\text{×}{10}^{-07}$	3.96	$1.14\text{×}{10}^{-07}$	3.96	$2.48\text{×}{10}^{-07}$	4.04

Table 4. Table 4: Shock states for d = 2 , 3 d=2,3 in Example 1

	$\stackrel{tilde}{s}$	$p_{-}$	$v_{-}$	$p_{+}$	$v_{+}$
$d=2$	$0.45503$	$15.75505$	$0$	$5.71869$	$-0.41629$
$d=3$	$0.52314$	$25.56463$	$0$	$17.16524$	$-0.17106$

Equations

$$g_{\alpha \beta }=\{\begin{array}{ll}+1,& \alpha =\beta =0,\\ -1,& \alpha =\beta =1,2,3,\\ 0,& \alpha \ne \beta .\end{array}$$

$$\text{sumop}{\text{slimits@}}_{\beta =0}^3\frac{\partial T_{\alpha \beta }}{\partial x_{\beta }}=0$$

$$T_{\alpha \beta }=-p{g}_{\alpha \beta }+4p{u}_{\alpha }{u}_{\beta }.$$

$$\mathbf{u}=\frac{\Gamma }{c}\mathbf{v}$$

$$\Gamma =\frac{1}{\sqrt{1-\frac{|\mathbf{v}{|}^2}{c^2}}}.$$

$$e=3p.$$

$$\frac{\partial}{\partial t}\left(3p+4p|\mathbf{u}|^{2}\right)+\sumop\slimits@\limits_{k=1}^{d}\frac{\partial}{\partial x_{k}}\left(4pu_{k}\sqrt{1+|\mathbf{u}|^{2}}\right)=0,$$

$$\frac{\partial}{\partial t}\left(4pu_{j}\sqrt{1+|\mathbf{u}|^{2}}\right)+\sumop\slimits@\limits_{k=1}^{d}\frac{\partial}{\partial x_{k}}\left(p\delta_{jk}+4pu_{j}u_{k}\right)=0,\quad j=1,\dots,d.$$

$$w_j:=4p{u}_j\sqrt{1+|\mathbf{u}{|}^2},j=1,\dots ,d,$$

$$w_{d+1}:=p(3+4|\mathbf{u}{|}^2).$$

$$u_j=\frac{w_j}{\sqrt{4p(w_{d+1}+p)}},j=1,\dots ,d,$$

$$p=\frac{1}{3}\left(-w_{d+1}+\sqrt{-3|\overline{\mathbf{w}}{|}^2+4w_{d+1}^2}\right)\ge 0.$$

$$\eta =p^{\frac{3}{4}}\sqrt{1+|\mathbf{u}{|}^2}$$

$$\mathbf{q}=p^{\frac{3}{4}}\mathbf{u}.$$

$${\nabla }_{\mathbf{w}}{q}_k^T={\nabla }_{\mathbf{w}}{\eta }^T\cdot {\mathbf{D}}_{\mathbf{w}}{\mathbf{F}}_k$$

$$h_j:={\partial }_{w_j}\eta =\{\begin{array}{ll}-\frac{1}{4}\frac{\eta }{p}\frac{u_j}{\sqrt{1+|\mathbf{u}{|}^2}},& j=1,\dots ,d,\\ \frac{1}{4}\frac{\eta }{p},& j=d+1.\end{array}$$

$${\psi }_k={\nabla }_{\mathbf{W}}\eta \cdot {\mathbf{F}}_k-q_k=-\frac{1}{4}\eta \frac{u_k}{\sqrt{1+|\mathbf{u}{|}^2}}\iff \mathbf{\psi }=-\frac{1}{4}\eta \frac{\mathbf{u}}{\sqrt{1+|\mathbf{u}{|}^2}}.$$

$$\left\llbracket\mathbf{h}\right\rrbracket\cdot\tilde{\mathbf{F}}_{k}-\left\llbracket\psi_{k}\right\rrbracket=0,$$

$${\mathbf{F}}_k=\left(\begin{array}{c}p{\delta }_{1k}+4p{u}_1{u}_k\\ \text{vbox}...\\ p{\delta }_{dk}+4p{u}_d{u}_k\\ 4p{u}_k\sqrt{1+|\mathbf{u}{|}^2}\end{array}\right)=\left(\begin{array}{c}p{\delta }_{1k}+\frac{w_1{w}_k}{w_{d+1}+p}\\ \text{vbox}...\\ p{\delta }_{dk}+\frac{w_d{w}_k}{w_{d+1}+p}\\ w_k\end{array}\right)$$

$${\partial }_t\mathbf{w}+\text{sumop}{\text{slimits@}}_{k=1}^d{\partial }_{x_k}{\mathbf{F}}_k(\mathbf{w})=0.$$

$$\text{resizebox}0.0pt0.0pt\{\{a\text{resizebox}0.0pt0.0pt\}\}=\frac{a^{+}+a^{-}}{2}.$$

$$\displaystyle\left\llbracket ab\right\rrbracket$$

$$\displaystyle\left\llbracket a^{2}\right\rrbracket$$

$$\left\llbracket\sqrt{a}\right\rrbracket=\frac{\left\llbracket a\right\rrbracket}{2\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\{\mkern-6.0mu\{$}}}\sqrt{a}\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\}\mkern-6.0mu\}$}}}}.$$

$$-(a_{-}\sqrt{a_{+}}+a_{+}\sqrt{a_{-}})\left\llbracket a^{-1/2}\right\rrbracket=-(a_{-}+\sqrt{a_{+}}\sqrt{a_{-}}-\sqrt{a_{+}}\sqrt{a_{-}}-a_{+})=\left\llbracket a\right\rrbracket,$$

$$\frac{\left\llbracket a^{-1/2}\right\rrbracket}{\left\llbracket a\right\rrbracket}=-\frac{1}{a_{-}\sqrt{a_{+}}+a_{+}\sqrt{a_{-}}}.$$

$$\left\llbracket h_{i}\right\rrbracket=-\frac{1}{4}\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$[\mkern-3.0mu[$}}}u_{i}p^{-1/4}\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$]\mkern-3.0mu]$}}},\quad i\in\{1,\dots,d\}.$$

$$\displaystyle\left\llbracket h_{d+1}\right\rrbracket$$

$$\displaystyle=\frac{1}{8\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\{\mkern-6.0mu\{$}}}\sqrt{p^{-1/2}+\left(|\mathbf{u}|p^{-1/4}\right)^{2}}\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\}\mkern-6.0mu\}$}}}}\left(\frac{\left\llbracket p^{-1/2}\right\rrbracket}{\left\llbracket p\right\rrbracket}\left\llbracket p\right\rrbracket+2\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\{\mkern-6.0mu\{$}}}u_{i}p^{-1/4}\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\}\mkern-6.0mu\}$}}}\left\llbracket u_{i}p^{-1/4}\right\rrbracket\right),$$

$$\left\llbracket\psi_{k}\right\rrbracket=-\frac{1}{4}\left\llbracket p^{3/4}u_{k}\right\rrbracket=-\frac{1}{4}\left\llbracket p\,u_{k}p^{-1/4}\right\rrbracket=-\frac{1}{4}\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\{\mkern-6.0mu\{$}}}p\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\}\mkern-6.0mu\}$}}}\left\llbracket u_{k}p^{-1/4}\right\rrbracket-\frac{1}{4}\left\llbracket p\right\rrbracket\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\{\mkern-6.0mu\{$}}}u_{k}p^{-1/4}\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\}\mkern-6.0mu\}$}}}.$$

$$\displaystyle 0=\left\llbracket\mathbf{h}\right\rrbracket\cdot\tilde{\mathbf{F}}_{k}-\left\llbracket\psi_{k}\right\rrbracket$$

$$\displaystyle\quad+\left(\frac{1}{8\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\{\mkern-6.0mu\{$}}}\sqrt{p^{-1/2}+\left(|\mathbf{u}|p^{-1/4}\right)^{2}}\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\}\mkern-6.0mu\}$}}}}\frac{\left\llbracket p^{-1/2}\right\rrbracket}{\left\llbracket p\right\rrbracket}\tilde{\mathbf{F}}_{k}^{(d+1)}+\frac{1}{4}\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\{\mkern-6.0mu\{$}}}u_{k}p^{-1/4}\raisebox{0.0pt}{\resizebox{0.0pt}{0.0pt}{\raisebox{2.5pt}{$\}\mkern-6.0mu\}$}}}\right)\left\llbracket p\right\rrbracket.$$

$${\tilde{\mathbf{F}}}_k^{(d+1)}$$

Full text

**AMS subject classification. 65M06, 65M20, 65M70 **

1 Introduction

Relativistic fluid dynamics is an area of ongoing study in recent research literature, covering various scales from ions to neutron stars [undefb]. The commonality is that the energy or velocity of the fluid in question is so significant that it surpasses the confines of classical Newtonian mechanics. For a recent review of relativistic fluid dynamics, we refer to [undefp]. In this work, we examine relativistic fluid dynamics in the context of hyperbolic partial differential equations (PDEs). The study of the structural properties of the governing equations, based on their symmetric structure, can be found in [undefaam, undefaal, undefaak, undefay]. For works exploring diverse aspects of the relativistic Euler equations, we refer to [undefaas, undefaar, undefj]. Additional recent results can be found in [undefz, undefaaq, undefar, undefu, undefv]. Here, we are particularly interested in the ultra-relativistic Euler equations, i.e., the scenario where thermal energy predominates. Results on this specific system can be located in the aforementioned literature. For a general introduction to the mathematical theory of hyperbolic conservation laws, see Smoller [undefaap] and Dafermos [undefm].

In this work, we focus on the numerical aspects of the relativistic flows under consideration. For a recent treatment of numerical methods for relativistic fluid dynamics, we refer to [undefaaab, undefw, undefr] and the references therein. Previous results on the numerical treatment of the ultra-relativistic Euler equations are presented by Abdelrahman et al. [undef], who propose a front tracking scheme; for kinetic schemes, we cite Kunik et al. [undefaq, undefap, undefao]. In Kunik et al. [undefan], the ultra-relativistic Euler equations are studied, and a one-dimensional radially symmetric solver that can provide non-trivial reference solutions for genuinely three-dimensional calculations is developed. Recently, in [undefam], this work continues, and multi–dimensional benchmark problems are introduced. These benchmarks are solved using full-dimensional simulations, and the results are compared to the specific solver provided in [undefan, undefam]. In this work, we extend the previous works by developing an entropy-stable method for the ultra-relativistic equations and applying it not only to two-dimensional but also to more challenging three-dimensional problems. The entropy-stable numerical methods developed in this paper are based on the framework pioneered by Tadmor [undefaaw, undefaav]. A key ingredient of the method is an entropy-conservative two-point flux, which we derive for the first time for the ultra-relativistic Euler equations. This flux is not based on integral averages like the first entropy-conservative numerical fluxes proposed by Tadmor, but rather on differential averages similar to the affordable numerical fluxes developed for the compressible Euler equations of an ideal gas in the non-relativistic setting [undefaf, undefi, undefaab]. The second-order method is extended to a higher-order scheme using flux differencing [undefat, undeft, undefaab, undefk] combined with summation-by-parts (SBP) operators [undefaau, undefs]. While many numerical schemes, including finite differences [undefak, undefaat], finite volumes [undefav, undefaw], continuous finite elements [undefad, undefac, undefa], flux reconstruction [undefae, undefaay, undefaag], active flux [undefd], and meshless methods [undefab] can be formulated using SBP operators, we focus specifically on discontinuous Galerkin (DG) methods [undefx, undefg, undefh]. Finally, we combine the entropy-stable baseline schemes with shock capturing using finite volume subcells [undefaa] and adaptive mesh refinement implemented in the open-source software Trixi.jl [undefaai, undefaan]. In particular, we provide all Julia source code and data needed to reproduce the numerical results presented in this paper in the accompanying reproducibility repository [undefaax].

The outline of this paper is as follows. In Section 2, the governing equations relevant to this study are provided. In Section 3, an entropy-conservative numerical flux is derived. Section 4 summarizes the numerical methods used. In Section 5, several numerical tests are provided, and the simulation results are presented. Concluding remarks are found in Section 6. Appendix A offers detailed and supplementary calculations for the readers’ convenience.

2 Conservative Formulations of the Equations

In this section we specify the governing equations precisely and provide necessary details on the conservative, primitive as well as entropy variables. Here we consider the flat Minkowski metric given by

[TABLE]

It is of further interest to study other space–times which then would result in curvi–linear coordinates, see for example [undefc, undefn]. A precise study of such cases is beyond the scope of this work and will be considered in the future. As mentioned before we study the ultra-relativistic equations for a perfect fluid in Minkowski space-time $\mathbf{x}=(x_{1},x_{2},x_{3})$, $t=x_{0}$, namely

[TABLE]

with $\alpha,\beta\in\{0,1,2,3\}$ and the energy-momentum tensor

[TABLE]

Here $p$ represents the pressure, $\mathbf{u}\in\mathbb{R}^{3}$ is the spatial part of the four-velocity vector $(u_{0},u_{1},u_{2},u_{3})=(\sqrt{1+|\mathbf{u}|^{2}},\mathbf{u})$. Note that in this work we consider the dimensionless velocity $\mathbf{u}$ and we have the following relation to the usual velocity $\mathbf{v}$

[TABLE]

where $c$ is the speed of light and we have the Lorentz factor

[TABLE]

Using the dimensionless velocity we hence have $\Gamma=\sqrt{1+|\mathbf{u}|^{2}}$ and we further scale to $c=1$. For the ideal ultra-relativistic gas the thermal energy dominates and thus the local energy density $e$ at rest is linked to the pressure $p$ by

[TABLE]

Note that (2.4) is universal for an ideal gas in the ultra-relativistic regime. For the physical background we refer to Weinberg [undefaaaa, Part I, pp 47-52], further details can be found in Kunik [undefal, Chapter 3.9] and in the works by Ruggeri and co-authors [undefay, undefaal, undefaak].

The conservative multi–dimensional $d+1$ form of the ultra–relativistic Euler equations (2.2) is presented in the following. We consider either $d=2$ or $d=3$ space dimensions. Then the unknown quantities $p$ and $\mathbf{u}=(u_{1},\ldots,u_{d})\in\mathbb{R}^{d}$ satisfying (2.2) depend on time $t\geq 0$ and position $\mathbf{x}=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$.

Putting $\alpha=0$ in (2.2) gives the conservation of energy

[TABLE]

whereas for $\alpha=j=1,\dots,d$ we obtain the conservation of momentum

[TABLE]

We want to rewrite the conserved variables $\mathbf{w}\mathrel{\mathop{\ordinarycolon}}=(w_{1},\ldots,w_{d},w_{d+1})=(\overline{\mathbf{w}},w_{d+1})\in\mathbb{R}^{d+1}$, $\overline{\mathbf{w}}\mathrel{\mathop{\ordinarycolon}}=(w_{1},\ldots,w_{d})\in\mathbb{R}^{d}$ in terms of the primitive variables determined by the velocity $\mathbf{u}\mathrel{\mathop{\ordinarycolon}}=(u_{1},\ldots,u_{d})\in\mathbb{R}^{d}$ and the pressure $p$, i.e.,

[TABLE]

Reversely, the primitive variables can be rewritten in terms of the conserved variables as

[TABLE]

Concerning the entropy for this system we want to remark the following. First, it should be noted that the equation for the particle number density $n$ decouples from the system in the ultra–relativistic regime, cf. [undefaal, undefaaaa, undefal]. This provides the opportunity to find an entropy, additionally to the one from the full system, particularly for the system under consideration. According to [undefal] the system is equipped with a concave physical entropy given by

[TABLE]

and the corresponding flux reads

[TABLE]

To justify this entropy one can use arguments from [undefq] as well as from Weinberg [undefaaaa, Chap. 10, p. 51]. In order to be self–contained we provide the proof that the presented entropy is a mathematical entropy for this system in the Appendix A. In particular we show concavity for the physical entropy and the relation of the entropy and the entropy flux, i.e.,

[TABLE]

where $\mathbf{D}_{\mathbf{w}}\mathbf{F}_{k}$ denotes the Jacobian of the flux with respect to the conservative variables and $\mathbf{F}_{k}$ denotes the flux of the system in $x_{k}$ direction for $k=1,\dots,d$. The entropy variables $\mathbf{h}=(h_{1},\dots,h_{d+1})^{T}$ which are also known as main field, see [undefaal], are given by

[TABLE]

The main field of the present system may be obtained from the main field of the relativistic Euler equations since the ultra–relativistic Euler equations are a subsystem of the latter. We refer the interested reader to [undefaak, undefaal] and to the recent survey [undefaaz]. To construct the entropy-stable numerical flux, cf. [undefaaw, undefk], we need to calculate the entropy - flux potential. The corresponding flux is given by

[TABLE]

3 Derivation of an entropy-conservative numerical flux

In this section we present the detailed derivation of the entropy-conservative numerical flux. Following the seminal works of Tadmor [undefaaw, undefaav], the condition for an entropy-conservative numerical flux reads

[TABLE]

where $\mathbf{h}$ are the entropy variables given by (2.12), $\tilde{\mathbf{F}}_{k}$ is the numerical flux in space direction $k$, $\psi_{k}$ is the flux potential in direction $k$ defined in (2.13), and \left\llbracket a\right\rrbracket=a^{+}-a^{-} denotes the jump of a quantity $a$. The numerical fluxes $\tilde{\mathbf{F}}_{k}$ should be consistent approximations of the fluxes

[TABLE]

of the ultra-relativistic Euler equations (2.5) and (2.6) which form a conservative system

[TABLE]

We follow the algorithmic approach of [undefaab] to derive an entropy-conservative numerical flux for the ultra-relativistic Euler equations. Thus, we choose a set of variables and express all jumps in terms of these variables. We then solve the resulting system of equations in order to obtain entropy-conservative numerical fluxes which are “affordable” in the sense of [undefaf].

To express jumps of general expressions, we use mean values such as the arithmetic mean

[TABLE]

Applying straight-forward algebraic manipulations we obtain the analogues of well-known rules for differentiation, i.e.,

[TABLE]

and for positive $a$ we also have

[TABLE]

Moreover, we will use

[TABLE]

i.e.,

[TABLE]

Because of the structure of the entropy variables $\mathbf{h}$, we choose the variables $u_{i}p^{-1/4}$ and $p$ to express all jumps. Then, we get

[TABLE]

The jump of the entropy variable $h_{d+1}$ is given by

[TABLE]

where (3.6) has been used. Moreover, the ratio of the jumps of the pressure terms can be expressed using (3.8). For the jump of the flux potential we get

[TABLE]

Thus, the condition for the entropy-conservative numerical flux in the present situation reads

[TABLE]

This equation is satisfied if we choose the numerical flux $\tilde{\mathbf{F}}_{k}$ so that each term in parentheses vanishes. First, we obtain the numerical energy flux

[TABLE]

This is clearly a consistent approximation of the energy flux $\mathbf{F}_{k}^{(d+1)}=4pu_{k}\sqrt{1+|\mathbf{u}|^{2}}$. Inserting this numerical energy flux into the other term in parentheses, we get

[TABLE]

This term vanishes if we choose the numerical momentum fluxes

[TABLE]

This is a consistent approximation of the momentum flux $\mathbf{F}_{k}^{(i)}=4pu_{i}u_{k}+\delta_{ik}p$ of the ultra-relativistic Euler equations. Thus, we have proven

Theorem 3.1**.**

The numerical fluxes

[TABLE]

are symmetric, consistent, and entropy-conservative two-point fluxes for the ultra-relativistic Euler equations (2.6) and (2.5).

Remark 3.2**.**

The pressure in the momentum flux of (3.16) is approximated using an arithmetic mean. For the classical compressible Euler equations, this is required for the stronger version of kinetic energy preservation [undefag, undefi] discussed in [undefaac, undefaad], see also [undefaaf, undefaaa]. Even without provable entropic properties, kinetic energy-preserving split forms have been widely used for compressible computational fluid dynamics and turbulence simulations because of the increased robustness [undefy, undefaao, undefah].

4 Numerical Methods

We want to compare numerical results for radially symmetric solutions of the original multi–dimensional initial-value problem (IVP) for the ultra-relativistic Euler equations (2.5)–(2.6). To avoid confusion of the pressure in the multi–dimensional case and the pressure in the radially symmetric case we denote the pressure in the multi–dimensional case by $\tilde{p}$ from now on. For $j=1,\dots,d$ the IVP reads

[TABLE]

We require radially symmetric solutions of this system, i.e., for $t\geq 0$ and $x=|\mathbf{x}|>0$ the restrictions for pressure and velocity are given by

[TABLE]

For $t=0$ we obtain given initial data $\tilde{p}(0,\mathbf{x})=p(0,x)=p_{0}(x)$ and

[TABLE]

and for $\mathbf{x}=\mathbf{0}\in\mathbb{R}^{d}$ we may also put $\mathbf{u}(0,\mathbf{0})=\mathbf{0}$ and $v_{0}(0)=0$.

4.1 Multi–dimensional discontinuous Galerkin methods

We use the numerical simulation framework Trixi.jl [undefaai, undefaan] implemented in Julia [undefe] to solve the multi–dimensional IVP (4.1) using discontinuous Galerkin methods. We follow the standard approach to first partition the spatial domain into quadrilateral/hexahedral elements with disjoint interiors. In each element, the numerical solution is represented by its values at tensor product Gauss-Lobatto-Legendre nodes able to represent polynomials of degree three in each coordinate direction. The associated quadrature rule and differentiation matrices are used to compute integrals and derivatives, respectively, i.e., we use a collocation approach. Thus, the resulting scheme is a high-order accurate nodal discontinuous Galerkin spectral element method (DGSEM) [undefai].

DG schemes are often derived from a weak form of the equations. First, multiply the conservation law $\partial_{t}\mathbf{w}+\sumop\slimits@_{k=1}^{d}\partial_{x_{k}}\mathbf{F}_{k}(\mathbf{w})=0$ by a test function $\varphi$ and integrate over an element $E$ to get

[TABLE]

Next, integrate the volume term by parts to get

[TABLE]

where $n_{k}$ is the $k$-th component of the outward normal vector on the boundary $\partial E$ of the element $E$. We replace the flux in the surface integral by a numerical flux $\hat{\mathbf{F}}$ to get

[TABLE]

where $\mathbf{w}_{-}$ and $\mathbf{w}_{+}$ are the values of the solution on the interior and exterior of the element, respectively, and $\mathbf{n}$ is the outward normal vector on the boundary $\partial E$ of the element $E$. Concretely, we use a local Lax-Friedrichs/Rusanov flux $\hat{\mathbf{F}}$.

Flux differencing [undefat, undeft, undefaab, undefk] is more similar to the strong form DG scheme. Thus, we integrate the volume term once more by parts to get

[TABLE]

A classical strong-form DGSEM is obtained by replacing the volume flux $\mathbf{F}_{k}(\mathbf{w})$ by the polynomial approximation obtained by collocation at the Gauss-Lobatto-Legendre nodes, replacing integrals by the associated quadrature rule [undefai], and requiring the resulting equations to hold for all test functions $\varphi$ in the same polynomial space. The resulting volume term can be written as

[TABLE]

where the indices $i$ and $j$ run over the Gauss-Lobatto-Legendre nodes in the element and $D_{k}$ is the discrete (polynomial) differentiation matrix in the $k$-th coordinate direction.

Flux differencing methods [undefat, undeft, undefaab, undefk] replace the volume term above by

[TABLE]

where $\tilde{\mathbf{F}}_{k}$ is a two-point numerical flux. An efficient implementation of such flux differencing volume terms is described in [undefaah]. We use the newly developed entropy-conservative two-point flux given in Theorem 3.1 for the volume term. Thus, the flux differencing theory [undefat, undeft, undefaab, undefk] guarantees that we obtain an entropy-stable high-order method.

To ensure robustness also in the presence of discontinuities, we use the finite volume subcell shock capturing approach of [undefaa]. There, the high-order flux differencing volume terms are interpreted as high-order subcell finite volume methods. For robustness, these high-order methods are blended with a first-order finite volume method using a local Lax-Friedrichs/Rusanov flux. The blending factor determining the convex combination of high- and low-order methods is based on a shock indicator which is applied to the variable $\tilde{p}\sqrt{1+|\mathbf{u}|^{2}}$. For efficiency, we use adaptive mesh refinement (AMR) based on p4est [undeff] and an adaptation of the indicator of [undefau] applied to the same indicator variable used for shock capturing.

The resulting system of ordinary differential equations is integrated in time using the third-order, four-stage explicit strong stability preserving Runge-Kutta method of [undefaj] with embedded method of [undefl] and adaptive step size controller of [undefaae] implemented in OrdinaryDiffEq.jl [undefaz]. Thus, the step size is adapted automatically to the maximal resolution of the mesh and the local wave speeds, see also [undefaaj] for further details on the error control and step size adaptation. For the self-similar expansion (Example 2 below), we also apply the positivity-preserving limiter of Zhang and Shu for the pressure [undefaaac].

All data and the Julia source code required to reproduce the numerical results presented in this paper are available online in our reproducibility repository [undefaax].

4.2 One-dimensional radially symmetric solver RadSymS

For a detailed presentation of the scheme RadSymS used to compute the solution to the quasi one-dimensional radially symmetric form of the ultra-relativistic Euler equations we refer to [undefam] and [undefan]. The solution of the radially symmetric quasi one-dimensional problem satisfies

[TABLE]

for $t>0$ and $x>0$ representing the radial direction and with the initial conditions

[TABLE]

Here $d=2$ or $d=3$ denotes the space dimension and the variables $(a,b)$ are linked to $(p,u)$ through a one-to-one transformation on the sets

[TABLE]

The mapping $\Theta\mathrel{\mathop{\ordinarycolon}}\mathcal{S}\to\tilde{\mathcal{S}}$ and its inverse are given by

[TABLE]

Using the transformation (4.6), its inverse (4.7) and the velocity

[TABLE]

we replace the state variables $a$ and $b$ by $p$ and $v$, respectively. The variable $c=c(a,b)$ is given by

[TABLE]

Observe, that the variables $a$ and $b$ correspond to the conservative quantities in radial direction, while $c(a,b)=p(1+4u^{2})$.

We are looking for weak solutions $a=a(t,x)$, $b=b(t,x)$ with $|b|<a$, i.e.

[TABLE]

Here $\Omega\subset Q$ is a convex domain with piecewise smooth boundary $\partial\Omega$ with $Q=\{(t,x)\,|\,t>0,\;x>0\}$ and the integral formulation is based on [undefax, undefal]. Note that the integrals along $\partial\Omega$ are line integrals over a 1-form written in their local coordinates, cf. [undefas, undefo].

Based on this reformulation of the problem the one-dimensional numerical scheme (RadSymS) for the initial value problem of the radially symmetric ultra-relativistic Euler equations in two and three space dimensions is developed. It further can be shown that this scheme preserves positive pressure, see [undefam, undefan]. The method of contour-integration for the formulation of the balance laws (4.10) is used to construct a function called “Euler”. This function enables us to obtain the time evolution of the numerical solution on a staggered grid. More precisely it allows us to construct the solution $(a^{\prime},b^{\prime})$ at the next time step from the solution $(a_{\pm},b_{\pm})$ in two neighboring grid points at the previous time step according to Figure 1. First we determine the computational domain and define some quantities which are needed for its discretization.

Given are $t_{*},x_{*}>0$ in order to calculate a numerical solution of the initial value problem (4.10), with the initial conditions given in (4.4), in the time range $[0,t_{*}]$ and the spatial range $[0,x_{*}]$. 2. 2)

We want to use a staggered grid scheme. Any given number $N\in\mathbb{N}$ with $N\cdot x_{*}\geq t_{*}$ determines the time step size $\Delta t=t_{*}/2N$ and the time steps are $t_{n}=(n-1)\Delta t,n=1,\ldots,2N+1$. 3. 3)

Put

[TABLE]

then the spatial mesh size is $\Delta x=x_{*}/M$ with the spatial grid points $x_{j}=(j-1)\Delta x\,,j=1,\ldots,N+M+1$. Note that our scheme uses a trapezoidal computational domain $\mathcal{D}$ defined below that includes the target domain $[0,t_{*}]\times[0,x_{*}]$. Thereby, we can use all initial data that influence the solution on the target domain. In this way we avoid using a numerical boundary condition at $x_{*}$. 4. 4)

The number $\lambda=\Delta x/(2\Delta t)\geq 1$ is used to satisfy the CFL-condition and to define the computational domain

[TABLE]

For the numerical discretization of the integral balance laws (4.10) we choose the triangular balance domain depicted in Figure 1. We assume that the midpoints $P_{-}=(\overline{t},\overline{x}-\Delta x/2)$, $P_{+}=(\overline{t},\overline{x}+\Delta x/2)$ and $P^{\prime}=(\overline{t}+\Delta t,\overline{x})$ of the cords of $\partial\Omega$ are numerical grid points for the computational domain $\mathcal{D}$. Let the numerical solution $(a_{\pm},b_{\pm})$ be given at the grid points $P_{\pm}$. We have to require $|b_{\pm}|<a_{\pm}$ for the numerical solution in the actual time step $\overline{t}=t_{n}$ with $n=1,\ldots,2N$. The major task is to calculate the numerical solution $(a^{\prime},b^{\prime})$ for the next time step $\overline{t}+\Delta t=t_{n+1}$ at its grid point $P^{\prime}$, see Figure 1.

The spatial value $\overline{x}\geq 0$ is given. We have to determine a function

[TABLE]

for the calculation of $(a^{\prime},b^{\prime})$. This leads to the structure of a staggered grid scheme. Note that at the boundary the balance region may have parts outside $\mathcal{D}$, e.g. points below the half-space $x\geq 0$. In the latter case we employ a simple reflection principle for the numerical solution in order to use the function Euler as well for the boundary points with $\overline{x}=0$.

Next we make use of the fact that the points $P_{\pm}$ with the numerical values $(a_{\pm},b_{\pm})$ and $P^{\prime}$ with the unknown value $(a^{\prime},b^{\prime})$ are the midpoints of the three boundary cords of the balance region .

Theorem 4.1** (Numerical solution $(a^{\prime},b^{\prime})$ for the balance region ).**

Given are real quantities $\overline{x}\geq 0$ and $a_{\pm}$, $b_{\pm}$. Assume that $|b_{\pm}|<a_{\pm}$. We recall $\lambda\geq 1$ defined in terms of $\Delta t$ and $\Delta x$ and put $c_{\pm}=c(a_{\pm},b_{\pm})$ in (4.9).

[TABLE]

Then we have $|b^{\prime}|<a^{\prime}$ in both cases.

Definition 4.2** (The function Euler).**

The state $(a^{\prime},b^{\prime})$ from Theorem 4.1 defines the function Euler in (4.11).

Now we formulate the numerical scheme for the solution of the initial-boundary value problem (4.10). We construct staggered grid points in the computational domain $\mathcal{D}$ and compute the numerical solution at these grid points. Using the function Euler we obtain the evolution of the numerical solution in time, i.e., it allows us to construct the solution at time $t=t_{n+1}$ from the solution which is already calculated in the grid points at the former time step $t=t_{n}$.

1. (I)

The staggered grid points are $(t_{n},x_{n,j})\in\mathcal{D}$ for $t_{n}=(n-1)\Delta t$, $n=1,\ldots,2N+1$ and $j=1,\ldots,M+N-\lfloor(n-1)/2\rfloor$ with

[TABLE]

We want to calculate the numerical solution $(a_{n,j},b_{n,j})$ at $(t_{n},x_{n,j})$. 2. (II)

For $j=1,\ldots,M+N$ we calculate the numerical solution $(a_{1,j},b_{1,j})$ at the grid point $(t_{1},x_{1,j})=(0,(x_{j}+x_{j+1})/2)$ from the given initial data by

[TABLE]

This corresponds to taking the integral average of the initial data on $(x_{j},x_{j+1})$ and using the midpoint rule as quadrature. 3. (III)

Assume that for a fixed odd index $n\in\{1,\ldots,2N\}$ we have already determined the numerical solution $(a_{n,j},b_{n,j})$ at the grid points $(t_{n},x_{n,j})$, $j=1,\ldots,M+N-(n-1)/2$.

First we determine the solution $(a_{n+1,1},b_{n+1,1})$ at the boundary point $(t_{n+1},x_{n+1,1})=(t_{n+1},0)$ according to Thm. 4.1. For this purpose we put $a_{+}=a_{n,1}$, $b_{+}=b_{n,1}$, $a_{-}=a_{n,1}$, $b_{-}=-b_{n,1}$ and have

[TABLE]

Next we put $a_{-}=a_{n,j-1}$, $b_{-}=b_{n,j-1}$ and $a_{+}=a_{n,j}$, $b_{+}=b_{n,j}$ for $j=2,\ldots,M+N-(n-1)/2$ and determine the values $a_{n+1,j}$, $b_{n+1,j}$ at time $t_{n+1}$ and position $\overline{x}=x_{n+1,j}=x_{j}$ from

[TABLE] 4. (IV)

Assume that for a fixed even index $n\in\{1,\ldots,2N\}$ we have already determined the numerical solution $(a_{n,j},b_{n,j})$ at the grid points $(t_{n},x_{n,j})$, $j=1,\ldots,M+N-n/2+1$.

We put $a_{-}=a_{n,j}$, $b_{-}=b_{n,j}$ and $a_{+}=a_{n,j+1}$, $b_{+}=b_{n,j+1}$ for $j=1,\ldots,M+N-n/2$ and determine the values $a_{n+1,j}$, $b_{n+1,j}$ at time $t_{n+1}$ and position $\overline{x}=x_{n+1,j}=(x_{j}+x_{j+1})/2$ from

[TABLE]

The quasi one-dimensional problem (4.4) is approximately solved by the one-dimensional radially symmetric scheme with $N=5000$. We prescribe the initial pressure $p_{0}=p(0,\cdot)$ as well as the initial velocity $v_{0}=v(0,\cdot)$. The variable $v$ is chosen, because the restriction $|v|<1$ leads to better color plots (the variable $u$ is unbounded). Note that for constant pressure $p_{0}>0$ and $v_{0}=0$, we obtain a stationary solution. Such a steady part is contained in the solutions to the examples presented below.

5 Numerical Results

In the following we provide the results for the benchmarks proposed and studied in [undefam]. We computed numerical solutions in two and three space dimensions using the genuine multi–dimensional solver Trixi.jl and compare the obtained results with the one-dimensional radially symmetric scheme RadSymS presented in [undefan, undefam]. Note that all of the test cases are radially symmetric problems where two of them provide self-similar solutions, see Example 1 and 2.

In Figs. 2–19 for Examples 1–5, respectively, we present 2D and 3D results in the $t$–$x$ plane and a line plot of the solution in radial direction at the final time $t_{end}$. We compare the results by means of the pressure $p$ and the velocity $v$ (4.8). Additionally, for Examples 1 and 2 we compare the solutions also with the reference solution that can be determined for self-similar problems by solving the ODE (5.2) in the smooth part and using the jump conditions (5.4) to determine the stationary part in Example 1.

For all examples presented subsequently, we summarize the computational domain $\Omega\subset\mathbb{R}^{d}$ and the maximum refinement level $L\in\mathbb{N}$ in Table 1. The initial cell size $\Delta\mathbf{x}_{0}\in\mathbb{R}^{d}$ at the coarsest refinement level is always the whole computational domain . Due to the dyadic grid hierarchy in Trixi.jl the cell size at refinement level $l$ is $\Delta\mathbf{x}_{l}=2^{-l}\Delta\mathbf{x}_{0}$, for $l=0,\dots,L$.

We want to remark that the refinement level in each 3D computation is less than the level used for the 2D computation of the same example, respectively. Nevertheless, the computational time increases dramatically with the space dimension. Therefore, the final results for the pressure in examples 3, 4, and 5 in 3D cannot match the excellent agreement of the corresponding 2D results due to the computational costs and complexity. Here we decided to provide results that can be obtained within hours of computational time on a standard workstation rather than using HPC resources.

Convergence test. To verify the implementation and demonstrate the high-order accoracy of the multi–dimensional solver Trixi.jl, we perform a convergence test for a smooth manufactured solution of the ultra-relativistic Euler equations. We choose the solution

[TABLE]

where $\alpha_{2}=\sqrt{10}{4}$ and $\alpha_{3}=\sqrt{15}{4}$. The resulting source terms for the momentum equations are

[TABLE]

and no source term has to be added to the energy equation. We use a uniformly refined mesh with periodic boundaries and compute the solution until the final time $t_{end}=1$. The time integration tolerance is set to $10^{-8}$ and the time step is adapted automatically by the step size controller of [undefaae] to ensure that the temporal error is negligible compared to the spatial error. The resulting discrete $L^{2}$ errors and the corresponding experimental order of convergence (EOC) are presented in Tables 2 and 3 for the 2D and 3D version of Trixi.jl, respectively. As expected for polynomials of degree three, the EOC is close to four.

Example 1: Solutions including a shock and a stationary part. Following [undefar], self-similar solutions can be constructed solving an ODE system in radial direction $x$. These solutions are in particular constant along rays $\xi=x/t$ or, equivalently, $\vartheta=t/x$.

We consider constant initial data with pressure $p_{0}=1$ and radial velocity $v_{0}\in(-1,0)$. Due to [undefar, Section 2.3] there is a solution $p(t,x)=P(\vartheta)$ and $v(t,x)=V(\vartheta)$ depending only on $\vartheta=t/x$ for $t,x>0$, with a single straight line shock emanating from the zero point. Let $\tilde{s}\in(0,1)$ be the unknown constant shock speed. Then we put $v_{-}=0$ and can find an unknown pressure $p_{-}>0$ with

[TABLE]

Due to the Rankine-Hugoniot shock conditions introduced in [undefal, Section 4.4] we obtain from [undefal, page 82] for a so called 3-shock after a lengthy calculation the algebraic shock conditions

[TABLE]

[TABLE]

Here $p_{-},v_{-}$ are the values of pressure and velocity left to the 3-shock, and $p_{+},v_{+}$ are the values of pressure and velocity right to the 3-shock. Due to Lai [undefar, Section 2.3] the solution $p=P(\vartheta)$, $v=V(\vartheta)$ satisfies the IVP of ordinary differential equations

[TABLE]

Moreover, it is shown in the appendix of [undefam] that Lai’s results guarantee a unique solution for $0\leq\vartheta<\vartheta_{max}$ with a certain value $\vartheta_{max}>\sqrt{3}$ and a unique value $\sqrt{3}<\tilde{\vartheta}<\sqrt{v_{0}^{2}+3}-v_{0}<3$, see [undefam, (A.10)] such that

[TABLE]

After having determined $\tilde{\vartheta}$ we finally obtain

[TABLE]

For the numerical simulation we prescribe a constant initial pressure $p_{0}=1$ and a constant initial velocity $v_{0}=-1/\sqrt{2}$. This corresponds to the following initial data for the original IVP (4.1)

[TABLE]

The radially symmetric solution of (4.1) is determined by the solution of the ODE system (5.2) which gives a shock wave moving at constant speed $\tilde{s}$. We want to point out that these ODEs are written in terms of $\vartheta=t/x$. Thus, the initial data $P(0)$ and $V(0)$ for (5.2) prescribe the solution for (4.1) at infinity given by $p_{0}=P(0)$ and $v_{0}=V(0)$. The ODE system (5.2) is solved by applying the classical fourth-order RK scheme with step size $h=10^{-6}$. For the computation of $\tilde{\vartheta}$ the ODE solver is run until (5.3) is satisfied with a tolerance of $\varepsilon=10^{-9}$. The solution values for the shock speed, the states ahead and behind of the shock, respectively, are presented in Table 4.

The results reported in Fig. 2 (2D) and Fig. 3 (3D) clearly show an excellent agreement of both numerical methods with the reference solution provided by an ODE solver for the final time $t_{end}=1$ in radial direction. Due to the self-similar structure of the solution we omit the presentation of results in the $t$–$x$ plane.

In order to verify the entropy conservation property of the presented scheme we provide a numerical test for slightly modified initial data of the present example, i.e., we use a constant initial pressure $p_{0}=1$ and a constant initial velocity $v_{0}=-1/\sqrt{26}$. This corresponds to the following initial data for the original IVP (4.1)

[TABLE]

To verify the entropy conservation property, we compute the entropy rate, i.e., the discrete version of

[TABLE]

where the integral is replaced by the quadrature associated to the DGSEM, and the time derivative is given by the entroy-conservative semidiscretization. The results for the entropy rate over time in 2D and 3D are given in Fig. 4. Clearly entropy is conserved up to machine precision. Checking discrete entropy conservation like this is a demanding test for the derivation of the entropy-conservative flux in Theorem 3.1 as well as for its implementation.

Example 2: Self-Similar Expansion. In this example we prescribe initial data leading to a smooth self-similar solution by applying Lai’s approach as in Example 1. Note that for $d=3$ the solution will expand into vacuum at the speed of light when the initial value for $v_{0}$ is close enough to one. We do not consider this case here and refer to [undefar] for more details. For the numerical simulation we consider a constant initial pressure $p_{0}=1$ and a constant initial velocity $v_{0}=1/\sqrt{5}$. This corresponds to the following initial data for the original IVP (4.1)

[TABLE]

The solution consists of a rarefaction wave that is determined by the ODE (5.2). In contrast to Example 1 there is no shock.

The results reported in Fig. 5 for 2D and Fig. 6 for 3D, respectively, clearly show an excellent agreement of both numerical methods with the reference solution provided by an ODE solver for the final time $t_{end}=1$ in radial direction. Due to the self-similar structure of the solution we omit the presentation of results in the $t$–$x$ plane.

Example 3: Expansion of a Spherical Bubble. Next we consider the expansion of a spherical bubble with the following initial data

[TABLE]

These initial values correspond to the following initial values for the original IVP (4.1)

[TABLE]

Initially, the pressure inside the bubble is ten times larger than outside, which leads to a fast expansion of the bubble into the outer low pressure area. This in turn gives rise to the formation of another low pressure area. We observe the formation of a shock wave, running towards the new low pressure area and reaching the zero point around time $t=5.032$ for $d=2$ and $t=4.165$ for $d=3$. The formation of this new shock wave is a peculiar nonlinear phenomenon. Shortly before the shock reaches the zero point the pressure takes very low values, but its reflection from the zero point depicted in Fig. 8 and Fig. 9 indicates a blow-up of the pressure in a very small time-space range near the boundary. This is similar to [undefan, Section 5, Example 4], but here the blow-up is weaker than in three space dimensions, and its illustration requires a higher numerical resolution. A similar behavior for the corresponding solution of the explosion problem with the classical Euler equations should be expected.

Before discussing the numerical results in details, we show the evolution of the total entropy $\intslop\ilimits@\eta(\mathbf{w}(t,\mathbf{x}))\,\textup{d}\mathbf{x}$ (computed numerically via the quadrature rule associated to the method) over time for the present example in Fig. 7. In contrast to the entropy results shown earlier in Fig. 4, we use the the local Lax-Friedrichs/Rusanov flux and do not show the discrete entropy rate, but the total entropy itself. Since the entropy $\eta$ is concave, entropy stability means that the total entropy increases over time. Clearly, the numerical results shown in Fig. 7 confirm the entropy stability of the presented scheme, as expected based on the theory of flux differencing and our construction of the entropy-conservative flux in Theorem 3.1.

From the depicted solution at time $t=6$ we observe the reflected shock curve at radius around $x=0.55$ for $d=2$ and near $x=1.21$ for $d=3$, respectively. The results in the $t$–$x$ plane given in Fig. 8 and Fig. 9 show the expansion of the bubble and the formation of the new shock which focuses at the origin and then is reflected again. At the focus point the pressure rises drastically in a small vicinity of the origin and away from that the pressure values do not vary a lot. Therefore, we present a zoom plot of the pressure. We observe an excellent agreement of both numerical methods for the velocity in Figs. 8 and 9 (b). For the pressure it is visible in 2D that Trixi.jl is capable of resolving the solution close to the reflection point of the shock at the origin rather well, see Fig. 8 (a). However, in 3D the reflection point is not captured well, which is due to the computational efforts needed in 3D as stated before. Indeed, the maximum pressure value computed by Trixi.jl is $p_{Trixi}\approx 9.6562$ versus $p_{RadSymS}\approx 289.2772$ for the case in 3D, see Fig. 9 (a). This is due to the lower resolution we for the 3D results compared to the 2D results, see Table 1. Matching the maximum pressure value also in 3D would require at least the same resolution as in 2D, which was not feasible for the workstation we could use for this project. Further note that this pressure peak occurs in the very vicinity of the origin and is very much localized in time. Thus we present zoom plot for the pressure to visualize it. For the velocity this is not needed since we have proper scaling properties due to $|v|<1$.

Additionally, we compare both methods at final time $t_{end}=6$ where the solutions coincide very well, see Figs. 10 and 11.

Example 4: Collapse of a Spherical Bubble. Next we study the collapse of a spherical bubble with the following initial data

[TABLE]

These initial values correspond to the following initial values for the original IVP (4.1)

[TABLE]

The results in the $t$–$x$ plane, see Figs. 13 and 14, show the collapse of the bubble with a focus point at the origin which then is reflected again. At the focus point the pressure rises drastically in a small vicinity of the origin and away from that the pressure values do not vary a lot. Thus, we again present a zoom plot of the pressure for this example. We observe for both, $d=2$ and $d=3$, an excellent agreement of both numerical methods for the velocity in Figs. 13 and 14, respectively. For the pressure it is again visible that in 2D the results agree well whereas due to the computational complexity in three dimensions the pressure peak is not very well resolved. Indeed, the maximum pressure value computed by Trixi.jl is $p_{Trixi}\approx 390$ versus $p_{RadSymS}\approx 3550$ for the case in 3D, see Fig. 14 (a). This is again caused by the lower resolution we for the 3D results compared to the 2D results, see Table 1. Matching the maximum pressure value also in 3D would require at least the same resolution as in 2D, which was not feasible for the workstation we could use for this project.

Additionally, we compare both methods at time $t_{end}=6$ where the solutions coincide very well, see Figs. 12 and 15.

Example 5: Initially Sine–Shaped Radial Velocity. Finally, we study initial data with a velocity in radial direction described by a sine function, i.e.,

[TABLE]

These correspond to the following initial data for the multi–dimensional simulation

[TABLE]

The results in the $t$–$x$ plane, see Figs. 16 and 17, show a complex wave structure for the velocity and a high pressure focus at around time $t=0.77$ for $d=2$ and $t=0.71$ for $d=3$, respectively. Again, we present a zoom plot of the pressure for this example. As for the previous examples it is again visible near the reflection point that in 2D the results agree well whereas due to the computational complexity in three dimensions the pressure peak is not very well resolved. Indeed, the maximum pressure value computed by Trixi.jl is $p_{Trixi}\approx 503$ versus $p_{RadSymS}\approx 5957$ for the case in 3D, see Fig. 17 (a). This is again due to the lower resolution we for the 3D results compared to the 2D results, see Table 1. Matching the maximum pressure value also in 3D would require at least the same resolution as in 2D, which was not feasible for the workstation we could use for this project.

We observe an excellent agreement of both numerical methods for the velocity in Figs. 16 and 17 (b). Additionally, we compare both methods at final time $t_{end}=6$ where the solutions coincide very well, see Figs. 18 and 19.

6 Conclusion

We have investigated the system of the ultra–relativistic Euler equations in two and three space dimensions. For the first time an entropy-conservative flux has been derived and applied for numerical computations. We studied five different benchmarks in two and three dimensions, where the full-dimensional system is treated using Trixi.jl and the results are compared to numerical reference solutions. It should be emphasized that this is the first time that genuine three-dimensional results for these benchmarks have been obtained and presented. The obtained numerical results coincide well, although they emphasize the need to reduce the computational complexity for the three-dimensional case. Despite the interest in the particular system under consideration, these benchmarks may serve as challenging benchmarks for other multi–dimensional solvers. It will be interesting to investigate the extension of the present work to different space–times.

Acknowledgments

Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through SPP 2410 Hyperbolic Balance Laws in Fluid Mechanics: Complexity, Scales, Randomness (CoScaRa) within the Projects 525939417 (FT) and 526031774 (HR). HR additionally acknowledges the DFG individual research grants 513301895 and 528753982. FT thanks Matthias Kunik for introducing him to the interesting topic of relativistic Euler equations and fruitful discussions. HR thanks Joshua Lampert for reviewing HR’s contributions to Trixi.jl that improved the efficiency of radially symmetric output data.

Data Availability Statement

All Julia source code and data needed to reproduce the numerical results presented in this paper are available in the accompanying reproducibility repository [undefaax].

Appendix A Derivatives and Supplementary Calculations for the Entropy

In this section we want to provide further details on the derivatives. Furthermore, we show the concavity of the entropy and prove relation (2.11). For simplicity, we will express the derivatives in terms of the primitive variables, which is more convenient to work with.

Using (2.8b) and (2.7) we can calculate the derivatives of the pressure

[TABLE]

The derivatives of the velocity components $u_{i},i=1,\dots,d$ can be obtained using (2.8a) and (2.7)

[TABLE]

We yield for the derivatives of the entropy with respect to the primitive variables

[TABLE]

Note that corresponding to the flux potential (2.13) we can also calculate the entropy potential. Sometimes it is also called generating potential; we refer to the recent survey [undefaaz] and the references therein. We obtain for the entropy potential

[TABLE]

Next, we provide the second order derivatives of $\eta$ with respect to the conservative variables, i.e.,

[TABLE]

For simplicity, we are going to show the convexity of $\tilde{\eta}=-\eta$ and the Hessian can be written as follows

[TABLE]

with

[TABLE]

To prove convexity we now want to show that

[TABLE]

Therefore, we first note that all elements $(v_{1},\dots,v_{d},0)^{T}\in\mathbb{R}^{d+1}$ belong to $\ker(\mathbf{A}_{3})$ and conversely all elements $(0,\dots,0,v_{d+1})^{T}\in\mathbb{R}^{d+1}$ belong to $\ker(\mathbf{A}_{1})$ as well as to $\ker(\mathbf{A}_{2})$. Clearly, $\mathbf{A}_{1}>0$ for all $(v_{1},\dots,v_{d},0)^{T}\in\mathbb{R}^{d+1}\setminus\{\mathbf{0}\}$. For $\mathbf{A}_{2}$ we first consider the non–zero block matrix and verify that the rows are pairwise linearly dependent, since the $i$-th row is written as $u_{i}\mathbf{u}$. Thus, we have that $\mathbf{A}_{2}\geq 0$ for all $(v_{1},\dots,v_{d},0)^{T}\in\mathbb{R}^{d+1}\setminus\{\mathbf{0}\}$. For $\mathbf{A}_{3}$ we yield when multiplying with elements $\tilde{\mathbf{v}}=(0,\dots,0,v_{d+1})^{T}\in\mathbb{R}^{d+1}\setminus\{\mathbf{0}\}$

[TABLE]

and thus $\mathbf{A}_{3}>0$ for all $(0,\dots,0,v_{d+1})^{T}\in\mathbb{R}^{d+1}\setminus\{\mathbf{0}\}$. Altogether we hence have shown the convexity of $\tilde{\eta}=-\eta$.

We now want to verify relation (2.11). Therefore, we first calculate the derivatives of the entropy flux with respect to the primitive quantities. These are given by

[TABLE]

From these we calculate the derivatives with respect to the conservative variables

[TABLE]

We further need the Jacobians of the conservative fluxes $\mathbf{F}_{k}$ in the normal directions w.r.t. $x_{k}$. To this end we refer to [undefam] where the Jacobians are given for any normal direction, see [undefam, (B.5)]. Here the situation simplifies, since we only consider the normal directions $\mathbf{n}=\mathbf{e}_{k},k=1,\dots,d$. The flux Jacobian is determined by

[TABLE]

with

[TABLE]

If (2.11) is written component wise we yield

[TABLE]

Considering $i=1,\dots,d,i\neq k$ yield

[TABLE]

Now we deal with $i=k=1,\dots,d$ and obtain

[TABLE]

We now consider the final case and get

[TABLE]

Thus, (2.11) is verified.