High order wellbalanced finite volume methods for multidimensional systems of hyperbolic balance laws
Abstract
We introduce a general framework for the construction of wellbalanced finite volume methods for hyperbolic balance laws. The phrase wellbalancing is used in a wider sense, since the method can be applied to exactly follow any solution of any system of hyperbolic balance laws in multiple spatial dimensions. The solution has to be known a priori, either as an analytical expression or as discrete data. The proposed framework modifies the standard finite volume approach such that the wellbalancing property is obtained. The potentially high order of accuracy of the method is maintained under the modification. We show numerical tests for the compressible Euler equations with and without gravity source term and with different equations of state, and for the equations of compressible ideal magnetohydrodynamics. Different grid geometries and reconstruction methods are used. We demonstrate high order convergence numerically.
keywords:
finitevolume methods, wellbalancing, hyperbolic balance laws, compressible Euler equations with gravity, ideal magnetohydrodynamics1 Introduction
Several problems in engineering and science are modeled by conservation properties. A common approach is to describe their evolution using hyperbolic conservation laws. However, only in few cases can these equations be solved analytically. Numerical methods to approximate solutions are used to forecast the behavior of the system. One successful approach amongst these is the finite volume method, based on Godunov’s method Godunov (1959). Finite volume methods introduce discretization errors. Typically, they are only exact on the constant state, some polynomials (in the case of higher order methods), and maybe other numerical stationary states, which coincide with the stationary states of the equations in very few cases.
As soon as external forces enter the modeled system, a source term has to be added to the hyperbolic conservation laws turning these into hyperbolic balance laws. The constant state is not a solution of these anymore. Instead, other static or stationary solutions take its place. However, standard methods are not able to maintain these solutions accurately but introduce discretization errors. This gives rise to the need for the development of socalled wellbalanced methods, i.e. methods which are designed to be exact on special stationary solutions of the system.
In the wellknown shallow water equations with nonflat bottom topography, the most widely considered static state, which is the lakeatrest solution, can be formulated in a closed form. This favors the construction of wellbalanced methods for this system. There is a rich literature about wellbalanced methods for shallow water equations (Audusse et al. (2004); Bermudez and Vázquez (1994); LeVeque (1998) and references therein) and related systems like the Ripa model (Desveaux et al. (2016); Touma and Klingenberg (2015) and references therein). This includes higher order methods for static Noelle et al. (2006) and nonstatic stationary states Noelle et al. (2007). The relevance for methods for nonstatic stationary states has been pointed out in Xing et al. (2011). For tsunami modeling applications high order methods for shallow water equations on nonflat manifolds have been developed e.g. in Castro and Semplice (2018) considering the earth’s surface geometry. For the Euler equations with gravitational potential, on the other hand, static solutions have to be found by solving a differential equation for density and pressure. The second equation relating these quantities is the equation of state (EoS). This makes the construction of wellbalanced methods much more delicate and typically restricts the resulting method to special cases. Many methods have been developed for some classes of hydrostatic states assuming an ideal gas EoS. Examples are given in Chandrashekar and Klingenberg (2015); Cargo and LeRoux (1994); Chertock et al. (2018); Ghosh and Constantinescu (2016); LeVeque and Bale (1999); LeVeque et al. (2011); Touma et al. (2016). There are also higher order methods, see e.g. Klingenberg et al. (2018); Xing and Shu (2013).
While the hydrostatic equation for compressible Euler equations with gravity is basically onedimensional, the spatial structure of this relation is much richer for compressible ideal magnetohydrodynamics (MHD) equations with gravity since it includes offdiagonal terms. In Fuchs et al. (2010) a wellbalanced method for MHD is derived to compute waves on the stationary background. This method is designed to balance isothermal hydrostatic states of the Euler equations together with a magnetic field, which satisfies certain stationarity conditions and is known a priori. Part of this method, namely considering deviations to a background magnetic field, goes back to to Tanaka Tanaka (1994) and is also used by Powell et al. Powell et al. (1999). To do so, the background magnetic field is assumed to be static and free of as well rotation as divergence.
There are different approaches to obtain the wellbalanced property. Some methods are based on a relaxation approach, in which the hydrostatic equation is included in the relaxation system Desveaux et al. (2014, 2016); Thomann et al. (2018). Another widespread idea is the hydrostatic reconstruction, i.e. reconstruction of variables which are constant if the system is in the considered stationary state. An early example of this method is Audusse et al. (2004) for the shallow water system. For Euler equations this approach has been used in Berberich et al. (2018, 2019); Chandrashekar and Klingenberg (2015); Chertock et al. (2018); Ghosh and Constantinescu (2016); Klingenberg et al. (2018); Touma et al. (2016).
The methods for Euler equations mentioned before are restricted to a certain EoS and certain classes of hydrostatic solutions. For astrophysical applications, for example, this restriction is a severe limitation. Equations of state describing physics in the stellar interior are much more complex than the EoS of an ideal gas.
More general methods have been developed in Käppeli and Mishra (2014, 2016); Varma and Chandrashekar (2019); GrosheintzLaval and Käppeli (2019). The wellbalanced methods introduced in these publications can be applied for any EoS. They are exact on certain hydrostatic solutions, in all other cases they are exact on a second order approximation of the considered hydrostatic solution. In Berberich et al. (2019) a second order wellbalanced method for Euler equations with gravity is introduced. This method can be applied for any hydrostatic solution of Euler equations with any EoS if the hydrostatic solution is known. The method is then exact up to machine precision. The method is extended to higher order in Klingenberg et al. (2018). Notably, there is also an extension to stationary states with nonzero velocity in the same article.
The method we present in this paper is designed in the manner of the method in Berberich et al. (2019). It uses the idea of hydrostatic reconstruction and a modification of the source term discretization to obtain the wellbalanced property. The main point in which it differs from all of the methods mentioned above is that our method is not restricted to a certain system of hyperbolic balance laws. Instead, we present a general framework which modifies finite volume methods for any hyperbolic conservation or balance laws such that they obtain the wellbalancing property. Also, the method can be used to balance any solution if it is known before and can be given in an analytical expression or as discrete data. Due to this property we use the phrase wellbalancing in a broader sense than it is typically used in literature. This reference solution, which is chosen to be wellbalanced, can even depend on time, as we will show in this paper.
Our method is also general in the possibility of combination with other modules of a finite volume scheme: It can be applied on any grid system, with any numerical flux function, reconstruction routine, source term discretization, and ODE solver for timediscretization. It allows for higher order in the sense that, if all these components are high order accurate, the resulting method is also high order accurate.
There are several applications in which the hydrostatic solution is given. Consider the application of stellar astrophysics; there, the EoS is often given in the form of a table since there is no analytical expression. Consequently, hydrostatic solutions which depend on the EoS can only be found numerically and given in the form of discrete data. While methods which incorporate analytical expressions are not able to exactly maintain these hydrostatic solutions, it is very well possible with the methods in Berberich et al. (2019); Klingenberg et al. (2018) and the method we present in this paper. Especially, if we consider the better approximation of stellar structure which is given by a stationary state including rotation, our method can be applied to maintain this stationary solytion. Another example from astrophysical application are rotating Keplerian disks. These are twodimensional disks of matter which follows Newton’s laws of motion in the gravitational field of a massive attractor. One way to describe this disk is a stationary solution of Euler equations with gravity including nonzero velocities. Since this is not a hydrostatic solution, conventional wellbalanced methods can not preserve this solution. A special method designed for this application is introduced in Gaburro (2018). In this paper we will show that our method is also able to preserve this solution on different grids. Besides the applicability to any system of hyperbolic balance laws, the balancing of moving and timedepending solutions is one of the key features of our method.
The rest of the paper is structured as follows: In Section 2 we introduce the standard finite volume framework for systems of hyperbolic conservation laws in three spatial dimensions on arbitrary grids. In Section 3, we introduce our general wellbalanced modification for this framework. The wellbalanced property we claim for our method is then shown in Section 4. The validity of the wellbalanced property also depends on a consistent choice of boundary conditions. Therefore, we add a discussion about wellbalanced boundary conditions in Section 6. In Section 7 we add some remarks regarding our method concerning structure, accuracy, and the range of application. To emphasize how simple it is to add our method to an existing finite volume code, we comment on the implementation of the method in Section 8. Finally, in Section 9, we show a variety of numerical tests. These range from applications on Euler equations to ideal magentohydrodynamics (MHD) equations. They include classical wellbalanced tests on the balance laws and also test on the homogeneous hyperbolic conservation laws. Different equations of state are used for the Euler equations. We include a test in which the wellbalanced solution is not analytically known but has been obtained numerically. Also, we present tests in which the wellbalanced solution depends on time. We verify higher order accuracy for solutions close to and far away from the wellbalanced solution numerically. To show the efficiency of the method, we present CPU time comparisons of simulations with and without the wellbalanced modification in Section 10.
2 A standard finite volume method
In this section we present the standard higher order finite volume framework for threedimensional hyperbolic balance laws. There is a rich literature on these methods, e.g. LeVeque (1992); Toro (2009).
Consider the 3d system of hyperbolic balance laws
(1) 
with , where is the flux in direction. Using any mesh we divide the domain into control volumes. For the th control volume () we define the cellaverage
(2) 
where is the control cell volume. Integrating Eq. 1 over and applying the divergence theorem yields the evolution equation for
(3) 
An equivalent formulation which is useful for the following discretization is
(4) 
where is the set of indexes of all control volumes sharing an interface with . For the discretization of the interface fluxes we use a numerical flux function consistent with . The consistency conditions are Lipschitz continuity in the first two arguments and the relation for all normalized vectors . We apply this discretization to Eq. 4 and obtain
(5) 
where the reconstructed functions are obtained using a consistent conservative reconstruction routine on the cell average values . Examples for popular consistent conservative reconstruction routines can be found in Harten et al. (1987); Liu et al. (1994); Levy et al. (2000); Toro (2009) and simple methods in the appendix. In the next step we use numerical quadrature rules for the interface flux integral and a discretization of the source term integral. The semidiscrete method is then
(6) 
is the number of quadrature points at the interfaces, are the quadrature points at the interface and are the corresponding weights. The symbol denotes a consistent discretization of the integral over the argument in the domain . The quadrature rules and source term discretizations we use in our tests are given in the appendix.
The semidiscrete scheme Eq. 6 is th order accurate if the applied reconstruction routine, interface flux quadrature and source term discretization are all at least th order accurate. It can then be evolved in time using a th order accurate ODE solver to obtain a th order accurate fully discrete scheme.
3 The wellbalanced modification of the standard finite volume method
In this section we will introduce a wellbalanced modification for the threedimensional finite volume method presented in Section 2. Reducing it to one or two spatial dimensions is straight forward.
Let be a given continuous solution of Eq. 1. Plugging this reference solution into Eq. 3 we get
(7) 
where is the average of the reference solution in the th control volume. In the next step we subtract Eq. 7 from Eq. 3 to obtain
(8) 
Now, let us rewrite Eq. 8 in terms of the deviation from the reference solution
(9) 
This yields
(10) 
where is the set of indexes of all control volumes sharing an interface with . At this point, we start to discretize. For that we define a numerical flux difference approximation
(11) 
where is a numerical flux function consistent with . We apply this discretization to Eq. 10 and obtain
(12) 
where the reconstructed functions are obtained using a consistent conservative reconstruction routine on the cell average values . In the next step we use numerical quadrature rules for the interface flux integral and a discretization of the source term integral. The semidiscrete method is then
(13) 
where the notations used are as in Eq. 6.
As in the standard method, this semidiscrete scheme Eq. 13 is th order accurate if the applied reconstruction routine, interface flux quadrature and source term discretization are all at least th order accurate. It can then be evolved in time using a th order accurate ODE solver to obtain a th order accurate fully discrete scheme.
4 Proof of the wellbalanced property
In this section we show the wellbalanced property of our method.
Theorem 4.1.
The modified finite volume method introduced in Section 3 satisfies the following property: If
(14) 
at initial time, then this holds for all .
Proof: Let for all . The consistency of the applied reconstruction leads to at all flux quadrature points. The flux consistency then yields
(15) 
Now, consider the contribution from the source term: With the source term discretization in Eq. 13 reduces to
(16) 
We have shown that the right hand side in Eq. 13 vanishes and thus the initial data are conserved for all time. ∎
In creftype 4.1 the formulation of the wellbalanced property is quite simple and maybe not intuitive. If we formulate the result in terms of the actual solution, it might read like this:
Corollary 4.2.
If the initial condition , , equals the cell averages of the reference solution , , the computed solution equals the reference solution for all time.
5 The treatment of discrete reference solutions
It is not necessary that the reference solution in the wellbalanced modification introduced in Section 3 has to be known analytically. It can also be given in the form of discrete data. For consistency with the system of balance laws Eq. 1 it is important that the discrete data which are used for the reference solution converge to a solution of Eq. 1 when the computational grid is refined. To ensure the high order of our method, this convergence should also be of high order. When only discrete data are given we need a method to compute values at the interface quadrature points and cellaveraged values for the grid on which we use our wellbalanced method. In this section we will describe how this can be done. Depending on the form in which the discrete reference data are given, the smoothness of the data, and the required order of accuracy, different methods for this reconstruction have to be used. Here, we will give two examples.
5.1 Example 1: pointwise 1d data on a fine grid
One application of our method is wellbalancing hydrostatic solutions of Euler equations. Especially in physical applications with complex EoS hydrostatic solutions have to be obtained by numerical methods. Even for multidimensional simulations, the underlying hydrostatic solution can be onedimensional in its nature (see cases (a) and (b) below). Now assume such a numerically approximated hydrostatic solution in one spatial dimension is given in the form of point values , (in conservative variables) on a fine equidistant grid. Assume it is supposed to be used in a twodimensional third order accurate modified finite volume method as introduced in Section 3 on a Cartesian grid. For that, in a first step, we use a cubic spline interpolation to construct a continuous function (e.g. Press et al. (1992)). This function is then extended to two spatial dimensions. How this is done depends on the symmetry of the 2d problem which allowed the reduction to a 1d hydrostatic solution. We consider two different cases:

Suppose we have an essentially 1d hydrostatic solution where gravitational force is at an angle to axis. Then we extend the onedimensional hydrostatic solution to a twodimensional solution via .

Suppose we have a radial hydrostatic solution with a gravity vector pointing towards the center . In that case we extend the hydrostatic solution back to 2d by setting .
The values for the reference solution at interface quadrature points can then be evaluated pointwise as . The cell average values of the reference solution are computed using the third order accurate 2d quadrature rule Eq. 46. Method (a) is used in the test in Section 9.1.3.
5.2 Example 2: reference solution from a highly resolved finite volume simulation
Consider a reference solution given as numerical solution of a finite volume simulation on a twodimensional structured static grid with curvilinear coordinates as described in A.1. In this case the data are given as cell averages , instead of point values. Assume this reference solution is supposed to be used in a second order accurate twodimensional modified finite volume method as described in Section 3. The grid to be used is a coarser version of the grid on which the reference solution has been computed, such that all interfaces on the coarse grid coincide with interfaces of the fine grid. For each time step of the stored reference data we map the fine grid on the coarse grid with
(17) 
where all quantities with correspond to the coarse grid and all quantities with correspond to the fine grid. The values of at intermediate times are obtained via interpolation in time. Note, that the interpolation has to be conservative in the sense that for any if . This has to hold for any subset of the set of grid indexes. We choose linear interpolation in time, i.e.
(18) 
which satisfies this property. Note that this conservation property poses a severe restriction for higher order in time interpolation. For a second order method only one interface quadrature point at the center of the interface is needed. For that we can in the computational coordinates just interpolate linearly from the cellaverages interpreted as second order approximation to the cellcentered values. Thus we have and .
6 Boundary conditions
In the previous sections (including the proof of the wellbalanced property) we omitted to include boundary conditions in the discussion. Yet, the validity of the wellbalanced property also depends on the correct choice of boundary conditions. In this section we will describe boundary conditions which have all of the following properties:

They are compatible with the wellbalancing property.

They support the potentially high order of the scheme.

They have relevance for actual application.
Some of the proposed numerical boundary conditions require knowledge of the reference solution outside the domain. If this is not given, one can simply extrapolate the reference solution to the ghost cells. This will not affect the wellbalanced property nor order of accuracy, if the extrapolation is done with a sufficiently high order.
Extrapolation boundary conditions: One way to treat boundaries is the extrapolation of data in the domain to ghost cells. This can be done with high order to support the high order of the applied scheme. In our method we extrapolate the deviations . In the case that holds for all control volumes in the domain, this will also be true for the extrapolated states. Hence, the wellbalanced property also holds at the boundary. Extrapolation boundary conditions for one and two spatial dimensions with different orders of accuracy can be found in A.2.
Wall boundary conditions: To mimic the effect of a solid wall in a finite volume simulation, the values in the cells next to the boundary are copied to ghost cells outside the boundary. The normal velocity component switches the sign. Assuming that the reference solution is compliant to the boundary conditions, we can simply apply wall boundary conditions to the deviation average states in our wellbalanced method.
Periodic boundary conditions: In the case of our method the deviations to the reference states are communicated to the opposite side of the domain. Periodic boundary conditions are descibed in A.2.
7 Remarks on the wellbalanced method
Remark 7.1.
If a stationary solution is chosen as reference solution (which is the case for classical wellbalancing applications), the time derivative of the reference solution vanishes by definition. This leads to . The described method can then also be used to directly evolve the in time instead of . For that the scheme can be adapted by just substituting all terms with the corresponding terms. The reconstruction has then to be applied on the states .
Remark 7.2.
Note that
(19) 
if the source term in Eq. 1 is linear in the first argument. Due to the linearity of the corresponding source term discretizations, this relation then also holds for the discretized source terms. Naming examples, this is the case for the gravitational source term in Euler or ideal MHD equations and the bottom topography source term in the shallow water equations.
Remark 7.3.
One argument for high order accuracy in the finite volume framework is the polynomial accuracy. We can say that a method is th order, if reconstruction, interface flux quadrature, source term discretization, and boundary condition are exact on th order polynomials. Our method has this property for the deviation to the reference solution. Since the reference solution need not to be a polynomial, this does not translate to . However, because of
(20) 
for any discrete norm the method is also th order accurate in (assuming that as well as are sufficiently smooth). The notations and denote the cellaverage of the exact analytical solutions opposed to the cellaverages of the numerically approximated solutions and .
The accuracy will also be demonstrated numerically in Section 9. In the following two Remarks (7.4 and 7.5) we discuss the scope of application of our method.
Remark 7.4.
Our wellbalanced method can even be beneficially applied if there is no source term. Applications could include stationary solutions based on vorticity in multidimensional simulations. Corresponding numerical tests are presented in this paper.
Remark 7.5.
We want to emphasize that the reference solution in our wellbalanced method can very well be timedependent. This does not change the consistency, accuracy, or the wellbalanced property formulated in creftype 4.1. Thus, we use the phrase “wellbalancing” in a wider sense than it is typically used.
8 Notes on the implementation
We have seen that our wellbalanced method can be applied for a wide range of problems. In this section we will discuss another property useful in application: The method is also easy to implement. This holds especially if there is an existing finite volume code which shall be modified to obtain a wellbalancing capability. In a typical finite volume code the following changes have to be applied:

Implement a function returning or an array carrying the cell averages and values at flux quadrature points if is timeindependent. In the second case can just be set alongside the initial condition.

After the program received the initial data for the , it should transform to the deviations .

In the routine evaluating the numerical flux, compute instead of . Basically, this just means subtracting the exact flux after evaluation of the numerical flux.

This step is only necessary if the source term is not linear in (see creftype 7.2): Evaluate the source term at the states and . The difference of these source terms is added to the residual.

Do not forget to correct the output routine such that the states are written to the output.
Let us remind of creftype 7.1 and point out that an alternative implementation could also evolve instead of in time if is timeindependent. In some codes this can be the easier way.
9 Numerical tests of the scheme
We test the method we proposed using a simple python finite volume code. It is build in a modular way, such that different methods can be applied.
Grids: Twodimensional discretization is realized using a structured grid. In some tests we use curvilinear grids. The implementation of these grids restricts the total method to second order accuracy. Note that this is not a general statement. This restriction is due to the special implementation of the grid in our code. Higher order can only be achieved with a Cartesian grid in our code.
Numerical flux function: As numerical flux function we use the Lax–Friedrichs flux LeVeque (1992), since it is simple and can be applied for any hyperbolic system. In some tests we use the Roe’s approximate Riemann solver for Euler equations Roe (1981) to obtain more accurate results.
First order method: To formally obtain a first order method we use constant reconstruction to obtain the interface values. The numerical fluxes are computed at the center of the interfaces. The source term is evaluated at the cellcenter. For the gravity source term used in our tests, we need the gradient of the given gravitational potential. This gradient is computed exactly at the cellcenter.
Second order method: To formally obtain a second order method we use a conservative linear reconstruction to obtain the interface values. This is the only difference to the first order method
Third order method: To formally obtain a third order method we use a conservative parabolic reconstruction to obtain the interface values. In the twodimensional case this is realized using the polynomial from Levy et al. (2000). In the twodimensional case, the numerical fluxes are evaluated at the Gauß–Legendre quadrature points. The source term is evaluated by exactly integrating a parabolic interpolation polynomial. We interpolate from cellcentered states which are obtained via conservative parabolic reconstruction from the cellaverage states.
Seventh order method: To formally obtain a onedimensional seventh order method we use a conservative polynomial reconstruction to obtain the interface and cellcentered values. On the cellcentered values multiplied with the pointwise gravity acceleration we interpolate the source term using a polynomial and integrate exactly.
Boundary conditions: If the setup has periodic character we use periodic boundary conditions. Otherwise we extrapolate the states to ghost cells with a sufficiently high spatial order. If we use the third order method, for example, we extrapolate using parabolas.
Timestepping: We evolve the semidiscrete scheme in time using an explicit third order four stage Runge–Kutta method from Kraaijevanger (1991).
For the seventh order method we use the explicit tenth order 17 stage Runge–Kutta method from Feagin (2007).
For brevity we only give a short description of the methods in this place. The interested reader can find details of the methods in A. Note, that we use simple and straight forward methods for the tests in this section. However, our wellbalanced method is not restricted to these. Instead, one can also use more sophisticated methods, like e.g. WENOtype methods for reconstruction or numerical flux functions designed for special problems (e.g. a low Mach number compliant method for Euler equations like in Barsukow et al. (2017)).
9.1 Euler equations with gravity
The 2d compressible Euler equations which model the balance laws of mass, momentum, and energy under the influence of gravity are given by
(21) 
where the conserved variables, fluxes and source terms are
(22) 
with . Moreover, is the total energy per unit volume with the velocity and specific internal energy . The scalar function is a given gravitational potential. An additional relation between density, pressure, and specific internal energy is given in the form of an equation of state (EoS). In our tests we will use the ideal gas EoS
(23) 
with , although our wellbalanced method can be applied for Euler equations with any EoS.
The 2d Euler equations can be reduced to 1d Euler equations by setting and removing the equation. It can be reduced to homogeneous Euler equations by setting .
9.1.1 1d isothermal hydrostatic solution
We consider an isothermal hydrostatic solution of the 1d compressible Euler equations with gravitational source term and the ideal gas equation of state given by
(24) 
We set these data on a 1d grid with 128 grid cells on the domain . These initial data are evolved up to the final time using the first, second, and third order method with the standard method and the wellbalanced method each. In the wellbalanced method we set the initial data Eq. 24 as timeindependent reference solution. The errors at final time compared to the initial grid can be seen in Table 1. We see that there is no error when the wellbalanced method is applied.
Method  error  error  error 

standard, first order  1.19156249e01  2.17623832e02  1.63898613e01 
wellbalanced, first order  0.00000000e+00  0.00000000e+00  0.00000000e+00 
standard, second order  5.38535901e04  3.05814377e05  7.59110453e04 
wellbalanced, second order  0.00000000e+00  0.00000000e+00  0.00000000e+00 
standard, third order  3.57840305e04  2.24036794e05  4.93235762e04 
wellbalanced, third order  0.00000000e+00  0.00000000e+00  0.00000000e+00 
standard, seventh order  1.49009288e08  1.22670251e10  1.49174239e08 
wellbalanced, seventh order  0.00000000e+00  0.00000000e+00  0.00000000e+00 
Remark 9.1.
To the reader, especially the reader with experience in wellbalanced methods, it might seem unusual and unconvincing that the error is actually zero. Most implementations of wellbalanced methods still produce a machine error, even if the discretization error is eliminated in the wellbalanced method. One reason for that is that typically the source term is balanced against the fluxes. Since the fluxes and source terms are computed in an inherently different way, the difference of these is not exactly zero but some value close to machine precision. In our method, we balance fluxes against fluxes and source term against source term. Thus, the differences can cancel out exactly. If this is the case depends on the specific implementation.
9.1.2 1d isothermal hydrostatic solution with perturbation
We add a perturbation to the pressure such that our initial conditions are
(25) 
in the domain . We choose to test the convergence of our method. We evolve this initial setup up to time using our wellbalanced method (first to third order and seventh order). The results and convergence rates are shown in Table 2. As a reference solution we use a numerical solution computed with the seventh order standard scheme on a grid with cells. All convergence rates match our expectations. The convergence rate for the seventh order scheme drops in the last step, since the error approaches machine precision. In Figs. 2 and 1, density deviations at time for the test with are shown. In Fig. 1 we see, that the discretization error on the hydrostatic background dominates the total error when the second order standard method is used with a low resolution. For higher resolutions, the perturbation can be resolved correctly. In Fig. 2 on the other hand, it gets evident that the second order wellbalanced method is capable of correctly resolving the perturbation on a coarse grid. This is due to the fact that there are no discretization errors on the hydrostatic background as we have already seen in Section 9.1.1.
9.1.3 2d numerically approximated hydrostatic solution
In stellar astrophysical applications, the hydrostatic state of the star can often be given in a discrete form. In this test we will show that our wellbalanced method can be used if the reference solution is given in the form of discrete data in a table.
We assume the following given data: Let the gravitational potential be and the temperature . We assume an ideal gas with radiation pressure, which satisfies the EoS Chandrasekhar (1958)
(26) 
where the temperature is defined implicitly via
(27) 
Using Chebfun Driscoll et al. (2014) in the numerical software MATLAB we solve the 1d hydrostatic equation and EoS for density and pressure corresponding to the given temperature profile. The result is shown in Fig. 3.The data are stored as point values on a fine grid (10,000 data points). The data are set on the 2d grid using the procedure (a) from Section 5.1 We use a grid to evolve the hydrostatic initial condition to the final time . For the conversion between pressure and internal energy we use Newton’s method to solve for the temperature. The errors at final time are shown in Table 3. In all tests using the wellbalanced modification, there is no error at the final time.
Method  error  error  error  error 

standard, first order  4.78e03  1.56e03  1.56e03  5.10e03 
wellbalanced, first order  0.00e+00  0.00e+00  0.00e+00  0.00e+00 
standard, second order  2.63e06  6.19e07  6.19e07  6.34e06 
wellbalanced, second order  0.00e+00  0.00e+00  0.00e+00  0.00e+00 
standard, third order  2.98e06  6.66e08  6.66e08  6.84e06 
wellbalanced, third order  0.00e+00  0.00e+00  0.00e+00  0.00e+00 
9.1.4 Double Gresho vortex
In this test we use a vortex for homogeneous 2d Euler equations first introduced in Gresho and Chan (1990). The pressure and angular velocity of this vortex in dependence of the distance to the center are given by
(28) 
The radial velocity is zero and the density is . In our test we set up the domain with two Gresho vortices centered at and respectively. The vortexes are advected with the velocity and the boundaries are periodic. At time the exact solution of this initial data equals the initial setup. We apply our first order wellbalanced method on a grid to evolve the initial condition up to final time . We use Roe’s numerical flux functions and a linear reconstruction. Only the vortex initially (and finally) centered at is included in the reference solution. The result is illustrated in Fig. 4.
9.1.5 2d Euler wave in gravitational field
method  grid  error  error  error  error 

standard, first order  Cartesian  1.47e02  1.90e02  1.90e02  1.39e01 
wellbalanced, first order  Cartesian  0.00e+00  0.00e+00  0.00e+00  0.00e+00 
standard, second order  Cartesian  7.12e04  7.82e04  7.82e04  2.50e03 
wellbalanced, second order  Cartesian  0.00e+00  0.00e+00  0.00e+00  0.00e+00 
standard, third order  Cartesian  8.68e05  9.98e05  9.98e05  4.50e04 
wellbalanced, third order  Cartesian  0.00e+00  0.00e+00  0.00e+00  0.00e+00 
standard, second order  polar  2.45e04  2.33e04  2.33e04  3.43e04 
wellbalanced, second order  polar  0.00e+00  0.00e+00  0.00e+00  0.00e+00 
To demonstrate that we can follow timedependent solutions exactly with our method we use a problem from Xing and Shu (2013) and Chandrashekar and Zenk (2017) which involves a known exact solution of the 2d Euler equations with gravity given by
(29) 
The gravitational potential is , the EoS is the ideal gas EoS. In accordance to Xing and Shu (2013) and Chandrashekar and Zenk (2017) we choose , on the domain . We use the standard method and the wellbalanced method, first, second, and third order each, to evolve the initial data with to a final time . The grid resolution is on a Cartesian grid. We also include tests with the second order methods on a polar grid with the same grid resolution. The errors at final time are presented in Table 4. If the wellbalanced method is used, there is no error.
9.1.6 Perturbation on the 2d Euler wave in gravitational field
In this test we want to verify the order of accuracy for perturbations to timedependent reference solutions if the wellbalanced method is used. For this we use the initial setup from Eq. 29 and add a pressure perturbation:
(30) 
We evolve these initial data to time using the third order standard and wellbalanced method with . The errors and corresponding convergence rates are presented in Table 5. As reference solution for determining the error we use a numerically approximated solution computed using the third order standard method on a grid. In this test we use exact boundary conditions for the standard method, which means that we evaluate the states in the ghost cells at any time from Eq. 29. We see third order convergence for both methods. However, there seems to be no significant benefit from using the wellbalanced method in this test. The choice of leads to a large discretization error in the perturbation which seems to dominate the total error. This is necessary since we use a solution computed from the standard method as a reference to compute the errors. For smaller perturbations the standard method fails to produce a sufficiently accurate solution. To show the improved accuracy of the wellbalanced modification we add a convergence test with a small perturbation of . To produce a sufficiently accurate reference solution we have to use the third order wellbalanced method on a grid. The errors and convergence rates can be seen in Table 6. Comparing the error of the wellbalanced method on the grid in Table 6 to the error of the third order standard method on the unperturbed solution in Table 4, we see that the wellbalanced method is significantly more accurate.
9.1.7 2d Keplerian disk
Consider a stationary solution given by Gaburro (2018)
(31) 
with the gravitational potential and , , . We use the initial conditions
(32) 
and on the domain . We chose the domain such that we omit the singularity in the velocity at . In Fig. 5 results of numerical tests are illustrated. The second order standard and wellbalanced methods are applied on a polar grid with cells and a Cartesian grid with cells. In the Cartesian grid we take out the center with using Dirichlet boundary conditions. We also use Dirichlet boundary conditions at all outer boundaries. Since there is a discontinuity in the initial setup we apply a minmod slope limiter. The density at time for each simulation is shown in Fig. 5. The exact solution of this problem can be seen in Gaburro (2018). Since this is an purely advective problem and there is no radial component to the velocity, the quantity , which describes the average distance of the density perturbation to the center, is conserved for all time in the exact solution. For our simulations we measure the quantity as a measure of the quality of the numerical solutions. For the exact solution we have for all time. The values of are shown in the captions in Fig. 5.
In the tests with the standard method (Figs. (a)a and (c)c) we see discretization errors in the Keplerian disk solution Eq. 31. This introduces radial velocities, the advection of the spot of increased density has a component towards the center. In the tests using our wellbalanced methods (Figs. (b)b and (d)d), the result is free of discretization errors in the Keplerian disk solution Eq. 31. The advection is more accurate, the only errors are diffusion errors. The polar grid is more suitable for this test problem, since it is adapted to the radial geometry. The test using our wellbalanced method on the Cartesian grid (Fig. (b)b) is more diffusive than the one on the polar grid (Fig. (d)d), yet we see that the wellbalanced modification improves the result significantly.
9.2 Homogeneous compressible ideal magnetohydrodynamics
The 2d compressible ideal magnetohydrodynamics equations which model the conservation of mass, momentum, magnetic field, and energy are given by
(33) 
The conserved variables and fluxes are
(34) 
where , are the  and component of the magnetic field. The total energy is . All other quatities are defined as for the Euler equations. We use the same EoS as for the Euler equations.
One can also define 2d compressible ideal MHD equations such that they include the and components. This is in principle reasonable due to the genuine threedimensional interactions between velocity and magnetic field. In our tests we set and to zero and there is no difference if we omit the corresponding equations.
9.2.1 Stationary Vortex  Long Time
We consider the following exact solution of the homogeneous 2d ideal MHD equations:
(35) 
This setup describes a stationary vortex which is advected through the domain with the velocity . The domain is . One vortex turnovertime is . In a first test we set , and run the test up to on a grid. We use the wellbalanced method and the reference solution equals the initial data. The numerical error at final time compared to the initial setup is exactly zero in all conservative variables.
9.2.2 Stationary vortex  order of accuracy
In a second test with the vortex from Section 9.2.1 we want to see if the wellbalanced method converges as expected, even if the reference solution which is set deviates from the actual solution over time. For that we set , in the initial condition. As reference solution we use the same vortex but with . In Table 7 the errors and rates at final time are presented for the formally first, second, and third order accurate wellbalanced method. We omitted the errors for and . Due to the symmetry of the setup these errors equal the errors in and respectively. We see that even if the reference solution moves away from the actual solution over time the method is consistent with the expected order of accuracy.
10 Computational cost of the modification
Wellbalanced methods are constructed to improve the accuracy with which solutions of balance laws are approximated. In the previous section we have shown that usage of our wellbalanced modification can improve the accuracy of a simulation significantly. On the other hand, an increase of computational effort can countervail the gain in accuracy, if it is too high. In this section we will compare the computation times of simulations using our wellbalanced modification to simulations using the corresponding standard method and show that the increase in CPU time is moderate.
10.1 The procedure
To compare the methods, we will run tests with different setups and grid resolutions using a standard method and the corresponding wellbalanced method. We use the simple python code described in Section 9 on a single CPU. Each test is repeated 20 times and the wall clock times are measured. We compute the average wallclock time and standard deviations of the single runs for every test. We compute and present the ratio of average wall clock time for the wellbalanced compared to the standard method. The errors are computed by adding the relative standard deviations of the wellbalanced and standard method.
10.2 The tests
To compare the runtimes we use test setups described in Section 9. In a first test, we use the perturbed onedimensional isothermal solution described in Section 9.1.2 with a final time . The first, second, third, and seventh order 1d methods are applied. The ratio of wall clock times for the tests with and without the wellbalanced modification can be seen in Fig. 6. Roe’s approximate Riemann solver has been used in the tests. The reference solution used in the wellbalanced modification is constant in time in this test case. It is hence computed once and stored in an array. From the figure we see that we can expect an increase in CPU time of about .
As a second test setup we choose the Keplerian disk from Section 9.1.7. We evolve it to time with a first and second order method on a polar grid. As in the previous test, the reference solution is timeindependent and thus only computed once. We see an increase in CPU time of about when using the wellbalanced method.
To also test the increase of CPU time consumption for a simulation in which the reference solution is time dependent, we use the Euler wave in a gravitational field with perturbation from Section 9.1.6. In this test, the reference solution is computed from a function every time it is used (which happens in every intermediate step). The result of these tests can be seen in Fig. 8. We see an increase in CPU time of about if the wellbalanced method is used.
11 Summary and conclusions
We introduced a new general framework for the construction of wellbalanced finite volume methods for hyperbolic balance laws. A standard finite volume method is modified such that it evolves the deviation from a reference solution instead of the actual solution. This makes the scheme exact on the reference solution. The finite volume method can include any consistent reconstruction, numerical flux function, interface quadrature, source term discretization, and ODE solver for time discretization. Thus it can have arbitrarily high order. The method can be defined on any computational grid geometry. One can view our method as a high order extension of Berberich et al. (2019) and Klingenberg et al. (2018) to all known solutions of all hyperbolic balance laws.
In numerical tests with Euler and MHD equations on different grids we could verify that the method can successfully be applied to exactly maintain static and stationary solutions or even follow timedependent solutions. For that, the solution has to be known either analytically or in the form of discrete data. The latter case is especially interesting for complex applications like stellar astrophysics, where static states of the Euler equations with gravity can be obtained numerically but only in few cases analytically. Also, for the case of differentially rotating stars, our method can be applied for wellbalancing since it can include nonzero velocities in stationary states.
High order accuracy has been confirmed in numerical experiments. For practical applications we emphasize that our wellbalanced method can be easily implemented in existing finite volume codes with minimal effort. Practical advice is given. Also, in a series of numerical tests we have shown that the increase in computational time is moderate. There are applications in which the wellbalanced solution is not known beforehand. In this case, one of the existing wellbalanced methods for the considered balance law has to be applied, if there is one. In all cases, in which the wellbalanced solution is known, our simple frame work can be applied to obtain the wellbalanced property.
Acknowledgments
The research of Jonas Berberich is supported by the Klaus Tschira Foundation. Praveen Chandrashekar would like to acknowledge the support received from Airbus Foundation Chair on Mathematics of Complex Systems at TIFRCAM, Bangalore.
References
 Godunov (1959) S. K. Godunov, Finite difference method for numerical computation of discontinous solution of the equations of fluid dynamics, Matematicheskii Sbornik 47 (1959) 271.
 Audusse et al. (2004) E. Audusse, F. Bouchut, M.O. Bristeau, R. Klein, B. Perthame, A fast and stable wellbalanced scheme with hydrostatic reconstruction for shallow water flows, SIAM Journal on Scientific Computing 25 (2004) 2050–2065.
 Bermudez and Vázquez (1994) A. Bermudez, M. E. Vázquez, Upwind methods for hyperbolic conservation laws with source terms, Computers & Fluids 23 (1994) 1049–1071.
 LeVeque (1998) R. J. LeVeque, Balancing source terms and flux gradients in highresolution godunov methods: the quasisteady wavepropagation algorithm, Journal of computational physics 146 (1998) 346–365.
 Desveaux et al. (2016) V. Desveaux, M. Zenk, C. Berthon, C. Klingenberg, Wellbalanced schemes to capture nonexplicit steady states: Ripa model, Mathematics of Computation 85 (2016) 1571–1602.
 Touma and Klingenberg (2015) R. Touma, C. Klingenberg, Wellbalanced central finite volume methods for the ripa system, Applied Numerical Mathematics 97 (2015) 42–68.
 Noelle et al. (2006) S. Noelle, N. Pankratz, G. Puppo, J. R. Natvig, Wellbalanced finite volume schemes of arbitrary order of accuracy for shallow water flows, Journal of Computational Physics 213 (2006) 474–499.
 Noelle et al. (2007) S. Noelle, Y. Xing, C.W. Shu, Highorder wellbalanced finite volume weno schemes for shallow water equation with moving water, Journal of Computational Physics 226 (2007) 29–58.
 Xing et al. (2011) Y. Xing, C.W. Shu, S. Noelle, On the advantage of wellbalanced schemes for movingwater equilibria of the shallow water equations, Journal of scientific computing 48 (2011) 339–349.
 Castro and Semplice (2018) M. J. Castro, M. Semplice, Thirdand fourthorder wellbalanced schemes for the shallow water equations based on the cweno reconstruction, International Journal for Numerical Methods in Fluids (2018).
 Chandrashekar and Klingenberg (2015) P. Chandrashekar, C. Klingenberg, A second order wellbalanced finite volume scheme for euler equations with gravity, SIAM Journal on Scientific Computing 37 (2015) B382–B402.
 Cargo and LeRoux (1994) P. Cargo, A. LeRoux, A well balanced scheme for a model of atmosphere with gravity, COMPTES RENDUS DE L ACADEMIE DES SCIENCES SERIE IMATHEMATIQUE 318 (1994) 73–76.
 Chertock et al. (2018) A. Chertock, S. Cui, A. Kurganov, Ş. N. Özcan, E. Tadmor, Wellbalanced schemes for the euler equations with gravitation: Conservative formulation using global fluxes, Journal of Computational Physics (2018).
 Ghosh and Constantinescu (2016) D. Ghosh, E. M. Constantinescu, Wellbalanced, conservative finite difference algorithm for atmospheric flows, AIAA Journal (2016).
 LeVeque and Bale (1999) R. LeVeque, D. Bale, Wave propagation methods for conservation laws with source terms., Birkhauser Basel (1999) pp. 609–618.
 LeVeque et al. (2011) R. J. LeVeque, D. L. George, M. J. Berger, Tsunami modelling with adaptively refined finite volume methods, Acta Numerica 20 (2011) 211–289.
 Touma et al. (2016) R. Touma, U. Koley, C. Klingenberg, Wellbalanced unstaggered central schemes for the euler equations with gravitation, SIAM Journal on Scientific Computing 38 (2016) B773–B807.
 Klingenberg et al. (2018) C. Klingenberg, G. Puppo, M. Semplice, Arbitrary order finite volume wellbalanced schemes for the euler equations with gravity, arXiv preprint arXiv:1807.02341 (2018).
 Xing and Shu (2013) Y. Xing, C.W. Shu, High order wellbalanced weno scheme for the gas dynamics equations under gravitational fields, Journal of Scientific Computing 54 (2013) 645–662.
 Fuchs et al. (2010) F. G. Fuchs, A. McMurry, S. Mishra, N. H. Risebro, K. Waagan, High order wellbalanced finite volume schemes for simulating wave propagation in stratified magnetic atmospheres, Journal of Computational Physics 229 (2010) 4033–4058.
 Tanaka (1994) T. Tanaka, Finite volume tvd scheme on an unstructured grid system for threedimensional mhd simulation of inhomogeneous systems including strong background potential fields, Journal of Computational Physics 111 (1994) 381–389.
 Powell et al. (1999) K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, D. L. de Zeeuw, A SolutionAdaptive Upwind Scheme for Ideal Magnetohydrodynamics, Journal of Computational Physics 154 (1999) 284–309.
 Desveaux et al. (2014) V. Desveaux, M. Zenk, C. Berthon, C. Klingenberg, A wellbalanced scheme for the euler equation with a gravitational potential, in: Finite Volumes for Complex Applications VIIMethods and Theoretical Aspects, Springer, 2014, pp. 217–226.
 Desveaux et al. (2016) V. Desveaux, M. Zenk, C. Berthon, C. Klingenberg, A wellbalanced scheme to capture nonexplicit steady states in the euler equations with gravity, International Journal for Numerical Methods in Fluids 81 (2016) 104–127.
 Thomann et al. (2018) A. Thomann, M. Zenk, C. Klingenberg, A second order positivity preserving wellbalanced finite volume scheme for euler equations with gravity for arbitrary hydrostatic equilibria, arXiv preprint arXiv:1804.09965 (2018).
 Berberich et al. (2018) J. P. Berberich, P. Chandrashekar, C. Klingenberg, A general wellbalanced finite volume scheme for euler equations with gravity, in: C. Klingenberg, M. Westdickenberg (Eds.), Theory, Numerics and Applications of Hyperbolic Problems I, Springer Proceedings in Mathematics & Statistics 236, 2018.
 Berberich et al. (2019) J. P. Berberich, P. Chandrashekar, C. Klingenberg, F. K. Röpke, Second order finite volume scheme for euler equations with gravity which is wellbalanced for general equations of state and grid systems, accepted for publication in Communications in Computational Physics (2019).
 Käppeli and Mishra (2014) R. Käppeli, S. Mishra, Wellbalanced schemes for the euler equations with gravitation, Journal of Computational Physics 259 (2014) 199–219.
 Käppeli and Mishra (2016) R. Käppeli, S. Mishra, A wellbalanced finite volume scheme for the euler equations with gravitationthe exact preservation of hydrostatic equilibrium with arbitrary entropy stratification, Astronomy & Astrophysics 587 (2016) A94.
 Varma and Chandrashekar (2019) D. Varma, P. Chandrashekar, A secondorder, discretely wellbalanced finite volume scheme for euler equations with gravity, Computers & Fluids (2019).
 GrosheintzLaval and Käppeli (2019) L. GrosheintzLaval, R. Käppeli, Highorder wellbalanced finite volume schemes for the euler equations with gravitation, Journal of Computational Physics 378 (2019) 324–343.
 Gaburro (2018) E. Gaburro, Well balanced ArbitraryLagrangianEulerian Finite Volume schemes on moving nonconforming meshes for nonconservative Hyperbolic systems, Ph.D. thesis, University of Trento, 2018.
 LeVeque (1992) R. J. LeVeque, Numerical Methods for Conservation Laws, 2 ed., BirkhÃ¤user, Basel, 1992.
 Toro (2009) E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, Springer, Berlin Heidelberg, 2009.
 Harten et al. (1987) A. Harten, B. Engquist, S. Osher, S. R. Chakravarthy, Uniformly high order accurate essentially nonoscillatory schemes, iii, in: Upwind and highresolution schemes, Springer, 1987, pp. 218–290.
 Liu et al. (1994) X.D. Liu, S. Osher, T. Chan, Weighted Essentially Nonoscillatory Schemes, Journal of Computational Physics 115 (1994) 200–212.
 Levy et al. (2000) D. Levy, G. Puppo, G. Russo, Compact central weno schemes for multidimensional conservation laws, SIAM Journal on Scientific Computing 22 (2000) 656–672.
 Press et al. (1992) W. H. Press, S. Teukolsky, W. Vetterling, B. Flannery, Numerical Recipes in C: the art of scientific computing, Second Edition, Cambridge Univ. Press, New York, 1992.
 Roe (1981) P. L. Roe, Approximate riemann solvers, parameter vectors, and difference schemes, Journal of Computational Physics 43 (1981) 357 – 372.
 Kraaijevanger (1991) J. F. B. M. Kraaijevanger, Contractivity of rungekutta methods, BIT Numerical Mathematics 31 (1991) 482–528.
 Feagin (2007) T. Feagin, A tenthorder rungekutta method with error estimate, in: Proceedings of the IAENG Conference on Scientific Computing, 2007.
 Barsukow et al. (2017) W. Barsukow, P. V. F. Edelmann, C. Klingenberg, F. Miczek, F. K. Röpke, A numerical scheme for the compressible lowmach number regime of ideal fluid dynamicsroe, Journal of Scientific Computing (2017) 1–24.
 Chandrasekhar (1958) S. Chandrasekhar, An introduction to the study of stellar structure, volume 2, Courier Corporation, 1958.
 Driscoll et al. (2014) T. A. Driscoll, N. Hale, L. N. Trefethen, Chebfun Guide, Pafnuty Publications, Oxford, 2014.
 Gresho and Chan (1990) P. M. Gresho, S. T. Chan, On the theory of semiimplicit projection methods for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix. part 2: Implementation, International Journal for Numerical Methods in Fluids 11 (1990) 621–659.
 Chandrashekar and Zenk (2017) P. Chandrashekar, M. Zenk, Wellbalanced nodal discontinuous galerkin method for euler equations with gravity, Journal of Scientific Computing (2017) 1–32.
Appendix A Details of the applied finite volume schemes
In all numerical tests we use structured grids. Hence, in the description of the details we restrict to structured grids. Some parts of the scheme, such as the reconstruction methods, are applied to in the standard method, to in the wellbalanced method. We will denote the states with . Depending on the method we have or . Analogously, we will use to denote or .
a.1 Curvilinear coordinates
We define a 2d curvilinear coordinate system. The coordinates in physical space are , the coordinates in computational space are . The th cell is denoted in the physical space and by in the computational space. We can rewrite Eq. 21 in the computational coordinates as
(36) 
where