Performance of high-order Godunov-type methods in simulations of astrophysical low Mach number flows

TL;DR

This study systematically benchmarks 18 combinations of spatial reconstruction schemes and approximate Riemann solvers within the fully compressible SLH code for subsonic astrophysical flows ($10^{-3} \lesssim \mathcal{M} \lesssim 10^{-1}$), focusing on a 2D Kelvin–Helmholtz instability and a 3D turbulent convection with internal gravity waves. It demonstrates that a low-dissipation Riemann solver (LHLLC) paired with a high-order reconstruction (notably PSH) yields the best balance between accuracy and computational cost, achieving near-reference results at lower grid resolution than less dissipative pairings. The work also reveals that unlimited high-order reconstructions improve spectral resolution and reduce dissipation but can generate artificial acoustic noise or boundary undershoots unless carefully managed with limiters or hybrid approaches. Overall, the results provide practical guidance for modeling low-Mach astrophysical flows, suggesting that LHLLC-based schemes with hybrid reconstructions offer substantial efficiency gains, while recognizing trade-offs in oscillations, noise, and boundary behavior. The findings have implications for simulations of stellar interiors, accretion flows, and geophysical-like convection in astrophysical contexts, where both accuracy and efficiency are crucial.

Abstract

High-order 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 Godunov-type methods, these two steps can be treated independently, so that many different schemes can in principle be built from the same numerical framework. In this work, we use our fully compressible Seven-League Hydro (SLH) code to test the accuracy of six reconstruction methods and three approximate Riemann solvers on two- and three-dimensional (2D and 3D) problems involving subsonic flows only. We consider Mach numbers in the range from $10^{-3}$ to $10^{-1}$ in a well-posed, 2D, Kelvin--Helmholtz instability problem and a 3D turbulent convection zone that excites internal gravity waves in an overlying stable layer. 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 "low-dissipation" Riemann solver and a sextic reconstruction scheme, (ii) the low-dissipation 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 slope-limiting techniques tend to generate artificial, high-frequency acoustic waves during the evolution of the flow, (v) unlimited reconstruction methods introduce oscillations in the thermal stratification near the convective boundary, where the entropy gradient is steep.

Figures (27)

• Figure 1: Space-time diagram showing the wave structure of the HLLC solver for a Riemann problem of gas dynamics. $S_{\mathrm{L},\mathrm{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, $\bm{U}^*_{\mathrm{L},\mathrm{R}}$.
• Figure 2: Reference solution to the Kelvin--Helmholtz problem with initial Mach number $\mathcal{M}_0 = 10^{-2}$. The solution was computed using PSH reconstruction and the LHLLC flux function on an $8192 \,{\times}\, 4096$ grid. The mass fraction $X$ of the passive scalar is shown at four points in time: late in the linear growth of the instability ($t = 0.2\,\mathcal{M}_0^{\,-1}$), at an early stage of non-linear evolution ($t = 0.4\,\mathcal{M}_0^{\,-1}$), at a stage when the primary vortices have fully formed ($t = 0.8\,\mathcal{M}_0^{\,-1}$; the final time for all of our other Kelvin--Helmholtz simulations), and at a late stage when fine threads have formed inside the primary vortices ($t = 1.6\,\mathcal{M}_0^{\,-1}$). We use the same colour scale as in Figs. \ref{['fig:kh2d-machx-1.000e-02-128x64-ps']}, \ref{['fig:kh2d-machx-1.000e-03-128x64-ps']}, and \ref{['fig:kh2d-machx-1.000e-01-128x64-ps']}, although $0 \leq X \leq 1$ in the reference solution.
• Figure 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 $\mathcal{M}_0 = 10^{-2}$ and we use PSH reconstruction and the LHLLC flux function in this series of simulations. Once the steepest gradients in $X$ become resolved the minimum scale height starts to follow the linear scaling relations shown.
• Figure 4: Distributions of the mass fraction $X$ of the passive scalar in simulations of the Kelvin--Helmholtz instability with the initial Mach number $\mathcal{M}_0 = 10^{-2}$ on a $128 \times 64$ grid. The results are shown at $t = 0.8\,\mathcal{M}_0^{\,-1}$. Rows and columns show different reconstruction schemes and numerical flux functions, respectively. The colour 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.
• Figure 5: Relative $L_1$ errors for different variables (rows), flux functions (columns), and reconstruction methods (legend) as functions of grid resolution in the simulations of the Kelvin--Helmholtz instability with the initial Mach number $\mathcal{M}_0 = 10^{-2}$. The dashed and dotted lines, which are at the same locations in all of the panels, show the $1$st- and $2$nd-order scalings to guide the eye.