R Eymard

T Gallouët

R Herbin

J.-C Latché

FINITE VOLUME SCHEMES AND LAX-WENDROFF CONSISTENCY

Keywords: 2010 AMS Subject Classification. Primary 65M08, 76N15 ; Secondary 65M12 Finite-volume schemes, Lax-Wendroff theorem

We present a (partial) historical summary of the mathematical analysis of finite difference and finite volume methods, paying special attention to the Lax-Richtmyer and Lax-Wendroff theorems. We then state a Lax-Wendroff consistency result for convection operators on staggered grids (often used in fluid flow simulations), which illustrates a recent generalization of the flux consistency notion designed to cope with general discrete functions.

The finite difference method and the Lax-Richtmyer theorem

In classical numerical analysis textbooks, we are taught that, in order to show the convergence of a finite difference (FD) scheme, one should show its stability and its consistency. The founding result in this regard is the Lax-Richtmyer theorem [START_REF] Lax | Survey of the stability of linear finite difference equations[END_REF] due to P.D. Lax, who exposed it in a seminar at NYU in 1954. The so-called Lax equivalence theorem can be summarized as follows (see for example [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations[END_REF]Theorem 1.5.1]): Theorem 1.1 (Lax-Richtmyer). Consider a linear partial differential equation for which the initial value problem is well posed, and a finite difference scheme consistent for its approximation; then this scheme is convergent if and only if it is stable.

In some articles and textbooks (see for instance [27, p. 142]), maybe because of the name "equivalence theorem", and because uniform meshes are considered, the Lax-Richtmyer theorem is replaced by the equivalence "consistency + stability ⇐⇒ convergence"

However, Theorem 1.1 does not state the equivalence [START_REF] Arakawa | A potential enstrophy and energy conserving scheme for the shallow water equations[END_REF], which, in fact does not hold in the general case of a PDE which is discretized with a non constant space step, even in the linear case. As an example, let us consider the approximation on R×]0, T [, where T > 0 is the final time of the study, of the linear transport equation

∂ t u(x, t) + ∂ x (au)(x, t) = 0, x ∈ R, t ∈]0, T [, (2)
u(x, 0) = u ini (x), x ∈ R, (3)
for some given a > 0, and initial data u ini ∈ C ∞ c (R, R), with ∂ t (resp. ∂ x ) the time (resp. space) partial derivative. Let ū(x, t) = u ini (x -at) be the exact (unique) solution of this problem. The FDM applied to (3) is classically defined by choosing a strictly increasing sequence (x i ) i∈Z of real numbers, such that h := max i∈Z (x i+1 -x i ) < ∞ and h = min i∈Z (x i+1 -x i ) > 0, and a time step δt = T /N , for N ∈ N with N > 1 (here a constant time step is considered for simplicity). The initial data is discretized by defining the following quantities that depend on h implicitly:

u 0 i = u ini (x i ) for any i ∈ Z. (4)
Since some upwinding is indeed necessary for stability purposes, the approximation of ∂ x u at the point x i is upwinded, so that the FD scheme reads

u n+1 i -u n i δt + a u n i -u n i-1 h i-1/2 = 0 for any i ∈ Z and n ∈ 0, N -1 , (5)
with h i-1/2 = x i -x i-1 . In this context, the definition of the terms consistency, stability and convergence is the following:

• Consistency: it requires two conditions:

the consistency of the discretization of the initial condition, that is

lim h→0 max i∈Z |u 0 i -u ini (x i )| = 0. (6)
the consistency of the discretization of the PDE (2); setting t n = nδt for n ∈ N, it reads:

lim h→0 δt→0 max i∈Z,n∈ 0,N -1 u(x i , t n+1 ) -u(x i , t n ) δt + a u(x i , t n ) -u(x n i-1 , t n ) h i-1/2 = 0. (7)
• Stability: There exists C depending only on u ini (and thus not on h nor on δt) such that max i∈Z,n∈ 0,N -1

|u n i | ≤ C. (8)
• Convergence:

max i∈Z,n∈ 0,N -1 |u n i -u(x i , t n )| → 0 as h → 0.
Clearly, condition (4) ensures the consistency condition (6) on the initial condition. Moreover, the consistency condition [START_REF] Diperna | Measure-valued solutions to conservation laws[END_REF] on the PDE can be obtained by Taylor expansions. Finally, the stability condition (8) can be shown under the CFL condition aδt ≤ h. The conditions for the Lax-Richtmyer to hold are therefore satisfied, and the convergence of the scheme (4)-( 5) is thus proven.

Consider now a variant of the scheme (5) obtained by keeping (4), but replacing

h i-1/2 by h i = x i+1 -x i-1 2 in (5): u n+1 i -u n i δt + a u n i -u n i-1 h i = 0 for any i ∈ Z and n ∈ 0, N -1 , (9)
In the specific case where x 2k+1 -x 2k = h/2 and x 2k+2 -x 2k+1 = h for all k ∈ Z, we get that h i = 3 4 h for all i ∈ Z and the consistency property no longer holds. Therefore, if the equivalence (1) were true, the scheme (4)- [START_REF] Eymard | Existence and uniqueness of the entropy solution to a nonlinear hyperbolic equation[END_REF] would not be convergent in the above sense. However, let us show that this scheme is in fact convergent in the same sense. Let us write the finite difference scheme obtained at the points x i = xi+xi+1 2 instead of the points x i , defining the values u n i :

u 0 i = u ini ( x i ) for any i ∈ Z, (10)
and

u n+1 i -u n i δt + a u n i -u n i-1 h i = 0 for any i ∈ Z and n ∈ 0, N -1 . (11)
Then the scheme ( 10)-( 11) again satisfies the consistency condition [START_REF] Crandall | Monotone difference approximations for scalar conservation laws[END_REF], and satisfies the consistency condition (7) since h i = x i -x i-1 . The stability condition [START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF] can again be shown under the CFL condition aδt ≤ h ≤ 3 4 h. Hence the Lax-Richtmyer theorem yields the convergence of the values

u n i to u( x i , t n ). Using | u 0 i -u 0 i | ≤ h max |u ′ ini |
, the maximum principle applied to the difference u n i -u n i , solution to the equation obtained by subtracting [START_REF] Eymard | Existence and uniqueness of the entropy solution to a nonlinear hyperbolic equation[END_REF] to [START_REF] Faille | Des mathématiciens découvrent les volumes finis[END_REF], shows that max i∈Z,n∈ 0,N -1

u n i -u n i ≤ h max |u ′ ini |, which implies that max i∈Z,n∈ 0,N -1 |u n i -u(x i , t n )| ≤ 2h max |u ′ ini | + max i∈Z,n∈ 0,N -1 | u n i -u( x i , t n )|,
leading to the convergence of the scheme in the same sense as above. This example shows that on the one hand the direction ⇐ of ( 1) is not true, and that on the other hand the Lax-Richtmyer theorem cannot be applied directly to obtain the convergence of such a scheme, since the scheme is not consistent in the above defined sense. In fact, the scheme ( 9) is the 1D upwind FV scheme with the control volumes

] xi-1+xi 2 , xi+xi+1 2
[, as shown in the next section. In the scheme (5), the partial derivative ∂ x u is upwinded. In the scheme (9), the unknown u itself is upwinded.

Note that the equivalence (1) also does not hold in the case of an elliptic operator, see the example of a non uniform 1D mesh in [START_REF] Faille | Des mathématiciens découvrent les volumes finis[END_REF] or [10, section 5.2]. For the analysis of such schemes, other notions must be introduced, and it seems that Peter D. Lax first identified them, in collaboration with B. Wendroff.

The finite volume method and the Lax-Wendroff theorem

In a very famous article of 1960, P.D. Lax and B. Wendroff [START_REF] Lax | Systems of conservation laws[END_REF] consider discretization schemes for nonlinear hyperbolic systems of conservation laws and show that if a conservative scheme with consistent fluxes, in the sense that they define and which is stated below (see Section 2.2), converges a.e. and boundedly towards a limit, then this limit is necessarily a weak solution of the system. We call this property Lax-Wendroff (LW) consistency. The notions of flux consistency and flux conservativity highlighted in [START_REF] Lax | Systems of conservation laws[END_REF] are truly fundamental for the convergence analysis of the FVM for the hyperbolic equations considered in [START_REF] Lax | Systems of conservation laws[END_REF], as are their extensions to elliptic and parabolic conservation equations. In order to explain these terms, we consider, again in the 1D case, a general differential form of a conservation law, written on the whole space R and on a time interval ]0, T [ where 0 < T < +∞ is the final time:

∂ t u(x, t) + ∂ x F(x, t) = 0, (12)
stating the conservation of the quantity u at each point x ∈ R and each time t ∈]0, T [, with F a vector function depending only on x and t through the unknown u. In addition to this equation, an initial condition is assumed to be given on u. Simple examples of such a conservation law include • the transport equation F(x, t) = au(x, t) (linear hyperbolic equation, introduced in the previous section), • the heat equation F(x, t) = -∂ x u(x, t) (linear parabolic equation),

• Burgers' equation F(x, t) = u 2 (x, t) (nonlinear hyperbolic equation),

• the porous media equation F(x, t) = -∂ x u p (x, t), p > 1 (nonlinear parabolic equation). The FVM consists in approximating the integral form of the conservation law, that is to say the balance on the time-space cuboid ]x, x + δx[×]t, t + δt[ (for given δx > 0 and δt > 0), rather than the PDE itself (this corresponds to the way such an equation is derived from physical conservation principles). Note that with a non zero source term on the right hand-side of (12), the equation is then usually renamed "balance law", and this modification poses no problem for a FV discretization. The integral form relative to the differential form [START_REF] Gallouët | Mesure,intégration, probabilités[END_REF] of a conservation law reads

x+δx x (u(x, t + δt) -u(x, t)) dx + t+δt t F(x + δx, t) -F(x, t) dt = 0. (13)
Let (]x i-1/2 , x i+1/2 [) i∈Z be a family of intervals of R (also called control volumes or grid cells), with x i-1/2 < x i+1/2 and such that ∪ i∈Z ]x i-1/2 , x i+1/2 [= R (or Ω ⊂ R), and let δt = T /N , N > 1 (the time step could be taken non constant). Let x i be some chosen point in the cell ]x i-1/2 , x i+1/2 [ (the choice of this point is constrained by the flux consistency property in the case of elliptic or parabolic problems, but not in the hyperbolic case, contrarily to the FD consistency in the sense of the previous section). The discrete unknowns, {u n i , i ∈ Z, n ∈ {0, . . . , N }} are expected to be approximations of u(x i , t n ) with t n = nδt. The integral form ( 13) is written on each control volume ]x i-1/2 , x i+1/2 [ and time interval ]t n , t n+1 [, leading to the following FV scheme (with an explicit time scheme for the sake of simplicity):

h i u n+1 i -u n i δt + F n i+ 1 2 -F n i-1 2 = 0 for any i ∈ Z and n ∈ 0, N -1 , (14)
where

h i = x i+ 1 2 -x i+ 1 2 , and F n i+ 1 2
is the numerical flux, which will be expressed in terms of the discrete unknowns to yield a numerical approximation of F(x i+1/2 , t n ). Note that F n i+ 1 2 is the numerical flux outgoing ]x i-1/2 , x i+1/2 [ to the right and its opposite is the numerical flux outgoing ]x i+1/2 , x i+3/2 [ to the left: this copies the situation of the exact flux F(x i+1/2 , t n ). This is the well known "local conservativity" or "flux conservativity" property, which is important in physical applications, but also fundamental in the mathematical analysis of the FVM. Indeed, it is thanks to this property that we may hope to prove some convergence properties of the method, both for elliptic or parabolic type and for hyperbolic equations despite the loss of the consistency in the FD sense, as presented in the first section of this paper.

Flux conservativity

Writing the FV scheme in one dimension naturally ensures the numerical flux conservativity, since only one flux is defined: F i+ 1 2 at the interface x i+ 1 2 . In a multi-dimensional framework, (d = 2 or 3) the PDE ( 12) is written as ∂ t u + divF = 0, where F is a vector function of x and t and div the space divergence operator. The scheme ( 14) is now written for a control volume K:

|K| u n+1 K -u n K δt + σ⊂∂K |σ|F n K,σ = 0, n ∈ 0, N -1 , (15)
where |K| (resp. |σ|) is the volume or surface of K (resp. the surface or length of |σ|) and |σ|F n K,σ is the numerical outgoing flux from K through the face σ; it is an approximation of the ougoing normal flux

σ F(x, t n ) • n K,σ (where n K,σ is the normal vector to σ outward K), which is expressed in terms of the discrete unknowns (u n M ) M∈M : F n K,σ = F n K,σ ((u n M ) M∈M ).
In 2 or 3D, the numerical flux is defined on either side of the interface σ. Suppose that the interface σ separates the control volumes K and L, which we write as σ = K|L, then the flux conservativity reads

F n K,σ = -F n L,σ . (16)

Flux consistency

Let us turn back to the 1D case, for ease of notations. The numerical flux is said to be consistent if for a sequence of time and space discretizations, indexed by m and such that h (m) → 0 and δt (m) → 0 as m → +∞, one has max i∈Z,n∈ 0,N -1

|F(x i+ 1 2 , t n ) -F n i+ 1 2 | → 0 as m → +∞, where F n i+ 1 2
is the quantity obtained from F n i+ 1 2 when replacing the discrete unknowns by the values of a regular function u:

F n i+ 1 2 = F n i+ 1 2 (u(x i , t n )) i∈Z .
In the context of nonlinear hyperbolic equations,

F(x i+ 1 2 , t n ) = f (u(x i+ 1 2 , t n )) with f ∈ C(R, R
), this definition of consistency is equivalent to the usual notion of consistency introduced by Lax (written here for a two-point scheme):

F n i+ 1 2 (u, u) = f (u) if one assumes F n i+ 1 2
to be Lipschitz continuous, or at least, "lip-diag" see [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF]Remark 5.2]. In the context of the heat equation (F (u) = -∇u), the numerical flux

F n i+ 1 2 = F n i+ 1 2 (u n i , u n i+1 ) = - u n i+1 -u n i h i+ 1 2 with h i+ 1 2 = x i+1 -x i , (17)
is consistent in the above sense, since, for a regular function u,

F(x i+ 1 2 , t n ) = -∂ x u(x i+ 1 2 , t n ) and F n i+ 1 2 (u(x i , t n ), u(x i+1 , t n )) = - u(x i+1 , t n ) -u(x i , t n ) h i+ 1 2 = -∂ x u(x i+ 1 2 , t n ) -h i+ 1 2 ∂ 2 xx u(c, t n ), c ∈]x i , x i+1 [.
We notice that in this case the flux depends on the choice of the points x i .

Stability, compactness, convergence

Stability.

There are different notions of stability of a numerical scheme. The notion which is of interest in the context of the convergence of a numerical scheme for a general, possibly non linear, PDE is an estimate on the approximate solutions, independently of the mesh. For instance, the L ∞ stability of a linear FD scheme for a linear elliptic equation may be obtained by writing the scheme in matrix form and by obtaining a bound of the infinity norm of the inverse of this matrix. Even though the FV approximate solutions are piecewise constant, the estimates are obtained in a norm which is in close relation with that which one uses for the estimates on the solutions of the continuous problem, and which depends of course on the considered problem; the notion of stability in the FVM is therefore linked to the stability of the continuous problem.

Let us take two examples:

• the heat equation on [0, 1] with homogeneous Dirichlet conditions: the natural norm for the continuous problem is L 2 (H 1 0 ), and the associated discrete norm corresponding to the choice ( 17) is

L 2 (H 1 0,d ) with u H 1 0,d = M i=0 h i+1/2 u i+1 -u i h i+1/2 2 1 2 ,
where M is the number of control volumes and u 0 = u M+1 = 0 and with h i+1/2 = x i+1 -x i . We denote by u the weak solution of the heat equation and u app the solution of the time-implicit scheme (implicit schemes are a natural choice for parabolic equations to avoid a condition of type δt ≤ Ch 2 ). The function u app is a piecewise constant function which is equal to u n i on the cuboid ]x i-1/2 , x i+1/2 [×]t n , t n+1 [. In the continuous PDE setting, an L 2 (H 1 0 ) estimate on u is obtained by taking u as a test function in the weak formulation of the heat equation, and integrating by parts. Similarly, the L 2 (H 1 0,d ) estimate on u app is obtained by multiplying the i-th discrete equation by δt u n i , summing over i and n and performing discrete integrations by parts (obtained by changes of indices in the sums). An L 2 (L 2 ) estimate on the approximate solutions then follows with a discrete Poincaré inequality [10, Lemma 9.1].

• the transport equation on R n , with initial condition u ini ∈ L ∞ : the natural norm for the continuous problem is L ∞ . It is classical and easy to show that the scheme [START_REF] Gallouët | Weak consistency of finite volume schemes for systems of non linear conservation laws: Extension to staggered schemes[END_REF], with the upwind flux F i+1/2 = a u i , is stable, see for example [START_REF] Eymard | Finite Volume Methods[END_REF]Lemma 20.1], under a CFL (for Courant Friedrichs Lewy, [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF]) condition.

3.2. The linear case: stability + conservativity + flux consistency =⇒ convergence

Convergence for linear operators are often obtained through error estimates (the same technique as in the proof of stability is applied to the error between the approximate solution and the exact one); the compactness analysis of sequences of approximate solutions is another means, which also gives the existence of a solution (see [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF] for a seminal paper on this type of proof). Let us detail this second means, which extends to non linear operators. Consider a sequence of approximate solutions, on meshes whose space and time steps tend to 0. The general principle of proof of convergence is then the following. Thanks to the stability, a uniform (with respect to the discretization parameters) estimate on the approximate solutions holds in a Lebesgue space, and so there exists a sub-sequence of this sequence which converges weakly (or ⋆-weakly) in that same Lebesgue space. Each equation of the scheme is then multiplied by the interpolation of a regular test function and by the time step; the summation of all these equations over the time index and space indices, and an integration by parts (using the conservativity of the fluxes) are performed so as to shift the discrete derivatives from the discrete unknowns to the regular test function (the flux consistency which is used on the test functions is thus that given by the dual discrete operator). Since the problem is linear, the weak convergence of the approximate solutions suffices to pass to the limit in all terms; however the parabolic and hyperbolic cases exhibit some different difficulties.

-Elliptic or parabolic case. The passage to the limit can be performed thanks to the consistency of the flux, see e.g. [10, Theorem 8.1] for the 1D case. In the multidimensional case, the consistency of the flux for the heat equation (Laplace operator) is obtained for meshes which respect an orthogonality condition [10, Definition 9.1]. In this case, the resulting matrix is symmetric and the flux consistency is identical to that given by the dual discrete operator. Anisotropic operators and general meshes have been the object of several different works in the last decades, we refer to [START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF] for a review of these methods, several of them leading to non symmetric matrices. In this latter case, one can either prove convergence by an error estimate, using the flux consistency of the primal discrete operator, or a compactness method, using the flux consistency of the dual discrete operator, i.e.on the test functions. In both cases the main difficulty is to establish the stability of the scheme.

-Hyperbolic case. For stability reasons, the numerical flux is upwinded and introduces an error term whose convergence to 0 must be shown. This term is a sum of products of differences of the values of the solutions in neighbouring cells by discrete derivatives of the test function. On 1D or multiD Cartesian uniform or non uniform meshes, if u ini ∈ L ∞ ∩ BV (an additional argument is used to handle the case u ini ∈ L ∞ ) this term is shown to tend to 0 thanks to a uniform BV estimate on the approximate solutions (and on the continuous solution), see [START_REF] Crandall | Monotone difference approximations for scalar conservation laws[END_REF][START_REF] Kuznecov | Monotone difference approximations for a first order quasilinear equation[END_REF] who consider the non linear case; this proof uses the "TVD" character (total variation diminishing) of the monotone schemes [START_REF] Harten | On a class of high resolution total-variation-stable finite-difference schemes[END_REF]. Unfortunately, in the case of an unstructured mesh in multiD, even for a linear equation, even if u ini ∈ L ∞ ∩ BV , it can be shown that the upwind scheme is not TVD; a counter-example is given in [START_REF] Champier | Convergence d'un schéma décentré amont sur un maillage triangulaire pour un problème hyperbolique linéaire[END_REF]. A "weak BV" estimate based on the numerical diffusion of the scheme is established therein in order to prove convergence.

Let us notice that in both the parabolic and hyperbolic cases, the Lax-Wendroff theorem is not used directly:

• in the parabolic case, because the continuous flux function cannot be applied to the approximate solutions, • in the hyperbolic case, because the approximate solutions only converge weakly, whereas the Lax-Wendroff theorem supposes a strong convergence.

Nevertheless in both cases, the fundamental notions introduced in [START_REF] Lax | Systems of conservation laws[END_REF] are used, namely:

• flux conservativity: it is the property that leads to a weak form of the FV scheme by shifting the discrete derivatives of the discrete unknown to the discrete derivatives of the interpolant of the test function,

• flux consistency: it is the property that is used to prove the fact that a limit of approximate solutions is a weak solution thanks to the convergence of the discrete derivatives of the interpolant of the test function to the exact derivatives of the same test function.

The nonlinear case: more compactness needed

In the nonlinear case, weak convergence is not sufficient; indeed, if a sequence (u n ) n∈N converges only weakly in a Lebesgue space to a limit u, there is no reason why the sequence (f (u n )) n∈N should converge to f (u), even weakly.

-Elliptic and parabolic equations. In the elliptic or parabolic setting, one of the fundamental tools to obtain more compactness is Kolmogorov's compactness theorem, a consequence of which is that any bounded sequence of L p , 1 ≤ p < +∞ which is "equicontinuous on average" admits a convergent sub-sequence (see for example [START_REF] Gallouët | Mesure,intégration, probabilités[END_REF]Theorem 8.16]). Equicontinuity in mean amounts to showing that the difference between the function and its translates in time and space converges to 0 in the L p norm, uniformly with respect to the time and space step, see for example [START_REF] Eymard | Finite Volume Methods[END_REF]Lemma 18.3] in the case of a nonlinear parabolic equation of Stefan type. Once the compactness of the sequence of approximate solutions in L 2 (L 2 ) is proven, we can then exhibit a sub-sequence tending to ū in L 2 (L 2 ).

By passing to the limit in the "very weak" form of the scheme (that is with the discrete divergence of the discrete normal derivatives of the test functions) we can then show, as in the linear case, that each term (in time and space) converges to the corresponding term in the "very weak" formulation of the continuous problem. In the case of the Stefan problem, namely ∂ t u -∆ϕ(u) = 0, this gives the convergence of the approximate solutions to the exact solution if ϕ is an increasing Lipschitz continuous fucntion. There is however an additional difficulty if ϕ is only nondecreasing. In this case, it is possible to prove compactness in L 2 (L 2 ) of ϕ(u app ) (where u app is the approximate solution) but not of u app for which only an L 2 (L 2 ) bound holds. It is possible however to conclude using the Minty trick (see [START_REF] Eymard | Finite Volume Methods[END_REF]Chapter 4]).

-Hyperbolic equations. Now consider a nonlinear hyperbolic conservation law of the form

∂ t u(x, t) + div(f (u(x, t))) = 0, x ∈ R n , ( 18
) u(x, 0) = u ini (x). ( 19
)
If f ∈ C 1 (R, R n ) and u ini ∈ L ∞ , there exists a unique entropy solution of this problem [START_REF] Kružkov | First order quasilinear equations with several independent variables[END_REF]. In order to show that a scheme approximates this entropy solution, it is first shown that it satisfies a discrete entropy equation.

The case of Cartesian meshes has been studied independently by Kuznetsov [START_REF] Kuznecov | Monotone difference approximations for a first order quasilinear equation[END_REF] and Crandall and Majda [START_REF] Crandall | Monotone difference approximations for scalar conservation laws[END_REF]. As in the linear case, if u ini ∈ BV , a BV estimate on the approximate solutions holds, uniformly with respect to the space and time step; Helly's lemma, which is itself a direct consequence of Kolmogorov's compactness theorem may then be invoked to obtain the convergence of a subsequence of the approximate solutions in L 1 (L 1 ), and one can then use the Lax-Wendroff theorem (which generalizes easily to the entropy formulation) [START_REF] Crandall | Monotone difference approximations for scalar conservation laws[END_REF] [10, section 21.5]). It is also possible to handle the case u ini ∈ L ∞ , using a contraction principle in L 1 for the exact solution and for the approximate solution, see [START_REF] Crandall | Monotone difference approximations for scalar conservation laws[END_REF] for instance.

In the general case of a non-Cartesian mesh, a suitable BV estimate seems out of reach, and the proof of convergence is performed with the following steps.

• Consider a sequence of space discretizations, indexed by m, and of time steps δt (m) ; for any m, the mesh size is defined as the maximum diameter of the cells and denoted by h (m) . Assume that h (m) → 0 and δt (m) → 0 as m → +∞. • Owing to the L ∞ estimate on the sequence of approximate solutions (u (m) ) m∈N , there exists a subsequence which converges in a "non linear weak sense", i.e. there exists

µ ∈ L ∞ (R × R + ×]0, 1[) such that for any function ψ ∈ C(R, R), R R+ ψ(u (m) )ϕ → R R+ µ ψ ϕ, ∀ϕ ∈ L 1 (L 1 ) with µ ψ = 1 0 ψ µ(x, t, α) dα.
This notion of convergence is equivalent to the convergence to a Young measure [START_REF] Diperna | Measure-valued solutions to conservation laws[END_REF]; it may seem a little simpler to handle in the sense that it involves a function, µ, rather than a measure ; however this function depends on an additional parameter α ∈]0, 1[, which we will have to get rid of in order to reach convergence to an entropy weak solution. • Using the numerical diffusion of the scheme, we get a uniform weak BV estimate on the sequence of approximate solutions; this estimate is called "weak" for two reasons: on the one hand it involves the differences |(f

(u K ) -f (u L ))
• n KL | and not the differences |u K -u L |, where (K, L) denotes a pair of control volumes sharing a common interface ; on the other hand, it only requires that the sum of these differences does not blow up too fast: the difference between the discrete gradient of the interpolated test function and the gradient of the test function itself behaves like the mesh size, and to pass to the limit, we only need that the sum of the differences involving the discrete unknowns be bounded by a term in C/h 1-ε with ε > 0 and C > 0 independent of h. • Using the nonlinear weak convergence and the weak BV estimate, we pass to the limit on the weak form of the discrete entropy and obtain a so-called "process solution" which is an entropy solution up to an integral with respect to the additional parameter α. • Starting from the discrete entropy inequalities that are verified by the scheme, the process solutions are shown to satisfy an entropy inequality. A uniqueness result on the process solutions can then be obtained thanks to a variable doubling technique "à la Kruskov" [START_REF] Eymard | Existence and uniqueness of the entropy solution to a nonlinear hyperbolic equation[END_REF]; this result differs from Di Perna's [START_REF] Diperna | Measure-valued solutions to conservation laws[END_REF] in that it takes into account the initial condition in the weak entropy formulation, which allows to avoid the more restrictive conditions on the mesh [START_REF] Cockburn | Convergence of the finite volume method for multidimensional conservation laws[END_REF]. The uniqueness of the process solution entails that a process solution is the unique entropy weak solution. It also yields the (strong) convergence, in L p -spaces, of the approximate solution to the exact solution. The proof of convergence of the FVM has been obtained for several other problems than the one considered here. However, for systems of PDEs, it is often difficult to obtain compactness results, and LW-consistency then seems an interesting way to make sure that an eventual limit of the scheme is indeed a weak solution of the system. In the next section, we show how this is feasible even on staggered grids, which are often used in the numerical simulation of fluid flows.

LW-consistency and staggered grids

Rectangular staggered grids have been used since the sixties in fluid mechanics (in the wellknown MAC method [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF]) including environmental flows [START_REF] Arakawa | A potential enstrophy and energy conserving scheme for the shallow water equations[END_REF], see also [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF]. The mathematical analysis of the MAC scheme has been the object of several recent works, see e.g. [START_REF] Gallouët | Convergence of the Marker-And-Cell scheme for the incompressible Navier-Stokes equations on non-uniform grids[END_REF][START_REF] Gallouët | Error estimates for the implicit MAC scheme for the compressible Navier-Stokes equations[END_REF]. Systems of partial differential equations for which no existence or uniqueness is known have also been discretized on staggered grids. Let us mention in particular the compressible Euler equations: one of the advantages of using a staggered grid is to produce a scheme that can be mathematically proven to be asymptotically stable to the incompressible limit, as shown for the isentropic case in [START_REF] Herbin | Low Mach number limit of some staggered schemes for compressible barotropic flows[END_REF]. For such systems, a Lax-Wendroff type theorem is very useful; indeed, although true convergence cannot be proven for lack of compactness properties, such a theorem allows to state that if the scheme converges, and provided some bounds on the approximate solutions are satisfied (these bounds are generally not attainable mathematically but can be verified numerically), the limit of the scheme is an entropy weak solution of the system, see [START_REF] Herbin | Consistent segregated staggered schemes with explicit steps for the isentropic and full euler equations[END_REF], [START_REF] Herbin | Conservativity and weak consistency of a class of staggered finite volume methods for the Euler equations[END_REF]. The proof of this result may be obtained thanks to the generalization of the Lax-Wendroff theorem to general grids, which include staggered grids such as the MAC grid [START_REF] Gallouët | Weak consistency of finite volume schemes for systems of non linear conservation laws: Extension to staggered schemes[END_REF]. One of the major additional difficulties for staggered grids is that the discrete unknowns are piecewise constant on different grids.

Consider a rectangular domain Ω ⊂ R 2 (the 3D case can be tackled in the same way), and a possibly non-uniform rectangular grid. We denote by F the set of edges of the mesh, and the internal edge separating the cells K and L is denoted by σ = K|L (see Figure 1). This mesh will be referred to in the following as the primal mesh, and denoted by P. Two dual meshes (three in 3D) are considered, each consisting in a partition of Ω indexed by the vertical and horizontal elements of F, i.e. Ω = ∪ σ∈F (i) D σ , i = 1, 2, where F (1) (resp. F (2) ) denotes the set of vertical (resp. horizontal) edges. The cells (D σ ) σ∈F are referred to as the dual cells. A half dual cell D K,σ is half of the rectangle K with side σ (see Figure 1). For an internal edge σ = K|L, the dual cell D σ is the subset of K ∪ L defined as D σ = D K,σ ∪ D L,σ ; for an external edge σ of a cell K, D σ is the subset D K,σ of K.

M N K L σ = M|L σ ′ = N|K τ = M|N τ ′ = K|L D σ D ′ σ M N K L σ ′ σ D τ D τ ′ M N K L τ τ ′ Figure 1.
Primal and dual meshes and associated notations for the MAC case.

-Left: the primal cells; the edges σ and σ ′ belong to F (1) and the edges τ and τ ′ to F (2) .

-Center: the dual cells associated to F (1) .

-Right: the dual cells associated to F (2) .

To illustrate the use of the generalized Lax-Wendroff theorem proven in [START_REF] Gallouët | Weak consistency of finite volume schemes for systems of non linear conservation laws: Extension to staggered schemes[END_REF], let us consider as a simple example the mass equation of, say, the compressible Euler equation:

∂ t ρ(x, t) + div ρu (x, t) = 0 (x, t) ∈ Ω×]0, T [, (20)
where ∂ t ρ denotes the time derivative of the density ρ, and div the space divergence. The scalar unknown ρ is associated to the primal cells:

ρ(x, t) = ρ n K for x ∈ K, K ∈ P, t ∈ [t n , t n+1 [, n ∈ 0, N -1 .
The unknowns associated to the i-th component of u are located at the center of the edges of the i-th dual mesh. The associated approximate vector function thus reads: u(x, t) = (u 1 (x, t), u 2 (x, t)) t where, for i = 1, 2, u i (x, t) = u n σ , for x ∈ D σ , σ ∈ F (i) and t ∈ [t n , t n+1 [, n ∈ 0, N -1 . Let e (i) denote the i-th unit vector; the discretization of (20) reads:

C(ρ, u) n K = (ð t ρ) n K + 1 |K| σ∈F(K) |σ| F n σ • n K,σ = 0, with (ð t ρ) n K = ρ n+1 K -ρ n K t n+1 -t n and F n σ = ρ n σ u n σ , where u n σ is defined as u n σ e (i) for σ ∈ F (i) , i = 1 or 2,
and, for σ = K|L, ρ n σ stands for a convex combination of ρ n K and ρ n L (for instance the upwind or a MUSCL choice with respect to u n σ ). The initial value for the scalar unknown ρ is defined by

ρ 0 K = 1 |K| K ρ 0 (x) dx. (21)
For i = 1, 2, let h(i) = max{|σ|, σ ∈ F (i) } and h (i) = min{|σ|, σ ∈ F (i) }. We define the space step by h(P) = max( h(1) , h(2) ), and the time step by δt = max n∈ 0,N -1 (t n+1 -t n ). A sequence of grids is said quasi-uniform if the quotients h(1) /h (2) and h(2) /h (1) are bounded by a constant independent of the grid. Lemma 4.1 (Lax-Wendroff consistency for the mass equation, MAC grid). Let a sequence of quasi-uniform MAC grids (P (m) ) m∈N and of time discretizations be given, with h(P (m) ) and δt (m) tending to zero; let (ρ (m) , u (m) ) m∈N be the associated sequence of discrete functions.

We suppose that the sequences (ρ (m) ) m∈N and (u (m) ) m∈N are bounded in L ∞ (Ω×]0, T [) and L ∞ (Ω×]0, T [) 2 respectively, and that, when m tends to +∞, they converge in L p (Ω×]0, T [) and L p (Ω×]0, T [) 2 , 1 ≤ p < +∞, to ρ ∈ L ∞ (Ω×]0, T [) and u ∈ L ∞ (Ω×]0, T [) 2 respectively. Then (ρ, u) is a weak solution of [START_REF] Herbin | Low Mach number limit of some staggered schemes for compressible barotropic flows[END_REF], in the sense that, for any function ϕ ∈ C ∞ c (Ω × [0, T [), ρ(x, t) ∂ t ϕ(x, t) + ρu (x, t) • ∇ϕ(x, t) dx dt as m → +∞, where ϕ n K stands for the mean value of ϕ(x, t n ) over the cell K, and so the right hand-side of this assertion vanishes.

The proof of Lemma 4.1 is given in [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF] and quite simple, especially using the ad hoc tools developed in [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF][START_REF] Gallouët | Weak consistency of finite volume schemes for systems of non linear conservation laws: Extension to staggered schemes[END_REF]. However, the developed arguments may be extended to more complex operators, to deal for instance with the momentum and energy balance equations of the compressible Euler equations, even if the proofs are more tricky.

Conclusion

In this paper, we have presented some concepts for the mathematical analysis of FV schemes, paying a special attention to the flux consistency issue and its most direct consequence, i.e. the LW-consistency of the scheme. We emphasize that we gave here a very partial picture of the mathematical world of finite volumes. In many problems (some of them evoked here), the convergence of the discrete solution to a limit follows by compactness arguments in norms strong enough to "feed" the consistency study; in fine, this yields a stronger result than just the scheme consistency, namely the convergence (up to a subsequence if the uniqueness of the solution of the continuous problem is not known) of the numerical solutions to the (a) continuous one (see e.g. [START_REF] Gallouët | Convergence of the Marker-And-Cell scheme for the incompressible Navier-Stokes equations on non-uniform grids[END_REF]). Many parabolic equations enter this framework, including, focusing on fluid flow simulations, incompressible, possibly variable density, or steady barotropic Navier-Stokes equations. In some problems, the continuous solution may be reasonably supposed (or even proven) to be regular, and an error analysis is possible.

The LW-consistency issue is especially important for practical applications in fluid flow simulations. indeed, in many cases of interest, stronger results are out of reach, and this property is the only one left to mathematically support the design of schemes. This is for instance the case for multi-dimensional flows governed by hyperbolic systems, as shallow-water equations, Euler equations or models for multi-phase flows. For instance, the study evoked in Section 4 is motivated by such a situation: in the last ten years, a class of staggered schemes has been designed for hyperbolic flow problems [START_REF] Herbin | Conservativity and weak consistency of a class of staggered finite volume methods for the Euler equations[END_REF][START_REF] Herbin | Consistent segregated staggered schemes with explicit steps for the isentropic and full euler equations[END_REF][START_REF] Herbin | A decoupled staggered scheme for the shallow water equations[END_REF], and implemented in the open-source software CALIF 3 S developed at IRSN [2]; they are now routinely used for industrial safety applications as hydrogen explosion problems, supposing inviscid or at least vanishing viscosity flows. The accuracy of the numerical schemes involved here is essentially supported by LW-consistency studies [START_REF] Herbin | Conservativity and weak consistency of a class of staggered finite volume methods for the Euler equations[END_REF].