© The Authors 2025

1 Introduction

Turbulent magnetised gas is ubiquitous in astrophysics. In a single system (e.g. stellar and planetary atmospheres, the interstellar medium, etc.), physical quantities such as density, gas pressure, entropy, and temperature can vary by orders of magnitude. We now take our closest star, the Sun, as an example. It has several distinct regions. In the core, energy is generated by nuclear fusion and diffused outward by radiation (gamma and X-rays). When the temperature drops sufficiently for ions to recombine, the opacity of the plasma increases and the radiative diffusion can no longer efficiently transport the energy outward, which causes the plasma to start to simmer. Here, the convection zone begins. The convection zone spans roughly 200 Mm and makes up the outer 30% of the solar interior. The mean free path of the photons increases as the density of the plasma decreases with increasing radius, eventually allowing the photons to escape. This is the optical surface that we can observe: the photosphere. Above the photosphere, in the chromosphere, the thermal energy of the gas no longer dominates the total energy budget; plasma β drops below unity, indicating that magnetic pressure becomes higher than the gas pressure; and the gas flows must follow the magnetic-field configuration. Even higher up, the gas temperature suddenly goes back up to millions of degrees Kelvin. The physical quantities, such as density, gas pressure, and temperature, vary by many orders of magnitude from the deep interior all the way up to the outermost parts of the Sun. However, in the convection zone, radially averaged entropy per unit mass is nearly constant. This suggests using entropy as one of the main quantities in solar and stellar simulations is necessary. Numerically maintaining hydrostatic equilibrium is much easier when a quantity varies very little with depth, unlike temperature or internal or total energy. Existing studies have used entropy as a primary thermodynamic variable – for example, works by Popovas et al. (2018, 2019) and Hotta et al. (2022); however, the partial differential equation (PDE) solvers in these works are either very diffuse in low-velocity motions (i.e. Popovas et al. 2018, 2019) or treat entropy as a conservative quantity and do not account for entropy increase in shocks.

With the advent of exascale computing, we are able to run numerical simulations of increasing complexity, either covering a wider range of physical and temporal scales, or including more physics, such as radiative heat transfer, conductivity, chemistry, non-local thermal equilibrium, the generalised Ohm’s law, and so on. A very popular group of tools to solve the equations of magnetohydrodynamics (MHD) is the Godunov-type Riemann solvers. These solvers are efficient in conditions where detailed turbulent structures ought to be preserved, but due to numerical precision cannot accurately evolve thermal energy when kinetic and magnetic energies dominate the energy budget. This can create negative pressure spots in, for example, solar corona; cause inaccurate convective flow velocities; or create areas with negative temperatures in a giant molecular cloud simulation. This is a well-known problem with HLL- and Roe-type Riemann solvers (e.g. Teyssier 2015, Sect. 6.6). Many attempts have been made to mitigate the issue (e.g. Ryu et al. 1993; Bryan et al. 1995; Ismail & Roe 2009; Winters & Gassner 2016; Herbin et al. 2020; Gallice et al. 2022); these did not necessarily resolve the problems, but they at least make them negligible under certain conditions. One alternative to these fixes would be using logarithmic density and pressure.

In this paper, we present a prescription for an approximate entropy-based Riemann solver, which for simplicity we call HLLS. It is based on the Harten-Lax-van Leer discontinuities (HLLD) approximate Riemann solver in RAMSES code (Fromang et al. 2006); see also Miyoshi & Kusano (2005) for a detailed description of an HLLD solver). The HLLS solver is not entropy conserving – shocks and magneto-acoustic waves generate entropy by converting kinetic and magnetic energy to heat. It is an irreversible process. In a closed system, energy cannot be created or destroyed; it can only change its form. This is a fundamental requirement and is mathematically¹ fulfilled in total energy-based Riemann solvers. As we are converting the solver to an entropy-based one, we must take into account the entropy increase in shocks and magneto-acoustic waves. In principle, the generated entropy S_gen depends on the state of the Riemann fan and is characterised by the change in primitive quantities – density, velocity, and total pressure – when crossing a given discontinuity. However, in practise, computing entropy production via discontinuities is insufficient. Riemann solvers are diffusive due to the use of slope limiters as well as finite numerical precision. Section 2 briefly covers the key aspects of the HLLD Riemann solver, which are relevant to the discussion of the changes and underlying considerations when converting to an entropy-based solver in Sect. 4. We briefly discuss the equation of state (EOS) where entropy is one of the primary variables and how to compute one in Sect. 3. We continue by showing a few numerical tests in Sect. 5. In Sect. 6, we discuss the current limitations of the solver and the applications, where this solver would be the most beneficial. We conclude by discussing our future development work.

2 The total energy-based HLLD Riemann solver

The ideal MHD equations are a classical system of equations used to study the dynamics of conducting fluids. In traditional Riemann solvers, the fluid equations for the total mass density, momentum density, total energy density, and the magnetic field of the system can be written as a system of conservation laws: $$\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left[\begin{array}{c}\rho \\ \rho \mathit{u}\\ E_{\mathrm{t}\mathrm{o}\mathrm{t}}\\ \mathit{B}\end{array}\right]+\mathrm{\nabla }\cdot \left[\begin{array}{c}\rho \mathit{u}\\ \rho (\mathit{u}\otimes \mathit{u})+P_{\mathrm{t}\mathrm{o}\mathrm{t}}-\mathit{B}\otimes \mathit{B}\\ \mathit{u}(E_{\mathrm{t}\mathrm{o}\mathrm{t}}+P_{\mathrm{t}\mathrm{o}\mathrm{t}})-\mathit{B}(\mathit{u}\cdot \mathit{B})\\ \mathit{B}\otimes \mathit{u}-\mathit{u}\otimes \mathit{B}\end{array}\right]=0,$$(1)

with $$\mathrm{\nabla }\cdot \mathit{B}=0,$$(2)

where ρ, u, E_tot, P_tot, and B are density, velocity, total energy, total pressure, and magnetic field, respectively. Total energy is a sum of internal, kinetic, and magnetic energies: $$E_{\mathrm{t}\mathrm{o}\mathrm{t}}=\rho \epsilon +\frac{\rho {\mathit{u}}^2}{2}+\frac{|\mathit{B}{|}^2}{2},$$(3)

where ε is internal energy per unit mass and the total pressure, P_tot, is a sum of gas and magnetic pressures: $$P_{\mathrm{t}\mathrm{o}\mathrm{t}}=P_{\mathrm{g}\mathrm{a}\mathrm{s}}+\frac{|\mathit{B}{|}^2}{2}.$$(4)

The system is closed by an EOS. In this work, we used a simple gamma-law (ideal gas) EOS: $$P=\rho (\gamma -1)\epsilon ,$$(5)

where γ is the adiabatic index.

The beauty of System (1) lies in the total energy – when one type of energy is converted to another, one does not have to explicitly compute the conversion via, for example, Maxwell stress tensors, Lorentz force, diagonal Reynolds stress, and so on. It saves time and retains very high precision levels, at least in theory. In practice, however, there are complications, mostly due to numerical floating point precision. If magnetic or kinetic energy strongly dominates over thermal energy, the precision of the latter drops to an extent that the expressed pressure or temperature can become negative.

3 The entropy-based equation of state

To define the temperature and entropy, we used an ideal gas law: $$P=\frac{\rho k_{B}T}{\mu m_u},$$(6)

where k_B is the Boltzmann constant, μ is the mean molecular weight, T is the gas temperature, and m_u is the atomic mass unit. Entropy per unit mass is defined as $$S=S_0+c_v\ln (P{\rho }^{-\gamma }),$$(7)

where $c_v={\left(\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }T}\right)}_{\rho }$ is the heat capacity at constant volume. Here, we introduce the entropy null-point S₀. It can be considered as an integration constant or the starting point of the summation of the microstates. Having S₀ is not necessary (i.e. S₀ = 0), and mathematically the entropy value can be both positive and negative, but physically the entropy should not go below 0². From a computational point of view, it is best to keep S as close to zero as possible; then, small gradients in entropy retain the precision.

It is really convenient to express Eq. (7) as $$S=S_0+\frac{k_B}{(\gamma -1)\mu m_u}\ln (P{\rho }^{-\gamma })=S_0+\frac{\mathrm{\wp }}{\gamma -1}\ln (P{\rho }^{-\gamma });$$(8)

here, $\mathrm{\wp }=\frac{k_B}{\mu m_u}=R$, where R is the gas constant. Now we can notice that ℘ repeats itself in both the entropy definition and the ideal gas law. Thus, we can use it as a scaling factor for both temperature and entropy, making the modified entropy, $$\mathfrak{S}(\gamma -1)=S_0+ln(P{\rho }^{-\gamma }),$$(9)

a dimensionless quantity, which is also much closer to unity; thus, numerical operations using it are more accurate. Converting to code units is a common practice in numerical simulations.

Entropy in Eq. (7) represents the thermodynamic concept of entropy, which is a measure of the disorder or randomness of a system. For more realistic gases, entropy per mole can be expressed from total partition sums (e.g. Popovas & Jørgensen 2016, Eq. (43)): $$S=R\ln Q_{\mathrm{t}\mathrm{o}\mathrm{t}}+\frac{RT}{Q_{\mathrm{t}\mathrm{o}\mathrm{t}}}\left(\frac{\mathrm{\partial }{Q}_{\mathrm{t}\mathrm{o}\mathrm{t}}}{\mathrm{\partial }T}\right),$$(10)

where Q_tot is the total partition sum.

For an arbitrary tabular EOS, entropy can be derived from other quantities by integrating $$dS=\left[\frac{1}{T}{\left(\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }T}\right)}_{\rho }\right]dT+[\frac{1}{T}{\left(\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }\ln \rho }\right)}_T-\frac{P}{T\rho }]d\ln \rho .$$(11)

The first term in the second brackets on the RHS of Eq. (11) is generally very small, and for ideal gas it vanishes completely, leading to $$dS=\left[\frac{1}{T}{\left(\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }T}\right)}_{\rho }\right]dT-\left[\frac{P}{T\rho }\right]d\ln \rho .$$(12)

It is, however, very important to check that the new entropy-based table satisfies Maxwell’s thermodynamic identity, $${\left(\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }\rho }\right)}_T=\frac{P}{{\rho }^2}-\frac{T}{{\rho }^2}{\left(\frac{\mathrm{\partial }P}{\mathrm{\partial }T}\right)}_{\rho }=\frac{P}{{\rho }^2}\left[{\left(\frac{\mathrm{\partial }\ln P}{\mathrm{\partial }\ln T}\right)}_{\rho }\right],$$(13)

and gives exact results to numerical precision for requested quantities, when compared to the original table. In fact, any EOS must satisfy the thermodynamic identity (Eq. 13).

Lastly, sometimes entropy is adopted to be S = Pρ^−γ (e.g. Ryu et al. 1993; Canivete Cuissa & Teyssier 2022); this is called modified entropy, which should be considered as a different kind of polytropic EOS and does not have a real physical meaning. Although this quantity could be used as a probe for thermal energy, it does not relate to thermodynamic entropy. On the other hand, it works really well in maintaining thermal energy stability in ideal gases, as shown by Ryu et al. (1993). It is not surprising, since the quantity is advected as a passive scalar, that it will be conserved to numerical precision despite not giving accurate values as to where entropy (or thermal energy in this case) should increase.

4 The conversion to an entropy-based HLLD solver

In the DISPATCH implementation (Nordlund et al. 2018) of this new solver, we used the MUSCL-Hancock algorithm with constrained transport (Evans & Hawley 1988, CT) for the induction equation, as well as a positivity-preserving three-dimensional (3D) unsplit TVD slope limiter (see Fromang et al. (2006) for details). For convenience, do not give full details of the Godunov scheme and we use standard notation. However, we summarize the modified numerical method in Appendix A with the changes depicting the conversion from total energy to entropy. For a deeper overview of HLL, HLLC, and HLLD Riemann solvers, we refer the reader to Miyoshi & Kusano (2005) and Toro (2019), and references therein. The main differences compared to typical HLLD solvers are the parts where total energy is replaced with entropy, with additional fluxes computed to calculate the entropy production (details below).

In the system of equations (1), we replace the total energy term with entropy per unit mass: $$\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left[\begin{array}{c}\rho \\ \rho \mathit{u}\\ \rho S\\ \mathit{B}\end{array}\right]+\mathrm{\nabla }\cdot \left[\begin{array}{c}\rho \mathit{u}\\ \rho (\mathit{u}\otimes \mathit{u})+P_{\mathrm{t}\mathrm{o}\mathrm{t}}-\mathit{B}\otimes \mathit{B}\\ \rho \mathit{u}S\\ \mathit{B}\otimes \mathit{u}-\mathit{u}\otimes \mathit{B}\end{array}\right]=\left[\begin{array}{c}0\\ \mathbf{\Phi }\\ ({\mathit{Q}}_{\mathrm{e}\mathrm{x}\mathrm{t}}+{\mathit{Q}}_S)/T\\ 0\end{array}\right].$$(14)

For completeness, we also show the additional source terms here, such as Φ, which is the force per unit volume and includes the force of gravity, Coriolis force, and so on. Q_ext is the external heating per unit volume, which includes Newton cooling, radiative heat transfer, and so on. Q_S is entropy production from converting kinetic and magnetic energy to heat. The source terms from gravity, Coriolis forces, Newton cooling, and so on are added during both the predictor and the corrector steps.

The system in Eq. (14) can be written in vectorial form, similarly to, for example, Londrillo & Del Zanna (2000): $$\frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\mathbf{F}}{\mathrm{\partial }x}+\frac{\mathrm{\partial }\mathbf{G}}{\mathrm{\partial }y}+\frac{\mathrm{\partial }\mathbf{H}}{\mathrm{\partial }z}=\mathrm{\Psi },$$(15)

where $$\mathbf{U}={(\rho ,\rho {\mathrm{u}}_{\mathrm{x}},\rho {\mathrm{u}}_{\mathrm{y}},\rho {\mathrm{u}}_{\mathrm{z}},\rho \mathrm{S},{\mathrm{B}}_{\mathrm{x}},{\mathrm{B}}_{\mathrm{y}},{\mathrm{B}}_{\mathrm{c}})}^{\mathrm{T}},$$(16) $$\mathrm{\Psi }={(0,{\mathrm{\Phi }}_x,{\mathrm{\Phi }}_y,{\mathrm{\Phi }}_z,Q/T+Q_S/T,0,0,0)}^T,$$(17)

and $$\mathbf{F}=\left(\begin{array}{c}\rho u_x\\ \rho u_x^2+P_{tot}-B_x^2\\ \rho u_x{u}_y-B_x{B}_y\\ \rho u_x{u}_z-B_x{B}_z\\ \rho u_{x}S\\ 0\\ u_x{B}_y-u_y{B}_x\\ u_x{B}_z-u_z{B}_x\end{array}\right)$$(18)

is the flux function. The expressions for the terms G and H are completely analogous. These equations have seven eigenvalues, corresponding to four magneto-acoustic waves (two slow and two fast), two Alfvén waves, and an entropy wave: $${\lambda }_{1,7}=u\mp c_f,{\lambda }_{3,5}=u\mp c_s,{\lambda }_{2,6}=u\mp c_a,{\lambda }_4=u,$$(19)

where $$c_{f,s}^2=d\pm \sqrt{d^2-a{B}_x^2/\rho },c_a=\frac{\left|B_x\right|}{\sqrt{\rho }}$$(20)

with $$d=\frac{a^2+|\mathit{B}{|}^2/\rho }{2},$$(21)

and a is the speed of sound. More often than not some eigenvalues in Eq. (19) coincide: $${\lambda }_1\le {\lambda }_2\le {\lambda }_3\le {\lambda }_4\le {\lambda }_5\le {\lambda }_6,$$(22)

depending on the direction and strength of the magnetic field. A solution to the Riemann problem may thus be composed not only of ordinary shock and rarefaction waves, but also compound waves (Brio & Wu 1988; Miyoshi & Kusano 2005). In the set of equations (18), we do not show the additional fluxes that will be needed to compute the entropy generation, which is discussed in detail below.

4.1 Entropy production

The entropy production in shocks from dissipation of kinetic energy has long been a topic of discussion (Salas & Iollo 1995). Kelvin and Rayleigh questioned the validity of shock discontinuity as it violated the conservation of entropy. A very curious result was put forward by Morduchow & Libby (1949) that equilibrium entropy has a maximum inside a Navier-Stokes shock profile, indicating that entropy was decreasing after passing the shock, strengthening the doubts by Stokes, Kelvin, and Rayleigh. Salas & Iollo (1995) explored this curiosity and concluded that because of this phenomenon, the entropy propagation equation cannot be used as a conservation law. Using a single jump condition for the propagation of entropy is not adequate, as going from equations with two jumps, [P] and [u], to an equation with a single jump, [S], information is lost (see Sect. 5 in Salas & Iollo (1995) for details). Shocks in the infinitesimally small area cannot be considered equilibrium, as the adiabatic index changes locally because of the rapid compression. Rankine–Hugoniot jump conditions are only applicable when elements are far enough from the discontinuity, allowing to assume that the gas is in a local equilibrium on each of the sides. The condition also implies the discontinuity is not resolved. Margolin (2017) studied non-equilibrium entropy in a shock and showed that it increases monotonically inside the shock and a certain modification can be applied to the equilibrium formulation so that it follows the non-equilibrium formulation better. Thornber et al. (2008) derived analytical formulae for the rate of increase of entropy at arbitrary jumps in primitive variables for Godunov methods. It is used later for total energy corrections.

Instead of correcting the total energy, we rewrote the Riemann solver to work in terms of entropy as a primary thermodynamic variable. We can split the evolution of entropy into two parts: the advection of the conserved entropy and new entropy production via the dissipation of kinetic and magnetic energy. The conserved entropy is advected through the Riemann fan as a passive scalar (see Appendix A for details). For the second part, we start by looking at the evolution of total, kinetic, magnetic, and thermal energies. The total energy can be split into its constituent thermal (E_th), kinetic (E_kin) and magnetic (E_mag) energies: $$E_{\mathrm{t}\mathrm{o}\mathrm{t}}=E_{\mathrm{t}\mathrm{h}}+E_{\mathrm{k}\mathrm{i}\mathrm{n}}+E_{\mathrm{m}\mathrm{a}\mathrm{g}},$$(23)

while its evolution is thus $$\frac{d{E}_{\mathrm{t}\mathrm{o}\mathrm{t}}}{dt}=\frac{d{E}_{\mathrm{t}\mathrm{h}}}{dt}+\frac{d{E}_{\mathrm{k}\mathrm{i}\mathrm{n}}}{dt}+\frac{d{E}_{\mathrm{m}\mathrm{a}\mathrm{g}}}{dt}.$$(24)

Since the total energy is conserved, any local change in one form of energy must be accounted for by corresponding fluxes or conversions into other forms. To properly track entropy production, we examine the contributions from kinetic, magnetic, and thermal energies separately: $$\frac{d{E}_{\mathrm{k}\mathrm{i}\mathrm{n}}}{dt}=-\mathbf{\nabla }\cdot {\mathbf{F}}_{\mathrm{k}\mathrm{i}\mathrm{n}}+W_{\mathrm{g}\mathrm{a}\mathrm{s}}-Q_{\mathrm{k}\mathrm{i}\mathrm{n}}+{\mathrm{\Theta }}_{\mathrm{k}\mathrm{i}\mathrm{n}},$$(25) $$\frac{d{E}_{\mathrm{t}\mathrm{h}}}{dt}=-\mathbf{\nabla }\cdot {\mathbf{F}}_{\mathrm{t}\mathrm{h}}-W_{\mathrm{g}\mathrm{a}\mathrm{s}}+Q_{\mathrm{k}\mathrm{i}\mathrm{n}}+{\mathrm{\Theta }}_{\mathrm{t}\mathrm{h}},$$(26) $$\frac{d{E}_{\mathrm{m}\mathrm{a}\mathrm{g}}}{dt}=-\mathbf{\nabla }\cdot {\mathbf{F}}_{\mathrm{m}\mathrm{a}\mathrm{g}}-{\mathrm{\Theta }}_{\mathrm{k}\mathrm{i}\mathrm{n}}-{\mathrm{\Theta }}_{\mathrm{t}\mathrm{h}},$$(27)

where ∇ · F_kin, ∇ · F_th, and ∇ · F_mag are the kinetic, thermal, and magnetic energy fluxes, respectively; $\frac{d{E}_{\text{kin }}}{dt},\frac{d{E}_{\text{th }}}{dt}$, and $\frac{d{E}_{\text{mag }}}{dt}$ are the change of kinetic, thermal, and magnetic energies in a given cell in a time-step dt. Q_kin is heating from kinetic energy dissipation, and Θ_kin and Θ_th are conversion of magnetic energy to kinetic and thermal energies, respectively (e.g. reconnection, Ohmic dissipation, etc.). Lastly, W_gas = P(∇ · u) is the pressure work. From Eqs. (25)–(27), we can derive the heating from the irreversible energy dissipation Q_S = Q_kin + Θ_th: $${\mathbf{Q}}_{\mathrm{S}}=-\mathrm{\nabla }\cdot {\mathbf{F}}_{\mathrm{k}\mathrm{i}\mathrm{n}}-\mathrm{\nabla }\cdot {\mathbf{F}}_{\mathrm{m}\mathrm{a}\mathrm{g}}-\frac{d{E}_{\mathrm{k}\mathrm{i}\mathrm{n}}}{dt}-\frac{d{E}_{\mathrm{m}\mathrm{a}\mathrm{g}}}{dt}+W_{\mathrm{g}\mathrm{a}\mathrm{s}}.$$(28)

Therefore, to compute Q_S we need additional quantities from the corrector step: the kinetic and magnetic energy fluxes, as well as the normal velocities at the cell faces to compute ∇ · u. These velocities can be extracted from the corrector step by passing the normal velocity along with the fluxes (see Appendix A for the complete description of the corrector step). To accurately compute Q_S, a consistent evaluation of the kinetic and magnetic energy time derivatives is needed.

The time derivatives of kinetic and magnetic energy, $\frac{d{E}_{\text{kin }}}{dt}$ and $\frac{d{E}_{\text{mag }}}{dt}$, are computed by evaluating these quantities at both the old (at t) and new (at t + 1) time steps. Specifically, E_kin and E_mag are first computed from the conserved variables at the beginning of the time step. Then, after updating the density, momentum, and magnetic fields via flux divergences, the new values of E_kin and E_mag are computed. The time derivatives are then approximated as simple finite differences between these two states. This approach ensures consistency with the numerical fluxes and avoids introducing additional diffusive or dispersive errors. Furthermore, this discretisation is fully consistent with the MUSCL-Hancock predictor-corrector scheme used in the solver. If a different time-integration scheme, such as a higher order Runge–Kutta method, were employed, the computation of $\frac{d{E}_{\text{kin }}}{dt}$ and $\frac{d{E}_{\text{mag }}}{dt}$ would need to be adjusted accordingly to maintain consistency with the chosen temporal discretisation.

This entropy-generation calculation recipe fully accounts for entropy production in shocks and shearing, as well as for numerical dissipation of kinetic and magnetic energy via slope limiters, for example, thus fully closing the system. Lastly, to fully adhere to the second law of thermodynamics, we require Q_S to be positive-definite: $$\frac{\mathrm{\partial }(\rho S)}{\mathrm{\partial }t}=-\mathbf{\nabla }\cdot (\rho \mathbf{u}S)+max(0,\frac{{\mathit{Q}}_{\mathbf{S}}}{T}).$$(29)

4.2 Summary

To reformulate a total-energy-based HLLD Riemann solver, we need to replace the total energy with entropy in both predictor and corrector steps, if the Godunov method is used. Here are the steps:

1. Replace thermal-energy-based EOSs (Eq. (5)) with entropy-based EOSs (Eq. (8)).

2. In the predictor step, it is sufficient to advect S as a passive scalar to obtain the source term σ, (e.g. for the right state): $${\sigma }_S=\frac{\mathrm{\Delta }t}{2}(-\mathit{u}\frac{\mathrm{\partial }s}{\mathrm{\partial }t})\pm \frac{\mathrm{\Delta }\xi }{2}(-\mathit{u}\frac{\mathrm{\partial }s}{\mathrm{\partial }\xi }),$$(30)

   where ξ = [x, y, z] is the spatial dimension.

3. Compute the Godunov flux for entropy instead of total energy.

4. Compute the Godunov flux for kinetic and magnetic energies (see Appendix A).

5. Compute Q_S using Eq. (28).

6. Compute the $S_{\text{gen }}=max(0,\frac{Q_{\mathrm{S}}}{T})$.

7. Apply Godunov flux, heating $\frac{Q_{\text{ext }}}{T}$ and S_gen to obtain S^t+1.

In the next section, we show a number of tests we put this new solver through.

5 The numerical tests

To check the validity of the HLLS solver, we conducted a series of experiments. They are in 1D, 2D, and 3D, and both hydrodynamic and MHD. We ran 1D tests in all principal directions and 2D tests in all three planes to make sure the solver is well balanced in all dimensions. We note that some tests were run on low resolution intentionally, as resolution-starved experiments may reveal problems that are inversely proportional to resolution. The experiments were run in single floating-point precision.

Fig. 1 Convergence in norm of L1 error vector for contact, fast, slow, and Alfvén waves after propagating a distance of one wavelength. The convergence rate P is given in the legend.

5.1 Linear wave convergence

This test shows that the HLLS formulation converges with second-order accuracy for linear amplitude waves. The test is nearly identical to that of Sect. 8.2 in Stone et al. (2008), with a uniform medium where the generic quantities are $(\rho ,\gamma ,P,\mathbf{u},B_x,B_y,B_z)=[1,5/3,3/5,0,1,\sqrt{2},1/2]$. The exact eigenfunctions for fast and slow magnetosonic, Alfvén, and contact waves for the background state are given in Gardiner & Stone (2005). Since we usually operate in single floating-point precision, the original wave amplitude of A = 10⁻⁶ is too close to the round-off error and is impractical to conduct an informative test. Thus, the only difference compared to Stone et al. (2008) is that we set the wave amplitude value to be A = 10⁻³. The resulting waves are still weak and do not steepen into shocks. Figure 1 shows the norm of the L₁ error vector for each wave family as a function of the numerical resolution from 48 to 1024 cells. The figure shows that we achieve second-order accuracy in both space and time for smooth solutions.

5.2 Entropy wave

Riemann solvers are generally very good at preserving shocks. However, very slowly moving waves are sometimes more difficult. To test the diffusion and dispersion error, we launched a very slowly moving discontinuity (entropy wave) in one dimension. Assuming the experiment is carried out in the x direction, the generic quantities are (ρ, P) = [0.9, 1.0] when x < 0 and [1.1, 1.0] when x ≥ 0. All other quantities are set to 0, with $\gamma =\frac{5}{3}$ and $S=ln(\frac{P}{{\rho }^{\gamma }})/(\gamma -1)$ in the entire experiment. Boundaries are periodic. We used a rather low resolution of 100 cells in the direction of the wave propagation to exaggerate the dispersion error. The velocity u_x of the whole box was set as constant to [1.0, 0.1, 0.01, 0.001, 0.001], and the runtime of the experiment was set to, respectively, t_end = [1, 10, 100, 1000, 1000] time units. This means the number of updates increased tenfold with each run. If the solver is diffusive, the wave form should smear out. Figure 2 shows the test results at the end time, when the wave crosses one period. Naturally, the corners are slightly rounded from slope limiters, but we can clearly see that with lower velocities the increase in diffusion is negligible. The diffusion is resolution-sensitive and depends on the slope limiters. The black curve is at time=0, the blue one is at u_x = 1, the red one is at u_x = 0.1, the dashed green one is at u_x = 0.01, and the dotted magenta one is at u_x = 0.001. At the lowest velocity, we observe noticeable distortion. This is expected, as many Riemann solvers that are not asymptotic-preserving become excessively diffusive at low Mach numbers. The issue arises because numerical dissipation is often linked to jumps in the velocity component that is normal to the cell interfaces, which do not scale correctly with Mach number. As a result, even small spurious variations in this velocity component can lead to disproportionately large diffusion, particularly when M ≤ 0.2.

Fig. 2 Entropy wave test. Density profiles at starting time (black) and for different wave speeds ux: blue – 1.0, red – 0.1, green dashed – 0.01, magenta dotted – 0.001. The overlap of wave profiles is perfect, with expected slopes present.

5.3 Shu & Osher shock tube

This test is a 1D Mach=3 shock interacting with sine waves in density (Shu & Osher 1989). It tests the solver’s ability to capture both shocks and small-scale smooth flow. The computational domain is nine length units long and is split into two regions with different conditions in each of them. Assuming the experiment is carried out in the x direction, the quantities are (ρ, u_x, P) = [3.857143, 2.629369, 10.33333] when x < −4 and [1 + 0.2 sin (5x), 0,1] when x ≥ −4. All the other quantities are set to 0, with γ = 1.4 and $S=ln(\frac{P}{{\rho }^{\gamma }})/(\gamma -1)$ in the whole experiment. Figure 3 shows the test results at time = 1.8 time units (red curve). The reference run is done with the HLLD solver (black curve, perfectly covered by the red curve). Both the test and the reference runs had 1500 cells. We additionally show an HLLS run with 500 cells (dots).

5.4 Brio & Wu shock tube

This is a classical test of an MHD shock tube, described by Brio & Wu (1988, Sect. V), where the right and left states are initialised to different values. The left state is initialised as (ρ, u_x, u_y, u_z, B_y, B_z, P) = [1, 0, 0, 1, 0, 1], and the right state is [0.125, 0, 0, −1, 0, 0.1]. B_x = 0.75 and γ = 2. This test shows whether the solver can accurately represent the shocks, rarefaction waves, compound structures, and contact discontinuities. Figure 4 shows the test results at time = 1.0. The black curve is for a reference HLLD run with 1200 cells. The red curve is for the HLLS run. The blue dots show the HLLS run with 400 cells. Both high- and low-resolution runs overlap each other nearly perfectly and follow the profiles in the literature very accurately. From left to right we can identify a fast rarefaction wave, compound wave, contact discontinuity, slow shock, and a fast rarefaction wave again.

Fig. 3 Shu & Osher shock-tube test. Density profile at time=1.8; black: reference HLLD solver with 1500 cells; red: HLLS solver with 1500 cells; blue dots: HLLS solver with 500 cells.

Fig. 4 Brio & Wu shock-tube test. Density profile (left) and y component of magnetic field (right) at time=1.0; black: reference HLLD solver with 1200 cells; red: HLLS solver with 1200 cells; blue dots: HLLS solver with 400 cells. Black and red overlap each other nearly perfectly.

5.5 Kelvin–Helmholtz instability

Kelvin-Helmholtz instability occurs when velocity shear is present within a continuous fluid or across fluid boundaries. We conducted the test as described in McNally et al. (2012). For convenience, here we summarize the setup. The domain is one unit by one unit in the x – and y-directions, with resolutions of 256 × 256 cells. All boundaries are periodic. The density is given by $$\rho =\{\begin{array}{ll}{\rho }_1-{\rho }_m{e}^{\frac{y-1/4}{L}}& \text{ if }1/4>\mathrm{y}\ge 0\\ {\rho }_2+{\rho }_m{e}^{\frac{1/4-y}{L}}& \text{ if }1/2>\mathrm{y}\ge 1/4\\ {\rho }_2+{\rho }_m{e}^{-\frac{3/4-y}{L}}& \text{ if }3/4>\mathrm{y}\ge 1/2\\ {\rho }_1-{\rho }_m{e}^{-\frac{y-3/4}{L}}& \text{ if }1>\mathrm{y}\ge 3/4\end{array},$$(31)

where ρ_m = (ρ₁ − ρ₂)/2, ρ₁ = 1.0, ρ₂ = 2.0 and L = 0.025. The x-direction velocity is given by $$u_x=\{\begin{array}{ll}{u}_1-u_m{e}^{\frac{y-1/4}{L}}& \text{ if }1/4>\mathrm{y}\ge 0\\ u_2+u_m{e}^{\frac{1/4-y}{L}}& \text{ if }1/2>\mathrm{y}\ge 1/4\\ u_2+u_m{e}^{-\frac{3/4-y}{L}}& \text{ if }3/4>\mathrm{y}\ge 1/2\\ u_1-u_m{e}^{-\frac{y-3/4}{L}}& \text{ if }1>\mathrm{y}\ge 3/4\end{array},$$(32)

where u_m = (u₁ − u₂)/2, u₁ = 0.5 and u₂ = −0.5. The background shear is perturbed by adding velocity in the y–direction: $$u_y=0.01\sin (4\pi x).$$(33)

An ideal EOS with γ = 5/3 is used. The initial gas pressure P_gas = 2.5. We ran the simulation until time t = 10. To make an exact comparison to codes in McNally et al. (2012), we show the gas density at t = 1.5 and the maximal value of vertical kinetic energy evolution in time in Figure 5. We note that we used 256 × 256 resolution, but the result is directly comparable to the 512² and 4096² runs in McNally et al. (2012). In the right panel of the figure, we show the maximum value of vertical kinetic energy in the simulation for three resolutions. From this figure, it can be seen that primary instability begins at exactly the same time. Secondary instabilities occur around the same time as well, indicating that we obtain a converged solution even at the low resolution of 256 × 256.

Fig. 5 HD Kelvin-Helmholtz instability. Density profile (left) and time evolution of kinetic energy (right).

5.6 Rayleigh–Taylor instability

Another classic test of a code’s ability to handle subsonic perturbations is the Rayleigh–Taylor instability, which has been described in a number of studies (e.g. Stone & Gardiner 2007; Abel 2011; Hopkins 2015). In this test, a layer of heavier fluid is placed on top of a layer of lighter fluid. With gravitational source term added to the forces, we test two things: whether the explicit addition of force is correct, and the solver’s ability to preserve instabilities while keeping things symmetric at the right time; that is, during the linear phase. During the nonlinear phase, in a sufficiently high resolution the symmetry is expected to break by construction. We used the initial setup similar to those of Abel (2011) and Hopkins (2015). Here, we recap the setup for convenience. In a two-dimensional domain with 0 ≤ x ≤ 0.5 (128 cells for low-resolution runs and 768 for high-resolution runs) and 0 ≤ y ≤ 1 (256 and 1536 cells for low- and high-resolution runs, respectively), we used periodic boundary conditions in the x direction and constant boundary conditions in the y direction. In this test, we used γ = 1.4, and the density profile was initialised as ρ(y) = ρ₁ + (ρ₂ − ρ₁)/(1 + e^{−(y−0.5)/Δ}), where ρ₁ = 1 and ρ₂ = 2 are the density below and above the contact discontinuity (y_c, respectively, with Δ = 0.025 being its width). The pressure gradient is in hydrostatic equilibrium with a uniform gravitational acceleration of g = −0.5 in the y direction, P = ρ₂/γ + g ρ(y) (y − y_c), and then entropy S = (ln(P) − ln(ρ)^γ)/(γ − 1). An initial y-velocity perturbation v_y = δv_y(1 + cos(8π(x + 0.25)))(1 + cos (5π(y − 0.5))) with δv_y = 0.025 is applied in the 0.3 ≤ y ≤ 0.7 range (Hopkins 2015).

Figure 6 shows the evolution of the instability at different times. The initial velocity grows, and the heavier fluid starts to sink. We note the single-cell resolution of contact discontinuities and mixing. Both blobs are perfectly symmetric during the linear phase. Soon enough, Kelvin-Helmholtz instabilities develop at the shear surface between the fluids, and the symmetry is harder to maintain due to low relative velocities. However, the low-resolution run maintained perfect symmetry throughout the whole simulation. Figure 7 shows the high-resolution run. Secondary instabilities are more pronounced, and thus the symmetry is harder to maintain. The figure can be directly compared to Fig. 22 in Hopkins (2015). Our high-resolution (768 × 1536) run shows very similar features to those of Hopkins (2015), although we maintain much better symmetries until the blobs reach the bottom.

Fig. 6 Rayleigh–Taylor instability in low (128 × 256) resolution. We plot density at different times.

Fig. 7 Same as Fig. 6, but in 768 × 1536 resolution.

5.7 ‘Hot-bubble’ experiment

To evaluate the entropy-preserving properties of our numerical scheme in a strongly stratified, low Mach-number regime, we performed the hot-bubble test following the setup of Edelmann et al. (2021). This experiment simulates the buoyant rise of a localised entropy perturbation embedded in an isentropic but highly stratified layer. In stellar interiors, entropy fluctuations drive convection; however, outside the heating or cooling layers these fluctuations should remain conserved, except for mixing effects. Maintaining accurate entropy advection in such conditions is usually numerically challenging due to the large variations in pressure and density. The experimental setup is identical to the original, with resolution of 128 × 192 cells in the x- and y-directions. The gravitational acceleration is then $$g_z=g_0\sin (k_{z}z),$$(34)

where g₀ = −1.099044 × 10⁵ cm and k_z = 2π/z_max, z_max = 1.5 × 10⁶ cm, while x_max is 1 × 10⁶ cm. γ = 5/3, at z = 0 P₀ = 10⁶ Ba, with temperature T₀ = 300 K. The stratification is isentropic; that is, $$\rho (z)={\left(\frac{P(z)}{A}\right)}^{1/\gamma },$$(35)

where A = A₀ = const everywhere outside of the bubble. As the A₀ value is not given in Edelmann et al. (2021), we used the ideal gas EOS to obtain ρ₀ = P/R_gasT = 4.00906 × 10⁻⁵ and consequently A0 = 2.129425 × 10¹³. A = P/ρ^γ is not a true entropy expression, as discussed in Sect. 3, but a rather good probe for it. Since ideal gas entropy can be simplified to S = ln(P/ρ^γ)/(γ − 1), we can express S = S₀ = ln(A)/(γ − 1) = 46.03419. The hydrostatic pressure stratification is given by $$P(z)={[P_0^{1-\frac{1}{\gamma }}+(1-\frac{1}{\gamma })\frac{g_0}{A_0^{\frac{1}{\gamma }}{k}_z}[1-\cos \left(k_{z}z\right)]]}^{\frac{\gamma }{\gamma -1}},$$(36)

which is not perturbed and has a low-to-high ratio of 100. The perturbation is only applied for density via $$A=A_0[1+{\left(\frac{\mathrm{\Delta }A}{A}\right)}_{t=0}\cos {\left(\frac{\pi r}{2r_0}\right)}^2],$$(37)

where (ΔA/A)_t=0 = 10⁻³ is the bubble’s amplitude and r = [(x − x₀)² + (y − y₀)²] is the distance from the bubble’s center. The radius of the bubble is r₀ = 1.25 × 10⁵ cm. This low perturbation makes the bubble rise at low Mach numbers: a few times 10⁻² (Edelmann et al. 2021). To improve the numerical precision, we rescale the length scale, l_sc = 10⁵, and density scale, d_sc = 4 × 10⁻⁵ leaving the timescale intact. This results in velocities and density, and entropy closer to unity, which in turn leads to the highest numerical accuracy. Since we considered the normalised fluctuations in A, ΔA/A, the scaling does not change the result. The simulation was run from time t = 0 to t = 300 s. Two time slices at t=150 and 300 s are shown in Fig. 8. This figure can be directly compared to Fig. 6 in Edelmann et al. (2021). The HLLS solver performs almost as well as a well-balanced scheme and noticeably better than simulations with energy as the primary thermodynamic variable without any well-balanced schemes. We note that the main challenge for us is the perturbation being very close to the numerical noise limit, as we normally operate with single precision. The adherence to the second law of thermodynamics is ensured to the numerical precision via Eq. (29); thus, the negative entropy and, consequently, A values are from numerical round-off errors.

Fig. 8 ‘Hot-bubble’ experiment at time=150 and 300 s. Both panels show fluctuations of ΔA/A, with the min/max values given in the floating box.

5.8 MHD blast

This is a very popular test, and various papers have presented results with slightly different problem setups. It is a very good test of the code’s ability to handle the evolution of strong MHD waves and look for directional biases. We chose to follow the setups of Ramsey et al. (2012), Clarke (2010), and Stone et al. (2008). A rectangular domain with −0.5 ≤ [x, y] ≤ 0.5, 512 × 512 cell resolution was used. All boundaries are periodic. The MHD quantities are set to (ρ, u, B_x, B_y, B_z) = $(1,0,5\sqrt{2},5\sqrt{2},0)$. The ambient pressure P_gas = 1 with P_gas = 100 within the radius r = 0.125. Figure 9 shows the density and magnetic pressure in the experiment at time t = 0.021 and can be directly compared to Ramsey et al. (2012), Clarke (2010), and Stone et al. (2008). The shock front is sharp, a fast magneto-acoustic wave is travelling at the correct speed, and the blast is symmetric.

Fig. 9 MHD blast experiment. Density (left) and magnetic pressure (right) at time=0.021. We note the sharp shock-fronts.

5.9 Orszag-Tang vortex

This test has become a staple for MHD codes. The initial setup is identical to those of Stone et al. (2008), and Ramsey et al. (2012). A rectangular domain with −0.5 ≤ [x, y] ≤ 0.5, 1024 × 1024 cells for a high-resolution run is used. All boundaries are periodic. Initially, the pressure and density were constant: P_gas = 5/12π and ρ = 25/36π. The ratio of specific heats γ = 5/3. The initial velocity (u_x, u_y, u_z) = (− sin (2πy), sin (2πx), 0), and the magnetic field was set through the vector potential $$A_z=\frac{B_0}{4\pi }\cos (4\pi x)+\frac{B_0}{2\pi }\cos (4\pi y),$$

where $B_0=1/\sqrt{4\pi }$. Figures 10 and 11 show density, entropy per unit mass, and magnetic pressure at two times: t = 0.5 and t = 0.75. The t = 0.5 is the typical time this test is shown in literature. When compared to Stone et al. (2008) and Ramsey et al. (2012), we can recognize similar features. We note the very sharp features and perfect symmetry between sides. When advancing the simulation further, the vortex starts producing plasmoids. Here, it is very difficult to maintain symmetry, but the lower panels show that plasmoids are released rather symmetrically.

This experiment is very effective in testing the solver stability. In high resolution, discontinuities and rarefactions are more severe, and Riemann solvers can crash from negative pressure or thermal energy values if they cannot deal with them. Since we have an HLLD solver that can also use both total and thermal energies as an available primary thermodynamic variable, we can investigate the conservation of different quantities in different formulations. Figure 12 shows the time evolution of horizontally averaged kinetic (top left), magnetic (bottom right), internal (top middle), and total (top middle) energies, as well as density (bottom left) and entropy (bottom middle). The resolution was identical in all the runs (256 × 256). HLLS represents the run with the solver (presented in this paper), and HLLS-2 allows Q_S to be negative, breaking the second law of thermodynamics, but more accurately following the energy dissipation and numerical dispersion (see discussion in Sect. 6.1); HLLD is a classical HLLD solver where total energy is conserved and magnetic fields are updated using the CT method, but no adjustments are made to the total energy. While HLLD CT is a modification of said solver, where the magnetic-energy contribution in the total energy is updated with the values from CT solver, making it more stable in areas where magnetic energy strongly dominates over kinetic energy. This correction makes the solver not conserve the total energy in favor of conserving all three energies separately.

In the figure, we can see that the HLLS solver does not conserve the total energy, by construction. However, if we allow Q_S < 0, it tracks the total energy-conserving HLLD solver rather closely, indicating where the second law of thermodynamics is broken in the HLLD solver. The primary reason for this non-conservation is the discrepancy of magnetic fields between the 1D Riemann solver and the CT 2D solver. The former is more diffusive and assumes a normal magnetic field is constant across the interface, while the latter is much less diffusive and preserves ∇ · B = 0 better. Secondary effects stem from numerical dispersion.

The discrepancy between the solvers diminishes with higher resolution. In realistic simulations, open boundary conditions, external forces, and a higher desired resolution may deem the discrepancies in total energy negligible.

Fig. 10 Orszag-Tang vortex test at time = 0.5. Density (left), entropy per unit mass (middle), and magnetic pressure (right).

Fig. 11 Orszag-Tang vortex test at time = 0.75. Density (left), entropy per unit mass (middle), and magnetic pressure (right). We note the plasmoids forming.

Fig. 12 Current-sheet test at different times. Density (top) and magnetic pressure (bottom).

Fig. 13 Current-sheet test at different times. Density (top) and magnetic pressure (bottom).

5.10 Current sheet

The current sheet test problem is designed to determine what an algorithm will do with a perturbed current sheet. Although an analytic solution for this test problem is not available, it is still suitable for testing the robustness of the algorithm. The experimental setup is similar to that of Hawley & Stone (1995). The test is run in two dimensions using a periodic square grid with −0.5 ≤ [x, y] ≤ 0.5 and 256 × 256 cell resolution. There is a uniform magnetic field that discontinuously reverses direction at a certain point. In the whole box (ρ, P_gas) = (1, β/2, where β is an input parameter. We set B_y: $$B_y=\{\begin{array}{ll}1/\sqrt{(4\pi )}& \text{ if }|\mathrm{x}|>0.25\\ -1/\sqrt{(4\pi )}& \text{ if }|\mathrm{x}|\le 0.25\end{array}.$$(38)

For |x| > 0.25, velocities (u_x, u_y) = (A sin (2πy), 0), where A is an amplitude. We ran the test with a range of β and A values, but here we present it only with the “standard” setup, where β = A = 0.1. In Fig. 13 we can see the density and magnetic pressure at different times. Initially linearly polarised Alfvén waves propagate along the field in the y-direction and quickly start generating magneto-acoustic waves since the magnetic pressure does not remain constant. As there are two current sheets in the setup (at x = ±0.25), reconnection inevitably occurs. If β < 1, this reconnection drives strong over-pressurised regions that launch magneto-acoustic waves transverse to the field. Moreover, as reconnection changes the topology of the field lines, magnetic islands form, grow, and merge. The point of the test is to make sure the algorithm can follow this evolution for as long as possible without crashing. We ran tests for 0.1 ≤ (β, A) ≤ 10, and the HLLS solver did not encounter any numerical problems.

Fig. 14 Magnetic-loop advection experiment. The initial magnetic pressure (left) and after the loop has been advected around the grid twice.

5.11 Magnetic-field loop advection

This is a very powerful test to check whether the scheme preserves ∇ · B = 0. The experiment is similar to those of Tóth & Odstrčil (1996) and Stone et al. (2008) with a two-dimensional domain, 0 ≤ x ≤ 2 (128 cells), and 0 ≤ y ≤ 1 (64 cells). Boundaries are periodic in the whole domain. In this test, γ = 5/3 and both density and gas pressure are constant throughout the box: ρ = P_gas = 1.0. The magnetic field is initialised using a vector potential A_z = max [A * (r₀ − r), 0], with A = 10⁻³ and r₀ = 0.3, and r is the radial distance from the domain center. The flow velocity (u_x, u_y) = (2, 1); thus, the problem is essentially an advection test for the vector potential.

Although this test does not test the HLLS solver directly, the CT scheme is nevertheless an integral part of the solver. Moreover, we found this test to be a particularly useful tool to find errors in the magnetic-field predictor step. Figure 14 shows the initial magnetic pressure and that after the magnetic loop has been advected twice through the domain. The shape of the loop is rather well preserved.

Fig. 15 Gresho vortex test. The initial setup (left) and after two rotation periods with different Mach numbers: Ma = 0.1, 0.01, 0.001. The rightmost panel is for a higher resolution run (see text for more details).

Fig. 16 Gas density (left) and pressure (right) in the Gresho vortex test with Ma = 0.01. We note that the quantities are within the noise limit.

5.12 Gresho vortex

The Gresho vortex is a time-independent rotation pattern. Angular velocity depends only on the radius, and centrifugal force is balanced by the pressure gradient. We used the slightly modified initial condition, which permits the variation of the Mach number; it can be found in Happenhofer et al. (2013) and Grimm-Strele et al. (2014). For convenience, the setup is summarised here. The simulation is in two dimensions − 0 ≤ x ≤ 1 (48 cells) and 0 ≤ y ≤ 1 (48 cells) – with periodic boundary conditions everywhere. The low resolution is deliberate; from our tests, we see that with higher resolution (e.g. Grimm-Strele et al. 2014) the experiment becomes much less challenging to the solver. In the whole domain, ρ = 1, γ = 5/3, $$u_{\varphi }=\{\begin{array}{ll}5r& \text{ if }0\le \mathrm{r}<0.2\\ 2-5r& \text{ if }0.2\le \mathrm{r}\le 0.4,and\\ 0& \text{ if }0.4<\mathrm{r}\end{array}$$(39) $$P_{\text{gas }}=\{\begin{array}{l}{P}_0+\frac{25}{2}{r}^2\text{ if }0\le \mathrm{r}<0.2\\ P_0+\frac{25}{2}{r}^2+4(1-5r-\ln (0.2r)\text{ if }0.2\le \mathrm{r}\le 0.4,\\ P_0-2+4\ln (2)\text{ if }0.4<\mathrm{r}\end{array}$$(40)

where $r=\sqrt{(x^2+y^2)}$ is the radial distance from the center of the domain, u_ϕ is the angular velocity in terms of the polar angle ϕ = atan2(y, x), and $$P_0=\frac{\rho }{\gamma {\mathrm{M}\mathrm{a}}_{\mathrm{r}\mathrm{e}\mathrm{f}}^2},$$(41)

with Ma_ref being a reference Mach number, which is the highest Mach number in the resulting flow. Four runs were executed, with Ma_ref = [0.1, 0.01, 0.001]. The last Ma_ref is repeated in two runs at the nominal resolution and a higher resolution, with 128 × 128 cells. Figure 15 shows the results. We deliberately did not run the experiment with higher Mach numbers, as those are just too easy to maintain and the final result looks identical to the initial condition. In the figure, we can see that the vortex is maintained very well down to Ma=0.01. With Ma=0.001, we have the expected outcome; that is, the Godunov-type Riemann solvers deal with very low Mach numbers poorly, even though the Ma=0.01 result is still very good and especially since we used single floating point precision. Figure 16 shows the gas density and pressure with Ma=0.01. It can be clearly seen that the vortex itself is barely distinguishable inside the noise of the box. On the other hand, with increased resolution, the result with Ma=0.001 (as expected) has significantly improved, and the vortex can be once again identified.

5.13 Magnetic rotor

The magnetic rotor tests the propagation of strong torsional Alfén waves as well as solver behavior in moderately high Mach numbers (Ma = 20). A dense disc of fluid rotates within a static fluid background, with a gradual velocity-tapering layer between the disc edge and the ambient fluid. An initially uniform magnetic field is present, which is twisted with the disc rotation. The magnetic field is strong enough that it wraps around the rotor, diminishing the rotor’s angular momentum. The increased magnetic pressure around the rotor compresses the fluid in the rotor, giving it an oblong shape. The experimental setup is introduced by Balsara & Spicer (1999), but we used a more stringent variant of it from Guillet et al. (2019), which is summarised below.

The simulation is in a two-dimensional square, 0 ≤ x, y ≤ 1 (512 × 512 cells), with periodic boundary conditions everywhere. The gas pressure and magnetic fields are uniform in the whole domain, with P_gas = 1, γ = 1.4 and $\mathbf{B}=(5/\sqrt{4\pi },0,0)$. The gas density is $$\rho =\{\begin{array}{ll}10& \text{ if }\mathrm{r}<{\mathrm{r}}_0\\ 1+9f& \text{ if }{\mathrm{r}}_0\le \mathrm{r}\le {\mathrm{r}}_1\\ 1& \text{ if }{\mathrm{r}}_1<\mathrm{r}\end{array},$$(42)

where r₀ = 0.1, r₁ = 0.115, r is the radial distance from the center of the box c, and f = (r₁ − r)/(r₁ − r₀) is the tapering function for the taper region between the disc and the background. The velocities are $$(u_x,u_y)=\{\begin{array}{ll}(\frac{v_0(c-y)}{r_0},\frac{v_0(x-c)}{r_0})& \text{ if }\mathrm{r}<{\mathrm{r}}_0\\ (f\frac{v_0(c-y)}{r_0},f\frac{v_0(x-c)}{r_0})& \text{ if }{\mathrm{r}}_0\le \mathrm{r}\le {\mathrm{r}}_1\\ (0,0)& \text{ if }{\mathrm{r}}_1<\mathrm{r}\end{array},$$(43)

with v₀ = 2. The experiment was run until time t = 0.15, by which point the torsional Alfvén waves have almost reached the boundary. Figure 17 shows the density ρ, magnetic pressure P_B, Mach number, and normalised magnetic-field divergence $\frac{(\mathrm{\nabla }\cdot \mathbf{B})\mathrm{\Delta }x}{|\mathbf{B}|}$. We can note the very sharp details with no distortions outside the now almond-shaped disc. In the compressed areas, the Mach number becomes very high, but it does not pose any issues. The magnetic field is divergence-free within the noise level.

Fig. 17 Magnetic rotor test. From left to right: Density, magnetic pressure, Mach number, and $\frac{\mathrm{\nabla }\cdot \mathbf{B}}{|\mathbf{B}|}$ (see text for more details).

6 Discussion and conclusions

In this paper, we present a new approximate entropy-based HLLD Riemann solver. It works well in both sub- and supersonic regimes and preserves positive temperature and gas pressure. The numerical tests are very encouraging and indicate that the HLLS solver can be readily used in a wide range of physical conditions and experimental setups.

6.1 Considerations about Q_S being positive-definite

Even though we require Q_S to be positive-definite, in reality some numerical dispersion and diffusion might artificially increase kinetic or magnetic energy, which in turn breaks total energy conservation if a corresponding decrease in thermal energy does not occur. Indeed, this is one of the main reasons many Riemann solvers fail in conditions where kinetic or magnetic energy strongly dominate over thermal energy: the numerical precision of the latter is reduced if it is several orders of magnitude smaller than the other total energy components, and it can lead to negative thermal energy. In the Orszag-Tang test, described in the Sect. 5.9, we did check the difference; that is, whether the requirement of Q_S being positive-definite is upheld or not. Until time ≈0.3, the difference is non-existent (indicating that numerical dispersion would not artificially decrease the thermal energy in a total energy solver); however, later we do observe considerable differences. This difference could sometimes be decisive of whether the experiment crashes or not when total or thermal energy is used as a primary thermodynamic variable.

6.2 Very low-Mach-number regimes

The numerical tests show that the HLLS solver can easily handle Mach numbers as low as 0.01, and with lower values it becomes rather diffusive. Of course, the perturbations at such Mach numbers are on the order of numerical precision (see the Gresho vortex in Sect. 5.12 for details) and can be absolutely indistinguishable from the background. This is very encouraging, as Godunov- and Roe-type Riemann solvers with singular precision usually struggle with Mach numbers below 0.1. Various good attempts have been made to modify solvers to reach very low Mach numbers, including those adding correction terms to star states – for example, Shima & Kitamura (2011), Dellacherie et al. (2016), Minoshima & Miyoshi (2021), and Sheng Chen et al. (2022) – or by using well-balanced schemes. In the latter, a hydrostatic equilibrium is imposed directly in the set of dynamic equations, separating primitive variables into equilibrium (stationary) states and dynamical perturbations, as was done in, for example, Greenberg & Leroux (1996), Hotta et al. (2022), and Canivete Cuissa & Teyssier (2022). This approach is less successful and is not as flexible as the first option, although it is indispensable in modelling deep atmospheres as it helps with hydrostatic equilibrium. The numerical precision is better when only perturbations are considered, and there is nothing really preventing gas in incompressible state from acting as compressible gas and ignoring the rotational flows. However, these methods are not mutually exclusive. In fact, combining low-dissipation solvers with well-balancing schemes, as shown by Edelmann et al. (2021), can be particularly effective for modelling highly subsonic flows in strongly stratified media (such as deep stellar convection) while keeping grid sizes reasonably small. Our initial attempts to use a simple ad hoc modification of u^⋆ and $P_{\text{tot }}^{\star }$ to imitate the transition to incompressible gas at around Mach = 0.2, as in Minoshima & Miyoshi (2021), did not result in any significant improvement; on the contrary, the test experiments produced highly undesirable results. For example, the Gresho vortex, presented in Sect. 5.12, had background gas rotating in the opposite direction. Clearly we need to do more work to balance out the equations if we want to use this method and later expand the formulation to use some form of well-balanced scheme.

Lastly, developing a fully 3D wave model for strongly interacting states, akin to the approach in Balsara (2015), is an avenue we would like to explore in the future. While our current HLLS scheme is inherently three-dimensional and unsplit, it still relies on one-dimensional Riemann problems along coordinate-aligned directions. Extending this to a truly multi-dimensional Riemann solver, incorporating a strongly interacting state that evolves self-similarly and enforces consistency with the full three-dimensional conservation laws would enable a more isotropic propagation of flow features and improved treatment of entropy generation from irreversible kinetic and magnetic energy dissipation. This would be particularly beneficial in MHD turbulence and low-Mach-number flows, where current schemes may suffer from excessive numerical dissipation or directional bias. However, significant work is required to implement such a model, including designing a suitable three-dimensional wave structure and computing numerical fluxes that maintain consistency across multiple interacting Riemann problems. We will continue our work on very low Mach-number regimes, as this is rather crucial to the applications we intend to use the solver for.

6.3 Applications

The solver was developed with the primary intent to use it in the context of solar, stellar, and planetary atmospheres. From the outermost regions, towards the cores, density, temperature, and pressure vary by orders of magnitude, but entropy per unit mass stays almost constant. Additionally, HLLS is the preferred solver in situations where magnetic energy and kinetic energy strongly dominate over thermal energy (e.g. the solar chromosphere and corona).

In strongly stratified media, evolving volumetric entropy rather than volumetric energy as a primary variable ensures that entropy gradients, which drive convection, are captured with more accuracy. This is particularly beneficial in nearly isentropic convection zones, where traditional total-energy-based schemes struggle due to the small thermal energy contrast relative to the total energy. However, even in sub-adiabatic regions, where entropy gradients are stable and convection is suppressed, the proposed method remains reliable. Since entropy is a logarithmic quantity by nature, its variability in sub-adiabatic regions does not span orders of magnitude, making it a well-conditioned variable for numerical evolution. This helps maintain accuracy in regions where the entropy contrast is small. For internal gravity waves, entropy-based formulations should naturally capture buoyancy oscillations, as entropy perturbations directly relate to restoring forces in a stratified medium. However, the very low Mach numbers typical of sub-adiabatic regions introduce additional numerical challenges, as discussed above.

We are already employing this solver to simulate the whole solar convective region from 0.655 to 0.998 R_⊙ (Popovas et al. in prep.), and we can see that the HLLS solver can maintain the hydrostatic equilibrium much better than an HLLD solver without introducing a well-balanced scheme, as in Edelmann et al. (2021) and Canivete Cuissa & Teyssier (2022), for example.

6.4 Future work

We will keep developing the solver so that it can eventually cover a larger range of Mach numbers and generally perform better. We will consider the following:

• solver modifications for very low Mach-number regimes;

• implementing a well-balanced scheme consistent with the solver;

• solver modifications for relativistic regimes;

• heavy optimisation to improve memory alignment and vectorisation that would be also suitable for modern GPUs;

• non-ideal MHD effects introducing entropy generation in Ohmic dissipation and ambipolar diffusion.

Acknowledgements

This research was supported by the Research Council of Norway through its Centres of Excellence scheme, project number 262622, and through grants of computing time from the Programme for Supercomputing, as well as through the Synergy Grant number 810218 (ERC-2018-SyG) of the European Research Council. Lastly, we are grateful to the anonymous referee for in-depth review, detailed and insightful comments, and useful suggestions, that significantly improved the manuscript.

Appendix A The numerical method

In Godunov-type schemes the volume-averaged conservative physical quantities at a next time-step are given by integrating an approximate solution to the Riemann problem with left and right states, U_L and U_R at the cell interface. An HLL Rieman solver (Harten et al. 1983) is constructed by assuming an average intermediate state between the fastest and slowest waves. The information is lost, as slow magneto-acoustic waves are merged together with Alfvén and entropy waves. An HLLC Toro et al. (1994) solver expands and estimates the middle wave of the contact surface. Lastly, HLLD Miyoshi & Kusano (2005) expands into four intermediate states. We start with computing the HLL wave speed. Usually the wave speeds are notated by S_L and S_R, but to avoid the confusion with entropy, we will note them with Σ_L and Σ_R: $$\begin{array}{rl}{\mathrm{\Sigma }}_L& =min(u_l,u_r)-max(c_{f,l},c_{f,r})\\ {\mathrm{\Sigma }}_R& =max(u_l,u_r)+max(c_{f,l},c_{f,r})\end{array}$$(A.1)

For convenience, we compute (u · B), which to avoid confusion we denote ϕ, kinetic energy E_kin and magnetic energy E_mag for left and right states as well: $${\varphi }_l=u_{l}A+v_l{B}_l+w_l{C}_l,$$(A.2) $$E_{\mathrm{k}\mathrm{i}\mathrm{n},l}=\frac{{\rho }_l}{2}(u_l^2+v_l^2+w_l^2),$$(A.3) $$E_{\mathrm{m}\mathrm{a}\mathrm{g},l}=\frac{1}{2}(A^2+B_l^2+C_l^2),$$(A.4)

and $${\varphi }_r=u_{r}A+v_r{B}_r+w_r{C}_r,$$(A.5) $$E_{kin,r}=\frac{{\rho }_r}{2}(u_r^2+v_r^2+w_r^2),$$(A.6) $$E_{mag,r}=\frac{1}{2}(A^2+B_r^2+C_r^2).$$(A.7)

We define the Lagrangian sound speed, $$\begin{array}{rl}{v}_L& =u_l-{\mathrm{\Sigma }}_L\\ v_R& ={\mathrm{\Sigma }}_R-u_r\end{array}$$(A.8)

and then compute the acoustic star state: $$u^{\star }=\frac{{\rho }_r{u}_r{v}_R+{\rho }_l{u}_l{v}_L+(P_{\mathrm{t}\mathrm{o}\mathrm{t},l}-P_{\mathrm{t}\mathrm{o}\mathrm{t},r})}{{\rho }_l{v}_L+{\rho }_r{v}_R}$$(A.9) $$P_{\mathrm{t}\mathrm{o}\mathrm{t}}^{\star }=\frac{{\rho }_r{v}_R{P}_{\mathrm{t}\mathrm{o}\mathrm{t},l}+{\rho }_l{v}_l{P}_{\mathrm{t}\mathrm{o}\mathrm{t},r}+{\rho }_r{\rho }_l{v}_R{v}_L(u_l-u_r)}{{\rho }_l{v}_L+{\rho }_r{v}_R}.$$(A.10)

The left star region variables are $${\rho }_L^{\star }=-\frac{{\rho }_l{v}_L}{{\mathrm{\Sigma }}_L-u^{\star }},$$(A.11) $$v_l^{\star }=\frac{v_l-A{B}_l(u^{\star }-u_l)}{e_l^{\star }},$$(A.12) $$w_l^{\star }=\frac{w_l-A{C}_l(u^{\star }-u_l)}{e_l^{\star }},$$(A.13) $$B_l^{\star }=\frac{B_l{e}_l}{e_l^{\star }},$$(A.14) $$C_l^{\star }=\frac{C_l{e}_l}{e_l^{\star }},$$(A.15)

where v_l and w_l are left tangential velocities, A = A_l = A_r is the normal component of the magnetic field, B_l and C_l are the left states of the tangential components of the magnetic field, $$e_l^{\star }=-\frac{{\rho }_l{v}_L}{{\mathrm{\Sigma }}_L-u^{\star }}-A^2,$$(A.16)

and $$e_l=-\frac{{\rho }_l{v}_L}{{\mathrm{\Sigma }}_L-u_l}-A^2.$$(A.17)

Entropy is updated as a passive scalar: $$S_l^{\star }=S_l\frac{{\rho }_l}{{\rho }_l^{\star }}\frac{-v_L}{{\mathrm{\Sigma }}_L-u^{\star }}.$$(A.18)

Left state Alfvén wave speed is $${\mathrm{\Sigma }}_{a,l}=u^{\star }+\frac{|A|}{\sqrt{{\rho }_l^{\star }}}.$$(A.19)

Left star region (u · B), E_kin and E_mag are $${\varphi }_l^{\star }=u^{\star }A+v_l^{\star }{B}_l^{\star }+w_l^{\star }{C}_l^{\star },$$(A.20) $$E_{\mathrm{k}\mathrm{i}\mathrm{n},l}^{\star }=\frac{{\rho }_l^{\star }}{2}(u^{\star 2}+v_l^{\star 2}+w_l^{\star 2}),$$(A.21) $$E_{\mathrm{m}\mathrm{a}\mathrm{g},l}^{\star }=\frac{1}{2}(A^2+B_l^{\star 2}+C_l^{\star 2}).$$(A.22)

Correspondingly, the right star region variables are: $${\rho }_r^{\star }=-\frac{{\rho }_r{v}_r}{{\mathrm{\Sigma }}_R-u^{\star }},$$(A.23) $$v_r^{\star }=\frac{v_r-A{B}_r(u^{\star }-u_r)}{e_r^{\star }},$$(A.24) $$w_r^{\star }=\frac{w_r-A{C}_r(u^{\star }-u_r)}{e_r^{\star }},$$(A.25) $$B_r^{\star }=\frac{B_r{e}_r}{e_r^{\star }},$$(A.26) $$C_r^{\star }=\frac{C_r{e}_r}{e_r^{\star }},$$(A.27)

where v_r and w_r are right tangential velocities, A is the normal component of the magnetic field, B_r and C_r are the tangential components of the magnetic field, $$e_r^{\star }=-\frac{{\rho }_r{v}_R}{{\mathrm{\Sigma }}_R-u^{\star }}-A^2,$$(A.28)

and $$e_r=-\frac{{\rho }_r{v}_R}{{\mathrm{\Sigma }}_R-u_l}-A^2.$$(A.29) $$S_r^{\star }=S_r\frac{{\rho }_r}{{\rho }_r^{\star }}\frac{-v_r}{{\mathrm{\Sigma }}_R-u^{\star }}.$$(A.30)

Right state Alfvén wave speed is $${\mathrm{\Sigma }}_{a,r}=u^{\star }+\frac{|A|}{\sqrt{{\rho }_r^{\star }}}.$$(A.31)

Right star region (u · B), E_kin and E_mag are $${\varphi }_r^{\star }=u^{\star }A+v_r^{\star }{B}_r^{\star }+w_r^{\star }{C}_r^{\star },$$(A.32) $$E_{\mathrm{k}\mathrm{i}\mathrm{n},r}^{\star }=\frac{{\rho }_r^{\star }}{2}(u^{\star 2}+v_r^{\star 2}+w_r^{\star 2}),$$(A.33) $$E_{\mathrm{m}\mathrm{a}\mathrm{g},r}^{\star }=\frac{1}{2}(A^2+B_r^{\star 2}+C_r^{\star 2}).$$(A.34)

Lastly, the double star region variables are: $$v^{\star \star }=\frac{v_l^{\star }\sqrt{{\rho }_l^{\star }}+v_r^{\star }\sqrt{{\rho }_r^{\star }}+\mathrm{\Xi }(B_r^{\star }-B_l^{\star })}{\sqrt{{\rho }_l^{\star }}+\sqrt{{\rho }_r^{\star }}},$$(A.35) $$w^{\star \star }=\frac{w_l^{\star }\sqrt{{\rho }_l^{\star }}+w_r^{\star }\sqrt{{\rho }_r^{\star }}+\mathrm{\Xi }(C_r^{\star }-C_l^{\star })}{\sqrt{{\rho }_l^{\star }}+\sqrt{{\rho }_r^{\star }}},$$(A.36) $$B^{\star \star }=\frac{B_r^{\star }\sqrt{{\rho }_l^{\star }}+B_l^{\star }\sqrt{{\rho }_r^{\star }}+\mathrm{\Xi }\sqrt{{\rho }_l^{\star }}\sqrt{{\rho }_r^{\star }}(v_r^{\star }-v_l^{\star })}{\sqrt{{\rho }_l^{\star }}+\sqrt{{\rho }_r^{\star }}},$$(A.37) $$C^{\star \star }=\frac{C_r^{\star }\sqrt{{\rho }_l^{\star }}+C_l^{\star }\sqrt{{\rho }_r^{\star }}+\mathrm{\Xi }\sqrt{{\rho }_l^{\star }}\sqrt{{\rho }_r^{\star }}(w_r^{\star }-w_l^{\star })}{\sqrt{{\rho }_l^{\star }}+\sqrt{{\rho }_r^{\star }}},$$(A.38)

where Ξ is a sign function, Ξ = sign(A). Note, that there is no double star variable for entropy. Then double star region (u · B) and E_mag are $${\varphi }^{\star \star }=u^{\star }A+v^{\star \star }{B}^{\star \star }+w^{\star \star }{C}^{\star \star },$$(A.39) $$E_{\text{mag }}^{\star \star }=\frac{1}{2}(A^2+B^{\star \star 2}+C^{\star \star 2}).$$(A.40)

while left and right double star regions for E_kin are correspondingly $$E_{\mathrm{k}\mathrm{i}\mathrm{n},l}^{\star \star }=\frac{{\rho }_l^{\star }}{2}(u^{\star 2}+v^{\star \star 2}+w^{\star \star 2})$$(A.41)

and $$E_{\mathrm{k}\mathrm{i}\mathrm{n},r}^{\star \star }=\frac{{\rho }_r^{\star }}{2}(u^{\star 2}+v^{\star \star 2}+w^{\star \star 2}).$$(A.42)

A.1 Godunov fluxes

The fluxes are then given by $${\mathit{F}}_{\mathrm{H}\mathrm{L}\mathrm{L}\mathrm{S}}=\{\begin{array}{ll}{\mathit{F}}_L,& \text{ if }{\mathrm{\Sigma }}_{\mathrm{L}}>0\\ {\mathit{F}}_L^{\star },& \text{ if }{\mathrm{\Sigma }}_{\mathrm{L}}\le 0\le {\mathrm{\Sigma }}_{\mathrm{a},1}\\ {\mathit{F}}_L^{\star \star },& \text{ if }{\mathrm{\Sigma }}_{\mathrm{a},1}\le 0\text{ and }{\mathrm{u}}^{\star }>0\\ {\mathit{F}}_R^{\star \star },& \text{ if }{\mathrm{u}}^{\star }\le 0\text{ and }{\mathrm{\Sigma }}_{\mathrm{a},\mathrm{r}}>0\\ {\mathit{F}}_R^{\star },& \text{ if }{\mathrm{\Sigma }}_{\mathrm{a},\mathrm{r}}\le 0\le {\mathrm{\Sigma }}_{\mathrm{R}}\\ {\mathit{F}}_R,& \text{ if }{\mathrm{\Sigma }}_{\mathrm{R}}<0\end{array},$$(A.43)

where $${\mathit{F}}_L=\left(\begin{array}{l}{\rho }_l{u}_l\\ {\rho }_l{u}_l^2+P_{\text{tot, }l}-A^2\\ {\rho }_l{u}_l{v}_l-A{B}_l\\ {\rho }_l{u}_l{w}_l-A{C}_l\\ {\rho }_l{u}_l{S}_l\\ 0\\ u_l{B}_l-A{v}_l\\ u_l{C}_l-A{w}_l\\ u_l(E_{\text{kin },l}+P_{\text{tot, }l})\\ u_l{E}_{\text{mag },l}-A{\varphi }_l\\ u_l\end{array}\right)$$(A.44) $${\mathit{F}}_L^{\star }=\left(\begin{array}{l}{\rho }_l^{\star }{u}^{\star }\\ {\rho }_l^{\star }{u}^{\star 2}+P_{tot}^{\star }-A^2\\ {\rho }_l^{\star }{u}^{\star }{v}_l^{\star }-A{B}_l^{\star }\\ {\rho }_l^{\star }{u}^{\star }{w}_l^{\star }-A{C}_l^{\star }\\ {\rho }_l^{\star }{u}^{\star }{S}_l^{\star }\\ 0\\ u^{\star }{B}_l^{\star }-A{v}_l^{\star }\\ u^{\star }{C}_l^{\star }-A{w}_l^{\star }\\ u^{\star }(E_{kin,l}^{\star }+P_{tot}^{\star })\\ u^{\star }{E}_{mag,l}^{\star }-A{\varphi }_l^{\star }\\ u^{\star }\end{array}\right)$$(A.45) $$F_L^{\star \star }=\left(\begin{array}{l}{\rho }_l^{\star }{u}^{\star }\\ {\rho }_l^{\star }{u}^{\star 2}+P_{tot}^{\star }-A^2\\ {\rho }_l^{\star }{u}^{\star }{u}^{\star \star }-A{B}^{\star \star }\\ {\rho }_l^{\star }{u}^{\star }{w}^{\star \star }-A{C}^{\star \star }\\ {\rho }_l^{\star }{u}^{\star }{S}_l^{\star }\\ 0\\ u^{\star }{B}^{\star \star }-A{v}^{\star \star }\\ u^{\star }{C}^{\star \star }-A{w}^{\star \star }\\ u^{\star }(E_{kin,l}^{\star \star }+P_{tot}^{\star })\\ u^{\star }{E}_{mag}^{\star \star }-A{\varphi }^{\star \star }\\ u^{\star }\end{array}\right)$$(A.46) $$F_R^{\star \star }=\left(\begin{array}{l}{\rho }_r^{\star }{u}^{\star }\\ {\rho }_r^{\star }{u}^{\star 2}+P_{tot}^{\star }-A^2\\ {\rho }_r^{\star }{u}^{\star }{u}^{\star \star }-A{B}^{\star \star }\\ {\rho }_r^{\star }{u}^{\star }{w}^{\star \star }-A{C}^{\star \star }\\ {\rho }_r^{\star }{u}^{\star }{S}_r^{\star }\\ 0\\ u^{\star }{B}^{\star \star }-A{v}^{\star \star }\\ u^{\star }{C}^{\star \star }-A{w}^{\star \star }\\ u^{\star }(E_{kin,r}^{\star \star }+P_{tot}^{\star })\\ u^{\star }{E}_{mag}^{\star \star }-A{\varphi }^{\star \star }\\ u^{\star }\end{array}\right)$$(A.47) $${\mathit{F}}_R^{\star }=\left(\begin{array}{l}{\rho }_r^{\star }{u}^{\star }\\ {\rho }_r^{\star }{u}^{\star 2}+P_{tot}^{\star }-A^2\\ {\rho }_r^{\star }{u}^{\star }{v}_r^{\star }-A{B}_r^{\star }\\ {\rho }_r^{\star }{u}^{\star }{w}_r^{\star }-A{C}_r^{\star }\\ {\rho }_r^{\star }{u}^{\star }{S}_r^{\star }\\ 0\\ u^{\star }{B}_r^{\star }-A{v}_r^{\star }\\ u^{\star }{C}_r^{\star }-A{w}_r^{\star }\\ u^{\star }(E_{kin,r}^{\star }+P_{tot}^{\star })\\ u^{\star }{E}_{mag,r}^{\star }-A{\varphi }_r^{\star }\\ u^{\star }\end{array}\right)$$(A.48) $${\mathit{F}}_R=\left(\begin{array}{l}{\rho }_r{u}_r\\ {\rho }_r{u}_r^2+P_{tot,r}-A^2\\ {\rho }_r{u}_r{v}_r-A{B}_r\\ {\rho }_r{u}_r{w}_r-A{C}_r\\ {\rho }_r{u}_r{S}_r\\ 0\\ u_r{B}_r-A{v}_r\\ u_r{C}_r-A{w}_r\\ u_r(E_{kin,r}+P_{tot,r})\\ u_r{E}_{mag,r}-A{\varphi }_r\\ u_r\end{array}\right)$$(A.49)

In the brackets above there are these Godunov fluxes:

1. density flux

2. normal momentum flux

3. first tangential velocity flux

4. second tangential velocity flux

5. entropy flux

6. normal magnetic field component flux

7. first tangential magnetic field component flux

8. second tangential magnetic field component flux

9. kinetic energy flux

10. magnetic energy flux

11. normal velocity

Magnetic field component fluxes are used to update the magnetic fields in order to get the $\frac{d{E}_{\text{mag }}}{dt}$. If CT method is used to update the magnetic fields, this step is still necessary for consistency (temporarily storing $\frac{d{E}_{\text{mag }}}{dt}$ separately), as the actual magnetic field update is done through the electromotive force from the cell edges. The kinetic and magnetic energy fluxes are used to compute ∇ · F_kin and ∇ · F_mag respectively, while normal velocity (last item) is used to compute P (∇ · u) to get the gas pressure work W_gas.

References

All Figures

	Fig. 1 Convergence in norm of L1 error vector for contact, fast, slow, and Alfvén waves after propagating a distance of one wavelength. The convergence rate P is given in the legend.
In the text

	Fig. 2 Entropy wave test. Density profiles at starting time (black) and for different wave speeds ux: blue – 1.0, red – 0.1, green dashed – 0.01, magenta dotted – 0.001. The overlap of wave profiles is perfect, with expected slopes present.
In the text

	Fig. 3 Shu & Osher shock-tube test. Density profile at time=1.8; black: reference HLLD solver with 1500 cells; red: HLLS solver with 1500 cells; blue dots: HLLS solver with 500 cells.
In the text

	Fig. 4 Brio & Wu shock-tube test. Density profile (left) and y component of magnetic field (right) at time=1.0; black: reference HLLD solver with 1200 cells; red: HLLS solver with 1200 cells; blue dots: HLLS solver with 400 cells. Black and red overlap each other nearly perfectly.
In the text

	Fig. 5 HD Kelvin-Helmholtz instability. Density profile (left) and time evolution of kinetic energy (right).
In the text

	Fig. 6 Rayleigh–Taylor instability in low (128 × 256) resolution. We plot density at different times.
In the text

	Fig. 7 Same as Fig. 6, but in 768 × 1536 resolution.
In the text

	Fig. 8 ‘Hot-bubble’ experiment at time=150 and 300 s. Both panels show fluctuations of ΔA/A, with the min/max values given in the floating box.
In the text

	Fig. 9 MHD blast experiment. Density (left) and magnetic pressure (right) at time=0.021. We note the sharp shock-fronts.
In the text

	Fig. 10 Orszag-Tang vortex test at time = 0.5. Density (left), entropy per unit mass (middle), and magnetic pressure (right).
In the text

	Fig. 11 Orszag-Tang vortex test at time = 0.75. Density (left), entropy per unit mass (middle), and magnetic pressure (right). We note the plasmoids forming.
In the text

	Fig. 12 Current-sheet test at different times. Density (top) and magnetic pressure (bottom).
In the text

	Fig. 13 Current-sheet test at different times. Density (top) and magnetic pressure (bottom).
In the text

	Fig. 14 Magnetic-loop advection experiment. The initial magnetic pressure (left) and after the loop has been advected around the grid twice.
In the text

	Fig. 15 Gresho vortex test. The initial setup (left) and after two rotation periods with different Mach numbers: Ma = 0.1, 0.01, 0.001. The rightmost panel is for a higher resolution run (see text for more details).
In the text

	Fig. 16 Gas density (left) and pressure (right) in the Gresho vortex test with Ma = 0.01. We note that the quantities are within the noise limit.
In the text

	Fig. 17 Magnetic rotor test. From left to right: Density, magnetic pressure, Mach number, and $\frac{\mathrm{\nabla }\cdot \mathbf{B}}{|\mathbf{B}|}$ (see text for more details).
In the text