© The Authors 2024

1 Introduction

High-resolution schemes for gas dynamics (see, e.g., van Leer 1979; Colella & Woodward 1984; Harten et al. 1987; Colella 1990; Liu et al. 1994; Jiang & Shu 1996; Colella & Sekora 2008; Toro 2009; Balsara 2017) are routinely used for modeling a broad variety of astrophysical flow phenomena. Their popularity derives from their conservation properties and robustness, which allow them to accurately capture both smooth and discontinuous solutions on the same computational grid without sacrificing numerical stability.

These schemes are based on higher-order extensions of the original first-order accurate method of Godunov (1959) and their time-integration algorithm is typically carried out in three steps. First, a pair of Riemann states is reconstructed at each grid cell interface by applying high-order monotonic interpolants to a set of cell-averaged hydrodynamic quantities. Second, the resulting Riemann problems are solved (either exactly or approximately) to obtain fluxes across every cell boundary. Finally, the cell surface integral of the fluxes is evaluated, allowing the cell-volume-averaged state quantities to be advanced in time¹.

In most high-order Godunov schemes, the solution strategy of the Riemann problem is independent of the spatial interpolant used for reconstructing the Riemann states. Therefore, many different schemes can be built from the same numerical framework. The choices made in the construction of a particular scheme, however, do have a strong effect on its accuracy, that is the difference between the numerical (U) and the true (u) solution, $$\parallel U-u\parallel =\mathcal{O}\left(\mathrm{\Delta }{x}^m\right)+\mathcal{O}\left(\mathrm{\Delta }{t}^n\right),$$(1)

computed in some norm ∥ ⋅ ∥ (see, e.g., LeVeque 2002). Here, Δx is the width of the grid cell and Δt is the time step. Although the formal order of the spatial and temporal accuracy of a Godunov-like scheme, “m” and “n” in Eq. (1), can be derived for smooth flows, they do not give any information about the magnitude of the numerical errors generated on a given grid, which is problem-dependent. The convergence rates can also be significantly lower than the formal order of accuracy of the scheme for problems that admit nonsmooth solutions. Consequently, formally higher-than-second-order interpolants do not always outperform simpler linear spatial reconstruction schemes when large gradients or discontinuities are present in the flow (Greenough & Rider 2003). Moreover, if the flow is stochastic or chaotic, as in the case of turbulence, convergence may not be achieved in the sense of Eq. (1), but rather the quality of the numerical results can only be judged in terms of the global or ensemble-averaged quantities that characterize the flow and its evolution.

Given these considerations, it is impossible to generalize the convergence properties of a certain combination of numerical methods in a Godunov-type scheme, so they have to be explored by running numerical tests. For such tests to be significant, they have to be challenging enough and close to the actual application case. The performance of the numerical scheme is another crucial aspect to be considered alongside its accuracy, especially in astrophysical simulations, which often test out the limits of what is feasible nowadays with available computational resources. Therefore, the question arises of what combination of different ingredients in a Godunov-type scheme should be used to produce the desired solution at a minimal computational cost.

Several comparison studies have been presented in the literature with the aim of shedding light on the behavior of different high-resolution schemes in simulations of complex hydrodynamic phenomena, such as forced turbulence (Klingenberg et al. 2007; Kritsuk et al. 2011; San & Kara 2015; Radice et al. 2015; Seo & Ryu 2023), convection (Müller 2020), jet evolution (Beckwith & Stone 2011; Musoke et al. 2020), magneto-rotational instabilities in accretion disks (Flock et al. 2010), Richtmyer-Meshkov instabilities (Latini et al. 2007), and shear instabilities (McNally et al. 2012; Lecoanet et al. 2017). These numerical experiments focused on supersonic or mildly subsonic flow regimes, for which Godunov-type methods are highly optimized (see, e.g., LeVeque 2002; Toro 2009). Nonetheless, highresolution schemes have also been proven to be a powerful tool for modeling regimes of low Mach numbers ($\mathcal{M}$: = ∣V∣/c ≲ 0.1, where V is the fluid velocity and c is the sound speed), especially in simulations of terrestrial (see, e.g., Day & Bell 2000; Klein 2009; Dumbser et al. 2009; Motheau et al. 2018) and stellar (see, e.g., Meakin & Arnett 2007; Muthsam et al. 2010; Woodward et al. 2014; Goffrey et al. 2017; Müller 2020; Horst et al. 2021; Canivete Cuissa & Teyssier 2022) flows. To our knowledge, other than idealized tests, no extensive work along the lines of the aforementioned comparison studies has been presented for low-Mach-number flows yet. Only a few studies evaluated the impact of the order of the spatial and temporal discretization in the numerical scheme on the properties of highly subsonic turbulent flows, but they kept the Riemann solver fixed (Wongwathanarat et al. 2016; Teissier & Müller 2023).

In this work, we used our fully compressible Seven-League Hydro (SLH) code to test 18 different combinations of spatial reconstruction schemes and Riemann solvers on two test problems in which flows are highly subsonic (10⁻³ ≲ $\mathcal{M}$ ≲ 10⁻¹). In particular, we considered a two-dimensional (2D) Kelvin–Helmholtz instability with smooth initial conditions and a 3D, turbulent convection zone that entrains material from an upper, stably stratified layer, where internal waves are free to propagate. The initial conditions of the latter test were adopted from the work of Andrassy et al. (2022). Here, we opted to reduce the strength of the heat source driving the convection in order to achieve lower convective speeds than those obtained by Andrassy et al. (2022). Also, contrary to that work, we provided performance measurements for all the methods tested in our study. In both tests, we conducted a resolution study to analyze the convergence of the numerical results obtained by each method.

The paper is structured as follows: in Sec. 2, we provide a detailed description of the equations solved and the numerical methods included in this study. In Sec. 3, we measure the convergence properties of different Godunov-type schemes for the 2D, Kelvin–Helmholtz instability test problem (see Sec. 3.1) and for the 3D setup involving turbulent convection, convective boundary mixing, and wave excitation (see Sec. 3.2). In Sec. 4, we use the kinetic energy spectrum of the convective flows simulated in the latter test, which is close to a real astrophysical application, to provide measurements of the computational cost per fixed accuracy for each method. Finally, in Sect. 5, we summarize the main results and we give some guidance on which methods to use for specific applications.

2 Methods

2.1 Governing equations

We solved the fully compressible, inviscid Euler equations with a source term S in the integral form $$\frac{1}{|\mathrm{\Omega }|}\frac{\mathrm{\partial }}{\mathrm{\partial }t}\left({\int }_{\mathrm{\Omega }}\mathit{U}\mathrm{d}\mathrm{\Omega }\right)+\frac{1}{|\mathrm{\Omega }|}{\oint }_{\mathrm{\partial }\mathrm{\Omega }}\mathbb{T}\cdot \mathit{n}\mathrm{d}A=\frac{1}{|\mathrm{\Omega }|}{\int }_{\mathrm{\Omega }}\mathit{S}\mathrm{d}\mathrm{\Omega },$$(2)

where $$\mathit{U}=\left[\begin{array}{c}\rho \\ \rho u\\ \rho v\\ \rho w\\ \rho e_{\text{tot }}\\ \rho X\end{array}\right]$$(3)

is the set of conserved quantities, ∣Ω∣ is the volume of a fluid element Ω enclosed by a surface ∂Ω, n is the outward normal vector to the surface, and $\mathbb{T}$ = [F∣G∣H] is a tensor defined by the flux vectors $$\mathit{F}=\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uv\\ \rho uw\\ (\rho e_{\mathrm{t}\mathrm{o}\mathrm{t}}+p)u\\ \rho Xu\end{array}\right],\mathit{G}=\left[\begin{array}{c}\rho v\\ \rho uv\\ \rho v^2+p\\ \rho vw\\ (\rho e_{\mathrm{t}\mathrm{o}\mathrm{t}}+p)v\\ \rho Xv\end{array}\right],\mathit{H}=\left[\begin{array}{c}\rho w\\ \rho uw\\ \rho vw\\ \rho w^2+p\\ (\rho e_{\mathrm{t}\mathrm{o}\mathrm{t}}+p)w\\ \rho Xw\end{array}\right].$$(4)

Here, ρ denotes the mass density, V = (u, v, w) the velocity field, $e_{\mathrm{t}\mathrm{o}\mathrm{t}}=e_{\mathrm{i}\mathrm{n}\mathrm{t}}+\frac{1}{2}|\mathit{V}{|}^2$ the total energy per unit mass, e_int the specific internal energy, and X the mass fraction of a passive scalar used as a tracer advected with the fluid. The system represented by Eq. (2) is closed by an equation of state (EoS), which gives the gas pressure as a function of the density and internal energy, $$p=p(\rho ,e_{\text{int }}).$$(5)

In this work, we only considered a perfect gas with a given adiabatic index γ, for which $$p(\rho ,e_{\text{int }})=(\gamma -1)\rho e_{\text{int }}.$$(6)

For modeling the test problem introduced in Sec. 3.2, which involves the presence of a gravitational field, g = (g_x, g_y, g_z), the corresponding source term, $${\mathit{S}}_{\text{grav }}=\left[\begin{array}{c}0\\ \rho g_x\\ \rho g_y\\ \rho g_z\\ \rho \mathit{g}\cdot \mathit{V}\\ 0\end{array}\right],$$(7)

must be included in the right-hand-side term of Eq. (2). Because we assumed g to be time-independent, we opted to solve an equivalent form of Eq. (2), in which the gravitational potential ϕ was directly added to e_tot, eliminating the ρg ⋅ V source term from the energy equation, $$e_{\text{tot }}\mapsto e_{\text{tot }}+\varphi ,{\mathit{S}}_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{v}}=\left[\begin{array}{c}0\\ \rho g_x\\ \rho g_y\\ \rho g_z\\ \rho \mathit{g}\cdot \mathit{V}\\ 0\end{array}\right]\mapsto \left[\begin{array}{c}0\\ \rho g_x\\ \rho g_y\\ \rho g_z\\ 0\\ 0\end{array}\right].$$(8)

Numerical schemes that solve this form of Eq. (2) are capable of conserving the total energy of the system over time exactly. Such a property is particularly important in low-Mach-number hydrodynamics, where even small energy conservation errors can become comparable to the kinetic energy content of the flows (Müller 2020; Edelmann et al. 2021).

2.2 Seven-League Hydro code

In our study, Eq. (2) was solved numerically using the Seven-League Hydro code (SLH, Miczek 2013; Edelmann 2014), which was originally developed to model the broad variety of hydrodynamic processes that characterize the deep interiors of stars, such as shear instabilities (Edelmann et al. 2017), excitation of internal waves (Horst et al. 2020), convective boundary mixing (Horst et al. 2021; Andrassy et al. 2022, 2024), and turbulent dynamos (Leidi et al. 2022, 2023). SLH makes use of the finite-volume discretization on an arbitrarily curvilinear, but logically rectangular, Eulerian grid to retain the conservation properties of the fluid-dynamics equations. Numerical solutions to Eq. (2) are computed by means of Godunov-type methods based on the definition of Riemann problems at cell interfaces. The code is parallelized using the Message Passing Interface (MPI) and it has been proven to scale up to several hundred thousand processes (Edelmann & Röpke 2016).

SLH allows the user to choose among many different numerical options at compile time, which makes this code perfectly suited to run the comparison study introduced in Sect. 1. In particular, in addition to the well-known Rusanov (Rusanov 1962), Roe (Roe 1981), and Harten-Lax-van Leer-Contact (HLLC, Toro et al. 1994) approximate Riemann solvers, SLH adopts special low-Mach-number methods (Liou 2006; Miczek et al. 2015; Minoshima & Miyoshi 2021) to reduce the excessive numerical dissipation introduced by shock-capturing schemes at low Mach numbers (see Sec. 2.4). A wide spectrum of spatial reconstruction methods can be used to generate a pair or Riemann states at each grid cell interface, ranging from (first-order accurate) constant reconstruction to very high order methods, some of which are described in Sec. 2.3. In problems involving the presence of a gravitational field, the deviation well-balancing method (Berberich et al. 2021; Edelmann et al. 2021) is used to preserve hydrostatic solutions and to reduce the strength of spurious flows generated by grid discretization errors in strongly stratified media.

In SLH, the cell-volume-averaged source term and the cell-surface-averaged fluxes are approximated using the midpoint method, making the code at best second-order accurate in space. On a 3D, evenly spaced Cartesian grid, the final expressions for these integrals read $$\begin{array}{rl}& \frac{1}{\mathrm{\Delta }V}{\int }_{{\mathrm{\Omega }}_{i,jk}}\mathit{S}\mathrm{d}x\mathrm{d}y\mathrm{d}z={\mathcal{S}}_{i,j,k}+\mathcal{O}(\mathrm{\Delta }{x}^2),\end{array}$$(9) $$\begin{array}{rl}& \frac{1}{\mathrm{\Delta }A}{\int }_{\mathrm{\Delta }{A}_{i+1/2,jk}}\mathit{F}\mathrm{d}y\mathrm{d}z={\mathcal{F}}_{i+1/2,j,k}+\mathcal{O}(\mathrm{\Delta }{x}^2),\end{array}$$(10) $$\begin{array}{rl}& \frac{1}{\mathrm{\Delta }A}{\int }_{\mathrm{\Delta }{A}_{i,j+1/2,k}}\mathit{G}\mathrm{d}x\mathrm{d}z={\mathcal{G}}_{i,j+1/2,k}+\mathcal{O}(\mathrm{\Delta }{x}^2),\end{array}$$(11) $$\begin{array}{rl}& \frac{1}{\mathrm{\Delta }A}{\int }_{\mathrm{\Delta }{A}_{i,jk+1/2}}\mathit{H}\mathrm{d}x\mathrm{d}y={\mathcal{H}}_{i,j,k+1/2}+\mathcal{O}(\mathrm{\Delta }{x}^2),\end{array}$$(12)

where S_i,j,k is the point value of the source term S at the center of the cell represented by the set of indices (i, j, k) and vector quantities such as $\mathcal{F}$_i+1/2,j,k refer to the face-centered value of the flux at the boundary between two adjacent cells, in this case (i, j, k) and (i + 1 , j, k). The volume of the cell and the surface of a cell face are ΔV = Δx³ and ΔA = Δx², respectively. The discretized source term and fluxes in Eqs. (9)–(12) are used to build a semidiscrete version of Eq. (2) which leaves the problem continuous in time, following the method of lines (see, e.g., LeVeque 2002), $$\begin{array}{rl}\frac{\mathrm{\partial }{\overline{\mathit{U}}}_{i,j,k}}{\mathrm{\partial }t}=& -\frac{1}{\mathrm{\Delta }x}({\mathcal{F}}_{i+1/2,j,k}-{\mathcal{F}}_{i-1/2,j,k}\\ & +{\mathcal{G}}_{i,j+1/2,k}-{\mathcal{G}}_{i,j-1/2,k}\\ & +{\mathcal{H}}_{i,j,k+1/2}-{\mathcal{H}}_{i,j,k-1/2})\\ & +{\mathit{S}}_{i,j,k}.\end{array}$$(13)

The time update on the cell-volume-averaged conserved variables, ${\overline{\mathit{U}}}_{i,j,k}$, is then carried out in a dimensionally unsplit fashion using explicit or implicit time stepping. Here, we only considered a limited set of all of the numerical methods available in SLH to avoid constructing a too large parameter space. In the following sections, we provide a detailed description of the algorithms that we used for running the tests presented in Sects. 3.1 and 3.2.

2.3 Spatial reconstruction methods

We used six reconstruction methods as summarized in Table 1. We did not attempt to be exhaustive in the choice of our methods, which is why, for example, essentially nonoscillatory (ENO) and weighted ENO (WENO) schemes (see, e.g., Liu et al. 1994; Jiang & Shu 1996; Shu 2009) were left out from our comparison study². We also excluded multidimensional reconstruction methods. However, the methods we did include cover a wide range in complexity and formal order of accuracy and several of them are often used in stellar hydrodynamics.

The simplest are unlimited linear (LIN) and parabolic (PAR) methods, which have proved to be well-behaved in implicit simulations of slow flows using our SLH code (e.g., Horst et al. 2020, 2021; Andrassy et al. 2024). Their main disadvantage – the generation of artificial oscillations around steep gradients – can be eliminated using slope limiters. Out of a wide spectrum of limiters available, we have decided to include the popular van Leer limiter (van Leer 1974) in combination with linear reconstruction (LIN+VL). This limiter has the total-variation-diminishing (TVD) property and eliminates the oscillations completely³. As examples of higher-order methods with limiters, we included two versions of the widely used piecewise-parabolic method (PPM: Colella & Woodward 1984, CW84 hereinafter and Colella & Sekora 2008, CS08 hereinafter), for which we introduce the acronyms PPM84 and PPM08, respectively. We also constructed a hybrid method, which we named piecewise sextic hybrid (PSH). It combines unlimited sextic reconstruction for dynamic variables (ρ, V, p) with the PPM08 method for passive scalars. We provide lower-order alternatives to the PSH method in Appendix A, although we did not test them in this study.

In the remainder of this section, we provide detailed descriptions of all the six reconstruction methods in unified notation. We do this to (i) maximize the reproducibility of our work, (ii) simplify the methods’ more general original forms (e.g., for nonuniform grids), (iii) make clear what choices we made if several options were available, and (iv) to point out a few typographic mistakes in the original works.

We refer to the variable being reconstructed as a in Sects. 2.3.1–2.3.6. Although reconstruction can be performed in physical coordinates, we chose to do it in logical coordinates defined by the cell index i, which leads to much simpler expressions. This choice did not influence our results in any way since we used uniform Cartesian grids. In some cases, we introduce the continuous linear logical coordinate ζ such that the coordinate of the center of cell i is ζ_i = i and the left and right cell interfaces are located at ${\zeta }_{i-1/2}=i-\frac{1}{2}$ and ${\zeta }_{i+1/2}=i+\frac{1}{2}$, respectively. The reconstruction is usually discontinuous at the interfaces. We use the notation a_i+1/2,L and a_i+1/2,R to denote the reconstructed states on the left and right side of the interface at ζ_i+1/2, respectively. Reconstruction is performed along each spatial axis separately, that is to say the reconstruction procedure is always one-dimensional. To highlight the difference between cell averages and point values, we use the notation ;ā_i for the average value of a within cell i.

In this work, we always reconstructed the set of cell-volume-averaged primitive variables $\overline{q}$, where q = (ρ, V, p, X), because it helps reducing oscillations near discontinuities as compared to reconstructing cell-volume-averaged conserved quantities $\overline{\mathit{U}}$. In SLH, transformations between $\overline{q}$ and $\overline{\mathit{U}}$ are performed using 2nd-order approximations, $$\begin{array}{rl}\overline{\mathit{q}}& =\mathit{\xi }(\overline{\mathit{U}})+\mathcal{O}\left(\mathrm{\Delta }{x}^2\right),\\ \overline{\mathit{U}}& ={\mathit{\xi }}^{-1}(\overline{\mathit{q}})+\mathcal{O}\left(\mathrm{\Delta }{x}^2\right),\end{array}$$(14)

where ξ is a nonlinear, invertible transformation, $$\mathit{\xi }(\mathit{U})=\left[\begin{array}{c}{U}_1\\ U_2/U_1\\ U_3/U_1\\ U_4/U_1\\ (\gamma -1)(U_5-\frac{1}{2U_1}(U_2^2+U_3^2+U_4^2))\\ U_6/U_1\end{array}\right],$$(15)

and U_i is the i th component of U.

2.3.1 The LIN method

The unlimited linear reconstruction method is based on a linear approximation to the underlying function a(ζ). Its slope δ_i (in cell-index coordinates) is estimated using the central difference (the Fromm method) $${\delta }_i=\frac{{\overline{a}}_{i+1}-{\overline{a}}_{i-1}}{2}={\mathrm{\Delta }x{\mathrm{\partial }}_{x}a|}_{x_i}+O\left(\mathrm{\Delta }{x}^3\right),$$(16)

using which we obtain the reconstructed states $$a_{i-1/2,\mathrm{R}}={\overline{a}}_i-\frac{{\delta }_i}{2}+O\left(\mathrm{\Delta }{x}^2\right),$$(17) $$a_{i+1/2,\mathrm{L}}={\overline{a}}_i+\frac{{\delta }_i}{2}+O\left(\mathrm{\Delta }{x}^2\right).$$(18)

The LIN method is exact wherever a(ζ) is locally linear, 2nd-order accurate for general but smooth functions a(ζ), and it requires two ghost cells at domain boundaries⁴.

2.3.2 The LIN+VL method

The van-Leer-limited linear reconstruction method is particularly easy to describe in terms of the 2nd-order-accurate, interface-centered slopes $${\delta }_{i-1/2}={\overline{a}}_i-{\overline{a}}_{i-1},$$(19) $${\delta }_{i+1/2}={\overline{a}}_{i+1}-{\overline{a}}_i.$$(20)

The final slope δ_i,lim is then obtained by applying the limiter of van Leer (1974), $${\delta }_{i,lim}=\{\begin{array}{ll}\frac{2{\delta }_{i-1/2}{\delta }_{i+1/2}}{{\delta }_{i-1/2}+{\delta }_{i+1/2}}& \text{ if }{\delta }_{i-1/2}{\delta }_{i+1/2}>0,\\ 0& \text{ otherwise, }\end{array}$$(21)

which gives the reconstructed states $$a_{i-1/2,\mathrm{R}}={\overline{a}}_i-\frac{{\delta }_{i,\mathrm{l}\mathrm{i}\mathrm{m}}}{2}+O\left(\mathrm{\Delta }{x}^2\right),$$(22) $$a_{i+1/2,\mathrm{L}}={\overline{a}}_i+\frac{{\delta }_{i,\mathrm{l}\mathrm{i}\mathrm{m}}}{2}+O\left(\mathrm{\Delta }{x}^2\right).$$(23)

The limiter makes the slope less steep where a(ζ) is strongly curved and flat at local extrema. With a smooth and monotonic function a(ζ), the effect of the limiter weakens upon grid refinement and the left and right states converge to those provided by the LIN method. This makes the LIN+VL method formally 2nd-order accurate away from any extrema. Two ghost cells are required at domain boundaries.

2.3.3 The PAR method

The unlimited parabolic method assumes that a(ζ) can within cell i be described using the parabola $$a(\zeta )=\sum _{n=0}^2{c}_n{(\zeta -{\zeta }_i)}^n.$$(24)

The three coefficients c_n are uniquely determined by the requirement that the averages of a(ζ) in cells i − 1, i, and i + 1 equal ā_i-1, ā_i, and ā_i+1, respectively. The reconstructed states are then obtained by evaluating Eq. (24) at ζ_i−1/2 and ζ_i+1/2, respectively. The resulting expressions are $$a_{i-1/2,\mathrm{R}}=\frac{1}{6}(2{\overline{a}}_{i-1}+5{\overline{a}}_i-{\overline{a}}_{i+1})+O\left(\mathrm{\Delta }{x}^3\right),$$(25) $$a_{i+1/2,\mathrm{L}}=\frac{1}{6}(-{\overline{a}}_{i-1}+5{\overline{a}}_i+2{\overline{a}}_{i+1})+O\left(\mathrm{\Delta }{x}^3\right).$$(26)

The PAR method is exact wherever a(ζ) is locally parabolic, 3rd-order accurate for general but smooth functions a(ζ), and it requires two ghost cells at domain boundaries.

2.3.4 The PPM84 method

The piecewise parabolic reconstruction of CW84 is a two-step process. The first is based on the 4th-order-accurate interpolation formula $$a_{i+1/2}=\frac{1}{2}({\overline{a}}_i+{\overline{a}}_{i+1})-\frac{1}{6}({\delta }_{i+1}-{\delta }_i)={a|}_{x_{i+1/2}}+O\left(\mathrm{\Delta }{x}^4\right).$$(27)

This expression is the equivalent of Eq. (1.6) of CW84 in the special case of a uniform grid. CW84 then replace δ_i by the limited value $${\delta }_{i,lim}=\{\begin{array}{ll}\mathrm{sgn}\left({\delta }_i\right)min(\left|{\delta }_i\right|,2\left|{\delta }_{i-1/2}\right|,2\left|{\delta }_{i+1/2}\right|)& \text{ if }{\delta }_{\mathrm{i}-1/2}{\delta }_{i+1/2}>0\\ 0& \text{ otherwise }\end{array}$$(28)

which is the monotonized central limiter of van Leer (1977). In our implementation, we did not set δ_i,lim = 0 if δ_i−1/2 δ_i+1/2 ≤ 0 (i.e., at local extrema). We have tested that, thanks to the presence of another limiter in PPM84 (see below), this modification has essentially no influence on the results. However, it makes the code faster because it removes three conditional expressions per reconstruction step (we need δ_i−1,lim, δ_i,lim, and δ_i+1,lim to obtain a_i−1/2 and a_i+1/2).

The interpolated value a_i+1/2 is initially assigned to both a_i+1/2,L and a_i+1/2,R, meaning there is no discontinuity at the interface. However, CW84 approximate the distribution of variable a in cell i by the parabola uniquely defined by a_i−1/2,R, a_i-1/2,L, and the cell average a_i. This parabola may in some cases take on values outside of the interval defined by a_i−1/2,R and a_i−1/2,L, that is overshoots may appear within the cell. To prevent this, a second limiting step is introduced. We express it in terms of the differences $$a_i^{-}=a_{i-1/2,\mathrm{R}}-{\overline{a}}_i,$$(29) $$a_i^{+}=a_{i+1/2,\mathrm{L}}-{\overline{a}}_i.$$(30)

The limiter is then defined by the variable replacements $$\begin{array}{rll}{a}_i^{-}\mapsto 0,a_i^{+}\mapsto 0& \text{ if }& a_i^{-}{a}_i^{+}\ge 0\\ a_i^{-}\mapsto -2a_i^{+}& \text{ if }& |a_i^{-}|>2|a_i^{+}|,\\ a_i^{+}\mapsto -2a_i^{-}& \text{ if }& |a_i^{+}|>2|a_i^{-}|.\end{array}$$(31)

The reconstructed states a_i−1/2,R and a_i−1/2,L are recovered using Eqs. (29) and (30). Equation (31) is equivalent to Eq. (1.10) of CW84. This second limiter introduces discontinuities at cell interfaces in regions where the gradient of a(ζ) changes rapidly and it flattens the assumed parabola at local extrema.

Unlike CW84, we did not use the parabolic model of a(ζ) inside the cell in any way. We only needed the states at the two sides of each interface to construct Riemann problems and time integration was done using a Runge-Kutta scheme (see Sec. 2.5 for details). However, the second limiter (Eq. (31)), which is based on the parabolic model, was still needed to remove oscillations and to introduce dissipation where necessary.

Although the interpolation formula that PPM84 starts with is 4th-order accurate, the slope flattening introduced at all extrema reduces the practically attainable order of accuracy substantially for nonmonotonic solutions even if they are smooth. Our 1D experiment in Appendix B.1 gives the order of 2.3 whereas Colella & Sekora (2008) reach the order of 2.6 in a similar advection experiment with a Gaussian-shaped profile and PPM84 reconstruction. The PPM84 method requires three ghost cells at domain boundaries.

2.3.5 The PPM08 method

The piecewise parabolic method of CS08 is based on ideas similar to those of CW84 in constructing the PPM84 scheme and PPM08 also contains two limiters. However, the limiters are modified such that PPM08 models smooth extrema instead of flattening them.

In PPM08, the first estimate of a_i+1/2 is obtained using the 6th-order-accurate interpolation formula (cf. Eq. (17) of CS08) $$\begin{array}{rl}{a}_{i+1/2}& =\frac{37}{60}({\overline{a}}_i+{\overline{a}}_{i+1})-\frac{2}{15}({\overline{a}}_{i-1}+{\overline{a}}_{i+2})+\frac{1}{60}({\overline{a}}_{i-2}+{\overline{a}}_{i+3})\end{array}$$(32) $$={a|}_{x_{i+1/2}}+O\left(\mathrm{\Delta }{x}^6\right).$$(33)

If a_i+1/2 does not satisfy the condition (cf. Eq. (13) of CS08) $$min({\overline{a}}_i,{\overline{a}}_{i+1})\le a_{i+1/2}\le max({\overline{a}}_i,{\overline{a}}_{i+1})$$(34)

a limiter is applied. It is based on three 2nd-order-accurate second derivatives, $$\begin{array}{rl}{\left({\mathrm{D}}^{2}a\right)}_{i+1/2}& =3({\overline{a}}_i-2a_{i+1/2}+{\overline{a}}_{i+1})={{\mathrm{\partial }}_x^{2}a|}_{x_{i+1/2}}\mathrm{\Delta }{x}^2+O\left(\mathrm{\Delta }{x}^4\right),\end{array}$$(35) $$\begin{array}{rl}{\left({\mathrm{D}}^{2}a\right)}_i=& {\overline{a}}_{i-1}-2{\overline{a}}_i+{\overline{a}}_{i+1}={{\mathrm{\partial }}_x^{2}a|}_{x_i}\mathrm{\Delta }{x}^2+O\left(\mathrm{\Delta }{x}^4\right),\end{array}$$(36) $$\begin{array}{rl}{\left({\mathrm{D}}^{2}a\right)}_{i+1}& ={\overline{a}}_i-2{\overline{a}}_{i+1}+{\overline{a}}_{i+2}={{\mathrm{\partial }}_x^{2}a|}_{x_{i+1}}\mathrm{\Delta }{x}^2+O\left(\mathrm{\Delta }{x}^4\right).\end{array}$$(37)

In Eq. (35), the more common finite difference formula with a prefactor 4 would arise if the involved quantities were of the same kind, meaning all three point values, or all three averages. These derivatives are combined to obtain a limited derivative (D²a)_i+1/2,lim such that $$\begin{array}{rl}{\left({\mathrm{D}}^{2}a\right)}_{i+1/2,\mathrm{l}\mathrm{i}\mathrm{m}}=& \mathrm{sgn}\left[{\left({\mathrm{D}}^{2}a\right)}_{i+1/2}\right]min(\left|{\left({\mathrm{D}}^{2}a\right)}_{i+1/2}\right|,\\ & C\left|{\left({\mathrm{D}}^{2}a\right)}_i\right|,C\left|{\left({\mathrm{D}}^{2}a\right)}_{i+1}\right|)\end{array}$$(38)

if all three derivatives have the same sign and $${\left({\mathrm{D}}^{2}a\right)}_{i+1/2,\text{ lim }}=0$$(39)

otherwise. We used C = 1.25 in Eq. (38). The limited derivative is then used to modify the value of a_i+1/2, $$a_{i+1/2}\mapsto \frac{1}{2}({\overline{a}}_i+{\overline{a}}_{i+1})-\frac{1}{6}{\left({\mathrm{D}}^{2}a\right)}_{i+1/2,\mathrm{l}\mathrm{i}\mathrm{m}}.$$(40)

Equation (19) of CS08 is equivalent to our Eq. (40) except for the factor in front of the second term, which is $\frac{1}{3}$ in CS08. This is likely a typographical error because the factor of $\frac{1}{6}$ is needed to obtain the original, unlimited value of a_i+1/2 when the solution is so smooth that the limiter does not change the second derivative significantly.

Just as in the PPM84 method, the interpolated value a_i+1/2 is initially assigned to both a_i+1/2,L and a_i+1/2,R, meaning there is no discontinuity at the interface. The second limiting step depends on whether cell i is in the vicinity of a local extremum or not. If (c.f. Eq. (20) of CS08) $$\begin{array}{c}({\overline{a}}_i-a_{i-1/2,\mathrm{R}})(a_{i+1/2,\mathrm{L}}-{\overline{a}}_i)\le 0\text{ or }\\ ({\overline{a}}_i-{\overline{a}}_{i-1})({\overline{a}}_{i+1}-{\overline{a}}_i)\le 0,\end{array}$$(41)

cell i is close to a local extremum, which should be preserved if smooth enough. The second derivatives $${\left({\mathrm{D}}^{2}a\right)}_i^{\ast }=6(a_{i-1/2,\mathrm{R}}-2{\overline{a}}_i+a_{i+1/2,\mathrm{L}}),$$(42) $${\left({\mathrm{D}}^{2}a\right)}_{i-1}={\overline{a}}_{i-2}-2{\overline{a}}_{i-1}+{\overline{a}}_i,$$(43) $${\left({\mathrm{D}}^{2}a\right)}_i={\overline{a}}_{i-1}-2{\overline{a}}_i+{\overline{a}}_{i+1},$$(44) $${\left({\mathrm{D}}^{2}a\right)}_{i+1}={\overline{a}}_i-2{\overline{a}}_{i+1}+{\overline{a}}_{i+2}.$$(45)

are combined to judge the solution’s smoothness (CS08 have one wrong index in their equivalent of our Eq. (43), cf. their Eq. (21)). We then set $$\begin{array}{rl}{\left({\mathrm{D}}^{2}a\right)}_{i,\mathrm{l}\mathrm{i}\mathrm{m}}^{\ast }=& \mathrm{sgn}\left[{\left({\mathrm{D}}^{2}a\right)}_i^{\ast }\right]min\left(\right|{\left({\mathrm{D}}^{2}a\right)}_i^{\ast }|,\\ & C\left|{\left({\mathrm{D}}^{2}a\right)}_{i-1}\right|,C\left|{\left({\mathrm{D}}^{2}a\right)}_i\right|,C\left|{\left({\mathrm{D}}^{2}a\right)}_{i+1}\right|)\end{array}$$(46)

if all four second derivatives have the same sign and $${\left({\mathrm{D}}^{2}a\right)}_{i,\text{ lim }}^{\ast }=0$$(47)

otherwise. Finally, the reconstructed states are updated, $$a_{i-1/2,\mathrm{R}}\mapsto {\overline{a}}_i+(a_{i-1/2,\mathrm{R}}-{\overline{a}}_i)\frac{{\left({\mathrm{D}}^{2}a\right)}_{i,\mathrm{l}\mathrm{i}\mathrm{m}}^{\ast }}{{\left({\mathrm{D}}^{2}a\right)}_i^{\ast }},$$(48) $$a_{i+1/2,\mathrm{L}}\mapsto {\overline{a}}_i+(a_{i+1/2,\mathrm{L}}-{\overline{a}}_i)\frac{{\left({\mathrm{D}}^{2}a\right)}_{i,\mathrm{l}\mathrm{i}\mathrm{m}}^{\ast }}{{\left({\mathrm{D}}^{2}a\right)}_i^{\ast }}.$$(49)

If ∣(D²a)_i^*∣ < ϵ = 10⁻¹² we do not modify the reconstructed states in this step. Equations (48) and (49) are equivalent to Eq. (23) of CS08.

If the condition in Eq. (41) is not satisfied, meaning cell i is not in the vicinity of a local extremum, we use Eqs. (29)–(31) instead of Eqs. (48) and (49) to limit the reconstructed states. CS08 propose to use a slightly less restrictive limiter away from extrema (their Eq. (26)), but that limiter produces oscillations with our time-discretization scheme and we did not use it.

Just as we did in the case of PPM84, we only used the reconstructed states and not the assumed parabolic model of a(ζ) within the cell, see Sect. 2.3.4 for details. The PPM08 method is 6th-order accurate for smooth functions a(ζ) even if they are not monotonic. Our 1D experiment in Appendix B.1 confirms this. The PPM08 method requires four ghost cells at domain boundaries.

2.3.6 The PSH method

Whereas all of the previous reconstruction methods are applied to all variables in the same way, the PSH method is hybrid: it is an unlimited piecewise-sextic method for all dynamic variables combined with PPM08 for passive scalars. This allows us to eliminate certain issues that occur with methods containing limiters when applied to slow flows (see Sec. 3.1 for details) while essentially eliminating oscillations in the passive scalars, which could represent mass fractions. The piecewise sextic reconstruction assumes that within cell i a(ζ) can be described by the sextic polynomial $$a(\zeta )=\sum _{n=0}^6{c}_n{(\zeta -{\zeta }_i)}^n.$$(50)

The seven coefficients c_n are uniquely determined by the requirement that the averages of a(ζ) in cells i − 3 + n equal ā_i−3+n for n = 0, 1,..., 6. The reconstructed states are then obtained by evaluating Eq. (50) at ζ_i−1/2 and ζ_i+1/2, respectively. The resulting expressions are $$\begin{array}{rl}{a}_{i-1/2,\mathrm{R}}=& \frac{1}{420}(4{\overline{a}}_{i-3}-38{\overline{a}}_{i-2}+214{\overline{a}}_{i-1}+319{\overline{a}}_i\\ & -101{\overline{a}}_{i+1}+25{\overline{a}}_{i+2}-3{\overline{a}}_{i+3})+O\left(\mathrm{\Delta }{x}^7\right),\end{array}$$(51) $$\begin{array}{rl}{a}_{i+1/2,\mathrm{L}}=& \frac{1}{420}(-3{\overline{a}}_{i-3}+25{\overline{a}}_{i-2}-101{\overline{a}}_{i-1}+319{\overline{a}}_i\\ & +214{\overline{a}}_{i+1}-38{\overline{a}}_{i+2}+4{\overline{a}}_{i+3})+O\left(\mathrm{\Delta }{x}^7\right),\end{array}$$(52)

The method is exact wherever a(ζ) is locally a sextic polynomial, 7th-order accurate for general but smooth functions a(ζ), and it requires four ghost cells at domain boundaries.

2.4 Approximate Riemann solvers

The reconstructed pair of primitive state quantities, q_i+1/2,L,R, defines a Riemann problem at the cell interface i + 1/2, which SLH solves by means of 1D approximate Riemann solvers to obtain the face-centered value of the fluxes $\mathcal{F}$_i+1/2. We ran the tests presented in Sec. 3 using two widely popular flux functions, namely the RUSANOV and HLLC solvers. Because in this work we only focused on simulations of subsonic flows, for comparison we also built a low-dissipation version of HLLC following the approach of Minoshima & Miyoshi (2021), who modified the Harten-Lax-van Leer-Discontinuities (HLLD, Miyoshi & Kusano 2005) scheme for magnetohydrodynamics to diminish the magnitude of the numerical dissipation for low-Mach-number flows. The authors called this low-dissipation flux “LHLLD,” so, for consistency, we refer to the low-dissipation HLLC solver as “LHLLC” throughout the text⁵. In the rest of this section, we summarize the main aspects of each of these solvers and provide the implementation details whenever several choices can be made for specifying the value of a certain quantity that is needed to evaluate the numerical flux.

2.4.1 RUSANOV

The RUSANOV flux is computed by adding an upwind, numerical diffusive term proportional to the maximum wave speed at the cell interface, S_max, to every component of the central flux. The final expression for the numerical flux reads⁶ $$\mathcal{F}({\mathit{U}}_{\mathrm{L}},{\mathit{U}}_{\mathrm{R}})=\frac{1}{2}[\mathit{F}\left({\mathit{U}}_{\mathrm{L}}\right)+\mathit{F}\left({\mathit{U}}_{\mathrm{R}}\right)]-\frac{1}{2}{S}_{max}({\mathit{U}}_{\mathrm{R}}-{\mathit{U}}_{\mathrm{L}}),$$(53)

where U_L,R are the left and right sets of conserved quantities, respectively. In SLH, S_max is estimated as $$S_{\text{max }}=max(\left|u_{\mathrm{L}}\right|+c_{\mathrm{L}},\left|u_{\mathrm{R}}\right|+c_{\mathrm{R}}),$$(54)

where c = (γp/ρ)^1/2 is the sound speed. The diffusive term in Eq. (53) scales with the Mach number of the flow $\mathcal{M}$ and allows the scheme to achieve numerical stability by smearing out any discontinuity that may arise in the vector of state quantities U.

The RUSANOV solver is one of the simplest schemes that can be used to approximate the fluxes at grid cell interfaces, which makes it very efficient in terms of Floating Point Operations per Second. However, it does not take into account the complex structure of the solution arising from the Riemann problem of gas dynamics (see, e.g., Toro 2009), so the states between the two outer waves in the Riemann fan are averaged out. For this reason, this flux function is particularly diffusive for transporting contact and shear waves, which lack the self-steepening property of sound waves.

Fig. 1 Space-time diagram showing the wave structure of the HLLC solver for a Riemann problem of gas dynamics. The Sl,r are the two outer sonic waves (solid lines), while S* is the speed of the linearly degenerate contact and shear waves (dashed line) that separate the two intermediate “star” states, U*L,r.

2.4.2 HLLC

Different from RUSANOV, the HLLC solver restores the linearly degenerate contact and shear waves back to the structure of the solution of the Riemann problem (see Fig. 1). In this method, the numerical flux is chosen according to the sign of the wave speeds in the Riemann fan, $$\mathcal{F}({\mathit{U}}_{\mathrm{L}},{\mathit{U}}_{\mathrm{R}})=\{\begin{array}{l}{\mathit{F}}_{\mathrm{L}}\text{ if }0\le S_{\mathrm{L}},\\ {\mathit{F}}_{\mathrm{L}}^{\ast }\text{ if }{S}_{\mathrm{L}}<0\le S^{\ast },\\ {\mathit{F}}_{\mathrm{R}}^{\ast }\text{ if }{S}^{\ast }<0\le S_{\mathrm{R}},\\ {\mathit{F}}_{\mathrm{R}}\text{ if }{S}_{\mathrm{R}}<0.\end{array}$$(55)

While the computation of the physical fluxes F_L = F(U_L) and F_R = F(U_R) is trivial, the fluxes in the intermediate regions, F*_L and F^*_R, are obtained by solving the Rankine–Hugoniot jump conditions across the two outer sonic waves, $$S_{\mathrm{L}}{\mathit{U}}_{\mathrm{L}}^{\ast }-{\mathit{F}}_{\mathrm{L}}^{\ast }=S_{\mathrm{L}}{\mathit{U}}_{\mathrm{L}}-{\mathit{F}}_{\mathrm{L}},$$(56) $$S_{\mathrm{R}}{\mathit{U}}_{\mathrm{R}}^{\ast }-{\mathit{F}}_{\mathrm{R}}^{\ast }=S_{\mathrm{R}}{\mathit{U}}_{\mathrm{R}}-{\mathit{F}}_{\mathrm{R}}.$$(57)

Here, ${\mathit{U}}_{\mathrm{L}}^{\ast }$ and ${\mathit{U}}_{\mathrm{R}}^{\ast }$ represent the state quantities in the star regions. In order to solve Eqs. (56)–(57), proper estimates of the wave speeds S_L and S_R must be provided beforehand. SLH computes these wave speeds as $$S_{\mathrm{L}}=min(u_{\mathrm{L}},u_{\mathrm{R}})-max(c_{\mathrm{L}},c_{\mathrm{R}}),$$(58) $$S_{\mathrm{R}}=max(u_{\mathrm{L}},u_{\mathrm{R}})+max(c_{\mathrm{L}},c_{\mathrm{R}}).$$(59)

After some assumptions and algebraic manipulations (see Toro 2009), these estimates allow Eqs. (56)–(57) to be solved for the star states ${\mathit{U}}_{\mathrm{L},\mathrm{R}}^{\ast }$, $$(\rho {)}_{\mathrm{L},\mathrm{R}}^{\ast }={\tilde{\alpha }}_{\mathrm{L},\mathrm{R}},$$(60) $$\begin{array}{rl}(\rho u{)}_{\mathrm{L},\mathrm{R}}^{\ast }=& {\tilde{\alpha }}_{\mathrm{L},\mathrm{R}}{S}^{\ast },\end{array}$$(61) $$\begin{array}{rl}(\rho v{)}_{\mathrm{L},\mathrm{R}}^{\ast }=& {\tilde{\alpha }}_{\mathrm{L},\mathrm{R}}{v}_{\mathrm{L},\mathrm{R}},\end{array}$$(62) $$\begin{array}{rl}(\rho w{)}_{\mathrm{L},\mathrm{R}}^{\ast }=& {\tilde{\alpha }}_{\mathrm{L},\mathrm{R}}{w}_{\mathrm{L},\mathrm{R}},\end{array}$$(63) $$\begin{array}{rl}{\left(\rho e_{\mathrm{t}\mathrm{o}\mathrm{t}}^{\ast }\right)}_{\mathrm{L},\mathrm{R}}^{\ast }=& {\tilde{\alpha }}_{\mathrm{L},\mathrm{R}}\frac{{\left(\rho e_{\mathrm{t}\mathrm{o}\mathrm{t}}\right)}_{\mathrm{L},\mathrm{R}}}{{\rho }_{\mathrm{L},\mathrm{R}}}\\ & +{\tilde{\alpha }}_{\mathrm{L},\mathrm{R}}(S^{\ast }-u_{\mathrm{L},\mathrm{R}})[S^{\ast }+\frac{p_{\mathrm{L},\mathrm{R}}}{{\rho }_{\mathrm{L},\mathrm{R}}(S_{\mathrm{L},\mathrm{R}}-u_{\mathrm{L},\mathrm{R}})}],\end{array}$$(64) $$\begin{array}{rl}(\rho X{)}_{\mathrm{L},\mathrm{R}}^{\ast }=& {\tilde{\alpha }}_{\mathrm{L},\mathrm{R}}{X}_{\mathrm{L},\mathrm{R}},\end{array}$$(65)

with $${\tilde{\alpha }}_{\mathrm{L},\mathrm{R}}={\rho }_{\mathrm{L},\mathrm{R}}\left(\frac{S_{\mathrm{L},\mathrm{R}}-u_{\mathrm{L},\mathrm{R}}}{S_{\mathrm{L},\mathrm{R}}-S^{\ast }}\right),$$(66)

and S^* being the speed of the intermediate wave, $$S^{\ast }=\frac{p_{\mathrm{R}}-p_{\mathrm{L}}+{\rho }_{\mathrm{L}}{u}_{\mathrm{L}}(S_{\mathrm{L}}-u_{\mathrm{L}})-{\rho }_{\mathrm{R}}{u}_{\mathrm{R}}(S_{\mathrm{R}}-u_{\mathrm{R}})}{{\rho }_{\mathrm{L}}(S_{\mathrm{L}}-u_{\mathrm{L}})-{\rho }_{\mathrm{R}}(S_{\mathrm{R}}-u_{\mathrm{R}})}.$$(67)

The gas pressure is preserved across the middle wave and takes the value $$\begin{array}{rl}{p}^{\ast }=& \frac{1}{2}[p_{\mathrm{L}}+p_{\mathrm{R}}+{\rho }_{\mathrm{L}}(S_{\mathrm{L}}-u_{\mathrm{L}})(S^{\ast }-u_{\mathrm{L}})\\ & +{\rho }_{\mathrm{R}}(S_{\mathrm{R}}-u_{\mathrm{R}})(S^{\ast }-u_{\mathrm{R}})].\end{array}$$(68)

These expressions are then inserted back into Eqs. (56)–(57) to compute ${\mathit{F}}_{\mathrm{L}}^{\ast }$ and ${\mathit{F}}_{\mathrm{R}}^{\ast }$. Finally, the interface flux is chosen according to Eq. (55).

In our study, we used a variant of the original HLLC solver of Toro et al. (1994), which allowed the low-Mach correction presented in the next section to be implemented trivially into the solver. In particular, we directly evaluated the physical fluxes in the selected state of the Riemann fan, $$\mathcal{F}({\mathit{U}}_{\mathrm{L}},{\mathit{U}}_{\mathrm{R}})=\{\begin{array}{ll}{\mathit{F}}_{\mathrm{L}}& \text{ if }0\le S_{\mathrm{L}},\\ \mathit{F}\left({\mathit{U}}_{\mathrm{L}}^{\ast }\right)& \text{ if }{S}_{\mathrm{L}}<0\le S^{\ast },\\ \mathit{F}\left({\mathit{U}}_{\mathrm{R}}^{\ast }\right)& \text{ if }{S}^{\ast }<0\le S_{\mathrm{R}},\\ {\mathit{F}}_{\mathrm{R}}& \text{ if }{S}_{\mathrm{R}}<0,\end{array}$$(69)

and we computed p^* using a linearized Riemann solver for the equations of gas dynamics (see, e.g., Toro 1991), $$p^{\ast }=\frac{1}{2}(p_{\mathrm{L}}+p_{\mathrm{R}})-\frac{1}{2}\tilde{\rho }\tilde{c}(u_{\mathrm{R}}-u_{\mathrm{L}}),$$(70)

where $\tilde{c}=(c_{\mathrm{L}}+c_{\mathrm{R}})/2$ and $\tilde{\rho }=({\rho }_{\mathrm{L}}+{\rho }_{\mathrm{R}})/2$ The system in Eq. (69) is consistent with the physical fluxes in the sense that $$\mathcal{F}(\mathit{U},\mathit{U})=\mathit{F}(\mathit{U})$$(71)

and it satisfies the Rankine–Hugoniot jump conditions across the contact wave S^* as the original solver, $$S^{\ast }{\mathit{U}}_{\mathrm{R}}^{\ast }-{\mathit{F}}_{\mathrm{R}}^{\ast }=S^{\ast }{\mathit{U}}_{\mathrm{L}}^{\ast }-{\mathit{F}}_{\mathrm{L}}^{\ast }.$$(72)

Diagnostic tests of the Kelvin–Helmholtz instability problem, described in Sec. 3.1, showed that the numerical solutions computed with our modified version and the original solver of Toro et al. (1994) are virtually indistinguishable for subsonic flows. We stress, however, that the fluxes in Eq. (69) do not satisfy the jump conditions across the sonic waves S_L and S_R (see Eqs. (56) and (57)). Therefore, there is no guarantee that the resulting scheme preserves positivity of density and internal energy when the flow is nearly transonic, in which case effects of compressibility and nonlinearities can become dominant. Such a flow regime, however, was not considered in this study.

2.4.3 LHLLC

As discussed in Sect. 2.4.2, HLLC restores the intermediate, linearly degenerate waves, so it is generally more accurate than two-wave solvers such as RUSANOV or HLL (Harten et al. 1983) in simulations involving the presence of material interfaces or the propagation of entropy waves. However, the effects of the numerical dissipation introduced by HLLC on the evolution of the flow become progressively more dominant as $\mathcal{M}$ → 0, thus producing unnecessarily large diffusive errors in highly subsonic velocity regimes (see, e.g., Fleischmann et al. 2020). In our variant of HLLC, this behavior is caused by the upwind term in the expression for p^* (see Eq. (70)), $$D\left(p^{\ast }\right)=-\frac{1}{2}\tilde{\rho }\tilde{c}(u_{\mathrm{R}}-u_{\mathrm{L}}).$$(73)

This term scales with $\mathcal{M}$, which is inconsistent with the scaling of pressure fluctuations dynamically generated by subsonic flows. In fact, in the asymptotic limit $\mathcal{M}$ → 0, the solution to the compressible Euler equations approaches the incompressible regime (Guillard & Viozat 1999), in which the gas pressure is homogeneous in space except for fluctuations proportional to $\mathcal{M}$². At low Mach numbers, the numerical term in Eq. (73) can eventually become larger than the physical pressure fluctuation at the cell interface, thus yielding highly inaccurate numerical results.

In order to correct for the flawed scaling of the numerical dissipation introduced by HLLC-like methods, here we followed the approach described in Minoshima & Miyoshi (2021), who proposed to multiply the diffusive term in Eq. (73) by a factor φ proportional to the local Mach number of the flow⁷. Such a correction was originally applied to the magnetohydrodynamic solver HLLD, but it can easily be used in HLLC by setting all magnetic field components to zero, resulting in $$\varphi =\chi (2-\chi ),$$(74)

with $$\chi =min\{1,max(\frac{\left|{\mathit{V}}_{\mathrm{L}}\right|}{c_{\mathrm{L}}},\frac{\left|{\mathit{V}}_{\mathrm{R}}\right|}{c_{\mathrm{R}}})\}.$$(75)

The final expression for p* then reads $$p^{\ast }=\frac{1}{2}(p_{\mathrm{L}}+p_{\mathrm{R}})-\varphi \frac{1}{2}\tilde{\rho }\tilde{c}(u_{\mathrm{R}}-u_{\mathrm{L}}).$$(76)

The resulting upwind term in this “low-dissipation” version of the HLLC flux (LHLLC) scales with $\mathcal{M}$² when the flow is subsonic, so the ratio of the numerical diffusive term to the amplitude of pressure fluctuations is independent of $\mathcal{M}$.

We note that the same fix cannot equally be applied to the RUSANOV flux without sacrificing numerical stability. In particular, a diffusive coefficient proportional to $\mathcal{M}$² would result in too little dissipation for sound waves. This is not the case for LHLLC, in which the complex upwinding performed in Eq. (69) guarantees that the scheme remains stable for the propagation of sound waves (see also Appendix B.1).

2.5 Time discretization

Because the acoustic Courant-Friedrichs-Lewy (CFL, Courant et al. 1928) criterion on the time step becomes excessively strict in regimes of very low Mach numbers, implicit time discretization techniques are typically better suited for simulating the evolution of such slow flows (see, e.g., Viallet et al. 2011; Miczek et al. 2015; Dumbser et al. 2019). However, we recognize that most hydrodynamic codes nowadays do not have time-implicit integration capabilities, whose implementation requires a considerable effort from code developers. Thus, to make our study easily reproducible, we decided to target in our test setups, see Sects. 3.1 and 3.2, flows with Mach numbers in the range 10⁻³ ≲ $\mathcal{M}$ ≲ 10⁻¹, where simple time-explicit marching schemes are still competitive with implicit ones. In this work, explicit time integration was performed in a semi-discrete fashion, in which the cell-surface integral of the fluxes and the cell-volume integral of the source terms in Eq. (2) were first separately discretized in space whilst the system was left continuous in time according to the method of lines (see Sec. 2.2). The resulting system of ordinary differential equations (see Eq. (13)), was then solved numerically to advance the cell-volume-averaged state quantities in time. To solve Eq. (13), we used the third-order accurate, strong stability preserving (SSP) RK3 method of Shu & Osher (1988), in which the update on ${\overline{\mathit{U}}}_{i,j,k}^{(n)}$ from time t_n to time t_n+1 = t_n + Δt is performed in three stages, $$\begin{array}{rl}{\overline{\mathit{U}}}_{i,j,k}^{(1)}& ={\overline{\mathit{U}}}_{i,j,k}^{(n)}-\mathcal{R}\left({\overline{\mathit{U}}}_{i,j,k}^{(n)}\right)\mathrm{\Delta }t,\end{array}$$(77) $$\begin{array}{rl}{\overline{\mathit{U}}}_{i,j,k}^{(2)}& =\frac{3}{4}{\overline{\mathit{U}}}_{i,j,k}^{(n)}+\frac{1}{4}{\overline{\mathit{U}}}_{i,j,k}^{(1)}-\frac{1}{4}\mathcal{R}\left({\overline{\mathit{U}}}_{i,j,k}^{(1)}\right)\mathrm{\Delta }t,\end{array}$$(78) $$\begin{array}{rl}{\overline{\mathit{U}}}_{i,j,k}^{(n+1)}& =\frac{1}{3}{\overline{\mathit{U}}}_{i,j,k}^{(n)}+\frac{2}{3}{\overline{\mathit{U}}}_{i,j,k}^{(2)}-\frac{2}{3}\mathcal{R}\left({\overline{\mathit{U}}}_{i,j,k}^{(2)}\right)\mathrm{\Delta }t.\end{array}$$(79)

In particular, we computed the spatial residuals at stage s, $$\begin{array}{rl}\mathcal{R}\left({\overline{\mathit{U}}}_{i,j,k}^{(s)}\right)=& \frac{1}{\mathrm{\Delta }x}({\mathcal{F}}_{i+1/2,j,k}^{(s)}-{\mathcal{F}}_{i-1/2,j,k}^{(s)}\\ & +{\mathcal{G}}_{i,j+1/2,k}^{(s)}-{\mathcal{G}}_{i,j-1/2,k}^{(s)}\\ & +{\mathcal{H}}_{i,j,k+1/2}^{(s)}-{\mathcal{H}}_{i,j,k-1/2}^{(s)})\\ & -{\mathcal{S}}_{i,j,k}^{(s)},\end{array}$$(80)

using the numerical techniques described in Sects. 2.3 and 2.4. Finally, in order to achieve numerical stability, we limited the time step according to $$\mathrm{\Delta }t=\frac{\mathrm{C}\mathrm{F}\mathrm{L}}{N_{\text{dim }}}\underset{i,j,k}{min}\left(\frac{\mathrm{\Delta }x}{|{\mathit{V}|}_{i,j,k}+c_{i,j,k}}\right),$$(81)

where N_dim is the number of spatial dimensions. In all the tests presented in Sec. 3, we always adopted CFL = 0.8. We preferred to use a third-order accurate time stepper over less computationally expensive (but more inaccurate) methods, such as the “Midpoint rule” or SSP-RK2 (Shu & Osher 1988), so that the spatial instead of the temporal discretization would contribute most to the building up of global truncation errors.

3 Convergence properties of different Godunov-type methods

In this section, we check if the methods included in our study (described in Sec. 2) converge to the same numerical solution for several physical quantities of interest. In particular, we test a Kelvin–Helmholtz instability and a more complex setup characterized by the presence of turbulent convective flows, turbulent entrainment, and wave excitation. The results of the latter set of simulations allow us to estimate the computational cost per fixed accuracy for any given scheme, which we show in Sec. 4.

3.1 Kelvin–Helmholtz instability

We first test all of our 18 combinations of numerical schemes as described in Sects. 2.3 and 2.4 on a 2D Kelvin–Helmholtz problem with the initial condition $$\rho =\gamma ,$$(82) $$u={\mathcal{M}}_0[1-2\eta (y)],$$(83) $$v=\frac{{\mathcal{M}}_0}{10}\sin (2\pi x),$$(84) $$p=1,$$(85) $$X=\eta (y),$$(86)

where λ = 1.4 and $$\eta (y)=\{\begin{array}{ll}\frac{1}{2}\{1+\sin [16\pi (y+0.25)]\},& \text{ for }y>-\frac{9}{32}\text{ and }y<-\frac{7}{32},\\ 1,& \text{ for }y\ge -\frac{7}{32}\text{ and }y\le \frac{7}{32},\\ \frac{1}{2}\{1-\sin [16\pi (y-0.25)]\},& \text{ for }y>\frac{7}{32}\text{ and }y<\frac{9}{32},\\ 0,& \text{ otherwise. }\end{array}$$(87)

The smooth function η(y) provides a resolvable transition between layers moving in opposite horizontal directions. The initial speed of sound is unity, so $\mathcal{M}$₀ is a tunable initial Mach number of the shear flow. We discuss solutions with $\mathcal{M}$₀ = 10⁻¹, 10⁻², and 10⁻³. In Eq. (84), a smooth initial perturbation with an amplitude of $\mathcal{M}$₀/10 is included as a velocity component perpendicular to the shear flow. The computational domain, assumed to be periodic in both x and y, spans 0 ≲ x ≲ 2, −0.5 ≲ y ≲ 0.5.

The fact that the transition function η(y) between the shearing layers is smooth⁸ allows us to compute numerically converged solutions even in the absence of physical viscosity as long as the simulations are stopped before the flow field becomes chaotic (see also Robertson et al. 2010; McNally et al. 2012; Lecoanet et al. 2017; Berlok & Pfrommer 2019). Each of the two transitions spans only ¹/16 of the domain height and is poorly resolved on the coarser grids used in our tests. Therefore, we improve the accuracy of the initial cell averages that involve η(y) by averaging η(y) over 100 points uniformly distributed in the y-range covered by each cell. We measure numerical errors with respect to a reference solution computed using PSH reconstruction and the LHLLC flux function on a 8192 × 4096 grid. The solution for $\mathcal{M}$₀ = 10⁻² is shown in Fig. 2 at four points in time⁹. As the instability grows in amplitude, the sinusoidal initial perturbation is rolled up into a series of vortices. Parts of the initial shear layers are stretched and become trapped in the centers of the vortices. Other parts of the shear layers become substantially narrower. We quantify this phenomenon by computing the minimum scale height min(H_X) ≡ 1/ max(∣∇ X∣) of the passive scalar X. Figure 3 shows that this quantity drops by as much as a factor of 28 between t = 0 and t = 0.8$\mathcal{M}$₀⁻¹. At the latter point in time, the minimum scale height is only 5.75 computational cells on the 8192 × 4096 reference grid. Extremely thin and difficult-to-resolve filaments appear at even later times (see Fig. 2). Therefore, we compare the solutions at $t=0.8{\mathcal{M}}_0^{-1}$, making the problem nonlinear and challenging enough but not computationally prohibitive. The maximum Mach number in the flow field is 1.8$\mathcal{M}$₀ at this point in time.

Figure 4 compares the distributions of the passive scalar in simulations with $\mathcal{M}$₀ = 10⁻² computed on a 128 × 64 grid. The steepest gradients are strongly under-resolved on this grid (see Fig. 3), which increases the amplitude and visibility of small-scale artifacts produced by different methods. All of the six reconstruction functions lead to extremely diffusive solutions with the RUSANOV flux, although high-order methods with limiters (PPM84, PPM08, PSH) still preserve steep gradients at some places. Numerical diffusion is strongly suppressed with the HLLC flux function owing to the explicit treatment of the contact wave in HLLC. The three highest-order methods (PPM84, PPM08, PSH) reproduce the structure of the primary vortices (cf. Fig. 2) much more closely than the lower-order methods (LIN, LIN+VL, PAR). However, the two PPM methods develop secondary instabilities around the primary vortices. This effect, not present in the reference solution, occurs also with the LHLLC flux function (see also McNally et al. 2012). The secondary instabilities become the dominant source of numerical errors. We find that these instabilities tend to grow when excessive velocity shear is generated at the grid scale and there is not enough numerical dissipation to suppress their growth. Therefore, the excitation of artificial, short wavelength Kelvin–Helmholtz instabilities is favored on coarser grids, which generate larger shear at the grid scale across the poorly resolved slip line, and by less dissipative Riemann solvers. For the same reason, simulations run with the HLLC solver are more prone to developing secondary instabilities at higher rather than lower Mach numbers, as HLLC introduces less numerical diffusion into the system when modeling faster flows (see, e.g., the panels for PPM84+HLLC in Fig. 4 and Fig. D.1). As the grid is refined, the shear layers are progressively better resolved thus reducing grid scale shear and suppressing the growth of the secondary instabilities. Not surprisingly, the smooth interiors of the primary vortices are best represented with the highest-order method PSH (cf. Fig. 2). Figures D.1 and D.2, respectively, show that the differences between the three flux functions become smaller with $\mathcal{M}$₀ = 0⁻¹ and much larger with $\mathcal{M}$₀ = 10⁻³. This is expected because the amount of numerical dissipation (relative to the flow of interest) introduced by the RUSANOV and HLLC flux functions increases with decreasing Mach number of the flow, see Sect. 2.4. With $\mathcal{M}$₀ = 10⁻³ only the PSH method reproduces the basic structure of the primary vortices when combined with the HLLC flux. Comparing Fig. 4 with Fig. D.2, we see that the LHLLC flux produces results independent of the initial Mach number $\mathcal{M}$₀ because all of the flows are considerably subsonic. All of the methods we test converge to the same flow pattern with $\mathcal{M}$₀ = 10⁻², see Fig. D.7 for solutions computed on the 4096 × 2048 grid¹⁰.

The magnitude of the largest overshoot (or undershoot) in the passive scalar is shown in the insets in Figs. 4, D.1, and D.2. Thanks to their use of limiters, the LIN+VL, PPM84, PPM08, and PSH methods largely eliminate the overshoots, making the methods useful for the advection of mass fractions. We observe only some accumulation of round-off-level overshoots (up to ≈10⁻¹⁴) with LIN+VL, independently of the initial Mach number or grid resolution. The PPM84, PPM08, and PSH methods produce small but finite overshoots. At t = 0.8$\mathcal{M}$₀⁻¹, the magnitude of the largest overshoots is ≈10⁻⁴ with PPM84, ≈10⁻² with PPM08, and ≈10⁻⁶ with PSH across all grids and all initial Mach numbers. However, they only occur in a few cells or groups of cells scattered across the computational grid, so the total mass affected by this effect is negligibly small. Moreover, the amplitude of such overshoots, drops to the round-off level in many of our simulations as the resolution is increased. Unlimited methods produce substantially larger overshoots, reaching 8% and 17% for LIN and PAR, respectively. They are the result of the dispersion errors characteristic of linear schemes of 2nd or higher order, so they are nearly always present around sharp structures and affect a much larger amount of mass than the sporadic and isolated overshoots produced by PPM84, PPM08, and PSH. The benefits of the hybrid approach in PSH thus become obvious – secondary instabilities are suppressed (as compared with PPM84 and PPM08) thanks to the absence of limiters for dynamic variables while mass fractions, for which limiters are used, are well behaved.

We quantify numerical convergence of all of our methods by measuring L₁ errors with respect to the reference solution in the density ρ, kinetic energy E_k,y associated with motions perpendicular to the initial shear flow, and mass fraction X of the passive scalar. We define the relative L₁ error for any quantity q as $$L_1=\frac{1}{{\sigma }^{\text{ref }}}\frac{\sum _{i=1}^{N_x}\sum _{j=1}^{N_y}|q_{i,j}-q_{i,j}^{\text{ref }}|}{N_x{N}_y},$$(88)

where q^ref is the reference solution and the sums run over the whole computational grid of N_x × N_y cells. The error is normalised using the standard deviation σ^ref of σ^ref. This choice is motivated by the fact that density fluctuations are much smaller than the mean density in our setup. The reference solution is always re-binned from its original, 8192 × 4096 grid, by repeatedly averaging groups of 2 × 2 neighboring cells until the desired grid resolution is achieved.

The L₁ errors for the set of simulations with the initial Mach number $\mathcal{M}$₀ = 10⁻² are shown in Fig. 5. We first focus on the kinetic energy E_k,y and mass fraction X. The solutions computed using the LHLLC flux approach 2nd-order convergence, as expected. Only the PPM84 reconstruction method gives a slight decrease in the convergence rate in E_k,y on the finest of our grids. The same effect is observed in the sets of simulations with $\mathcal{M}$₀ = 10⁻¹ (Fig. D.3) and $\mathcal{M}$₀ = 10⁻³ (Fig. D.4). The convergence rate of PPM08 also decreases on very fine grids in the latter case. Apart from this, the convergence curves are nearly independent of the initial Mach number with the LHLLC flux (cf. Figs. 5, D.3, and D.4), confirming its low-Mach property. The highest-order reconstruction methods usually (but not always) produce the smallest errors on a given and sufficiently fine grid. The errors differ by up to one order of magnitude, although the overall 2D scheme is 2nd-order accurate in all of the cases.

The magnitude of numerical errors strongly increases when we decrease initial Mach number with the non-low-Mach flux functions RUSANOV and HLLC. As described above, the spatial structure of errors produced by linear methods differs from that produced by methods with limiters. Figures 5, D.3, and D.4 show that the linear methods LIN, PAR, and for dynamic variables also PSH approach 2nd-order convergence on sufficiently fine grids even if the magnitude of the errors is much larger than what we obtain with the LHLLC flux function.

The two PPM methods converge poorly in ρ at low Mach numbers (Figs. 5 and D.4). Because the initial density is constant everywhere, all of the density fluctuations are the integrated effects of the divergence (or convergence) of the velocity field in the continuity equation. The magnitude of the velocity divergence in simulations computed on the 512 × 256 grid is shown in Figs. 6, D.5, and D.6. The solutions computed using the LIN, PAR, and PSH methods clearly show the structure of the primary vortices with some oscillations in the shear layers and a background of relatively weak, large-scale sound waves. On the other hand, the solutions computed using the PPM84 and PPM08 methods show a large amount of small-scale numerical noise. After inspecting the time dependence of these artificial structures¹¹, we conclude that some of these structures travel with the flow while others have the character of small-scale sound waves. We believe that both originate from cumulative effects of the limiters switching their local state many times over the time span of the simulation¹². The switching can result in rapid changes in the magnitude of the discontinuities at cell faces and, consequently, in the amount of dissipation applied in the Riemann solver (i.e., the numerical flux function). This effect is weakest with the low-dissipation flux function LHLLC but it is still clearly present on the 4096 × 2048 grid, see Fig. D.8. The structures disappear when we disable the limiters in PPM84 and PPM08 (not shown in the figures). We do not observe the fast-propagating, small-scale sound waves with LIN+VL but the method does produce thin structures in the velocity divergence around the primary vortices, which follow the flow.

Overall, the PSH method produces by far the smallest errors with the RUSANOV and HLLC flux functions and initial Mach numbers $\mathcal{M}$₀ ≤ 10⁻² (Figs. 5 and D.4). The only exception is the poor convergence (or even divergence) in the passive tracer X observed with the PSH+RUSANOV combination in simulations with $\mathcal{M}$₀ ≤ 10⁻² performed on very fine grids. This effect is also likely related to the limiters switching their local state many times over the time span of the simulation but we did not investigate it further.

Fig. 2 Reference solution to the Kelvin–Helmholtz problem with initial Mach number $\mathcal{M}$0 = 10−2. The solution was computed using PSH reconstruction and the LHLLC flux function on an 8192 × 4096 grid. The mass fraction X of the passive scalar is shown at four points in time: late in the linear growth of the instability (t = 0.2 $\mathcal{M}$0−1), at an early stage of nonlinear evolution (t = 0.4 $\mathcal{M}$0−1), at a stage when the primary vortices have fully formed (t = 0.8 $\mathcal{M}$0−1; the final time for all of our other Kelvin–Helmholtz simulations), and at a late stage when fine threads have formed inside the primary vortices (t = 1.6 $\mathcal{M}$0−1). We use the same color scale as in Figs. 4, D.2, and D.1, although 0 ≲ X & ≲ 1 in the reference solution.

Fig. 3 Resolution dependence of the minimum scale height min(Hx) of the passive scalar in the Kelvin–Helmholtz problem expressed in units of the computational cell width and shown at five points in time. Grid resolution is given by the number Nx of computational cells along the x axis. The initial Mach number is $\mathcal{M}$0 = 10−2 and we use PSH reconstruction and the LHLLC flux function in this series of simulations. Once the steepest gradients in X become resolved the minimum scale height starts to follow the linear scaling relations shown.

Fig. 4 Distributions of the mass fraction X of the passive scalar in simulations of the Kelvin–Helmholtz instability with the initial Mach number $\mathcal{M}$0 = 10−2 on a 128 × 64 grid. The results are shown at t = 0.8 $\mathcal{M}$0−1. Rows and columns show different reconstruction schemes and numerical flux functions, respectively. The color scale intentionally shows values X < 0 and X > 1 to highlight the overshoots that some of the schemes produce. The absolute value of the largest overshoot outside of this range is given in each panel.

3.2 3D simulations of convection, turbulent entrainment, and wave excitation

In this section, we consider a test problem in which a variety of complex hydrodynamic phenomena can be captured on the same computational grid. These include convective transport of energy, turbulent mixing at a convective boundary, and the propagation of internal waves in a stably stratified layer, all of which are often encountered in simulations of geophysical and stellar hydrodynamics. The initial conditions for this test are adopted from the work of Andrassy et al. (2022) and they represent a hydrostatic solution of Eq. (2). In particular, the thermodynamic conditions in this setup are similar to those found in an oxygen-burning shell of a massive star. We map the hydrostatic stratification (see Fig. 7) on an evenly spaced, 3D, Cartesian grid defined by (x, y, z) ∈ [−1, 1] × [1, 3] × [−1, 1]. The gravitational acceleration points downward in the y-direction, $$g_y=g_0{f}_{\mathrm{g}}(y)y^{-5/4},$$(89)

where g₀ = −1.414870 and $$f_{\mathrm{g}}(y)=\{\begin{array}{ll}\frac{1}{2}\{1+\sin \left[16\pi (y-\frac{1}{32})\right]\},& \text{ for }1\le y<1+\frac{1}{16},\\ 1,& \text{ for }1+\frac{1}{16}\le y\le 3-\frac{1}{16}\\ \frac{1}{2}\{1+\sin \left[16\pi (y-\frac{1}{32})\right]\},& \text{ for }3-\frac{1}{16}<y\le 3.\end{array}$$(90)

As in Andrassy et al. (2022), we work with rescaled quantities such that the speed of sound and the density at the base of the box are taken as units of velocity and density, and the thickness of the convective oxygen shell as unit of length. A detailed list of units with the associated values can be found in Table 1 of Andrassy et al. (2022).

The initially hydrostatic stratification is described by a piecewise-polytropic relation in the form $$\frac{\mathrm{\partial }\ln p}{\mathrm{\partial }\ln \rho }=\{\begin{array}{ll}{\gamma }_0,& \text{ for }1\le y<2-\frac{1}{16},\\ {\gamma }_0+\eta (y)({\gamma }_1-{\gamma }_0),& \text{ for }2-\frac{1}{16}\le y\le 2+\frac{1}{16},\\ {\gamma }_1,& \text{ for }2+\frac{1}{16}<y\le 3,\end{array}$$(91)

where γ₀ = 5/3, γ₁ = 1.3, and η(y) is a smooth function, $$\eta (y)=\{\begin{array}{ll}0,& \text{ for }1\le y<2-\frac{1}{16},\\ \frac{1}{2}[1+\sin (8\pi y)],& \text{ for }2-\frac{1}{16}\le y\le 2+\frac{1}{16},\\ 1,& \text{ for }2+\frac{1}{16}<y\le 3.\end{array}$$(92)