© ESO 2019

1. Introduction

Many astrophysical processes involving plasma flows can be modeled using equations of magnetohydrodynamics (MHD). These equations arise from the coupling between equations of fluid mechanics and Maxwell’s equations, and they allow us to model the plasma motion under certain conditions that can be assumed in multiple astrophysical phenomena. Magnetohydrodynamics equations do not have a general exact solution, therefore they have to be solved using alternative techniques. Thus, numerical simulations play a very important role in obtaining solutions for many processes that are modeled by these equations.

The dynamics of the solar atmosphere is often analyzed using MHD numerical simulations to study multiple phenomena that occur there. Some of these phenomena evolve in the context of a very gravitationally stratified medium with strong pressure and temperature gradients (e.g. Aschwanden 2005; Mei et al. 2012; Murawski et al. 2013; Krause et al. 2018). The numerical codes used in these simulations must be able to maintain the hydrostatic equilibrium of the background medium at some precision level so as not to introduce spurious results that affects the actual physical behavior. Sometimes a certain level of spurious accelerations is acceptable so that conventional numerical methods can be used in relatively fine grids, but in other situations it is necessary to minimize the presence of nonphysical velocities, which could mask or significantly alter the studied processes. For this reason, special modifications to the standard numerical schemes have been proposed in order to satisfy the discrete equilibrium in hydrostatic conditions exactly. We emphasize that hydrostatic equilibrium preservation does not only concern stationary problems, but can be a crucial aspect for many dynamic processes that occur in the context of an hydrostatically balanced environment. In these cases, a numerically induced velocity can strongly affect or even completely overlap the studied dynamics leading to erroneous results. Some examples of these phenomena are wave propagation in stratified magnetic atmospheres (Fuchs et al. 2010), low-amplitude wave propagation in helioseismology (Kobanov et al. 2008), or the interaction between coronal and chromospheric waves Krause et al. (2018).

The capability of a scheme to numerically satisfy an equilibrium condition is known as a well-balanced property after the work by Greenberg & Leroux (1996), who introduced this terminology to characterize those schemes that are able to exactly satisfy, in a numerical discrete sense, a defined stationary state. In the particular case of the hydrostatic equilibrium, some well-balanced schemes have been developed to model numerically the variants of that condition (e.g. Botta et al. 2004; Audusse et al. 2004; Felipe et al. 2010; Fuchs et al. 2010; Xing & Shu 2013; Käppeli & Mishra 2014, 2016). Among these, schemes that use a local hydrostatic reconstruction are the simplest since this technique does not require a modification of the system of governing equations nor the central aspects of the standard numerical schemes. This allows for a relatively simple implementation of the method in existing codes with the conventional Godunov-type schemes mostly used in the astrophysical community (e.g. Fryxell et al. 2000; Pen et al. 2003; Ziegler 2004; Mignone et al. 2007).

In this work we choose the local hydrostatic reconstruction approach following the ideas of works by Fuchs et al. (2010), Käppeli & Mishra (2014, 2016). We propose some contributions that allow us to simulate equilibrium atmospheres in more general conditions than the isothermal or isentropic cases, for both constant and variable gravity, without the need to define new variables nor solve additional equations.

The paper is organized as follows. In Sect. 2 we present the governing equations of the MHD problem and we show the analytic conditions of hydrostatic equilibrium. In Sect. 3 we briefly summarize the numerical methods used in this work and discuss the problem of hydrostatic equilibrium preservation for these types of numerical schemes. In Sect. 4 we present an extension of a local hydrostatic reconstruction technique and in Sect. 5 we detail its implementation in the FLASH¹ code. Finally, in Sect. 6 we show some results to verify our formulation and in Sect. 7 we present the main conclusions of this work.

2. Governing equations

Considering an ideal plasma (inviscid and without magnetic resistivity), the ideal MHD equations in conservative form and in CGS units are written as

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathbf{\nabla }·\left(\rho \mathit{u}\right)=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial \rho \mathit{u}}{\partial t}+\mathrm{\nabla }·{\displaystyle (\rho \mathit{u}\mathit{u}-\frac{1}{4\pi }\mathit{B}\mathit{B})}+\mathrm{\nabla }p+\frac{1}{8\pi }\mathbf{\nabla }{B}^2=\rho \mathit{g},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial \mathit{B}}{\partial t}+\mathbf{\nabla }·(\mathit{u}\mathit{B}-\mathit{B}\mathit{u})=\mathbf{0}\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial E}{\partial t}+\mathrm{\nabla }·[(E+p+\frac{B^2}{8\pi })\mathit{u}-\frac{1}{4\pi }(\mathit{u}·\mathit{B})\mathit{B}]=\rho \mathit{g}·\mathit{u},\hfill \end{array}$$(4)

where ρ is the plasma density, u = (u, v, w) is the velocity, B = (B_x, B_y, B_z) is the magnetic field, p is the thermal pressure, and E is the total energy (per unit volume), i.e.,

$$\begin{array}{c}\hfill E=\rho ϵ+\frac{1}{2}\rho (u^2+v^2+w^2)+\frac{B^2}{8\pi },\end{array}$$

where ϵ is the specific internal energy. The vector g = (g_x, g_y, g_z) is the gravitational acceleration.

The set of Eqs. (1)–(4) must be closed by means of an equation of state that links thermal pressure with other thermodynamic variables (p = f(ρ, ϵ)), In addition, the divergence-free condition of the magnetic field (∇ ⋅ B = 0) must be fulfilled everywhere.

The conservative MHD equations constitute a conservation law that expresses the balance between the variation, flux and source term of each conservative variable,

$$\begin{array}{c}\hfill \frac{\partial \mathit{U}}{\partial t}+\frac{\partial \mathit{F}}{\partial x}+\frac{\partial \mathit{G}}{\partial y}+\frac{\partial \mathit{H}}{\partial z}=\mathit{S},\end{array}$$(5)

where

$$\begin{array}{c}\hfill \mathit{U}=\left[\begin{array}{c}\rho \\ \rho u\\ \rho v\\ \rho w\\ B_x\\ B_y\\ B_z\\ E\end{array}\right]\end{array}$$(6)

is the vector of conservative variables, and

$$\begin{array}{cc}& \mathit{F}=\left[\begin{array}{c}\rho u\\ \rho u^2+p+\frac{B^2}{8\pi }-\frac{B_x^2}{4\pi }\\ \rho uv-\frac{1}{4\pi }{B}_y{B}_x\\ \rho uw-\frac{1}{4\pi }{B}_z{B}_x\\ 0\\ u{B}_y-v{B}_x\\ u{B}_z-w{B}_x\\ (E+p+\frac{B^2}{8\pi })u-\frac{1}{4\pi }{B}_x(u{B}_x+v{B}_y+w{B}_z)\end{array}\right],\hfill \\ \hfill & \mathit{G}=\left[\begin{array}{c}\rho v\\ \rho vu-\frac{1}{4\pi }{B}_x{B}_y\\ \rho v^2+p+\frac{B^2}{8\pi }-\frac{B_y^2}{4\pi }\\ \rho vw-\frac{1}{4\pi }{B}_z{B}_y\\ v{B}_x-u{B}_y\\ 0\\ v{B}_z-w{B}_y\\ (E+p+\frac{B^2}{8\pi })v-\frac{1}{4\pi }{B}_y(u{B}_x+v{B}_y+w{B}_z)\end{array}\right],\hfill \\ \hfill & \mathit{H}=\left[\begin{array}{c}\rho w\\ \rho wu-\frac{1}{4\pi }{B}_x{B}_z\\ \rho wv-\frac{1}{4\pi }{B}_y{B}_z\\ \rho w^2+p+\frac{B^2}{8\pi }-\frac{B_z^2}{4\pi }\\ w{B}_x-u{B}_z\\ w{B}_y-v{B}_z\\ 0\\ (E+p+\frac{B^2}{8\pi })v-\frac{1}{4\pi }{B}_z(u{B}_x+v{B}_y+w{B}_z)\end{array}\right]\hfill \end{array}$$(7)

are the flux vectors in x, y, and z direction in a Cartesian coordinate system. In this case, the source term is only given by the components of the gravitational acceleration

$$\begin{array}{c}\hfill \mathit{S}=\left[\begin{array}{c}0\\ \rho g_x\\ \rho g_y\\ \rho g_z\\ 0\\ 0\\ 0\\ \rho (g_{x}u+g_{y}v+g_{z}w)\end{array}\right]·\end{array}$$(8)

There are several numerical methods that were specifically developed to solve MHD equations (e.g. Colella & Woodward 1984; Balsara & Spicer 1999; Powell et al. 1999; Gardiner & Stone 2005), which were implemented in many computational codes with different features and capabilities. In general, these codes use some extension of the Godunov method, which is based on the finite volume method (FVM). In this type of scheme the physic domain is subdivided into small control volumes (cells) on each of which the conservation laws of mass, momentum, magnetic flux, and energy must be satisfied, i.e., the balance between the inbound and outbound fluxes of each conservative property, its source term, and its rate of change must be fulfilled on each cell. This integral approach allows us to capture discontinuous solutions such as shocks waves and contact discontinuities that can be present in high energy plasma flows, as regularly occurs in the solar atmosphere (blasts, coronal mass ejections, etc.).

Usually, the main challenge for numerically solving MHD equations lies in the evolution of the magnetic field, which has to satisfy the divergence-free constraint (∇ ⋅ B = 0) at all times. This requirement does not allow the use of direct extensions of fluid codes to solve the MHD equations, since the discrete evaluation of magnetic field gradients introduces artificial divergence even when the initial magnetic field is solenoidal, which violates the law of nonexistence of magnetic monopoles. The non-preservation of the divergence-free condition causes discrepancies between the conservative and primitive forms of the MHD equations. In addition, the Lorentz force is not normal to the magnetic field (Brackbill & Barnes 1980) and the jump conditions through MHD shocks can be incorrect under certain conditions (Tóth 2000). For these reasons, many efforts have been carried out to obtain numerical schemes that hold the magnetic field divergence in very small levels throughout the simulation. In the Godunov-type schemes three methods to attack this problem can be distinguished: “clean” the magnetic field of any spurious divergence after each time step using a Hodge projection (Gardiner & Stone 2005); expand the set of governing equations, by incorporating an equation for the divergence evolution designed to eliminate the error accumulation (Powell et al. 1999); and write the finite difference equations for the magnetic field such that the magnetic flux is explicitly conserved to preserve the divergence-free condition at the discretization level (e.g. Evans & Hawley 1988; Lee & Deane 2009). With some of these techniques and some Riemann solver developed for MHD equations (e.g. Powell 1994; Cargo & Gallice 1997; Gurski 2004; Miyoshi & Kusano 2005), the numerical solution of Eqs. (1)–(4) is obtained in a similar way as the gas dynamic equations. In our study, we are interested in the “gas dynamic part” of the MHD problem, since we are proposing a methodology to satisfy hydrostatic equilibrium in atmospheres at rest, which is a fluid problem. Therefore, even though we only consider the gas dynamic variables in the MHD equations, this approach can be applied on any Godunov-type MHD code. In addition, for simplicity we considered a 2D domain taking into account that the extension to 3D problems is trivial.

Then, for our analysis we used the 2D Euler equations, which can be viewed as a particular case of Eqs. (1)–(4) when the magnetic field is zero (B = 0). The resulting equations after eliminating the terms associated with the magnetic field and the components in the z-direction are written as

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathbf{\nabla }·\left(\rho \mathit{u}\right)=0,\hfill \end{array}$$(9)

$$\begin{array}{cc}& \frac{\partial \rho \mathit{u}}{\partial t}+\mathrm{\nabla }·\rho \mathit{u}\mathit{u}+\mathrm{\nabla }p=\rho \mathit{g},\hfill \end{array}$$(10)

$$\begin{array}{cc}& \frac{\partial E}{\partial t}+\mathrm{\nabla }·(E+p)\mathit{u}=\rho \mathit{g}·\mathit{u}.\hfill \end{array}$$(11)

The vector of conservative variables and the vectors of fluxes for the 2D Euler equations are written as

$$\begin{array}{c}\hfill \mathit{U}=\left[\begin{array}{c}\rho \\ \rho u\\ \rho v\\ E\end{array}\right],\mathit{F}=\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho vu\\ (E+p)u\end{array}\right],\mathit{G}=\left[\begin{array}{c}\rho v\\ \rho uv\\ \rho v^2+p\\ (E+p)v\end{array}\right]·\end{array}$$(12)

In order to simplify the analysis, the coordinate system is oriented so that the gravitational acceleration direction is parallel to the y coordinate pointing downward. Thus, for a gravity vector g of magnitude g, the gravitational source term is

$$\begin{array}{c}\hfill \mathit{S}=\left[\begin{array}{c}0\\ 0\\ -\rho g\\ -\rho gv\end{array}\right]·\end{array}$$(13)

The nonconservative form of the system (9)–(11) is written in terms of the Jacobian matrices A_x and A_y, and the vector of primitive variables V as

$$\begin{array}{c}\hfill \frac{\partial \mathit{V}}{\partial t}+{\mathbf{A}}_x\frac{\partial \mathit{V}}{\partial x}+{\mathbf{A}}_y\frac{\partial \mathit{V}}{\partial y}={\mathit{S}}_V,\end{array}$$(14)

where

$$\begin{array}{cc}& \mathit{V}=\left[\begin{array}{c}\rho \\ u\\ v\\ p\end{array}\right],{\mathbf{A}}_x={\left(\frac{\partial \mathit{U}}{\partial \mathit{V}}\right)}^{-1}\frac{\partial \mathit{F}}{\partial \mathit{V}},{\mathbf{A}}_y={\left(\frac{\partial \mathit{U}}{\partial \mathit{V}}\right)}^{-1}\frac{\partial \mathit{G}}{\partial \mathit{V}},\hfill \\ \hfill & {\mathit{S}}_V=\left[\begin{array}{c}0\\ 0\\ -g\\ \rho gv(\gamma -1)\end{array}\right]·\hfill \end{array}$$(15)

In the hydrostatic equilibrium condition the gravitational forces are balanced by the pressure forces so that the fluid remains at rest if there are no perturbations. This condition implies that the velocity in the gravity direction is zero (v = 0 in this case) and pressure and density vary only in that direction according to the height equation

$$\begin{array}{c}\hfill \frac{\partial p}{\partial y}=-\rho g.\end{array}$$(16)

The solution of the height equation is given by the prescribed distribution of one of the thermodynamic variables that has to be known to obtain a unique solution for the equilibrium. In the modeling of stellar atmospheres it is usual to prescribe the radial temperature distribution, since this profile can generally be estimated. In the particular case of the Sun, there are several models for the temperature variation in the different layers of the solar atmosphere (Aschwanden 2005). The link between the temperature and height Eq. (16) is given by the equation of state. Taking into account that the assumption of an ideal gas atmosphere is widely adopted in the study of the solar atmosphere, we focused our analysis on this model since it is a simple and reasonable supposition that allowed us to obtain explicit formulations for the numerical scheme. Therefore, for the equation of state of ideal gases we have

$$\begin{array}{c}\hfill p=\rho RT=(\gamma -1)\rho ϵ,\end{array}$$(17)

where R is the particular gas constant (relation between the universal gas constant and the gas molar mass) and γ is the specific heats relation. With this assumption, Eq. (16) reads

$$\begin{array}{c}\hfill \frac{\mathrm{d}p}{\mathrm{d}y}=-\frac{\mathit{pg}}{\mathit{RT}},\end{array}$$(18)

where the distribution T(y) is known. Integrating the last expression gives

$$\begin{array}{c}\hfill p(y)=p_0\exp (-{{\displaystyle \int }}_{y_0}^y\frac{g(y^{\prime })}{RT(y^{\prime })}\mathrm{d}{y}^{\prime }),\end{array}$$(19)

where p₀ is a reference pressure that is known at the reference height y₀.

Equation (19) gives the hydrostatic pressure distribution for a given temperature profile T(y) for an ideal gas atmosphere. In those regions where the velocity is zero and the thermodynamic variables satisfy Eq. (19), the fluid remains at rest if there are no perturbations. Under certain conditions such as strong pressure gradients, large height scales, or in long-lasting phenomena, the conventional Godunov-type numerical schemes fail to preserve the hydrostatic equilibrium; this is because the exponential variation of pressure in the gravity direction cannot be captured by the classical reconstruction-evolution methods even when they are of high order. As a consequence, unbalanced momentum fluxes arise in the gravity direction, which induce spurious vertical velocities that constantly increase with the advance of time. Although the magnitude of these spurious velocities can be reduced, thereby enhancing the grid resolution (provided that the implemented method is convergent), in some problems this strategy is not possible because of the strong penalty in computational cost that implies the excessive increment of computing cells.

3. Godunov-type numerical schemes

Co-located FVMs assume that inside each control subvolume, i.e., inside each cell, the variables are constant and equal to an averaged quantity. By definition, inside a cell identified by the subindices i, j in a structured Cartesian grid of rectangles of dimension Δx × Δy (see Fig. 1), the discrete vector of conservative variables for a given time tⁿ is written as

$$\begin{array}{c}\hfill {\mathit{U}}_{i,j}^n=\frac{1}{\mathrm{\Delta }x\mathrm{\Delta }y}{{\displaystyle \int }}_{x_i-\frac{\mathrm{\Delta }x}{2}}^{x_i+\frac{\mathrm{\Delta }x}{2}}{{\displaystyle \int }}_{y_j-\frac{\mathrm{\Delta }y}{2}}^{y_j+\frac{\mathrm{\Delta }y}{2}}\mathit{U}(x,y,t^n)\mathrm{d}y\mathrm{d}x,\end{array}$$(20)

Fig. 1. Definitions for a generic cell i, j in a structured Cartesian grid.

where x_i and y_j are the coordinates of the cell center.

The evolution of the discrete variable U_i, j from the time tⁿ to the time t^n + 1 = tⁿ + Δt is obtained by considering the integral form of the conservation law (5) applied on the control volume defined by the corner points $x_{i\mp \frac{1}{2}}=x_i\mp \mathrm{\Delta }x/2$ and ${\mathit{y}}_{j\mp \frac{1}{2}}={\mathit{y}}_j\mp \mathrm{\Delta }\mathit{y}/2$ as

$$\begin{array}{cc}\hfill {\mathit{U}}_{i,j}^{n+1}=& {\mathit{U}}_{i,j}^n+\mathrm{\Delta }t{\mathit{S}}_{i,j}^{n+\frac{1}{2}}-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}({\mathit{F}}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}-{\mathit{F}}_{i-\frac{1}{2},j}^{n+\frac{1}{2}})-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }y}({\mathit{G}}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}-{\mathit{G}}_{i,j-\frac{1}{2}}^{n+\frac{1}{2}}),\hfill \end{array}$$(21)

where the discrete fluxes are averaged over the cell faces of coordinates $x_{i\mp \frac{1}{2}}$ and ${\mathit{y}}_{j\mp \frac{1}{2}}$ and through the time interval tⁿ ≤ t ≤ tⁿ + Δt, i.e.,

$$\begin{array}{cc}& {\mathit{F}}_{i\mp \frac{1}{2},j}^{n+\frac{1}{2}}=\frac{1}{\mathrm{\Delta }t\mathrm{\Delta }y}{{\displaystyle \int }}_{t^n}^{t^n+\mathrm{\Delta }t}{{\displaystyle \int }}_{y_j-\frac{\mathrm{\Delta }y}{2}}^{y_j+\frac{\mathrm{\Delta }y}{2}}\mathit{F}(x_i\mp \frac{\mathrm{\Delta }x}{2},y,t)\mathrm{d}y\mathrm{d}t,\hfill \\ \hfill & {\mathit{G}}_{i,j\mp \frac{1}{2}}^{n+\frac{1}{2}}=\frac{1}{\mathrm{\Delta }t\mathrm{\Delta }x}{{\displaystyle \int }}_{t^n}^{t^n+\mathrm{\Delta }t}{{\displaystyle \int }}_{x_i-\frac{\mathrm{\Delta }x}{2}}^{x_i+\frac{\mathrm{\Delta }x}{2}}\mathit{G}(x,y_j\mp \frac{\mathrm{\Delta }y}{2},t)\mathrm{d}x\mathrm{d}t.\hfill \end{array}$$(22)

3.1. Numerical flux

Typically, the averaged fluxes are approximated by numerical fluxes $\stackrel{\mathbf{\sim }}{\mathit{F}}$ and $\stackrel{\mathbf{\sim }}{\mathit{G}}$, which must satisfy the conservation property. In the Godunov-type schemes the numerical fluxes are evaluated by constructing Riemann problems (RP) at the interfaces of adjacent cells, where discontinuities can exist owing to the difference between the averaged quantities of each cell. The variable values on each side of the intercell boundary, namely, the RP states, are defined by some reconstruction technique that is used to extrapolate the cell-centered values of the variables to the cell faces. With the solution of the RP the numerical fluxes are evaluated by means of the physic fluxes (e.g. Toro 2009)

$$\begin{array}{cc}& {\stackrel{\mathbf{\sim }}{\mathit{F}}}_{i\mp \frac{1}{2},j}^{n+\frac{1}{2}}=\mathit{F}\left({\stackrel{\mathbf{\sim }}{\mathit{V}}}_{i\mp \frac{1}{2},j}^{n+\frac{1}{2}}\right),\hfill \\ \hfill & {\stackrel{\mathbf{\sim }}{\mathit{G}}}_{i,j\mp \frac{1}{2}}^{n+\frac{1}{2}}=\mathit{G}\left({\stackrel{\mathbf{\sim }}{\mathit{V}}}_{i,j\mp \frac{1}{2}}^{n+\frac{1}{2}}\right),\hfill \end{array}$$(23)

where ${\stackrel{\mathbf{\sim }}{\mathit{V}}}_{i\mp \frac{1}{2},j}^{n+\frac{1}{2}}$ and ${\stackrel{\mathbf{\sim }}{\mathit{V}}}_{i,j\mp \frac{1}{2}}^{n+\frac{1}{2}}$ are the solutions of the RP on faces normal to the x-direction and the y-direction, respectively, i.e.,

$$\begin{array}{cc}\hfill {\stackrel{\mathbf{\sim }}{\mathit{V}}}_{i-\frac{1}{2},j}^{n+\frac{1}{2}}& =\mathrm{RP}\{{\mathit{V}}_{i-1,j,E}^{n+\frac{1}{2}},{\mathit{V}}_{i,j,W}^{n+\frac{1}{2}}\},{\stackrel{\mathbf{\sim }}{\mathit{V}}}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}=\mathrm{RP}\{({\mathit{V}}_{i,j,E}^{n+\frac{1}{2}},{\mathit{V}}_{i+1,j,W}^{n+\frac{1}{2}}\},\hfill \\ \hfill {\stackrel{\mathbf{\sim }}{\mathit{V}}}_{i,j-\frac{1}{2}}^{n+\frac{1}{2}}& =\mathrm{RP}\{{\mathit{V}}_{i,j-1,N}^{n+\frac{1}{2}},{\mathit{V}}_{i,j,S}^{n+\frac{1}{2}}\},{\stackrel{\mathbf{\sim }}{\mathit{V}}}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}=\mathrm{RP}\{{\mathit{V}}_{i,j,N}^{n+\frac{1}{2}},{\mathit{V}}_{i,j+1,S}^{n+\frac{1}{2}}\}.\hfill \end{array}$$

The variable averaging over the time interval (step $n+\frac{1}{2}$) depends on the implemented temporal integration scheme. The subindices E, W, N, and S (east, west, north, and south) denote the relative positions of the cell faces with respect to the cell center, according to the nomenclature usually used (see Fig. 1).

The way in which we extrapolated the cell-centered values to the faces to construct the RPs defines the spatial approximation order of a Godunov-type scheme. The MUSCL (monotone upwind scheme for conservation laws) schemes, initially proposed by van Leer (1979) reconstruct the variable inside the cell using the values at neighboring cells to calculate an interpolating polynomial, whose order is given by the number of surrounding cell values used in the calculation. This polynomial, whose integration over the cell must satisfy the definition (20) to verify the conservation property, allows us to evaluate the reconstructed variables at face positions. In addition, a flux limiter (FL) function has to be used to preserve the monotonicity of the solution (Toro 2009).

The solution of the RPs at each interface is obtained using some approximate Riemann solver, which must satisfy the consistency property that requires that in a constant state through the interface, the numerical flux has to be equal to the physics flux. This means that the solution of the Riemann solver for a constant state has to be equal to the constant values. On the other hand, in many applications it is also necessary that the Riemann solver can capture contact discontinuities, which implies that the numerical flux for a constant pressure state with zero velocity normal to the face only produces a flux momentum component that is normal to the face and equal to the pressure, i.e.,

$$\begin{array}{cc}& \mathit{F}\left(\mathrm{RP}\{\left[\begin{array}{c}{\rho }_{E^{-}}\\ 0\\ 0\\ p\end{array}\right],\left[\begin{array}{c}{\rho }_{W^{+}}\\ 0\\ 0\\ p\end{array}\right]\}\right)=\left[\begin{array}{c}0\\ p\\ 0\\ 0\end{array}\right],\hfill \\ \hfill & \mathit{G}\left(\mathrm{RP}\{\left[\begin{array}{c}{\rho }_{N^{-}}\\ 0\\ 0\\ p\end{array}\right],\left[\begin{array}{c}{\rho }_{S^{+}}\\ 0\\ 0\\ p\end{array}\right]\}\right)=\left[\begin{array}{c}0\\ 0\\ p\\ 0\end{array}\right]·\hfill \end{array}$$(24)

In other words, in a constant pressure state with zero velocity, the numerical flux has to be equal to the physical flux independent of the density values on each side of the interface. In general, this property is very important to improve the hydrostatic equilibrium preservation of conventional Godunov-type schemes.

3.2. Source terms

The numerical evaluation of the source term is also made by means of an averaging over the cell and the time interval, i.e.,

$$\begin{array}{c}\hfill {\mathit{S}}_{i,j}^{n+\frac{1}{2}}=\frac{1}{\mathrm{\Delta }t\mathrm{\Delta }x\mathrm{\Delta }y}{{\displaystyle \int }}_{t^n}^{t^n+\mathrm{\Delta }t}{{\displaystyle \int }}_{x_i-\frac{\mathrm{\Delta }x}{2}}^{x_i+\frac{\mathrm{\Delta }x}{2}}{{\displaystyle \int }}_{y_j-\frac{\mathrm{\Delta }y}{2}}^{y_j+\frac{\mathrm{\Delta }y}{2}}\mathit{S}(x,y,t)\mathrm{d}y\mathrm{d}x\mathrm{d}t.\end{array}$$(25)

For the particular case in which the source term is only due to gravity, which does not depend on time, the averaging is usually replaced by the cell-centered evaluation of the gravitational acceleration, i.e.,

$$\begin{array}{c}\hfill {\mathit{S}}_{i,j}^{n+\frac{1}{2}}=\left[\begin{array}{c}0\\ 0\\ -{\rho }_{i,j}^{n+\frac{1}{2}}g(y_j)\\ -{(\rho v)}_{i,j}^{n+\frac{1}{2}}g(y_j)\end{array}\right]·\end{array}$$(26)

This expression is the same as Eq. (25) when gravity is constant and is approximate to that average for variable gravity.

3.3. Hydrostatic equilibrium preservation

In this section we describe some problems concerning the hydrostatic equilibrium preservation in standard Godunov-type schemes. For the analysis, we start from a hydrostatic configuration where the pressure and density distributions satisfy Eq. (16) and the velocity is zero everywhere (u = 0). Naturally, this is a stationary problem, but we want to know how explicit Godunov-type methods, such as those usually implemented in astrophysical simulation, behave under this situation.

Numerically, the hydrostatic condition manifests as follows: Taking into account that there are no source terms acting in that x-direction and that hydrostatic pressure and density are constant with respect to x in our model with gravity parallel to the y-direction, the condition of zero velocity implies that the net flux $({\stackrel{\mathbf{\sim }}{\mathit{F}}}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}-{\stackrel{\mathbf{\sim }}{\mathit{F}}}_{i-\frac{1}{2},j}^{n+\frac{1}{2}})$ is zero; this happens in any trivial problem with constant variable distributions without source terms. On the other hand, in the y-direction pressure and density are not constant but present an exponential decay with height, thus the net flux $({\stackrel{\mathbf{\sim }}{\mathit{G}}}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}-{\stackrel{\mathbf{\sim }}{\mathit{G}}}_{i,j-\frac{1}{2}}^{n+\frac{1}{2}})$ is not zero for u = 0. The hydrostatic configuration remains at rest only if the gravitational source terms exactly cancel the net flux, i.e.,

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Delta }y}({\stackrel{\mathbf{\sim }}{\mathit{G}}}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}-{\stackrel{\mathbf{\sim }}{\mathit{G}}}_{i,j-\frac{1}{2}}^{n+\frac{1}{2}})={\mathit{S}}_{i,j}^{n+\frac{1}{2}}.\end{array}$$(27)

For the hydrostatic equilibrium to be numerically preserved, Eq. (27) should be reduced to a discrete form of the height Eq. (16) so that the velocity remains zero and the thermodynamic variables are unchanged. This implies that the y-momentum component of the net numerical flux must be canceled by the corresponding source term, i.e.,

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Delta }y}({\stackrel{\sim }{p}}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}-{\stackrel{\sim }{p}}_{i,j-\frac{1}{2}}^{n+\frac{1}{2}})=-{(\rho g)}_{i,j}^{n+\frac{1}{2}},\end{array}$$(28)

where ${\stackrel{\sim }{p}}_{i,j\pm \frac{1}{2}}^{n+\frac{1}{2}}$ should be the only non-zero components of the numerical fluxes ${\stackrel{\mathbf{\sim }}{\mathit{G}}}_{i,j\pm \frac{1}{2}}^{n+\frac{1}{2}}$.

Equation (28) shows the reason by which standard Godunov-type schemes cannot preserve a hydrostatic equilibrium condition. Firstly, we see that the y-momentum component of the numerical flux $\stackrel{\mathbf{\sim }}{\mathit{G}}$ must only be due to pressure since velocity has to be zero. This requires that the RP solution on faces normal to y should not produce non-zero velocity components in that direction (i.e., ${\stackrel{\sim }{\mathit{v}}}_{i,j\mp \frac{1}{2}}^{n+\frac{1}{2}}=0$) in order for Eq. (27) to reduces to Eq. (28). Secondly, the numerical evaluation of the source term −ρg at each cell must be equal to the y-momentum component of the net numerical flux divided by Δy.

For the first aspect, we need to construct a RP with a constant pressure state through the cell face when the hydrostatic equilibrium is satisfied, with which we can obtain a solution of the form of Eq. (24) for flux $\stackrel{\mathbf{\sim }}{\mathit{G}}$. This means that, under a hydrostatic equilibrium condition, the extrapolation in the y-direction of the pressure cell-centered values to a face should give the same value either using the below or above cell (i.e., $p_{i,j-1,N}^{n+\frac{1}{2}}=p_{i,j,S}^{n+\frac{1}{2}}$). Conventional polynomial MUSCL reconstruction techniques do not allow us to produce such a constant pressure state because they are not able to capture the exponential decay of the variable. Consequently, a pressure jump is generated at the interface (see Fig. 2 for the examples of a first and second order reconstruction), which induces a non-zero velocity. We must note that an increase of the reconstruction order using higher order polynomials, even though it reduces the difference between the extrapolated pressures on each side of the interface, always presents errors with respect to the actual exponential distribution. The strategy usually used to reduce the spurious velocities caused by these errors in conventional codes is to increase the grid resolution, but sometimes this is not possible owing to the strong impact on computational cost that is generated under certain conditions (for instance, very large pressure gradients in small regions of the domain).

Fig. 2. First (red lines) and second (blue lines) order standard reconstruction in a hydrostatic pressure distribution. The solid black line represents the equilibrium hydrostatic pressure p0(y).

The second aspect, related to the numerical evaluation of the gravitational source term, also introduces errors in the hydrostatic equilibrium preservation of standard schemes. It is clear that the y-momentum component of the source term evaluated by Eq. (26) is not equal to the right-hand side of Eq. (28) in the general case. Hence, even if we could satisfy the required conditions for the left-hand side of the discrete height equation, the numerical scheme would not preserve the hydrostatic equilibrium if the source term $-{(\rho g)}_{i,j}^{n+\frac{1}{2}}$ is not correctly evaluated. Naturally, these errors can be reduced by improving the grid resolution but, as mentioned, this option is not always available.

As expressed in the previous paragraphs, it is clear that the conventional Godunov-type finite volume schemes are not able to reproduce a discrete hydrostatic equilibrium condition exactly. A good analysis referring to this issue is found in Zingale et al. (2002), who studied the effect of the scheme type, the treatment of boundary conditions, and the influence of the grid resolution in the modeling of gravitational stratified mediums.

As we mentioned in the Introduction, a well-balanced numerical scheme is able to numerically reproduce a given stationary solution. For the particular case of hydrostatic equilibrium preservation, several techniques with different levels of difficulty can be applied. An option is to apply a splitting on the variables, which are separated in an “equilibrium part”, which remains unchanged and corresponds to the background medium, and in a “perturbed part”, which is evolved (Botta et al. 2004; Felipe et al. 2010; Tarr et al. 2016). Although this technique would allow us to obtain well-balanced schemes for general configurations, it requires the modification of the governing equations and the Riemann solvers when nonlinear approximate methods are used (Miyoshi et al. 2010). These requirements make the implementation of the splitting variables methodology difficult in existing codes. A simpler technique makes use of a local hydrostatic reconstruction by which the extrapolation of the cell-centered pressure values to construct the RPs at each intercell face is calculated considering the gravity action (Audusse et al. 2004; Fuchs et al. 2010; Xing & Shu 2013; Käppeli & Mishra 2014, 2016). This approach only requires the modification of the form in which the face-centered values and the gravitational source terms are calculated without the need to change the main aspects of the code, such as the Riemann solvers or the numerical integrators. The local hydrostatic reconstruction technique allows us to obtain Riemann states with constant pressure when a hydrostatic equilibrium exists, with which is possible to satisfy Eq. (28) at round-off orders when a consistent evaluation of the gravitational source term is made. In the next section we present a new local hydrostatic reconstruction formulation based on the ideas by Fuchs et al. (2010), Käppeli & Mishra (2016).

4. A new well-balanced Godunov-type scheme

The goal of a well-balanced scheme is to exactly satisfy (to numerical level) a given equilibrium or stationary condition. This means that possible spurious results remain to round-off orders and do not propagate or increase with the simulation progression. In addition, it is expected that a well-balanced scheme presents a better behavior than a conventional method similar to the equilibrium conditions for which it was formulated.

In Godunov-type schemes, a well-balanced scheme for gravitationally stratified atmospheres can be obtained using a local hydrostatic reconstruction. This consists in considering the gravity action in the extrapolation of cell-centered pressure values to the faces normal to the gravity direction. Recalling our assumption of an ideal gas atmosphere, this calculation is made using Eq. (19), choosing as the reference pressure the known cell-centered pressure value p_i, j, i.e. inside the i, j cell we have

$$\begin{array}{c}\hfill {p_0}_{i,j}(y)=p_{i,j}\exp (-{{\displaystyle \int }}_{y_j}^y\frac{g(y^{\prime })}{RT(y^{\prime })}\mathrm{d}{y}^{\prime }),y_j-\frac{\mathrm{\Delta }y}{2}\le y\le y_j+\frac{\mathrm{\Delta }y}{2}·\end{array}$$(29)

This equation requires that the temperature distribution inside the cell be known. For this estimation we propose a linear reconstruction as shown below.

As is known, the gravity acceleration can be constant or variable depending on the height scale of the problem. For the second case, we have the gravitation law

$$\begin{array}{c}\hfill g(y)=\frac{\mathit{GM}}{{(R_0+y)}^2},\end{array}$$(30)

where G is the gravitational constant, M is the mass of the Sun, and R₀ is the radial distance between the center of the Sun and the origin of the coordinate system (y = 0).

With the pressure reconstruction given by Eq. (29) it could be possible to construct RPs of the form of Eq. (24) when the variables are in hydrostatic conditions, provided that the temperature distribution used for the calculation is the actual distribution. In that case, we obtained p_{i, j − 1, N} = p_i, j, S for all j because the global pressure distribution responds to the same decay law than the local pressure, thus any contact discontinuity capturing Riemann solver gives the expected solution for the numerical fluxes. In this way, with the local hydrostatic reconstruction we can reach the first requirement established by Eq. (28). The second requirement for hydrostatic equilibrium preservation is related to the numerical evaluation of the y-momentum component of the gravitational source term. For this, we consider the height Eq. (16) written in finite differences as follows:

$$\begin{array}{c}\hfill \frac{{(\mathrm{\Delta }{p}_0)}_{i,j}}{\mathrm{\Delta }y}=-\frac{1}{\mathrm{\Delta }y}{{\displaystyle \int }}_{y_j-\frac{\mathrm{\Delta }y}{2}}^{y_j+\frac{\mathrm{\Delta }y}{2}}\rho g\mathrm{d}y.\end{array}$$

We used the notation Δp₀ because this relation makes reference to the hydrostatic pressure variation, then using Eq. (29) we write

$$\begin{array}{c}\hfill -\frac{1}{\mathrm{\Delta }y}{{\displaystyle \int }}_{y_j-\frac{\mathrm{\Delta }y}{2}}^{y_j+\frac{\mathrm{\Delta }y}{2}}\rho g\mathrm{d}y=\frac{{p_0}_{i,j}(y_j+\frac{\mathrm{\Delta }y}{2})-{p_0}_{i,j}(y_j-\frac{\mathrm{\Delta }y}{2})}{\mathrm{\Delta }y},\end{array}$$(31)

where ${p_0}_{i,j}({\mathit{y}}_j\pm \frac{\mathrm{\Delta }\mathit{y}}{2})$ are the hydrostatically extrapolated pressures at each cell face (i.e., using Eq. (29)). The gravitational source term vector then results to be

$$\begin{array}{c}\hfill {\mathit{S}}_{i,j}=\left[\begin{array}{c}0\\ 0\\ \frac{1}{\mathrm{\Delta }y}[{p_0}_{i,j}(y_j+\frac{\mathrm{\Delta }y}{2})-{p_0}_{i,j}(y_j-\frac{\mathrm{\Delta }y}{2})]\\ \frac{{(\rho v)}_{i,j}}{{\rho }_{i,j}}\frac{1}{\mathrm{\Delta }y}[{p_0}_{i.j}(y_j+\frac{\mathrm{\Delta }y}{2})-{p_0}_{i,j}(y_j-\frac{\mathrm{\Delta }y}{2})]\end{array}\right]·\end{array}$$(32)

The local hydrostatic reconstruction schemes proposed by others authors present different options in the modeling of gravity and temperature, which defines the conditions under which the scheme is well-balanced. In the works by Zingale et al. (2002) and Fuchs et al. (2010) a constant temperature is assumed a inside the cell with constant gravity; in the work by Käppeli & Mishra (2014) an isentropic condition is assumed with a formulation that allows for modeling self-gravity or non-explicit gravitational terms; in other work by the same authors (Käppeli & Mishra 2016) the well-balanced scheme does not assume thermal equilibrium, but it is needed to solve a system to obtain the discrete equilibrium values to set the initial condition. In this work we assume a linear temperature distribution at each cell to construct a second order scheme with local hydrostatic reconstruction, which is well-balanced for hydrostatic conditions with constant (isothermal) or linear temperature distribution and with constant or variable gravity. In addition, this formulation improve the results for hydrostatic conditions with arbitrary temperature distributions thanks to the second order reconstruction proposed for temperature.

Taking into account that the scheme presented in this work does not require the solution of additional equations nor the definition of new variables, it results to be simpler than the scheme proposed by Käppeli & Mishra (2014, 2016). On the other hand, unlike what happens with the logarithmic reconstruction used by Fuchs et al. (2010), our formulation allows us to satisfy the well-balanced property for more general configurations than the isothermal hydrostatic state with constant gravity. In addition, we consider that this scheme is relatively easy to implement in existing codes with conventional MUSCL schemes, since it only requires us to intervene in the calculation of the pressure Riemann states and in the evaluation of the gravitational source term, independent of the temporal integration scheme.

4.1. Temperature distribution

As we pointed out before, this formulation is based on the assumption that the ideal gas equation of state is valid. This is a key issue of the method, since it allows us to obtain explicit expressions. Although other models of equations of state could be used, this would lead to implicit equations that should be computed by numerical integration, which would impact the efficiency of the scheme. For an ideal gas, we can use the averaged pressure and density values corresponding to each cell to found the product RT needed to evaluate Eq. (29), since

$$\begin{array}{c}\hfill R{T}_{i,j}=\frac{p_{i,j}}{{\rho }_{i,j}}·\end{array}$$

We assume a linear temperature distribution inside each cell i, j, whose cell-centered value matches the actual value, then

$$\begin{array}{c}\hfill {T_0}_{i,j}(y)=T_{i,j}+\frac{{\mathrm{\Delta }T}_{i,j}}{\mathrm{\Delta }y}(y-y_j),y_j-\frac{\mathrm{\Delta }y}{2}\le y\le y_j+\frac{\mathrm{\Delta }y}{2}·\end{array}$$(33)

The slope ΔT_i, j/Δy is estimated by means of a conventional linear reconstruction and a slope limiter function (FL), i.e., using the temperature values at neighboring cells and using a limiter to avoid the appearance of new extrema. Then, the variation ΔT_i, j is obtained by

$$\begin{array}{c}\hfill {\mathrm{\Delta }T}_{i,j}=\mathrm{FL}({\mathrm{\Delta }T}_{i,j}^{-},{\mathrm{\Delta }T}_{i,j}^{+}),\end{array}$$(34)

where

$$\begin{array}{c}\hfill {\mathrm{\Delta }T}_{i,j}^{-}=T_{i,j}-T_{i,j-1},{\mathrm{\Delta }T}_{i,j}^{+}=T_{i,j+1}-T_{i,j}.\end{array}$$

Throughout this work we use the min-mod limiter function, which is defined as follows:

$$\begin{array}{cc}& \mathrm{m}in\text{-}mod({\mathrm{\Delta }T}^{-},{\mathrm{\Delta }T}^{+})=\hfill \\ \hfill & \{\begin{array}{cc}0,\hfill & \mathrm{if}{\mathrm{\Delta }T}^{-}{\mathrm{\Delta }T}^{+}\le 0,\hfill \\ {\mathrm{\Delta }T}^{-},\hfill & \mathrm{if}\left|{\mathrm{\Delta }T}^{-}\right|\le \left|{\mathrm{\Delta }T}^{+}\right|\mathrm{and}{\mathrm{\Delta }T}^{-}{\mathrm{\Delta }T}^{+}>0,\hfill \\ {\mathrm{\Delta }T}^{+},\hfill & \mathrm{if}\left|{\mathrm{\Delta }T}^{+}\right|<\left|{\mathrm{\Delta }T}^{-}\right|\mathrm{and}{\mathrm{\Delta }T}^{-}{\mathrm{\Delta }T}^{+}>0.\hfill \end{array}\hfill \end{array}$$

Naturally, the case of ΔT_i, j = 0 means that temperature is constant inside the cell and for ΔT_i, j ≠ 0 we have a linear temperature distribution. Therefore, we find the following solutions for Eq. (29) depending on the gravity model:

• Constant gravity and constant temperature

  $$\begin{array}{c}\hfill {p_0}_{i,j}(y)=p_{i,j}\exp [-\frac{g}{R{T}_{i,j}}(y-y_j)]·\end{array}$$(35)

• Variable gravity and constant temperature

  $$\begin{array}{c}\hfill {p_0}_{i,j}(y)=p_{i,j}\exp [-\frac{\mathit{GM}}{R{T}_{i,j}}\frac{y-y_j}{(R_0+y)(R_0+y_j)}]·\end{array}$$(36)

• Constant gravity and linear distribution temperature

  $$\begin{array}{c}\hfill {p_0}_{i,j}(y)=p_{i,j}{[1+\frac{{\mathrm{\Delta }T}_{i,j}}{T_{i,j}\mathrm{\Delta }y}(y-y_j)]}^{-g/({R\mathrm{\Delta }T}_{i,j}/\mathrm{\Delta }y)}.\end{array}$$(37)

• Variable gravity and linear distribution temperature

  $$\begin{array}{cc}\hfill {p_0}_{i,j}(y)=& p_{i,j}\exp \{\frac{-GM}{R[T_{i,j}-(R_0+y_j){\mathrm{\Delta }T}_{i,j}/\mathrm{\Delta }y]}\hfill \\ \hfill & \times [\frac{y-y_j}{(R_0+y)(R_0+y_j)}-\frac{{\mathrm{\Delta }T}_{i,j}}{\mathrm{\Delta }y}\frac{ln\left[\frac{T_{i,j}(R_0+y)}{[T_{i,j}+\frac{{\mathrm{\Delta }T}_{i,j}}{\mathrm{\Delta }y}(y-y_j)](R_0+y_j)}\right]}{T_{i,j}-\frac{{\mathrm{\Delta }T}_{i,j}}{\mathrm{\Delta }y}(R_0+y_j)}]\}·\hfill \end{array}$$(38)

Having imposed the temperature variation and knowing the hydrostatic pressure distribution, the hydrostatic density profile inside the cell is determined by means of the equation of state as follows:

$$\begin{array}{c}\hfill {{\rho }_0}_{i,j}(y)=\frac{{p_0}_{i,j}(y)}{R[T_{i,j}+\frac{{\mathrm{\Delta }T}_{i,j}}{\mathrm{\Delta }y}(y-y_j)]}·\end{array}$$(39)

With these equations we can calculate the Riemann states in the y-direction for both faces on each cell together with the corresponding gravitational source term. The final values p_i, j, S, p_i, j, N, ρ_i, j, S, ρ_i, j, N are obtained as functions of the approximation order required for the simulation.

4.2. First order approximation

In the standard first order Godunov scheme the extrapolation of the cell-centered values to the faces is simply made assuming a constant distribution of the variables inside each cell without considering what happens in neighboring cells. In the same way, for the first order approximation of the local hydrostatic reconstruction scheme only the values of the current cell are considered and the hydrostatic pressure and density values at face positions directly constitute the Riemann states. Then, for the normal to the gravity direction faces (S and N) we have

$$\begin{array}{cc}& \left[\begin{array}{c}{\rho }_{i,j,S}\\ u_{i,j,S}\\ v_{i,j,S}\\ p_{i,j,S}\end{array}\right]=\left[\begin{array}{c}{{\rho }_0}_{i,j}(y_j-\frac{\mathrm{\Delta }y}{2})\\ u_{i,j}\\ v_{i,j}\\ {p_0}_{i,j}(y-\frac{\mathrm{\Delta }y}{2})\end{array}\right]={{\mathit{V}}_0}_{i,j,S},\hfill \\ \hfill & \left[\begin{array}{c}{\rho }_{i,j,N}\\ u_{i,j,N}\\ v_{i,j,N}\\ p_{i,j,N}\end{array}\right]=\left[\begin{array}{c}{{\rho }_0}_{i,j}(y_j+\frac{\mathrm{\Delta }y}{2})\\ u_{i,j}\\ v_{i,j}\\ {p_0}_{i,j}(y+\frac{\mathrm{\Delta }y}{2})\end{array}\right]={{\mathit{V}}_0}_{i,j,N},\hfill \end{array}$$(40)

where we used the notation V_0i, j, X (with X = S, N) to indicate that the primitive variables vector corresponds to hydrostatically extrapolated pressure and density.

The source term is calculated using Eq. (32). With these definitions we obtain a well-balanced scheme for a hydrostatic equilibrium condition with constant or variable gravity and constant or linearly variable temperature. We should note that, although this scheme numerically satisfies a hydrostatic equilibrium with the mentioned conditions, it is only first order in space when the variables deviate from a hydrostatic distribution.

4.3. Second order approximation

In the MUSCL-type reconstruction schemes the second order precision in space is reached using a linear function to approximate the variables inside each cell, i.e., the distribution of a generic variable ϕ in the y-direction inside the i, j cell is assumed to be

$$\begin{array}{c}\hfill {\overline{\varphi }}_{i,j}(y)={\varphi }_{i,j}+\frac{{\mathrm{\Delta }\varphi }_{i,j}}{\mathrm{\Delta }y}(y-y_j),y-\frac{\mathrm{\Delta }y}{2}\le y\le y_j+\frac{\mathrm{\Delta }y}{2},\end{array}$$(41)

where Δϕ_i, j is a “limited” variation obtained by some limiter function to avoid the appearance of new extrema and then preserve the solution monotonicity.

Clearly, a standard linear reconstruction like this for pressure leads to a conventional scheme that does not preserve a hydrostatic equilibrium condition. However, this concept can be combined with the local hydrostatic reconstruction to obtain a second order precision scheme in space, as proposed by Käppeli & Mishra (2014, 2016). This approach consists in splitting the pressure in a hydrostatic term p₀ and in a equilibrium deviation term p′, which can be large such that in the i, j cell we have

$$\begin{array}{c}\hfill p_{i,j}(y)={p_0}_{i,j}(y)+p_{i,j}^{\prime }(y).\end{array}$$(42)

By definition we set p_0i, j(y_j) = p_i, j (the cell-centered value) and p′_i, j(y_j) = 0, since the hydrostatic profile is calculated choosing p_i, j at y = y_j as the reference pressure.

The second order local hydrostatic reconstruction is reached using the equilibrium deviations ${p^{\prime }}_{i,j}^{\mp }(y\mp \Delta y)$ in the neighboring cells i, j ∓ 1 to make a linear approximation of the local equilibrium deviation inside the cell. The equilibrium deviation in the neighboring cells are measured with respect to the local hydrostatic equilibrium in the current cell and are calculated by comparing the actual pressure values p_{i, j ∓ 1} at i, j ∓ 1 cells with the local hydrostatic profile p_0i, j(y) projected to the centroids of coordinates (x_i, y_j ∓ 1). From Eq. (42) we write

$$\begin{array}{cc}& p_{i,j}^{\prime -}=p_{i,j-1}-{p_0}_{i,j}(y_{j-1}),\hfill \\ \hfill & p_{i,j}^{\prime +}=p_{i,j+1}-{p_0}_{i,j}(y_{j+1}).\hfill \end{array}$$(43)

Taking into account that, by definition, p′_i, j(y_j) = 0, the local pressure equilibrium deviations at each neighboring cell are written as

$$\begin{array}{cc}& {\mathrm{\Delta }{p}^{\prime }}_{i,j}^{-}=-p_{i,j}^{\prime -},\hfill \\ \hfill & {\mathrm{\Delta }{p}^{\prime }}_{i,j}^{+}=p_{i,j}^{\prime +}.\hfill \end{array}$$(44)

With these values and a limiter function we obtain the variation to define the slope of the linear approximation for the deviation

$$\begin{array}{c}\hfill p_{i,j}^{\prime }(y)=\frac{{\mathrm{\Delta }{p}^{\prime }}_{i,j}}{\mathrm{\Delta }y}(y-y_j),y-\frac{\mathrm{\Delta }y}{2}\le y\le y+\frac{\mathrm{\Delta }y}{2},\end{array}$$

where

$$\begin{array}{c}\hfill {\mathrm{\Delta }{p}^{\prime }}_{i,j}=\mathrm{FL}\{{\mathrm{\Delta }{p}^{\prime }}_{i,j}^{-},{\mathrm{\Delta }{p}^{\prime }}_{i,j}^{+}\}.\end{array}$$(45)

Finally, the reconstructed pressure inside the cell is the sum of the local hydrostatic profile p_0i, j(y) and the linear approximation of the deviation equilibrium pressure p′_i, j(y), i.e.,

$$\begin{array}{c}\hfill {\overline{p}}_{i,j}(y)={p_0}_{i,j}(y)+\frac{{\mathrm{\Delta }{p}^{\prime }}_{i,j}}{\mathrm{\Delta }y}(y-y_j),y-\frac{\mathrm{\Delta }y}{2}\le y\le y+\frac{\mathrm{\Delta }y}{2}·\end{array}$$(46)

We emphasize that for a global hydrostatic condition with constant or linearly variable temperature, the deviation pressure is zero and the hydrostatically reconstructed pressure results to be equal to the global distribution.

The density reconstruction is made in a similar way, using in this case the density equilibrium deviations ${{\rho }^{\prime }}_{i,j}^{\mp }$, which are defined analogously to Eq. (44). The reconstructed density inside the cell is written as

$$\begin{array}{c}\hfill {\overline{\rho }}_{i,j}(y)={{\rho }_0}_{i,j}(y)+\frac{{\mathrm{\Delta }{\rho }^{\prime }}_{i,j}}{\mathrm{\Delta }y}(y-y_j),y-\frac{\mathrm{\Delta }y}{2}\le y\le y+\frac{\mathrm{\Delta }y}{2},\end{array}$$(47)

where, as before, the variation Δρ′_i, j is obtained with the limiter function

$$\begin{array}{c}\hfill {\mathrm{\Delta }{\rho }^{\prime }}_{i,j}=\mathrm{FL}\{{\mathrm{\Delta }{\rho }^{\prime }}_{i,j}^{-},{\mathrm{\Delta }{\rho }^{\prime }}_{i,j}^{+}\},\end{array}$$(48)

where

$$\begin{array}{cc}& {\mathrm{\Delta }{\rho }^{\prime }}_{i,j}^{-}=-{\rho }_{i,j}^{\prime -}={\rho }_{i,j-1}-{{\rho }_0}_{i,j}(y_{j-1}),\hfill \\ \hfill & {\mathrm{\Delta }{\rho }^{\prime }}_{i,j}^{+}={\rho }_{i,j}^{\prime +}={\rho }_{i,j+1}-{{\rho }_0}_{i,j}(y_{j+1}).\hfill \end{array}$$(49)

The pressure and density Riemann states for each cell face normal to the gravity direction are obtained evaluating Eqs. (46) and (47) at face positions. The velocity components for the RPs are calculated using the standard linear reconstruction. Then, we have for the south and north faces of the i, j cell,

$$\begin{array}{cc}& \left[\begin{array}{c}{\rho }_{i,j,S}\\ u_{i,j,S}\\ v_{i,j,S}\\ p_{i,j,S}\end{array}\right]=\left[\begin{array}{c}{{\rho }_0}_{i,j}(y_{j-\frac{1}{2}})-\frac{1}{2}{\mathrm{\Delta }{\rho }^{\prime }}_{i,j}\\ u_{i,j}-\frac{1}{2}{\mathrm{\Delta }u}_{i,j}\\ v_{i,j}-\frac{1}{2}{\mathrm{\Delta }v}_{i,j}\\ {p_0}_{i,j}(y_{j-\frac{1}{2}})-\frac{1}{2}{\mathrm{\Delta }{p}^{\prime }}_{i,j}\end{array}\right]={{\mathit{V}}_0}_{i,j,S}-\frac{1}{2}{{\mathbf{\Delta }\mathit{V}}^{\prime }}_{i,j},\hfill \\ & \left[\begin{array}{c}{\rho }_{i,j,N}\\ u_{i,j,N}\\ v_{i,j,N}\\ p_{i,j,N}\end{array}\right]=\left[\begin{array}{c}{{\rho }_0}_{i,j}(y_{j-\frac{1}{2}})+\frac{1}{2}{\mathrm{\Delta }{\rho }^{\prime }}_{i,j}\\ u_{i,j}+\frac{1}{2}{\mathrm{\Delta }u}_{i,j}\\ v_{i,j}+\frac{1}{2}{\mathrm{\Delta }v}_{i,j}\\ {p_0}_{i,j}(y_{j-\frac{1}{2}})+\frac{1}{2}{\mathrm{\Delta }{p}^{\prime }}_{i,j}\end{array}\right]={{\mathit{V}}_0}_{i,j,N}+\frac{1}{2}{{\mathbf{\Delta }\mathit{V}}^{\prime }}_{i,j},\hfill \end{array}$$(50)

where the variations Δu_i, j and Δv_i, j are limited values.

As in the first order scheme, the evaluation of the gravitational source term is made using Eq. (32). With these definitions, the second order scheme is well-balanced for a hydrostatic equilibrium condition with constant or variable gravity and constant or linearly variable temperature, since under this condition the equilibrium deviations values p′ and ρ′ are zero, leading to constant states RPs at each intercell.

We point out that this second order well-balanced scheme differs from a conventional MUSCL scheme in the form of reconstructing pressure and density and evaluating the gravitational source term. In this way, the implementation of the local hydrostatic reconstruction scheme in an existing code with MUSCL-type reconstruction is relatively simple. The additional computational cost of this implementation lies in the calculation of integral (29), which in the case of the equation of state of ideal gases can be explicitly performed.

4.4. Boundary conditions

As known, the boundary conditions treatment is a very important aspect for any gas dynamic numerical simulation, but for the specific case of hydrostatic equilibrium problems this aspect becomes crucial. If the boundary conditions are not correctly modeled in a hydrostatic problem, spurious mass fluxes can enter or exit the domain in the gravity direction, with velocities that increase in time and can reach very high values, even if we are using a well-balanced scheme. Naturally, these nonphysical fluxes can strongly affect the solution.

The key for modeling boundary conditions in contours normal to the gravity direction is to consider the gravitational acceleration in the outside extrapolation of pressure and density. For an outflow (zero gradient for velocity), reflecting (or symmetry) or any Neumann boundary condition applied on velocity parallel to the gravity direction, the pressure extrapolation has to be made considering a hydrostatic profile where the reference pressure is the cell-centered pressure in the last cell adjacent to the boundary, i.e.,

$$\begin{array}{ccc}& p_{i,-J}=p_{i,1}\exp (-{{\displaystyle \int }}_{y_1}^{y_{-J}}\frac{g(y^{\prime })}{RT(y^{\prime })}\mathrm{d}{y}^{\prime })\hfill & \hfill \text{(bottom boundary)},\\ \hfill & p_{i,M+J}=p_{i,M}\exp (-{{\displaystyle \int }}_{y_M}^{y_{M+J}}\frac{g(y^{\prime })}{RT(y^{\prime })}\mathrm{d}{y}^{\prime })\hfill & \hfill \text{(top boundary)},\end{array}$$(51)

where J = 1, …, n_g, where n_g is the number of ghost cells. The subindices 1 and M correspond to the first and the last cells of the physic domain in the gravity direction.

To evaluate (51) a temperature distribution has to be imposed, which can be made in different ways. The simplest temperature distribution is to assume a constant temperature in the ghost cells equal to the temperature in the last interior cell, with which Eq. (51) result to be equal to Eqs. (35) or (36) depending on the gravity model. The assumption of constant temperature works very well when the gradient temperature near the boundary is not too large. In case of the temperature variation in the last interior cells is large, it can be assumed to be a linear profile for temperature in the ghost cells, whose slope is calculated with the values at the last cells, i.e.,

$$\begin{array}{cc}& T(y)=T_{i,1}+\frac{\mathrm{\Delta }{T}_g^{-}}{\mathrm{\Delta }y}(y_1-y)\text{(bottom boundary)},\hfill \\ \hfill & T(y)=T_{i,M}+\frac{\mathrm{\Delta }{T}_g^{+}}{\mathrm{\Delta }y}(y_M-y)\text{(top boundary)}.\hfill \end{array}$$(52)

With this assumption, Eq. (51) take the form of Eqs. (37) or (38) according to the gravity model.

After calculating the pressure at the ghost cells, density is obtained using the equation of state. The remaining variables (velocity and magnetic field components) are computed considering the particular boundary condition for each of them.

5. Implementation in the FLASH code

The FLASH cCode is a publicly available high performance multiphysics open source software that can be used in multiple astrophysical problems (Fryxell et al. 2000). For the MHD simulation unit, FLASH uses Godunov’s method based on numerical schemes with adaptive mesh refinement (AMR) capabilities and allows the treatment of high energy compressible plasma flows. The version 4.5 of this code has two available numerical schemes: one with a directional split formulation based on the Powell method (Powell et al. 1999) and another with an unsplit formulation based on the scheme proposed by Lee et al. (2009). This second scheme, called the USM method for unsplit staggered mesh, presents a more consistent formulation for magnetic field treatment, since it uses a combination of the constrained transport (CT) method and the corner transport upwind (CTU) method to evolve the magnetic field in a staggered mesh, avoiding the appearance of numerical magnetic field divergence (Gardiner & Stone 2005, 2008).

In this section we detail the implementation of the local hydrostatic reconstruction scheme proposed in the previous section in the MHD USM unit of the FLASH Code. For the second order approximation, the USM algorithm uses the MUSCL-Hancock scheme with characteristic limiting option in the prediction step. For this implementation we only considered the gas dynamic terms, since the magnetic field terms are not affected by the local hydrostatic reconstruction. We accentuate that, taking into account that the coupling of the magnetic field in the MHD equations is only given through velocity and that the hydrostatic reconstruction is only applied on density and pressure with the velocity reconstructed in the standard form, it turns out that the magnetic field is not affected by the hydrostatic scheme.

The basic second order reconstruction-evolution scheme of the USM method is obtained starting from the quasi-linear form of the primitive MHD equations (Lee & Deane 2009)

$$\begin{array}{cc}\hfill {\mathit{V}}_{i,j,W,E}^{n+\frac{1}{2}}=& {\mathit{V}}_{i,j}^n\mp \frac{1}{2}{\mathbf{\Delta }\mathit{V}}_{x,i,j}-\frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}{\mathbf{A}}_x^n\left({\mathit{V}}_{i,j}^n\right){\mathbf{\Delta }\mathit{V}}_{x,i,j}\hfill \\ \hfill & -\frac{\mathrm{\Delta }t}{2\mathrm{\Delta }y}{\mathbf{A}}_y\left({\mathit{V}}_{i,j}^n\right){\mathbf{\Delta }\mathit{V}}_{y,i,j},\hfill \\ \hfill {\mathit{V}}_{i,j,S,N}^{n+\frac{1}{2}}=& {\mathit{V}}_{i,j}^n\mp \frac{1}{2}{\mathbf{\Delta }\mathit{V}}_{y,i,j}-\frac{\mathrm{\Delta }t}{2\mathrm{\Delta }y}{\mathbf{A}}_y\left({\mathit{V}}_{i,j}^n\right){\mathbf{\Delta }\mathit{V}}_{y,i,j}\hfill \\ \hfill & -\frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}{\mathbf{A}}_x\left({\mathit{V}}_{i,j}^n\right){\mathbf{\Delta }\mathit{V}}_{x,i,j},\hfill \end{array}$$(53)

where the plus and minus signs correspond to the W, E (for x-direction) and S, N (for y-direction), respectively. The last term in each equation represents the transverse flux contribution that arises from the unsplit formulation and that plays a crucial role in the numerical stability of the CTU method (Bell et al. 1988). We should note that the reconstruction-evolution step and the computation of the source term are the only parts of the numerical algorithm that need to be modified to obtain the well-balanced property for hydrostatic equilibrium. We emphasize that the slope limiters, Riemann solvers, numerical fluxes computation, and the advance of the variables to the next time step remain in their original form, which allows us to keep all the FLASH options for the simulations.

Firstly, we treated the parallel to the gravity direction reconstruction variables, which is indicated by ${\widehat{\mathit{V}}}_{\alpha ,i,j}^{n+\frac{1}{2}}$ and is equal to the reconstruction-evolution equation for S, N removing the last term. Because the gravity acceleration is parallel to the y-direction, the reconstruction in the x-direction (W and E faces) is not affected with respect to the original scheme, except for the transverse flux contribution term, which is treated below.

For the characteristic limiting of the reconstruction, the relationships that link the characteristic variables vector W with the primitive variables V are

$$\begin{array}{c}\hfill \mathit{V}={\mathit{R}}_y{\mathit{W}}_y,{\mathit{L}}_y{\mathbf{A}}_y{\mathit{R}}_y={\mathbf{\Lambda }}_y,\end{array}$$(54)

where R_y and ${\mathit{L}}_{\mathit{y}}={\mathit{R}}_{\mathit{y}}^{-1}$ are the right and left eigenvectors matrices in y-direction, respectively, and Λ_y is the diagonal eigenvalues matrix of the Jacobian matrix A_y.

With these relationships, and taking into account that we only considered the normal to the face reconstruction, the reconstruction-evolution equation in the gravity direction reads

$$\begin{array}{c}\hfill {\widehat{\mathit{V}}}_{i,j,S,N}^{n+\frac{1}{2}}={\mathit{V}}_{i,j}^n\mp \frac{1}{2}{\mathit{R}}_{y,i,j}^n{\mathbf{\Delta }\mathit{W}}_{y,i,j}-\frac{\mathrm{\Delta }t}{2\mathrm{\Delta }y}{\mathit{R}}_{y,i,j}^n{\mathbf{\Lambda }}_{y,i,j}^n{\mathbf{\Delta }\mathit{W}}_{y,i,j},\end{array}$$(55)

where ΔW_y, i, j is the characteristic limited vector given by

$$\begin{array}{c}\hfill {\mathbf{\Delta }\mathit{W}}_{y,i,j}=\mathrm{FL}\{{\mathit{L}}_{y,i,j}^n·{\mathbf{\Delta }\mathit{V}}_{y,i,j}^{-},{\mathit{L}}_{y,i,j}^n·{\mathbf{\Delta }\mathit{V}}_{y,i,j}^{+}\}.\end{array}$$(56)

For the hydrostatic equilibrium conditions considered in our well-balanced scheme, we expect that ΔW_y, i, j = 0, the pressure and density values projected to a face from adjacent cells are equal and the gravitational source term cancels the pressure gradient inside the cell, so that the hydrostatic equilibrium is numerically preserved. Therefore, we needed to reformulate the reconstruction-evolution Eq. (55) to ensure the USM algorithm of FLASH satisfies the well-balanced property. For this purpose, we used the hydrostatically extrapolated values proposed in the first order approximation (40) to set the face values when ΔW_y, i, j = 0. In addition, we defined a new limiting of the characteristic variables to avoid taking into account the changes due to the hydrostatic stratification. This was made using Eqs. (44) and (49) for the variations ${\mathrm{\Delta }p}_{\mathit{y},i,j}^{\mp }$ and ${\mathrm{\Delta }\rho }_{\mathit{y},i,j}^{\mp }$ in Eq. (56). With these definitions, Eq. (55) becomes

$$\begin{array}{c}\hfill {\widehat{\mathit{V}}}_{i,j,S,N}^{n+\frac{1}{2}}={{\mathit{V}}_0}_{i,j,S,N}^n\mp \frac{1}{2}{\mathit{R}}_{y,i,j}^n{\mathbf{\Delta }\mathit{W}}_{y,i,j}^{\prime }-\frac{\mathrm{\Delta }t}{2\mathrm{\Delta }y}{\mathit{R}}_{y,i,j}^n{\mathbf{\Lambda }}_{y,i,j}^n{\mathbf{\Delta }\mathit{W}}_{y,i,j}^{\prime },\end{array}$$(57)