High-order well-balanced finite volume schemes for the Euler equations with gravitation
HHigh-order well-balanced finite volume schemes for the Euler equationswith gravitation
L. Grosheintz-Laval a, ∗ , R. K¨appeli a a Seminar for Applied Mathematics (SAM), Department of Mathematics, ETH Z¨urich, CH-8092 Z¨urich, Switzerland
Abstract
A high-order well-balanced scheme for the Euler equations with gravitation is presented. The scheme is ableto preserve a spatially high-order accurate discrete representation of a large class of hydrostatic equilibria.It is based on a novel local hydrostatic reconstruction, which, in combination with any standard high-orderaccurate reconstruction procedure, achieves genuine high-order accuracy for smooth solutions close or awayfrom equilibrium. The resulting scheme is very simple and can be implemented into any existing finitevolume code with minimal e ff ort. Moreover, the scheme is not tied to any particular form of the equation ofstate, which is crucial for example in astrophysical applications. Several numerical experiments demonstratethe robustness and high-order accuracy of the scheme nearby and out of hydrostatic equilibrium. Keywords:
Numerical methods, Hydrodynamics, Source terms, Well-balanced schemes
1. Introduction
A multitude of interesting physical phenomena are modeled by the Euler equations with gravitationalsource terms. Applications range from the study of atmospheric phenomena, such as numerical weatherprediction and climate modeling, to the numerical simulation of the climate of exoplanets, convection instars and core-collapse supernova explosions. The Euler equations with gravitational source terms expressthe conservation of mass, momentum and energy: ∂ρ∂ t + ∇ · ( ρ v ) = ∂ρ v ∂ t + ∇ · ( ρ v ⊗ v ) + ∇ p = − ρ ∇ φ (1.2) ∂ E ∂ t + ∇ · ( v ( E + p )) = − ρ v ∇ φ. (1.3)Here ρ is the mass density, v the velocity and E = ρ e + ρ v (1.4) ∗ Corresponding author
Email address: [email protected] (L. Grosheintz-Laval )
Preprint submitted to Elsevier July 12, 2018 a r X i v : . [ m a t h . NA ] J u l he total fluid energy density being the sum of internal and kinetic energy densities. The pressure p is relatedto the density and specific internal energy through an equation of state p = p ( ρ, e ).The source terms on the right-hand side of the momentum and energy equations model the e ff ect of thegravitational forces on the fluid. They are dictated by the variation of the gravitational potential φ , whichcan either be a given function or, in the case of self-gravity, be determined by the Poisson equation ∇ φ = π G ρ, (1.5)where G is the gravitational constant.In many physically relevant applications, such as the ones named above, (parts of) the flow of interestmay be realized close to hydrostatic equilibrium ∇ p = − ρ ∇ φ. (1.6)As a matter of fact, the numerical simulation of near equilibrium flows is challenging for standard finitevolume methods. The reason for this is that these methods may in general not satisfy a discrete equivalent ofthe equilibrium. Thus such states are not preserved exactly but are solely approximated with an error propor-tional to the truncation error of the scheme. So if the interest relies in the simulation of small perturbationson top of a hydrostatic equilibrium, the numerical resolution has to be increased to the point that the trun-cation errors do not obscure these small perturbations. This may result in prohibitively high computationalcosts, especially in several space dimensions.A design principle to overcome the challenge was introduced by Greenberg and Leroux [1] leading tothe concept of so-called well-balanced schemes. In these schemes, a discrete equivalent of the equilibrium isexactly satisfied. Therefore, they possess the ability to maintain discrete equilibrium states down to machineprecision and are capable of resolving small equilibrium perturbations e ff ectively. Many well-balancedschemes have been designed, especially for the shallow water equations with non-trivial bottom topography,see e.g. [2, 3, 4] and references therein. An extensive review on well-balanced schemes for many di ff erentapplications is also given in the book by Gosse [5].Well-balanced schemes for the Euler equations with gravitation have received a considerable amount ofattention in the recent literature. First, LeVeque and Bale [6] have applied the quasi-steady wave-propagationalgorithm [2] to the Euler equations with gravity. Few years later, Botta et al. [7] designed a well-balancedfinite volume scheme for numerical weather prediction applications. More recently, several well-balancedfinite volume [8, 9, 10, 11, 12, 13, 14, 15, 16] , finite di ff erence [17, 18] and discontinuous Galerkin [19, 20,21] schemes have been presented. Magnetohydrostatic steady state preserving well-balanced finite volumeschemes were devised in [22]. To the best of our knowledge, many of the mentioned schemes are at mostsecond-order accurate and only [17, 19, 20, 18, 21] go to higher orders. However, with the notable exceptionof [20], it appears that these schemes need the equilibrium to be predetermined.In fact, equation (1.6) only specifies a mechanical equilibrium. In order to fully characterize the equilib-rium a thermal variable, such as the specific entropy s or the temperature T , needs to be supplemented. As aconcrete astrophysically relevant example of a stationary state we consider the case of constant entropy. Therelevant thermodynamic relation for isentropic hydrostatic equilibrium isd h = T d s + d p ρ , (1.7)where h is the specific enthalpy h = e + p ρ , (1.8)2 the temperature and s the specific entropy. Then we can write (1.6) for the isentropic case (d s =
0) as1 ρ ∇ p = ∇ h = −∇ φ. (1.9)The last equation can then be trivially integrated to obtain h + φ = const . (1.10)In [9] this equilibrium was used to build a second-order accurate well-balanced finite volume scheme. Alongthe same lines, well-balanced schemes for isothermal hydrostatic equilibrium can be constructed [15]. In thelatter case, the relevant thermodynamic potential is the Gibbs free energy.In this paper, we extend the well-balanced finite volume schemes [9] beyond second-order accuracy. Thescheme possesses the following novel features: • An arbitrarily high-order accurate local hydrostatic profile is constructed based on the equilibrium(1.10). • An arbitrarly high-order equilibrium preserving reconstruction is designed on the basis of any standardhigh-order reconstruction procedure. • A well-balanced source term discretization is built from the equilibrium preserving reconstruction. • It is well-balanced for any consistent numerical flux, which allows a straightforward implementationwithin any standard finite volume method. • It is well-balanced for multi-dimensional hydrostatic equilibria. • It is not tied to any particular equation of state such as the ideal gas law. This is important, especiallyfor astrophysical applications.The rest of the paper is structured as follows: the well-balanced finite volume scheme is presented insection 2. Extensive numerical results are presented in section 3 and conclusions are provided in section 4.
2. Numerical Method
We first consider the Euler equations with gravitation (1.1–3) in one space dimension and write them inthe following compact form ∂ u ∂ t + ∂ f ∂ x = s (2.1)with u = ρρ v x E , f = ρ v x ρ v x + p ( E + p ) v x and s = − ρρ v x ∂φ∂ x , (2.2)where u , f and s are the vectors of conserved variables, fluxes and source terms. An equation of state (EoS) p = p ( ρ, e ) relates the pressure to the density ρ and specific internal energy e (or any other thermodynamicquantity such as specific entropy s or temperature T ). For example, a simple EoS is provided by the idealgas law p = ρ e ( γ − , (2.3)3here γ is the ratio of specific heats. We stress that the well-balanced scheme derived below is not tied toany particular form of EoS, which is crucial especially in astrophysical applications.In the next section we will briefly describe a standard high-order finite-volume discretization and it’s corecomponents in order to fix the notation. The following sections will then describe our novel well-balancedscheme in detail. For the numerical approximation of (2.1), the spatial domain of interest is discretized by a number ofcells or finite volumes I i = [ x i − / , x i + / ]. Here x i ± / denotes the left and right cell interface, respectively,and x i = ( x i − / + x i + / ) / I i . For ease of presentation, we assume a regular cell size ∆ x = x i + / − x i − / . Nevertheless, varying cell sizes can easily be accommodated for.A one-dimensional semi-discrete finite volume scheme is then given byd U i d t = L ( U ) = − ∆ x (cid:0) F i + / − F i − / (cid:1) + S i , (2.4)where U i = U i ( t ) denotes the approximate cell average of the conserved variables in cell I i at time t . Itapproximates the exact cell average u i = u i ( t ) of the true solution u ( t , x ) at time t : U i ( t ) ≈ u i ( t ) = ∆ x (cid:90) I i u ( t , x ) d x . (2.5)In the following, a quantity with an overbar indicates a cell average while a quantity without indicates apoint value. By S i ( t ) is denoted the approximate cell average of the true source terms at time t : S i ( t ) ≈ s i ( t ) = ∆ x (cid:90) I i s ( u , ∂φ∂ x ) d x . (2.6)Note that we have suppressed the time dependence of the gravitational potential since we are mainly con-cerned with flows close to hydrostatic equilibrium and for ease of notation. Numerical flux.
The numerical flux is obtained by solving (approximately) the Riemann problem at cellinterfaces F i + / = F ( U i + / − , U i + / + ) , (2.7)where the point values U i + / ∓ are the cell interface extrapolated conserved variables and F is a consistent,i.e. F ( u , u ) = f ( u ), and Lipschitz continuous numerical flux function.Below, we will make use of the HLLC approximate Riemann solver with simple wave speed estimatesfrom [23, 24]. Though, our well-balanced scheme is independent of this particular choice. Reconstruction.
The purpose of a reconstruction procedure R is to compute accurate point values of theapproximate solution U i ( t , x ) within each cell from the cell averages U . We denote such a reconstructionprocedure, which recovers a r -th order accurate point value of a quantity c at location x within cell I i fromthe cell averages c , by c i ( x ) = R ( x ; { c k } k ∈ S i ) . (2.8)Here S i is the stencil for the reconstruction procedure for cell I i , i.e. S i is a finite set of neighbors of I i .The values of the conserved variables extrapolated to the interface are then given by U i + / − = U i ( t , x i + / ) = R (cid:16) x i + / ; { U k } k ∈ S i (cid:17) and U i + / + = U i + ( t , x i + / ) = R (cid:16) x i + / ; { U k } k ∈ S i + (cid:17) . Source term discretization.
The approximate cell average of the source term S i is obtained by numericalintegration. Let Q i denote a q -th order accurate quadrature rule over cell I i . Then the cell average of thesource term is approximated by S i = ∆ x Q i (cid:32) s ( U , ∂φ∂ x ) (cid:33) = ∆ x N q (cid:88) α = ω α s (cid:32) U i ( t , x i ,α ) , ∂φ∂ x ( x i ,α ) (cid:33) , (2.9)where the x i ,α ∈ I i and ω α denote the N q quadrature nodes and weights of Q i , respectively. For example,the two-point Gauss-Legendre quadrature rule can be used, which is the choice we will make below. Thepoint values of the conserved variables at the quadrature nodes U i ( t , x i ,α ) are obtained by the reconstructionprocedure: U i ( t , x i ,α ) = R (cid:16) x i ,α ; { U k } k ∈ S i (cid:17) . (2.10)If the gravitational potential is known analytically, it can be evaluated directly at the quadrature nodes. If itis not, then a suitable interpolation has to be applied. Temporal discretization.
The temporal domain of interest [0 , T ] is discretized into time steps ∆ t = t n + − t n ,where the superscript n labels the di ff erent time levels. For the temporal integration, the high-order strongstability-preserving Runge-Kutta (SSP-RK) schemes [31] can be used. In particular, we use the third-orderSSP-RK method for the numerical results presented in this paper U (1) i = U ni + ∆ t L ( U n ) U (2) i = U ni + (cid:18) U (1) i + ∆ t L ( U (1) ) (cid:19) U n + i = U ni + (cid:18) U (2) i + ∆ t L ( U (2) ) (cid:19) , (2.11)where L denotes the spatial discretization operator from (2.4). Furthermore, the time step ∆ t has to fulfill acertain CFL condition.This concludes the description of a standard high-order finite volume scheme for the Euler equations. Werefer to the excellent books available in the literature for detailed derivations, e.g. [32, 33, 34, 2]. However,a standard reconstruction procedure and source term discretization will in general not preserve a discreteequivalent of hydrostatic equilibrium. In order to achieve this, we need the ingredients presented in thefollowing two sections 2.1.2 and 2.1.3. The local hydrostatic reconstruction consists of two parts. First, within each cell a high-order accurateequilibrium profile that is consistent with the cell-averaged conserved variables is determined. Second, the5ell’s equilibrium profile is extrapolated to neighboring cells to perform a high-order accurate reconstructionof the equilibrium perturbation.We begin by describing how the local high-order accurate equilibrium profile is determined. Within the i -th cell I i , we define a subcell equilibrium reconstruction of the specific enthalpy h eq , i ( x ) by assuming (1.10)as h eq , i ( x ) = h , i + φ i − φ ( x ) . (2.12)Here h , i = h eq , i ( x i ) and φ i = φ ( x i ) are point values of the specific enthalpy and the gravitational potentialat the cell center, respectively. In the following, we assume that the gravitational potential can be evaluatedanywhere, either because it is a given function or obtained by a suitable interpolation.In combination with the (assumed constant) equilibrium entropy s , i in cell I i , the equilibrium density ρ eq , i ( x ) and internal energy density ρ e eq , i ( x ) profiles can be computed through the EoS: ρ eq , i ( x ) = ρ ( h eq , i ( x ) , s , i ) and ρ e eq , i ( x ) = ρ e ( h eq , i ( x ) , s , i ) . The computational complexity of this computation depends strongly on the functional form of the EoS. Forthe ideal gas case, explicit expressions are given in Appendix A.We note that the equilibrium specific enthalpy h , i and entropy s , i are not specified so far. In order tofix h , i and s , i , we demand that the equilibrium density and internal energy density profiles agree up to thedesired order of accuracy with their respective cell average in cell I i . Hence, we seek h , i and s , i such that ρ i = ∆ x Q i ( ρ eq , i ) = ∆ x N q (cid:88) α = ω α ρ ( h eq , i ( x i ,α ) , s , i ) ρ e i = ∆ x Q i ( ρ e eq , i ) = ∆ x N q (cid:88) α = ω α ρ e ( h eq , i ( x i ,α ) , s , i ) , (2.13)where Q i denotes the previously introduced q -th order accurate quadrature rule over cell I i . In the aboveexpression, an estimate of the cell average of the internal energy density ρ e i is needed. We simply estimateit directly from the cell-averaged conserved variables by ρ e i = E i − ρ v x , i ρ i , (2.14)which is exact at equilibrium ( v x ≡ h , i and the (constant) specific entropy s , i . This system must be solved iteratively, e.g.with Newton’s method. In practice, the iterative process is started from the specific entropy and enthalpycomputed from the cell-averaged conserved variables U i . The cost of this iterative process is mitigated bythe fact that it is local to each cell and the initial guess is a spatially second order accurate estimate, i.e.a very small two-by-two system of equations must be solved, independently, in every cell starting from agood initial guess. For the ideal gas law, the system can be reduced to a single nonlinear equation for whichexistence and uniqueness of the solution can be guaranteed under very weak requirements. This is shown inAppendix A.Once h , i and s , i have been fixed, we have the following high-order accurate representation of the equi-librium in cell I i : U eq , i ( x ) = ρ eq , i ( x )0 ρ e eq , i ( x ) . (2.15)6ext we develop the high-order equilibrium preserving reconstruction procedure. The idea is to decom-pose the solution into an equilibrium and a (possibly large) perturbation part. Within cell I i , the equilibriumpart is simply given by the previously derived equilibrium profile U eq , i ( x ). The perturbation part is obtainedby applying the standard reconstruction R procedure on the equilibrium perturbation cell averages δ U i ( x ) = R (cid:16) x ; { U k − Q k ( U eq , i ) } k ∈ S i (cid:17) , (2.16)which results in a min( q , r )-th order accurate representation of the equilibrium perturbation in cell I i . Notethat the equilibrium perturbation cell average in cell I k is obtained by taking the di ff erence between theactual cell average U k in cell I k and the cell average of the equilibrium profile U eq , i ( x ) in cell I k . The latter isevaluated by applying the I k cell’s quadrature rule Q k to U eq , i ( x ).The full equilibrium preserving reconstruction W is then obtained by simply adding the equilibriumprofile to the perturbation U i ( x ) = W ( x ; { U k } k ∈ S i ) = U eq , i ( x ) + δ U i ( x ) . (2.17)We observe that, by construction, this reconstruction will preserve any equilibrium of the form (1.10), sincethe perturbation δ U i ( x ) vanishes under these conditions. Remark 2.1
Any function can be written as some other function plus the di ff erence. Clearly, this di ff erencecan be reconstructed from the cell-averages of the di ff erence. Therefore, the well-balanced reconstructionprocedure (2.17) is high-order accurate, for any smooth function U eq , i ( x ) . In particular, the choice of anonly second order accurate estimate of ρ e i does not a ff ect the overall order of the reconstruction.2.1.3. Well-balanced source term discretization For the momentum source discretization, we use the previous splitting of the cell I i ’s density ρ i ( x ) intoequilibrium ρ eq , i ( x ) and perturbation δρ i ( x ) as S ρ v , i ( x ) = − ρ i ( x ) ∂φ∂ x ( x ) = − (cid:16) ρ eq , i ( x ) + δρ i ( x ) (cid:17) ∂φ∂ x ( x ) = − ρ eq , i ( x ) ∂φ∂ x ( x ) − δρ i ( x ) ∂φ∂ x ( x ) , which is clearly a pointwise min( q , r )-th order accurate approximation of the true source term. However, astraightforward numerical integration will not result in a well-balanced scheme. Instead, we use the fact thatfor the equilibrium profiles we have ∂ p eq , i ∂ x = − ρ eq , i ∂φ∂ x by construction. As a result, the equilibrium part of the momentum source term can be trivially integratedand numerical integration is only applied to the perturbation part: S ρ v , i = p eq , i ( x i + / ) − p eq , i ( x i − / ) ∆ x − ∆ x Q i (cid:32) δρ i ∂φ∂ x (cid:33) . (2.18)Since we are only concerned with stationary equilibria, the energy source term S E , i discretization is leftunchanged from (2.9).We summarize the developed high-order well-balanced finite volume scheme in the following theorem: Theorem 2.2
Consider the scheme (2.4) with a consistent and Lipschitz continuous numerical flux F , a r-thorder accurate spatial reconstruction procedure R , a q-th order accurate quadrature rule Q , the hydrostaticreconstruction W (2.17) and the gravitational source term S (2.9) (with (2.18) ).This scheme has the following properties: The scheme is consistent with (2.1) and it is min( q , r ) -th order accurate in space (for smooth solutions). (ii) The scheme is well-balanced and preserves the discrete hydrostatic equilibrium given by (1.10) andv x = exactly. Proof (i) The consistency and formal order of accuracy of the scheme is straightforward.(ii) Let the hydrostatic equilibrium (1.10) be characterized by the constant specific entropy s and specificenthalpy profile h eq ( x ). The equilibrium conserved variables are then given by u eq ( x ) = [ ρ ( h eq ( x ) , s ) , , ρ e ( h eq ( x ) , s ] T and let U i (0) = ∆ x Q i (cid:16) u eq (cid:17) be the discrete initial conditions. Then the iterative process for solving (2.13)will, in each cell, find the local equilibrium h , i = h eq ( x i ) and s , i = s . We prove this fact for ideal gasesin Appendix A. Therefore, in every cell δ U i ( x ) = R ( x ; { } k ∈ S i ) =
0. Hence, we have U i + / − = U i + / + and by consistency of the numerical flux F i + / = f ( U i + / − ) = [0 , p eq ( x i + / ) , T . Likewise, by definition(2.18) the cell-averaged source term becomes S i = ∆ x [0 , p eq ( x i + / ) − p eq ( x i − / ) , T . By plugging the aboveexpressions for the numerical flux and source term into the semi-discrete finite volume scheme (2.4) we getd U i d t = L ( U ) = − ∆ x (cid:0) F i + / − F i − / (cid:1) + S i = (cid:4) Remark 2.3
The presented scheme reduces to the second-order accurate scheme presented in [9] by settingthe quadrature rule Q to the midpoint rule and the reconstruction procedure R to piecewise linear.2.2. Extension to several space dimensions We now describe the extension of our well-balanced scheme for hydrostatic equilibrium to the multi-dimensional case. For ease of presentation, we describe it for two dimensions and the extension to threedimensions is straightforward. As in the one-dimensional case, we briefly introduce a standard high-orderfinite volume scheme and then detail the well-balanced scheme.The two-dimensional Euler equations with gravity in Cartesian coordinates are given by ∂ u ∂ t + ∂ f ∂ x + ∂ g ∂ y = s (2.19)with u = ρρ v x ρ v y E , f = ρ v x ρ v x + p ρ v y v x ( E + p ) v x , g = ρ v y ρ v x v y ρ v y + p ( E + p ) v y and s = s x + s y = − ρ − ρ v x ∂φ∂ x + − ρ − ρ v y ∂φ∂ y , (2.20)where u is the vector of conserved variables, f and g the fluxes in x - and y -direction, and s the gravitationalsource terms.We consider a rectangular spatial domain Ω = [ x min , x max ] × [ y min , y max ] discretized uniformly (for easeof presentation) by N x and N y cells or finite volumes in x - and y -direction, respectively. The cells are labeledby I i , j = I i × I j = [ x i − / , x i + / ] × [ y j − / , y j + / ] and the constant cell sizes by ∆ x = x i + / − x i − / and ∆ y = y j + / − y j − / . We denote the cell centers by x i = ( x i − / + x i + / ) / y j = ( y j − / + y j + / ) / c over the cell faces are approximated by q -th order accurate quadrature rules as Q i ± / , j ( c ) = N q (cid:88) β = ω β c ( x i ± / , y j ,β ) ≈ (cid:90) I j c ( x i ± / , y ) d xQ i , j ± / ( c ) = N q (cid:88) α = ω α c ( x i ,α , y j ± / ) ≈ (cid:90) I i c ( x , y i ± / ) d y , (2.21)where the x i ,α ∈ I i , y j ,β ∈ I j and ω α , ω β denote the N q quadrature nodes and weights, respectively. Likewise,integrals over the cells are approximated by Q i , j ( c ) = N q (cid:88) α = N q (cid:88) β = ω α ω β c ( x i ,α , y j ,β ) ≈ (cid:90) I i , j c ( x , y ) d x d y . (2.22)A semi-discrete finite volume scheme for the numerical approximation of (2.19) then takes the followingform d U i , j d t = L ( U ) = − ∆ x (cid:16) F i + / , j − F i − / , j (cid:17) − ∆ y (cid:16) G i , j + / − G i , j − / (cid:17) + S i , j , (2.23)where U i , j denotes the approximate cell averages of the conserved variables, F i ± / , j and G i , j ± / the facialaverages of the fluxes through the cell boundary and S i , j the cell averages of the source term. The fluxes areobtained by applying the above quadrature rules along the cell boundary to the numerical flux formulas F and G in respective direction: F i + / , j = ∆ y Q i + / , j (cid:16) F ( U i , j , U i + , j ) (cid:17) G i , j + / = ∆ x Q i , j + / (cid:16) G ( U i , j , U i , j + ) (cid:17) , (2.24)where U i , j = U i , j ( x , y ) is a suitable reconstruction to be defined in detail at a later point. Similarly, the sourceterm is obtained by quadrature over the cell S i , j = ∆ x ∆ y Q i , j ( s ( U , ∇ φ )) . (2.25)In the evaluation of the quadrature rules, a r -th reconstruction procedure R is used to obtain pointwiserepresentations of the solution from the cell-averaged conserved variables: U i , j ( x , y ) = R (cid:18) x , y ; (cid:110) U k , l (cid:111) ( k , l ) ∈ S i , j (cid:19) . (2.26)Here S i , j is the stencil of the reconstruction for cell I i , j . Many such reconstruction procedures have beendeveloped in the literature and we refer to the references previously mentioned in section 2.1.1.As in the one-dimensional case, we need two ingredients to construct our well-balanced scheme. Thefirst is a high-order equilibrium preserving reconstruction and the second is a well-balanced discretizationof the momentum source terms.Let us begin with the description of the first ingredient and consider cell I i , j . Then the high-order equi-librium preserving reconstruction W takes the following form U i , j ( x , y ) = W (cid:18) x , y ; (cid:110) U k , l (cid:111) ( k , l ) ∈ S i , j (cid:19) = U eq , i , j ( x , y ) + δ U i , j ( x , y ) , (2.27)9hich again separates the solution into an equilibrium U eq , i , j and a (possibly large) perturbation δ U i , j .The equilibrium profile is built from (1.10), which is indeed also valid in more than one dimensions.Hence, we construct the local equilibrium profile in cell I i , j by h eq , i , j ( x , y ) = h , i , j + φ i , j − φ ( x , y ) , (2.28)where h , i , j = h eq , i , j ( x i , y j ) and φ i , j = φ ( x i , y j ) are the point values of the specific enthalpy and the gravita-tional potential at cell center, respectively. Given the (constant) equilibrium entropy s , i , j , the equilibriumprofiles of density ρ eq , i , j and internal energy density ρ e eq , i , j can be computed through the EoS.The equilibrium enthalpy at cell center h , i , j and the (constant) entropy s , i , j are again fixed by demandingagreement with the local cell averages up to the desired order of accuracy: ρ i , j = ∆ x ∆ y Q i , j ( ρ eq , i , j ) ρ e i , j = ∆ x ∆ y Q i , j ( ρ e eq , i , j ) . (2.29)Here ρ e i , j is the cell average of the internal energy density, which we estimate simply from the cell-averagedconserved variables by ρ e i , j = E i , j − ρ i , j (cid:16) ρ v x , i , j + ρ v y , i , j (cid:17) , (2.30)The latter estimate is again exact at equilibrium. As in the one-dimensional case, these equations represent,in general, a nonlinear system of two equations in the equilibrium specific enthalpy at cell center h , i , j andthe (constant) specific entropy s , i , j . Their resolution proceeds as in the one-dimensional case. In the end,we have the following equilibrium profile U eq , i , j ( x , y ) = ρ eq , i , j ( x , y )00 ρ e eq , i , j ( x , y ) . (2.31)The perturbation part is reconstructed as in the one-dimensional case by δ U i , j ( x , y ) = R (cid:18) x , y ; (cid:110) U k , l − Q k , l ( U eq , i , j ) (cid:111) ( k , l ) ∈ S i , j (cid:19) . (2.32)This simply extrapolates the cell’s local equilibrium profile, computes equilibrium cell averages by numer-ical integration, and uses the standard reconstruction procedure to obtain a high-order representation of theperturbation.We observe that the reconstruction procedure (2.27) preserves the equilibrium by construction, since δ U i , j vanishes, and it is min( q , r )-th order accurate in and away from equilibrium (for su ffi ciently smoothsolutions).Like in the one-dimensional case, only the momentum source terms need to be modified. The well-balanced momentum source terms are simply obtained on a dimension-by-dimension basis S ρ v x , i , j = ∆ x (cid:16) Q i + / , j ( p eq , i , j ) − Q i − / , j ( p eq , i , j ) (cid:17) − ∆ x ∆ y Q i , j (cid:32) δρ i , j ∂φ∂ x (cid:33) S ρ v y , i , j = ∆ y (cid:16) Q i , j + / ( p eq , i , j ) − Q i , j − / ( p eq , i , j ) (cid:17) − ∆ x ∆ y Q i , j (cid:32) δρ i , j ∂φ∂ y (cid:33) . (2.33)10his completes the description of the two-dimensional well-balanced scheme for hydrostatic equilibriumand its properties are summarized in the corollary below: Corollary 2.4
Consider the scheme (2.23) with consistent and Lipschitz continuous numerical fluxes F and G , a r-th order accurate spatial reconstruction procedure R , a q-th order accurate quadrature rule Q , thehydrostatic reconstruction W (2.27) and the gravitational source term S (2.25) (with (2.33) ).This scheme has the following properties: (i) The scheme is consistent with (2.19) and it is min( q , r ) -th order accurate in space (for smooth solu-tions). (ii) The scheme is well-balanced and preserves the discrete hydrostatic equilibrium given by (1.10) andv x = v y = exactly. Proof
The proof follows directly by applying theorem 2.2 dimension-by-dimension. (cid:4)
3. Numerical Experiments
In this section we assess the performance of our well-balanced scheme on a series of numerical exper-iments. For comparison, we also present results obtained with a standard (unbalanced) base scheme. Thefully-discrete finite volume base scheme consists of • the temporally third-order accurate SSP-RK scheme for time integration (see [31]), • the spatially third-order accurate CWENO3 [35] reconstruction procedure R , • the spatially fourth-order accurate two-point Gauss-Legendre quadrature rule for Q .Overall the scheme is third-order accurate in space and time. This scheme is conditionally stable under theusual CFL condition. We use a CFL number of C CFL = .
85. In the following, we will refer to this schemeas the unbalanced scheme. The well-balanced scheme is built with the same base components, but uses thewell-balanced reconstruction procedure and source term computation as outlined in the previous section.Below, all the initial conditions will be given in functional form u ( x ). The discrete initial conditions areobtained simply by quadrature, i.e. U i = Q i ( u ) , U i , j = Q i , j ( u ) (3.1)in the one- and two-dimensional case, respectively. It is important to notice that the well-balancing onlyrequires that the initial conditions are obtained by the exact same quadrature rule used in the numericalscheme. Therefore, showing that the initial conditions are well-balanced will immediately imply that thepreserved discrete state is a high-order accurate approximation of the exact equilibrium, simply because thediscrete initial conditions are nothing else than a high-order quadrature of the exact equilibrium.We will be using three distinct notions of “error”. The first error is the usual L -error err ( q ) : = (cid:88) i ∆ x | q i − q re f , i | , (3.2)where q is any scalar variable of interest, e.g. q = ρ, p , v , . . . . Furthermore, q re f , i is computed by down-sampling a high-resolution reference solution or, where available, a highly accurate quadrature of an analyticsolution. A subtlety is that even in a well-balanced scheme the err of a discrete preserved state is not, ingeneral, zero. 11o answer the question of how big the error of a perturbation δ q from equilibrium is, we define the L -error of δ q as err ( δ q ) : = (cid:88) i ∆ x | ( q i − Q i ( q eq )) − δ q re f , i | (3.3)where q eq is the background equilibrium profile and δ q re f , i is the cell-average of the perturbation in a ref-erence solution. This measures the error of the perturbation from numerical equilibrium. This is subtlydi ff erent than the error of the perturbation from the exact equilibrium. The di ff erence is that err ( δ q ) conve-niently uses the quadrature rule used in the finite volume method to compute the average of the equilibriumprofile, i.e. Q i , j q eq . Therefore, for equilibria, a well-balanced scheme should have zero err ( δ q ), but mayhave non-zero err ( q ).When computing the err ( δ q ) for hydrostatic equilibria, the reference solution Q i , j ( q eq ) is known exactly,it’s simply the initial condition. Therefore, err ( δ q ) can be computed at a greatly reduced computational costby err eq , ( q ) : = (cid:88) i ∆ x | q i − Q i ( q eq ) | (3.4)In order to be clear about how the errors where computed we will make the distinction throughout thenumerical experiments. Moreover, the above error measures readily generalize to the two-dimensional case.To characterize a time scale on which a model reacts to perturbations of its equilibrium, we define thesound crossing time τ sound τ sound = (cid:90) c − s d x , (3.5)where c s denotes the speed of sound and the integral has to be taken over the extent of the stationary state ofinterest. The sound crossing time is basically the time in which a sound wave travels back and forth throughthe equilibrium.We begin by several simple one- and two-dimensional numerical experiments employing the ideal gasEoS. The interested reader may readily reproduce these experiments in order to check his or her implementa-tion. Finally, we demonstrate the performance of the scheme on a problem involving a complex multiphysicsEoS. We consider an isentropic hydrostatic atmosphere in a constant gravitational field. The gravitationalpotential is a linear function φ ( x ) = gx where g is the constant gravitational acceleration. The initial densityand pressure profiles are then given by ρ ( x ) = (cid:32) K γ − γ ( h − gx ) (cid:33) /γ − , p ( x ) = K ρ ( x ) γ + A exp (cid:32) − ( x − / . (cid:33) . (3.6)with the constants g = . γ = . h = .
75 and K =
1. The atmosphere’s pressure is perturbed by aGaussian bump of amplitude A . The velocity is set to zero everywhere.The computational domain is set to [0 , L ] with L = N cells, i.e. we set thecell size ∆ x = L / N , the cell interfaces x i + / = i ∆ x and the cell centers x i = ( x i − / + x i + / ) / i = , ..., N .The following resolutions are used N = , , , , , U i = Q i ( U eq , ) for i < U i = Q i ( U eq , N ) for i > N . (3.7) We first verify the well-balanced property of our scheme. For this we evolve the isentropic atmospherein hydrostatic equilibrium without pressure perturbation, A =
0, up to time t =
10. This corresponds toroughly 6 sound crossing time ( τ sound = . weno weno err eq , ( ρ ) rate err eq , ( ρ ) rate32 8 . × − – 9 . × − –64 1 . × − . × − -2.01128 1 . × − . × − . × − . × − -0.80512 1 . × − . × − -0.501024 1 . × − . × − -0.89 Table 1: Convergence data for the one-dimensional test Section 3.1.1 without perturbation. We show err eq , ( ρ ) at t = .
0. On the leftwe show the error for the unbalanced scheme. It converges at slightly higher rates than expected. This is likely due to the fact that theinitial conditions are a fourth order approximation of the exact cell-averages. Clearly, the numerical solution isn’t stationary and thetruncation error has accumulated over time. The right hand side column shows the errors for the well-balanced scheme. Note that theerrors are at the level of round-o ff . This result also implies that the preserved discrete state is a fourth order approximation of the exactequilibrium. Next we add a small pressure perturbation to the isentropic atmosphere in order to examine the schemesability to propagate small waves. The amplitude of the pressure perturbation is set to A = − , whichgenerates one smooth wave propagating upwards and one downwards through atmosphere. As the wavespropagate, they are modified by the density and pressure stratification of the atmosphere. We evolve thesetup until time t = .
2, shortly before the waves reach the boundaries.The errors of the density perturbation err ( δρ ) are shown in Table 2. The density perturbation is thedensity at the final time minus the density of the unperturbed atmosphere. These errors were obtained onthe basis of a reference solution computed by the unbalanced scheme with a high resolution N =
32 768.We observe that the errors of the well-balanced scheme are roughly four orders of magnitude smaller thanthe unbalanced scheme. The convergence rate of both the unbalanced and well-balanced reach the expectedrate of three. The somewhat irregular convergence rates of the unbalanced scheme can be explained by the13cheme still being (heavily) pre-asymtotic at the lower resolutions. The slow convergence rate of the well-balanced scheme is a feature of the well-balanced scheme. Since it was designed to have very small errorsclose to equilibrium.In Figure 1 the profile of the velocity and the pressure perturbation are shown at the final time for boththe unbalanced (blue crosses) and well-balanced (red circle) schemes. The well-balanced solution is shownfor N =
64. Even at this relatively low resolution the well-balanced scheme resolves the perturbation well.The errors of the unbalanced scheme for N =
64 are too small to be shown on the same plot. Instead we plotthe solution of the unbalanced scheme at N = weno weno err ( δρ ) rate err ( δρ ) rate32 5 . × − – 3 . × − –64 6 . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − Table 2: Convergence data for the one-dimensional test Section 3.1.2 with the small perturbation A = − . We show the error of thedensity perturbation err ( δρ ) at t = .
2. On the left we show the errors for the unbalanced scheme. The first observation is that the errorsare large and even exceed the size of the perturbation, the second is that the convergence rates are not quite as expected. The reason forthe second observation is that given the very small perturbation we are trying to resolve, the scheme is likely still pre-asymptotic. Theerrors of the perturbation for the well-balanced scheme is given in the right hand side column. The overall error is much smaller, andless than the amplitude of the perturbation. Furthermore, the scheme converges at the expected rate.Figure 1: Snapshot of the smooth test case, see Section 3.1.2. The simulation is performed with N =
256 and N =
64 cells for theunbalanced and well-balanced scheme respectively. It is run up to time t = .
2. On the left we show the velocity, on the right thepressure perturbation. Even with N =
64 cells, the well-balanced scheme can resolve the small perturbation cleanly. The error in theunbalanced scheme on the other hand is too big to be shown on the plot. The reference solution, plotted in black, is obtained by ahigh-resolution N =
32 768 simulation using the unbalanced scheme. .1.3. Large pressure perturbation propagation For the purpose of testing that the well-balanced reconstruction does not destroy the robustness of theshock-capturing base scheme, we increase the pressure perturbation by several orders of magnitude to A =
10. This generates two strong waves quickly steepening into shock waves. The setup is evolved until time t = . ff ected the performanceof the scheme away from equilibrium. Figure 2: Snapshot of one-dimensional test with large perturbation A =
10, see Section 3.1.3. The simulation is performed with N =
64 cells and run up to time t = .
06. On the left the velocity is shown, on the right the pressure. The well-balanced scheme isnearly indistinguishable from the unbalanced scheme. Clearly, the well-balancing has not a ff ected the quality of the numerical solutionaway from equilibrium. Furthermore, no spurious oscillations are observed. The reference solution, plotted in black, is obtained by ahigh-resolution N =
32 768 simulation using the unbalanced scheme. .2. Two-dimensional polytrope The following numerical experiment is a two-dimensional version of the one in [9]. This experimentsimulates a so-called polytrope, which is a static configuration of an adiabatic gaseous sphere held togetherby self-gravitation. These model stars are constructed in spherical symmetry from hydrostatic equilibrium,Poisson’s equation and the polytropic relation p = K ρ γ , which can be combined into the so-called Lane-Emden equation (see e.g. [36]). The latter equation can be solved analytically for three values of the ratio ofspecific heats ( γ = / , , ∞ ).As in [9] we use γ =
2. Then the density and pressure profiles are given by ρ ( r ) = ρ C sin( α r ) α r , p ( r ) = K ρ ( r ) γ (3.8)where r is the radial coordinate, ρ C is the central density of the polytrope and α = (cid:114) π G K . (3.9)The gravitational potential is given by φ ( r ) = − K ρ C sin( α r ) α r . (3.10)In the following we set K = G = ρ C =
1. Note that the polytrope (obviously) fulfills the equilibrium (1.10)for any r ≥ − . , . by N uniform cells for N = , , , , , u ( x , y ) = [ ρ ( r ) , , , p ( r ) / ( γ − T where the radial coordinate is given by r = x + y . Notethat the velocity is set to zero in the whole domain.The boundary conditions are applied along the coordinates axes as in Section 3.1. In the corner bound-aries (needed by the reconstruction procedure), we extrapolate the equilibrium from the relevant corner cellin the computational domain. For example, the ghost cells in the upper right corner are set as follows U i , j = Q i , j ( U eq , N , N ) for N < i , j . (3.11)The gravitational potential is simply given by the above analytical expression. We begin by evolving the polytrope with the well-balanced and unbalanced scheme until time t = τ sound ≈ . L norm. Furthermore, italso works for equilibria which are not grid aligned. The unbalanced scheme, however, su ff ers from largespurious deviations. Next we add a perturbation to the equilibrium pressure of the polytrope as p ( r ) = (cid:16) + A exp( − r / . ) (cid:17) p ( r ) (3.12)16 igure 3: Scatter plots of the two-dimensional polytrope, see Section 3.2.1 at time t = .
2. The panel on the left shows the velocityfor A = − , the one on the right for A = − . The resolution is generally N = , however, for the smallest perturbation, theerrors of the unbalanced scheme at N =
128 exceed the limits of the plot. Therefore, for A = − , we plot the unbalanced scheme at N = N =
32 768 cells.Figure 4: Scatter plots of the two-dimensional polytrope, see Section 3.2.1. The simulation is performed with N = cells and runup to time t = .
2. The panel on the left shows the velocity for A = − , the one on the right for A =
10. For these larger perturbationswe see that the unbalanced and well-balanced scheme agree very well. In fact, Table 4 shows that the errors are on the same orderof magnitude for A = − and di ff er only about a percent for A =
10. Clearly the well-balancing does not a ff ect the quality of theapproximate solution away from equilibrium. Furthermore, A =
10 shows that the well-balanced scheme resolves discontinuities just aswell as the unbalanced scheme. The reference solution was computed with a one-dimensional finite volume code assuming cylindricalsymmetry on N =
32 768 cells.
17 C weno weno err eq , ( ρ ) rate err eq , ( ρ ) rate32 1 . × − – 2 . × − –64 1 . × − . × − . × − . × − . × − . × − -1.00512 3 . × − . × − -0.951024 3 . × − . × − -0.95 Table 3: Convergence data for the polytrope at rest, see Section 3.2.1. We show err eq , ( ρ ) for the unbalanced scheme (left) and thewell-balanced scheme (right). The unbalanced scheme converges at the expected rate, but even with 1024 cells, it has not reachedround o ff . The well-balanced scheme is again shown to be in fact well-balanced. Like in the one-dimensional test, this implies that thepreserved discrete state is fourth order accurate. This shows that the scheme also works in two-dimensions, even in cases where thegravity is non-constant and not grid-aligned. with three di ff erent amplitudes A = − , − , − ,
10. The setup is evolved up to time t = . N =
32 768.For perturbations of size A = − the well-balanced scheme clearly outperforms the unbalanced scheme(by at least four orders of magnitude). Scatter plots of the velocity and pressure perturbation are shown inFigure 3. At N = the well-balanced scheme resolves the perturbation well and with no discerniblescatter. Which is not trivial, since the radially symmetric solution is approximated on a uniform Cartesiangrid which does not respect the radial symmetry.At the next larger perturbation, A = − the well-balanced scheme still outperforms the unbalancedscheme by a factor 10. Unlike the unbalanced scheme, the well-balanced scheme shows no scatter, as canbe seen in Figure 3. Furthermore, once the perturbation has traveled away from the center of the domain,the solution returns back to equilibrium in the well-balanced scheme, but not in the unbalanced one. Bothschemes converge at the expected rate.The second largest perturbation, A = − , was chosen such that the perturbation is very well approxi-mated by the unbalanced scheme, yet small enough to not develop any discontinuities. The aim is to showthat away from equilibrium the well-balancing does not have a negative impact on the quality of the solution.This is confirmed in Figure 4 and Table 4.For A =
10 both schemes perform equally well, both converge at first order and the errors di ff er byapproximately one percent. Therefore, well-balancing has not a ff ected the quality of the solution away fromequilibrium. The scatter plot of the velocity and pressure is shown in Figure 4. Neither scheme shows anysign of spurious oscillations. In order to further verify that our well-balanced scheme does not deteriorate the robustness and shock-capturing properties of the base scheme, we add to the polytrope several localized high pressure regions. Tothis end, we add the following pressure perturbation to the equilibrium polytrope δ p ( x ) = (cid:88) i = B ( x i , r ) ( x ) , (3.13)18 C weno weno err ( δρ ) rate err ( δρ ) rate32 5 . × − – 5 . × − –64 6 . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − weno weno err ( δρ ) rate err ( δρ ) rate32 5 . × − – 5 . × − –64 6 . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − weno weno err ( δρ ) rate err ( δρ ) rate32 8 . × − – 5 . × − –64 2 . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − weno weno err ( δρ ) rate err ( δρ ) rate32 2 . × − – 2 . × − –64 1 . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − Table 4: Convergence data for the polytrope with perturbation. We show err ( δρ ) at t = . A = − . Clearly, the well-balanced scheme outperformsthe unbalanced scheme. Furthermore, the expected rate is observed until round o ff sets in at N = . The second table shows theerrors for the medium perturbation, A = − . The well-balanced scheme is still slightly better than the unbalanced one. However thereal benefit can be seen much more clearly in the scatter plots, c.f. Figure 3. The third table is for A = − . We clearly see that thewell-balancing had no negative a ff ect on the quality of the solution, even though the solution is no longer near equilibrium. Finally,the fourth table shows the case for A =
10. The smooth perturbation turns into a discontinuity and only first order convergence can beexpected. It’s interesting to see that the error of the unbalanced and well-balanced scheme di ff er by only about one percent. B x , R = { x (cid:48) ∈ R : (cid:107) x (cid:48) − x (cid:107) < R } denotes the open ball of radius R centered on x and B the indicatorfunction for the set B , i.e. B ( x ) = x ∈ B , R = .
05 and centers x = [ − . , . T , x = [ − . , . T , x = [0 . , . T , x = [0 . , . T , x = x = [0 . , − . T . The velocity is set to zero everywhere. The initial conditions are shown in the upper panel of Figure 5.We evolve the setup until time t = .
02 with the well-balanced and unbalanced scheme at resolution N = We show a snapshot at t = .
02 in Figure 5. Even under these much more extreme conditionswith non-trivial wave interactions, the well-balanced scheme is stable and by eye indistinguishable from theunbalanced scheme.
The final numerical experiment assesses the performance of our well-balanced scheme on a astrophys-ically relevant problem involving a complex multiphysics EoS. We simulate the equilibrium and some per-turbations of a model white dwarf. A white dwarf is the final evolutionary state of a star not massive enoughto go through the final nuclear burning stages and become a neutron star or a black hole (see e.g. [37]).This numerical experiment is a two-dimensional version of the one presented in [9]. Likewise, we usethe publicly available Helmholtz EoS (see [38] for a detailed descriptions and [39]). This EoS includescontributions of (photon) radiation, nuclei, electrons and positrons and is well adapted to large range ofstellar environments. The radiation is treated as a blackbody in local thermal equilibrium and the nucleiare modeled by the ideal gas law. For computational e ffi ciency, the electrons and positrons are treated in atabular manner with a thermodynamically consistent interpolation procedure.The white dwarf model is fully characterized by specifying the central density, the chemical compositionand the thermodynamic equilibrium. We set the central density ρ = × g / cm and temperature T = × K. We assume a constant specific entropy and set the composition to half carbon C and halfoxygen O. Then the model can be constructed by simple numerical integration of the self-gravitatinghydrostatic equilibrium equations in spherical symmetry. We refer to [9] for the detailed procedure.The one-dimensional white dwarf profile is then mapped onto the two-dimensional computational do-main [ − L , L ] with L = × cm. The velocity is set to zero. The same hydrostatic extrapolation boundaryconditions are used as in Section 3.2. We evolve the hydrostatic equilibrium without perturbation on a grid with N = cells until time t = err eq , ( ρ ) = . × g / cm and err eq , ( E ) = . × erg / cm .The well-balanced scheme is confirmed to be exact up to machine precision, with errors of err eq , ( ρ ) = . × − g / cm and err eq , ( E ) = . × erg / cm . This shows that equation (2.29) can be solvednumerically and the solution is e ff ectively unique. If the iterative procedure where to find a di ff erent equi-librium in any cell, it would be very unlikely that the resulting scheme would be well-balanced.20 igure 5: Snapshot of the two dimensional blast waves, see Section 3.2.3. The upper image shows the initial total energy (on a linearscale). The remaining three plots show the total energy at t = .
02 (on a logarithmic scale) for the unbalanced (bottom, left) andwell-balanced (bottom, right) scheme. Even for these extreme initial conditions the unbalanced and well-balanced schemes performequally well. .3.2. Wave propagation Next we add a small Gaussian pressure perturbation at the origin, i.e. p ( x ) = (1 + A exp( −| x | / b )) p eq ( x ) , (3.14)with A = − and b = × cm. The solution is evolved to t = . × − s on N = . A scatterplot of the solution is shown in Figure 6. The scatter in the well-balanced scheme is significantly reducedcompared to the unbalanced scheme. Unlike the unbalanced solution, the well-balanced solution remainsconstant ahead of the perturbation and returns to rest after the perturbation has passed.The reference solution is computed using a one-dimensional, cylindrically symmetric, well-balancedfinite volume code with a resolution of N = Figure 6: Snapshots of the two-dimensional white-dwarf, see Section 3.3.2 with a small perturbation. The simulation is performed with N = cells and run up to time t = . × − s. The radial velocity is shown on the left, the pressure perturbation on the right.The size of the initial perturbation was chosen such that the main feature of the perturbation is resolved similarly well in both solvers.However, the well-balanced scheme has significantly less scatter and shows no deviation from equilibrium ahead of the perturbation.Furthermore, the well-balanced scheme returns to rest after the wave has moved away from the center of the domain, i.e. r =
4. Conclusion
We presented a novel well-balanced, high-order finite volume scheme for Euler equations with gravity.We are able to well-balance a large class of astrophysically relevant hydrostatic equilibria without imposingthe exact equilibrium apriori. Rather, we only assume some thermodynamic information about the equilib-rium, e.g. constant entropy, and solve for the equilibrium in every time step. Since the equilibrium definedby ∇ p = − ρ ∇ φ (4.1)is only a mechanical equilibrium, it seems natural that some additional information about the thermodynamicnature of the equilibrium must always be imposed.The important features of the proposed scheme are: • Its independence of a particular form of equation of state. This scheme can handle arbitrary equationsof state including tabulated ones, as shown in the final numerical experiment.22
Its independence of a particular gravitational potential. The only requirement is that the gravitationpotential and its gradient can be evaluated at apriori known locations in the computation domain. Infact the gravitation source term does not need to be constant in time. Therefore this scheme can also beapplied to simulations which include self-gravity. Such simulations may benefit from well-balancingif the initial conditions are at rest and perturbed by some other means, such as a heating source term. • Its modular nature. The scheme clearly describes how any reconstruction procedure can be made well-balanced. Therefore, the proposed scheme can be extended to arbitrary orders in a straightforwardmanner. • Its local nature. The well-balancing is local to each cell. In particular it does not change the stencilrequired to update the cell.The numerical experiments have shown that the scheme is high-order accurate for flows both near andfar away from hydrostatic equilibrium. In fact, the numerical results suggest that the scheme is no worse onlarge perturbations than the equivalent unbalanced scheme. We have also shown that the schemes are stablein the presence of shocks. Furthermore, the smooth tests show that the well-balanced scheme preservesradial symmetry much better than the unbalanced scheme. Furthermore, the well-balanced solutions do notcause any changes in the part of the domain the perturbation has not reached yet. Additionally, the well-balanced scheme returns to rest after the perturbation has passed over some region in the domain, while theunbalanced scheme does neither. These tests were performed under a variety of di ff erent conditions, i.e. inone dimension for constant gravity, in two dimensions for non-grid aligned gravity with both the ideal gaslaw and a complex multiphysics equation of state.Our scheme only a ff ects the reconstruction procedure and the numerical source term. Therefore, largeparts of an existing finite volume code would remain una ff ected by adding our well-balancing. By reusingthe existing unbalanced reconstruction procedure for the perturbation the cost of implementing our schemeis further reduced. These very localized and modular changes ensure that our method can be used to well-balance a variety of di ff erent existing finite volume schemes with minimal e ff ort. Acknowledgments
The work was supported by the Swiss National Science Foundation (SNSF) under grant 200021-169631.The authors also acknowledge the use of computational resources provided by the Swiss SuperComputingCenter (CSCS), under the allocation grant s661, s665, s667 and s744. We acknowledge the computationalresources provided by the EULER cluster of ETHZ.
References [1] J. Greenberg, A. Leroux, A well-balanced scheme for the numerical processing of source terms inhyperbolic equations, SIAM Journal on Numerical Analysis 33 (1) (1996) 1–16. arXiv:http://epubs.siam.org/doi/pdf/10.1137/0733001 , doi:10.1137/0733001 .URL http://epubs.siam.org/doi/abs/10.1137/0733001 [2] R. J. LeVeque, Balancing source terms and flux gradients in high-resolution godunov methods: Thequasi-steady wave-propagation algorithm, J. Comput. Phys. 146 (1) (1998) 346–365.URL doi:10.1137/S1064827503431090 .URL http://link.aip.org/link/?SCE/25/2050/1 [4] S. Noelle, N. Pankratz, G. Puppo, J. R. Natvig, Well-balanced finite volume schemes of arbitrary orderof accuracy for shallow water flows, J. Comput. Phys. 213 (2) (2006) 474–499. doi:10.1016/j.jcp.2005.08.019 .URL http://dx.doi.org/10.1016/j.jcp.2005.08.019 [5] L. Gosse, Computing Qualitatively Correct Approximations of Balance Laws, Springer Milan, 2013. doi:10.1007/978-88-470-2892-0 .[6] R. J. LeVeque, D. S. Bale, Wave propagation methods for conservation laws with source terms, in:Hyperbolic problems: theory, numerics, applications, Vol. II (Z¨urich, 1998), Vol. 130 of Internat. Ser.Numer. Math., Birkh¨auser, Basel, 1999, pp. 609–618.[7] N. Botta, R. Klein, S. Langenberg, S. L¨utzenkirchen, Well balanced finite volume methodsfor nearly hydrostatic flows, Journal of Computational Physics 196 (2) (2004) 539 – 565. doi:DOI:10.1016/j.jcp.2003.11.008 .URL [8] R. J. LeVeque, A well-balanced path-integral f-wave method for hyperbolic problems with sourceterms, Journal of Scientific Computing 48 (2011) 209–226. doi:10.1007/s10915-010-9411-0 .[9] R. K¨appeli, S. Mishra, Well-balanced schemes for the euler equations with gravitation, Journal of Com-putational Physics 259 (0) (2014) 199 – 219. doi:http://dx.doi.org/10.1016/j.jcp.2013.11.028 .URL [10] V. Desveaux, M. Zenk, C. Berthon, C. Klingenberg, A well-balanced scheme for the Euler equationwith a gravitational potential, Springer Proceedings in Mathematics & Statistics (2014) 217’ ¨A`ı226 doi:10.1007/978-3-319-05684-5_20 .URL http://dx.doi.org/10.1007/978-3-319-05684-5_20 [11] P. Chandrashekar, C. Klingenberg, A second order well-balanced finite volume scheme for euler equa-tions with gravity, SIAM Journal on Scientific Computing 37 (3) (2015) B382–B402. arXiv:http://dx.doi.org/10.1137/140984373 , doi:10.1137/140984373 .URL http://dx.doi.org/10.1137/140984373 [12] R. K¨appeli, S. Mishra, A well-balanced finite volume scheme for the Euler equations with gravitation.The exact preservation of hydrostatic equilibrium with arbitrary entropy stratification, Astronomy andAstrophysics 587 (2016) A94. doi:10.1051/0004-6361/201527815 .[13] R. Touma, U. Koley, C. Klingenberg, Well-balanced unstaggered central schemes for the euler equa-tions with gravitation, SIAM Journal on Scientific Computing 38 (5) (2016) B773–B807. doi:10.1137/140992667 .[14] G. Li, Y. Xing, High order finite volume WENO schemes for the euler equations under gravitationalfields, Journal of Computational Physics 316 (2016) 145–163. doi:10.1016/j.jcp.2016.04.015 .2415] R. K¨appeli, A well-balanced scheme for the euler equations with gravitation, in: Innovative Al-gorithms and Analysis, Springer International Publishing, 2017, pp. 229–241. doi:10.1007/978-3-319-49262-9_8 .[16] A. Chertock, S. Cui, A. Kurganov, S¸ . N. ¨Ozcan, E. Tadmor, Well-balanced schemes for the eulerequations with gravitation: Conservative formulation using global fluxes, Journal of ComputationalPhysics 358 (2018) 36–52. doi:10.1016/j.jcp.2017.12.026 .[17] Y. Xing, C.-W. Shu, High order well-balanced WENO scheme for the gas dynamics equations undergravitational fields, J. Sci. Comput. 54 (2-3) (2013) 645–662. doi:10.1007/s10915-012-9585-8 .URL http://dx.doi.org/10.1007/s10915-012-9585-8 [18] G. Li, Y. Xing, Well-balanced finite di ff erence weighted essentially non-oscillatory schemes for theeuler equations with static gravitational fields, Computers & Mathematics with Applications 75 (6)(2018) 2071–2085. doi:10.1016/j.camwa.2017.10.015 .[19] G. Li, Y. Xing, Well-balanced discontinuous galerkin methods for the euler equations under gravita-tional fields, Journal of Scientific Computing (2015) 1–21 doi:10.1007/s10915-015-0093-5 .URL http://dx.doi.org/10.1007/s10915-015-0093-5 [20] P. Chandrashekar, M. Zenk, Well-balanced nodal discontinuous galerkin method for euler equa-tions with gravity, Journal of Scientific Computing 71 (3) (2017) 1062–1093. doi:10.1007/s10915-016-0339-x .[21] G. Li, Y. Xing, Well-balanced discontinuous galerkin methods with hydrostatic reconstruction for theeuler equations with gravitation, Journal of Computational Physics 352 (2018) 445–462. doi:10.1016/j.jcp.2017.09.063 .[22] F. Fuchs, A. McMurry, S. Mishra, N. Risebro, K. Waagan, High order well-balanced finite volumeschemes for simulating wave propagation in stratified magnetic atmospheres, Journal of ComputationalPhysics 229 (11) (2010) 4033 – 4058. doi:DOI:10.1016/j.jcp.2010.01.038 .URL [23] P. Batten, N. Clarke, C. Lambert, D. M. Causon, On the choice of wavespeeds for the hllc rie-mann solver, SIAM Journal on Scientific Computing 18 (6) (1997) 1553–1570. doi:10.1137/S1064827593260140 .URL http://link.aip.org/link/?SCE/18/1553/1 [24] E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction,Springer-Verlag GmbH, 1997.[25] B. van Leer, Towards the ultimate conservative di ff erence scheme. V. A second-order sequel togodunov’s method, Journal of Computational Physics 32 (1) (1979) 101–136. doi:10.1016/0021-9991(79)90145-1 .[26] A. Harten, High resolution schemes for hyperbolic conservation laws, Journal of ComputationalPhysics 49 (3) (1983) 357–393. doi:10.1016/0021-9991(83)90136-5 .2527] P. Colella, P. R. Woodward, The piecewise parabolic method (ppm) for gas-dynamical simulations, J.Comput. Phys. 54 (1) (1984) 174–201.URL [28] A. Harten, B. Engquist, S. Osher, S. R. Chakravarthy, Uniformly high order accurate essentially non-oscillatory schemes, III, Journal of Computational Physics 71 (2) (1987) 231–303. doi:10.1016/0021-9991(87)90031-3 .[29] C.-W. Shu, High order weighted essentially nonoscillatory schemes for convection dominated prob-lems, SIAM Review 51 (1) (2009) 82–126. doi:10.1137/070679065 .[30] I. Cravero, G. Puppo, M. Semplice, G. Visconti, CWENO: Uniformly accurate reconstructions forbalance laws, Mathematics of Computation (2017) 1 doi:10.1090/mcom/3273 .[31] S. Gottlieb, C.-W. Shu, E. Tadmor, Strong stability-preserving high-order time discretization methods,SIAM Review 43 (1) (2001) 89–112. doi:10.1137/S003614450036757X .URL http://link.aip.org/link/?SIR/43/89/1 [32] E. Godlewski, P.-A. Raviart, Numerical approximation of hyperbolic systems of conservation laws,Applied Mathematical Sciences doi:10.1007/978-1-4612-0713-9 .URL http://dx.doi.org/10.1007/978-1-4612-0713-9 [33] C. Hirsch, Numerical Computation of Internal and External Flows: The Fundamentals of Computa-tional Fluid Dynamics: The Fundamentals of Computational Fluid Dynamics, Vol. 1, Butterworth-Heinemann, 2007.[34] C. B. Laney, Computational Gasdynamics, Cambridge University Press, 1998. doi:10.1017/cbo9780511605604 .URL http://dx.doi.org/10.1017/CBO9780511605604 [35] D. Levy, G. Puppo, G. Russo, Compact central weno schemes for multidimensional conservation laws,SIAM Journal on Scientific Computing 22 (2) (2000) 656–672. arXiv:https://doi.org/10.1137/S1064827599359461 , doi:10.1137/S1064827599359461 .URL https://doi.org/10.1137/S1064827599359461 [36] S. Chandrasekhar, An introduction to the study of stellar structure, New York: Dover, 1967.[37] S. L. Shapiro, S. A. Teukolsky, Black Holes, White Dwarfs, and Neutron Stars, Wiley-VCH VerlagGmbH, 2007. doi:10.1002/9783527617661.fmatter .[38] F. X. Timmes, F. D. Swesty, The accuracy, consistency, and speed of an electron-positron equation ofstate based on table interpolation of the helmholtz free energy, The Astrophysical Journal SupplementSeries 126 (2) (2000) 501.URL http://stacks.iop.org/0067-0049/126/i=2/a=501 [39] F. X. Timmes, http://cococubed.asu.edu/code_pages/eos.shtml (2013). [link].URL http://cococubed.asu.edu/code_pages/eos.shtml ppendix A. Equilibrium reconstruction for the ideal gas law In the ideal gas case, it can be shown that a unique equilibrium exists which matches the cell-averages,i.e. satisfies (2.33) (in one dimension) and (2.13) (in two dimensions).In a first step the system is reduced to a single nonlinear equation in one unknown. To this end, we writethe ideal gas law in the polytropic form p = p ( K , ρ ) = K ρ γ , (A.1)where K = K ( s ) is a function of entropy s alone and γ is the ratio of specific heats. Then the equilibriumdensity and internal energy density can be expressed as functions of the constant K , i and the enthalpy at cellcenter h , i : ρ eq , i ( x ) = (cid:32) K , i γ − γ h eq , i ( x ) (cid:33) γ − ρ e eq , i ( x ) = γ − (cid:32) K , i (cid:33) γ − (cid:32) γ − γ h eq , i ( x ) (cid:33) γγ − . (A.2)By plugging the latter into (2.13), one obtains a single equation for h , i ρ e i = ρ i γ − (cid:80) N q j = w j (cid:16) γ − γ (cid:16) h , i + φ i − φ ( x j ) (cid:17)(cid:17) γγ − (cid:80) N q j = w j (cid:16) γ − γ (cid:16) h , i + φ i − φ ( x j ) (cid:17)(cid:17) γ − = : f ( h , i ) (A.3)and the constant K , i is simply given by K , i = ∆ x ρ i N q (cid:88) j = w j (cid:32) γ − γ (cid:16) h , i + φ i − φ ( x j ) (cid:17)(cid:33) γ − γ − . (A.4)To show that (A.3) has a unique solution we show that it is monotone. Therefore, we di ff erentiate f and find f (cid:48) ( h , i ) = ρ i γ − − γ (cid:80) N q j = w j (cid:16) h , i + φ i − φ ( x j ) (cid:17) γγ − · (cid:80) N q j = w j (cid:16) h , i + φ i − φ ( x j ) (cid:17) − γγ − (cid:18)(cid:80) N q j = w j (cid:16) h , i + φ i − φ ( x j ) (cid:17) γ − (cid:19) . (A.5)Clearly, the second term is positive, and if it where less than γ the derivative of f would be positive, every-where, and therefore f would be a strictly monotone function. If within every cell φ does not vary too much,this turns out to be true and can be proven as follows.Let h max , i = h , i + max x ∈ [ x i − / , x i + / ] φ i − φ ( x ) (A.6) h min , i = h , i + min x ∈ [ x i − / , x i + / ] φ i − φ ( x ) (A.7)27hen for 1 < γ ≤ (cid:80) N q j = w j (cid:16) h , i + φ i − φ ( x j ) (cid:17) γγ − · (cid:80) N q j = w j (cid:16) h , i + φ i − φ ( x j ) (cid:17) − γγ − (cid:18)(cid:80) N q j = w j (cid:16) h , i + φ i − φ ( x j ) (cid:17) γ − (cid:19) (A.8) ≤ (cid:80) N q j = w j h γγ − max , i · (cid:80) N q j = w j h − γγ − max , i (cid:18)(cid:80) N q j = w j h γ − min , i (cid:19) = h max , i h min , i . (A.9)Therefore, under the condition that h max , i h min , i < γ / (A.10) f has a unique solution. By a very similar estimate we find that for γ > h max , i h min , i < γ γ − γ ..