Well-balanced high-order finite difference methods for systems of balance laws

01/27/2020
by   Carlos Parés, et al.
0

In this paper, high order well-balanced finite difference weighted essentially non-oscillatory methods to solve general systems of balance laws are presented. Two different families are introduced: while the methods in the first one preserve every stationary solution, those in the second family only preserve a particular stationary solution. The accuracy, well-balancedness, and conservation properties of the methods are discussed, as well as their application to systems with singular source terms. The strategy is applied to derive third and fifth order well-balanced methods for a linear scalar balance law, Burgers' equation with a nonlinear source term, and for the shallow water model.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 2

page 3

page 4

07/02/2020

An order-adaptive compact approximation Taylor method for systems of conservation laws

We present a new family of high-order shock-capturing finite difference ...
03/23/2020

High-order accurate entropy stable finite difference schemes for the shallow water magnetohydrodynamics

This paper develops the high-order accurate entropy stable (ES) finite d...
07/18/2020

An iterative scaling function procedure for solving scalar non-linear hyperbolic balance laws

The scaling of the exact solution of a hyperbolic balance law generates ...
05/18/2020

Entropy conservation property and entropy stabilization of high-order continuous Galerkin approximations to scalar conservation laws

This paper addresses the design of linear and nonlinear stabilization pr...
04/11/2020

High-order well-balanced finite volume schemes for hydrodynamic equations with nonlocal free energy

We propose high-order well-balanced finite-volume schemes for a broad cl...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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ázquez-Cendón introduced in [2] 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, [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 well-balanced high-order 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 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 [5] to develop well-balanced high-order finite volume numerical methods: see [11] 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 [6].

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 [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 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 [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 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 [23], [35], [36]. Uniform meshes of constant step will be considered and the following notation will be also used for the intercells

Once the reconstruction operator has been chosen, a possible semi-discrete numerical method for solving (1) can be written as follows:

(3)

where are the flux reconstructions:

(4)

The semidiscrete method (3) will be discretized in time by using TVD-RK methods: see [21].

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 Lax-Friedrichs flux-splitting:

where is the maximum of the absolute value of the eigenvalues of , this maximum being taken over either local (WENO-LLF) or global (WENO-LF): 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 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:

(13)

Once the reconstruction operator has been chosen, the following semi-discrete numerical method is proposed:

(14)

where are computed as follows:

  1. Look for the stationary solution such that:

    (15)
  2. Define

  3. Compute

The following result holds:

Proposition 1.

If the numerical method (3) is well-defined and the stationary solutions of (1) are smooth, then it has order of accuracy .

Proof.

Let be a smooth solution of (1). Given a time and an index , the reconstruction procedure is applied to to obtain , where represents the solution of (17) that satisfies

(16)

One has:

where the facts that is a solution of (1) and a stationary solution satisfying (16) have been used. ∎

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

(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 Well-balanced property

The numerical method (14) is well-balanced in the sense given by the following

Proposition 2.

Given a stationary solution of (1

), the vector of its point values

is an equilibrium of the ODE system given by the semi-discrete method (14).

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 flux-splitting 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 .

  • Define

    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:

  • Define

  • Compute

  • Define

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.

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:

  1. Apply the reconstruction operator to the point values

    to obtain

  2. 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 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 [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 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:

  1. Look for the solution of (21) such that:

    (23)

    where .

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 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 [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, 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:

(31)

where is the 11th degree polynomial

such that

see Figure 1.

Figure 1: Test 5.1.1: initial condition (left). Exact solution and numerical solution obtained with WBWENO3 and WBWENO5 at time using a mesh of 200 cells

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.

WENO3 WBWENO3 WENO5 WBWENO5
Cells Error Order Error Order Error Order Error Order
100 1.000E-1 - 1.023E-1 - 4.0902E-2 - 4.0910E-2 -
200 2.053E-2 2.28 2.084E-2 2.29 2.4404E-3 4.06 2.4407E-3 4.06
400 2.978E-3 2.78 3.019E-3 2.78 9.1307E-5 4.74 9.1315E-5 4.74
800 3.815E-4 2.96 3.867E-4 2.96 3.0118E-6 4.92 3.0121E-6 4.92
1600 4.788E-5 2.99 4.855E-5 2.99 9.4849E-8 4.98 9.4857E-8 4.98
Table 1: Test 5.1.1. Errors in norm and convergence rates for WB and WBWENO, at time .

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.

Figure 2: Test 5.1.2: exact solution and numerical solutions obtained with WBWENO3 and WENO3 using a mesh of 200 cells, (left); zoom of the differences between the numerical and the exact solutions (right)
Figure 3: Test 5.1.2: exact solution and numerical solutions obtained with WBWENO5 and WENO5 using a mesh of 100 cells, (left); zoom of the differences between the numerical and the exact solutions (right)

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

Figure 4: Test 5.2.1: zoom of the differences between the numerical solutions at time and the stationary solution using a 200-cell mesh. Left: WBWENO3. Right: WENO3
Figure 5: Test 5.2.1: zoom of the differences between the numerical solutions at time and the stationary solution using a 200-cell mesh. Left: WBWENO5. Right: WENO5
WENO3 WBWENO3
Cells Error Order Error
100 1.9044E-06 - 8.9928E-17
200 2.4762E-07 2.94 1.4543E-16
400 3.1550E-08 2.97 1.5304E-14
800 3.9817E-09 2.98 1.6560E-14
Table 2: Test 5.2.1. Errors in norm and convergence rates for WB3 and WBWENO3 at time .
WENO5 WBWENO5
Cells Error Order Error
20 7.7695E-07 - 2.2759e-16
40 3.5170E-09 7.78 1.5543e-16
80 2.0005E-10 4.13 1.1657e-16
160 1.0352E-11 4.27 2.9559e-16
Table 3: Test 5.2.1. Errors in norm and convergence rates for WENO5 and WBWENO5 at time .

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.

Figure 6: Test 5.2.2: graph of the function (left) and stationary solution (right)
Figure 7: Test 5.2.2: exact solution and numerical solutions obtained at time and a mesh of 100 cells. Left: WBWENO3, WENO3. Right: WBWENO5, WENO5

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

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

Figure 8: Test 5.2.3: initial condition. Left: graph. Right: difference with the stationary solution
Figure 9: Test 5.2.3: reference and numerical solutions obtained with WBWENO3 and WBWENO5 at time and a mesh of 100 cells. Left: graphs. Right: difference with the stationary solutions

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.4984E-15 5.3790E-14
200 1.5432E-16 1.7763E-16
300 2.1464E-16 4.1611E-15
Table 4: Test 5.2.4. Errors in norm for WBWENO3 and WBWENO5 at time .
Figure 10: Test 5.2.4: graph of the function (left) and stationary solution (right)

For WENO methods, the non-conservative 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 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 [7].

Figure 11: Test 5.2.3: exact and numerical solutions obtained at time with WBWENO3,WENO3-UPW1, WENO3-UPW2,WENO3-LF, WENO5-LF. Left: global view. Rigth: zoom close to the discontinuity

5.2.5 Perturbation of a stationary solution with piecewise continuous

We consider again (32) with (35) and an initial condition that is the stationary solution approximated in the previous test with a small perturbation