1 Introduction
We consider 1d systems of balance laws of the form
(1) 
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:
(2) 
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ázquezCendón introduced in [2] the condition called Cproperty: a scheme is said to satisfy this condition if it solves correctly the steadystate solutions corresponding to water at rest. Since then, the design of highorder wellbalanced numerical methods for systems of balance laws has been a very active front of research: see, for instance, [1], [4], [5], [8], [9], [10], [11], [12], [13], [16] [17], [19], [22], [24], [28], [29], [30], [31], [34], [39] [40], [41], [42]…
We focus here on finite difference methods. Previous wellbalanced highorder finite difference types have been based in techniques like the one proposed in [20] 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, [4], [16]). In [39][40], 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
one has
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 highorder 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 wellbalanced reconstruction introduced in [5] to develop wellbalanced highorder finite volume numerical methods: see [11] for a recent followup. 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 [6].
While the finite volume methods based on wellbalanced 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 [15], 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 RankineHugoniot 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 [11] 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 [32]), 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 wellbalanced property. Next, methods that preserve only one stationary solution are introduced and the conservation property is discussed. The implementation of wellbalanced 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 highorder 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. Finitedifference ENO or WENO reconstructions are examples of such reconstruction operators: see [23], [35], [36]. 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 :
(5)  
(6) 
Once these reconstructions have been computed, an upwind criterion can be chosen to define . For instance, for scalar problems
(7) 
an approximate value of is chosen and the numerical flux is defined then as follows:
where
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
(8) 
where
(9) 
Here is the diagonal matrix whose coefficients are
and
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 :
(10)  
(11) 
and finally,
(12) 
A standard choice is the LaxFriedrichs fluxsplitting:
where is the maximum of the absolute value of the eigenvalues of , this maximum being taken over either local (WENOLLF) or global (WENOLF): see [23], [36].
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 Wellbalanced highorder 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:
(13) 
Once the reconstruction operator has been chosen, the following semidiscrete numerical method is proposed:
(14) 
where are computed as follows:

Look for the stationary solution such that:
(15) 
Define

Compute
The following result holds:
Proposition 1.
Proof.
Observe that the method is welldefined if the first step of the reconstruction procedure can be always performed, i.e. if for every the Cauchy problem
(17) 
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.
Remark 1.
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
(18) 
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 .
Remark 2.
If the Jacobian of the flux function is regular, the ODE system can be written in normal form:
(19) 
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 Wellbalanced property
The numerical method (14) is wellbalanced in the sense given by the following
Proposition 2.
Proof.
Observe that, in this case
and thus
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 fluxsplitting approach.
3.3.1 Upwind approach
The implementation in this case is as follows: once the solution has been computed:

Define

Compute

Choose intermediate matrices .
3.3.2 Fluxsplitting approach
The implementation of WENO with splitting approach will be as follows: once the solution has been computed:

Define

Compute

Define
In the particular case of the LaxFriedrichs 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.
Remark 3.
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:
(20) 
where are computed now as follows:

Apply the reconstruction operator to the point values
to obtain

Apply the reconstruction operator to the point values
to obtain
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.
Another important difference between numerical methods (14) and (20) lies in their conservation properties: let us assume that there exists such that:
i.e. the first equations of (1) are conservation laws. Then, (20) is conservative for these equations: in effect, the first components reduce to the conservative expressions:
where and represent respectively the th component of and . This is not the case for since inequality (18) is expected to hold at every component for (14) due to the fact that, in general, .
4 Wellbalanced highorder 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 [15]. Like in [11], 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
(21) satisfying
(22)
This admissibility criterion can be understood as a particular choice of paths within the theory developed by DalMasso, LeFloch, and Murat in [15] that allows one to give a sense to the nonconservative product as a Borel measure. In [11] 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 [11] 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 wellbalanced 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:
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 :
(24) 
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 [33], [7]) 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
(25) 
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:
(26)  
where
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:
(27)  
with
Clearly, if for the righthand 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 LaxFriedrichs fluxsplitting approach is used for WENO implementation and the third order TVDRK3 method is applied for the time discretization: see [21]. The CFL parameter is set to 0.5.
5.1 A linear problem
We consider the linear scalar problem
(28) 
In this case, (21) reduces to
(29) 
whose solutions are
The stationary solutions of (28) for any given are thus given by:
(30) 
The solution of (23) is thus
Therefore, wellbalanced 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:
(31) 
where is the 11th degree polynomial
such that
see Figure 1.
We solve (28) with initial condition (31) with wellbalanced and non wellbalanced 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.
WENO3  WBWENO3  WENO5  WBWENO5  

Cells  Error  Order  Error  Order  Error  Order  Error  Order 
100  1.000E1    1.023E1    4.0902E2    4.0910E2   
200  2.053E2  2.28  2.084E2  2.29  2.4404E3  4.06  2.4407E3  4.06 
400  2.978E3  2.78  3.019E3  2.78  9.1307E5  4.74  9.1315E5  4.74 
800  3.815E4  2.96  3.867E4  2.96  3.0118E6  4.92  3.0121E6  4.92 
1600  4.788E5  2.99  4.855E5  2.99  9.4849E8  4.98  9.4857E8  4.98 
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 wellbalanced 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
(32) 
with
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
(33) 
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 200cell mesh: the wellbalanced 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, .
WENO3  WBWENO3  

Cells  Error  Order  Error 
100  1.9044E06    8.9928E17 
200  2.4762E07  2.94  1.4543E16 
400  3.1550E08  2.97  1.5304E14 
800  3.9817E09  2.98  1.6560E14 
WENO5  WBWENO5  

Cells  Error  Order  Error 
20  7.7695E07    2.2759e16 
40  3.5170E09  7.78  1.5543e16 
80  2.0005E10  4.13  1.1657e16 
160  1.0352E11  4.27  2.9559e16 
5.2.2 Preservation of a stationary solution with oscillatory smooth
Let us consider now (32) with a function that has an oscillatory behavior:
(34) 
(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 100cell 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 wellbalanced methods preserve the stationary solution with machine precision, the non wellbalanced 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
We consider again Burgers’ equation (32) with as in (34), and an initial condition that is the stationary solution approximated in the previous test with a small perturbation
(see Figure 8). If a mesh with 100 cells is used, the non wellbalanced 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 wellbalanced 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 :
(35) 
(see Figure 10). We consider again the interval and we take as initial condition the stationary solution
(36) 
(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.
Cells  WBWENO3  WBWENO5 

100  7.4984E15  5.3790E14 
200  1.5432E16  1.7763E16 
300  2.1464E16  4.1611E15 
For WENO methods, the nonconservative product appearing at the source term is discretized by (26) with two different definitions of : a centered one
(37) 
or an upwind one
(38) 
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 300cell 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:

WENO3UPW1: WENO3 with upwind implementation and (38);

WENO3UPW2: WENO3 with upwind implementation and (37);

WENO3LF: WENO3 with LF implementation and (38);

WENO5LF: 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 [7].
Comments
There are no comments yet.