High order well-balanced finite volume methods for multi-dimensional systems of hyperbolic balance laws

TL;DR

This paper presents a flexible finite volume framework that ensures exact preservation of known solutions for multi-dimensional hyperbolic balance laws, maintaining high order accuracy and applicable to various physical systems.

Contribution

A novel general framework for constructing well-balanced finite volume methods that exactly follow known solutions in multiple dimensions, applicable to diverse hyperbolic systems.

Findings

Successfully applied to compressible Euler equations with gravity

Maintains high order accuracy while being well-balanced

Effective for magnetohydrodynamics equations

Abstract

We introduce a general framework for the construction of well-balanced finite volume methods for hyperbolic balance laws. We use the phrase well-balancing in a broader sense, since our proposed method can be applied to exactly follow any solution of any system of hyperbolic balance laws in multiple spatial dimensions and not only time independent solutions. 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 and in case the method is high order accurate, this is maintained under our modification. We present 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.

Tables

Table 1. Table 1: L 1 superscript 𝐿 1 L^{1} -errors for an isothermal hydrostatic solution of the Euler equations after time t = 2 𝑡 2 t=2 computed using the standard (Std) and well-balanced (WB) methods with different orders of accuracy (O m 𝑚 m for m 𝑚 m -th order). The setup is described in Section   8.3 .

error in	Std-O1	WB-O1	Std-O2	WB-O2	Std-O3	WB-O3	Std-O7	WB-O7
$\rho$	1.19e-01	0.00e+00	4.60e-04	0.00e+00	9.72e-05	0.00e+00	1.28e-09	0.00e+00
$\rho u$	2.18e-02	0.00e+00	7.24e-04	0.00e+00	1.50e-04	0.00e+00	1.40e-09	0.00e+00
$E$	1.64e-01	0.00e+00	3.20e-03	0.00e+00	3.92e-04	0.00e+00	3.99e-09	0.00e+00

Table 2. Table 2: L 1 superscript 𝐿 1 L^{1} -errors and convergence rates in total energy for a small pressure perturbation on an isothermal hydrostatic solution of the Euler equations after time t = 0.2 𝑡 0.2 t=0.2 . Different well-balanced methods are used. The setup is described in Section   8.4 .

$N$	WB-O1		WB-O2		WB-O3		WB-O7
	$E$ error	rate	$E$ error	rate	$E$ error	rate	$E$ error	rate
256	5.73e-03	0.9	5.98e-05	2.0	6.93e-05	2.6	1.32e-09	6.5
512	3.08e-03	0.9	1.49e-05	2.0	1.18e-05	2.7	1.46e-11	6.6
1024	1.60e-03	1.0	3.73e-06	2.0	1.78e-06	2.8	1.46e-13	4.1
2048	8.15e-04		9.36e-07		2.53e-07		8.41e-15

Table 3. Table 3: L 1 superscript 𝐿 1 L^{1} -errors in total energy for the numerically integrated hydrostatic solution of the Euler equations with radiation pressure after time t = 2 𝑡 2 t=2 computed using different methods. The setup is described in Section   8.6 .

$N$	Std-O1	WB-O1	Std-O2	WB-O2	Std-O3	WB-O3
64	5.10e-03	0.00e+00	4.38e-05	0.00e+00	1.36e-07	0.00e+00

Table 4. Table 4: L 1 superscript 𝐿 1 L^{1} -errors and convergence rates in total energy for different pressure perturbations ( η = 0.1 , 10 − 5 𝜂 0.1 superscript 10 5 \eta=0.1,10^{-5} ) on the wave in a gravitational field solution of the 2-D Euler equations after time t = 0.1 𝑡 0.1 t=0.1 . The third order standard and well-balanced method are used. The setup is described in Section   8.9 .

grid cells	$\eta =0.1$				$\eta ={10}^{-5}$
	Std-O3		WB-O3		Std-O3		WB-O3
	$E$ error	rate	$E$ error	rate	$E$ error	rate	$E$ error	rate
$64\times 64$	3.96e-04	2.8	3.84e-04	2.8	4.38e-05	2.7	9.38e-08	2.7
$128\times 128$	5.81e-05	3.0	5.65e-05	3.0	6.68e-06	3.0	1.49e-08	2.9
$256\times 256$	7.50e-06	3.1	7.30e-06	3.1	7.72e-07	3.1	1.93e-09	3.2
$512\times 512$	8.50e-07		8.30e-07		5.46e-08		2.17e-10

Table 5. Table 5: L 1 superscript 𝐿 1 L^{1} -errors and convergence rates in total energy for a 2-D MHD stationary vortex after time t = 0.2 𝑡 0.2 t=0.2 . Different methods are used. In the well-balanced method, a target solution is chosen, which deviates from the actual solution over time. The setup is described in  Section   8.12 .

grid cells	WB-O1		WB-O2		WB-O3		Std-03
	$E$ error	rate	$E$ error	rate	$E$ error	rate	$E$ error	rate
$128\times 128$	2.42e-02	1.0	1.87e-03	1.9	8.46e-05	3.0	8.46e-06	3.3
$256\times 256$	1.23e-02	1.0	5.11e-04	1.9	1.08e-05	3.0	8.38e-07	3.1
$512\times 512$	6.17e-03		1.34e-04		1.36e-06		9.67e-08

Equations

$${\partial }_t\mathbf{q}(\mathbf{x},t)+\nabla \cdot \mathcal{F}(\mathbf{q}(\mathbf{x},t))=\mathbf{s}(\mathbf{q}(\mathbf{x},t),\mathbf{x},t)$$

$${\mathbf{Q}}_i(t):=\frac{1}{V_i}{\int }_{{\Omega }_i}\mathbf{q}(\mathbf{x},t)d\mathbf{x},$$

$$\frac{d}{dt}{\mathbf{Q}}_i(t)+\frac{1}{V_i}{\int }_{\partial {\Omega }_i}\mathcal{F}(\mathbf{q}(\mathbf{x},t))\cdot \mathbf{n}(\mathbf{x})d\sigma =\frac{1}{V_i}{\int }_{{\Omega }_i}\mathbf{s}(\mathbf{q}(\mathbf{x},t),\mathbf{x})d\mathbf{x},$$

$$\frac{d}{dt}{\mathbf{Q}}_i(t)=-\frac{1}{V_i}\sum _{k\in N(i)}{\int }_{\partial {\Omega }_{ik}}\mathcal{F}(\mathbf{q}(\mathbf{x},t))\cdot \mathbf{n}(\mathbf{x})d\sigma +\frac{1}{V_i}{\int }_{{\Omega }_i}\mathbf{s}(\mathbf{q}(\mathbf{x},t),\mathbf{x},t)d\mathbf{x},$$

$$\frac{d}{dt}{\mathbf{Q}}_i(t)=-\frac{1}{V_i}\sum _{k\in N(i)}{\int }_{\partial {\Omega }_{ik}}\mathbf{F}\left({\mathbf{Q}}_i^{\text{rec}}(\mathbf{x},t),{\mathbf{Q}}_k^{\text{rec}}(\mathbf{x},t),\mathbf{n}(\mathbf{x})\right)d\sigma +\frac{1}{V_i}{\int }_{{\Omega }_i}\mathbf{s}\left({\mathbf{Q}}_i^{\text{rec}}(\mathbf{x},t),\mathbf{x},t\right)d\mathbf{x},$$

$$\frac{d}{dt}{\mathbf{Q}}_i(t)=-\frac{1}{V_i}\sum _{k\in N(i)}\left(\sum _{j=1}^M{\omega }_j\mathbf{F}\left({\mathbf{Q}}_i^{\text{rec}}({\mathbf{x}}_{ikj},t),{\mathbf{Q}}_k^{\text{rec}}({\mathbf{x}}_{ikj},t),\mathbf{n}({\mathbf{x}}_{ikj})\right)\right)+\frac{1}{V_i}{I}_{\mathbf{x}\in {\Omega }_i}\left[\mathbf{s}\left({\mathbf{Q}}_i^{\text{rec}},\mathbf{x},t\right)\right].$$

$$I_{x\in \Omega }[\psi ]\approx {\int }_{\Omega }\psi (\mathbf{x})d\mathbf{x}$$

$$\frac{d}{dt}{\tilde{\mathbf{Q}}}_i(t)+\frac{1}{V_i}{\int }_{\partial {\Omega }_i}\mathcal{F}(\tilde{\mathbf{q}}(\mathbf{x},t))\cdot \mathbf{n}(\mathbf{x})d\sigma =\frac{1}{V_i}{\int }_{{\Omega }_i}\mathbf{s}(\tilde{\mathbf{q}}(\mathbf{x},t),\mathbf{x})d\mathbf{x},$$

$$\frac{d}{dt}{\mathbf{Q}}_i(t)-\frac{d}{dt}{\tilde{\mathbf{Q}}}_i(t)+\frac{1}{V_i}{\int }_{\partial {\Omega }_i}(\mathcal{F}(\mathbf{q}(\mathbf{x},t))-\mathcal{F}(\tilde{\mathbf{q}}(\mathbf{x},t)))\cdot \mathbf{n}(\mathbf{x})d\sigma =\frac{1}{V_i}{\int }_{{\Omega }_i}\mathbf{s}(\mathbf{q}(\mathbf{x},t),\mathbf{x},t)-\mathbf{s}(\tilde{\mathbf{q}}(\mathbf{x},t),\mathbf{x},t)d\mathbf{x}.$$

$$\Delta \mathbf{q}$$

$$\frac{d}{dt}(\Delta {\mathbf{Q}}_i(t))=$$

$$+\frac{1}{V_i}{\int }_{{\Omega }_i}\mathbf{s}((\Delta \mathbf{q}+\tilde{\mathbf{q}})(\mathbf{x},t),\mathbf{x},t)-\mathbf{s}(\tilde{\mathbf{q}}(\mathbf{x},t),\mathbf{x},t)d\mathbf{x},$$

$$\Delta \hat{\mathcal{F}}\left(\Delta {\mathbf{Q}}^L,\Delta {\mathbf{Q}}^R,\tilde{\mathbf{q}},\mathbf{n}\right):=\mathbf{F}(\Delta {\mathbf{Q}}^L+\tilde{\mathbf{q}},\Delta {\mathbf{Q}}^R+\tilde{\mathbf{q}},\mathbf{n})-\mathbf{n}\cdot \mathcal{F}(\tilde{\mathbf{q}})\approx \mathbf{n}\cdot (\mathcal{F}(\Delta \mathbf{q}+\tilde{\mathbf{q}})-\mathcal{F}(\tilde{\mathbf{q}})),$$

$$\frac{d}{dt}(\Delta {\mathbf{Q}}_i(t))=$$

$$+\frac{1}{V_i}{\int }_{{\Omega }_i}\mathbf{s}((\Delta {\mathbf{Q}}_i^{\text{rec}}+\tilde{\mathbf{q}})(\mathbf{x},t),\mathbf{x},t)-\mathbf{s}(\tilde{\mathbf{q}}(\mathbf{x},t),\mathbf{x},t)d\mathbf{x},$$

$$\frac{d}{dt}(\Delta {\mathbf{Q}}_i(t))=$$

$$+\frac{1}{V_i}{I}_{\mathbf{x}\in {\Omega }_i}\left[\mathbf{s}((\Delta {\mathbf{Q}}_i^{\text{rec}}+\tilde{\mathbf{q}})(\mathbf{x},t),\mathbf{x},t)\right]-\frac{1}{V_i}{I}_{\mathbf{x}\in {\Omega }_i}\left[\mathbf{s}(\tilde{\mathbf{q}}(\mathbf{x},t),\mathbf{x},t)\right].$$

$$\mathbf{s}((\tilde{\mathbf{q}}+\Delta \mathbf{q})\left(t,\mathbf{x}\right),\mathbf{x},t)-\mathbf{s}(\tilde{\mathbf{q}}\left(t,\mathbf{x}\right),\mathbf{x},t)=\mathbf{s}(\Delta \mathbf{q}\left(t,\mathbf{x}\right),\mathbf{x},t).$$

$$\Delta \hat{\mathcal{F}}\left(\Delta {\mathbf{Q}}^L,\Delta {\mathbf{Q}}^R,{\tilde{\mathbf{q}}}^L,{\tilde{\mathbf{q}}}^R,\mathbf{n}\right):=\mathbf{F}(\Delta {\mathbf{Q}}^L+{\tilde{\mathbf{q}}}^L,\Delta {\mathbf{Q}}^R+{\tilde{\mathbf{q}}}^R,\mathbf{n})-\mathbf{F}({\tilde{\mathbf{q}}}^L,{\tilde{\mathbf{q}}}^R,\mathbf{n}),$$

$$\Delta {\mathbf{Q}}_i=0\forall i\in \{1,\dots ,N\}$$

$$\Delta \hat{\mathcal{F}}\left(\Delta {\mathbf{Q}}^L,\Delta {\mathbf{Q}}^R,\tilde{\mathbf{q}},\mathbf{n}\right)=\Delta \hat{\mathcal{F}}\left(0,0,\tilde{\mathbf{q}},\mathbf{n}\right)=\mathbf{F}(\tilde{\mathbf{q}},\tilde{\mathbf{q}},\mathbf{n})-\mathbf{n}\cdot \mathcal{F}(\tilde{\mathbf{q}})=\mathbf{n}\cdot \mathcal{F}(\tilde{\mathbf{q}})-\mathbf{n}\cdot \mathcal{F}(\tilde{\mathbf{q}})=0.$$

$$I_{\mathbf{x}\in {\Omega }_i}\left[\mathbf{s}((\Delta {\mathbf{Q}}_i^{\text{rec}}+\tilde{\mathbf{q}})(\mathbf{x},t),\mathbf{x},t)\right]-I_{\mathbf{x}\in {\Omega }_i}\left[\mathbf{s}(\tilde{\mathbf{q}}(\mathbf{x},t),\mathbf{x},t)\right]=I_{\mathbf{x}\in {\Omega }_i}\left[\mathbf{s}(\tilde{\mathbf{q}}(\mathbf{x},t),\mathbf{x},t)\right]-I_{\mathbf{x}\in {\Omega }_i}\left[\mathbf{s}(\tilde{\mathbf{q}}(\mathbf{x},t),\mathbf{x},t)\right]=0.$$

$${\tilde{\mathbf{Q}}}_{ij}(t^n)=\frac{1}{{\tilde{V}}_{ij}}\sum _{\hat{}{\mathbf{x}}_{kl}\in {\tilde{\Omega }}_{ij}}{\hat{V}}_{kl}{\hat{\mathbf{Q}}}_{kl}^n,$$

$${\partial }_t\mathbf{q}+{\partial }_x\mathbf{f}+{\partial }_y\mathbf{g}=\mathbf{s},$$

$$\mathbf{q}$$

$$p=(\gamma -1)\rho \epsilon$$

$${\partial }_t\mathbf{q}+{\partial }_x\mathbf{f}+{\partial }_y\mathbf{g}=0.$$

$$\mathbf{q}$$

$$\varphi (x)$$

$$\rho (x)$$

$$\tilde{\rho }(x)$$

Full text

††journal:

1 Introduction

Several problems in engineering and science are modeled by conservation principles and lead to non-linear partial differential equation which are hyperbolic. These equations can rarely be solved exactly and we must resort to some form of numerical approximation. One successful numerical approach for solving hyperbolic PDE is the finite volume method based on Godunov’s idea [1]. Finite volume methods introduce discretization errors such that they are in general not exact on non-trivial solutions. 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. Although numerical methods for hyperbolic balance laws might admit some discrete stationary states, they are in general grid dependent and different from the stationary states of the PDEs. Standard numerical methods are not able to accurately maintain the latter solutions for long times, since they introduce discretization errors which are especially large on coarse grids. This gives rise to the need to develop 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 (e.g. [2, 3, 4, 5] and references therein) and related systems like the Ripa model ([6, 7] and references therein). This includes high order methods for static [8] and non-static stationary states [9]. The relevance for methods for non-static stationary states has been pointed out in [10]. For tsunami modeling applications, high order methods for shallow water equations on non-flat manifolds have been developed e.g. in [11] 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 together with an equation of state (EoS). This makes the construction of well-balanced methods much more difficult and typically restricts the resulting method to some special cases. Many methods have been developed for some classes of hydrostatic states assuming an ideal gas EoS  (e.g. [12, 13, 14, 15, 16, 17, 18, 19, 20] and references therein) and there are also high order methods, see e.g. [21, 22, 23].

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 [24], 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 [25] and is also used by Powell et al. [26]. To do so, the background magnetic field is assumed to be static and free of rotation as well 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 [27, 28, 29]. Path-conservative methods are introduced in [30, 31, 32]. Another widespread idea is the hydrostatic reconstruction, i.e. the reconstruction of variables which are constant if the system is in the considered stationary state. An early example of this method is [3] for the shallow water system; for Euler equations, this approach has e.g. been used in [33, 34, 15, 19, 16, 22, 17]. The methods for Euler equations are restricted to a certain EoS and certain classes of hydrostatic solutions. For astrophysical applications, for example, this restriction is a severe limitation since the 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 [35, 36, 37, 38]. The well-balanced methods introduced in these publications can be applied for any EoS. They are exact on certain hydrostatic solutions, and in all other cases they are exact for a second order approximation of the considered hydrostatic solution. In [34], 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 and it has been extended to high order in [22]. Notably, there is also an extension to stationary states with non-zero velocity in the same article. Similar techniques can be found in the context of numerical atmospheric modeling (e.g. [39, 40, 16]). Those well-balanced schemes strongly rely on the structure of the discretized equations or the static solutions to be well-balanced.

The method we present in this paper is designed in the manner of the method in [34]. It uses the idea of hydrostatic reconstruction and a modification of the source term discretization to obtain the well-balanced property. The main point where 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 known solution which is either given by an analytical expression or as discrete data. Additionally, unlike all the other well-balanced methods mentioned above which deal with time independent hydrostatic solutions, our well-balanced method can be used to follow known time-dependent solutions exactly, as we will show in this paper. We will refer to the solution that is to be exactly captured by the scheme as the target solution. Depending on the application, the target solution can be a time-independent hydrostatic solution, or some general time-dependent solution of interest. Our method is also general in that it possible to combine it 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 high 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 known a priori, e.g., stellar astrophysics, but the EoS is often given in the form of a table. Consequently, hydrostatic solutions can only be found numerically and are available 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 [34, 22] 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 solution. 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 the velocity is not zero in such a solution, conventional well-balanced methods cannot preserve this solution. A special method designed for this application is presented in [41]. 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, but the numerical results will be presented only for 1-D and 2-D test cases. 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. In Section 5 the treatment of discrete target solutions is discussed. 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. 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 7. Finally, in Section 8, we show a variety of numerical tests. The range of applications goes from Euler equations to ideal magnetohydrodynamics (MHD) equations. They include classical well-balanced tests on the balance laws and also tests 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 high order accuracy for solutions close to and far away from the well-balanced solution numerically. A simple example for using a target solution which is obtained via numerical simulation is given. The robustness of our approach is validated in a shock tube on a hydrostatic solution for Euler equations with gravity. To show the efficiency of the method, we present CPU time comparisons of simulations with and without the well-balanced modification in Section 9.

2 A standard finite volume method

In this section we present the standard high order finite volume framework for three-dimensional hyperbolic balance laws [42, 43]. Consider the 3-d system of hyperbolic balance laws

[TABLE]

with $\mathcal{F}=({\bf f}_{1},{\bf f}_{2},{\bf f}_{3})$, where ${\bf f}_{l}$ is the flux in $l$-direction. The spatial domain is partitioned by a mesh consisting of $N$ non-overlapping control volumes. For the $i$-th control volume $\Omega_{i}$ ($i\in\{1,\dots,N\}$), we define the cell-average

[TABLE]

where $V_{i}=|\Omega_{i}|$ is the control cell volume. Integrating Eq. 1 over $\Omega_{i}$ and applying the divergence theorem yields an ordinary differential equation for ${\bf Q}_{i}$

[TABLE]

which can be written as

[TABLE]

where $N(i)$ is the set of indices of all control volumes sharing an interface with $\Omega_{i}$ and $\partial\Omega_{ik}$ denotes the interface between $\Omega_{i}$ and $\Omega_{k}$. For the discretization of the interface fluxes we use a numerical flux function ${\bf F}(\cdot,\cdot,{\bf n})$ consistent with ${\bf n}\cdot\mathcal{F}$. The consistency conditions are Lipschitz continuity in the first two arguments and the relation ${\bf F}({\bf q},{\bf q},{\bf n})={\bf n}\cdot\mathcal{F}({\bf q})$ for all unit vectors ${\bf n}$. We apply this discretization to Eq. 4 and obtain

[TABLE]

where the reconstructed functions ${\bf Q}_{i}^{\text{rec}},{\bf Q}_{k}^{\text{rec}}$ are obtained using a consistent conservative reconstruction routine on the cell average values ${\bf Q}$. Examples for popular consistent conservative reconstruction routines can be found in [44, 45, 46, 43]. 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

[TABLE]

$M$ is the number of quadrature points at the interfaces, ${\bf x}_{ikj}$ are the $M$ quadrature points at the $ik$ interface and $\omega_{j}$ are the corresponding weights. The symbol $I_{x\in\Omega}[\cdot]$ denotes a consistent discretization of the integral over the argument in the domain $\Omega$,

[TABLE]

The quadrature rules and source term discretizations we use in our tests are explained in the Appendix.

The semi-discrete scheme Eq. 6 is $k$-th order accurate if the applied reconstruction routine, interface flux quadrature and source term discretization are all at least $k$-th order accurate. It can then be evolved in time using a $k$-th order accurate ODE solver to obtain a $k$-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.

Let $\tilde{{\bf q}}$ be a given continuous and sufficiently smooth solution of Eq. 1. Plugging this target solution $\tilde{\bf q}$ into Eq. 3 we get

[TABLE]

where $\tilde{{\bf Q}}_{i}$ is the average of the target solution $\tilde{\bf q}$ in the $i$-th control volume. In the next step we subtract Eq. 7 from Eq. 3 to obtain

[TABLE]

Now, let us rewrite Eq. 8 in terms of the deviation from the target solution

[TABLE]

This yields

[TABLE]

where $N(i)$ is the set of indices of all control volumes sharing an interface with $\Omega_{i}$. At this point, we start to discretize. For that we define a numerical flux difference approximation

[TABLE]

where ${\bf F}(\cdot,\cdot,{\bf n})$ is a numerical flux function consistent with ${\bf n}\cdot\mathcal{F}$. We apply this discretization to Eq. 10 and obtain

[TABLE]

where the reconstructed functions $\Delta{\bf Q}_{i}^{\text{rec}},\Delta{\bf Q}_{k}^{\text{rec}}$ are obtained using a consistent conservative reconstruction routine on the cell average values $\Delta{\bf Q}$. Note that using a consistent reconstruction on the deviations $\Delta{\bf Q}$ is equivalent to the application of a hydrostatic reconstruction as e.g. in [3, 47]. A typical hydrostatic reconstruction consists of a transformation to hydrostatic variables (i.e. a set of variables that is constant in the hydrostatic case), a consistent reconstruction, and a transformation back to conservative variables. In our description of the method the deviations $\Delta{\bf Q}$ correspond to the hydrostatic states. Hence, no transformation is required. 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

[TABLE]

where the notations used are as in Eq. 6.

If the source term ${\bf s}$ is linear in the first argument, we can use the following relation

[TABLE]

Due to the linearity of the corresponding source term discretizations, this relation then also holds for the discretized source terms, which further simplifies the scheme. For example, 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.

As in the standard method, this semi-discrete scheme Eq. 13 is $k$-th order accurate if the applied reconstruction routine, interface flux quadrature and source term discretization are all at least $k$-th order accurate. It can then be evolved in time using an at least $k$-th order accurate ODE solver to obtain a $k$-th order accurate fully discrete scheme.

Remark 3.1**.**

In the description of the method, we assume the target solution $\tilde{\bf q}$ to be smooth. In the case of discontinuous $\tilde{\bf q}$, the two values $\tilde{\bf q}^{L}$ and $\tilde{\bf q}^{R}$ which are different if a discontinuity is present at the interface, have to given instead of only one value. Equation 11 has then to be modified to

[TABLE]

i.e. the numerical flux function is also applied to the target solution. If $\tilde{\bf q}$ is continuous, Eq. 15 reduces to Eq. 11 due to the consistency of the numerical flux function ${\bf F}$. However, in the case of a discontinuous target solution $\tilde{\bf q}$ no high order convergence can be expected and the computational cost of the well-balanced modification increases.

4 Proof of the well-balanced property

In this section we show the well-balanced property of our method.

Theorem 4.1**.**

*The modified finite volume method introduced in Section 3 satisfies the following property: If *

[TABLE]

at initial time, then this holds for all $t>0$. Consequently, if the initial condition ${\bf Q}_{i}(t=0)$, $i=1,\dots,N$, equals the cell averages of the target solution $\tilde{\bf Q}_{i}(t=0)$, $i=1,\dots,N$, the computed solution equals the target solution for all time.

Proof: Let $\Delta{\bf Q}_{i}=0$ for all $i\in\{1,\dots,N\}$. The consistency of the applied reconstruction leads to $\Delta{\bf Q}_{i}^{\text{rec}}\equiv 0$ at all flux quadrature points. The flux consistency then yields

[TABLE]

Now, consider the contribution from the source term: With $\Delta{\bf Q}_{i}=0$ the source term discretization in Eq. 13 reduces to

[TABLE]

We have shown that the right hand side in Eq. 13 vanishes and thus the initial data $\Delta{\bf Q}_{i}=0$ are conserved for all time. The second part of the theorem follows easily. ∎

Remark 4.2**.**

If a stationary solution is chosen as target solution (which is the case for classical well-balancing applications), the time derivative of the target solution vanishes by definition. This leads to $\frac{d}{dt}{\bf Q}_{i}=\frac{d}{dt}(\Delta{\bf Q}_{i})$. The described method can then also be used to directly evolve the ${\bf Q}_{i}$ in time instead of $\Delta{\bf Q}_{i}$.

5 The treatment of discrete target solutions

The target solution, introduced in Section 3 for the well-balanced modification, has not to be known analytically, it can also be given only 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 target 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. Note that the target solution obtained from the discrete data are in general only approximate solutions of the hyperbolic system, but we assume that it is a sufficiently accurate approximation to the true solution of the PDE. The well-balanced modification introduced in this article is then exact on this approximate solution.

Depending on the form in which the discrete target data are given, e.g. point values or cell averages, and the required order of accuracy, different methods for this reconstruction have to be used. Also, in the case of a discontinuous target solution (see 3.1), two interface values of $\tilde{\bf q}$ have to be determined for each interface quadrature point instead of only one. In the following we give two examples how cell average and interface values of a target solution can be obtained from discrete data.

5.1 Example 1: pointwise 1-D data on a fine grid

An important application of our method is in computing hydrostatic solutions of Euler equations and perturbations around such solutions. 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 ${\bf q}^{hs}_{i}$, $i=1,2,\dots,L$ (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 ${\bf q}^{hs}_{1-D}(x)$ (e.g. [48]). 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:

• (a)

Suppose we have an essentially 1-D hydrostatic solution where the gravitational force is at an angle $\alpha\in[0,2\pi)$ to the $x$-axis. Then we extend the one-dimensional hydrostatic solution ${\bf q}^{hs}_{1-D}$ to a two-dimensional solution via ${\bf q}^{hs}({\bf x}):={\bf q}^{hs}_{1-D}(x\cos(\alpha)+y\sin(\alpha))$.

• (b)

Suppose we have a radial hydrostatic solution with a gravity vector pointing towards the center ${\bf x}_{\text{center}}$. In that case we extend the hydrostatic solution back to 2-D by setting ${\bf q}^{hs}({\bf x}):={\bf q}^{hs}_{1-D}(\|{\bf x}-{\bf x}_{\text{center}}\|)$.

The values for the target solution at interface quadrature points can then be evaluated pointwise as $\tilde{\bf q}({\bf x},t)={\bf q}^{hs}({\bf x})$. The cell average values of the target solution are computed using a third order accurate 2-D Gauß–Legendre quadrature rule. Method (a) is used in the test in Section 8.6.

5.2 Example 2: target solution from a highly resolved finite volume simulation

Consider a target 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 $\hat{{\bf Q}}^{n}_{ij}$, instead of point values. Assume this target solution is supposed to be used in a $k$-th 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 target solution has been computed, such that all interfaces on the coarse grid coincide with interfaces of the fine grid. For each time step $t^{n}$ of the stored target data we map the fine grid on the coarse grid with

[TABLE]

where all quantities with $\tilde{\cdot}$ correspond to the coarse grid and all quantities with $\hat{\cdot}$ correspond to the fine grid. The values of $\tilde{\bf Q}_{ij}$ at intermediate times are obtained via a $k$-th order accurate interpolation in time.

The value of the target solution at all quadrature points required in the scheme are obtained using a $k$-th order accurate interpolation on the cell-centered point values $\tilde{\bf q}^{\text{rec}}({\bf x}_{ij})$. Those cell-centered values are obtained using a $k$-th order accurate conservative reconstruction on the cell-averages of the cell-average values $\tilde{{\bf Q}}_{ij}$. This method is applied in a numerical test in Section 8.13 for $k=3$.

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 some boundary conditions which are compatible with the well-balancing property and support the potentially high order accuracy of the scheme. Some of the proposed numerical boundary conditions require knowledge of the target solution outside the domain. If this is not available, one can simply extrapolate the target solution to the ghost cells using a sufficiently high order accurate extrapolation. The well-balanced scheme will in that case exactly balance the target solution including the extrapolated values. Consequently, this will not affect the well-balanced property nor order of accuracy.

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 $\Delta{\bf Q}$. In the case that $\Delta{\bf Q}_{i}=0$ 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 and periodic boundary conditions Periodic boundary conditions and wall boundary conditions as e.g. described in [49] can be applied to the deviations $\Delta{\bf Q}$. In the case of $\Delta{\bf Q}=0$ for each cell in the domain, this implies $\Delta{\bf Q}=0$ for each ghost cell. Hence, these boundary conditions are consistent with the well-balanced property.

7 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 for applications: The fact that the method is also easy to be implemented. This holds especially if there is an existing finite volume code evolving Eq. 6. It can be easily modified to evolve Eq. 13 and hence obtain a well-balancing capability. The changes that should be introduced in a typical finite volume scheme are the following:

Implement a function able to return the target solution $\tilde{\bf q}$ in any point $({\bf x},t)$, or a structure containing the cell averages and all the values needed at each quadrature point if $\tilde{\bf q}$ is time-independent. 2. 2.

When the main routine receives the initial cell-averages ${\bf Q}_{i}$, it should transform it to the deviations $\Delta{\bf Q}_{i}={\bf Q}_{i}-\tilde{{\bf Q}_{i}}$ and work directly with the deviation. 3. 3.

In the routine evaluating the numerical flux, it has to be computed $\Delta\hat{\mathcal{F}}\left(\Delta{\bf Q}^{L},\Delta{\bf Q}^{R},\tilde{\bf q},{\bf n}\right)$ instead of ${\bf n}\cdot\mathcal{F}(\Delta{\bf q}+\tilde{{\bf q}})$. Basically, this just means subtracting the exact flux after the evaluation of the standard numerical flux. 4. 4.

This step is only necessary if the source term is not linear in ${\bf q}$: Evaluate the source term at the states $\tilde{\bf Q}_{i}+\Delta{\bf Q}_{i}$ and $\tilde{{\bf Q}_{i}}$. The difference of these source terms is computed and used to evolve the approximate solution as described in Eq. 13.

Let us remind of 4.2 and point out that an alternative implementation could also evolve ${\bf Q}$ instead of $\Delta{\bf Q}$ in time if $\tilde{{\bf q}}$ is time-independent. In the computer implementation, this can be the easier and computationally efficient approach than evolving the deviations.

8 Numerical tests of the scheme

8.1 Hyperbolic systems used in the tests

Since the well-balancing procedure introduced in this article is applicable to any hyperbolic balance law, we present numerical experiments for two different hyperbolic systems. They are introduced in the following.

8.1.1 Compressible Euler equations with gravity source term

The 2-D compressible Euler equations which model the balance laws of mass, momentum, and energy under the influence of gravity are given by

[TABLE]

where the conserved variables, fluxes and source terms are

[TABLE]

with $\rho,p>0$. Moreover, $E=\rho\varepsilon+\tfrac{1}{2}\rho|{\bf v}|^{2}+\rho\phi$ is the total energy per unit volume with the velocity ${\bf v}=(u,v)^{T}$ and specific internal energy $\varepsilon$. The scalar function $\phi$ 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

[TABLE]

with $\gamma=1.4$, 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 ${\bf g}=0$ and removing the $\rho v$ equation. It can be reduced to homogeneous Euler equations by setting ${\bf s}=0$.

8.1.2 Homogeneous compressible ideal magnetohydrodynamics

The 2-D compressible ideal magnetohydrodynamics (MHD) equations which model the conservation of mass, momentum, magnetic field, and energy are given by

[TABLE]

The conserved variables and fluxes are

[TABLE]

where $B_{x}$, $B_{y}$ are the $x$- and $y$-component of the magnetic field. The total energy is $E=\rho\varepsilon+\tfrac{1}{2}\rho|{\bf v}|^{2}+\frac{1}{2}(B_{x}^{2}+B_{y}^{2})$. All other quantities 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 $\rho w$ and $B_{z}$ components. This is in principle reasonable due to the genuine three-dimensional interactions between velocity and magnetic field. In our tests we set $\rho w$ and $B_{z}$ to zero and there is no difference if we omit the corresponding equations.

8.2 Code and numerical methods

We test the methods proposed in this paper using a finite volume code implemented in Python. The code is built in a modular way, such that different schemes can be easily implemented. For brevity, we only give a short description of the methods in this section and the interested reader can find details in the A. Note that our well-balanced method is not restricted to the methods we choose to use in the tests and one can for example also use other reconstruction methods or quadrature formulae. One can also apply numerical flux functions designed for special problems (e.g. a low Mach number compliant method for Euler equations like in [50, 51]).

Grids: The domain in two-dimensional problems is discretized using a structured grid and in some tests, we use curvilinear grids. The implementation of these grids restricts the overall method to only second order accuracy. Note that this is not a general statement about the proposed method, but this restriction is due to the special implementation of the grid in our code. Higher order accuracy can only be achieved with a Cartesian grid in our code.

Numerical flux function: As numerical flux function we use the local Lax–Friedrichs flux (e.g. [42]), 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 [52] 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 and the source term is evaluated at the cell-center. For the gravity source term used in our tests, we need the gradient $\nabla\phi$ of the given gravitational potential which is evaluated at the cell center using analytical differentiation.

Second order method: To formally obtain a second order method, we use a conservative linear reconstruction (e.g. [43]) with a minmod limiter (e.g. [43]) 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 CWENO3 ([46] for 1-D, [53] for 2-D) reconstruction to obtain the interface values. In the two-dimensional case, the numerical fluxes are evaluated at the Gauß–Legendre quadrature points. To compute the source term we multiply the CWENO3 polynomials in momentum with the interpolation polynomial of $\partial_{x}\phi$ or $\partial_{y}\phi$ respectively. This results in a polynomial function of the spatial variables whose cell average can be computed by analytical integration.

Seventh order method: To formally obtain a one-dimensional seventh order method, we use a conservative CWENO7 reconstruction [54] to obtain the interface values. To compute the source term we multiply the CWENO7 polynomial in momentum with the interpolation polynomial of $\partial_{x}\phi$. The resulting source term polynomial is cell-averaged using analytical integration.

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: The first order accurate scheme is evolved in time using the explicit forward Euler method. For the second and third order accurate semi-discrete scheme we use the explicit third order, four stage Runge–Kutta method from [55] and the explicit tenth order, 17 stage Runge–Kutta method from [56] is used in the seventh order accurate method.

8.3 1-D isothermal hydrostatic solution of the Euler system with gravity

We consider an isothermal hydrostatic solution of the 1-D compressible Euler equations with the ideal gas law and the gravitational source term given by

[TABLE]

We set these data on a 1-D grid with 128 grid cells on the domain $[0,1]$. These initial data are evolved up to the final time $t=2$ 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. 25 as time-independent target solution. The $L^{1}$-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.

Remark 8.1**.**

Most other well-balanced methods balance fluxes against the source term which leads to machine errors. In our method we balance fluxes against fluxes and source term against source term. Thus, the differences can cancel out exactly and the error can be exactly zero.

8.4 Perturbed 1-D isothermal hydrostatic solution of the Euler system with gravity

We add a perturbation to the pressure such that our initial conditions are

[TABLE]

in the domain $[0,1]$. We choose $\eta=0.1$ to test the convergence of our method. We evolve this initial setup up to time $t=0.2$ 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 $4096$ 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. 1 and 2, density deviations at time $t=0.2$ for the test with $\eta=10^{-5}$ 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 8.3.

8.5 Riemann problem on a 1-D isothermal hydrostatic solution of the Euler system with gravity

To test the robustness of our well-balanced methods in combination with CWENO reconstruction we use the initial data

[TABLE]

with $\phi(x):=-10x$. Eq. 27 describes a piecewise isothermal hydrostatic solution with a jump, which includes all three waves of the Euler equations. We set these initial data on the domain $[0,0.25]$ and evolve them to the final time $t=0.02$ using our third and seventh order well-balanced method with CWENO reconstruction and Roe’s approximate Riemann solver on $128$ grid cells. As target solution for the well-balanced methods we choose

[TABLE]

The results at final time are presented in Fig. 3. As a reference we use a result obtained with a first order standard method on 8192 grid cells (top right panel). Neither the third order method (bottom left panel) nor the seventh order method (bottom right panel) show significant oscillations. The wave structure is captured correctly by both methods.

8.6 2-D numerically approximated hydrostatic solution of the Euler system with gravity

In stellar astrophysics 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 target solution is given in the form of discrete data in a table.

The thermodynamical quantities shall be related by the EoS for an ideal gas with radiation pressure, which is given by[57]

[TABLE]

where the temperature $T$ is defined implicitly via

[TABLE]

We assume the following data is given. Let the gravitational potential be $\phi({\bf x})=\phi(x,y)=x+y$ and the hydrostatic temperature profile is $\bar{T}({\bf x})=1-0.1\phi({\bf x})$.

Using Chebfun [58] 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. 4.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 $64\times 64$ grid to evolve the hydrostatic initial condition to the final time $t=2$. For the conversion between pressure and internal energy we use Newton’s method to solve for the temperature. The $L^{1}$-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.

8.7 Double Gresho vortex

In this test we use a vortex for homogeneous 2-D Euler equations first introduced in [59]. The pressure and angular velocity of this vortex in dependence of the distance to the center are given by

[TABLE]

The radial velocity is zero and the density is $\rho\equiv 1$. In our test we set up the domain $[0,1]\times[0,2]$ with two Gresho vortices centered at $(0.5,0.5)$ and $(0.5,1.5)$ respectively. The vortices are advected with the velocity ${\bf v}_{0}=(u_{0},v_{0})^{T}=(0.2,0.4)^{T}$ and the boundaries are periodic. At time $t=5$ the exact solution of this initial data equals the initial setup. We apply our well-balanced method on a $64\times 128$ grid to evolve the initial condition up to final time $t=5$. We use Roe’s numerical flux functions and a linear reconstruction. Only the vortex initially (and finally) centered at $(0.5,0.5)$ is included in the target solution. The result is illustrated in Fig. 5.

8.8 2-D Euler wave in gravitational field

To demonstrate that we can follow time-dependent solutions exactly with our method we use a problem from [21] and [60] which involves a known exact solution of the 2-D Euler equations with gravity given by

[TABLE]

The gravitational potential is $\phi({\bf x})=x+y$, the EoS is the ideal gas EoS. In accordance to [21] and [60] we choose $u_{0}=v_{0}=1$, $p_{0}=4.5$ on the domain $[0,1]^{2}$. We use the first, second, and third order accurate well-balanced method to evolve the initial data with $t=0$ to a final time $t=0.1$ on a $64\times 64$ Cartesian grid and the second order well-balanced method on a $64\times 64$ polar grid. The $L^{1}$-error in every component of the state vector is exactly zero in each of the tests. We omit showing a table since it does not provide additional insight.

8.9 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 target solutions if the well-balanced method is used. For this we use the initial setup from Eq. 32 and add a pressure perturbation:

[TABLE]

We evolve these initial data to time $t=0.1$ using the third order standard and well-balanced method with $\eta=0.1$. The $L^{1}$ errors and corresponding convergence rates are presented in Table 4. As reference solution for determining the error we use a numerically approximated solution computed using the third order standard method on a $1024^{2}$ 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. 32. 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 $\eta=0.1$ leads to a large discretization error in the perturbation which seems to dominate the total error. Choosing a large perturbation 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 provide a sufficiently accurate reference solution. To yet show the improved accuracy of the well-balanced modification we add a convergence test with a small perturbation of $\eta=10^{-5}$ for which a sufficiently accurate reference solution is produced using the third order well-balanced method on a $1024^{2}$ grid. The errors and convergence rates for the third order accurate standard and well-balanced methods can also be seen in Table 4 and it gets evident that the well-balanced method is significantly more accurate on the small perturbation.

8.10 2-D Keplerian disk

Consider a stationary solution given by [20]

[TABLE]

with the gravitational potential $\phi(r)=-\frac{Gm_{s}}{r}$ and $r=\sqrt{x^{2}+y^{2}}$, $\alpha=\arctan(\frac{y}{x})$, $G=m_{s}=1$. We use the initial conditions

[TABLE]

and $u=\tilde{u},v=\tilde{v},p=\tilde{p}$ on the domain $[-2,2]\times[-2,2]$. We chose the domain such that we omit the singularity in the velocity at $(x,y)=(0,0)$. In Fig. 6 results of numerical tests are illustrated. The second order standard and well-balanced methods are applied on a polar grid with $32\times 256$ cells and a Cartesian grid with $128\times 128$ cells. In the Cartesian grid we take out the center with $r<1$ 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 to the linear reconstruction of the conserved variables in the standard method or the deviations in conserved variables in the well-balanced method. The density at time $t=2.5$ for each simulation is shown in Fig. 6 together with the exact solution. Since this is an purely advective problem and there is no radial component to the velocity, the quantity $\|(\rho-1)r\|_{1}$, 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 $d=\|(\rho(t=2.5)-1)r\|_{1}/\|(\rho(t=0)-1)r\|_{1}$ as a measure of the quality of the numerical solutions. For the exact solution we have $d=1$ for all time. The values of $d$ are shown in the center of the plots in Fig. 6.

In the tests with the standard method we see discretization errors in the Keplerian disk solution Eq. 34. 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 the result is free of discretization errors in the Keplerian disk solution Eq. 34. 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 is more diffusive than the one on the polar grid, yet we see that the well-balanced modification improves the result significantly on both grids.

8.11 Stationary MHD vortex - long time

We consider the following exact solution of the homogeneous 2-D ideal MHD equations:

[TABLE]

This setup describes a stationary vortex which is advected through the domain with the velocity $(u_{0},v_{0})$. The domain is $[-5,5]\times[-5,5]$. One vortex turnover-time is $t_{\text{turnover}}=\frac{2\pi}{\sqrt{e}k_{p}}\approx\frac{3.81}{k_{p}}$. In a first test we set $m_{p}=k_{p}=0.1$, $u_{0}=v_{0}=0$ and run the test up to $t=100t_{\text{turnover}}$ on a $32\times 32$ grid. We use the well-balanced method and the target solution equals the initial data. The numerical error at final time compared to the initial setup is exactly zero in all conservative variables.

8.12 Stationary MHD vortex - order of accuracy

In a second test with the vortex from Section 8.11 we want to see if the well-balanced method converges as expected, even if the target solution which is set deviates from the actual solution over time. For that we set $m_{p}=k_{p}=0.1$, $u_{0}=v_{0}=0$ in the initial condition. As target solution we use the same vortex but with $u_{0}=v_{0}=1$. In Table 5 the $L^{1}$ errors and rates in energy at final time $t=0.2$ are presented for the formally first, second, and third order accurate well-balanced method. We see that even if the target solution moves away from the actual solution over time the method is consistent with the expected order of accuracy.

8.13 Stationary MHD vortex - numerical target solution

In this test we present a simple application of the method described in Section 5.2. Again, we use the stationary MHD vortex test case described in Section 8.11. The parameters are $k_{p}=m_{p}=0.1$ and $u_{0}=v_{0}=0.1$, the final time is $t_{\text{final}}=5$. First, we compute a target solution using our third order non-well-balanced method with $128\times 128$ grid cells with parabolic extrapolation boundary conditions. Every time-step is stored. This is then used in our well-balanced method as described in Section 5.2 using third order interpolation in time and in space. The resulting pressure for well-balanced and standard methods on different Cartesian meshes is shown in Fig. 7. All methods use CWENO3 reconstruction and parabolic extrapolation boundary conditions. On the $128\times 128$ grid the solutions for the well-balanced and non-well-balanced method are exactly the same, since the solution from the standard method is used as target solution in the well-balanced method. For smaller resolutions the standard method is too diffusive to resolve the vortex. The quality of the results obtained with the well-balanced method is the same for all resolutions, since all of them use the same $128\times 128$ simulation as target solution.

9 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.

9.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 8 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. The ratio of average wall clock time for the well-balanced compared to the standard method is visualized depending on the grid resolution.

Note that we use final times that are significantly smaller than the final times used in the corresponding tests above. The reason for this is that after some time the solutions obtain with and without well-balancing differ. This results in different sizes for the time steps which significantly influences the wall-clock time necessary to reach the final time. Since we aim to compare the efficiency of the methods without taking the quality of the solution into account, we use reduced final times.

9.2 The tests

To compare the runtimes we use test setups described in Section 8 using exactly the same methods. In a first test, we use the perturbed one-dimensional isothermal solution described in Section 8.4 with a final time $t=0.02$. 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 the top right panel of Fig. 8. The target solution $\tilde{q}$ 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 $20\%$ for the first and second order accurate methods. The difference in runtime reduces significantly for methods with higher order of accuracy. For the seventh order methods the increase in wall-clock time is only about $5\%$.

As a second test setup we choose the Keplerian disk from Section 8.10. We evolve it to time $t=0.01$ with a first and second order method on a polar grid. As in the previous test, the target solution is time-independent and thus only computed once. The wall-clock time ratios are visualized in the bottom left panel of Fig. 8. We observe an increase in CPU time of less than $20\%$ when using the first order well-balanced method and less than $10\%$ when using the second order well-balanced method.

To also test the increase of CPU time consumption for a simulation in which the target solution is time dependent, we use the Euler wave in a gravitational field with perturbation from Section 8.9. We evolve the solution up to the final time $t=0.01$ with each method. In this test, the target 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 the bottom right panel of Fig. 8. We see an increase in CPU time of less than $30\%$ if the well-balanced method is used. Note that, again, the wall-clock time ratio is smaller for methods with higher order accuracy. For the third order well-balanced method we only observe an increase of about $15\%$ in wall-clock time.

10 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 target solution instead of the actual solution. This makes the scheme exact on the target 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 be arbitrarily high order accurate and the method can be defined on any computational grid geometry. One can view our method as a high order extension of [34] and [22] 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, i.e. stars that are in a rotating stationary state with an angular velocity depending on the distance to the axis of rotation, 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. Also, in a series of numerical tests we have shown that the increase in computational time is moderate.

The proposed well-balanced method can be easily implemented in existing finite volume codes with minimal effort. However, there are applications in which the well-balanced solution is not known beforehand and the current method cannot be applied. In that case, another well-balanced methods existing in the literature for the balance law under consideration has to be applied, in case such a scheme exists. In all other cases, in which the well-balanced solution is known, our simple framework can be applied to obtain the well-balanced 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 TIFR-CAM, Bangalore and Department of Atomic Energy, Government of India, under project no. 12-R&D-TFR-5.01-0520.

Appendix A Details of the applied finite volume schemes

We use structured grids in all the numerical tests and hence in the description of the details, we restrict ourselves to structured grids. Some parts of the scheme, such as the reconstruction methods, are applied to ${\bf Q}$ in the standard method but to $\Delta{\bf Q}$ in the well-balanced method. We will denote the states with ${\bf U}$; depending on the method we have ${\bf U}={\bf Q}$ or ${\bf U}=\Delta{\bf Q}$. Analogously, we will use ${\bf u}$ to denote ${\bf q}$ or $\Delta{\bf q}$.

A.1 Curvilinear coordinates

We define a 2-D curvilinear coordinate system. The coordinates in physical space are ${\bf x}=(x,y)$, the coordinates in computational space are $\bm{\xi}=(\xi,\eta)$. The $(i,j)$-th cell is denoted $\Omega_{i,j}$ in the physical space and by $\hat{\Omega}_{i,j}$ in the computational space. We can rewrite Eq. 20 in the computational coordinates as

[TABLE]

where

[TABLE]

To solve Eq. 20 on the curvilinear physical grid, we can now solve Eq. 37 on a Cartesian grid. We construct the grid from the nodes and approximate the derivatives of the coordinate transformation using central differences on the nodal coordinates. This implementation of curvilinear grids restricts the scheme to second order accuracy. We can achieve higher order accuracy only on Cartesian grids. More details on the finite volume method on a curvilinear mesh can be found in [34].

Polar grid: The polar grid can be defined by the function

[TABLE]

for $\xi>0$, $\eta\in\left[0,2\pi\right)$. Note, that this functions can not be inverted at $\xi=0$, i.e. ${\bf x}=0$. Hence, the origin in physical coordinates has to be omitted, when this grid is used.

A.2 Implementation of boundary conditions

Boundary conditions are implemented on the structured grid by using ghost cells to artificially extend the computational domain. In this section we assume the physical domain is in the cells $\Omega_{ij}$, $(i,j)\in\{0,\dots,N\}\times\{0,\dots,M\}$. The necessary amount of ghost cells depends on the stencil of the method, especially on the quadrature and reconstruction (plus one row of ghost cells for the flux computation at the boundary). The values in the ghost cells are set after each intermediate Runge–Kutta step. In the following, the used boundary conditions are presented for 2-D grids. If there is no description for the 1-D method, the method is reduced to 1-D in the trivial way.

Constant extrapolation: The constant extrapolation boundary conditions are suitable to be used with a first oder accurate method. They are obtained by setting

[TABLE]

in the ghost cells.

Linear extrapolation: In combination with a second order accurate method on a Cartesian grid, the following linear extrapolation to the ghost cells can be used as boundary conditions: We set

[TABLE]

for $(i,j)\in\{0,\dots,N\}\times\{0,\dots,M\}$, $k=1,2$. The diagonal ghost cells ($\Omega_{N+1,M+1}$, etc.) are not needed for the second order scheme used in our tests.

Parabolic extrapolation: Parabolic extrapolation to the ghost cells is suitable to be used in combination with a third order accurate method. On a 1-D equidistant grid we set the values in the ghost cells to

[TABLE]

On a two-dimensional grid we have to use a genuine 2-Dimensional parabolic extrapolation to obtain third order accuracy. The basis for that is a 2-D parabolic reconstruction parabola. We use the $P_{\text{opt}}$ parabola from [46] which is given by

[TABLE]

for $(x,y)\in\Omega_{ij}$. The values in the ghost cells are then computed by integrating the closest possible reconstruction parabola (computed from only values inside the domain) over the ghost cell. This yields

[TABLE]

and correspondingly for $\Omega_{-k,j}$, $\Omega_{i,M+k}$, and $\Omega_{i,-k}$ with $k=1,2,..$ and $(i,j)\in\{1,N-1\}\times\{1,M-1\}$. At the upper right edge we obtain

[TABLE]

for a ghost cell $\Omega_{N+k,M+l}$($k,l\in\{0,1,2,..\}^{2}\setminus\{0\}^{2}$) and correspondingly for the other edges.

A.3 Interface quadrature

For the two-dimensional third order method it is necessary to apply a quadrature rule to compute the interface flux. For that we use the Gauß–Legendre quadrature rule. For that, in Eq. 13, we use $M=3$ and the weights

[TABLE]

The quadrature points are

[TABLE]

at the interface between the $\Omega_{ik}$ and the $\Omega_{i+i,k}$ control volumes and

[TABLE]

at the interface between the $\Omega_{ik}$ and the $\Omega_{i+i,k}$ control volumes. Note the change in notation from Eq. 13 to Eqs. 47 and 48. In Eq. 13 the indices $i$ and $k$ are indices of the $i$-th and $k$-tho control volumes. In Eqs. 47 and 48 we consider the case of a structured grid and the indices denote the position of the $\Omega_{ik}$ control volume in the grid. The interface is then denoted using half values.

A.4 Source term discretization

In some tests we use a gravity source term for Euler equations. The source term component in the momentum equation has to be approximated to sufficiently high order.

Second order source term discretization: For the first and second order method, we use the second order accurate source term discretization

[TABLE]

in the one-dimensional case and

[TABLE]

in the two-dimensional case in the momentum equation. The cell-averaged density is denoted $\bar{\rho}$ and the gravitational acceleration ${\bf g}=\nabla\phi$ is given exactly at the cell-center.

Third and seventh order source term discretization in 1-D: For the one-dimensional methods with CWENO reconstruction we define

[TABLE]

where $\rho^{\text{rec}}$ is the density polynomial obtained from the CWENO reconstruction and $g^{\text{int}}$ is the gravitational acceleration interpolated from the cell centered values to third or seventh order respectively. Since $s^{\rho u}_{h}$ is a polynomial, the cell-average

[TABLE]

can be computed exactly and is used as source term approximation in the momentum equation in the scheme.

Third order source term discretization in 2-D: Similar to the 1-D case we define the vector-valued function

[TABLE]

where $\rho^{\text{rec}}$ is obtained from the two-dimensional CWENO3 reconstruction and ${\bf g}^{\text{int}}$ is given by the interpolation

[TABLE]

for $(x,y)\in\Omega_{ij}$. This polynomial is constructed such that it satisfies ${\bf g}^{\text{int}}({\bf x}_{ij})={\bf g}_{ij}$, ${\bf g}^{\text{int}}({\bf x}_{i\pm 1,j\pm 1})={\bf g}_{i\pm 1,j\pm 1},$ and ${\bf g}^{\text{int}}({\bf x}_{i\pm 1,j\mp 1})={\bf g}_{i\pm 1,j\mp 1}$ for the cell-centered point values in coordinate direction. The diagonal values ${\bf g}_{i+1,j+1}$, ${\bf g}_{i-1,j+1}$, ${\bf g}_{i+1,j-1}$, and ${\bf g}_{i-1,j-1}$ are approximated in the least square sense.

Since $s^{\rho u}_{h}$ and $s^{\rho v}_{h}$ are polynomials, the cell-averages

[TABLE]

can be computed exactly and are used as source term approximations in the momentum equations in the scheme.