Issue 
A&A
Volume 686, June 2024



Article Number  A34  
Number of page(s)  38  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202348882  
Published online  28 May 2024 
Performance of highorder Godunovtype methods in simulations of astrophysical low Mach number flows
^{1}
Heidelberger Institut für Theoretische Studien,
SchlossWolfsbrunnenweg 35,
69118
Heidelberg,
Germany
email: giovanni.leidi@hits.org
^{2}
Bordeaux Institute of Mathematics, Bordeaux University and CNRS/UMR5251,
33405
Talence,
France
^{3}
HighPerformance Computing Center Stuttgart,
Nobelstraße 19,
70569
Stuttgart,
Germany
^{4}
Computer, Computational and Statistical Sciences (CCS) Division and Center for Theoretical Astrophysics (CTA), Los Alamos National Laboratory,
PO Box 1663,
Los Alamos,
NM 87545,
USA
^{5}
Zentrum für Astronomie der Universität Heidelberg, Institut für Theoretische Astrophysik,
Philosophenweg 12,
69120
Heidelberg,
Germany
Received:
7
December
2023
Accepted:
23
February
2024
Highorder Godunov methods for gas dynamics have become a standard tool for simulating different classes of astrophysical flows. Their accuracy is mostly determined by the spatial interpolant used to reconstruct the pair of Riemann states at cell interfaces and by the Riemann solver that computes the interface fluxes. In most Godunovtype methods, these two steps can be treated independently, so that many different schemes can in principle be built from the same numerical framework. Because astrophysical simulations often test out the limits of what is feasible with the computational resources available, it is essential to find the scheme that produces the numerical solution with the desired accuracy at the lowest computational cost. However, establishing the best combination of numerical options in a Godunovtype method to be used for simulating a complex hydrodynamic problem is a nontrivial task. In fact, formally more accurate schemes do not always outperform simpler and more diffusive methods, especially if sharp gradients are present in the flow. For this work, we used our fully compressible SevenLeague Hydro (SLH) code to test the accuracy of six reconstruction methods and three approximate Riemann solvers on two and threedimensional (2D and 3D) problems involving subsonic flows only. We considered Mach numbers in the range from 10^{−3} to 10^{−1}, which are characteristic of many stellar and geophysical flows. In particular, we considered a wellposed, 2D, Kelvin–Helmholtz instability problem and a 3D turbulent convection zone that excites internal gravity waves in an overlying stable layer. Although the different combinations of numerical methods converge to the same solution with increasing grid resolution for most of the quantities analyzed here, we find that (i) there is a spread of almost four orders of magnitude in computational cost per fixed accuracy between the methods tested in this study, with the most performant method being a combination of a lowdissipation Riemann solver and a sextic reconstruction scheme; (ii) the lowdissipation solver always outperforms conventional Riemann solvers on a fixed grid when the reconstruction scheme is kept the same; (iii) in simulations of turbulent flows, increasing the order of spatial reconstruction reduces the characteristic dissipation length scale achieved on a given grid even if the overall scheme is only second order accurate; (iv) reconstruction methods based on slopelimiting techniques tend to generate artificial, highfrequency acoustic waves during the evolution of the flow; and (v) unlimited reconstruction methods introduce oscillations in the thermal stratification near the convective boundary, where the entropy gradient is steep.
Key words: convection / hydrodynamics / instabilities / turbulence / waves / methods: numerical
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
Highresolution 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 higherorder extensions of the original firstorder accurate method of Godunov (1959) and their timeintegration algorithm is typically carried out in three steps. First, a pair of Riemann states is reconstructed at each grid cell interface by applying highorder monotonic interpolants to a set of cellaveraged 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 cellvolumeaveraged state quantities to be advanced in time^{1}.
In most highorder 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, (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 Godunovlike 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 problemdependent. 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 higherthansecondorder 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 ensembleaveraged 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 Godunovtype 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 Godunovtype 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 highresolution 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), magnetorotational instabilities in accretion disks (Flock et al. 2010), RichtmyerMeshkov 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 Godunovtype 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 (: = ∣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 lowMachnumber 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 SevenLeague 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^{−3} ≲ ≲ 10^{−1}). In particular, we considered a twodimensional (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 Godunovtype 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 (2)
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 = [F∣G∣H] is a tensor defined by the flux vectors (4)
Here, ρ denotes the mass density, V = (u, v, w) the velocity field, 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, (5)
In this work, we only considered a perfect gas with a given adiabatic index γ, for which (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, (7)
must be included in the righthandside term of Eq. (2). Because we assumed g to be timeindependent, 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, (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 lowMachnumber 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 SevenLeague Hydro code
In our study, Eq. (2) was solved numerically using the SevenLeague 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 finitevolume discretization on an arbitrarily curvilinear, but logically rectangular, Eulerian grid to retain the conservation properties of the fluiddynamics equations. Numerical solutions to Eq. (2) are computed by means of Godunovtype 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 wellknown Rusanov (Rusanov 1962), Roe (Roe 1981), and HartenLaxvan LeerContact (HLLC, Toro et al. 1994) approximate Riemann solvers, SLH adopts special lowMachnumber methods (Liou 2006; Miczek et al. 2015; Minoshima & Miyoshi 2021) to reduce the excessive numerical dissipation introduced by shockcapturing 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 (firstorder 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 wellbalancing 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 cellvolumeaveraged source term and the cellsurfaceaveraged fluxes are approximated using the midpoint method, making the code at best secondorder accurate in space. On a 3D, evenly spaced Cartesian grid, the final expressions for these integrals read (9) (10) (11) (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 _{i}_{+1/2,j,k} refer to the facecentered 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^{3} and ΔA = Δx^{2}, 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), (13)
The time update on the cellvolumeaveraged conserved variables, , 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^{2}. 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 wellbehaved 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 totalvariationdiminishing (TVD) property and eliminates the oscillations completely^{3}. As examples of higherorder methods with limiters, we included two versions of the widely used piecewiseparabolic 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 lowerorder 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 and , 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 onedimensional. 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 cellvolumeaveraged primitive variables , where q = (ρ, V, p, X), because it helps reducing oscillations near discontinuities as compared to reconstructing cellvolumeaveraged conserved quantities . In SLH, transformations between and are performed using 2ndorder approximations, (14)
where ξ is a nonlinear, invertible transformation, (15)
and U_{i} is the i th component of U.
Overview of spatial reconstruction methods used in this work.
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 cellindex coordinates) is estimated using the central difference (the Fromm method) (16)
using which we obtain the reconstructed states (17) (18)
The LIN method is exact wherever a(ζ) is locally linear, 2ndorder accurate for general but smooth functions a(ζ), and it requires two ghost cells at domain boundaries^{4}.
2.3.2 The LIN+VL method
The vanLeerlimited linear reconstruction method is particularly easy to describe in terms of the 2ndorderaccurate, interfacecentered slopes (19) (20)
The final slope δ_{i}_{,lim} is then obtained by applying the limiter of van Leer (1974), (21)
which gives the reconstructed states (22) (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 2ndorder 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 (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 ā_{i1}, ā_{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 (25) (26)
The PAR method is exact wherever a(ζ) is locally parabolic, 3rdorder 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 twostep process. The first is based on the 4thorderaccurate interpolation formula (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 (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,li}_{m}, δ_{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_{i1/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 (29) (30)
The limiter is then defined by the variable replacements (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 RungeKutta 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 4thorder 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 Gaussianshaped 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 6thorderaccurate interpolation formula (cf. Eq. (17) of CS08) (32) (33)
If a_{i}_{+1/2} does not satisfy the condition (cf. Eq. (13) of CS08) (34)
a limiter is applied. It is based on three 2ndorderaccurate second derivatives, (35) (36) (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^{2}a)_{i}_{+1/2,lim} such that (38)
if all three derivatives have the same sign and (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}, (40)
Equation (19) of CS08 is equivalent to our Eq. (40) except for the factor in front of the second term, which is in CS08. This is likely a typographical error because the factor of 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) (41)
cell i is close to a local extremum, which should be preserved if smooth enough. The second derivatives (42) (43) (44) (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 (46)
if all four second derivatives have the same sign and (47)
otherwise. Finally, the reconstructed states are updated, (48) (49)
If ∣(D^{2}a)_{i}^{*}∣ < ϵ = 10^{−12} 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 timediscretization 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 6thorder 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 piecewisesextic 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 (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 (51) (52)
The method is exact wherever a(ζ) is locally a sextic polynomial, 7thorder 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 facecentered value of the fluxes _{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 lowdissipation version of HLLC following the approach of Minoshima & Miyoshi (2021), who modified the HartenLaxvan LeerDiscontinuities (HLLD, Miyoshi & Kusano 2005) scheme for magnetohydrodynamics to diminish the magnitude of the numerical dissipation for lowMachnumber flows. The authors called this lowdissipation flux “LHLLD,” so, for consistency, we refer to the lowdissipation HLLC solver as “LHLLC” throughout the text^{5}. 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^{6} (53)
where U_{L,R} are the left and right sets of conserved quantities, respectively. In SLH, S_{max} is estimated as (54)
where c = (γp/ρ)^{1/2} is the sound speed. The diffusive term in Eq. (53) scales with the Mach number of the flow 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 selfsteepening property of sound waves.
Fig. 1 Spacetime diagram showing the wave structure of the HLLC solver for a Riemann problem of gas dynamics. The S_{l,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, (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, (56) (57)
Here, and 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 (58) (59)
After some assumptions and algebraic manipulations (see Toro 2009), these estimates allow Eqs. (56)–(57) to be solved for the star states , (60) (61) (62) (63) (64) (65)
and S^{*} being the speed of the intermediate wave, (67)
The gas pressure is preserved across the middle wave and takes the value (68)
These expressions are then inserted back into Eqs. (56)–(57) to compute and . 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 lowMach 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, (69)
and we computed p^{*} using a linearized Riemann solver for the equations of gas dynamics (see, e.g., Toro 1991), (70)
where and The system in Eq. (69) is consistent with the physical fluxes in the sense that (71)
and it satisfies the Rankine–Hugoniot jump conditions across the contact wave S^{*} as the original solver, (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 twowave 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 → 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)), (73)
This term scales with , which is inconsistent with the scaling of pressure fluctuations dynamically generated by subsonic flows. In fact, in the asymptotic limit → 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 ^{2}. 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 HLLClike 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^{7}. 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 (74)
The final expression for p* then reads (76)
The resulting upwind term in this “lowdissipation” version of the HLLC flux (LHLLC) scales with ^{2} when the flow is subsonic, so the ratio of the numerical diffusive term to the amplitude of pressure fluctuations is independent of .
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 ^{2} 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 CourantFriedrichsLewy (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 timeimplicit 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^{−3} ≲ ≲ 10^{−1}, where simple timeexplicit marching schemes are still competitive with implicit ones. In this work, explicit time integration was performed in a semidiscrete fashion, in which the cellsurface integral of the fluxes and the cellvolume 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 cellvolumeaveraged state quantities in time. To solve Eq. (13), we used the thirdorder accurate, strong stability preserving (SSP) RK3 method of Shu & Osher (1988), in which the update on from time t_{n} to time t_{n}_{+1} = t_{n} + Δt is performed in three stages, (77) (78) (79)
In particular, we computed the spatial residuals at stage s, (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 (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 thirdorder accurate time stepper over less computationally expensive (but more inaccurate) methods, such as the “Midpoint rule” or SSPRK2 (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 Godunovtype 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 (82) (83) (84) (85) (86)
The smooth function η(y) provides a resolvable transition between layers moving in opposite horizontal directions. The initial speed of sound is unity, so _{0} is a tunable initial Mach number of the shear flow. We discuss solutions with _{0} = 10^{−1}, 10^{−2}, and 10^{−3}. In Eq. (84), a smooth initial perturbation with an amplitude of _{0}/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^{8} 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 ^{1}/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 yrange 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 _{0} = 10^{−2} is shown in Fig. 2 at four points in time^{9}. 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_{0}^{−1}. 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 difficulttoresolve filaments appear at even later times (see Fig. 2). Therefore, we compare the solutions at , making the problem nonlinear and challenging enough but not computationally prohibitive. The maximum Mach number in the flow field is 1.8_{0} at this point in time.
Figure 4 compares the distributions of the passive scalar in simulations with _{0} = 10^{−2} computed on a 128 × 64 grid. The steepest gradients are strongly underresolved on this grid (see Fig. 3), which increases the amplitude and visibility of smallscale artifacts produced by different methods. All of the six reconstruction functions lead to extremely diffusive solutions with the RUSANOV flux, although highorder 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 highestorder methods (PPM84, PPM08, PSH) reproduce the structure of the primary vortices (cf. Fig. 2) much more closely than the lowerorder 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 highestorder method PSH (cf. Fig. 2). Figures D.1 and D.2, respectively, show that the differences between the three flux functions become smaller with _{0} = 0^{−1} and much larger with _{0} = 10^{−3}. 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 _{0} = 10^{−3} 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 _{0} because all of the flows are considerably subsonic. All of the methods we test converge to the same flow pattern with _{0} = 10^{−2}, see Fig. D.7 for solutions computed on the 4096 × 2048 grid^{10}.
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 roundofflevel overshoots (up to ≈10^{−14}) 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_{0}^{−1}, the magnitude of the largest overshoots is ≈10^{−4} with PPM84, ≈10^{−2} with PPM08, and ≈10^{−6} 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 roundoff 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_{1} 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_{1} error for any quantity q as (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 rebinned from its original, 8192 × 4096 grid, by repeatedly averaging groups of 2 × 2 neighboring cells until the desired grid resolution is achieved.
The L_{1} errors for the set of simulations with the initial Mach number _{0} = 10^{−2} 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 2ndorder 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 _{0} = 10^{−1} (Fig. D.3) and _{0} = 10^{−3} (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 lowMach property. The highestorder 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 2ndorder accurate in all of the cases.
The magnitude of numerical errors strongly increases when we decrease initial Mach number with the nonlowMach 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 2ndorder 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, largescale sound waves. On the other hand, the solutions computed using the PPM84 and PPM08 methods show a large amount of smallscale numerical noise. After inspecting the time dependence of these artificial structures^{11}, we conclude that some of these structures travel with the flow while others have the character of smallscale 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^{12}. 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 lowdissipation 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 fastpropagating, smallscale 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 _{0} ≤ 10^{−2} (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 _{0} ≤ 10^{−2} 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 _{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 _{0}^{−1}), at an early stage of nonlinear evolution (t = 0.4 _{0}^{−1}), at a stage when the primary vortices have fully formed (t = 0.8 _{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 _{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(H_{x}) 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 N_{x} of computational cells along the x axis. The initial Mach number is _{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 _{0} = 10^{−2} on a 128 × 64 grid. The results are shown at t = 0.8 _{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 oxygenburning 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 ydirection, (89)
where g_{0} = −1.414870 and (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 piecewisepolytropic relation in the form (91)
where γ_{0} = 5/3, γ_{1} = 1.3, and η(y) is a smooth function, (92)
The adiabatic index is γ= γ_{0}. In this work, we assume periodic boundary conditions in the horizontal x and zdirection, whereas reflecting, stressfree boundaries are used in the vertical ydirection. This problem is set up such that turbulent convective flows develop in the lower half of the domain, which is initially adiabatic, while internal gravity waves are free to propagate in the upper, stably stratified layer. To keep track of the position of the upper convective boundary, at t = 0 we fill the stable layer with a passive scalar whose abundance smoothly turns to 0 across the upper boundary of the adiabatic region according to X = η(y).
In order to drive the convection, we include a timeindependent heat source that continuously injects energy into the system close to the base of the box. The rate of energy released per unit volume, (93)
is added to the righthand side of Eq. (2) as (94)
In the work of Andrassy et al. (2022), the amplitude of the heat source was = 3.795720 × 10^{−4}, which gave rise to convective flows with a rootmeansquare Mach number _{rms} ≈ 0.04. To make the problem more challenging, here we decrease by a factor of ten, so = 3.795720 × 10^{−5}. The lower heating rate, according to the wellestablished relation (see, e.g., Woodward et al. 2014; Käpylä 2021; Horst et al. 2021), should drive convection at _{rms} ≈ 0.02.
In this setup, internal gravity waves (IGWs) are excited by the interaction of the convective flows with the bottom boundary of the stably stratified layer. The wavelength of IGWs in the direction of gravity becomes shorter when the waves are excited at progressively lower temporal frequencies (Sutherland 2010). At the heating rate we consider, the most prominent IGWs that originate at the convective boundary are only barely spatially resolved on the coarsest of our grids with 128^{3} cells. Therefore, we decide not to decrease even further for these simulations because it would give rise to convective flows with lower characteristic frequencies and lead to the generation of unresolved IGWs in the stable layer. Due to the fully compressible nature of SLH, we also expect shortwavelength sound waves to be generated, although at much lower amplitudes than those of IGWs at the typical Mach numbers encountered in this test problem (Lighthill 1952).
As done for the test described in Sec. 3.1, here we run simulations for each of the 18 considered combinations of Riemann solvers and spatial reconstruction schemes. To judge the numerical convergence of our results, each combination of methods is run on grids with 128^{3} and 256^{3} cells. Additionally, we run a single simulation on a 512^{3} grid using the LHLLC Riemann solver and the PAR reconstruction scheme, which we consider the reference solution for this test problem. However, because of the chaotic nature of the turbulent flows that arise in the convective layer, convergence is not expected in the exact flow morphology, so we do not compute L_{1} errors as done in Sect. 3.1. Instead, we analyze the convergence of the numerical results in terms of ensembleaveraged quantities that are representative of the dynamical properties of the system, such as kinetic energy spectra computed in the convective and stable layers.
To break the initial symmetry, we add a perturbation to the hydrostatic density stratification in the form (95)
The subsequent evolution of the system and the development of convection is shown in Fig. 8. The density perturbation, alongside the action of the heat source, generate packets of fluid with higher entropy content than the adiabatic surroundings. The packets of hot and lowdensity material buoyantly rise in the adiabatic stratification until they reach y ≈ 2. At this height, the temperature stratification turns subadiabatic and the buoyant force acting on the rising plumes changes sign, forcing them to overturn. IGWs excited at the bottom boundary of the subadiabatic region propagate upward in the stratification (with characteristic Mach numbers in the range from 0.005 to 0.01) and are subsequently reflected at the top boundary of the domain. Shear instabilities break the largescale buoyant structures that arise in the adiabatic layer and initiate the cascade of kinetic energy toward smaller scales. Turbulent convection fully develops after approximately one convective turnover time scale, τ_{conv} = 133 time units, which we define according to (96)
In Eq. (96), the rootmeansquare convective speed is averaged over several convective turnover time scales and L_{conv} = 1 is taken as representative of the vertical extent of the convection zone. In our reference simulation, the rootmeansquare Mach number in the convection zone is^{13} _{rms} = 0.019 ± 0.001, which is in agreement with the value predicted by the scaling relation. All simulations are run until τ_{max} = 32τ_{conv} to have a proper coverage of the dynamical evolution of the system and to compute meaningful time averages needed for the following analysis.
Fig. 5 Relative L_{1} errors for different variables (rows), flux functions (columns), and reconstruction methods (legend) as functions of the grid resolution in the simulations of the Kelvin–Helmholtz instability with the initial Mach number _{0} = 10^{−2}. The dashed and dotted lines, which are at the same locations in all of the panels, show the 1st and 2ndorder scalings to guide the eye. 
Fig. 6 Magnitude of velocity divergence in simulations of the Kelvin–Helmholtz instability with the initial Mach number _{0} = 10^{−2} on a 512 × 256 grid. The results are shown at t = 0.8 _{0}^{−1}. Rows and columns show different reconstruction schemes and numerical flux functions, respectively. The smallscale structure produced by the LIN+VL, PPM84, and PPM08 schemes, which does not occur with the LIN, PAR, and PSH schemes, is caused by the use of limiters. 
Fig. 7 Hydrostatic stratification of density, gas pressure, pseudoentropy (p/ρ^{γ}), and gravity at t = 0 in the problem of turbulent convection and wave excitation at a convective boundary, based on the initial conditions provided in Andrassy et al. (2022). The distribution of the mass fractional abundance of a passive tracer (X) and the heat source () are also shown, here represented by a dashdotted pink and dashed black line, respectively. In this setup, the gravitational acceleration smoothly turns to zero at the bottom and top boundaries of the spatial domain to make the problem consistent with the usage of reflecting, stress free boundary conditions, which for the gas pressure imply ∂p/∂y = 0. 
3.2.1 Flow morphology and spatial kinetic energy spectra
In Fig. 9, we show snapshots of the Mach number taken at the vertical plane z = 0, obtained at the final time t = 32τ_{conv} with all the methods tested here. Although all panels show results computed using the same grid resolution (256^{3}), there is a vast difference between the methods in terms of effects induced by numerical dissipation. With the RUSANOV flux, convection mainly happens in the form of largescale plumes and IGWs are barely excited in the stable layer except when highorder methods such as PSH or PPM08 are used to reconstruct the Riemann states. Smallerscale structures typical of turbulent flows become progressively more visible in the convection zone with higherorder reconstruction methods and less dissipative solvers. Internal waves with progressively shorter wavelengths also tend to be more visible in the stable layer with less dissipative methods.
The effects of the numerical dissipation on the properties of the flow are better quantified in Fig. 10 and 11, where we show the power spectra of kinetic energy extracted from a horizontal plane in the middle of the convection zone (y = 1.5) and stable layer (y = 2.5) as functions of the horizontal wavenumber, . We compute k_{x} and k_{z} as (97) (98)
where ⌊.⌋ represents the floor function and N_{x} and N_{z} are the numbers of cells in the x– and z–direction, respectively. The spectra are averaged over the time interval t ∈ (10τ_{conv}, 32τ_{conv}). In the convection zone, all of the curves approximately agree with the Kolmogorov scaling law (k_{h}^{−5/3}) on some intermediate scales. We notice that the kinetic energy spectra shown in Fig. 10 present a pileup of kinetic energy at the bottom of the inertial range. This phenomenon, also known as “bottleneck effect,” is often observed both in hydrodynamic simulations (see, e.g., Dobler et al. 2003) and experimental studies (see, e.g., Küchler et al. 2019) of turbulent flows. The extent of the inertial range of the turbulent kinetic energy spectrum greatly differs from method to method. With a fixed spatial reconstruction scheme, the kinetic energy spectrum obtained with LHLLC enters the viscous subrange (where the kinetic energy of the turbulent eddies is dissipated into internal energy of the gas) at higher wavenumbers as compared to both HLLC and RUSANOV. Reconstruction schemes based on slope limiters are characterized by shallower spectra in the viscous subrange than those generated by unlimited schemes. In the stable layer (see Fig. 11), RUSANOV and HLLC generate much weaker IGWs than the reference run even on large scales, except when used in combination with PSH. Both in the convective and stable layers, all spectra converge to the reference solution if the order of the spatial reconstruction method or the grid resolution are increased.
Fig. 8 Pseudoentropy (A = p/ρ^{γ}) fluctuations with respect to the hydrostatic, background state computed at different times during the development of the convection, as indicated by the insets, in the problem of turbulent convection and wave excitation. Data are visualized in the z = 0 plane. This simulation was run using the PAR method and LHLLC on a 512^{3} grid and we consider it our reference solution. Convective plumes start exciting IGWs at the upper convective boundary by t ≈ τ_{conv} (see center bottom panel). 
3.2.2 Turbulent entrainment at the convective boundary
The overturning of the convective flows at the upper convective boundary gives rise to a variety of hydrodynamic processes (including shear instabilities, convective overshoot, and breaking of surface waves) that lead to the entrainment of material from the upper, stably stratified layer into the convection zone. The conglomeration of these processes is known in the literature as “convective boundary mixing” (see, e.g., Meakin & Arnett 2007; Andrassy et al. 2020; Horst et al. 2021 ; Anders & Pedersen 2023) and it increases the size of convective regions over time. As part of our comparison study, we quantify the impact of the choice of a specific combination of methods in Godunovtype schemes on the evolution of the convective boundary. At any given point in time, we assume that the vertical coordinate of the boundary between the convective and stable layer, y_{cb}, is the position at which the horizontal average of the passive scalar X takes the value (99)
where X_{min} and X_{max} are the minimum and maximum values of X, respectively. Such a choice for y_{cb} is justified by the fact that the initial abundance of the passive scalar, X, smoothly transitions from 0 to 1 at y ≈ 2, which is the boundary of the initially adiabatic layer where convection sets in first. The time evolution of y_{cb} is shown in Fig. 12 for all of our combinations of methods. We note that PPMbased methods used in combination with RUSANOV do not show signs of numerical converge to the reference solution. Also, the other reconstruction schemes used with RUSANOV seem to accelerate the mixing at the convective boundary with respect to the reference run, especially on the 128^{3} grid. In reality, this is an artifact of the method used to estimate the position of the convective boundary y_{cb}. In fact, because RUSANOV does not resolve shear or contact waves (see Sec. 2.4), the initial gradient in X at the convective boundary is further smoothed out by the effects of numerical dissipation. When convection fully develops, it mixes the passive scalar that has diffused inside the convection zone and homogenizes its abundance, thus shifting the formal position of X = X outward. This effect is further enhanced with the most diffusive reconstruction methods tested here. Therefore, with the RUSANOV solver, entrainment of material from the stable to the convective layer is mostly determined by numerical diffusion rather than turbulent mixing and the distribution of the passive tracer is not representative of y_{cb}. On the other hand, all of the other methods tested here quickly achieve numerical converge to the reference solution with increasing grid resolution, the order of the reconstruction scheme, or by using progressively less dissipative Riemann solvers. Among the data sets obtained with the six reconstruction schemes, the maximum relative deviation of y_{cb} from the reference solution at the final time is 20% and 15% with HLLC on the 128^{3} and 256^{3} grids, respectively, while with LHLLC it decreases from 12% to 5% with increasing grid resolution.
Fig. 9 Snapshots obtained at the final simulation time (t = 32τ_{conv}) showing the distribution of the Mach number at z = 0 in the simulations of turbulent convection and wave propagation in a 3D box. The grid resolution is 256^{3}. Each panel is cut between y = 1.3 and y = 2.7 for a better visualization of the flows in the proximity of the upper convective boundary. 
Fig. 10 Kinetic energy spectra computed on the horizontal plane y = 1.5 of the convective layer for all of our 18 combinations of Riemann solvers and spatial reconstruction schemes in the problem of turbulent convection and wave excitation. The spectra have been averaged over the time interval t ∈ (10τ_{conv}, 32τ_{conv}) and rescaled by the Kolmogorov law (k_{h}^{−5/3}) and by the value Ê_{0} of the spectral energy density of the reference run at k_{h} = 15. 
3.2.3 Shape of the convective boundary
The properties of the upper convective boundary are also studied by means of horizontal averages in the BruntVäisälä frequency, which determines both the spectrum of convectively driven IGWs (see, e.g., Sutherland 2010) and the rate of mass entrainment from the stable layer according to the bulkRichardsonnumber mixing model (see, e.g., Meakin & Arnett 2007; Andrassy et al. 2020; Rizzuti et al. 2023). Here, the BruntVäisälä frequency is computed as (100)
T = p/ρ is the temperature of the gas^{14}, H_{p} = − (d ln p/dy)^{−1} is the local pressure scale height, and ∇_{ad} = 0.4 is the adiabatic temperature gradient. We average the profiles of N_{BV} over t ∈ (30τ_{conv}, 32τ_{conv}) so that the change in the position of the convective boundary is smaller than its width. At the same time, the chosen averaging time interval is long enough to partly suppress the statistical fluctuations induced by the turbulent nature of the convective flows. The results are shown in Fig. 13. As convection sets in and entrains material from the upper, stably stratified layer, it steepens the vertical entropy gradient ds/dy across the convective boundary. Therefore, due to the proportionality (see, e.g., Maeder 2009), the vertical profile of the BruntVäisälä frequency develops a spike in the region close to the convective boundary and it quickly approaches in the rest of the stable layer. The fact that, in this setup, N_{BV} is characterized by large vertical gradients near the convective boundary (at least with HLLC and LHLLC), makes this quantity particularly suited to measure the amplitude of numerical undershoots or overshoots generated by the methods included in our study. In the regions right above the position of the spike in N_{BV}, limited reconstruction methods (i.e., LIN + VL, PPM08, and PPM08) do not generate undershoots below the curve except for dynamical fluctuations smaller than 2% induced by the flows in the stable layer. On the other hand, reconstruction schemes that do not use limiters for dynamic variables (i.e., LIN, PAR, and PSH) produce undershoots whose amplitude is considerably larger than that of the dynamical fluctuations. The values of the largest undershoots relative to are shown in the insets of Fig. 13. Overall, the amplitude of the undershoots increases with the order of the unlimited reconstruction method, and it ranges from ≈ 7% with LIN to almost 40% with PSH. The profiles of N_{BV} obtained on the 128^{3} grid, not shown here (but available on Zenodo^{15}), reveal that the amplitude of the undershoots slightly increases from the 128^{3} to the 256^{3} grid for LIN and PAR, although the undershoot region becomes narrower on progressively finer grids. We note that, despite the largeamplitude undershoots generated by some of our unlimited reconstruction methods, such numerical errors do not have any significant impact on the growth of the convection zone, as shown in Fig. 12. Nonetheless, they could still affect the frequency spectrum of IGWs in the stable layer, although the spatial spectra from the stable layer (see Fig. 11) do show any such effect and converge to the reference solution. If numerical undershoots are largely to be avoided whilst generating a steep convective boundary, PPM84 + LHLLC or PPM08 + LHLLC could be the methods of choice.
Fig. 11 Kinetic energy spectra computed on the horizontal plane y = 2.5 of the stable layer for all of our 18 combinations of Riemann solvers and spatial reconstruction schemes in the problem of turbulent convection and wave excitation. The spectra have been averaged over the time interval t ∈ (10τ_{conv}, 32τ_{conv}). The black solid line represents the reference run (REF), which is computed using the LHLLC Riemann solver and the PAR reconstruction scheme on a 512^{3} grid. 
Fig. 12 Time evolution of the position of the upper convective boundary (y_{cb}) in the problem of turbulent convection and entrainment described in Sec. 3.2 The top and bottom rows show results from 128^{3} and 256^{3} simulations, respectively. The black solid line represents the reference run (REF). 
3.3 Frequency power spectra
Finally, in Fig. 14, we show the frequency power spectrum of the vertical velocity component, υ, obtained in the middle of the stable layer for all of our 18 Godunovtype schemes. In particular, the frequency power spectra shown in Fig. 14 are Fourier projections of an array^{16} containing the value of υ at (x, y, z) = (−0.2, 2.5, −0.2) at each time step. There is a clear distinction between the frequency range corresponding to the regime of IGWs (ω < N_{BV}), where all spectra are almost flat, and the highfrequency range (ω > N_{BV}, where only sound waves are nonevanescent) in which the spectra are characterized by a steep decrease in power (although non monotonic) toward higher frequencies. The forest of lines visible at frequencies in the range from 1 to 50 correspond to the resonant pressure modes of the cavity as predicted by the linear wave theory in the Cowling approximation (Aerts 2021, see also Fig. 15). In the frequency regime of IGWs, the power spectrum seems to converge to the reference solution when increasing the grid resolution or when using less dissipative Riemann solvers. A more quantitative analysis of the IGW spectrum would require substantially longer simulations, which would benefit from more time averaging as well as reach higher frequency resolution.
On the other hand, the highfrequency domain is well resolved even in simulations run on the 128^{3} grid, and differences between the spectra are clearly visible. Overall, there is a large spread in power (by almost 20 orders of magnitude at the Nyquist frequency) among the different methods tested in our study, which does not decrease with grid resolution. Therefore, numerical convergence is not achieved at high frequencies. We note that the spectra with the highest power density in the highfrequency range are always those obtained with limited reconstruction methods, in particular PPM84 and PPM08, while unlimited reconstruction methods tend to generate much quieter spectra and that are closer to the reference solution. In the simulations of the Kelvin–Helmholtz instability shown in in Sec. 3.1, we find that a power excess in highfrequency, shortwavelength sound waves is generated by the complex limiting procedure performed in PPM84 and PPM08 to reconstruct the Riemann states at cell interfaces. In this setup, a similar phenomenon is likely responsible for generating an acoustic power excess close to the Nyquist frequency, where the spread between the spectra is maximum. The power excess, however, is still large down to frequencies similar to the that of the fundamental oscillation mode of the cavity ω_{0} = 1.1. Also, PPMbased methods used in combination with HLLC excite much stronger resonant lines than those generated with LHLLC, and their peak power density is almost as high as the flat part of the spectrum associated with IGWs. When used with PPM84, both RUSANOV and LHLLC produce a broad feature in the frequency spectrum whose peak is at ω_{peak} ≈ 50 and 100 on the 128^{3} and 256^{3} grids, respectively. If these were sound waves, their wavelength λ would be (102)
which is close to the Nyquist frequency in space. Therefore, such a feature in the spectrum may be caused by oddeven cell decoupling traveling at the local speed of sound. Additional peaks are observed at even higher frequencies than the frequency of the least resolvable sound wave on the grid (with a wavelength of two cells), so these are most likely numerical artifacts. Tests performed with lower CFL factors of 0.4 and 0.2 show that the amplitude of these peaks is slightly reduced when using shorter time steps in the simulation, but the level of the continuum in the power spectrum remains essentially unaltered. In light of these results, we advise against using PPMbased methods in simulations of sound generation by lowMachnumber turbulence.
Fig. 13 Vertical profiles of the BruntVäisälä frequency (N_{BV}) averaged over the time interval t ∈ (30τ_{conv}, 32τ_{conv}) in the problem of turbulent convection and wave excitation. Each panel shows the results of simulations run using the same reconstruction scheme but different Riemann solvers. Here, we only show results from the 256^{3} grid to avoid cluttering the figures. The reference run (REF) is represented by a black dashed line. The black dashdotted line is the profile of the BruntVäisälä frequency at t = 0, and the vertical black dotted lines represent the position of the convective boundary at the beginning and the end of the chosen averaging time interval in the reference run, y = 2.168 and y = 2.174, respectively. The percentages shown in the insets for simulations run with the HLLC and LHLLC solvers represent the amplitude of the largest undershoot in N_{Bv} relative to N_{BV,t=0}, in the spatial range y ∈ (2.15, 2.40). The amplitude of the largest undershoot in the reference run is 12%. 
Fig. 14 Frequency power spectrum of the vertical velocity component υ extracted in the middle of the stable layer at (x, y, z) = (−0.2, 2.5, −0.2) over the time series t ∈ (10τ_{conv}, 32τ_{conv}) in the problem of turbulent convection and wave excitation. The left and right panels show results obtained on the 128^{3} and the 256^{3} grid, respectively. Each row of panels shows the results of simulations run using 6 reconstruction schemes with the same Riemann solver. The reference solution (indicated with REF) is the black curve. The convective turnover frequency (ω_{conv} = 2π/τ_{conv}), the BruntVäisälä frequency at (x, y, z) = (−0.2, 2.5, −0.2) (N_{BV}), and the frequency of the fundamental oscillation mode of the cavity (ω_{0} = 1.1) are represented by the black dasheddotted, dashed, and dotted lines, respectively. 
4 Performance metrics
After proving that all of the methods tested in our study converge to the correct solution for most of the physical quantities of interest^{17}, we can now search for the most efficient way to generate a numerical solution at the desired accuracy. In principle, one could use the L_{1} errors computed from the simulations of the Kelvin–Helmholtz instability (see Sec. 3.1) to find the resolution of the grid on which the scheme achieves a given level of accuracy. Then, the computational cost of the simulation run on such a grid can be estimated if the wallclock time spent by the program to perform a single cell update is known. Here, we prefer to use the results from the test problem involving turbulent convective flows and wave excitation (see Sec. 3.2), which is more challenging than the Kelvin–Helmholtz instability test and much closer to a real application of stellar hydrodynamics. However, due to the chaotic nature of the turbulent flows that develop in the convection zone, convergence in the L_{1} error norm cannot be achieved. Therefore, rather than measuring the computational cost per fixed accuracy in the sense of L_{1} errors, we opt to measure the computational cost of a simulation run with a specific combination of numerical methods that achieves a given level of effective resolution of the turbulent flows, N_{eff}. We estimate the effective resolution obtained on a given grid with cells as (103)
where λ_{vis}(N_{x}) is a characteristic dissipation length scale in the turbulent kinetic energy spectrum. In highresolution schemes, the amount of kinematic viscosity introduced into the system is not fixed, but rather its value depends on the local dynamical properties of the flow and it is often found to be a steep function of the spatial wavenumber (see, e.g., Porter & Woodward 1994). Therefore, λ_{vis}(N_{x}) cannot be defined uniquely. Here, we choose (104)
as a representative value for the dissipation length scale, with k_{h}_{,10}(N_{x}) being the spatial wavenumber at which the kinetic power spectrum rescaled by the Kolmogorov law drops by one dex from its maximum. In this approximation, the effective resolution is simply given by (105)
To perform a more precise measurement of N_{eff}(N_{x}), both the array of wavenumbers and of the kinetic energy spectrum are linearly interpolated on a finer grid. The values of the effective resolution (rescaled by the grid resolution N_{x}) obtained in the simulations included in our study are collected in Table 2. There is a clear trend toward higher effective resolution when using progressively higherorder reconstruction methods or less dissipative Riemann solvers.
By combining the values of N_{eff}(N_{x}) and the average wallclock time δt(N_{x}) spent by the program to perform a single cell update, shown in Table 3, it is possible to estimate the computational effort required by each combination of numerical options to achieve the effective resolution of PPM08 + LHLLC^{18} as (106)
Equation (106) is rescaled such that the cost of the simulation run with PPM08 + LHLLC on any given grid is unity. All values of Θ(N_{x} = 128) and Θ(N_{x} = 256) are shown in Table 4. The spread of almost four orders of magnitude in the relative computational cost among the different methods is due to the steep dependence of Θ(N_{x}) on N_{eff}(N_{x}), which varies by as much as a factor of 10 (see Table 2). On the other hand, Θ(N_{x}) only scales linearly with δt(N_{x}), which in our simulations varies at most by a factor of ≈3 (with LIN + RUSANOV and PPM84 + LHLLC achieving the lowest and highest wallclock time per cellupdate, respectively). The most expensive combination of methods among those tested in our study (in terms of computational cost per fixed resolving power) is LIN + VL + RUSANOV, ~ 1000 times as expensive as PPM08 + LHLLC. The strong numerical dissipation generated by RUSANOV leads to very poor performance of the finitevolume scheme even when used in combination with PPMbased reconstruction methods. Only with the unlimited PSH method such a flux function is capable of achieving acceptable performance (Θ(N_{x}) ≈ 4). When the reconstruction scheme is kept the same, using progressively less dissipative Riemann solvers decreases the cost of reaching the same effective resolution. At the typical Mach numbers encountered in the convection zone ( ≈ 0.02), the combination of HLLC and secondorder reconstruction schemes or PAR is considerably more expensive than PPM08 + LHLLC (Θ(N_{x}) in the range from 10 to 30). The performance of the scheme increases when HLLC is used with PPMbased reconstruction schemes (Θ(N_{x}) in the range from 2 to 6) and it is even higher than that of PPM08 + LHLLC when coupled to PSH. PSH + LHLLC is the most performant method on both grids according to the chosen metric. Overall, the computational cost of the finitevolume scheme is considerably reduced when using the lowdissipation solver (by a factor from 2 to 10 with respect to a scheme using HLLC and the same spatial reconstruction method). The relative performance of the scheme only varies slightly when LHLLC is used in combination with reconstruction methods less accurate than PSH, with the worst performance being achieved by LIN + VL (Θ(N_{x}) ≈ 3).
We note that, for most methods, the value of Θ(N_{x}) sensibly increases (by as much as a factor of ≈2) from the 128^{3} to the 256^{3} grid. Such differences are due to the small but systematic decrease of the rescaled effective resolution N_{eff} (N_{x})/N_{x} with increasing the grid resolution (this behavior is also confirmed by the reference solution computed on the 512^{3} grid, see Table 2). Although the differences observed in N_{eff}(N_{x})/N_{x} among the two grids are at most 20%, they are significantly amplified after applying the steep scaling relation between Θ(N_{x}) and N_{eff}(N_{x}). However, this effect is negligible if compared to the large spread obtained in the values of Θ(N_{x}) on a given grid. Furthermore, we stress that our measure of Θ(N_{x}) is based on a crude approximation of effective resolution of the turbulent flows, so the values provided in Table 4 should only be taken as estimates of the relative performance of different Godunovtype methods in simulations of turbulent convection. The measurements of the absolute performance provided in Table 3, which are needed to estimate Θ(N_{x}), may also depend on the parallelization strategy and the number of cores used to run the simulations, especially for the methods that require many ghost cells (e.g., PPM08 and PSH) and are therefore characterized by higher communication costs.
Fig. 15 Frequency power spectrum of ν as in Fig. 14 but showing only the results of PPM08 + HLLC on the 256^{3} grid. The vertical black dotted lines represent the resonant modes of the cavity which we derived from the theory of linear oscillations in Cowling approximation (see, e.g., Aerts 2021). The pink line at ω = 1.1 is the frequency of the fundamental oscillation mode of the cavity. 
Mean wallclock time spent by the program to advance the solution by one time step with RK3 in a single cell of the computational grid, in units of µs.
5 Summary and conclusions
Highresolution, finitevolume schemes are popular methods for simulating the behaviour of astrophysical fluids. There is a wide range of spatialreconstruction schemes, numerical flux functions and timeintegration methods that can be combined into a numerically stable and robust, finitevolume scheme. Focusing on subsonic flows, we have analyzed the accuracy and computational cost of all possible combinations of six spatialreconstruction schemes and three numerical flux functions, meaning 18 methods in total. The numerical solution was marched in time with a semidiscrete scheme based on a thirdorder Runge–Kutta method. This choice was motivated by our focus on the spatial accuracy of the schemes and the related need to suppress timestepping errors.
We considered two main test problems. The first is a Kelvin–Helmholtz instability problem, in which the initial shear flows have Mach numbers of 10^{−1}, 10^{−2}, and 10^{−3}. We used a smooth initial condition to make it possible to obtain numerically convergeable solutions to the inviscid Euler equations at a fixed point in time in the nonlinear phase of the instability. We characterized the accuracy of the solutions by (i) performing a qualitative assessment of the spatial structure of the solutions and (ii) by measuring L_{1} errors with respect to a reference solution obtained on a fine grid of 8192 × 4096 cells.
The lowMach flux function LHLLC is found to be much less dissipative and much more accurate than the HLLC and RUSANOV flux functions at the two lowest Mach numbers considered. Even though the overall 2D numerical scheme is 2ndorder accurate, the errors produced by different spatial reconstruction schemes span as much as an order of magnitude. Unlimited reconstruction schemes of increasing order, up to 7th, provide progressively more accurate results when the solution is sufficiently well resolved. However, these schemes are of limited use for the advection of mass fractions due to the oscillations and overshoots they produce around discontinuities. This issue is largely eliminated by the use of limiters in the schemes LIN+VL, PPM84, and PPM08. However, we find that the switching behaviour of the limiters introduces spurious structures and smallscale sound waves to the solutions, which severely reduce the accuracy of the methods when applied to slow flows. Our hybrid method PSH, which combines unlimited, 7thorder reconstruction for dynamic variables with PPM08 for mass fractions, is found to provide the most accurate solutions in nearly all of our simulations of the Kelvin–Helmholtz instability.
Our second test problem, chosen to be as close as possible to practical applications in the dynamics of stellar interiors, involves stratified turbulent convection generating waves in and entraining mass from an overlying stably stratified layer. The typical Mach number of the convection is ≈0.02. The turbulent nature of the convective flow makes the solutions chaotic but space and timeaveraged quantities can still be meaningfully compared between different simulations. We ran the simulations on grids of 128^{3} and 256^{3} cells with one additional simulation on a 512^{3} grid serving as a reference solution.
Qualitatively speaking, the solutions match the trends seen in the simulations of the Kelvin–Helmholtz instability, that is to say combinations of the least dissipative flux functions with the highestorder reconstruction schemes provide the highest resolving power in both the convective and stable layers. Spatial spectra of kinetic energy computed in the convective and stable layers converge toward the reference solution with any of the 18 methods but the least dissipative combinations (especially PSH+LHLLC) are much closer to the reference solution than the most dissipative ones (especially LIN+VL+RUSANOV). The spatial spectrum in both layers obtained with PSH+LHLLC on a 256^{3} grid is essentially the same as that with PAR+LHLLC on a 512^{3} grid.
The growth rate of the convective layer converges upon grid refinement with all 6 reconstruction schemes and both the HLLC and LHLLC flux function but the convergence is significantly faster with the lowMach flux function LHLLC. The RUSANOV flux function produces extremely viscous flows and large deviations from the reference solution in terms of the position of the upper convective boundary.
We show that schemes that do not use limiting for dynamic variables (LIN, PAR, PSH), unlike those that do (LIN+VL, PPM84, PPM08), produce overshoots in the BruntVäisälä frequency at the relatively sharp convective boundary. The magnitude of the overshoots increases with increasing order of accuracy of the method, as expected. Neither the growth rate of the convective layer, nor the temporal and spatial spectra of internal gravity waves are affected by this phenomenon in our simulations but we recommend careful testing when such highorder methods are applied to other setups or used to derive quantities not investigated here.
Temporal spectra of the vertical component of velocity in the stable layer show that all six reconstruction schemes produce qualitatively similar, flat and featureless spectra in the regime of internal gravity waves. On the other hand, the spectra of pressure (i.e., high frequency) waves are sensitive to the choice of the reconstruction scheme even with the lowMach flux function LHLLC. The spectra also reinforce our conclusions based on the Kelvin–Helmholtz problem that methods that apply limiters to dynamic variables generate much more acoustic noise than methods that do not. This includes the popular methods PPM84 and PPM08. Although we do not know what the acoustic spectrum in our problem should look like, we show that the frequencies of vertical resonant modes agree with 1D linear wave theory.
Finally, we used the steep dropoff of the spatial kinetic energy spectra in the convective layer close to the grid scale to define a measure of effective resolution. Rescaling the wallclock time of our 128^{3} and 256^{3} simulations, we estimated the computational cost of matching the same effective resolution using our 18 methods. The cost is found to span almost four orders of magnitude. The largest cost reduction comes from choosing the least dissipative flux function, which is LHLLC in our study. We show that the reconstruction schemes PSH, PAR, and PPM08, in order of increasing cost, are the cheapest at the same effective resolution.
Our study demonstrates that it is generally advisable to use lowdissipation Riemann solvers in Godunovtype finitevolume schemes for simulating lowMachnumber flows. The lowMach fix φ in Eq. (76) can easily be implemented in any HLLCtype scheme without affecting the absolute performance of the code, and it reduces the amount of computation required to achieve the same accuracy as HLLC by a factor ranging from 2 to 10 at typical Mach numbers of ≈ 0.01. At Mach numbers of ≈10^{−4} or ≈10^{−3}, similar to those encountered in the convective cores of main sequence stars, the usage of a lowMach solver decreases the computational cost per fixed accuracy with respect to a conventional Riemann solver by even larger factors (Leidi et al. 2022). Although lowdissipation Riemann solvers such as LHLLC bring clear advantages in simulations of subsonic flows when used in Eulerian codes, their robustness and accuracy properties in quasiLagrangian, movingmesh schemes still has to be verified. We note that such schemes hold promise for a particularly significant improvement because mesh motions aim at reducing advection errors by minimizing the velocities relative to the cell interfaces. This naturally leads to lowMachnumber flows in the comoving frame, which are better modeled by lowMachnumber Riemann solvers. We are currently testing an implementation of the LHLLC/D solvers of Minoshima & Miyoshi (2021) in the movingmesh MHD code AREPO (Springel 2010), which we will report in a forthcoming study.
On the other hand, the choice of the spatial reconstruction scheme depends on the physical application of interest. Overall, secondorder reconstruction methods lead to the generation of considerably more inaccurate results than higherorder schemes when the Mach number of the flow is low. If numerical overshoots have to be suppressed whilst retaining high resolving power, PPM84 and PPM08 should be the methods of choice. However, higherorder unlimited reconstruction methods should be preferred over PPMbased methods in simulations of sound wave generation by subsonic turbulence because they generate much less numerical noise in the acoustic spectrum. Less oscillatory alternatives to the hybrid PSH method, such as the 3thorder piecewise parabolic hybrid (PPH) or the 5thorder piecewise quartic hybrid (PQH) methods described in Appendix A, could offer the best compromise between the complete removal of numerical overshoots and reduction of acoustic noise.
Finally, we note that there are complex astrophysical environments where both high and lowMachnumber flows can coexist. This scenario often occurs, for instance, in accretion disks (see, e.g., Klessen & Hennebelle 2010), star forming regions (see, e.g., Wang et al. 2024), or outer layers of stars (see, e.g., Wedemeyer et al. 2017). To model such diverse flow regimes, the simple use of unlimited higherorder reconstructions and lowdissipation solvers may give rise to numerical instabilities, especially in the proximity of shock fronts. To cure this problem, an alternative approach based on the local dynamical properties of the flow is needed (see, e.g., Mignone et al. 2011; Fleischmann et al. 2020). One possibility is to add a selection criterion to the function that computes the flux at a cell interface such that the LHLLC solver is used if the Mach number of the flow is lower than a conservative threshold (e.g., ≲ 0.6), whereas a proper shockcapturing method (e.g., HLL or HLLC) is chosen for modeling faster flows. A similar argument can be made for the choice of the spatial reconstruction scheme. To avoid introducing large oscillations in the state quantities near large, nonlinear discontinuities, schemes that are more dissipative and possibly close to being TVD (such as LIN+VL, PPM84, PPM08, or WENO) should be preferred over high order unlimited methods.
Acknowledgements
The work of GL and FKR is supported by the German Research Foundation (DFG) through the grant RO 3676/31. We acknowledge support by the Klaus Tschira Foundation. This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1 – 390900948 (the Heidelberg STRUCTURES Excellence Cluster). The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gausscentre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUCNG at Leibniz Supercomputing Centre (www.lrz.de). PVFE was supported by the U.S. Department of Energy through the Los Alamos National Laboratory (LANL). LANL is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy (Contract No. 89233218CNA000001). This work has been assigned a document release number LAUR2333257.
Appendix A The PPH and PQH methods
The idea of applying limiters to passive scalars (or mass fractions) only, which lead us to formulating the PSH method (Sect. 2.3.6), can be applied at any order of accuracy. In this section, we provide two methods that can be seen as lowerorder alternatives to PSH. Their potential advantages include a lower amplitude of overshoots in dynamic quantities and the need for fewer ghost cells.
We start with the unlimited parabolic method PAR (Sect. 2.3.3) and apply a twostep limiter to passive scalars. The first step is defined by the variable replacements (A.1)
These replacements are followed by the application of the PPM84 limiter defined by Eq. (31). We call the resulting method piecewise parabolic hybrid (PPH). The PPH method is exact wherever a(ζ) is locally parabolic, 3^{rd}order accurate for general but smooth functions a(ζ), and it requires two ghost cells at domain boundaries.
A more accurate method, which we call piecewise quartic hybrid (PQH), can be obtained by assuming that within cell i a(ζ) can be described by the quartic polynomial (A.2)
The five coefficients c_{n} are uniquely determined by the requirement that the averages of a(ζ) in cells i − 2 + n equal ā_{i}_{−2+n} for n = 0, 1, . . . , 4. The reconstructed states are then obtained by evaluating Eq. (A.2) at ζ_{i}_{−1/2} and ζ_{i}_{+1/2}, respectively. The resulting expressions are (A.3) (A.4)
The limiter we apply to passive scalars in the PQH method is the same as that in the PPH method, see above. The PQH method is exact wherever a(ζ) is locally a quartic polynomial, 5^{th}order accurate for general but smooth functions a(ζ), and it requires three ghost cells at domain boundaries.
We have only performed a small number of tests with these two methods. Specifically, we ran simulations of the KelvinHelmholtz problem (Sec. 3.1) with the initial Mach number of _{0} = 10^{−3} and the HLLC flux function on grids of 64 × 32 and 128 × 64 cells. There were no overshoots in the passive tracer in the results of these tests.
Appendix B 1D test cases
We use two simple 1D experiments – linear advection and the propagation of a linear sound wave – to compare the accuracy that the six reconstruction schemes can reach if not constrained by the 2^{nd}order accuracy limit imposed by our multidimensional scheme.
B.1 Linear advection
The initial conditions in this experiment correspond to a rightgoing contact wave: (B.1) (B.2) (B.3)
We consider the interval 0 ≤ x ≤ 1 with periodic boundary conditions and the equation of state of an ideal gas with the ratio of specific heats γ = 1.4. The average speed of sound is unity and it varies by 0.5% due to the density variation. The solution is sought at t = 10, when the sinusoid has been advected by one period and the analytic solution becomes identical to the initial condition. This makes the quantification of numerical errors trivial.
We use the same code for this experiment as we do for all the other experiments reported in this work. Since some of the schemes tested exceed 2^{nd} order of accuracy, cell averages cannot be approximated by sampling the initial condition at cell centers. We initialize the discrete density profile using the formula (B.4)
where Δx is the constant grid spacing and x_{i} the location of ith cell’s center. This formula is easily obtained by analytically averaging Eq. (B.1) in the interval x_{i} − Δx/2 ≤ x ≤ x_{i} + Δx/2. We suppress timestepping errors by using the RK3 scheme with CFL = 0.01. Because the Mach number of the flow is 0.1, this means that the wave is advected by only 10^{−3} Δx during every 3^{rd}orderaccurate time step. We use the LHLLC flux function for this experiment.
The L^{1} errors we obtain on grids with different numbers of cells N_{x} are shown in Fig. B.1. We test the LIN, LIN+VL, PAR, and PPM84 schemes on grids with 8 to 512 cells. We are forced to stop at grids of 64 cells for the most accurate schemes PPM08 and PSH, because the errors rapidly become dominated by the finite precision of floatingpoint arithmetic. The asymptotic orders of accuracy based on L^{1} errors measured on the three finest grids available are reported in Table 1. Both LIN and LIN+VL reach the expected 2^{nd} order of accuracy. The absolute errors produced by LIN+VL are somewhat larger as compared with LIN because of the presence of the van Leer limiter. The PPM84 scheme is formally 4^{th}order accurate but that only holds for monotonic solutions. PPM84 contains a limiter that flattens the slope at local extrema and that reduces the scheme’s asymptotic order to 2.3 in this test case. We have experimentally turned off both limiters present in PPM84, which increased the asymptotic order of accuracy to 4.0, as expected. The PAR, PPM08, and PSH methods reach the asymptotic orders of 3.0, 6.0, and 7.0, also matching theoretical expectations.
B.2 Sound wave
The initial conditions in this experiment correspond to a rightgoing sound wave: (B.5) (B.6) (B.7)
Fig. B.1 Relative L_{1} errors obtained by advecting a sinusoid by one period on 1D grids with N_{x} cells using different reconstruction schemes. A few scaling relations are shown to guide the eye. 
We use ρ_{0} = γ, p_{0} = 1 and the equation of state of an ideal gas with the ratio of specific heats γ = 1.4, so that the unperturbed speed of sound is c_{0} = 1. We set the amplitude Ψ_{0} = 10^{−10} to suppress nonlinear effects. This allows us to avoid accuracy constraints imposed by our 2^{nd}orderaccurate transformations between primitive and conserved variables. Roundoff errors are suppressed by using 128bit floating point numbers. In a way analogous to Eq. (B.4), we initialize the discrete wave such that the cell average of Ψ(x) in cell i is^{19} (B.9)
The simulations are stopped at t = 1, when the sound wave has propagated by a single wavelength around the periodic domain 0 ≤; x ≤ 1. Neglecting nonlinear effects, which are of order Ψ_{0}^{2}, = 10^{−20}, the evolved solution is expected to be identical to the initial condition. Timestepping errors are suppressed by using the RK3 scheme with CFL = 0.001, so the wave moves by only 10^{−3} Δx per time step. We use the LHLLC flux function for this experiment.
We show L^{1} errors in the velocity u on grids with 8 to 512 cells in Fig. B.2. The simulation series with the most accurate PPM08 and PSH schemes are stopped at grids of 64 cells, because the absolute errors rapidly approach 10^{−20}, that is the magnitude of the residual nonlinear effects. The asymptotic orders of accuracy, as defined by the L^{1} errors measured on the three finest grids available, are 2.0, 3.0, 6.0, and 7.0 for the LIN, PAR, PPM08, and PSH methods, respectively, matching theoretical expectations. The LIN+VL and PPM84 methods only reach the order of 1.6, which seems to be caused by the combination of the lowMach flux function LHLLC with the limiters contained in these two methods. We have checked that both methods reach the order of 2.0 with the RUSANOV and HLLC flux functions (PPM84 flattens the slope at local extrema, which reduces its order, see also Sect. B.1).
The L^{1} errors discussed so far do not distinguish between amplitude and phase errors. We quantify amplitude errors by computing the relative loss of the total kinetic energy from t = 0 to t = 1, see Fig. B.3. Some kinetic energy is lost in all of the simulations, which is a sign of stability. The amount of energy dissipated decreases with the 3^{rd} power of the grid spacing for LIN, LIN+VL, PAR, and PPM84. It is not immediately clear why this is the case given that the methods have different orders of accuracy, but we did not investigate this further. In case of PSH, the dissipation rate decreases with the 7^{th} power of the grid spacing, matching the method’s order of accuracy. PPM08 is a special case – the method preserves extrema and sinusoids resolved by 16 or more cells per wavelength turn out to be smooth enough not to trigger any of the method’s limiters. PPM08 reduces to a simple interpolation function in this special case, eliminating jumps at all cell interfaces, see Sect. 2.3.5. This, in turn, eliminates all explicit dissipative terms in the flux function. The dissipation rate drops by many orders of magnitude and timestepping errors start to dominate.
Fig. B.2 Relative L_{1} errors in velocity u obtained by propagating a sinusoidal sound wave by one period on 1D grids with N_{x} cells using different reconstruction schemes. A few scaling relations are shown to guide the eye. 
Fig. B.3 Relative decrease in kinetic energy after propagating a sinusoidal sound wave by one period on 1D grids with N_{x} cells using different reconstruction schemes. Two scaling relations are shown to guide the eye. 
Appendix C Sound wave generation with a compact third order WENO scheme
The results presented in Sects. 3.1 and 3.2 reveal that nonlinear, switching reconstruction schemes (such as LIN+VL, PPM84, and PPM08) generate excess acoustic power which is not present in simulations run with unlimited (linear) reconstruction methods. The question naturally arises whether nonlinear but nonswitching schemes (e.g., WENO) also generate artificial sound waves in simulations of lowMachnumber flows. In this Appendix we try to answer this question using the compact third order (CTO) WENO scheme of Kolb (2014). This scheme reconstructs the quantity a within the cell i using a weighted sum of three polynomials, a_{L,R,C}(ζ). The ultimate expression for the underlying function a(ζ) reads (C.1)
The weighting coefficients are computed as (C.2)
m ∈ {L, C, R}, c_{L} = c_{R} = 0.25, c_{C} = 0.5, and IS_{m} are the smoothness indicators of the scheme^{20}, (C.4) (C.5) (C.6)
In Eq. (C.3), we set p = 2. The functions a(ζ)_{L,R} are onesided linear reconstructions (C.7) (C.8)
whereas a_{C}(ζ) is defined such that (C.9)
is the unique parabola that conserves the cell volume averages ā_{i}_{−1}, ā_{i}, and ā_{i}_{+1} over cells i − 1, i, and i + 1, respectively. This constraint implies (C.10)
The smoothness indicators IS_{L,R,C} are such that the smoothest among the polynomials a(ζ)_{L.R,C} has the largest weight in Eq. (C.1). This feature allows CTOWENO to achieve thirdorder spatial accuracy in smooth parts of the flow while at the same time it remains robust near discontinuities. The parameter ε _{WENO} that appears in the denominator of α_{m} avoids division by zero in the case a smoothness indicator IS_{m} = 0. Thus, the oscillatory behavior of a(ζ) is also determined by the value of this parameter. In fact, when ε _{WENO} ≫ IS_{L,R,C}, the weighting coefficients become almost equal and the scheme is unlimited even near discontinuities, making it similar to our PAR method. On the other hand, if ε _{WENO} ≪ IS_{m}, the weighting coefficients are only determined by IS_{L,R,C} and the scheme becomes close to being TVD near discontinuities or in poorly resolved parts of the flow.
Although CTOWENO, unlike LIN+VL or the PPM methods, does not involve any conditional statements, the weighting coefficients that occur in a(ζ) can still abruptly change from time step to time step if a barely resolved wave or feature in the flow crosses that particular cell. Such a rapid change in the form of the polynomial a(ζ) can generate highfrequency perturbations in the state quantities on the grid scale and affect the propagation of sound waves. Therefore, we expect a version of CTOWENO that uses small values of ε _{WENO} to generate artificial acoustic noise, whereas large values of ε _{WENO} should generate results closer to our unlimited reconstruction methods.
To test this hypothesis, we run a series of simulations of the setup involving turbulent convective flows and excitation of internal waves described in Sec. 3.2. We fix the grid resolution to 128^{3} cells and the Riemann solver is LHLLC. We run one simulation for each value of ε_{WENO} ∈ (10^{−12}, 10^{−10}, 10^{−8}, 10^{−6}, 10^{−4}, 10^{2}), so that several intermediate cases between the two extreme behaviors of the scheme (close to being TVD and fully oscillatory) are considered. We extract the frequency spectrum of the vertical velocity component from the middle of the stable layer as done for the analysis in Sec. 3.3. The results are shown in Fig. C.1. As expected, the power stored in sound waves is considerably increased when very small values of ε_{WENO} are used. For ε_{WENO} = 10^{−12}, the power spectrum obtained using CTOWENO resembles that produced by LIN+VL. Decreasing the value of ε_{WENO} progressively reduces (in a monotonic way) the energy of the sound waves and eventually the continuum of the power spectrum converges for ε_{WENO} ≳ 10^{−4}. In the study performed by Kolb (2014), the value of ε_{WENO} that achieves the optimal order of accuracy lies within the range^{21} Δx^{3} ≲ ε_{WENO} ≲ Δx^{2}, which, for this setup, corresponds to 10^{−6} ≲ ε_{WENO} ≲ 10^{−4}. However, we note that, in this test, the values of the primitive variables reconstructed at the grid cell interfaces span several orders of magnitude. Therefore, using a unique value of ε_{WENO} could potentially result in different oscillatory properties of the scheme depending on which variable is being reconstructed. One way to avoid this problem is to rescale the reconstructed quantity a such that its mean value across the stencil is close to unity. However, we decide not to investigate this effect on the generation of artificial sound waves further.
Fig. C.1 Frequency power spectrum of the vertical velocity component υ extracted in the middle of the stable layer at (x, y, z) = (−0.2, 2.5, −0.2) over the time series t ∈ (10τ_{conv}, 32τ_{conv}) in the problem of turbulent convection and wave excitation described in Sec. 3.2. These results are obtained on a grid with 128^{3} cells, using the LHLLC Riemann solver and the CTOWENO scheme of Kolb (2014) for different values of the parameter ε _{WENO} occurring in the smoothness indicators of the scheme (see Sect. C), ranging from 10^{−12} to 10^{2}. As a reference, power spectra obtained with the numerical options LHLLC + PPM84 and LHLLC + LIN + VL are also shown. The convective turnover frequency (ω_{conv} = 2π/τ_{conv}), the BruntVäisälä frequency at the location of the point probe (N_{BV}), and the frequency of the fundamental oscillation mode of the cavity (ω_{0} = 1.1) are represented by the black dasheddotted, dashed, and dotted lines, respectively. 
Appendix D Additional plots for the Kelvin–Helmholtz problem
References
 Aerts, C. 2021, Rev. Mod. Phys., 93, 015001 [Google Scholar]
 Anders, E. H., & Pedersen, M. G. 2023, Galaxies, 11 [Google Scholar]
 Andrassy, R., Herwig, F., Woodward, P., & Ritter, C. 2020, MNRAS, 491, 972 [NASA ADS] [CrossRef] [Google Scholar]
 Andrassy, R., Higl, J., Mao, H., et al. 2022, A&A, 659, A193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Andrassy, R., Leidi, G., Higl, J., et al. 2024, A&A, 683, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Balsara, D. S. 2017, Living Rev. Computat. Astrophys., 3, 2 [CrossRef] [Google Scholar]
 Beckwith, K., & Stone, J. M. 2011, ApJS, 193, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Berberich, J. P., Chandrashekar, P., & Klingenberg, C. 2021, Comput. Fluids, 219, 104858 [Google Scholar]
 Berlok, T., & Pfrommer, C. 2019, MNRAS, 485, 908 [NASA ADS] [CrossRef] [Google Scholar]
 Canivete Cuissa, J. R., & Teyssier, R. 2022, A&A, 664, A24 [Google Scholar]
 Chen, S., Lin, B., Li, Y., & Yan, C. 2020, SIAM J. Sci. Comput., 42, B921 [CrossRef] [Google Scholar]
 Colella, P. 1990, J. Computat. Phys., 87, 171 [NASA ADS] [CrossRef] [Google Scholar]
 Colella, P., & Woodward, P. R. 1984, J. Computat. Phys., 54, 174 [NASA ADS] [CrossRef] [Google Scholar]
 Colella, P., & Sekora, M. D. 2008, J. Computat. Phys., 227, 7069 [NASA ADS] [CrossRef] [Google Scholar]
 Courant, R., Friedrichs, K., & Lewy, H. 1928, Math. Ann., 100, 32 [Google Scholar]
 Day, M. S., & Bell, J. B. 2000, Combust. Theory Model., 4, 535 [NASA ADS] [CrossRef] [Google Scholar]
 Dobler, W., Haugen, N. E., Yousef, T. A., & Brandenburg, A. 2003, Phys. Rev. E, 68, 026304 [NASA ADS] [CrossRef] [Google Scholar]
 Dumbser, M., Castro, M., Parés, C., & Toro, E. F. 2009, Comput. Fluids, 38, 1731 [CrossRef] [MathSciNet] [Google Scholar]
 Dumbser, M., Balsara, D. S., Tavelli, M., & Fambri, F. 2019, Int. J. Numer. Methods Fluids, 89, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Edelmann, P. V. F. 2014, Dissertation, Technische Universität München [Google Scholar]
 Edelmann, P. V. F., & Röpke, F. K. 2016, in JUQUEEN Extreme Scaling Workshop 2016, eds. D. Brömmel, W. Frings, & B. J. N. Wylie, JSC Internal Report No. FZJJSCIB201601, 63 [Google Scholar]
 Edelmann, P. V. F., Röpke, F. K., Hirschi, R., Georgy, C., & Jones, S. 2017, A&A, 604, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Edelmann, P. V. F., Horst, L., Berberich, J. P., et al. 2021, A&A, 652, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fleischmann, N., Adami, S., & Adams, N. A. 2020, J. Comput. Phys., 423, 109762 [NASA ADS] [CrossRef] [Google Scholar]
 Flock, M., Dzyurkevich, N., Klahr, H., & Mignone, A. 2010, A&A, 516, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Godunov, S. K. 1959, Mat. Sbornik, 89, 271 [Google Scholar]
 Goffrey, T., Pratt, J., Viallet, M., et al. 2017, A&A, 600, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Greenough, J. A., & Rider, W. J. 2003, J. Computat. Phys., 196 [Google Scholar]
 Guillard, H., & Viozat, C. 1999, Comput. Fluids, 28, 63 [Google Scholar]
 Harten, A., Lax, P. D., & Leer, B. V. 1983, SIAM Rev., 25, 35 [Google Scholar]
 Harten, A., Engquist, B., Osher, S., & Chakravarthy, S. R. 1987, J. Computat. Phys., 71, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Horst, L., Edelmann, P. V. F., Andrássy, R., et al. 2020, A&A, 641, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Horst, L., Hirschi, R., Edelmann, P. V. F., Andrassy, R., & Roepke, F. K. 2021, A&A, 653, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jiang, G.S., & Shu, C.W. 1996, J. Computat. Phys., 126, 202 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J. 2021, A&A, 651, A66 [Google Scholar]
 Klein, R. 2009, Theoret. Computat. Fluid Dyn., 23, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17 [NASA ADS] [CrossRef] [Google Scholar]
 Klingenberg, C., Schmidt, W., & Waagan, K. 2007, J. Computat. Phys., 227, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Kolb, O. 2014, SIAM J. Numer. Analys., 52, 2335 [CrossRef] [Google Scholar]
 Kritsuk, A. G., Åke Nordlund, Collins, D., et al. 2011, ApJ, 737, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Küchler, C., Bewley, G., & Bodenschatz, E. 2019, J. Stat. Phys., 175, 617 [CrossRef] [Google Scholar]
 Latini, M., Schilling, O., & Don, W. S. 2007, J. Computat. Phys., 221, 805 [NASA ADS] [CrossRef] [Google Scholar]
 Lecoanet, D., McCourt, M., Quataert, E., et al. 2017, MNRAS, 455, 4274 [Google Scholar]
 Leidi, G., Birke, C., Andrassy, R., et al. 2022, A&A, 668, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leidi, G., Andrassy, R., Higl, J., Edelmann, P. V. F., & Röpke, F. K. 2023, A&A, 679, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LeVeque, R. J. 2002, Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Applied Mathematics (Cambridge University Press) [CrossRef] [Google Scholar]
 Lighthill, M. J. 1952, Proc. Roy. Soc. Lond. Ser. A, 211, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Liou, M.S. 2006, J. Computat. Phys., 214, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, X.D., Osher, S., & Chan, T. 1994, J. Computat. Phys., 115, 200 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars [Google Scholar]
 McNally, C. P., Lyra, W., & Passy, J.C. 2012, ApJS, 201, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448 [NASA ADS] [CrossRef] [Google Scholar]
 Miczek, F. 2013, Dissertation, Technische Universität München [Google Scholar]
 Miczek, F., Röpke, F. K., & Edelmann, P. V. F. 2015, A&A, 576, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mignone, A., Zanni, C., Tzeferacos, P., et al. 2011, ApJS, 198, 7 [Google Scholar]
 Minoshima, T., & Miyoshi, T. 2021, J. Computat. Phys., 446, 110639 [NASA ADS] [CrossRef] [Google Scholar]
 Miyoshi, T., & Kusano, K. 2005, J. Computat. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Motheau, E., Duarte, M., Almgren, A., & Bell, J. B. 2018, J. Computat. Phys., 372, 1027 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, B. 2020, Living Rev. Computat. Astrophys., 6, 3 [CrossRef] [Google Scholar]
 Musoke, G., Young, A. J., & Birkinshaw, M. 2020, MNRAS, 498, 3870 [NASA ADS] [CrossRef] [Google Scholar]
 Muthsam, H. J., Kupka, F., LöwBaselli, B., et al. 2010, New Astron., 15, 460 [Google Scholar]
 Porter, D. H., & Woodward, P. R. 1994, ApJS, 93, 309 [NASA ADS] [CrossRef] [Google Scholar]
 Radice, D., Couch, S. M., & Ott, C. D. 2015, Computat. Astrophys. Cosmol., 2, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Rieper, F. 2011, J. Computat. Phys., 230, 5263 [NASA ADS] [CrossRef] [Google Scholar]
 Rizzuti, F., Hirschi, R., Arnett, W. D., et al. 2023, MNRAS, 523, 2317 [NASA ADS] [CrossRef] [Google Scholar]
 Robertson, B. E., Kravtsov, A. V., Gnedin, N. Y., Abel, T., & Rudd, D. H. 2010, MNRAS, 401, 2463 [NASA ADS] [CrossRef] [Google Scholar]
 Roe, P. L. 1981, J. Computat. Phys., 43, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Rusanov, V. 1962, USSR Computat. Math. Math. Phys., 1, 304 [CrossRef] [Google Scholar]
 San, O., & Kara, K. 2015, Comput. Fluids, 117, 24 [CrossRef] [Google Scholar]
 Seo, J., & Ryu, D. 2023, ApJ, 953, 39 [CrossRef] [Google Scholar]
 Shu, C.W. 2009, SIAM Rev., 51, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, C.W., & Osher, S. 1988, J. Computat. Phys., 77, 439 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V. 2010, MNRAS, 401, 791 [Google Scholar]
 Sutherland, B. R. 2010, Internal Gravity Waves (Cambridge University Press) [CrossRef] [Google Scholar]
 Teissier, J.M., & Müller, W.C. 2023, JCP, submitted [arXiv:2306.09856] [Google Scholar]
 Thornber, B., Mosedale, A., Drikakis, D., Youngs, D., & Williams, R. 2008, J. Computat. Phys., 227, 4873 [NASA ADS] [CrossRef] [Google Scholar]
 Toro, E. F. 1991, Proc. Roy. Soc. Lond. Ser. A, 434, 683 [NASA ADS] [CrossRef] [Google Scholar]
 Toro, E. F. 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction (Berlin, Heidelberg: Springer) [CrossRef] [Google Scholar]
 Toro, E. F., Spruce, M., & Speares, W. 1994, Shock Waves, 4, 25 [Google Scholar]
 van Leer, B. 1974, J. Computat. Phys., 14, 361 [NASA ADS] [CrossRef] [Google Scholar]
 van Leer, B. 1977, J. Computat. Phys., 23, 276 [NASA ADS] [CrossRef] [Google Scholar]
 van Leer, B. 1979, J. Computat. Phys., 32, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Viallet, M., Baraffe, I., & Walder, R. 2011, A&A, 531, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wang, C., Wang, K., Xu, F.W., et al. 2024, A&A, 681, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wedemeyer, S., Kučinskas, A., Klevas, J., & Ludwig, H.G. 2017, A&A, 606, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wongwathanarat, A., GrimmStrele, H., & Müller, E. 2016, A&A, 595, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Woodward, P. R., Herwig, F., & Lin, P.H. 2014, ApJ, 798, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Xie, W., Zhang, R., Lai, J., & Li, H. 2019, Int. J. Numer. Methods Fluids, 89, 430 [NASA ADS] [CrossRef] [Google Scholar]
Highresolution schemes for gas dynamics can be fully discrete, where the system of equations is discretized both in space and time, or semidiscrete, where spatial discretization is performed first while leaving the problem continuous in time. In the latter approach, state quantities are then advanced in time using any standard numerical solver for systems of ordinary differential equations.
The influence of the parameter that occurs in the smoothness indicators of the compact third order WENO scheme of Kolb (2014) on the generation of sound waves by turbulent flows is discussed in Appendix C.
This formally holds for linear advection in one spatial dimension. In practice, we do not observe oscillations even when the LIN+VL method is applied to a multidimensional system of conservation laws, see Sec. 3.1.
Although SLH was already equipped with several lowMach solvers, such as AUSM^{+}up (Liou 2006) and MiczekRoe (Miczek et al. 2015), we decided not to use them in this study because they all suffer from a more restrictive stability criterion on the time step than LHLLC when used in combination with explicit time steppers, whilst all of these fluxes provide very similar results in terms of accuracy.
Other lowMach corrections for the HLLC Riemann solver can be found, for instance, in Thornber et al. (2008); Rieper (2011); Xie et al. (2019); Chen et al. (2020); Fleischmann et al. (2020).
Similar plots for all initial Mach numbers and all computational grids are available on Zenodo (https://zenodo.org/doi/10.5281/zenodo.10280900).
See an animation available on Zenodo (https://zenodo.org/doi/10.5281/zenodo.10280900).
The animations show high values of the velocity divergence developing early on where the second derivative of the transition function η(y) (Eq. (87)) is discontinuous. This higherorder discontinuity may influence the behaviour of limiters.
A prominent example of a nonconvergent numerical solution is the frequency power spectrum of sound waves in the test problem of turbulent convection and mass entrainment shown in Sec. 3.2.
For a formal derivation of these indicators, see Sect. 2.1 of Kolb (2014).
All Tables
Mean wallclock time spent by the program to advance the solution by one time step with RK3 in a single cell of the computational grid, in units of µs.
All Figures
Fig. 1 Spacetime diagram showing the wave structure of the HLLC solver for a Riemann problem of gas dynamics. The S_{l,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}. 

In the text 
Fig. 2 Reference solution to the Kelvin–Helmholtz problem with initial Mach number _{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 _{0}^{−1}), at an early stage of nonlinear evolution (t = 0.4 _{0}^{−1}), at a stage when the primary vortices have fully formed (t = 0.8 _{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 _{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. 

In the text 
Fig. 3 Resolution dependence of the minimum scale height min(H_{x}) 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 N_{x} of computational cells along the x axis. The initial Mach number is _{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. 

In the text 
Fig. 4 Distributions of the mass fraction X of the passive scalar in simulations of the Kelvin–Helmholtz instability with the initial Mach number _{0} = 10^{−2} on a 128 × 64 grid. The results are shown at t = 0.8 _{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. 

In the text 
Fig. 5 Relative L_{1} errors for different variables (rows), flux functions (columns), and reconstruction methods (legend) as functions of the grid resolution in the simulations of the Kelvin–Helmholtz instability with the initial Mach number _{0} = 10^{−2}. The dashed and dotted lines, which are at the same locations in all of the panels, show the 1st and 2ndorder scalings to guide the eye. 

In the text 
Fig. 6 Magnitude of velocity divergence in simulations of the Kelvin–Helmholtz instability with the initial Mach number _{0} = 10^{−2} on a 512 × 256 grid. The results are shown at t = 0.8 _{0}^{−1}. Rows and columns show different reconstruction schemes and numerical flux functions, respectively. The smallscale structure produced by the LIN+VL, PPM84, and PPM08 schemes, which does not occur with the LIN, PAR, and PSH schemes, is caused by the use of limiters. 

In the text 
Fig. 7 Hydrostatic stratification of density, gas pressure, pseudoentropy (p/ρ^{γ}), and gravity at t = 0 in the problem of turbulent convection and wave excitation at a convective boundary, based on the initial conditions provided in Andrassy et al. (2022). The distribution of the mass fractional abundance of a passive tracer (X) and the heat source () are also shown, here represented by a dashdotted pink and dashed black line, respectively. In this setup, the gravitational acceleration smoothly turns to zero at the bottom and top boundaries of the spatial domain to make the problem consistent with the usage of reflecting, stress free boundary conditions, which for the gas pressure imply ∂p/∂y = 0. 

In the text 
Fig. 8 Pseudoentropy (A = p/ρ^{γ}) fluctuations with respect to the hydrostatic, background state computed at different times during the development of the convection, as indicated by the insets, in the problem of turbulent convection and wave excitation. Data are visualized in the z = 0 plane. This simulation was run using the PAR method and LHLLC on a 512^{3} grid and we consider it our reference solution. Convective plumes start exciting IGWs at the upper convective boundary by t ≈ τ_{conv} (see center bottom panel). 

In the text 
Fig. 9 Snapshots obtained at the final simulation time (t = 32τ_{conv}) showing the distribution of the Mach number at z = 0 in the simulations of turbulent convection and wave propagation in a 3D box. The grid resolution is 256^{3}. Each panel is cut between y = 1.3 and y = 2.7 for a better visualization of the flows in the proximity of the upper convective boundary. 

In the text 
Fig. 10 Kinetic energy spectra computed on the horizontal plane y = 1.5 of the convective layer for all of our 18 combinations of Riemann solvers and spatial reconstruction schemes in the problem of turbulent convection and wave excitation. The spectra have been averaged over the time interval t ∈ (10τ_{conv}, 32τ_{conv}) and rescaled by the Kolmogorov law (k_{h}^{−5/3}) and by the value Ê_{0} of the spectral energy density of the reference run at k_{h} = 15. 

In the text 
Fig. 11 Kinetic energy spectra computed on the horizontal plane y = 2.5 of the stable layer for all of our 18 combinations of Riemann solvers and spatial reconstruction schemes in the problem of turbulent convection and wave excitation. The spectra have been averaged over the time interval t ∈ (10τ_{conv}, 32τ_{conv}). The black solid line represents the reference run (REF), which is computed using the LHLLC Riemann solver and the PAR reconstruction scheme on a 512^{3} grid. 

In the text 
Fig. 12 Time evolution of the position of the upper convective boundary (y_{cb}) in the problem of turbulent convection and entrainment described in Sec. 3.2 The top and bottom rows show results from 128^{3} and 256^{3} simulations, respectively. The black solid line represents the reference run (REF). 

In the text 
Fig. 13 Vertical profiles of the BruntVäisälä frequency (N_{BV}) averaged over the time interval t ∈ (30τ_{conv}, 32τ_{conv}) in the problem of turbulent convection and wave excitation. Each panel shows the results of simulations run using the same reconstruction scheme but different Riemann solvers. Here, we only show results from the 256^{3} grid to avoid cluttering the figures. The reference run (REF) is represented by a black dashed line. The black dashdotted line is the profile of the BruntVäisälä frequency at t = 0, and the vertical black dotted lines represent the position of the convective boundary at the beginning and the end of the chosen averaging time interval in the reference run, y = 2.168 and y = 2.174, respectively. The percentages shown in the insets for simulations run with the HLLC and LHLLC solvers represent the amplitude of the largest undershoot in N_{Bv} relative to N_{BV,t=0}, in the spatial range y ∈ (2.15, 2.40). The amplitude of the largest undershoot in the reference run is 12%. 

In the text 
Fig. 14 Frequency power spectrum of the vertical velocity component υ extracted in the middle of the stable layer at (x, y, z) = (−0.2, 2.5, −0.2) over the time series t ∈ (10τ_{conv}, 32τ_{conv}) in the problem of turbulent convection and wave excitation. The left and right panels show results obtained on the 128^{3} and the 256^{3} grid, respectively. Each row of panels shows the results of simulations run using 6 reconstruction schemes with the same Riemann solver. The reference solution (indicated with REF) is the black curve. The convective turnover frequency (ω_{conv} = 2π/τ_{conv}), the BruntVäisälä frequency at (x, y, z) = (−0.2, 2.5, −0.2) (N_{BV}), and the frequency of the fundamental oscillation mode of the cavity (ω_{0} = 1.1) are represented by the black dasheddotted, dashed, and dotted lines, respectively. 

In the text 
Fig. 15 Frequency power spectrum of ν as in Fig. 14 but showing only the results of PPM08 + HLLC on the 256^{3} grid. The vertical black dotted lines represent the resonant modes of the cavity which we derived from the theory of linear oscillations in Cowling approximation (see, e.g., Aerts 2021). The pink line at ω = 1.1 is the frequency of the fundamental oscillation mode of the cavity. 

In the text 
Fig. B.1 Relative L_{1} errors obtained by advecting a sinusoid by one period on 1D grids with N_{x} cells using different reconstruction schemes. A few scaling relations are shown to guide the eye. 

In the text 
Fig. B.2 Relative L_{1} errors in velocity u obtained by propagating a sinusoidal sound wave by one period on 1D grids with N_{x} cells using different reconstruction schemes. A few scaling relations are shown to guide the eye. 

In the text 
Fig. B.3 Relative decrease in kinetic energy after propagating a sinusoidal sound wave by one period on 1D grids with N_{x} cells using different reconstruction schemes. Two scaling relations are shown to guide the eye. 

In the text 
Fig. C.1 Frequency power spectrum of the vertical velocity component υ extracted in the middle of the stable layer at (x, y, z) = (−0.2, 2.5, −0.2) over the time series t ∈ (10τ_{conv}, 32τ_{conv}) in the problem of turbulent convection and wave excitation described in Sec. 3.2. These results are obtained on a grid with 128^{3} cells, using the LHLLC Riemann solver and the CTOWENO scheme of Kolb (2014) for different values of the parameter ε _{WENO} occurring in the smoothness indicators of the scheme (see Sect. C), ranging from 10^{−12} to 10^{2}. As a reference, power spectra obtained with the numerical options LHLLC + PPM84 and LHLLC + LIN + VL are also shown. The convective turnover frequency (ω_{conv} = 2π/τ_{conv}), the BruntVäisälä frequency at the location of the point probe (N_{BV}), and the frequency of the fundamental oscillation mode of the cavity (ω_{0} = 1.1) are represented by the black dasheddotted, dashed, and dotted lines, respectively. 

In the text 
Fig. D.1 As Fig. 4 but with the initial Mach number _{0} = 10 ^{1}. 

In the text 
Fig. D.2 As Fig. 4 but with the initial Mach number _{0} = 10^{−3}. 

In the text 
Fig. D.3 As Fig. 5 but with the initial Mach number _{0} = 10^{−}^{1}. 

In the text 
Fig. D.4 As Fig. 5 but with the initial Mach number _{0} = 10^{−3}. 

In the text 
Fig. D.5 As Fig. 6 but with the initial Mach number _{0} = 10^{−}^{1}. 

In the text 
Fig. D.6 As Fig. 6 but with the initial Mach number _{0} = 10^{−3}. 

In the text 
Fig. D.7 As Fig. 4 but showing simulations computed on the 4096 × 2048 grid. 

In the text 
Fig. D.8 As Fig. 6 but showing simulations computed on the 4096 × 2048 grid. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.