Huy Duong

Brigitte Jaumard

David Coudert

Minimum Disturbance Rerouting to Optimize Bandwidth Usage

Keywords: Network reconfiguration, rerouting, reoptimization, make-before-break

Dynamic traffic leads to bandwidth fragmentation, which drastically reduces network performance, resulting in increased blocking rate and reduced bandwidth usage. When rerouting traffic flows at Layer 3 of an optical network, network operators are interested in minimizing the disturbances in order to satisfy their Service Level Agreements. Therefore, they turn to the Make-Before-Break (MBB) paradigm.

In this paper, we revisit MBB rerouting with the objective of identifying the reroute sequence planning that minimizes the number of reroutes in order to minimize the resource usage. We propose a Dantzig-Wolfe decomposition mathematical model to solve this complex rerouting problem. We instigate how multiple or parallel rerouting reduces the overall minimum number of rerouting events (shortest makespan), and achieve the best resource usage. Numerical results bring interesting insights on that question and show a computational time reduction by about one order of magnitude over the state of the art.

Introduction

Network resource fragmentation is incurred by dynamic traffic and reduces network efficiency [START_REF] Chatterjee | Fragmentation problems and management approaches in elastic optical networks: A survey[END_REF]. As part of the resource reallocation by rerouting, many studies have looked at optical path (lightpath) rerouting, both in the context of routing and wavelength allocation (RWA) [START_REF] Jaumard | Wavelength defragmentation for seamless migration[END_REF] and routing and spectrum allocation (RSA) [START_REF] Soumplis | Dynamic connection establishment and network re-optimization in flexible optical networks[END_REF], or IP / VPN rerouting to overcome failures. However, very few authors have studied rerouting at the IP layer (or Layer 3) for the sole purpose of overcoming fragmentation, optimizing the network resources and increasing the throughput. The latter problem was studied in the 90s in the context of MPLS rerouting, see, e.g., [START_REF] Klopfenstein | Rerouting tunnels for MPLS network resource optimization[END_REF] or [START_REF] Józsa | On the solution of reroute sequence planning problem in MPLS networks[END_REF], but then no recent study has been devoted to it, despite a few unanswered questions, and a renewed interest in the context of 5G network slicing. Indeed, network slicing provides a key functionality in 5G networks, and offers flexibility in creating customized virtual networks and supporting different services on a common physical infrastructure [6].

To reroute a connection, the network operator first finds a new route for this connection. To prevent the traffic disruption that occurs when the connection is switched to a new route after the termination of the current route, network operators favor the use of the make-beforebreak (MBB) paradigm. With MBB, the connection is switched to the new route before the termination of the current route, thus avoiding disruptions.

Recent advances for MBB rerouting [START_REF] Jaumard | Efficient realtime scalable make-before-break network re-routing[END_REF]8] enable to significantly improve upon the state of the art on the size of solved instances. However, these proposals did not consider the parallel and multiple MBB rerouting, i.e., the ability to reroute multiple connections in parallel or to reroute a connection more than once, which could help to further reduce the use of network resources. This approach has only been considered in few studies and the state of the art is the model proposed in [START_REF] Klopfenstein | Rerouting tunnels for MPLS network resource optimization[END_REF].

In this paper, we propose an optimization model that takes into account both parallel and multiple MBB rerouting for a capacitated network. We then present decomposition methods based on column generation algorithms to solve this model. The numerical results show that our algorithm is remarkably faster and provides better solutions than [8].

The paper is organized as follows. Section 2 contains the literature review. After a concise problem statement, we propose a decomposition model in Section 3, and its solution process in Section 4. Numerical results are presented in Section 5. Conclusions are drawn in the last section

Literature Review

Several studies have been devoted to network reconfiguration with the minimum number of disruptions, following the strategy of migrating from a legacy ineffective provisioning to a given pre-computed optimized/optimal one. As a result, it usually prevents the existence of a strategy using only MBB due to the presence of dependency cycles, also called deadlocks [START_REF] Maggi | Reroute backward to better break deadlocks[END_REF]. In order to find a rerouting strategy, authors have then proposed to use the Break-Before-Make (BBM) paradigm sparingly to allow temporary interruption of connection requests, and so to break dependency cycles. For instance, Jose and Somani [START_REF] Jose | Connection rerouting/network reconfiguration[END_REF] propose heuristics for minimizing the total number of BBMs used in the rerouting strategy, and Cohen et al. [START_REF] Cohen | Tradeoffs in process strategy games with application in the WDM reconfiguration problem[END_REF] and Solano [START_REF] Solano | Analyzing two conflicting objectives of the WDM lightpath reconfiguration problem[END_REF] provide scalable exact algorithms to minimize the concurrent number of BBMs and investigate tradeoffs between these two conflicting objectives.

To further reduce the total or concurrent number of BBMs, Kadohata et al. [START_REF] Kadohata | Wavelength path reconfiguration design in transparent optical WDM networks[END_REF] propose to use spare wavelengths to reroute a connection request to a temporary route rather than using a BBM. For example, assume that the current connection k needs to be rerouted from path p to path p , but such a rerouting cannot be MBB due to resource dependence. Then one unavoidable BBM rerouting is performed. However, using an intermediate reroute, it may be possible to reroute k under the MBB paradigm. For instance, assume that there exists a path p such that the reroutings from p to p and from p to p satisfy the MBB condition. In other words, one BBM can be avoided at the expense of performing two MBBs. This idea is similar to multiple rerouting for capacity reoptimization.

The idea of the second direction is to compute the best provisioning that is reachable from the legacy provisioning by a sequence of connection reroutings with no disruption, i.e., under the so-called MBB paradigm. While many studies have investigated the first direction, this second direction has received very little attention [START_REF] Klopfenstein | Rerouting tunnels for MPLS network resource optimization[END_REF][START_REF] Jaumard | Efficient realtime scalable make-before-break network re-routing[END_REF].

Problem Statement

Notations

We consider a network represented by a directed multi-graph G = (V, L), where V is the set of nodes (indexed by v) and L is the set of fiber links (indexed by ). Different links may exist between two nodes in order to model different logical links, with, e.g., different types of traffic. We denote ω -(v) (resp. ω + (v)) the set of incoming (resp. outgoing) links of node v ∈ V . Let C denote the transport capacity of link .

Let K be the set of connection requests (indexed by k). Connection request k ∈ K is characterized by its source s k , its destination d k , and its bandwidth requirement b k . In what follows, we call rerouting operation the action of rerouting a connection request k ∈ K, and rerouting event the action of either performing a single rerouting operation, or a set of parallel rerouting operations. A reoptimization event is an ordered sequence of rerouting events, and so of rerouting operations. Let T be the set of rerouting events, indexed by t.

A life-line configuration (or LL config. for short) for request k, denoted by γ, is defined by an ordered sequence of rerouting operations of a connection k, which is characterized by:

• λ γ t is 1 if k is rerouted at rerouting event t, 0 otherwise. • α γ t is 1 if k uses link at the end of rerouting event t. • δ γ t = α γ t-1 -α γ t
is the difference in usage of link between the begin of rerouting event t and its end.

• x γ t is 1 if k does not use link after rerouting event t -1 and uses link after rerouting event t.

Denote by Γ k the set of all feasible life-line configurations for connection k ∈ K with at most R mul rerouting operations.

Multiple Parallel MBB Reoptimization Model

Assuming that all feasible life-line configurations are enumerated, the multiple parallel MBB reoptimization model (mul_par_ro) is defined formally by the following decomposition mathematical model.

We first define its parameters:

• R = limit on the total number of rerouting operations.

• R = limit on the number of parallel rerouting operations at a rerouting event.

• R mul = limit on the number of rerouting operations for connection k.

• a init k = 1 if link ∈ L is used in the initial routing of connection request k ∈ K, 0 otherwise. • C init = k∈K b k a init k = initial load of link ∈ L.
It also uses the following variables:

• y γ = 1 if life-line γ ∈ Γ is selected, 0 otherwise. • C t = bandwidth usage on link ∈ L after rerouting event t ∈ T . Minimize ∈L C |T | (1)
subject to:

k∈K γ∈Γ k λ γ t y γ ≤ R t ∈ T (2) γ∈Γ ( t∈T λ γ t )y γ ≤ R (3) γ∈Γ k y γ ≤ 1 k ∈ K (4) C t ≤ C ∈ L, t ∈ T (5)
C init = k∈K b k a init k ∈ L (6)
C t-1 + k∈K γ∈Γ k b k λ γ t x γ t y γ ≤ C ∈ L, t ∈ T (7)
C t-1 + k∈K γ∈Γ k b k λ γ t δ γ t y γ ≤ C t ∈ L, t ∈ T ( 8
)
C t ≥ 0 ∈ L, t ∈ T (9)
y γ ∈ {0, 1} k ∈ K, γ ∈ Γ k (10)
The objective ( 1) is to minimize the capacity usage at the end of the reoptimization event. Constraints (2) prevent the selection of more than R rerouting operations at each rerouting event. Constraints (3) restrict the total number of rerouting operations of the whole reoptimization event. In case one does not want to restrict the number of parallel rerouting operations per rerouting event, R is set to infinite and Constraints (2) are removed from the model. Note that, in case R mul = ∞, |T | ≤ |K| is an upper bound on the number of rerouting operations.

Constraints (4) ensure that a connection request has at most one life-line. Constraints (5) make sure that transport capacities are never exceeded after any rerouting event.

Constraints (6) specifies the initial bandwidth usage of each link. Constraints [START_REF] Jaumard | Efficient realtime scalable make-before-break network re-routing[END_REF] ensure that the bandwidth which is needed on link for the "make" part does not exceed its capacity at rerouting event t. Note that if the old and new routes of a connection go through a same link, reserved capacity on that link is not duplicated. Constraints (8) update the bandwidth usage on link after rerouting event t. Constraints ( 9)-( 10) define the domain of the variables.

Solution Process

Life-Line Pricing Algorithmcpp

We use a column generation algorithm to solve the linear relaxation of Model ( 1)- [START_REF] Jose | Connection rerouting/network reconfiguration[END_REF] (see, e.g., [START_REF]Linear Programming[END_REF] if not familiar with column generation techniques). It results in decomposing the original problem into a Restricted Master Problem (RMP), i.e., Model (1)-( 10) with a very restricted number of variables, and one or several pricing problems (PPs), which are solved alternately. The process continues until the optimality condition is satisfied, i.e., all the so-called reduced costs (or equivalently, the optimal values of the objective functions of the pricing problems) are positive. An ε-optimal solution is derived by solving exactly the ILP model associated with the last RMP, with ε defined as follows:

ε = (z ilp -z lp ) /z lp , (11)
where z lp and zilp denote the optimal LP value and the optimal ILP value of the last RMP, respectively.

We now describe the pricing problem model. As the parameters of a life-line configuration become variables in the pricing problem, we do not change their names (slight abuse of notations).

Let u

(2) t ≤ 0, u (3)
t ≤ 0, u

k ≤ 0, u

t ≤ 0, and u

t ≤ 0 be the dual values of Constraints (2),( 3), ( 4), [START_REF] Jaumard | Efficient realtime scalable make-before-break network re-routing[END_REF], and (8) respectively.

Variables:

• α t is 1 if the routing of connection k uses link ∈ L at the end of rerouting event t, 0 otherwise.

• λ t is 1 if t is a rerouting event, i.e. if a rerouting operation occurs, 0 otherwise.

• x t is 1 if link ∈ L is used at the end of rerouting event t but was not used at t -1, 0 otherwise.

• δ t is 1 if link ∈ L is used at t but not at t -1, -1 if it was used at t -1 and is no longer used, and 0 if its usage is unchanged.

• p t v ≥ 0 encodes the relative position of vertex v ∈ V in the path at t. It is used for the elimination of subtours.

Objective (reduced cost):

[pp k ] ≡ min c pp k = - t∈T (u (2) t + u (3) )λ t -u (4) k - ∈L t∈T u (7) t b k λ t x t - ∈L t∈T b k u (8) t λ t δ t (12)
Subject to:

∈ω -(v) α t - ∈ω + (v) α t =      -1 if v = s k 1 if v = d k 0 otherwise t ∈ T, v ∈ V ( 13
) ∈ω + (v) α t ≤ 1 t ∈ T, v ∈ V (14)
α 0 = a init k ∈ L (15) t∈T λ t ≤ R mul (16) λ t ≥ α ,t-1 -α ,t t ∈ T, ∈ L (17)
x t ≥ α ,t -α ,t-1 t ∈ T, ∈ L (18) x t ≤ λ t t ∈ T, ∈ L (19) ∈L x t ≥ λ t t ∈ T ( 20
)
δ t = α ,t -α ,t-1 t ∈ T, ∈ L (21)
p t v ≥ p t u + 1 + |V |(α ,t -1) t ∈ T, (u, v) ∈ L (22)
Domains of the variables:

α t ∈ {0, 1} t ∈ T ∪ {0}, ∈ L (23) λ t ∈ {0, 1} t ∈ T (24) x t ∈ {0, 1}, t ∈ T, ∈ L (25) δ t ∈ {-1, 0, 1} t ∈ T, ∈ L (26)
p t v ≥ 0 t ∈ T, v ∈ V. (27)
The objective function ( 12) is quadratic. However, the first quadratic term λ t x t can be equivalently rewritten x t thanks to Constraints (19), and the second quadratic term λ t δ t can be equivalently rewritten δ t thanks to the combination of Constraints (17) and Constraints (21).

The overall set of Constraints ( 13)-( 26) describes the life-line of request k. Constraints ( 13) are flow conservation constraints defining a path for each connection after each rerouting event, while avoiding loops along the path thanks to Constraints [START_REF]Linear Programming[END_REF]. Constraints (15) specify the links that are used in the initial routing of request k. Constraint (16) limits the number of rerouting operations on the life-line of connection k. Constraints ( 17)-(20) identify rerouting operations and consequently rerouting events. Let us first assume that λ t = 0. Then, Constraints (17) ensure that the links used at t-1 are still used at t, and Constraints (18)-(19) prevent from using new links at t. Hence, the routing at t -1 and t are the same. Assume now that a rerouting operation occurs at t, that is λ t = 1. Constraints (20) ensure that at least one variable x t is set to 1, thus enabling (i) to use at t a link that was not used at t -1 with Constraints (19) (ii) to stop using at t some links that were used at t -1 with Constraints (17). On the other hand, if either a link is no longer used at t, or a link is used at t but was not used at t -1, Constraints (18)-( 19) identify that a rerouting event occurs at t and set variable λ t to 1. If no new link is used at t, Constraints (20) force variable λ t to 0, Constraints (17) force to continue using at t the links that were used at t -1, and so the paths at t -1 and t are the same. Constraints (21) encode in δ t the changes in link usage. Constraints (22) prevent the selection of cycles. Finally, Constraints (23)-( 27) define the domains of the variables.

Solution Process of Non-Multiple (Single) Rerouting per Connection

When R mul = 1, each connection k ∈ K can be rerouted at most once, and so, at most one λ t can be one. As in Duong et al. [8], we can decompose the compact pricing problem into a set of pricing problems pp sr kt in which λ t = 1. We then get:

c pp sr k = min t∈T c pp sr kt = -u (4) k -max t∈T u (2) t + u (3) + ∈L b k u (7) t x t + ∈L b k u (8) t δ t (28)
Assuming λ t = 1, Constraints ( 17) and ( 19) become redundant for the selected t, and can therefore be eliminated. The simplified pricing problem pp sr kt with single rerouting per connection can be written as follows.

[pp sr kt ] ≡ c pp sr kt ≡ min c pp sr kt = -u (2) t -u (3) -u (4) k + ∈L b k u (8) t a init k -b k ∈L (u (7)
t x t + u (8) t α t ) (29)
subject to:

∈ω -(v) α t - ∈ω + (v) α t =      -1 if v = s k 1 if v = d k 0 otherwise v ∈ V (30) ∈ω + (v) α t ≤ 1 v ∈ V (31) x t ≥ α t -0 ∈ L k : α ,t-1 = 0 (32) x t ≥ α t -1 ∈ L k : α ,t-1 = 1 (33) ∈L x t ≥ 1 ( 34
)
p v ≥ p u + 1 + |V |(α ,t -1) (u, v) = ∈ L (35) α t ∈ {0, 1}, x t ∈ {0, 1} ∈ L ( 36
)
p v ≥ 0 v ∈ V. ( 37
)
Note that in this pricing problem, when t -1 = 0, we have α ,t-1 = a init k .

Proposition 1. If we are looking for only negative optimal objective value (reduced cost), pricing problem pp sr kt can be reduced to a shortest path problem with non negative weights.

Proof. Firstly, we show that variables x can be eliminated. Let

u (7)(8) t = u (8) t if a init k, = 1 u (7) t + u (8) t if a init k, = 0 (38)
and consider the following model.

cc pp sr kt ≡ min cc pp sr kt = -u (2) t -u (3) -u (4) k + ∈L b k u (8) t a init k, - ∈L u (7)(8) t b k α ,t (39)
subject to Constraints (30)-(37)

We will prove that the model using Objective (39) is equivalent to the original pricing model of pp sr kt using Objective (29). To do that, we first prove the two following claims. (i) The feasible regions of these models are equal. This claim is trivial because these problems have the same sets of constraints.

(ii) The optimal solutions of these models are equal. This claim will be proved if the two following statements are true:

(ii-a) c pp sr kt ≥ cc pp sr kt ; (ii-b) ∃(x , α ) ∈ {(30) -(37)} : cc pp sr kt = c pp sr kt (x , α ).
To prove statement (ii-a), we must show that it is true for any valid assignment of the variables, that is:

∀(x , α ) c pp sr kt (x , α ) ≥ cc pp sr kt (x , α ). (40) Since u (7)
t ≤ 0 and u

t ≤ 0, this statement holds if

∀ ∈ L -u (7)
t x t -u (8)
t α t ≥ -u (7)(8) t α t
is true, and so equivalently, if the following statement is true.

∀ ∈ L u (7) t x t + u (8) t α t ≤ u (7)(8) t α t .
Indeed,

• If ∈ L is such that a 0 k, = 0, we have by Constraints (32) that x t ≥ α t . Since u

t ≤ 0, we get u (7)

t x t + u (8) t α t ≤ u (7) t α t + u (8)
t α t , and Equation (38) sets u

(7)(8) t = u (7) t + u (8)
t . • If ∈ L is such that a 0 k, = 1, we have by Constraints (33) that x t ≥ α t -1. Since x t ∈ {0, 1} and u [START_REF] Jaumard | Efficient realtime scalable make-before-break network re-routing[END_REF] t ≤ 0, we therefore have u (7)

t x t ≤ 0. Since Equation (38) sets u (7)(8) t = u (8)
t , we can conclude that u (7)

t x t + u (8) t α t ≤ u (7)(8) t α t .
We now prove statement (ii-b). Let x and α be the optimal solution corresponding to cc pp sr kt . We show that if α = α , and x such that

x t = 0 if a 0 k, = 1 α t if a 0 k, = 0 , (41)
then cc pp sr kt = cc pp sr kt (α , x ) = c pp sr kt (α , x ). First, we need to show that α and x are a feasible solution of (30)-(37). Clearly, α satisfies (30)-(31) because α = α . Furthermore, Constraints (32)-( 33) are satisfied by α and x because they hold in all cases of ∈ L. Indeed, we have

• If ∈ L is such that a init k, = 0, then we have x t = α t = α t ≥ α t -0. • If ∈ L is such that a init k, = 1, then we have x t = 0 ≥ α t -1 (as α t ∈ {0, 1}).
Note that we only consider the original pricing when cc pp sr kt ≤ 0. In that case, there has to exist one ∈ L such that a init k,l = 0 and α lt > 0. Otherwise, α t = a init k, for all ∈ L (the only feasible path in this case), then

cc pp sr kt = -u (2) t -u (3) -u (4) k > 0 =⇒ c pp sr kt ≥ cc pp sr kt > 0,
this contradicts the assumption that c pp sr kt ≤ 0. Overall, it concludes that, in this case, Constraint (34) is satisfied by x .

When cc pp sr kt > 0, it implies c pp sr kt > 0, so the pricing problem cannot generate improving configuration for the restricted master problem.

We now show that the objective values are equal. Note that cc pp sr kt does not depend on x , and so is computed with α only. Now, we show that

∀ ∈ L u (7)(8) t α t = u (7)
t x t + u (8) t α t . (42)
This holds trivially by the combination of (38) and (41).

To this end, we observe that variables x can be removed from the simplified problem and it becomes a shortest path problem with non negative weights.

Numerical Results

Data Sets

We use the instances resulting from the simulation performed in [8] on a network with 32 nodes and 250 directed links, which corresponds approximately to a Ciena customer network. Connection requests had Poisson arrivals based on a traffic matrix and random durations drawn from a common exponential distribution. Each connection had a Weibull distributed bandwidth with a coefficient of variation of 0.3. A load factor parameter was used to globally vary the connection arrival rates: the corresponding equilibrium connection states represent a range of congestion levels from light to heavy.

All experiments are run with Scientific Linux release 7.9 (Nitrogen), Intel(R) Xeon(R) CPU E7-4890 v2 @ 2.80GHz, 15 cores per socket, 2 threads per core, and 1056 GB of RAM.

Regular Multiple Parallel Rerouting CG Performance

We now evaluate the performances of the cpp algorithm in data instances with 12 nodes.

In Table 1, we report experiments when at most one rerouting operation is allowed per request (R mul = 1) while the number of rerouting operations per rerouting event is unlimited. The average gap of the cpp algorithm over 10 reoptimization events is 1.3%. It shows that the proposed algorithm is close to the optimal solution in this setting. We also observe that the time required to solve the last ILP problem is negligible in the overall resolution time of the cpp algorithm, meaning that the difficult task of the proposed algorithm is the generation of life-lines. In Table 2, we report experiments when two rerouting operations are allowed per request (i.e., R mul = 2) and the number of rerouting operations per rerouting event is unlimited. Before solving the last ILP problem of this case, we use the solutions of the experiments with R mul = 1 as additional columns to generated columns. The rows in Table 2 where the number of routing operations is larger than the number of selected life-line configurations (# ILP LLs) indicate that there are connections that are rerouted twice in the solution.

In Table 2, we compute the bandwidth usage gain when we can reroute a connection at most twice. Gain is not significant (less than 1%) when we allow 2 reroutings per connection for small networks (about 10 nodes), as connections are often initially routed on their shortest paths or close to shortest paths. So, there is very limited room for improvement. Note that, for the fifth rerouting event, the gain is negative. It is because the gap of the fifth rerouting event of Table 2 (0.2%) is worse than the corresponding event of Table 1 (0.1%).

When the number of allowed rerouting operations per request is increased, it leads to more difficult problem and takes more time to converge, as demonstrated by the significant difference of the number of generated life-lines ("# LP LLs" column) and total cpu time between Tables 1 and2. literature, i.e., with multiple rerouting operations per connection and parallel rerouting operations per rerouting event, while minimizing the overall number of rerouting in order to minimize the number of disturbance.

In addition, this modeling leads to an efficient algorithm to solve the specific case of one rerouting operation per connection: we can perform up to 400 rerouting operations in about 5 minutes, while the model of [8] was only scalable for 150 rerouting operations requiring about one hour.

Figure 1 :

 1 Figure 1: Reduction (%) of capacity requirement for load factors 0.5 and 1.

Table 1 :

 1 cpp Algo. Performance, T = 5, R = 15, R = ∞, R mul = 1, |V | = 12, |L| = 52

	Scenarios |K| ε (%) # LP LLs ILP cpu (s)	zilp	# rerouting op. # ILP LLs Init BW	z lp	cpu (s)
	1	168	1.9	637	0.1	165,765	11	11	176,879 162,721.0	65.3
	2	184	0.5	338	0.0	180,652	15	15	197,224 179,751.0	32.6
	3	190	1.3	381	0.0	179538	15	15	196,782 177,264.1	37.2
	4	191	3.7	653	0.1	189,029	12	12	200,844 182,311.4	55.4
	5	180	0.1	321	0.0	175989	15	15	196,493 175.751.4	48.7
	6	181	0.9	386	0.0	172824	9	9	182,751 171.357.0	39.1
	7	177	0.6	509	0.1	171347	15	15	188,211 170.345.0	44.2
	8	184	0.6	449	0.0	178,750	12	12	191574	177.726.0	52.4
	9	180	2.5	649	0.1	184,146	13	13	202182	179,693.9	62.7
	10	167	0.5	382	0.0	168583	15	15	190894	167677.6	43.5
	Average 180.2 1.3	470.5	0.1	176662.3	13.2	13.2	192383.4 174459.8 48.1

Table 2 :

 2 cpp Algo. Performance, T = 5, R = 15, R = ∞, R mul = 2, |V | = 12, |L| = 52 Scenarios |K| ε (%) # LP LLs ILP cpu (s)zilp Gain # rerouting op. # ILP LLs Init BW

	z lp	cpu (s)

Acknowledgments B. Jaumard has been supported by a Concordia University Research Chair (Tier I) and by an NSERC (Natural Sciences and Engineering Research Council of Canada) grant. H. Duong was supported by a MITACS & Ciena Converge Fellowship. D. Coudert has been supported by the French National Research Agency (ANR), through the UCA JEDI Investments in the Future project with the reference number ANR-15-IDEX-0001, and the Inria associated-team project EfDyNet.

Non-multiple Rerouting Algorithm

In this section, we evaluate the performance of the cpp algorithm when it is specialized for the particular case where each connection is allowed to be rerouted at most once, i.e., R mul = 1. As aforementioned in Section 4.2, the life-line generation pricing problem can be solved by decomposing it into a set of simple non-negative weighted shortest path problems, which can be solved using Dijkstra's algorithm.

In the first part of Table 3, we report average performance of specialized cpp algorithm over 10 reoptimization events for each traffic intensity (load) with T = 20, R = 400, R mul = 1 and R = 20. The reported performance parameters are as in Table 1 andTable 2. Besides, in this table, "BW reduction" column reports reduced portions of bandwidth requirements before and after rerouting events using proposed cpp algorithm (i.e., zilp-"Init. BW" "Init. BW"

). In the second part of Table 3, the similar results in [8] are recalled for the case where the total number of rerouting operations is 150.

We observe that the CPU times of non-multiple rerouting algorithm is about one order of magnitude smaller than the algorithm of [8], as we use a much smaller number of rerouting events (T = 20 instead of T = 150 in [8]), while preserving a comparable accuracy to that of [8].


T Gallouët

R Herbin

J.-C Latché

LAX-WENDROFF CONSISTENCY OF FINITE VOLUME SCHEMES FOR SYSTEMS OF NON LINEAR CONSERVATION LAWS: EXTENSION TO STAGGERED SCHEMES

Keywords: 2010 AMS Subject Classification. Primary 65M08, 76N15 ; Secondary 65M12, 76N19 Finite-volume schemes, convection, consistency

Lax-Wendroff consistency of finite volume schemes for systems of non linear conservation laws: extension to staggered schemes

Introduction

The well-known Lax-Wendroff theorem [START_REF] Lax | Systems of conservation laws[END_REF] states that, on uniform 1D grids, if the approximate solutions of a flux-consistent and conservative cell-centred finite volume (FV) scheme for a system of conservation laws converge a.e. and boundedly as the mesh and time steps tend to zero, then the limit is a weak solution of the conservation law; we call this property "Lax-Wendroff consistency" (or LW-consistency for short); it is also stated in a different form [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]Section 12.10], with a BV bound assumption on the scheme. The Lax-Wendroff theorem is an "if-theorem", which fails to solve the convergence issue for FV schemes since compactness is lacking. Nevertheless, it introduces two crucial tools for the analysis of FV schemes, namely the conservativity and consistency of the numerical fluxes; note that this analysis cannot be handled by the famous Lax-Richtmyer theorem, even in the linear case and even if the exact solution is assumed to be regular, as soon as the mesh is non uniform, see e.g. [START_REF] Eymard | Finite volume schemes and Lax-Wendroff consistency[END_REF] for more on this subject. Moreover, the Lax-Wendroff theorem remains a useful tool to check whether a particular scheme gives a reasonable approximation when no estimates on the approximate solutions are available to yield some compactness, such as in the case of general hyperbolic systems. The Lax-Wendroff theorem was generalised to non uniform 1D or Cartesian meshes in [START_REF] Eymard | Finite volume methods[END_REF]Theorem 21.2]. In a recent work [START_REF] Ben-Artzi | Consistency and convergence of finite volume approximations to nonlinear hyperbolic balance laws[END_REF], the Lax-Wendroff theorem is extended to obtain some error estimates for higher order schemes on uniform 1D meshes. The case of general (and, in particular, unstructured) discretisations has been also been tackled over the past decades: [START_REF] Kroner | A Lax-Wendroff type theorem for upwind finite volume schemes in 2-D[END_REF], [6, Section 4.2.2] [2], [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF]. In [START_REF] Elling | A Lax-Wendroff type theorem for unstructured quasi-uniform grids[END_REF], a quasi-uniformity assumption is required on the mesh, but the flux is only required to be continuous, while in [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF], there is no uniformity assumption on the mesh but the flux is supposed to be locally Lipschitz continuous or at least locally "Lipschitz-diagonal" continuous, see Section 3 below and [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF]Remark 5.2]. In all the above cited works, the scheme is supposed to be colocated, in the sense that the discrete unknowns are associated to the cells of the mesh, so these results may not be used directly to cope with staggered approximations, for instance.

The aim of this paper is to address more general approximations, including those of co-located or staggered type; indeed, we prove the LW-consistency of a generic finite volume convection operator acting on discrete functions that are possibly not piecewise-constant over the cells of the mesh and over the time steps; the result is otained under sufficient conditions which, in usual cases, turn to essentially boil down to a new flux consistency requirement; this LW-consistency result is stated in Theorem 2.1 below. The flux consistency constraint, formulated by Assertion [START_REF] Lax | Systems of conservation laws[END_REF], demands a control on the difference between the discrete flux at a face (or edge) and the mean value over the adjacent cell of the continuous flux function applied to the discrete unknown function. Theorem 2.1 is valid for general polygonal or polyhedral meshes without any supplementary assumptions on the mesh; as a by product of this work, we thus also obtain a consistency result for colocated schemes (i.e. schemes using only piecewise-constant per cell unknowns) with possibly relaxed assumptions for the mesh compared to [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF]. However, let us note that the proof that the assumption [START_REF] Lax | Systems of conservation laws[END_REF] is satisfied is usually based on the control of the difference between the numerical solution and its space or time translates, see [5, Section 4] and that these latter results may require some regularity assumptions on the mesh, see also Remark 2.2. This paper is organized as follows. We state and prove the general consistency result in Section 2. We then apply it in Section 3 to the colocated case and then, in Section 4, to a staggered discretisation; precisely speaking, we show the consistency of a finite volume discretisation of a nonlinear convection operator for a scalar variable ρ of the form ∂ t β(ρ) + div(g(ρ)u), where β and g are regular functions and u is a velocity field, and where we use a cell-centred approximation for ρ and a face-centred approximation of u.

The general LW-consistency result

The aim is to prove the LW-consistency of finite volume approximations of nonlinear convective terms which appear in most models of fluid flow. The general context is the following. Given a numerical scheme which yields some approximate solutions to the system of conservative partial differential equations, we assume that these approximate solutions converge to some functions strongly in L 1 , and we wish to show that the limit is indeed a solution to the system, at least in a weak sense. In order to do so, the usual idea is to mutiply the numerical scheme by an interpolate of a smooth function, sum over the cells of the mesh and over the time steps and show that passing to the limit, we get a weak formulation of the system of partial differential equations. The theorem that we prove below is a mean to prove that one may indeed pass to the limit in the terms that involve nonlinear convection operators. Let us begin with an example. Consider the barotropic Euler equations, which read:

∂ t ρ + div(ρū) = 0, (1a)
∂ t (ρū) + div(ρū ⊗ ū) + ∇p = 0, ( 1b
)
where ρ is the density, ū the velocity and p the pressure, which, for barotropic flows, is a function of ρ only: p = p(ρ). Here and in the remainder of the paper, we use overlined letters when referring to the solution of the continuous problem, while non overlined letters will be used for discrete unknowns. This system of equations is supplemented by an initial condition and suitable boundary conditions if Ω is bounded. An entropy weak solution of the system satisfies the equations (1) and also satisfies (in a weak sense, which includes the initial condition) the following entropy condition:

∂ t Ē + div(( Ē + p)ū) ≤ 0, with Ē = 1 2 ρ|ū| 2 + H(ρ) and H(s) = s p(s) s 2 ds. (2)
The weak consistency of staggered finite volume schemes for this system of equations discretised on multidimensional Cartesian or unstructured meshes has been the object of several recent papers, see e.g. [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].

The system (1)-( 2) may be written as

C1 (ρ, ū) = 0, (3a) C2 (ρ, ū) + ∇p = 0, (3b) C3 ( Ē, ū) + div(pū) ≤ 0, (3c)
with C1 (ρ, ū) = ∂ t ρ + div(ρū), C2 (ρ, ū) = ∂ t (ρū) + div(ρū ⊗ ū), and C3 ( Ē, ū) = ∂ t Ē + div( Ē ū). In the above cited works, the system is discretised with an explicit or implicit in time scheme, and the convection operators C 1 and C 2 by a first or second order finite volume scheme. In fact, the system of the barotropic equations can be discretised by different schemes: explicit or implicit, colocated meshes or staggered meshes, using a Riemann solver or using an equation-by-equation procedure. In all cases, the consistency study will have to deal with each of the discrete non linear convection operator C i associated to Ci . The present work aims at simplifying the proofs of consistency by giving a general result for any nonlinear convection term, discretised on a colocated or staggered mesh, thereby extending our previous result of [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF] to staggered meshes. Theorem 2.1 below is an efficient tool to this purpose; it may be used for any of the terms in (3), and is specifically useful to tackle the terms featuring discrete variables with different space approximations, as the operators C i of these equations in case of staggered discretizations. We emphasize that both implicit or explicit schemes may be addressed, since the proof deals separately with the discrete time operator and the discrete space divergence operator.

Let us then turn to the general setting; we suppose that:

Ω ⊂ R d , d = 1, 2, 3, T ∈ (0, +∞), p ∈ N * , β ∈ C 0 (R p , R), f ∈ C 0 (R p , R d ). (4)
We consider the conservative convection operator C( Ū ) acting on a vector Ū ∈ R p of functions, real-valued, and defined (in the distributional sense), for Ū ∈ L ∞ (Ω × (0, T ), R p ), by:

C( Ū ) : Ω × (0, T ) → R, (x, t) → ∂ t (β( Ū ))(x, t) + div(f ( Ū ))(x, t). (5)
Note that, here and throughout the paper, we use β( Ū ) (resp. f ( Ū )) to denote the function β • Ū obtained by composition of β and Ū (resp. f and Ū ), so, for instance, β( Ū )(x, t) stands for β( Ū (x, t)). In the above example of the barotropic Euler equations (1), we have, for i = 1, 2, Ci ( Ū ) = ∂ t (β i ( Ū )) + div(f i ( Ū )), with Ū = (ρ, ū), β 1 ( Ū ) = ρ, f 1 ( Ū ) = ρū, β 2 ( Ū ) = ρū, and f 2 ( Ū ) = ρū ⊗ ū (in fact, to match precisely the formalism of Equation ( 4), these last two functions have to be considered as d functions, one for each velocity component, associated to d convection operators, which has no consequence for the matter at hand).

Let us denote by P a mesh of the domain Ω, consisting of a set of disjoint open polyhedral or polygonal subsets of Ω, called cells, whose union of closures is Ω. To avoid cumbersome notations, we assume that any pair of adjacent cells shares a whole face (in 3D) or edge (in 2D), and not only a part of it; however this assumption is not necessary for the result of Theorem 2.1 to hold. Throughout the paper and when the space dimension is not specified, we use "face" to define the interface between two cells; for a face ζ, |ζ| stands for its (d -1)-dimensional measure in 2D and 3D, and we set |ζ| = 1 by convention in one space dimension. The notation |P | stands for the d-dimensional measure of a cell P . We denote by δ(P) the space step, defined by δ(P) = max

P ∈P diam(P ).
Let F denote the set of faces of the mesh, and F int denote the set of faces that are not located on the boundary ∂Ω; for a given cell P ∈ P, let F(P ) be the set of faces of P . Let t 0 = 0 < t 1 < . . . < t N = T be a partition of (0, T ), denoted by T; for such a partition T, we define the time step by δt = max {t n+1 -t n , n ∈ 0, N -1 }, where 0, N -1 denotes the set of integers n such that 0 ≤ n ≤ N -1.

The unknown is supposed to be represented by a function U ∈ L ∞ (Ω × (0, T ), R p ). For a colocated FV scheme, it is the piecewise constant function defined by U (x, t) = U n K for x ∈ K and t ∈]t n , t n+1 [. For a FV staggered scheme for a system of equations, each component of U is piecewise constant on each associated mesh. But U could also be a piecewise affine function, for instance if say, a DG scheme is used. We emphasize that for non colocated schemes, some unknowns are not piecewise-constant over the cells of the mesh and over the time steps. For instance, when using staggered discretisations in fluid flow simulations, if P is a primal cell, the velocity is possibly discontinuous along surfaces or lines included in P (see the example developed in Section 4). The discrete convection operator that we consider here takes the following form:

P Q ζ = P | Q n P,ζ
C(U ) : Ω × (0, T ) → R, (x, t) → C(U ) n P ,
β n+1 P -β n P t n+1 -t n , n ∈ 0, N -1 , (6)
F n ζ , ζ ∈ F, n ∈ 0, N - 1
is a family of real vectors of R d and n P,ζ stands for the normal vector to ζ pointing outward P . Note that this form of the flux implies that the scheme is conservative. Of course, the real numbers β n P , P ∈ P, n ∈ 0, N and F n ζ , ζ ∈ F, n ∈ 0, N -1 are related to the unknown U ; it is the object of Theorem 2.1 below to state precisely the assumptions that must be satisfied by these quantities to ensure the consistency of the discrete convection operator.

Theorem 2.1 (LW-consistency for a multi-dimensional conservative convection operator). Under the assumptions (4), let (P (m) , T (m) ) m∈N be a sequence of possibly non uniform space-time discretisations, with δ(P (m) ) and δt (m) tending to zero as m → +∞, and let (U (m) ) m∈N be the associated sequence of discrete functions. We suppose that the sequence (U (m) ) m∈N is bounded and converges to a limit:

∃ C u ∈ R * + s.t. U (m) ∞ ≤ C u , ∀m ∈ N, (7)
∃ Ū ∈ L ∞ (Ω × (0, T ), R p ) s.t. U (m) -Ū L 1 (Ω×(0,T ),R p ) → 0 as m → +∞. ( 8
)
We also assume that the family {(β (m) ) n P , P ∈ P (m) , n ∈ 0, N (m) , m ∈ N} is bounded. In addition, let U 0 ∈ L ∞ (Ω, R p ) and let us suppose that, as m → +∞,

P ∈P (m) int P (β (m) ) 0 P -β(U 0 )(x) ϕ(x) dx → 0, for any ϕ ∈ C ∞ c (Ω), (9)
N (m) n=1 P ∈P (m) int tn tn-1 P (β (m) ) n P -β(U (m) )(x, t) ϕ(x, t) dx dt → 0, for any ϕ ∈ C ∞ c Ω × [0, T ) , (10)
N (m) -1 n=0 P ∈P (m) int diam(P ) |P | ζ∈F(P ) |ζ| tn+1 tn P (F (m) ) n ζ -f (U m )(x, t) • n P,ζ dx dt → 0, ( 11
)
where

P (m)
int denotes the set of cells of P (m) that have no face on the boundary ∂Ω. Then, for any

ϕ ∈ C ∞ c (Ω × [0, T )), T 0 Ω C (m) (U (m) ) I (m) (ϕ)(x, t) dx dt → - Ω β(U 0 )(x) ϕ(x, 0) dx - T 0 Ω β( Ū )(x, t) ∂ t ϕ(x, t) + f ( Ū )(x, t) • ∇ϕ(x, t) dx dt as m → +∞, ( 12
)
where I (m) (ϕ) is an interpolate of ϕ defined a.e. by

I (m) (ϕ)(x, t) = ϕ n P for x ∈ P and t ∈ (t n , t n+1 ), with ϕ n P = 1 |P | P ϕ(x, t n ) dx, for P ∈ P and n ∈ 0, N -1 . ( 13
)
Before we give the proof of Theorem 2.1, let us first briefly comment on its assumptions.

Remark 2.2 (Flux consistency). The required flux consistency is stated by Assertion [START_REF] Lax | Systems of conservation laws[END_REF], which requires the flux (F (m) ) n ζ through a face ζ of a cell P to be close to the mean value over P of the actual flux function f applied to the unknown. For a scheme involving only cell unknowns, for instance, the quantity (F (m) ) n ζ is generally a function of the unknowns in the cell P and in the neighbouring cells, and checking the assumption [START_REF] Lax | Systems of conservation laws[END_REF] amounts to bound the difference between the unknowns and their translates. Note that, while Theorem 2.1 holds for very general meshes, as we have already mentioned in the introduction, some regularity assumptions on the sequence of meshes may be required at this step. To clarify this point, let us consider a simple one-dimensional problem for the scalar unknown u, with β(u) = f (u) = u, leading to the linear convection operator C(u) = ∂ t u+∂ x u, which we discretise with the firstorder explicit-in-time upwind scheme. Let us suppose that the discrete functions are defined by u(x, t) = u n P for x ∈ P and t ∈ (t n , t n+1 ). Then, for x ∈ P and t

∈ (t n , t n+1 ), |((F (m) ) n ζ -f (U m )(x, t))•n P,ζ | = |(u (m) ) n P -- (u (m) ) n P |
where P -is the left cell to P when ζ is its left face, and |((

F (m) ) n ζ -f (U m )(x, t))
•n P,ζ | = 0 otherwise (disregarding the boundary cells, according to the formulation of the theorem). Checking Assumption [START_REF] Lax | Systems of conservation laws[END_REF] thus consists in proving that the term R (m) defined by

R (m) = N (m) -1 n=0 (t n+1 -t n ) P ∈P (m) diam(P ) |u n P -u n P -|
tends to zero as m tends to +∞. This is implied by the convergence in L 1 (Ω × (0, T )) of the sequence of discrete solutions provided that the ratio |P |/|P -| is bounded independently of m for the sequence of meshes under consideration [5, Section 4]. More elaborate example of application, using colocated then staggered meshes, are provided below, in Sections 3 and 4 respectively. ) that boundary cells may be excluded in the sum in Assertion [START_REF] Lax | Systems of conservation laws[END_REF]. This is the reason why only the cells in P (m) int are considered in Assumption [START_REF] Lax | Systems of conservation laws[END_REF]. For numerical fluxes involving wider stencils (for instance in the case of higher order schemes), one could in fact reduce the set of involved cells furthermore.

Remark 2.4 (Regularity of β and f ). The proof of Theorem 2.1 holds if β and f are only continuous functions, which is the assumption made in the present section; however, to prove Assertions ( 10) and ( 11), a locally Lipschitz-diagonal continuity (see Definition 3.1 below) is often required, as in sections 3 and 4.

Remark 2.5 (Stronger convergence assumptions on {(β (m) ) m∈N }). In most situations, stronger convergence properties hold for (β (m) ) m∈N , in the sense that the LW-convergence assumptions ( 9) and ( 10) are implied by the following strong convergence asssumptions:

P ∈P (m) int P |(β (m) ) 0 P -β(U 0 (x))| dx → 0 as m → +∞, Nm-1 n=0 P ∈P (m) int tn+1 tn P |(β (m) ) n P -β(U (m) (x, t))| dx dt → 0 as m → +∞.
This is the case, for instance, for the convection operators considered in Sections 3 and 4 below. However, there are cases where the convergence of β (m) is only weak, see for instance the reconstructed kinetic energy for the full compressible Euler equations in [START_REF] Herbin | Conservativity and weak consistency of a class of staggered finite volume methods for the Euler equations[END_REF].

Remark 2.6 (On the interpolate of the test function). Note that in the definition (13) of I (m) (ϕ) in [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF], the quantities ϕ n P , n ∈ 0, N , may be also defined as

ϕ n P = 1 |P | P ϕ(x, t n+1 ) dx,
with minor changes in the arguments of the present section, essentially a slightly different assumption [START_REF] Kroner | A Lax-Wendroff type theorem for upwind finite volume schemes in 2-D[END_REF], which reads:

N (m) n=1 P ∈P (m) int tn tn-1 P (β (m) ) n-1 P -β(U (m) )(x, t) ϕ(x, t) dx dt → 0, for any ϕ ∈ C ∞ c Ω × [0, T ) .
For instance, for a scalar problem, if the discrete function is defined as u(x, t) = u n-1 P for x ∈ P and t ∈ [t n-1 , t n ) (choice often used in explicit schemes) and β n-1 P is defined in the scheme as β(u n-1 P ), this assumption is trivially satisfied, since (β (m) ) n-1 P = β(U (m) )(x, t) in P ×(t n-1 , t n ), while checking the original assumption [START_REF] Kroner | A Lax-Wendroff type theorem for upwind finite volume schemes in 2-D[END_REF] needs to bound the time translates of the discrete solution. This is however an easy task, under a very mild regularity assumption for the time discretisation (see Section 4 below). The opposite situation occurs (i.e. this is Assumption [START_REF] Kroner | A Lax-Wendroff type theorem for upwind finite volume schemes in 2-D[END_REF] which is now trivially satisfied) if the discrete function is defined as u(x, t) = u n P for x ∈ P and t ∈ [t n-1 , t n ), which is often done for implicit schemes.

Proof of Theorem 2.1. Theorem 2.1 is the consequence of the two following lemmas, which prove respectively the convergence of the time derivative part and the space derivative part. Let us decompose

T 0 Ω C (m) (U (m) ) I (m) (ϕ)(x, t) dx dt = X (m) 1 + X (m) 2
, with

X (m) 1 = N (m) -1 n=0 (t n+1 -t n ) P ∈P (m) |P |(ð t β (m) ) n P ϕ n P , (14)
X (m) 2 = N (m) -1 n=0 (t n+1 -t n ) P ∈P (m) ζ∈F(P ) |ζ| (F (m) ) n ζ • n P,ζ ϕ n P . (15)
Then, by Lemma 2.7 below,

X (m) 1 → - Ω β(U 0 )(x) ϕ(x, 0) dx - T 0 Ω β( Ū )(x, t) ∂ t ϕ(x, t) dx dt as m → +∞,
and by Lemma 2.8 below,

X (m) 2 → - T 0 Ω f ( Ū )(x, t) • ∇ϕ(x, t) dx dt as m → +∞,
which concludes the proof.

Lemma 2.7 (LW-consistency, time derivative). Let the sequence (X

) m∈N be defined by (14). Then, under the assumptions and notations of Theorem 2.1,

X (m) 1 → - Ω β(U 0 )(x) ϕ(x, 0) dx - T 0 Ω β( Ū )(x, t) ∂ t ϕ(x, t) dx dt as m → +∞.
Proof. By the definition (6) of ð n t β (m) P (x, t) and thanks to a discrete integration by parts, we get that

X (m) 1 = - P ∈P (m) |P | (β (m) ) 0 P ϕ 0 P - N (m) n=1 (t n -t n-1 ) P ∈P (m) |P | (β (m) ) n P ϕ n P -ϕ n-1 P t n -t n-1 .
Let us write the first term of the right-hand side as

- P ∈P (m) |P | (β (m) ) 0 P ϕ 0 P = - P ∈P (m) P (β (m) ) 0 P ϕ 0 P -ϕ(x, 0) dx - P ∈P (m) P (β (m) ) 0 P ϕ(x, 0) dx.
On the one hand, the piecewise-constant function equal to ϕ 0 P on each cell P ∈ P (m) converges to ϕ(x, 0) in L ∞ (Ω) as m tends to +∞, and (β (m) ) 0 is supposed to be bounded; the first integral at the right-hand side thus tends to zero. Hence, invoking Assumption (9) for the second integral,

- P ∈P (m) |P | (β (m) ) 0 P ϕ 0 P → - Ω β(U 0 )(x) ϕ(x, 0) dx as m → +∞.
Let the piecewise constant function ð

(m) t ϕ : Ω × (0, T ) → R d be defined by ð (m) t ϕ(x, t) = ϕ n+1 P -ϕ n P t n+1 -t n for (x, t) ∈ P × (t n , t n+1 ).
The function ð

(m) t ϕ converges uniformly to ∂ t ϕ in L ∞ (Ω × (0, T )). The second term of X (m) 1
may be decomposed as

- N (m) n=1 (t n -t n-1 ) P ∈P (m) |P | (β (m) ) n P ϕ n P -ϕ n-1 P t n -t n-1 = Y (m) 1 + Y (m) 2 with Y (m) 1 = - N (m) n=1 P ∈P (m) tn tn-1 P (β (m) ) n P ð (m) t ϕ(x, t) -∂ t ϕ(x, t) dx dt, Y (m) 2 = - N (m) n=1 P ∈P (m) tn tn-1 P (β (m) ) n P ∂ t ϕ(x, t) dx dt.
Since the family {(β (m) ) n P , P ∈ P (m) , n ∈ 0, N (m) , m ∈ N} is assumed to be bounded, the uniform convergence of ð ) m∈N tends to zero. Invoking the assumption [START_REF] Kroner | A Lax-Wendroff type theorem for upwind finite volume schemes in 2-D[END_REF], the continuity of β and the convergence of (U (m) ) m∈N to Ū , we get that

lim m→+∞ X (m) 1 = lim m→+∞ Y (m) 2 = - T 0 Ω β( Ū )(x, t) ∂ t ϕ(x, t) dx dt.
Lemma 2.8 (LW-consistency, space derivative). Let the sequence (X

) m∈N be defined by (15). Then, under the assumptions and notations of Theorem 2.1,

X (m) 2 → - T 0 Ω f ( Ū )(x, t) • ∇ϕ(x, t) dx dt as m → +∞.
Proof. Since ϕ is compactly supported and since δ(P (m) ) → 0 as m → 0, there exists M ∈ N such that for m ≥ M , ϕ n P = 0 for all x ∈ P (m) \ P (m) int . Moreover, since for a face ζ separating P and P ′ , one has n P,ζ = -n P ′ ,ζ , we get that

X (m) 2 = N (m) -1 n=0 (t n -t n-1 ) P ∈P (m) int ζ∈F(P ) |ζ| (F (m) ) n ζ • n P,ζ ϕ n P = N (m) -1 n=0 (t n -t n-1 ) Pint∈P (m)
A n P with

A n P = ζ∈F(P ) |ζ| (F (m) ) n ζ • n P,ζ (ϕ n P -ϕ n ζ ),
where ϕ n ζ denotes the mean value of ϕ(x, t n ) over ζ. For any x ∈ P , t ∈ [t n , t n+1 ), we decompose A n P as

A n P = B n P (x, t) + R n P (x, t) with B n (x, t) = ζ∈F(P ) |ζ| f (U (m) )(x, t) • n P,ζ (ϕ n P -ϕ n ζ ), R n P (x, t) = ζ∈F(P ) |ζ| (F (m) ) n ζ -f (U (m) )(x, t) • n P,ζ (ϕ n P -ϕ n ζ ). (16)
Since ζ∈F(P )

|ζ| n P,ζ = 0, we have

B n P (x, t) = - ζ∈F(P ) |ζ| f (U (m) )(x, t) • n P,ζ ϕ n ζ = -|P | f (U (m) )(x, t) • (∇ϕ) n P , with (∇ϕ) n P = 1 |P | ζ∈F(P ) |ζ| ϕ n ζ n P,ζ = 1 |P | P ∇ϕ(x, t n ) dx. ( 17
)
Note that the piecewise constant function ∇ (m) ϕ : Ω × (0, T ) → R d defined by

∇ (m) ϕ(x, t) = (∇ϕ) n P for (x, t) ∈ P × (t n , t n+1 ) converges uniformly to ∇ϕ in L ∞ (Ω × (0, T )) d .
Owing to (16), we have

A n P = 1 (t n+1 -t n ) |P
=

N (m) -1 n=0 P ∈P (m) int 1 |P | tn+1 tn P B n P (x, t) dx dt + tn+1 tn P R n P (x, t) dx dt = - T 0 Ω f (U (m) )(x, t) • ∇ (m) ϕ(x, t) dx dt + N (m) -1 n=0 P ∈P (m) int 1 |P | tn+1 tn P R n P (x, t) dx dt. ( 18
)
Then, thanks to the boundedness and convergence assumptions on U (m) and to the uniform convergence of ∇ (m) ϕ to ∇ϕ, the first term tends to -

T 0 Ω f ( Ū )(x, t) • ∇ϕ(x, t) dx dt as m → +∞. Since |ϕ n ζ -ϕ n P | ≤ C ϕ diam(P )
, with C ϕ depending only on ϕ, we get, for any x ∈ P and t ∈ (t n , t n+1 ),

|R n P (x, t)| ≤ C ϕ ζ∈F(P ) |ζ| (F (m) ) n ζ -f (U (m) )(x, t) • n P,ζ diam(P ).
The second term of the right-hand side of Relation (18) thus tends to 0 as m → +∞ thanks to the assumption [START_REF] Lax | Systems of conservation laws[END_REF], which concludes the proof.

Application to a colocated scheme

In this section, we apply Theorem 2.1 to a first specific case, namely the case of a convection operator for a single scalar unknown ū, with a piecewise-constant discretisation of the unknown. We are going to prove a general consistency result for multipoint schemes, assuming minimal regularity of the mesh.

The considered convection operator reads:

C(ū) : Ω × (0, T ) → R, (x, t) → ∂ t (β(ū))(x, t) + div(f (ū))(x, t), (19)
where

Ω ⊂ R d , d = 1, 2, 3, T ∈ (0, +∞), β ∈ C 0 (R, R), f ∈ C 0 (R, R d ).
The functions β and f are supposed to be locally-Lipschitz. The discrete unknowns are (u n P ) P ∈P, n∈ 0,N -1 and the discrete convection operator reads:

C(u) : Ω × (0, T ) → R, (x, t) → β(u n+1 P ) -β(u n P ) t n+1 -t n + 1 |P | ζ∈F(P ) |ζ| F n ζ • n P,ζ , (20)
for x ∈ P, P ∈ P, and t ∈ (t n , t n+1 ), n ∈ 0, N -1 .

For a face ζ of the mesh, we denote by S ζ a set of cells, and suppose that the flux F n ζ reads

F n ζ = F ζ,n (u n L ) L∈S ζ .
The set S ζ is often referred to as the stencil of the scheme. We denote by (a) L∈S ζ the family of real numbers whose cardinal is the same as S ζ and whose elements are all equal to a. The flux is supposed to satisfy the usual consistency assumption:

for a ∈ R, ∀ζ ∈ F, F ζ,n (a) L∈S ζ = f (a). (21)
In addition, we suppose the following local "Lip-diag" (for Lipschitz-diagonal) property of the flux:

F ζ,n (u n L ) L∈S ζ -f (u n P ) ≤ C I L∈S ζ |u n L -u n P |. (22)
Note that this condition is weaker than the local Lipschitz-continuity of the numerical flux function F ζ,n . For P ∈ P, we denote by N 1 (P ) the set of the neighbours of P , i.e. the cells of P sharing a face with P ; for ℓ > 1, we define N ℓ (P ) as the set of the cells sharing a face with a cell of N ℓ-1 (P ). We assume that there exists l ∈ N such that ∀ζ ∈ F, for any cell P adjacent to ζ, S ζ ⊂ Nl(P ).

The integer l characterizes the compactness of the stencil of the scheme. The initial value for the scalar unknowm u is defined by

u 0 P = 1 |P | P u 0 (x) dx, ∀P ∈ P, (24)
where u 0 is a given function of L 1 (Ω). Finally, we define the discrete function associated to the unknowns as:

u(x, t) = u n P for x ∈ P, P ∈ P, t ∈ [t n , t n+1 ), n ∈ 0, N -1 .
The consistency result for the discrete convection operator is given in the next lemma; it uses the following regularity parameters of the mesh: We also measure the regularity of the time discretisation by the parameter θ 3 (T) defined by

θ 3 (T) = max 1≤n≤N -1 t n+1 -t n t n -t n-1 , t n -t n-1 t n+1 -t n .
Lemma 3.2 (Consistency, colocated scheme). Consider a sequence of space and time discretisations (P (m) ) m∈N and (T (m) ) m∈N , with δ(P (m) ) and δt (m) tending to zero; let (u (m) ) m∈N be an associated sequence of discrete functions, and let C (m) (u (m) ) be the associated sequence of discrete convection operators defined by (20). We assume that for each m ∈ N, (21)-( 24) hold with, in (23), l independent of m. We also suppose that

∃ θ ∈ R such that max{θ 1 (P (m) ), θ 2 (P (m) ), θ 3 (T (m) ), m ∈ N} ≤ θ, (25)
and that the number of faces of each cell of the meshes P (m) is bounded independently of m. Finally, we suppose that the sequence (u (m) ) m∈N is bounded in L ∞ (Ω × (0, T )), and that, when m tends to +∞, it converges in

L 1 (Ω × (0, T )) to ū ∈ L ∞ (Ω × (0, T )). Then, for any function ϕ ∈ C ∞ c (Ω × [0, T )), T 0 Ω C (m) (U (m) )(x, t) I (m) (ϕ) dx dt → - Ω β(u 0 )(x) ϕ(x, 0) dx - T 0 Ω β(ū)(x, t) ∂ t ϕ(x, t) + f (ū) (x, t) • ∇ϕ(x, t) dx dt as m → +∞.
Proof. Let m u ∈ R and M u ∈ R be two real numbers such that

m u ≤ (u (m) ) n P ≤ M u , ∀P ∈ P (m) , n ∈ 0, N (m) , ∀m ∈ N,
and let C β be the Lipschitz modulus of β over the interval [m u , M u ]. We check the assumptions of Theorem 2.1. The consistency of the initialization with the initial condition (Assumption ( 9)) follows from its definition (24); indeed, for any ϕ ∈ C ∞ c (Ω),

P ∈P (m) P (β (m) ) 0 P -β(u 0 )(x) ϕ(x) dx ≤ C β ϕ L ∞ (Ω) P ∈P (m) P |u 0 (x) -u 0 P | dx,
and thus tends to zero for any function u 0 ∈ L 1 (Ω). Since (β (m) ) n P = β((u (m) ) n P ), the left-hand side of Assertion [START_REF] Kroner | A Lax-Wendroff type theorem for upwind finite volume schemes in 2-D[END_REF]

β (u (m) ) n P -β (u (m) ) n-1 P ϕ(x, t) dx dt.
We thus have

|R (m) t | ≤ C β ϕ L ∞ (Ω×[0,T )) N (m) n=1 (t n -t n-1 ) P ∈P (m) int |P | |(u (m) ) n P -(u (m) ) n-1 P |,
and thus R (m) t tends to zero thanks to the assumed regularity of the sequence of time discretisations, invoking the bound of the time-translates of a converging sequence of functions of L 1 (Ω × (0, T )) stated by Lemma A.1 in Appendix. We now check Assumption [START_REF] Lax | Systems of conservation laws[END_REF]. For n ∈ 0, N (m) , P ∈ P (m) int and ζ ∈ F(P ), let

R n P,ζ = 1 |P | tn+1 tn P (F (m) ) n ζ -f (u (m) )(x, t) • n P,ζ dx dt and let R (m) = N (m) -1 n=0 P ∈P (m) int diam(P ) ζ∈F |ζ| R n P,ζ .
By definition of the discrete flux and the discrete functions, we get:

R n P,ζ = 1 |P | tn+1 tn P F ζ,n (u (m) ) n L L∈S ζ -f (u (m) ) n P • n P,ζ dx dt.
Thanks to Assumption (22), we thus have

R n P,ζ ≤ (t n+1 -t n ) C [mu,Mu] L∈S ζ |u n L -u n P |.
The remainder term R (m) is thus a collection of differences between the values taken by the unknown in two different cells. In order to apply Lemma A.3 of the appendix, we need to evaluate, for P ∈ P (m) , the sum ω P of the weights multiplying the jumps where (u (m) ) n P appears. We first notice that, thanks to the assumed regularity of the mesh, for

Q ∈ P (m) and ζ a face of Q, diam(Q) |ζ| ≤ diam(Q) d ≤ θ |Q|.
For P to appear in a difference associated to a face of a cell Q, we need, by assumption, that P ∈ Nl(Q); this in turn requires that Q ∈ Nl(P ). The sum ω P thus satisfies:

ω P ≤ θ L∈Nl(P ) |L|,
and, invoking once again the regularity constraints on the mesh:

ω P ≤ θ l+1 card(Nl(P )) |P |.
The proof is thus complete thanks to Lemma A.3, since we have supposed that the number of faces of the cells is uniformly bounded, and then so is card(Nl(P )).

An example of application for staggered discretisations

The interest of Theorem 2.1 lies in the fact that it may deal with terms combining several variables, associated to different meshes and time discretisations. A typical exemple of a such a case is the balance equation for the entropy in barotropic compressible flows [START_REF] Elling | A Lax-Wendroff type theorem for unstructured quasi-uniform grids[END_REF], where the entropy E is a nonlinear function of the density ρ and the velocity u which, in staggered discretisation, are approximated on different meshes, and may also be evaluated at different time levels. Hence, Theorem 2.1 is a suitable tool to prove the consistency of this equation. In this section, we focus on a similar but simpler problem, namely a staggered discretisation of a convection operator combining the time derivative of the function of a single scalar variable and a space divergence term, with a flux obtained as the product of another function of this scalar variable with the velocity.

We suppose that Ω is an open bounded polygonal set of R 2 , and consider the following convection operator:

C( Ū ) : Ω × (0, T ) → R, (x, t) → ∂ t (β(q))(x, t) + div g(q) v (x, t), ( 26
) with Ū = (q, v) : Ω × (0, T ) → R × R 2 , f ( Ū ) = f (q, v) = g(q)
v, where β : R → R and g : R → R are locally Lipschitz-continuous real functions. Note that, for instance, the convection term of Equation (1a) may be written as (26) with Ū = (ρ, ū), β(s) = s and g(s) = s.

In order to discretise this convection operator, we consider two types of staggered arrangements. In both arrangements, the scalar unknowns are located at the centre of the cells. However, these arrangements differ in the use of the vector unknowns. The first discretisation uses the whole velocity vector unknown on each edge of the mesh; this corresponds to the Rannacher-Turek (RT) discrete unknowns in the finite element setting [START_REF] Rannacher | Simple nonconforming quadrilateral Stokes element[END_REF]. The second discretisation uses only the normal component of the velocity on each edge; this latter arrangement of the discrete unknowns is very often referred to as the Marker-and-Cell (MAC) scheme [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF]. Hence we will refer to the first arrangement as the RT case, and the second as the MAC case. Such discretisations are called staggered and are widely used in computational fluid dynamics; an example of the implementation of a staggered discretisation for the solution of the barotropic and full Euler equations may be found e.g. in [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]. We suppose that the mesh is composed either of general quadrangles (RT case), or of rectangles (MAC case). We recall that F stands for the set of edges of the mesh, and the internal edge separating the cells P and Q is denoted by ζ = P |Q. This mesh will be referred to in the following as the primal mesh.

We also introduce now one or two dual meshes, depending on the case.

-RT case -In this case, the (unique) dual mesh consists in a new partition of Ω indexed by the elements of F, i. and that the sub-cells connectivities (i.e. the way the half-dual cells share a common edge) is left unchanged. Note that the actual geometry of the dual cells does not need to be specified (and a dual cell may not be a polytope, a dual edge being possibly curved). -MAC case -In this case, two dual meshes are considered, each of them consisting in a partition of Ω indexed by the vertical and horizontal elements of The scalar unknown q is associated to the primal cells:

P Q P Q ζ = P | Q D Q ,ζ D P ,ζ
F, i.e. Ω = ∪ ζ∈F (i) D ζ , i = 1, 2, where F (1) (resp. F (
q(x, t) = q n P for x ∈ P, P ∈ P, t ∈ [t n , t n+1 ), n ∈ 0, N -1 .
The unknowns corresponding to the vector-valued unknown v are located at the centre of the edges in the RT case; in the MAC case, the unknowns associated to the i-th component of v are located at the centre of the edges of the i-th dual mesh. Hence the associated approximate vector function reads:

M N P Q ζ = M|Q ζ ′ = N|P κ = M|N κ ′ = P|Q D ζ D ′ ζ M N P Q ζ ′ ζ D κ D κ ′ M N P Q κ κ ′ Figure 3.
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) .

-RT case -the whole vector unknown is associated to each dual cell :

v(x, t) = v n ζ for x ∈ D ζ , ζ ∈ F, t ∈ [t n , t n+1 ), n ∈ 0, N -1 .
-MAC case -The i-th component of the vector unknown is associated to the cells of the i-th dual mesh, so that v(x, t) = (v 1 (x, t), v 2 (x, t)) t where, for i = 1, 2,

v i (x, t) = v 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 of the canonical basis of R 2 ; with the notations of the previous section, the considered discrete convection operator reads:

C P (q, v) n P = (ð t β) n P + 1 |P | ζ∈F(P ) |ζ| F n ζ • n P,ζ
, with β n P = β(q n P ) and

F n ζ = f (q n ζ , v n ζ ) = g(q n ζ ) v n ζ
where v n ζ is the vector of discrete unknowns in the RT case, defined as v n ζ e (i) for ζ ∈ F (i) , i = 1 or 2, in the MAC case, and, for ζ = P |Q, q n ζ stands for a convex combination of q n P and q n Q . For instance the upwind choice would be q n ζ = q n P if v n ζ ≥ 0 and q n ζ = q n Q otherwise. Note that for the LW-consistency result, any convex combination works, but this is not so for the stability of the scheme, for which some unpwinding is required.

The initial value for the scalar unknown q is defined by

q 0 P = 1 |P | P q 0 (x) dx. ( 27
)
The consistency result for the discrete convection operator is given in the next lemma; it uses the same regularity parameters of the mesh as in the colocated case, which we recall:

θ 1 (P) = max P ∈P diam(P ) 2 |P | , θ 2 (P) = max |P | |Q| , P and Q adjacent cells of P .
Note that in the MAC case (in fact, for a Cartesian grid), the regularity parameter θ 1 (P) controls the ratio between the two dimensions (i.e. the height and the width) of a cell. For a rectangular computational domain, we thus observe that the ratio

|ζ|/|ζ ′ |, for (ζ, ζ ′ ) ∈ (F (i) ) 2 , i = 1, 2
, is bounded by θ 1 (P) 2 , which is a quasi-uniformity property of the mesh. This also implies that θ 2 (P) ≤ θ 1 (P) 2 , and so the second regularity parameter is useless. It may easily be checked that similar relations holds for a general MAC scheme, i.e. a union of matching Cartesian grids, with powers of θ 1 (P) possibly higher than 2. Hence, the regularity of a MAC mesh (or of a Cartesian grid) may be equivalently measured by

θ(P) = max h(1) h (2) , h (2)
h (1) ,

with, for i = 1, 2, h(i) = max{|ζ|, ζ ∈ F (i) } and h (i) = min{|ζ|, ζ ∈ F (i) }.
We also measure the regularity of the time discretisation by the parameter θ 3 (T) defined by

θ 3 (T) = max 1≤n≤N -1 t n+1 -t n t n -t n-1 , t n -t n-1 t n+1 -t n .
Lemma 4.1 (Consistency, staggered scheme). Let a sequence of discretisations (P (m) ) m∈N and (T (m) ) m∈N be given, with δ(P (m) ) and δt (m) tending to zero, and let (q (m) , v (m) ) m∈N be the associated sequence of discrete functions. We suppose that

∃ θ ∈ R such that max{θ 1 (P (m) ), θ 2 (P (m) ), θ 3 (T (m) ), m ∈ N} ≤ θ. ( 28
)
We suppose that the sequences (q (m) ) m∈N and (v (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 q ∈ L ∞ (Ω × (0, T )) and v ∈ L ∞ (Ω × (0, T )) 2 respectively. Then, for any function

ϕ ∈ C ∞ c (Ω × [0, T )), T 0 Ω C (m) (U (m) )(x, t) I (m) (ϕ) dx dt → - Ω β(q 0 )(x) ϕ(x, 0) dx - T 0 Ω β(q)(x, t) ∂ t ϕ(x, t) + g(q) v) (x, t) • ∇ϕ(x, t) dx dt as m → +∞.
Proof. In this proof, we denote by C β and C g the Lipschitz modulus of β and g respectively on the interval [m q , M q ], where m q ∈ R and M q ∈ R are such that m q ≤ (q (m) ) n P ≤ M q , ∀P ∈ P (m) , n ∈ 0, N (m) , ∀m ∈ N.

The proof of this lemma relies on Theorem 2.1. The consistency of the initialization with the initial condition (Assumption ( 9)) follows from its definition (27); indeed, for any ϕ ∈ C ∞ c (Ω),

P ∈P (m) P (β (m) ) 0 P -β(q 0 )(x) ϕ(x) dx ≤ C β ϕ L ∞ (Ω) P ∈P (m) P |q 0 (x) -q 0 P |,
and thus tends to zero for any function q 0 ∈ L 1 (Ω). Since (β (m) ) n P = β((q (m) ) n P ), the left-hand side of Assertion [START_REF] Kroner | A Lax-Wendroff type theorem for upwind finite volume schemes in 2-D[END_REF] β (q (m) ) n P -β (q (m) ) n-1 P ϕ(x, t) dx dt.

We thus have

|R (m) t | ≤ C β ϕ L ∞ (Ω×[0,T )) N (m) n=1 (t n -t n-1 ) P ∈P (m) int |(q (m) ) n P -(q (m) ) n-1 P |,
and thus R (m) t tends to zero thanks to the assumed regularity of the sequence of time discretisations, invoking the bound of the time-translates of a converging sequence of functions of L 1 (Ω × (0, T )) stated by Lemma A.1 in Appendix.

We now check Assumption [START_REF] Lax | Systems of conservation laws[END_REF]. For n ∈ 0, N (m) , P ∈ P (m) int and ζ ∈ F(P ), let

R n P,ζ = 1 |P | tn+1 tn P (F (m) ) n ζ -f (q (m) , v (m) )(x, t) • n P,ζ dx dt and let R (m) = N (m) -1 n=0 P ∈P (m) int diam(P ) ζ∈F |ζ| R n P,ζ .
We now express R n P,ζ , for the RT and MAC discretisations successively. -RT case -In the case of general quadrangular meshes with the whole vector unknowns located on the edges, we have

(F (m) ) n ζ = g(q n ζ ) v n ζ and f (q (m) , v (m) )(x, t) = g(q n P ) v n ζ ′ for x ∈ D P,ζ ′ , ζ ′ ∈ F(P ).
We thus get, denoting by |a| the Euclidean norm of any vector

a ∈ R 2 , (F (m) ) n ζ -f (U (m) )(x, t) • n P,ζ ≤ g(q n ζ ) v n ζ -g(q n P ) v n ζ ′ for x ∈ D P,ζ ′ , ζ ′ ∈ F(P ).
Let Q be the primal cell such that ζ = P |Q. Since q n ζ is a convex combination of q n P and q n Q , we thus get, for x ∈ P , and t ∈ [t n t n+1 ),

(F (m) ) n ζ -f (U (m) )(x, t) • n P,ζ ≤ C |q n P -q n Q | + ζ ′ ∈F(P ) |v n ζ -v n ζ ′ | ,
where C only depends on q (m) L ∞ (Ω×(0,T )) , v (m) L ∞ (Ω×(0,T )) 2 and C g . Integrating over P × (t n , t n+1 ), we obtain

R n P,ζ ≤ C (t n+1 -t n ) |q n P -q n Q | + ζ ′ ∈F(P ) |v n ζ -v n ζ ′ | .
-MAC case -In this case, the velocity components are piecewise constant on different grids. Let i be the index such that ζ ∈ F (i) , and let ζ ′ be the other edge of P normal to e (i) , i.e. the opposite of ζ in P . We have

(F (m) ) n ζ • n P,ζ = g(q n ζ ) v n ζ δ ζ and f (q (m) , v (m) )(x, t) = g(q n P ) v n ζ δ ζ if x ∈ D P,ζ , g(q n P ) v n ζ ′ δ ζ if x ∈ D P,ζ ′ , with δ ζ = n P,ζ • e (i) . We thus get (F (m) ) n ζ -f (U (m) )(x, t) • n P,ζ = g(q n ζ ) v n ζ -g(q n P ) v n ζ if x ∈ D P,ζ , g(q n ζ ) v n ζ -g(q n P ) v n ζ ′ if x ∈ D P,ζ ′ ,
and hence, for x ∈ P , and t ∈ [t n t n+1 ), denoting by Q the primal cell such that ζ = P |Q, where C only depends on q (m) L ∞ (Ω×(0,T )) , v (m) L ∞ (Ω×(0,T )) 2 and C g . Therefore, integrating over P × (t n , t n+1 ), we finally get

(F (m) ) n ζ -f (U (m) )(x, t) • n P,ζ ≤ C |q n P -q n Q | + |v n ζ -v n ζ ′ | , ζ ζ ′ η = ζ | ζ ′ ζ ζ ′′ ζ ′ η = ζ | ζ ′ ′ η ′ = ζ ′ | ζ ′ ′
R n P,ζ ≤ C (t n+1 -t n ) |q n P -q n Q | + |v n ζ -v n ζ ′ | .
Note that, in these computations, we have not addressed the case where ζ is an external edge, taking benefit of the fact that, in the expression of R (m) , the sum is retricted to the internal cells.

From the definition of R (m) , we thus get that, for both cases, it satisfies the following inequality:

R (m) ≤ C R (m) 1 + R (m) 2 , with R (m) 1 = N (m) -1 n=0 (t n+1 -t n ) P ∈P (m)
diam(P )

ζ∈F(P ), ζ=P |Q |ζ| |q n P -q n Q |, and
R (m) 2 =                  N (m) -1 n=0 (t n+1 -t n ) P ∈P (m) diam(P ) {ζ,ζ ′ }⊂F(P ) 2 (|ζ| + |ζ ′ |) |v n ζ -v n ζ ′ | in the RT case, N (m) -1 n=0 (t n+1 -t n ) P ∈P (m) diam(P ) i=1,2, {ζ,ζ ′ }⊂F (i) (P ) 2 (|ζ| + |ζ ′ |) |v n ζ -v n ζ ′ | in the MAC case.
There only remains to prove that R

R (m) 1 = N (m) -1 n=0 (t n+1 -t n ) ζ∈F (m) int , ζ=P |Q ω ζ |q n P -q n Q |, with ω ζ = diam(P ) + diam(Q) |ζ|.
Lemma A.1 states that R 2 , this is easily obtained using Assumption (28).

As to the term R (m) 2 , let us start with the RT case. We have:

{ζ,ζ ′ }⊂F(P ) 2 (|ζ| + |ζ ′ |) |v n ζ -v n ζ ′ | ≤ 2 diam(P ) {ζ,ζ ′ }⊂F(P ) 2 |v n ζ -v n ζ ′ |,
We distinguish two cases for the subsets {ζ, ζ ′ }) ⊂ F(P ) 2 that appear in the summation: either the dual cells D ζ and D ′ ζ share a common (dual) edge η = ζ|ζ ′ ∈ F * , where F * denotes the set of edges of the dual mesh, or they are opposite edges in the quadrilateral cell P ; in this latter case, we may write that

|v n ζ -v n ζ ′ | ≤ |v n ζ -v n ζ ′′ | + |v n ζ ′′ -v n ζ ′ |,
where ζ ′′ ∈ F(P ) is such that the dual cell D ζ ′′ shares a common (dual) edge η (resp. η ′ ) ∈ F * with D ζ (resp. D ζ ′ ) as shown in Figure 4. There is one jump between two adjacent faces that appears directly in the summation over {ζ, ζ ′ }) ⊂ F(P ) 2 , and at most two coming from the decompositions of the jumps needed for pairs of opposite edges, so that altogether,

(ζ,ζ ′ )∈F(P ) 2 (|ζ| + |ζ ′ |) |v n ζ -v n ζ ′ | ≤ 6 diam(P ) η=ζ|ζ ′ ∈F * (P ) |v n ζ -v n ζ ′ |,
with F * (P ) the edges of the dual mesh included in P . We thus get

R (m) 2 ≤ 6 N (m) -1 n=0 (t n+1 -t n ) η=ζ|ζ ′ ∈F * diam(P η ) 2 |v n ζ -v n ζ ′ |,
where P η stands for the primal cell in which η is included. The right-hand side of this inequality is thus a collection of jumps across the dual edges, with, for an edge η, a weight given by Let us now turn to the MAC case, which is in fact simpler; indeed, the differences of velocities appearing in the expression of R (m) 2 are all jumps across dual edges, and we may thus recast R

ω η = 6 diam(P η ) 2 . Thanks to Lemma A.1, R (m)
(m) 2 as R (m) 2 = N (m) -1 n=0 (t n+1 -t n ) 2 i=1 η=ζ|ζ ′ ∈F (i, * ) diam(P η ) (|ζ| + |ζ ′ |) |v n ζ -v n ζ ′ |,
where, once again, P η is the primal cell in which lies η (note that this sum only involves a subset of the dual edges, which corresponds of the dual edges included in primal cell), and F (1, * ) (resp. F (2, * ) ) denotes the set of vertical (respectively horizontal) dual edges. We thus again have a collection of jumps across the dual edges, with, for an edge η included in a primal cell P η and separating the dual cells D The assumption θ 3 (T (m) ) ≤ θ, for m ∈ N, may be avoided thanks to a different choice of the interpolation of the test function (see Remark 2.6). However, this assumption is very mild (in fact, we do not have in mind any scheme where the ratio between two consecutive time-steps is likely to blow up when refining the discretisation).

Let M be a conforming mesh of the domain Ω of R d , d = 1, 2, 3, in polygonal or polyhedral subsets, and T = (t i ) i∈ 0,N be a time discretisation of the interval (0, T ), i.e. a sequence of real numbers such that 0 = t 0 < • • • < t n < . . . t N = T . We denote by δt T the time step, defined by δt T = max{t n+1 -t n , n ∈ 0, N -1 }. For u ∈ L 1 (Ω × (0, T )), K ∈ M and n such that n ∈ 0, N -1 , let u n K be the mean value of u over K × (t n , t n+1 ). We denote by E int the set of internal faces of the mesh and the face σ ∈ E int separating the cells K and L is denoted by σ = K|L. We define the following quantity:

T M,T u = N -1 n=0 (t n+1 -t n ) σ∈Eint σ=K|L ω σ |u n K -u n L | + N -2 n=0 δ n+1/2 K∈M |K| |u n+1 K -u n K |, (29)
where (ω σ ) σ∈Eint and (δ n+1/2 ) n∈ 0,N -2 are two sets of non-negative weights. We introduce the two following parameters:

θ M = max K∈M max σ∈Eint(K) ω σ |K| , θ T = max n∈ 0,N -2 δ n+1/2 t n+1 -t n , δ n+1/2 t n+2 -t n+1 , (30)
with E int (K) the set of internal faces of K. We denote by δ(M) the space step characterizing M, i.e. δ(M) = max K∈M diam(K). Then the following convergence result holds.

Lemma A.1. Let θ > 0 and (M (m) ) m∈N be a sequence of meshes and for each m ∈ N, θ M (m) be defined by (30). We assume that θ M (m) ≤ θ for all m ∈ N and lim m→+∞ δ(M (m) ) = 0. We suppose that the number of faces of a cell K ∈ M (m) is bounded independently from m ∈ N. For m ∈ N, we suppose given a time discretisation T (m) , and suppose that δt T (m) also tends to zero when m tends to +∞, and that θ T (m) ≤ θ for all m ∈ N. Let u ∈ L 1 (Ω × (0, T )) and (u p ) p∈N be a sequence of functions of L 1 (Ω × (0, T )) such that u p → u in L 1 (Ω × (0, T )) as p → +∞.

Then T M (m) ,T (m) u p defined by (29) tends to zero when m tends to +∞ uniformly with respect to p ∈ N.

Remark A.2. The difference between Lemma A.1 and the formulation of the same convergence result in [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF] lies in the definition of the weight of the jumps, which is more general in Lemma A.1. Indeed, the weight of the jumps through the faces featured in the definition of T M,T u are defined in [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF]Lemma 4.3] as a function of the volume of some dual cells associated to the faces, but a careful examination of the proof itself shows that the introduction of a dual mesh is in fact useless. Therefore, the proof of Lemma [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF]Lemma 4.3] readily extends to prove Lemma A.1.

This generalization is in most cases sufficient. However, we may go one step further, still with minor modifications of the proof of [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF], as follows. Let S x be a set of cardinal 2 -subsets of M, and S t be a set of cardinal 2 -subsets of 0, N -1 . Let T M,T u be defined by

T M,T u = N -1 n=0 (t n+1 -t n ) {K,L}∈Sx ω K,L |u n L -u n K | + {p,q}∈St δ p,q K∈M |K| |u p K -u q K |, (31)
where (ω K,L ) {K,L}∈Sx and (δ p,q ) {p,q}∈St are two sets of non-negative weights. We introduce the two following parameters:

θ M = max

Figure 1 .

 1 Figure 1. An example of a two-dimensional mesh and associated notations: P and Q are two generic cells, P, Q ∈ P with P the set of cells, ζ = P |Q is the face separating P and Q, ζ ∈ F, with F the set of faces and n P,ζ is the normal to ζ pointing outward P .

Remark 2 . 3 (

 23 Disregarding boundary cells in Assumption[START_REF] Lax | Systems of conservation laws[END_REF]). Since the support of the test function ϕ is compact in Ω × [0, T ), for δ(P (m) ) small enough, ϕ vanishes in the boundary cells. Consequently, it is clear from the proof of the theorem below (see the expression (15) of the term X (m) 2

  (m) t ϕ to ∂ t ϕ yields that the sequence (Y (m) 1

Definition 3 . 1 (

 31 Local Lipschitz-diagonal continuity). The numerical flux function F ζ,n is said to be locally Lipschitz-diagonal continuous, or Lip-diag, if for any bounded interval I ⊂ R, there exists C I ∈ R + depending only on I such that, for any n ∈ 0, N -1 , for any face ζ ∈ F, for any family (u L ) L∈S ζ ⊂ I, and for any P adjacent cell to ζ,

θ 1 (

 1 Q adjacent cells of P .

N

  reads, with ϕ ∈ C ∞ c (Ω × [0, T )): m) ) n P -β(u (m) )(x, t) ϕ(x, t) dx dt =

  e. Ω = ∪ ζ∈F D ζ . For an internal edge ζ = P |Q, the set D ζ is supposed to be a subset of P ∪ Q and we define D P,ζ = D ζ ∩ P , so that D ζ = D P,ζ ∪ D Q,ζ (see Figure 2); for an external edge ζ of a cell P , D ζ is a subset of P , and D ζ = D P,ζ . The cells (D ζ ) ζ∈F are referred to as the dual or diamond cells, and D P,ζ as half dual cells or half diamond cells. For a rectangular cell P , we define D P,ζ as the simplex having the mass centre of P as vertex and the edge ζ as basis; this definition is extended to general primal meshes by supposing that |D P,ζ | is still equal to |P |/4

Figure 2 .

 2 Figure 2. Primal and dual meshes and associated notations for the quadrangular mesh and Rannacher-Turek like unknowns.

2 )

 2 ) denotes the set of vertical (resp. horizontal) edges. The cells (D ζ ) ζ∈F are still referred to as the dual cells. They are no longer diamond shaped; indeed, a half dual cell D P,ζ is now half of the rectangle P with side ζ (see Figure 3). As in the former case, for an internal edge ζ = P |Q, the dual cell D ζ is the subset of P ∪ Q defined as D ζ = D P,ζ ∪ D Q,ζ ; for an external edge ζ of a cell P , D ζ is a subset of P , and D ζ = D P,ζ .

N

  reads, with ϕ ∈ C ∞ c (Ω × [0, T )): m) ) n P -β(U (m) )(x, t) ϕ(x, t) dx dt =

Figure 4 .

 4 Figure 4. Left: the primal edges ζ and ζ ′ are adjacent. Right: the primal edges ζ and ζ ′ are opposite.

  as m tends to +∞. Reordering the summation in R

  if the weight ω ζ is controlled by both |P | and |Q|; since we have ω ζ ≤ 2 max(diam(P ), diam(Q))

2

 2 tends to zero when m tends to +∞ if ω η is controlled by both |D ζ | and |D ζ ′ |; this is indeed the case thanks to Assumption (28), since |D ζ | ≥ |P η |/4 and |D ζ ′ | ≥ |P η |/4.

  ζ and D ζ ′ , a weight given by ω η = diam(P η ) |ζ| + |ζ ′ | . Thus, again thanks to Lemma A.1, R (m) 2 tends to zero when m tends to +∞ since, remarking that |D ζ | ≥ |P η |/2, |D ζ ′ | ≥ |P η |/2 and ω η ≤ 2diam(P η ) 2 , so the weight ω η is controlled by both |D ζ | and |D ζ ′ | thanks to Assumption (28).

Remark 4 . 2 (

 42 On the required regularity of the time discretisation).

  for x ∈ P, P ∈ P, and t ∈ (t n , t n+1 ), n ∈ 0, N -1 ,

	with
	C(U ) n P = (ð t β) n P +	1 |P |	ζ∈F(P )	|ζ| F n ζ • n P,ζ ,
	where β n P , P ∈ P, n ∈ 0, N is a family of real numbers,
	(ð t β) n P =

  |

	tn+1 tn	P	B n P (x, t) dx dt +	tn+1 tn	P	R n P (x, t) dx dt ,
	and, thanks to (17),
	X

  L} ∈ S x and {p, q} ∈ S t , let

	1 t n+1 -t n p∈ 0,N -1	δ n,p .	(32)
	{n,p}∈St
	For {K,

d({K, L}) = max (x,y)∈K×L |y -x|, d({p, q}) = t q+1 -t p if q > p, t p+1 -t q otherwise

Appendix A. Convergence of discrete functions in L 1 We recall a result proven in [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF]Lemma 4.3]. To facilitate its use in the proof of Lemma 4.1, it is rephrased here under more general forms than in [START_REF] Gallouët | On the weak consistency of finite volume schemes for conservation laws on general meshes[END_REF] (see Remark A.2 below for the differences). and let

Then the following convergence result holds.

Lemma A. [START_REF] Eymard | Finite volume methods[END_REF]. Let (M (m) ) m∈N and (T (m) ) m∈N be a given sequence of meshes and time discretisations. Let us suppose there exists θ > 0 such that θ M (m) ≤ θ and θ T (m) ≤ θ for all m ∈ N, with θ M (m) and θ T (m) given by Equation (32). Let us assume that d(M (m) ) and d(T (m) ) tend to zero when m tends to +∞. Let u ∈ L 1 (Ω×(0, T )) and (u p ) p∈N be a sequence of functions of L 1 (Ω×(0, T )) such that u p → u in L 1 (Ω×(0, T )) as p → +∞. Then T M (m) ,T (m) u p defined by (31) tends to zero when m tends to +∞ uniformly with respect to p ∈ N.