We consider 1d systems of balance laws of the form
where takes value in , is the flux function; ; and is a known function from (possibly the identity function ). The system is supposed to be hyperbolic, i.e. the Jacobian of the flux function is assumed to have
different real eigenvalues. PDE systems of this form appear in many fluid models in different contexts: shallow water models, multiphase flow models, gas dynamic, elastic wave equations, etc.
Systems of the form (1) have non trivial stationary solutions that satisfy the ODE system:
The objective of well balanced schemes is to preserve exactly some of these steady state solutions. In the context of shallow water equations Bermúdez and Vázquez-Cendón introduced in  the condition called C-property: a scheme is said to satisfy this condition if it solves correctly the steady-state solutions corresponding to water at rest. Since then, the design of high-order well-balanced numerical methods for systems of balance laws has been a very active front of research: see, for instance, , , , , , , , , ,  , , , , , , , , ,  , , …
We focus here on finite difference methods. Previous well-balanced high-order finite difference types have been based in techniques like the one proposed in  consisting on a formal transformation of (1) into a conservative system through the definition of a ‘combined flux’ formed by the flux and a primitive of the source term (see, for instance, , ). In -, finite difference WENO schemes where applied to problems in which can be written as a function of variables that are constants for the stationary solutions to be preserved.
The technique introduced here is different and it is based on a simple idea: if is the numerical approximation of the solution at the node at time , then if one can find a stationary solution satisfying
Therefore, the source term can be numerically computed with high order accuracy by applying the finite difference reconstruction operator to the point values . We will show that this simple strategy leads to high-order numerical methods that preserve either one particular stationary solution or all of them. Of course, the main difficulty comes from the computation of the stationary solution at every node at every time step.
This strategy has been inspired by the concept of well-balanced reconstruction introduced in  to develop well-balanced high-order finite volume numerical methods: see  for a recent follow-up. In that case, a stationary solution whose average is the numerical approximation at every cell has to be computed at every time step. From this point of view, the technique introduced here for finite difference methods is easier, since the local problems to be solved are standard Cauchy problems for ODE system (2). In all the numerical tests considered here, explicit or implicit expressions of the stationary solutions are available, which will allow us to find easily. For cases where this is not possible, a numerical method for ODE can be used to compute at the stencil of in the spirit of .
While the finite volume methods based on well-balanced reconstructions have the property of being conservative for any conservation law contained in (1) this is only the case for the finite difference methods that preserve one particular stationary solution to be described here but not for those that preserve every stationary solutions.
The case in which
has jump discontinuities will be also considered. At a discontinuity of, a solution is expected to be discontinuous too and the source term cannot be defined within the distributional framework: it becomes a nonconservative product whose meaning has to be specified. There are different mathematical theories that allow one to give a sense to nonconservative products. In the theory developed in , nonconservative products are interpreted as Borel measures whose definition depends on the choice of a family of paths that, in principle, is arbitrary. As the Rankine-Hugoniot conditions, and thus the definition of weak solution, depend on the selected family of paths, its choice has to be consistent with the physics of the problem. Although for general nonconservative systems the adequate selection of paths may be difficult, in the case of systems of balance laws with singular source term there is a natural choice in which the paths are related to the stationary solutions of a regularized system: the interested reader is addressed to  for a detailed discussion. Although in general finite volume methods can be more easily adapted to deal with nonconservative products than finite difference schemes (see ), it will be shown that the numerical methods that preserve every stationary solution can be easily adapted to deal properly with singular source terms.
The organization of the article is as follows: finite difference high order methods based on a reconstruction operator are recalled in Section 2 where the particular example of WENO methods is highlighted. In Section 2 the case in which is continuous and a.e. differentiable is considered. First, the numerical methods that preserve any stationary solution are introduced, together with the proofs of their accuracy and their well-balanced property. Next, methods that preserve only one stationary solution are introduced and the conservation property is discussed. The implementation of well-balanced WENO methods is also discussed. Section 3 is devoted to the more difficult case in which is a.e. differentiable and piecewise continuous with isolated jump discontinuities. The definition of the nonconservative product is briefly discussed and the numerical methods introduced in Section 2 are adapted to this case. The numerical methods are applied to the linear transport equation with linear source term, Burgers’ equation with a nonlinear source term, and the shallow water equations. In particular, it will be shown that the methods introduced here preserve all the stationary solutions for the shallow water, including cases in which the bottom has a step. Finally, some conclusions are drawn and further developments are discussed.
2 Numerical methods
2.1 General case
We will consider here high-order finite difference numerical methods for (1) based on a reconstruction operator, i.e. an operator which, given the point values of a function in the nodes of a mesh
provides a numerical flux function
such that the flux difference approximates the derivative to -th order accuracy
when is smooth enough. Finite-difference ENO or WENO reconstructions are examples of such reconstruction operators: see , , . Uniform meshes of constant step will be considered and the following notation will be also used for the intercells
2.2 WENO reconstructions
In the particular case of the WENO reconstruction of order , two flux reconstructions are computed using the values at the points :
Once these reconstructions have been computed, an upwind criterion can be chosen to define . For instance, for scalar problems
an approximate value of is chosen and the numerical flux is defined then as follows:
For systems, a matrix with real different eigenvalues
that approximates the Jacobian of the flux function has to be chosen and then the numerical flux can be defined by
Here is the diagonal matrix whose coefficients are
is a matrix whose columns are eigenvectors.
An alternative approach is to split the flux
in such a way that the eigenvalues of the Jacobian (resp. ) of (resp. ) are positive (resp. negative). Then, the reconstruction operator is applied to :
A standard choice is the Lax-Friedrichs flux-splitting:
In both cases (the upwind or the splitting implementations) the values used to compute are those at the points
so that, in the general notation, and .
3 Well-balanced high-order finite difference methods: continuous
3.1 Definition of the method: general case
In order to tackle the difficulties gradually, let us suppose first that is continuous and a.e. differentiable. In this case, the stationary solutions of (1) are the solutions of the ODE system:
Once the reconstruction operator has been chosen, the following semi-discrete numerical method is proposed:
where are computed as follows:
Look for the stationary solution such that:
The following result holds:
Observe that the method is well-defined if the first step of the reconstruction procedure can be always performed, i.e. if for every the Cauchy problem
has a unique solution whose interval of definition contains the stencil . When (17) has no solution, or when it has a unique one but the stencil is not contained in its interval of definition, the cell values , to which the reconstruction operator is applied, cannot be the point values of a stationary solution of (1). In this case, there is no stationary solution to preserve and the numerical method (3) is applied to compute . When (17) has more than one solution, an additional criterion is needed to select one of them in the first stage of the reconstruction procedure.
In the notation the index corresponds to the intercell and the index to the center of the cell where the initial condition of (17) is imposed. Therefore, in general
as one can expect due to the non conservative nature of the system of equations. Notice that two reconstructions have to be computed at every stencil: and .
If the Jacobian of the flux function is regular, the ODE system can be written in normal form:
and it is expected to have a unique maximal solution . If one of the eigenvalues of vanishes in the stencil (the problem is said to be resonant in this case), (17) may have no solution or have more than one.
3.2 Well-balanced property
The numerical method (14) is well-balanced in the sense given by the following
Observe that, in this case
The proof is trivial from this equality. ∎
3.3 Numerical method with WENO reconstructions
Let us discuss the implementation of the numerical method (14) in the particular case of WENO reconstructions with the upwind or the flux-splitting approach.
3.3.1 Upwind approach
The implementation in this case is as follows: once the solution has been computed:
Choose intermediate matrices .
where the projection matrices are given by (9).
3.3.2 Flux-splitting approach
The implementation of WENO with splitting approach will be as follows: once the solution has been computed:
In the particular case of the Lax-Friedrichs splitting, the reconstruction operators and will be applied to the values
in the corresponding stencil. Again is the local or global maximum eigenvalue of : although is not explicitly taken into account in the numerical viscosity, the numerical method has been shown to be stable under a CFL number of 1/2 in all the test cases considered here.
Observe that, while for a conservative system 2 reconstructions are computed at every stencil , 4 reconstructions have to be computed at every stencil: 2 using , 1 using , 1 using .
3.4 Numerical methods that preserve only one stationary solution
If there is only one known stationary solution to preserve (which is the case if, for instance, the initial condition is a perturbation of a given steady state), the following semidiscrete numerical method can be used:
where are computed now as follows:
Apply the reconstruction operator to the point values
Apply the reconstruction operator to the point values
It can be easily checked that this method preserves the stationary solution and, if this solution is smooth enough, the order of the method is equal to that of the reconstruction operator. Observe that the first stage of the reconstruction procedure, that consisted of the resolution of a Cauchy problem for a ODE system, is avoided; therefore the computational cost is expected to be lower than that of the methods that preserve any stationary solution.
4 Well-balanced high-order finite difference methods: discontinuous
Let us suppose now that is a.e. differentiable with finitely many isolated jump discontinuities. In this case, the definition of weak solutions (and, in particular, of stationary solutions) of (1) becomes more difficult: a solution is expected to be discontinuous at the discontinuities of and, in this case, the source term cannot be defined within the distributional framework. The source term becomes then a nonconservative product that can be defined in infinitely many different forms: see . Like in , the following criterion is assumed here to define the weak solutions of (1):
A pair of states can be the limits of an admissible weak solution of (1) to the left and to the right of a discontinuity point of , , if and only if there exists a solution of the ODE system
This admissibility criterion can be understood as a particular choice of paths within the theory developed by DalMasso, LeFloch, and Murat in  that allows one to give a sense to the nonconservative product as a Borel measure. In  it has been shown that the admissible jumps according to this criterion can be regularized through a smooth stationary solution in the following sense: is an admissible jump at a discontinuity point of if and only if for any given it is possible to find two smooth functions , defined in such that is monotone and:
Moreover, the jumps at the discontinuities of can be interpreted as stationary contact discontinuities of an extended system and the admissibility criterion (AC) is equivalent to assuming that Riemann invariants are preserved through these discontinuities (see  for details).
Let us remark that, if is a solution of (21), then is a stationary solution of (1): observe that (13) is trivially satisfied in smooth regions and (AC) is trivially satisfied in discontinuities. The well-balanced numerical method (14) can be easily adapted to this case by looking for stationary solutions of this form. More precisely, let us assume that the mesh has been designed so that the discontinuities of lie in intercells. Then, the numerical method (14) can be still applied with the only difference that now the first stage of the reconstruction step (15) is now as follows:
Look for the solution of (21) such that:
Every discontinuity point will be placed at an intercell . The reconstruction so defined leads to the following approximation of the Dirac delta produced by the discontinuity of :
where is a solution of (21): therefore, it is consistent with the assumed definition of the nonconservative products. Although consistency is not enough to guarantee the convergence to the right weak solution of nonconservative systems (see , ) the numerical tests in Section 5 show that the numerical methods capture the correct weak solutions.
In order to compare the behaviour of the different methods in the presence of a discontinuity of , for WENO methods (1) the Dirac delta issued from a discontinuity of will be approached by
where is some intermediate value of at the discontinuity. An upwind treatment of the singular source term is then used, so that the numerical method for the neighbor nodes writes as follows:
Here, are the projection matrices given by (9). It will be seen in Section 5 that these methods do not converge to the assumed weak solutions. Moreover, the behaviour of the numerical methods at the discontinuity of depends both on and the chosen intermediate state.
In the case of the methods that only preserve one stationary solution , the source term in the neighbor nodes of the discontinuity will be computed as follows:
Clearly, if for the right-hand sides vanish and the approximation of the Dirac mass is given again by (24) but, if it is not the case, a combination of (24) and (25) is used. Therefore, this method is only consistent with the definition of weak solution when the stationary solution is not perturbed at the neighbour nodes of the discontinuity, as it will be seen in Section 5
5 Numerical tests
In this section we apply the numerical methods introduced in Sections 2 and 3 to a number of test cases with continuous or discontinuous. Three families of methods based on WENO reconstructions of order will be compared:
WENO: methods of the form (3).
WBWENO: methods of the form (14) that preserve any stationary state.
WB1WENO: methods of the form (20) that preserve only one given stationary state.
Nevertheless, in many test cases the results obtained with WBWENO and WB1WENO are indistinguishable: in those cases, only the results corresponding to WBWENO will be shown. In all cases, the global Lax-Friedrichs flux-splitting approach is used for WENO implementation and the third order TVD-RK3 method is applied for the time discretization: see . The CFL parameter is set to 0.5.
5.1 A linear problem
We consider the linear scalar problem
In this case, (21) reduces to
whose solutions are
The stationary solutions of (28) for any given are thus given by:
The solution of (23) is thus
Therefore, well-balanced methods are based on the reconstructions of
5.1.1 Order test
Let us consider (28) with
It can be easily checked that the solution of (28) with initial condition:
is given by
Let us consider the initial condition:
where is the 11th degree polynomial
see Figure 1.
We solve (28) with initial condition (31) with well-balanced and non well-balanced third order methods in the interval . Free boundary conditions based on the use of ghost cells are used at both extremes. Table 1 shows the -errors and the empirical order of convergence corresponding to WENO and WBWENO, . As can be seen, both methods are of the expected order and the errors corresponding to methods of the same order are almost identical. In order to capture the expected order, the smooth indicators of the WENO reconstruction have been set to 0 and has been chosen for the fifth order methods.
5.1.2 A moving discontinuity linking two stationary solutions
Next, we consider (28) with again , and initial condition
The solution consists of a discontinuity linking two steady states that travels at speed 1:
Figure 2 shows the exact and the numerical solutions at time obtained with WBWENO3 and WENO3 (left) and a zoom of the differences of the numerical and the exact solution at the same time (right). It can be observed that the stationary states at both sides of the discontinuity are better captured with the well-balanced method. For fifth order methods the differences are lower, but still noticeable: see Figure 3. The results obtained with WBWENO, and the WB1WENO, that only preserve the stationary solution , are indistinguishable.
5.2 Burgers’ equation with source term
We consider next the scalar equation
The stationary solutions are also given by (30). The reconstruction operator has to be applied in this case to
5.2.1 Preservation of a stationary solution with smooth
In this test case we consider and the stationary solution
Let us solve (32) taking this stationary solution as initial condition in the interval . As boundary conditions, the value of the stationary solution is imposed at the ghost cells. Figures 4 and 5 show the differences between the stationary solution and the numerical solutions obtained at time with WENO and WBWENO, using a 200-cell mesh: the well-balanced methods capture the stationary solution with machine accuracy. This is confirmed by Tables 2 and 3 that show the -errors and the empirical order of convergence corresponding to WENO and WBWENO, .
5.2.2 Preservation of a stationary solution with oscillatory smooth
Let us consider now (32) with a function that has an oscillatory behavior:
(see Figure 6). We consider again the interval and we take as initial condition the stationary solution
(see Figure 6) and, as boundary conditions, the value of this stationary solution is also imposed at the ghost cells.
We consider a 100-cell mesh, so that the period of the oscillations is close to . Figure 7 shows the numerical solutions at time corresponding to WBWENO, WENO, : while the well-balanced methods preserve the stationary solution with machine precision, the non well-balanced methods give a wrong numerical solution. Of course, they give more accurate solutions if the mesh is refined: see next paragraph, where the reference solution is computed using WENO3.
5.2.3 Perturbation of a stationary solution with oscillatory smooth
(see Figure 8). If a mesh with 100 cells is used, the non well-balanced methods are unable to follow the evolution of the perturbation, since the numerical errors observed in the previous test are much larger than the perturbation. Let us see what happens when well-balanced methods are used: Figure 9 shows the numerical solutions obtained with WBWENO3 and WBWENO5 and a reference solution computed with WENO3 using a mesh of 5000 cells at time . As it can be seen, both methods are able to follow the evolution of the perturbation.
5.2.4 Preservation of a stationary solution with piecewise continuous
Let us consider now (32) with a piecewise continuous function :
(see Figure 10). We consider again the interval and we take as initial condition the stationary solution
(see Figure 10) and, as boundary conditions, the value of this stationary solution is also imposed at the ghost cells. WBWENO, preserve the stationary solution again with machine accuracy: see Table 4.
For WENO methods, the non-conservative product appearing at the source term is discretized by (26) with two different definitions of : a centered one
or an upwind one
Here, is the index of the intercell is located and is the arithmetic mean of and .
In Figure 11 we compare the exact solution with the numerical solutions obtained at time obtained with WBWENO3 using a 300-cell mesh (its graph and the one of the exact solution are identical at the scale of the figure) and with WENO3 or WENO5 using different implementations:
WENO3-UPW1: WENO3 with upwind implementation and (38);
WENO3-UPW2: WENO3 with upwind implementation and (37);
WENO3-LF: WENO3 with LF implementation and (38);
WENO5-LF: WENO5 with LF implementation and (38).
As it can be observed, the results of the WENO methods depends on the chosen implementation, on the numerical definition of the source term, and on the order. Moreover, these differences remain as tends to 0. This is in good agreement with the difficulties of convergence of finite difference methods to nonconservative systems: see .