High order well-balanced finite volume methods for multi-dimensional systems of hyperbolic balance laws
We introduce a general framework for the construction of well-balanced finite volume methods for hyperbolic balance laws. The phrase well-balancing 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 well-balancing 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:finite-volume methods, well-balancing, hyperbolic balance laws, compressible Euler equations with gravity, ideal magnetohydrodynamics
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 so-called well-balanced methods, i.e. methods which are designed to be exact on special stationary solutions of the system.
In the well-known shallow water equations with non-flat bottom topography, the most widely considered static state, which is the lake-at-rest solution, can be formulated in a closed form. This favors the construction of well-balanced methods for this system. There is a rich literature about well-balanced 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 non-static stationary states Noelle et al. (2007). The relevance for methods for non-static stationary states has been pointed out in Xing et al. (2011). For tsunami modeling applications high order methods for shallow water equations on non-flat 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 well-balanced 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 one-dimensional, the spatial structure of this relation is much richer for compressible ideal magnetohydrodynamics (MHD) equations with gravity since it includes off-diagonal terms. In Fuchs et al. (2010) a well-balanced 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 well-balanced 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); Grosheintz-Laval and Käppeli (2019). The well-balanced 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 well-balanced 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 non-zero 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 well-balanced 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 well-balancing 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 well-balancing in a broader sense than it is typically used in literature. This reference solution, which is chosen to be well-balanced, 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 time-discretization. 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 two-dimensional 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 non-zero velocities. Since this is not a hydrostatic solution, conventional well-balanced 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 time-depending 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 well-balanced modification for this framework. The well-balanced property we claim for our method is then shown in Section 4. The validity of the well-balanced property also depends on a consistent choice of boundary conditions. Therefore, we add a discussion about well-balanced 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 well-balanced 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 well-balanced solution is not analytically known but has been obtained numerically. Also, we present tests in which the well-balanced solution depends on time. We verify higher order accuracy for solutions close to and far away from the well-balanced solution numerically. To show the efficiency of the method, we present CPU time comparisons of simulations with and without the well-balanced modification in Section 10.
2 A standard finite volume method
In this section we present the standard higher order finite volume framework for three-dimensional hyperbolic balance laws. There is a rich literature on these methods, e.g. LeVeque (1992); Toro (2009).
Consider the 3-d system of hyperbolic balance laws
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 cell-average
where is the control cell volume. Integrating Eq. 1 over and applying the divergence theorem yields the evolution equation for
An equivalent formulation which is useful for the following discretization is
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
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 semi-discrete method is then
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 semi-discrete 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 well-balanced modification of the standard finite volume method
In this section we will introduce a well-balanced modification for the three-dimensional finite volume method presented in Section 2. Reducing it to one or two spatial dimensions is straight forward.
Now, let us rewrite Eq. 8 in terms of the deviation from the reference solution
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
where is a numerical flux function consistent with . We apply this discretization to Eq. 10 and obtain
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 semi-discrete method is then
where the notations used are as in Eq. 6.
As in the standard method, this semi-discrete 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 well-balanced property
In this section we show the well-balanced property of our method.
The modified finite volume method introduced in Section 3 satisfies the following property: If
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
Now, consider the contribution from the source term: With the source term discretization in Eq. 13 reduces to
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 well-balanced property is quite simple and maybe not intuitive. If we formulate the result in terms of the actual solution, it might read like this:
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 well-balanced 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 cell-averaged values for the grid on which we use our well-balanced 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 1-d data on a fine grid
One application of our method is well-balancing hydrostatic solutions of Euler equations. Especially in physical applications with complex EoS hydrostatic solutions have to be obtained by numerical methods. Even for multi-dimensional simulations, the underlying hydrostatic solution can be one-dimensional 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 two-dimensional 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 2-d problem which allowed the reduction to a 1-d hydrostatic solution. We consider two different cases:
Suppose we have an essentially 1-d hydrostatic solution where gravitational force is at an angle to -axis. Then we extend the one-dimensional hydrostatic solution to a two-dimensional 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 2-d 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 2-d 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 two-dimensional 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 two-dimensional 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
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.
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 cell-averages interpreted as second order approximation to the cell-centered values. Thus we have and .
6 Boundary conditions
In the previous sections (including the proof of the well-balanced property) we omitted to include boundary conditions in the discussion. Yet, the validity of the well-balanced 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 well-balancing 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 well-balanced 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 well-balanced 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 well-balanced 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 well-balanced method
If a stationary solution is chosen as reference solution (which is the case for classical well-balancing 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 .
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.
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
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 cell-average of the exact analytical solutions opposed to the cell-averages of the numerically approximated solutions and .
Our well-balanced method can even be beneficially applied if there is no source term. Applications could include stationary solutions based on vorticity in multi-dimensional simulations. Corresponding numerical tests are presented in this paper.
We want to emphasize that the reference solution in our well-balanced method can very well be time-dependent. This does not change the consistency, accuracy, or the well-balanced property formulated in creftype 4.1. Thus, we use the phrase “well-balancing” in a wider sense than it is typically used.
8 Notes on the implementation
We have seen that our well-balanced 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 well-balancing 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 time-independent. 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 time-independent. 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: Two-dimensional 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 cell-center. 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 cell-center.
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 two-dimensional case this is realized using the polynomial from Levy et al. (2000). In the two-dimensional 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 cell-centered states which are obtained via conservative parabolic reconstruction from the cell-average states.
Seventh order method: To formally obtain a one-dimensional seventh order method we use a conservative polynomial reconstruction to obtain the interface and cell-centered values. On the cell-centered values multiplied with the point-wise 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.
Time-stepping: We evolve the semi-discrete 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 well-balanced method is not restricted to these. Instead, one can also use more sophisticated methods, like e.g. WENO-type 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 2-d compressible Euler equations which model the balance laws of mass, momentum, and energy under the influence of gravity are given by
where the conserved variables, fluxes and source terms are
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
with , although our well-balanced method can be applied for Euler equations with any EoS.
The 2-d Euler equations can be reduced to 1-d Euler equations by setting and removing the equation. It can be reduced to homogeneous Euler equations by setting .
9.1.1 1-d isothermal hydrostatic solution
We consider an isothermal hydrostatic solution of the 1-d compressible Euler equations with gravitational source term and the ideal gas equation of state given by
We set these data on a 1-d 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 well-balanced method each. In the well-balanced method we set the initial data Eq. 24 as time-independent 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 well-balanced method is applied.
|standard, first order||1.19156249e-01||2.17623832e-02||1.63898613e-01|
|well-balanced, first order||0.00000000e+00||0.00000000e+00||0.00000000e+00|
|standard, second order||5.38535901e-04||3.05814377e-05||7.59110453e-04|
|well-balanced, second order||0.00000000e+00||0.00000000e+00||0.00000000e+00|
|standard, third order||3.57840305e-04||2.24036794e-05||4.93235762e-04|
|well-balanced, third order||0.00000000e+00||0.00000000e+00||0.00000000e+00|
|standard, seventh order||1.49009288e-08||1.22670251e-10||1.49174239e-08|
|well-balanced, seventh order||0.00000000e+00||0.00000000e+00||0.00000000e+00|
To the reader, especially the reader with experience in well-balanced methods, it might seem unusual and unconvincing that the error is actually zero. Most implementations of well-balanced methods still produce a machine error, even if the discretization error is eliminated in the well-balanced 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 1-d isothermal hydrostatic solution with perturbation
We add a perturbation to the pressure such that our initial conditions are
in the domain . We choose to test the convergence of our method. We evolve this initial setup up to time using our well-balanced 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 well-balanced 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 2-d 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 well-balanced 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)
where the temperature is defined implicitly via
Using Chebfun Driscoll et al. (2014) in the numerical software MATLAB we solve the 1-d 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 2-d 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 well-balanced modification, there is no error at the final time.
|standard, first order||4.78e-03||1.56e-03||1.56e-03||5.10e-03|
|well-balanced, first order||0.00e+00||0.00e+00||0.00e+00||0.00e+00|
|standard, second order||2.63e-06||6.19e-07||6.19e-07||6.34e-06|
|well-balanced, second order||0.00e+00||0.00e+00||0.00e+00||0.00e+00|
|standard, third order||2.98e-06||6.66e-08||6.66e-08||6.84e-06|
|well-balanced, 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 2-d 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
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 well-balanced 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 2-d Euler wave in gravitational field
|standard, first order||Cartesian||1.47e-02||1.90e-02||1.90e-02||1.39e-01|
|well-balanced, first order||Cartesian||0.00e+00||0.00e+00||0.00e+00||0.00e+00|
|standard, second order||Cartesian||7.12e-04||7.82e-04||7.82e-04||2.50e-03|
|well-balanced, second order||Cartesian||0.00e+00||0.00e+00||0.00e+00||0.00e+00|
|standard, third order||Cartesian||8.68e-05||9.98e-05||9.98e-05||4.50e-04|
|well-balanced, third order||Cartesian||0.00e+00||0.00e+00||0.00e+00||0.00e+00|
|standard, second order||polar||2.45e-04||2.33e-04||2.33e-04||3.43e-04|
|well-balanced, second order||polar||0.00e+00||0.00e+00||0.00e+00||0.00e+00|
To demonstrate that we can follow time-dependent 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 2-d Euler equations with gravity given by
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 well-balanced 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 well-balanced method is used, there is no error.
9.1.6 Perturbation on the 2-d Euler wave in gravitational field
In this test we want to verify the order of accuracy for perturbations to time-dependent reference solutions if the well-balanced method is used. For this we use the initial setup from Eq. 29 and add a pressure perturbation:
We evolve these initial data to time using the third order standard and well-balanced 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 well-balanced 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 well-balanced 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 well-balanced method on a grid. The errors and convergence rates can be seen in Table 6. Comparing the error of the well-balanced 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 well-balanced method is significantly more accurate.
9.1.7 2-d Keplerian disk
Consider a stationary solution given by Gaburro (2018)
with the gravitational potential and , , . We use the initial conditions
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 well-balanced 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 well-balanced 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 well-balanced 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 well-balanced modification improves the result significantly.
9.2 Homogeneous compressible ideal magnetohydrodynamics
The 2-d compressible ideal magnetohydrodynamics equations which model the conservation of mass, momentum, magnetic field, and energy are given by
The conserved variables and fluxes are
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 2-d compressible ideal MHD equations such that they include the and components. This is in principle reasonable due to the genuine three-dimensional 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 2-d ideal MHD equations:
This setup describes a stationary vortex which is advected through the domain with the velocity . The domain is . One vortex turnover-time is . In a first test we set , and run the test up to on a grid. We use the well-balanced 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 well-balanced 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 well-balanced 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
Well-balanced 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 well-balanced 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 well-balanced 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 well-balanced 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 wall-clock time and standard deviations of the single runs for every test. We compute and present the ratio of average wall clock time for the well-balanced compared to the standard method. The errors are computed by adding the relative standard deviations of the well-balanced 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 one-dimensional isothermal solution described in Section 9.1.2 with a final time . The first, second, third, and seventh order 1-d methods are applied. The ratio of wall clock times for the tests with and without the well-balanced modification can be seen in Fig. 6. Roe’s approximate Riemann solver has been used in the tests. The reference solution used in the well-balanced 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 time-independent and thus only computed once. We see an increase in CPU time of about when using the well-balanced 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 well-balanced method is used.
11 Summary and conclusions
We introduced a new general framework for the construction of well-balanced 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 time-dependent 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 well-balancing since it can include non-zero velocities in stationary states.
High order accuracy has been confirmed in numerical experiments. For practical applications we emphasize that our well-balanced 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 well-balanced solution is not known beforehand. In this case, one of the existing well-balanced methods for the considered balance law has to be applied, if there is one. In all cases, in which the well-balanced solution is known, our simple frame work can be applied to obtain the well-balanced property.
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 TIFR-CAM, Bangalore.
- 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 well-balanced 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 high-resolution godunov methods: the quasi-steady wave-propagation algorithm, Journal of computational physics 146 (1998) 346–365.
- Desveaux et al. (2016) V. Desveaux, M. Zenk, C. Berthon, C. Klingenberg, Well-balanced schemes to capture non-explicit steady states: Ripa model, Mathematics of Computation 85 (2016) 1571–1602.
- Touma and Klingenberg (2015) R. Touma, C. Klingenberg, Well-balanced 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, Well-balanced 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, High-order well-balanced 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 well-balanced schemes for moving-water equilibria of the shallow water equations, Journal of scientific computing 48 (2011) 339–349.
- Castro and Semplice (2018) M. J. Castro, M. Semplice, Third-and fourth-order well-balanced 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 well-balanced 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 I-MATHEMATIQUE 318 (1994) 73–76.
- Chertock et al. (2018) A. Chertock, S. Cui, A. Kurganov, Ş. N. Özcan, E. Tadmor, Well-balanced 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, Well-balanced, 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, Well-balanced 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 well-balanced schemes for the euler equations with gravity, arXiv preprint arXiv:1807.02341 (2018).
- Xing and Shu (2013) Y. Xing, C.-W. Shu, High order well-balanced 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 well-balanced 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 three-dimensional 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 Solution-Adaptive 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 well-balanced scheme for the euler equation with a gravitational potential, in: Finite Volumes for Complex Applications VII-Methods and Theoretical Aspects, Springer, 2014, pp. 217–226.
- Desveaux et al. (2016) V. Desveaux, M. Zenk, C. Berthon, C. Klingenberg, A well-balanced scheme to capture non-explicit 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 well-balanced 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 well-balanced 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 well-balanced 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, Well-balanced 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 well-balanced finite volume scheme for the euler equations with gravitation-the exact preservation of hydrostatic equilibrium with arbitrary entropy stratification, Astronomy & Astrophysics 587 (2016) A94.
- Varma and Chandrashekar (2019) D. Varma, P. Chandrashekar, A second-order, discretely well-balanced finite volume scheme for euler equations with gravity, Computers & Fluids (2019).
- Grosheintz-Laval and Käppeli (2019) L. Grosheintz-Laval, R. Käppeli, High-order well-balanced finite volume schemes for the euler equations with gravitation, Journal of Computational Physics 378 (2019) 324–343.
- Gaburro (2018) E. Gaburro, Well balanced Arbitrary-Lagrangian-Eulerian Finite Volume schemes on moving nonconforming meshes for non-conservative 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 non-oscillatory schemes, iii, in: Upwind and high-resolution schemes, Springer, 1987, pp. 218–290.
- Liu et al. (1994) X.-D. Liu, S. Osher, T. Chan, Weighted Essentially Non-oscillatory 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 runge-kutta methods, BIT Numerical Mathematics 31 (1991) 482–528.
- Feagin (2007) T. Feagin, A tenth-order runge-kutta 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 low-mach 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 semi-implicit 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, Well-balanced 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 well-balanced 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 2-d 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