A low-dissipation HLLD approximate Riemann solver for a very wide range of Mach numbers

We propose a new Harten-Lax-van Leer discontinuities (HLLD) approximate Riemann solver to improve the stability of shocks and the accuracy of low-speed flows in multidimensional magnetohydrodynamic (MHD) simulations. Stringent benchmark tests verify that the new solver is more robust against a numerical shock instability and is more accurate for low-speed, nearly incompressible flows than the original solver, whereas additional computational costs are quite low. The novel ability of the new solver enables us to tackle MHD systems, including both high and low Mach number flows.


Introduction
A magnetohydrodynamic (MHD) simulation is an indispensable tool for studying the macroscopic dynamics of laboratory, space, and astrophysical plasmas.For compressible MHD simulations, shock-capturing schemes have been developed based on the solution to the Riemann problem in onedimensional hyperbolic conservation laws, which allows us to tackle a situation including supersonic flows.In particular, the Harten-Lax-van Leer discontinuities (HLLD) approximate Riemann solver developed by Miyoshi and Kusano (2005) is extensively implemented in modern MHD simulation codes by virtue of its robustness and accuracy (e.g., Fromang et al., 2006;Mignone et al., 2007;Stone et al., 2008;Lee and Deane, 2009;Zenitani and Miyoshi, 2011;Matsumoto et al., 2019).
In practical multidimensional MHD simulations, however, familiar shockcapturing schemes may suffer from numerical difficulties, which include a numerical shock instability for high Mach number flows and a degradation of the solution accuracy for low Mach number flows.We proposed a multistate lowdissipation advection upstream splitting method (MLAU; Minoshima et al., 2020) to remedy these difficulties, which is an MHD extension of all-speed advection upstream splitting methods for hydrodynamics (Liou, 2006;Shima and Kitamura, 2011;Kitamura and Shima, 2013).The MLAU scheme provides highly robust and accurate solutions of stringent problems, including high and low Mach number flows, whereas familiar shock-capturing schemes fail to resolve them, and it preserves the MHD discontinuities comparable to the HLLD scheme.This novel ability allows reliable simulations of wideranging Mach number flows in magnetized plasma.For users of the HLLD scheme to easily enjoy this ability, we propose a new HLLD approximate Riemann solver that implements the techniques used in the MLAU scheme.

Low-dissipation HLLD approximate Riemann solver
We consider one-dimensional MHD equations written in the following conservative form: where U and F are the state vector of conservative variables and the corresponding flux vector, respectively, and ρ, u = (u, v, w), B = (B x , B y , B z ), and e are the mass density, velocity, magnetic field, and total energy density.
The gas pressure P is determined from the equation of state for the ideal gas, where γ is the specific heat ratio.The solenoidal condition of the magnetic field gives B x = constant in one dimension.Equation ( 1) is discretized on computational cells into a finite volume form as follows: where F i±1/2 is the numerical flux at the interfaces of a cell ].The quality of the numerical solutions largely relies on an evaluation of the numerical flux.
The HLLD scheme solves the MHD Riemann problem at a cell interface for the left-and right-side variables U L,R as an initial state by allowing five eigenmodes in the Riemann fan.To obtain four intermediate states in the Riemann fan, the HLLD scheme consistently assumes that the normal velocity S M and the total (gas + magnetic) pressure P * t are constant over the Riemann fan, where ) are the robust estimations of the minimum and maximum signal speeds taken from Equation (67) in Miyoshi and Kusano (2005), and c f,max = max(c f,L , c f,R ) is the maximum fast magnetosonic wave speed, Given S M and P * t , the intermediate states and the corresponding numerical fluxes are algebraically calculated from the jump conditions across the five waves (see Miyoshi and Kusano (2005) for details).In the following subsections, we modify Equations ( 6) and (7) using the techniques applied in the MLAU scheme.

Shock detection
Shock-capturing schemes that can preserve the contact discontinuity, such as the Roe (Brio and Wu, 1988) and HLLD schemes, tend to suffer from a numerical shock instability when a multidimensional shock is well aligned to the grid spacing, and lead to catastrophic solutions such as odd-even decoupling and the Carbuncle phenomena (Quirk, 1994).Liou (2000) and Kim et al. (2003) argued that the pressure difference term in the interface mass flux (numerical flux of the density) is a possible cause of the instability when the flux is nearly parallel to the shock surface.The MLAU scheme then introduces a shock-detecting factor θ to eliminate the pressure difference term only at regions dangerous to the instability.This technique is implemented into the HLLD scheme through the normal velocity as follows: where a > 0 is a free parameter to make the factor θ (multiplied by the third term in Equation ( 9)) quickly approaching zero when the shock surface is nearly parallel to the z − x or x − y plane.Our numerical tests suggest that a = 1 is insufficient to suppress the instability and a ≥ 2 give reasonable solutions.We empirically use a = 4.

Pressure correction
Familiar one-dimensional shock-capturing schemes (Roe and HLLD) build the numerical flux of the normal momentum to include the velocity differ-ence term with a scale of fast magnetosonic speed, ∼ −ρc f (u R − u L )/2, which acts as a diffusion to numerically stabilize the compressible flows.The term is inappropriate for low-speed flows in multiple dimensions because a finite velocity difference in one direction does not necessarily mean that the flow is compressible (i.e., a rotational flow).A family of all-speed advection upstream splitting methods for hydrodynamics improves the accuracy of low-speed flows by correcting the velocity difference term with a scale of the convection speed (Liou, 2006;Shima and Kitamura, 2011;Kitamura and Shima, 2013).The MLAU scheme adopts this strategy, but uses the scale of a modified fast magnetosonic speed c u instead of the convection speed in consideration of the inequalities of the MHD eigenmodes, The HLLD scheme calculates the numerical flux of the normal momentum as x , where the both ρS 2 M and P * t terms include the velocity difference.Whereas the velocity difference term from ρS 2 M scales with the convection speed (not explicitly shown here), the term from P * t scales with the fast magnetosonic speed (third term in Equation ( 7)), and thus it is corrected as follows: The factor φ multiplied by the third term should satisfy ∝ c u /c f (c u /c f 1) to improve the accuracy of low Mach number flows and 1 (c u /c f ≥ 1) to reduce to the original scheme at high Mach numbers.We adopt the function used in Shima and Kitamura (2011), Note that the corrected pressure given by Equation ( 15) does not necessarily reduce to the original one (Equation ( 7)) even though the flow is purely one dimensional.
An asymptotic analysis by Liou (2006) provides a physical interpretation of Equation ( 15).In the following, we omit the magnetic field for simplicity.The physical variables are normalized such that each nondimensional variable remains on the order of unity, i.e., the spatial coordinates by the characteristic length L, the density by the ambient density ρ 0 , the velocity by the reference velocity u 0 , and the gas pressure by ρ 0 a 2 0 , where a 0 is the characteristic speed of sound.When one considers compressible flows, the time should be normalized by L/a 0 , and the resulting momentum equation is as follows: where M = u 0 /a 0 is the reference Mach number, and the nondimensional variables are denoted by a tilde.Expanding the pressure with respect to M , Equation ( 17) implies that for M 1, When one considers nearly incompressible flows, however, it is appropriate to normalize the time by L/u 0 , and the resulting momentum equation is as follows: which implies that the pressure is constant in space up to the first order in M for M 1, Owing to the difference between Equations ( 18) and ( 20), compressible fluid simulations may overestimate the spatial fluctuation of pressure in low Mach number flows, and their solution deviates from a correct solution with decreasing Mach number.
To cover an incompressible range with the compressible scheme, we consider the modified momentum equation: where φ ∝ M (M 1) and φ = 1 (M ≥ 1).By assuming φ is uniform over the Riemann fan, the normal velocity is evaluated from the HLL average of the modified momentum equation (as is done in the HLLD scheme): Once the normal velocity is determined, the pressure in the Riemann fan is consistently derived from the modified momentum equation to satisfy the jump condition across the S α wave (α = L or R), which is identical to Equation ( 15).The third term in Equation ( 22) is divided by φ, and one needs a cutoff Mach number to avoid division by zero (Liou, 2006).Even though division by zero is avoided, this term severely restricts the CFL condition of an explicit scheme.We therefore do not adopt Equation ( 22) to the new scheme for simplicity.Relationship between the new scheme and the original scheme is discussed in the one-dimension case.Since the normal velocity is unchanged in one dimension, the only difference is the total pressure in the Riemann fan, which is rewritten as where ρ = 2ρ L ρ R /(ρ R + ρ L ), ∆u = u R − u L , and we approximate S R − u R = u L − S L = c f for subsonic flows.Substituting Equation (24) into Equation (3) and recalling the fact that the HLLD scheme satisfies the jump conditions to calculate the numerical flux F HLLD = F (U * HLLD ), we obtain This numerical flux is interpreted as the combination of the HLLD flux and the anti-diffusion terms for the normal momentum and kinetic energy (φ−1 ≤ 0); however, it does no longer satisfy the jump conditions.
To summarize, we propose the new HLLD scheme that adopts Equations ( 9) and (15) instead of Equations ( 6) and ( 7).These modifications do not violate the preservation of the MHD discontinuities inherent in the original scheme.Hereafter, this scheme is termed a low-dissipation HLLD (LHLLD) approximate Riemann solver.

Numerical experiments
Numerical experiments are conducted to assess the capability of the LH-LLD scheme.The design of the numerical code is the same as that used in Minoshima et al. (2020); physical variables are interpolated using the second-order MUSCL scheme with a minmod limiter (van Leer, 1979), and are integrated in time through the third-order strong stability preserving Runge-Kutta method (Shu and Osher, 1988).We do not use the characteristic decomposition for low-speed flows because it may cause unphysical oscillation similar to the MLAU scheme (see Section 3.4 in Minoshima et al. (2020) for details).The solenoidal condition of the magnetic field is preserved within the machine precision by the central upwind constrained transport method (Minoshima et al., 2019).A CFL number of 0.4 is used.
We confirm that solutions with the LHLLD scheme are indistinguishable from those with the HLLD scheme in standard benchmark tests such as onedimensional shock tube and Orszag-Tang vortex problems; thus, they are not shown here.This implies that the inconsistency with respect to the jump conditions (Equation ( 25)) has little impact on at least these tests.Furthermore, the LHLLD scheme is almost identical to the HLLD scheme in the case of strongly magnetized plasma without strong grid-aligned shocks because the factors θ and φ are almost in unity throughout the domain.We then present two stringent problems conducted in Minoshima et al. (2020) to focus on the difference between the LHLLD and HLLD schemes, which include extremely low and high Mach number flows in weakly magnetized plasma.
The first problem is the two-dimensional Kelvin-Helmholtz instability (KHI) in nearly incompressible flows.The initial condition has a velocity shear u = ((V 0 /2) tanh(y/λ), 0, 0), uniform density and pressure ρ = ρ 0 , P = P 0 , and a uniform magnetic field B = B 0 (cos θ, 0, sin θ), with ρ 0 = V 0 = B 0 = λ = 1.0 and θ = 71.565• .We use γ = 2.0.To initiate the instability, we impose a small (1%) perturbation into the y-component of the velocity around the shear layer with a wavelength of 20λ equal to the fastest growing mode.The computational domain ranging from 0 ≤ x < 20λ and −10λ ≤ y < 10λ is resolved by N × N cells.The boundary condition is periodic and symmetric in the x and y directions.Figure 1(a) shows the growth rate obtained using the HLLD and LHLLD schemes at N = 64, 128, and 256 and P 0 = 500.The linear growth is shown for verification, which is obtained by solving the linearized MHD equations (Equation (68) in Minoshima et al. (2020)) and is mostly independent of the pressure for nearly incompressible flows.The LHLLD scheme converges to the linear growth rate of γ = 0.051V 0 /λ for all runs, whereas the HLLD scheme requires N = 256 to converge.Figure 1(b) shows the growth rate at N = 64 and P 0 = 500, 5000, 50000, and 500000 (corresponding Mach numbers of 1.58 × 10 −2 , 5.0 × 10 −3 , 1.58 × 10 −3 , and 5.0 × 10 −4 ).The solutions with the HLLD scheme get worse upon increasing the pressure (decreasing the Mach number), which stems from the fact that the numerical diffusion for velocity scales with a fast magnetosonic speed (Section 2.2).By contrast, the solutions with the LHLLD scheme are mostly independent of the Mach number with the help of the pressure correction, although the growth rate and the saturation level at P 0 ≥ 5000 are slightly lower than those at P 0 = 500.
The upstream Mach number is M = |u|/a = 100 and the plasma beta is β = 2P/|B| 2 = 10 5 .We use γ = 5/3.A corrugated contact discontinuity is imposed in the upstream region, y cd = Y 0 + ψ 0 cos(2πx/λ), where Y 0 = 1.0, ψ 0 = 0.1 is the corrugation amplitude, and λ = 1.0 is the wavelength.We shift a frame moving with v = −0.625,which is the interface velocity after the collision with the incident shock so that the structure of the RMI remains at approximately y = 0 throughout the simulation run.The computational domain ranging from 0 ≤ x < λ and −40λ ≤ y < 40λ is resolved by N × 80N cells, where N = 128.The boundary condition is periodic in the x direction and is fixed to the initial state in the y direction.Figure 2 shows the density profile at t = 25 obtained using (a) the HLLD scheme, (b) the LHLLD scheme without the pressure correction (φ = 1), and (c) the LHLLD scheme.The grid-scale oscillation observed at the front of the transmitted shock (y = 4) in the HLLD scheme is successfully removed with the LHLLD scheme, which is shown in Figure 3 (a).A comparison between Figure 2(a) and (b) indicates that the oscillation affects the structure until y = 1, far beyond the shock.The vortices centered around (x, y) = (0.3, 0.3) and (0.7, 0.3) are the consequence of the nonlinear development of the RMI, in which the magnetic field is amplified (Sano et al., 2012).2019)).The rotational flow is faster in the run with the LHLLD scheme than in other runs owing to the reduction of the numerical diffusion by the pressure correction, which will affect the saturation level of the magnetic field amplification.The results of KHI and RMI are in good agreement with those obtained using the MLAU scheme (Sections 4.5 and 4.6 in Minoshima et al., 2020).

Conclusion
We proposed a new low-dissipation HLLD (LHLLD) approximate Riemann solver that implements the techniques used in the MLAU scheme (Minoshima et al., 2020).The LHLLD scheme modifies the normal velocity and the total pressure in the Riemann fan to avoid a numerical shock instability and improve the accuracy of low-speed flows (Equations ( 9) and ( 15)).Stringent benchmark tests verify the capability of the LHLLD scheme.The scheme is accurate for nearly incompressible flows as long as the flow is super Alfvénic because the velocity difference term in the pressure is corrected to scale with the modified fast magnetosonic speed c u ≥ c a in consideration of the inequalities of the MHD eigenmodes.Furthermore, the baseline HLLD scheme intrinsically includes the numerical diffusion for the tangential momentum with a scale of Alfvén speed (Minoshima et al., 2020).We call this capability as quasi-all speeds for super-Alfvénic flows.
The clear advantage of the LHLLD scheme over the MLAU scheme is its ease of use for current HLLD scheme users; one only needs to implement two factors θ and φ (Equations ( 10) and ( 16)).One can independently implement these factors depending on the purpose, as is demonstrated in Section 3.For example, one can use only the pressure correction φ when the flow is known to be subsonic in a whole domain.Additional computational costs are low.The ratio of the CPU time for the serial computation of a two-dimensional problem in our code is (HLLD:LHLLD:MLAU)=(1:1.03:1.13),although the computational efficiency depends on the design of other procedures (e.g., interpolation) and the level of optimization.The proposed techniques can be implemented to the HLLC scheme for hydrodynamics (Toro et al., 1994) as well, and will be applicable to the simulations with a general equation of state.We hope that current HLLD scheme users can easily enjoy the novel ability of the new scheme to tackle MHD systems, including both high and low Mach number flows.
The source code written in C programming language can be downloaded from the GitHub website1 .Figure 3: One-dimensional profiles of (a) the density along the x-direction at y = 4.0, and (b) the rotational velocity component along the y-direction at x = 0.3.The cyan, blue, and red lines are the solutions obtained using the HLLD scheme, the LHLLD scheme without the pressure correction, and the LHLLD scheme.

Figure 1 :
Figure 1: Time profile of the Fourier amplitude of the y-component of the velocity in the Kelvin-Helmholtz instability.(a) Dependence on grid resolution.(b) Dependence on the initial pressure.The cold-and warm-colored lines are the solutions obtained using the HLLD and LHLLD schemes.The dashed lines indicate the solution obtained with the linear theory.

Figure 2 :
Figure 2: Density profile in the Richtmyer-Meshkov instability at t = 25 obtained using (a) the HLLD scheme, (b) the LHLLD scheme without the pressure correction, and (c) the LHLLD scheme.The arrows represent the rotational velocity component.