Nina Aguillon

Christophe Chalons

Nondiffusive conservative schemes based on approximate Riemann solvers for Lagrangian gas dynamics

In this paper, we present a conservative finite volume scheme for the gas dynamics in Lagrangian coordinates, which is fast and nondiffusive. Fast, because it relies on an approximate Riemann solver, and hence the costly resolution of Riemann problems is avoided. Nondiffusive, because the solution is exact when the initial data is an admissible isolated shock, and discontinuities are sharply captured in general. The construction of the scheme uses two main tools: the extension to the barotropic Euler equations of the discontinuous reconstruction strategy presented in [Agu14], and the approximate Riemann solver of [CC14], which is exact on isolated admissible shocks. Numerical experiments in 1D and 2D are proposed.

Introduction

In this paper, we consider the barotropic gas dynamics equations in Lagrangian coordinates:

     ∂ t τ -∂ x u = 0, ∂ t u + ∂ x p(τ ) = 0, τ (0, x) = τ 0 (x), u(0, x) = u 0 (x), (1)
where the pressure law p is a strictly convex and strictly decreasing function. This system is well-known to be strictly hyperbolic, with eigenvalues ± -p ′ (τ ) and two genuinely nonlinear characteristic fields. The phase space Ω of this system is

Ω = {(τ, u) ∈ R 2 , τ > 0 }.
As usual, we supplement (1) with the validity of the so-called entropy inequality

∂ t U (τ, u) + ∂ x W(τ, u) ≤ 0, (2)
where (U , W) is the entropy-entropy flux pair given by U (τ, u) = u 2 2 + e(τ ), with e(τ ) = -τ p(y)dy 1 W(τ, u) = p(τ )u, and we consider entropy weak solutions of (1)-(2), see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] for example.

In order to motivate our study, let us temporarily assume that the initial data is a Riemann initial data

τ 0 (x) = τ L 1 x<0 + τ R 1 x>0 , u 0 (x) = u L 1 x<0 + u R 1 x>0 , (3)
where the left and right states (τ L , u L ) and (τ R , u R ) verify the Rankine-Hugoniot relations

u L -u R = -s(τ L -τ R ) and p(τ L ) -p(τ R ) = s(u L -u R ), (4)
for a given speed of propagation s, and the triple (τ L , u L ), (τ R , u R ), s satisfies (2) in the weak sense, i.e.

-s

(U (τ R , u R ) -U (τ L , u L )) + (W(τ R , u R ) -W(τ L , u L )) ≤ 0.
In that case the initial condition under consideration is such that the initial discontinuity will move at velocity s to form an admissible entropy weak solution of (1)-(2), see again [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF].

We say that we have an isolated shock wave.

On the other hand and from the numerical point of view, the Godunov method [God59] and Godunov-type methods based on exact or approximate Riemann solvers (see [START_REF] Eleuterio | Riemann solvers and numerical methods for fluid dynamics[END_REF] for a review) are certainly the most celebrated numerical schemes that provide good numerical approximations of (1)-(2) for general initial data. In particular, the Lax-Wendroff theorem (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]) ensures the convergence to the entropy weak solution of (1)-(2) provided that the method does converge and verifies a discrete version of (2).

When applied to the particular initial data leading to an isolated shock wave, such conservative finite volume schemes are well-known to introduce numerical diffusion, which means here intermediate values that do not correspond to values taken by the exact solution.

The basic motivation of this work is to propose a new conservative finite volume scheme which will be able to exactly capture isolated shock waves; and more generally entropy weak solutions of (1)-(2), with no numerical diffusion near shocks. By no numerical diffusion, we mean here that only one intermediate value will be present in shock profiles and importantly that this intermediate value will correspond to the average of the exact solution on a given cell for isolated shock wave. This is of course the best one can do.

To achieve this goal, let us first emphasize that the drawbacks of Godunov and Godunovtype methods are two fold.

• First, the numerical diffusion already discussed above leads to a loss of accuracy of the approximate solutions, especially around discontinuities.

• Second, the Godunov method uses the knowledge of the exact solution of (1)-( 2)-(3) with any given left and right states (τ L , u L ) and (τ R , u R ) in the phase space Ω. This exact Riemann solver is of course exact when (3) is an isolated shock, but it can be very expensive or difficult to compute the exact solution in general. When approximate Riemann solvers are used, it can be much less expensive but the property of being exact when (3) is an isolated shock is generally lost.

To avoid these issues, we follow the following strategy. Concerning the first point, i.e. to suppress the numerical diffusion, we use a discontinuous reconstruction strategy. The basic idea is to reconstruct entropy satisfying shocks inside some cells of the mesh in order to regain accuracy. The idea of reconstructed an initial data from the constant by cell function given by the scheme is not new. For example, in the MUSCL scheme [vL97], a linear by cell reconstruction of the solution is used. It allows to built schemes with high order on smooth regions. On the contrary, the discontinuous reconstruction strategy is not built to be precise in smooth areas, but yields to schemes that are exact on isolated shocks by the mean of allowing the reconstructed solution to have more discontinuities. This strategy was introduced in [START_REF] Boutin | Convergent and conservative schemes for nonclassical solutions based on kinetic relations[END_REF]. In the case of scalar advection, this scheme has an interpretation in terms of downhill decentering, and is originally presented that way in [START_REF] Després | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF]. The case of scalar conservation law has been deeply investigated in [START_REF] Lagoutière | Stability of reconstruction schemes for scalar hyperbolic conservation laws[END_REF], [LagXX], [START_REF] Bouchut | An antidiffusive entropy scheme for monotone scalar conservation laws[END_REF] and more recently in [START_REF] Chalons | A conservative scheme for non-classical solutions to a strongly coupled pde-ode problem a conservative scheme for non-classical solutions to a strongly coupled pde-ode problem[END_REF]. It has been extended to the capture of contact discontinuities in [START_REF] Després | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF] and [BFBC + 11]. In [START_REF] Alouges | The reservoir technique: a way to make Godunov-type schemes zero or very low diffuse. Application to Colella-Glaz solver[END_REF], [START_REF] Li | Numerically neither dissipative nor compressive scheme for linear advection equation and its application to the Euler system[END_REF], the reader can find other nondiffusive schemes based on discontinuous reconstructions. The ENO scheme with subcell resolution [START_REF] Harten | ENO schemes with subcell resolution[END_REF] is also based on a discontinuous reconstruction; this is the only other scheme we know which is also exact on isolated shocks. In [START_REF] Aguillon | A reconstruction scheme for the euler equations[END_REF], the discontinuous reconstruction strategy has been extended to the system of gas dynamics in Eulerien coordinates. Shocks are also sharply captured, and the proposed scheme is exact on isolated shock, but the exact Riemann solver is used. We follow this generalization to systems in the present work. Concerning the second point, i.e. to avoid the costly exact resolution of Riemann problems which is used in [START_REF] Aguillon | A reconstruction scheme for the euler equations[END_REF], we use the approximate solver proposed in [START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF]. Up to our knowledge, this is the first (and only) entropy consistent approximate Riemann solver which is exact for isolated shock waves. This is the key property an approximate Riemann solver should have to be coupled successfully to the discontinuous reconstruction strategy. This is exactly what we propose to do in the present paper.

To conclude this introduction, let us illustrate the fact that usual approximate Riemann solvers (abbreviated ARS in the sequel) are not exact for isolated shock waves. We consider here the Suliciu relaxation procedure [START_REF] Suliciu | On the thermodynamics of rate-type fluids and phase transitions. I. Rate-type fluids[END_REF] (see also [START_REF] Bouchut | Semi-discrete entropy satisfying approximate Riemann solvers. The case of the Suliciu relaxation approximation[END_REF]), which is also the starting point of the solver proposed in [START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF] and used in the present paper. The idea is to approach the solutions of (1) by the solutions of a larger but simpler system, namely

     ∂ t τ -∂ x u = 0, ∂ t u + ∂ x π = 0, ∂ t T = 1 ǫ (τ -T ), (5)
where π is related to the pressure p via the expansion

π = π(T ) = p(T ) + a 2 (T -τ ).
When ǫ tends to 0, we recover asymptotically the p-system (1) we are interested in, and for stability reason a must be chosen larger than the maximum speed of the wave appearing in the initial system (1), see the next section and [START_REF] Bouchut | A reduced stability condition for nonlinear relaxation to conservation laws[END_REF], [START_REF] Qiang Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF] and [START_REF] Shi | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF] for more details. This is the so-called subcharacteristic stability condition.

System (5) is actually simpler to solve than the original p-system (1) since the Riemann problem associated to (5), namely

               ∂ t τ -∂ x u = 0, ∂ t u + ∂ x π = 0, ∂ t T = 0, (τ, u, T )(t = 0, x) = (τ L , u L , T L )1 x<0 + (τ R , u R , T R )1 x≥0 , T L = τ L and T R = τ R , (6)
where the initial data is taken at equilibrium, can be explicitly solved. The characteristic fields are easily shown to be linearly degenerate, and the solution of (6) contains three contact discontinuities, propagating with velocities -a, 0, and a. We denote by (τ L, * , u L, * , T L, * ) the state on the left of the stationary wave, and by (τ R, * , u R, * , T R, * ) the state on its right (see Figure 1). Similarly we denote by The intermediate states can be obtained using the Rankine-Hugoniot relations across each contact discontinuities, and are given by

π L := π(T L ) = p(τ L ) and π R := π(T R ) = p(τ R ). u L τ L T L u R τ R T R u R * = u * u L * = u * T L * T R * τ L * τ R * a -a 0
                     u L * = u R * := u * = u L +u R 2 -π R -π L 2a , π L * = π R * := π * = π L +π R 2 -a 2 (u R -u L ), τ L * = τ L + u * -u L 2a , τ R * = τ R -u * -u R 2a , T L * = τ L , T R * = τ R . (7)
We refer once again to [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] for details. Importantly, since a is chosen larger than any speed of propagation of (1) for stability reasons, there is no hope to capture exactly isolated shocks with any Godunov-type method using this approximate Riemann solver. The idea of [START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF] is to modify the approximate Riemann solver given by (6) by introducing a new wave propagating with a velocity σ in order for the new ARS to be exact on isolated shock waves. The new wave pattern is depicted on Figure 2. Attached to this new wave is a parameter θ, which makes the link between the case of an isolated shock (θ = 1, the other three waves are trivial) and the classical solver (θ = 0, the new wave is trivial and the solver coincides with (7)). Moreover, this parameter can be chosen in such a way that the new ARS is entropy satisfying and exact on isolated shocks. The parameter θ can be seen as a detection parameter for shocks. Thus the shock detection step needed in [START_REF] Aguillon | A reconstruction scheme for the euler equations[END_REF] is already enclosed in the ARS of [START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF].

This paper is organized as follows. The derivation of the scheme is detailed in the first section. More precisely, we describe the approximate solver introduced in [CC14] and the discontinuous reconstruction strategy. Then we put together the two tools and present the numerical fluxes. We prove in Subsection 1.3 that the scheme is exact whenever the initial data is an isolated shock. In the Section 2 we present 1-dimensional simulations showing that the scheme has a very low numerical dissipation. In the third and last section, we extend the scheme to 2-dimensional situation by the mean of a directional splitting. Once again, discontinuities are sharply captured.

1 Construction of the scheme

The approximate Riemann solver

In this section we briefly describe the approximate Riemann solver introduced in [CC14] to solve (1)-( 2)-(3) and refer to this paper for details. The proposed approximate solution is the exact solution of the following system with Riemann initial data at equilibrium:

               ∂ t τ -∂ x u = 0, ∂ t u + ∂ x π(T ) = 0, ∂ t T = M(θ)δ x=σt , (τ, u, T )(t = 0, x) = (τ L , u L , T L )1 x<0 + (τ R , u R , T R )1 x>0 , T L = τ L and T R = τ R . (8)
System (8) is a modified version of (5), where the measure valued right hand side M(θ)δ x=σt allows the solution to jump along the line x = σt. The parameter σ is chosen such that if the initial datum is an isolated shock, σ is the exact speed of this shock. The Rankine-Hugoniot relations (4) yield

s 2 = - p L -p R τ L -τ R .
In [START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF] the authors then propose to set

σ = sign(τ R -τ L ) - p L -p R τ L -τ R , (9)
the role of the sign function is to distinguish in between entropy shocks of the first (resp. the second) eigenvalue -√ -p ′ (resp. √ -p ′ ), for which the shock has a negative velocity and τ R < τ L (resp. a positive velocity and τ R > τ L ). The parameter θ will be chosen so that the approximate Riemann solver is entropy satisfying. The solution of System (8) has four waves: the three usual waves having speeds -a, 0 and a, plus a wave at speed σ driven by the source term. According to [START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF], this σ-wave has to be understood as an approximation of the shock wave with the largest amplitude in the exact Riemann solution of (1)-( 2)-(3). This approximation turns out to be exact in the case of an isolated shock wave. The structure of the new ARS is depicted on Figure 2. To define the three intermediate states (and hence the nine unknowns), consistency relations and Rankine-Hugoniot relations are imposed by the authors in [START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF]. We note in particular that u and τ are constant through the stationary wave while the quantity denoted by

replacements u L u L τ L τ L π L π L u R u R τ R τ R π R π R u - u - u - u + u + u + π - π - π - π + π + π + τ 1 τ 1 τ 2 τ 2 τ 3 τ 3 a a -a -a 0 0 σ σ
I = π(T ) + a 2 τ = p(T ) + a 2 T (10)
is constant through the -a-and a-waves. We denote with a minus subscript (resp. a plus subscript) their values on the left (resp. on the right) of the σ-wave. The intermediate states are found by solving a linear system (in τ , u and I) and are given , when σ > 0, by

               u -= u * + σθ 2a (a -σ)(τ R -τ L ), u + = u * -σθ 2a (a + σ)(τ R -τ L ), τ 1 = τ L + 1 a (u --u L ), τ 3 = τ R + 1 a (u R -u + ), τ 2 = τ 3 -θ(τ R -τ L ), and
               π -= π * -σθ 2 (a -σ)(τ R -τ L ), π + = π * -σθ 2 (a + σ)(τ R -τ L ), I 1 = I L , I 2 = (1 -θ)I R + θI L , I 3 = I R .
(11) and when σ < 0, by

               u -= u * + σθ 2a (a -σ)(τ R -τ L ), u + = u * -σθ 2a (a + σ)(τ R -τ L ), τ 1 = τ L + 1 a (u --u L ), τ 3 = τ R + 1 a (u R -u + ), τ 2 = τ 1 + θ(τ R -τ L ),
and

               π -= π * -σθ 2 (a -σ)(τ R -τ L ), π + = π * -σθ 2 (a + σ)(τ R -τ L ), I 1 = I L , I 2 = (1 -θ)I L + θI R I 3 = I R . ( 12
)
The states denoted with a star subscript * correspond to intermediate states in the Suliciu ARS. Their expressions have been given in (7).

The key point of [START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF] is to pick θ is such a way that the solver satisfies a discrete entropy inequality and gives the exact solution when the initial data is an isolated shock propagating with speed s (in which case σ = s thanks to the choice (9)). Under the classical subcharacteristic condition

|a| > max τ ∈{τ L ,τ L * ,τ R * ,τ R } -p ′ (τ ), (13)
it can be achieved with the choice

σθ(I R -I L ) = max 0, min σ(I R -I L ), -2a(a 2 -σ 2 )A(v L , v R ), a 2 |σ|(a + |σ|) a + |σ|/2 (τ R * -τ L * )sign(σ) , (14)
where

A(v L , v R ) = τ R p(s) ds + π R 2a 2 - τ L p(s) ds -π L 2a 2 I R -I L - π * a 2 .
We recall below the main result of [START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF].

Theorem 1.1. Under the subcharacteristic condition (13), the approximate Riemann solver defined by (9), ( 11)-( 12) and ( 14) is conservative and entropy satisfying, preserves the phase space Ω, is Lipschitz continuous with respect to the initial Riemann data, and is exact on isolated shocks (which, again, means that if the initial condition is such that the exact solution is an isolated entropy shock wave, the proposed approximate solution coincides with the exact one).

Discontinuous reconstruction schemes

This approximate Riemann solver will be coupled with the discontinuous reconstruction strategy introduced in [BCLL08] on the scalar conservation law with increasing flux, and extended in [START_REF] Aguillon | A reconstruction scheme for the euler equations[END_REF] to the Euler equations using the exact Riemann solver. The aim of the discontinuous reconstruction strategy is to obtain nondiffusive conservative Godunov-type finite volume schemes, which moreover are exact for isolated and admissible shock waves. Let us first recall the formulation of finite volume schemes on a staggered grid. We denote by t 0 = 0 < t 1 < t 2 < • • • the time discretization and by ∆t n = t n+1 -t n the n-th time step. The real line is divided with cells that are always of length ∆x but the cells will move at each time step. In the sequel we will also propose a scheme on a fixed grid, but for the sake of clarity it seemed more natural to begin with a (fixed-size) staggered grid. We denote by x 0 j the centers of the cells at time t 0 and by x 0 j+1/2 = x 0 j + ∆x/2 their extremities. At the n-th iteration, a mesh velocity V n mesh is given, and the mesh moves from time to time according to the mesh velocity: x n+1 j = x n j + ∆t n V n mesh , and accordingly x n+1 j+1/2 = x n+1 j + ∆x 2 . We denote by U = (τ, u) the vector of conservative variables and by F = (-u, p(τ )) the flux. Integrating Equation (1) on the space-time slab

(t, x) ∈ R + × R : t n ≤ t < t n+1 , x n j-1/2 + V n mesh t ≤ x < x n j+1/2 + V n mesh t ,
we obtain the scheme

U n+1 j = U n j - ∆t n ∆x (F n j+1/2 -F n j-1/2 ), (15)
where U n j = (τ n j , u n j ) is supposed to be an approximation of the mean value of the exact solution at time t n on the j-th cell:

U n j ≈ 1 ∆x x n j+1/2 x n j-1/2 U (t n , x) dx,
and the numerical flux

F n j+1/2 = (F n,τ j+1/2 , F n,u j+1/2
) is an approximation of the exact flux along the line x n j+1/2 + V n mesh t:

F n j+1/2 ≈ 1 ∆t n t n+1 t n F (U (s, x n j+1/2 + V n mesh s)) -V n mesh U (s, x n j+1/2 + V n mesh s) ds.
The choice of a formula expressing F n j+1/2 as a function of the mean values (U n k ) k∈Z defines the finite volume scheme. In the discontinuous reconstruction scheme presented below, the numerical flux F n j+1/2 will depend on (U n k ) k∈{j-1,j,j+1} when V n mesh is nonpositive and on (U n k ) k∈{j,j+1,j+2} when V n mesh is nonnegative. The idea of the scheme is to reconstruct entropy satisfying shocks inside each cell of the mesh, using the neighboring cells, and let them evolve during a time ∆t n to compute the flux. The use of a staggered grid is crucial to avoid to deal with waves interactions, but as already stated, a natural version of the scheme on fixed grids will also be proposed. Now, we first address the CFL restriction of time step and then we describe precisely the derivation of the discontinuous reconstruction scheme based on the ARS described in the previous section.

Moving mesh and CFL condition

For all integer j, we solve the Riemann problem (R j ) for the augmented system (8), with initial data at equilibrium

(τ L , u L , T L ) = (τ n j-1 , u n j-1 , τ n j-1 ) and (τ R , u R , T R ) = (τ n j+1 , u n j+1 , τ n j+1 ), (16)
i.e. involving the neighboring cells [x n j-3/2 , x n j-1/2 ] and [x n j+1/2 , x n j+3/2 ]. We denote by ±a n j the speed of the extremal waves, such that (13) holds true, and by σ n j the speed of the additional wave (in accordance with the notation of Figure 2). We denote by (u n j,-, u n j,+ , π n j,-, π n j,+ , τ n j,1 , τ n j,2 , τ n j,3 ) the intermediate states appearing in the solution, defined by (11) or (12) (depending on the sign of σ n j ), with (16) for the left and right states. Then, we fix a mesh velocity V n mesh such that

|V n mesh | ≥ V n waves := max k∈Z a n k . ( 17
)
The time step ∆t is constrained by the CFL condition

∆t n ≤ ∆x (|V n mesh | + V n waves ) . ( 18
)
Conditions ( 17) and (18) insures the following facts, that will be widely used in the sequel.

• A wave emitted at time t n at point x n j+1/2 will not cross the interface x = x n j+1/2 + V n mesh t by (17). By (18), it will not cross the neighboring interfaces x = x n j-1/2 +V n mesh t and x = x n j+3/2 + V n mesh t within the time step either.

• More generally, a wave emitted somewhere in the cell [x n j-1/2 , x n j+1/2 ] at time t n can only cross the interface x = x n j+1/2 + V n mesh t if V mesh < 0 and the interface x = x n j-1/2 + V n mesh t if V mesh > 0. • If a wave emitted in the j-th cell interact with a wave emitted in the j + 1-th cell this interaction occurs inside a cell and the waves created by this interaction stay in that cell and do not interact with the interface x = x n j+1/2 + V n mesh t. This last point is true for scalar conservation law, because the velocity of the waves is smaller after the interaction than before. When considering systems of conservation law, this velocity can increase, so it is wise to overestimate V mesh at least for the first iterations in time.

In practice we take, for some µ in (0, 0.5) referred to as the Courant number,

∆t n = µ∆x V n waves ,
and V n mesh = (-1) n ∆x 2∆t n . Conditions (17) and ( 18) are fulfilled with this choice and we have

x n+1 j = x n j + (-1) n ∆x 2
. Note that the mesh moves alternatively with positive and negative speeds and has always cells of size ∆x.

Reconstruction of shocks

We denote by U n j,-= (τ n j,-, u n j,-) and U n j,+ = (τ n j,+ , u n j,+ ) the states around the σ n j -wave in the Riemann problem (R j ) (see Figure 2 for the link between τ ± and τ 1 , τ 2 , τ 3 in Formula (11) and ( 12)). The core of the discontinuous reconstruction strategy is to see the mean value U n j as the average of a shock between the state U n j,-and the state U n j,+ , located somewhere in the cell [x n j-1/2 , x n j+1/2 ] and propagating with velocity σ n j . This process and the notation are illustrated on Figure 3.

If we use the conservation of τ at time t n on the cell [x n j-1/2 , x n j+1/2 ] to place the reconstructed shock, we obtain that it should lie at a distance d n,τ j of x n j-1/2 such that

d n,τ j τ n j,-+ (∆x -d n,τ j )τ n j,+ = ∆xτ n j ,
which yields

d n,τ j = ∆x τ n j,+ -τ n j τ n j,+ -τ n j,- . (19)
Similarly, if we use the conservation of u, we obtain that it should lie at the distance are different in general, but we will see later on that the they coincide in the case of an isolated shock wave. If one of those distances is outside of the interval (0, ∆x), we consider that the mean value in cell j does not come from the average of a shock, and we do not perform any reconstruction within the cell [x n j-1/2 , x n j+1/2 ]. It yields the following definition.

d n,u j = ∆x u n j,+ -u n j u n j,+ -u n j,- (20)
x x x n j+1/2 x n j+1/2 x n j-1/2 x n j-1/2 τ τ n j τ n j,+ τ n j,- τ n j+1 τ n j-1 u u n j u n j,+ u n j,-
Definition 1.2. The left and right reconstructed states in the j-th cell at time t n , (τ n j,L , u n j,L ) and (τ n j,R , u n j,R ) , are defined as follow:

(τ n j,L , u n j,L ) = (τ n j,-, u n j,-) if d n,τ j ∈ [0, ∆x] and d n,u j ∈ [0, ∆x], (τ n j , u n j ) otherwise, (21)
and

(τ n j,R , u n j,R ) = (τ n j,+ , u n j,+ ) if d n,τ j ∈ [0, ∆x] and d n,u j ∈ [0, ∆x], (τ n j , u n j ) otherwise. ( 22
)

Computation of the flux

The fluxes are computed by letting the shocks reconstructed in the previous step evolve during the time ∆t n . To illustrate the interest of using a moving mesh, consider the case illustrated in Figure 4, where a 2-shock is reconstructed in the cell [x n j-1/2 , x n j+1/2 ] and a 1-shock is reconstructed in the cell [x n j+1/2 , x n j+3/2 ] (σ n j > 0 and σ n j+1 < 0). At time t n , both shocks are located near x n j+1/2 , and they will interact within the time step. As it can be seen on Figure 4, it is impossible to compute the flux on the interface x = x n j+1/2 (blue dotted line) without resolving the wave interaction. On the other hand, computing the flux along the interface x = x j+1/2 + V n mesh t (red dotted line) is much easier. Under Conditions (17) and (18), this flux is piecewise constant and its computation only requires the knowledge of the crossing time between the shock reconstructed in cell j when V n mesh < 0 (in cell j + 1 when V n mesh > 0) and the interface.

σ n j σ n j+1 t n t n+1
x n j+1/2

x n j+1/2

x n j+1/2

x n j+1/2

x n+1 j+1/2

x n j-1/2

x n j-1/2 U n j,R U n j,L U n j+1,R U n j+1,L U n j,R U n j,L U n j+1,R U n j+1,L
Figure 4: Top: Reconstruction of a 1-shock in the j-th cell (σ n j > 0) and of a 2-shock in the j + 1-th cell (σ n j+1 < 0). Bottom: in black, positions of the reconstructed shocks through times, and the waves emitted when they interact. It is easy to compute the flux through the red and green interfaces, but not through the blue one.

Let us focus on the case where V n mesh is negative. The two components of the flux are computed in the same way, except that the position of the reconstructed shock in the j-th cell at the beginning of the time step may be different for each conserved variable. More precisely, in order to compute the flux in τ , we consider that the shock reconstructed in cell [x n j-1/2 , x n j+1/2 ] is initially located at the distance d n,τ j of x n j+1/2 . We recall that this insures the conservation of τ inside the j-th cell. As V n mesh is negative, the reconstructed shock crosses the right interface x = x n j+1/2 + V n mesh t at the time T n,τ j that verifies (whatever the sign of σ n j is)

d n,τ j + σ n j T n,τ j = ∆x + V n mesh T n,τ j , which is T n,τ j = ∆x -d n,τ j σ n j -V n mesh . ( 23
)
As we already explained, under Conditions ( 17) and ( 18), the only wave that can cross the interface x = x n j+1/2 +V n mesh t during the time step is the shock reconstructed in the j-th cell. Thus the flux passing through this interface (in red on Figure 4) is simply -u j,R -V n mesh τ n j,R before the crossing time T n,τ j and -u j,L -V n mesh τ n j,L after it. Of course, T n,τ j can be larger than ∆t n , in which case the flux is -u j,R -V n mesh τ n j,R during the whole time step. Eventually the flux in τ writes

∆tF n,τ j+1/2 = -u n j,R -V n mesh τ n j,R min(∆t, T n,τ j ) + -u n j,L -V n mesh τ n j,L (∆t -min(∆t, T n,τ j )).
To compute the flux in the u variable, the only difference is that the shock reconstructed in the j-th cell is initially located, using the conservation of u instead of τ , at the distance d n,u j of x n j-1/2 . Thus the crossing time is

T n,u j = ∆x -d n,u j σ n j -V n mesh ,
and we obtain similarly the second component of the numerical flux. Finally we have

               ∆tF n,τ j+1/2 = -u n j,R -V n mesh τ n j,R min(∆t, T n,τ j ) + -u n j,L -V n mesh τ n j,L (∆t -min(∆t, T n,τ j )), ∆tF n,u j+1/2 = π n j,R -V n mesh u n j,R min(∆t, T n,u j ) + π n j,L -V n mesh u n j,L (∆t -min(∆t, T n,u j )). (24)
If V n mesh is positive, the only wave that can cross the interface x = x n j+1/2 + V n mesh t during the time step is the shock reconstructed in cell j + 1, which propagates with velocity σ n j+1 . The corresponding crossing times are

T n,τ j+1 = d n,τ j+1 V n mesh -σ n j+1 and T n,u j+1 = d n,u j+1 σ n j+1 -V n mesh .
In that case, the flux F n j+1/2 passing through the interface is once again piecewise constant, but this time is computed with U n j+1,L before the crossing time, and with U n j+1,R after it (green interface on Figure 4). It writes

               ∆tF n,τ j+1/2 = -u n j+1,L -V n mesh τ n j+1,L min(∆t, T n,τ j+1 ) + -u n j+1,R -V n mesh τ n j+1,R (∆t -min(∆t, T n,τ j+1 )), ∆tF n,u j+1/2 = π n j+1,L -V n mesh u n j+1,L min(∆t, T n,u j+1 ) + π n j+1,R -V n mesh u n j+1,R (∆t -min(∆t, T n,u j+1 )). (25)
Finally, the conservative variables are updated with (15). Remark 1.3. When no reconstruction is performed, i.e. when U n j,L = U n j,R = U n j , we have, whatever the values of d n j , σ n j and T n j are,

F n j+1/2 = F (U n j ) if V n mesh < 0, F n j-1/2 = F (U n j+1 ) if V n mesh > 0.
In other words in that case, the numerical fluxes ( 24) and (25) coincide with the staggered Lax-Friedrichs scheme, which is well-known to be stable under Conditions (17) and (18).

On a fixed grid

In this section we present a discontinuous reconstruction scheme on a fixed grid. In this section only we take x n+1 j+1/2 = x n j+1/2 and thus V n mesh = 0. There is at least three advantages in using a fixed grid:

• the Courant number can be taken in (0, 1), instead of in (0, 0.5) with a staggered grid.

Thus we can use time step twice as large;

• at equal Courant number, the scheme on a fixed grid is less diffusive than with a staggered grid;

• dealing with boundary conditions is easier.

The problem is to deal with interactions between reconstructed shocks near an interface. The basic idea is to cancel some reconstructions when it happens. Typically in the example of Figure 4, where two reconstructed shocks cross the blue interface within the time step, we decide to cancel those two reconstructions and to use a classical numerical flux. This will be sufficient to obtain the main property of the scheme, namely to be exact when the initial data is an isolated shock.

Introducing discussion

If a shock is reconstructed in the j-th cell and has a positive speed (σ n j > 0), it is tempting to use the flux Formula (24) which gives the flux passing through the j + 1/2 interface if the initial data were

   τ 0 (x) = τ n j,L × 1 x<x n j-1/2 +d n,τ j + τ n j,R × 1 x>x n j-1/2 +d n,τ j , u 0 (x) = u n j,L × 1 x<x n j-1/2 +d n,u j + u n j,R × 1 x>x n j-1/2 +d n,u j ,
(see what happens through the red interface of the bottom of Figure 4). However, if a shock is reconstructed in the j +1 cell and has a negative speed (σ n j+1 < 0), it is similarly tempting to use the flux Formula (25), which gives the flux passing through the j + 1/2 interface if the initial data were

   τ 0 (x) = τ n j+1,L × 1 x<x n j+1/2 +d n,τ j+1 + τ n j+1,R × 1 x>x n j+1/2 +d n,τ j+1 , u 0 (x) = u n j+1,L × 1 x<x n j+1/2 +d n,u j+1 + u n j+1,R × 1 x>x n j+1/2 +d n,u j+1 ,
(see what happens through the green interface of the bottom of Figure 4). In the variable τ , the two shocks interact at the time T n,inter,τ j+1/2 such that

d n,τ j + σ n j T n,inter,τ j+1/2 = ∆x + d n,τ j+1 + σ n j+1 T n,inter,τ j+1/2 ,
and we obtain a different interaction time in the variable u. We evaluate the time at which the two shocks interact by

T n,inter j+1/2 = min(T n,inter,τ j+1/2 , T n,inter,u j+1/2 ) = min

d n j+1 + ∆x -d n j s n j -s n j+1 ,
where the second minimum is taken on the two components of the vectors d n j and d n j+1 . The interaction can also occur when when σ n j and σ n j+1 have the same sign and σ n j > σ n j+1 . If this time is smaller than ∆t n , the shocks interact within the time step and the waves created by the resulting interaction are likely to meet the j + 1/2-th interface (which was not the case under Condition (17), see once again Figure 4). In that case, we decide to not to take into account the reconstructions in cells j and j + 1, and to simply use the Godunov-type flux (associated with (8)) instead of the reconstruction flux (24) or (25).

Definition of the fluxes

Let us now be more precise and introduce some notation. We denote by F GOD (U L , U R ) the Godunov-type scheme associated to the ARS (8):

F GOD (U L , U R ) = (-u -, π -) if σ > 0, (-u + , π + ) if σ ≤ 0, (26)
(see Figure 2). We also denote by F n,← j-1/2 the flux given by (25) and by F n,→ j+1/2 the flux given by (24). Moreover, we use the convention σ n j = 0 if no reconstruction is performed in cell j, and we extend the definition of the interaction time by

T n,inter j+1/2 =    min d n j+1 +∆x-d n j σ n j -σ n j+1 if σ n j > σ n j+1 , and σ n j+1 σ n j = 0 +∞ otherwise. ( 27
)
We propose the following fluxes:

F n j+1/2 =            F n,→ j+1/2 if σ n j > 0, T n,inter j+1/2 > ∆t and σ n j+1 ≥ 0, F n,← j+1/2 if σ n j+1 < 0, T n,inter j+1/2 > ∆t and σ n j ≤ 0, F GOD (U n j , U n j+1 ) if (σ n j > 0 and σ n j+1 < 0) or T n,inter j+1/2 ≤ ∆t, F GOD (U n j,R , U n j+1,L ) otherwise. (28)
In particular for simplicity, we cancel the reconstructions in cells j and j + 1 if σ n j > 0 and σ n j+1 < 0. We now justify the choice of flux F n,→ j+1/2 if σ n j > 0 and σ n j+1 = 0, which corresponds to the case where a shock with positive speed is reconstructed in the j-th cell but no shock is reconstructed in the j + 1-th cell. In that case the -a n j , 0 and σ n j -waves in the solution of (8) with the left state U n j,R and the right state U n j+1 are trivial, and only the a n j -wave is nontrivial. Thus, as a n j > σ n j and as depicted on Figure 5, no wave interaction occur at interface j + 1/2. The flux through the interface j + 1/2 is easily computed and given by F n,→ j+1/2 .

U n j-1 U n j+1 U n j+1 U n j,L U n j,L U n j,R U n j,R a n j a n j -a n j 0 σ n j σ n j t x x n j-1/2 x n j+1/2
Figure 5: Left: the structure of the ARS used to determine U n j,L and U n j,R . Right: the waves emitted in the reconstructed solution: there is only a wave at speed a n j between U n j,R and U n j+1 , which is faster than the reconstructed shock. Thus the flux through the interface j + 1/2 is easily computed and given by F n,→ j+1/2 .

We now prove a nontrivial property that states that if the initial data is an isolated shock wave, the discontinuous reconstruction scheme (on a staggered or on a fixed grid) is exact.

Exact approximation of isolated shocks

Let us start this section by the following proposition, which is easily verified.

Proposition 1.4. The numerical fluxes (24) and (25) (on a staggered grid) and (28) (on a fixed grid) are consistent: if

U n j-1 = U n j = U n j+1 = U n j+2 := U , then F n j+1/2 = f (U )-V n mesh U (we recall that on a fixed grid, V n mesh = 0).
The aim of this subsection is to prove the following theorem.

Theorem 1.5. The discontinuous reconstruction scheme is exact whenever the initial data is an isolated entropy satisfying shock. More precisely, suppose that the initial data is

τ 0 (x) = τ L 1 x<0 + τ R 1 x>0 u 0 (x) = u L 1 x<0 + u R 1 x>0
with u L > u R and (τ L , u L ) and (τ R , u R ) linked by the Rankine-Hugoniot relations (4). Then U n j is the average of the exact solution at time t n over the interval [x n j-1/2 , x n j+1/2 ] for all j in Z.

Remark 1.6. In particular, this theorem means that for such particular initial data, no spurious numerical diffusion is created by the scheme. For general initial data, the same behavior is observed, see Section 2 below.

Proof. We prove the result for a 1-shock (the reasoning is similar for a 2-shock), thus in the following we suppose that τ L τ R . Suppose that the property is verified for an integer n. Then, there exists an integer j n 0 and a distance d n j 0 ∈ [0, ∆x] such that

U n k =        U L if k < j n 0 , d j n 0 ∆x U L + ∆x-d j n 0 ∆x U R if k = j n 0 , U R if k > j n 0 .
By construction, the approximate Riemann solver of [START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF] is exact on (R j n 0 ). Therefore, the shock is correctly reconstructed in cell j n 0 , it by (9) has the correct speed

σ := u L -u R τ R -τ L .
The two distances d τ , and correspond to the exact position of the shock at time t n . On the other hand, no reconstruction is performed in cells j n 0 -1 and j n 0 + 1, because it is impossible to have d n,τ

j n 0 -1 or d n,τ j n 0 +1 in (0, ∆x).
Let us indeed explain in detail the case of the cell j n 0 -1. Consider the Riemann problem between the left state (τ L , u L ) and (τ R , ũR ) :

= (ατ R + (1 -α)τ L , αu R + (1 -α)u L ), where α = d n j n 0
∆x . We denote by σ the speed of the shock wave in this Riemann problem, which by ( 9) is given by

σ := p L -pR τR -τ L
, where pR = p(τ R ).

As τ + ≥ τ -, if τ -is larger than τ L , it is impossible to reconstruct τ in cell j n 0 -1 in a conservative manner. The difference τ --τ L expresses as a function of θ: τ -(θ) -τ L = τ + (θ) -θ(τ R -τ L ) -τ L = τR + 1 a (ũ R -u + (θ)) -θ(τ R -τ L ) -τ L = (1 -θ)(τ R -τ L ) + 1 a ũR -u * + σθ 2a (a + σ)(τ R -τ L ) = (τ R -τ L ) + ũR -u * a -θ(τ R -τ L ) 1 - σ 2a 2 (a + σ) .
This function is decreasing with respect to θ. Indeed, τR is larger than τ L , and as a is larger than σ we have

σ(a + σ) 2a 2 ≤ 1.
Thus for all θ in [0, 1],

τ -(θ) -τ L = τ -(θ = 1) -τ L = 1 a ũR -u * + σ 2a (a + σ)(τ R -τ L ) .
Let us prove that the quantity

Q := ũR -u * + σ 2a (a + σ)(τ R -τ L )
is positive. Replacing u * by its value, we obtain

Q = ũR -u L 2 + πR -π L 2a + σ 2a (a + σ)(τ R -τ L ).
Now, we divide by τR -τ L and use the definition of σ to write

2Q τR -τ L = -σ - σ2 a + σ a (a + σ) = σ -σ
We conclude by remarking that as the pressure is convex,

πR = p(ατ R + (1 -α)τ L ) ≤ απ R + (1 -α)π L and thus σ = π L -πR τR -τ L ≥ π L -π R τ R -τ L = σ.
Proving that no reconstruction is performed in cell j n 0 +1 follows the same lines. We now focus on the Riemann problem between the left state (τ L , ũL ) :

= (ατ L + (1 -α)τ R , αu L + (1 -α)u R ) and (τ R , u R ). We have u + -u R ≥ u * - σ 2a (a + σ)(τ R -τL ) -u R = ũL -u R 2 - π R -πL 2a - σ 2a (a + σ)(τ R -τL ) = τ R -τL 2 σ + σ2 a - σ a (a + σ) = τ R -τL 2 (σ -σ) . It follows that u + is larger than u R , because σ = πL -π R τ R -τL ≤ π L -π R τ R -τ L = σ. Thus, τ + = τ R + 1 a (u R -u + )
is smaller than τ L , and as τ -is smaller than τ + , it is impossible to reconstruct τ .

The rest of the proof consists in an elementary checking that ( 24) and (25) give the correct flux, it can be found in [START_REF] Aguillon | A reconstruction scheme for the euler equations[END_REF]. The case of a fixed grid follows because no reconstruction is cancelled.

2 One dimensional numerical simulations

Riemann problems

To begin with we test the discontinuous reconstruction scheme on a staggered grid. The space interval [-1, 1] is discretized with 100 cells and the pressure law is p(τ ) = τ -2 . We take ∆t n = 0.45 ∆x V n waves and V n mesh = (-1) n ∆x 2∆t n , thus both (17) and (18) hold. In each cases, we compare the discontinuous reconstruction scheme and the Godunov-type scheme, both based on the approximate Riemann solver of [START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF]. The flux of the "Godunov-type scheme" is computed on a fixed grid (V n mesh = 0) and its expression is given in (26).

The first three test cases are taken from [START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF], the fourth one is a variant of the fast shock of Arora and Roe [START_REF] Arora | On postshock oscillations due to shock capturing schemes in unsteady flows[END_REF]. On those test cases, the results given by the proposed schemes on staggered and fixed grids are very close to each other (the scheme on a fixed grid being less diffusive in general, but with almost no impact here). This is why we only show the results given by the scheme on a staggered grid.

Test 1: Isolated shock

The first Riemann problem is

τ L = 1, u L = 0, τ R = 0 and u R = - √ 3/2.
- It corresponds to an isolated shock. Theorem 1.5 is illustrated on Figure 6: the shock is perfectly advected by the reconstruction scheme, while it is diffused by the Godunov-type scheme. On the bottom of this figure, we can see that there is only one intermediate value in the shock profile given by the reconstruction scheme, which corresponds to the average of the exact solution on the cell.

Test 2: Rarefaction and shock

The second Riemann problem corresponds to τ L = 0.3, u L = 0, τ R = 0.6 and u R = 0.

It contains a 1-rarefaction wave and a 2-shock. On Figure 7, we can see that the shock is sharply captured by the reconstruction scheme.

- For this test case we have

τ L = 1, u L = 1, τ R = 52 and u R = -8,
which is a slight modification of the isolated fast shock considered in [START_REF] Arora | On postshock oscillations due to shock capturing schemes in unsteady flows[END_REF]. When computing a strong shock (also called fast shock), spurious oscillations appear in the velocity when using usual Godunov-type methods. They also appear clearly in the characteristic variable u ± 8/τ . Those two quantities are plotted on Figure 9. Theorem 1.5 states that the scheme is exact on isolated strong shocks. We can see on Figure 9 that the ability of the reconstruction scheme to approach correctly fast shocks is not lost when we introduce a perturbation (here, a second wave in the Riemann problem).

-1 This test case consists in reversing the time starting from a developed 1-rarefaction wave, to recover a discontinuity. This is a difficult test case. More precisely, let us denote

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 X U -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -9 -8 -7 -6 -5 -4 -3 -2 -1 X U-sqrt(8/Tau) Godunov Reconstruction Exact
(t, x) → U rar (x/t; U L , U R ) = (τ rar (x/t; U L , U R ), u rar (x/t; U L , U R ))
the self-similar rarefaction wave associated with the Riemann initial condition

U L 1 x<0 + U R 1 x≥0 , with U L = (τ L , u L ), U R = (τ R , u R ) and τ L = 0.5, u L = 10, τ R = 5 and u R = u L - √ 2(τ -1/2 R -τ -1/2 L ).
For this test case, the initial data is τ 0 (x) = τ rar (-x/T rar ) u 0 (x) = u rar (-x/T rar )

with T rar = 0.2. This initial condition is plotted on Figure 10. For times t < T rar , the exact solution is U (t, x) = U rar (-x/(T rar -t)), and for times t ≥ T rar , the solution coincides with the solution of the Riemann problem associated with left state U R and right state U L (at time t -T rar ), which contains a 1-shock and a 2rarefaction. In particular at time t = T rar the exact solution is nothing but a discontinuity between U R and U L and located at x = 0.

We compare the reconstruction scheme on a staggered grid (abbreviated RecSG in the legends) and on a fixed grid (abbreviated RecFG) at five different times. This is the only test case where the results are significantly different for the two schemes. On Figures 11 and12 we plot the solutions for t = 0.14, t = 0.17, t = 0.2, t = 0.23 and t = 0.26. We observe a spike near the discontinuity when the time t is close to T rar . It seems that the reconstruction scheme starts to reconstruct a little before t = T rar the expected shock that will appear in the exact solution for times t > T rar , and that will be eventually sharply computed by the proposed scheme. The maximal height of the spike is at t = T rar , and it diminishes rapidly afterwards, the shock being correctly approximated after T rar . The spike does not prevent the scheme from converging in L 1 , and its size is independent from the Courant number. On this test case, we observe that on the left of the solution, the results given on a fixed grid are much better than the ones given on a staggered grid. On the other hand, we note that the Godunov-type scheme does not create any spike, but introduces numerical diffusion on the shock wave for times t > T rar , as it is expected. converge toward the correct solution for the L 1 -norm. The order of convergence is plotted on Figure 14. Note also that even though the scheme is not entropy satisfying, it contains entropy information as the underlying approximate Riemann solver of [START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF] is entropy satisfying. It is thus likely that the problem comes from the discontinuous reconstruction strategy. More precisely, we believe that some reconstructions are not entropy satisfying and should be cancelled using a stronger constraint than the ones proposed in (21)-(22) (i.e. d n,τ j ∈ [0, ∆x] and d n,u j ∈ [0, ∆x]). For example in this case, we expect the reconstructed solution U n rec = (τ n rec , u n rec ) to be monotonous for times t < T rar , but as depicted on Figure 15, this is not the case at all in the course of the first iteration (n = 0).

In the scalar case, entropy satisfying versions of the anti-diffusive scheme of [DL01] (which can be reinterpreted in terms of discontinuous reconstruction) have been proposed in [START_REF] Bouchut | An antidiffusive entropy scheme for monotone scalar conservation laws[END_REF], but remain to be adapted to the present setting. It is a challenging and open problem at the stage of the present work.

Two dimensional simulation

In two space dimensions, the p-system writes

           ∂ t τ -∂ x u -∂ y v = 0, ∂ t u + ∂ x p(τ ) = 0, ∂ t v + ∂ y p(τ ) = 0, τ (
0, x) = τ 0 (x), u(0, x) = u 0 (x), v(0, x) = v 0 (x).

(29)

The velocity of the fluid is (u, v). We use a directional splitting to approach the solution of (29). The space domain is the square [-1, 1] 2 with rigid walls. We take p(τ ) = τ -2 , rec , i.e. the solution reconstructed from the initial data before any iteration in time is done. We take T rar = 0.015 and 500 cells. In blue, the piecewise constant function corresponding to (U 0 j ) j∈ Z which corresponds to the pressure law for the shallow-water equation written in Lagrangian coordinates. The initial data is τ 0 (x, y) = 1 √ x 2 +y 2 ≥0.5 + 2 * 1 √ x 2 +y 2 <0.5 , u 0 (x, y) = v 0 (x, y) = 0.

We compare the solution given by the discontinuous reconstruction scheme on a fixed grid and the Godunov-type scheme at different times. The specific volume τ is represented on Figures 16, 17 and 18. In that case, it corresponds to the inverse of the height of water. The discontinuities are sharply captured with the reconstruction scheme, with much more details captured with this method.

Figure 1 :

 1 Figure 1: The approximate Riemann solver of Suliciu.

Figure 2 :

 2 Figure 2: Structure of the Riemann solution depending on the sign of σ.

Figure 3 :

 3 Figure 3: Toeq: Reconstruction of a shock in the j-th cell.

  19) and (20) are both equal to d n j n 0

Figure 6 :

 6 Figure 6: Test 1 : specific volume (top left) and velocity (top right) at time t = 0.15. Bottom: close-up on the piecewise constant shock profile.

Figure 7 :Figure 8 :

 78 Figure 7: Test 2 : specific volume (top left) and velocity (top right) at time t = 0.04. Bottom: close-up on the piecewise constant shock profile.

Figure 9 :

 9 Figure 9: Test 4 : velocity (left) and u -8/τ (right) at time t = 0.3.

Figure 10 :

 10 Figure 10: Initial condition for test 5.

Figure 11 :Figure 12 :Figure 14 :

 111214 Figure11: Solution of test 5 with, from top to bottom, t = 0.14 and t = 0.17.

Figure 15 :

 15 Figure15: In black, the reconstructed solution U 0 rec , i.e. the solution reconstructed from the initial data before any iteration in time is done. We take T rar = 0.015 and 500 cells. In blue, the piecewise constant function corresponding to (U 0 j ) j∈ Z

Figure 16 :

 16 Figure 16: Comparison of τ for the discontinuous reconstruction scheme (left) and the Godunov-type scheme (right) at time t = 3.05 for meshes with (from top to bottom): 200 × 200, 400 × 400 and 800 × 800 cells

Figure 18 :

 18 Figure 18: Comparison of τ for the discontinuous reconstruction scheme (left) and the Godunov-type scheme (right) at time t = 4.58 for meshes with (from top to bottom): 200 × 200, 400 × 400 and 800 × 800 cells.

test case. More precisely, we plot the evolution of the quantity

where

are the reconstructed solution at time t n (we recall than on a fixed grid, V n-1 mesh = 0 for all n). On Figure 13, we can see that this quantity oscillates around 0 and often takes positive values until the time T rar = 0.2. It indicates that the scheme is not entropy satisfying in the strict sense, but only in a "weak" sense to be defined. Indeed, the scheme seems to


Vincenzo Zambianchi

Michel Kieffer

Francesca Bassi

Gianni Pasolini

Davide Dardari

Distributed SPS Algorithms for Non-Asymptotic Confidence Region Evaluation

In this paper, the distributed computation of confidence regions for parameter estimation is considered. Some information diffusion strategies are proposed and compared in terms of the required number of data exchanges to get the corresponding region. The effects of algorithms truncation is also addressed. As support for the theoretical part, numerical results are presented.

I. INTRODUCTION

A wireless sensor network (WSN) is composed of powerconstrained sensing devices, cooperating in processing local measurements to complete a common task [START_REF] Verdone | Wireless Sensor and Actuator Networks: technologies, analysis and design[END_REF]. The task can be simple environmental monitoring (e.g., temperature, pollution levels) or the evaluation of a complex function of the measurements (e.g. localization of a source of signal [START_REF] Matamoros | On the estimation of randomly sampled 2D spatial fields under bandwidth constraints[END_REF]). The reference scenario considered here is the estimation of a parameter vector p * from the set of local and noisy measurements at the nodes. To provide coherent estimates at all nodes, robustness against equipment failure, and scalability versus network size, a distributed protocol working in absence of a central processing unit is desirable.

In [START_REF] Schizas | Distributed LMS for consensus-based in-network adaptive processing[END_REF], [START_REF] Mateos | Distributed recursive least-squares for consensus-based in-network adaptive estimation[END_REF] this is accomplished via recursive weighted least square estimation, where a consensus-based algorithm allows to incorporate information from neighbouring nodes in the local estimate. A similar approach is taken for the Bayesian framework in [START_REF] Olfati-Saber | Kalman-consensus filter: Optimality, stability, and performance[END_REF], where consensus-based distributed Kalman filtering is proposed. The algorithms in [START_REF] Schizas | Distributed LMS for consensus-based in-network adaptive processing[END_REF]- [START_REF] Olfati-Saber | Kalman-consensus filter: Optimality, stability, and performance[END_REF] converge asymptotically in time, but it is not possible to predict the quality of the local estimates in the non-asymptotic scenario without statistical knowledge of noise distribution.

In some applications (e.g., in source localization) the characterization of the confidence region associated to the estimate can be more valuable than the evaluation of the point estimate itself. In this work we thus address the problem of the evaluation of non-asymptotic confidence regions for the estimation of p * . Non-asymptotic indicates that the estimate is evaluated from a finite amount of measurements (unlike classical Cramer-Rao bound-like asymptotic confidence regions).

In the presence of a central processing unit the characterization of the non-asymptotic confidence regions can be accomplished using the results in [START_REF] Campi | Guaranteed non-asymptotic confidence regions in system identification[END_REF]- [START_REF] Kieffer | Guaranteed characterization of exact nonasymptotic confidence regions as defined by {LSCR} and {SPS}[END_REF]. Specifically, the methods in [START_REF] Campi | Guaranteed non-asymptotic confidence regions in system identification[END_REF], [START_REF] Campi | Non-asymptotic confidence regions for model parameters in the presence of unmodelled dynamics[END_REF] allow to derive a lower bound on the probability that p * falls within the confidence region, whereas Vincenzo Zambianchi, Davide Dardari and Gianni Pasolini are with CNIT, University of Bologna, Italy, Michel Kieffer is with L2S, CNRS, SUPELEC, Univ. Paris-Sud, Gif-sur-Yvette, France. Francesca Bassi is with ESME-Sudria, Ivry-sur-Seine, France. This work has been supported by the European commission within the NoE Newcom#.

the exact probability can be obtained using the sign perturbed sums (SPS) algorithm in [START_REF] Csàji | Non-asymptotic confidence regions for the least-squares estimate[END_REF]. In [START_REF] Kieffer | Guaranteed characterization of exact nonasymptotic confidence regions as defined by {LSCR} and {SPS}[END_REF] efficient computation of the confidence regions is obtained using interval analysis techniques. These methods do not require precise statistical knowledge of the noise, and work under very mild assumptions on its distribution.

The aim of this work is to provide fully distributed solutions for the in-node computation of the non asymptotic confidence regions, based on the SPS algorithm [START_REF] Csàji | Non-asymptotic confidence regions for the least-squares estimate[END_REF]. The SPS method is easily amenable to a distributed form, which converges in time to the centralized algorithm. We consider three distributed approaches (data flooding and parallel in-node processing, distributed processing via average consensus, and mixed flood-ing+consensus) able to provide exact confidence regions at each node, even for early truncation of the communication protocol. The performance is evaluated, for a target confidence level, in terms of the size of the confidence regions expressed as a function of the amount of data to be exchanged, and hence the energetic cost. This allows to compare the energetic efficiency of the proposed solutions.

II. THE CENTRALIZED SPS ALGORITHM

Let us recall, first of all, the centralized SPS algorithm [START_REF] Csàji | Non-asymptotic confidence regions for the least-squares estimate[END_REF] for the computation of non-asymptotic confidence regions. Consider a network of N nodes taking different noisy observations {y i }, for i = 1, 2, . . . , N , of a parameter p * belonging to the parameter space P, having dimension n p . The following local linear scalar measurement model is considered for the ith node:

y i = ψ T i p * + w i , (1)
where ψ i is the regressor vector, supposed known only at the corresponding node i, and w i is the measurement noise. The only assumption we make on the measurement noise is that it is independent from node to node and its distribution is symmetric with respect to zero, this setting a quite general domain of applications.

The SPS algorithm returns the exact confidence region around the least squares estimate p of p * , obtained as the solution of the normal equations N i=1 ψ i y i -ψ T i p = 0. Specifically, define the unperturbed sum

s 0 (p) = N i=1 ψ i y i -ψ T i p (2)
and the m -1 sign perturbed sums

s j (p) = N i=1 α j,i ψ i y i -ψ T i p , (3)
where j = 1, . . . , m -1 and α j,i ∈ {±1} are independent, identically distributed and equiprobable random signs, drawn at each node i. Given the squared norm of all previous quantities

z j (p) = ||s j (p)|| 2 2 , j = 0, . . . , m -1, (4)
the confidence region Σ q for the true parameter value p * is defined as

Σ q = {p ∈ P|z 0 (p) is not among the q largest z j (p)} =    p ∈ P m-1 j=1 τ j (p) ≥ q    , (5)
where τ j (p) = 1 if z j (p) -z 0 (p) > 0 and 0 otherwise. The rationale is that, for the true parameter value, the unperturbed sum is closely related to the measurement noise and it is therefore small, while for parameter values that are far from the true one, the perturbed sums are larger than the unperturbed one. In [START_REF] Csàji | Non-asymptotic confidence regions for the least-squares estimate[END_REF] it was proven that the confidence level is

Prob(p * ∈ Σ q ) = 1 - q m . (6)
It is worth noting that this algorithm requires a central entity with full knowledge of regressors ψ i , measurements y i , and coefficients α j,i . In the following section we will investigate three different solutions to overcome this limitation and let each node capable of computing (5) locally.

III. INFORMATION DIFFUSION STRATEGIES

This section presents three different solutions for the distributed evaluation of confidence regions. In all cases no central processing unit is assumed, hence confidence regions are locally evaluated at each node.

The first solution here investigated is based on the flooding of all required quantities across the entire network and the subsequent in-node computation of the SPS algorithm. This will not only ensure that the computed confidence region is the same at all nodes, but also that it is exactly the same as in the centralized version of the algorithm.

In the second solution, the information needed for the in-node computation of the SPS algorithm are, instead, exchanged according to an average consensus algorithm [START_REF] Xiao | A scheme for robust distributed sensor fusion based on average consensus[END_REF].

The third solution, whose details will be provided in the following, is a combination of flooding and consensus.

Let us observe, finally, that for energy saving reasons, truncation of the information diffusion algorithms is often desirable. The effects of truncations in the considered algorithms is also investigated in subsection III-D. This constitutes, to the author knowledge, a relevant advance with respect to the current literature.

A. Information Diffusion Through Flooding

In a flooding algorithm, the generic node i initially broadcasts its own state x (0) i = ψ T i , y i at time 0. At the successive step, it collects the states from neighbors, then forwards a new data packet containing the aggregated state x

(1) i = ψ T i , y i , ψ T j , y j j∈Ni
, where ψ T j , y j j∈Ni denotes the set of regressors and measurements pairs collected from nodes j in the neighborhood N i of node i. To avoid network collapse due to cycles, duplicated states are discarded. This process is repeated until each node in the network has collected the state from all nodes. Afterwards, each node is able to compute the perturbed and unperturbed sums in ( 2) and (3) for any p, and hence derive the confidence region. No transmission of α j,i is necessary, provided that all nodes agree on seeds for their random generators. Note that each node has to transmit a packet containing D

B. Information Diffusion Through Consensus

Looking at (2) and ( 3) it is evident that the computation of s 0 (p) and s j (p) does not necessarily require the knowledge of each single quantity, but rather of aggregated values. This suggests the adoption of consensus strategies that are well suited for distributed computation of sums or averages.

A consensus scheme may be viewed as the following discrete time evolving system [10]- [START_REF] Olfati-Saber | Consensus and cooperation in networked multi-agent systems[END_REF] x (k+1) = W x (k) [START_REF] Campi | Non-asymptotic confidence regions for model parameters in the presence of unmodelled dynamics[END_REF] where

x (k) = x (k) 1 , x (k) 2 , . . . , x (k) N T
is the global state at time k, whose ith entry is the local state at node i. W is the system dynamic matrix that depends on the consensus algorithm as well as on network topology. Proper structures of W asymptotically lead to a global state whose entries are all equal to the average of the initial quantities x (0) [START_REF] Xiao | Fast linear iterations for distributed averaging[END_REF], [START_REF] Olfati-Saber | Consensus and cooperation in networked multi-agent systems[END_REF]. In particular, the convergence of the consensus algorithm is assured by the following three necessary and sufficient conditions:

1 T W = 1 T , W1 = 1, ρ(W -11 T /N ) < 1 (8)
where ρ denotes the spectral radius (see [START_REF] Xiao | Fast linear iterations for distributed averaging[END_REF]). We choose the Metropolis matrix, which allows a faster convergence than others. Each entry of the matrix is

w i,j =        1 1+max{di,dj } if (i, j) ∈ E, 1 - (i,q)∈E w i,q if i = j, 0 otherwise, (9)
where d i = |N i | denotes the degree of node i (i.e., the number of neighbors) and E is the set of edges in the network topology.

According to [START_REF] Campi | Non-asymptotic confidence regions for model parameters in the presence of unmodelled dynamics[END_REF], each node performs its local state update as follows

x (k+1) i = N j=1 w i,j x (k) j . (10)
To the purpose of applying the described consensus algorithm to the problem of distributed calculation of ( 2) and (3), define the following averages that do not depend on

p b 0 = 1 N N i=1 ψ i y i b j = 1 N N i=1 α j,i ψ i y i (11)
A 0 = 1 N N i=1 ψ i ψ T i A j = 1 N N i=1 α j,i ψ i ψ T i ( 12
)
for j = 1, 2, . . . , m -1.

The consensus algorithm is launched on all A j 's and b j 's (including b 0 and A 0 ). At step k = 0, the local state at node i is given by

x (0) i = (ψ i yi) T , ψ i ψ T i , (αj,iψ i yi) T j , αj,iψ i ψ T
i j , with j = 1, 2, . . . , m -1, that is, the single addends in [START_REF] Xiao | Fast linear iterations for distributed averaging[END_REF] and [START_REF] Olfati-Saber | Consensus and cooperation in networked multi-agent systems[END_REF]. At each successive step the ith node updates its own state according to [START_REF] Xiao | A scheme for robust distributed sensor fusion based on average consensus[END_REF]. Once a consensus on A j and b j is reached, each node is able to locally evaluate (2) and (3) for any value p in the parameter search space. We want to underline that no particular value of the parameter p has to be transmitted and that getting averages, instead of the true sums, does not affect the SPS algorithm since the comparison of rescaled norms or norms gives the same ordering in [START_REF] Olfati-Saber | Kalman-consensus filter: Optimality, stability, and performance[END_REF]. Therefore the algorithm works also without the knowledge of N . Note that the state is composed of D cons = m(3n p + n 2 p )/2 values, where symmetry of ψ i ψ T i is exploited. The state dimension is constant during the entire running of the algorithm but is larger than that initially required by flooding.

C. Mixed Flooding+Consensus Approach

As observed, the pure average consensus algorithm requires an amount of data to be transmitted that is initially strictly larger than that required by flooding. However, in few iterations the amount of data transmitted with flooding exceeds D cons due to data accumulation at nodes. This fact suggested the introduction of a third mixed strategy conceived as follows: flooding, as described in Sec. III-A, works until the amount of data flooded by nodes exceeds D cons . Afterwards, nodes switch to a consensus strategy as the one described in Sec. III-B. The correct initialization for the consensus state x (0) i , in this case, will be obtained as the average of the already received quantities during the initial flooding. The numerical results will show the benefits of this approach when the maximum amount of data exchanged is limited due to energy constraints.

D. Effects of Information Diffusion Truncation

When energy or traffic constraints are present, as typically happens in WSNs, whatever algorithm is used it has to be truncated after a certain number of steps or when the total amount of data exchanged has reached a maximum tolerable threshold. In this case, the confidence region is derived before collecting the whole set of measurements and regressors, in the flooding approach, or before reaching the consensus on the A j 's and b j 's, in the consensus approach. Our aim is to determine the characteristics of the confidence region defined by SPS when truncation is considered. Truncating the consensus algorithm entails that ( 2) and ( 3) are estimated taking into account only the data actually received at node k, that is:

s k,0 (p) = N i=1 c k,i ψ i y i -ψ T i p (13)
s k,j (p) = N i=1 c k,i α j,i ψ i y i -ψ T i p , (14)
where j = 1, . . . , m -1, and c k,i ∈ [0, 1]. Coefficients c k,i do not depend on index j since they depend only on the effectively covered communication links, that are the same for all j's.

As for the truncated flooding algorithm, the estimations of (2) and ( 3) have the same expressions ( 13) and ( 14) that hold for the consensus case, but with c k,i being exactly 1 or 0 if data corresponding to node i have, respectively, been received or not at node k.

In any case, ( 13) is the normal equation that would be obtained in a centralized context, considering a weighted least-squares estimator, with a diagonal weight matrix C k = diag (c k,1 , . . . , c k,N ). Similarly, ( 14) is the sign perturbed sum that would be obtained when considering weighted leastsquares. The confidence region, obtained considering ( 13) and ( 14) in ( 5), is a non-asymptotic confidence region associated to the weighted least-squares estimate. Upon reaching completion of the flooding or convergence of the consensus algorithm, the c k,i are either all equal to one (flooding) or equal to 1/N (consensus). As a consequence, in a distributed setup, the same asymptotic result as in the centralized approach can be obtained (asymptotic equivalence). This is no more true when the truncated flooding or consensus algorithms are adopted.

Fortunately, it is possible to assert that any time-truncation of the information-diffusing algorithms still produces confidence regions with the same level of confidence of the asymptotic one. Extended details for this will appear in a companion paper; for now it suffices to say that following a similar approach as in [START_REF] Csàji | Non-asymptotic confidence regions for the least-squares estimate[END_REF], the evaluation of ( 13) and ( 14)

for p * gives ŝk,0 (p * ) = N i=1 c k,i ψ i w i (15)
and

ŝk,j (p * ) = N i=1 c k,i α j,i ψ i w i , (16)
with j = 1, . . . , m -1. The truncation effect is therefore a deterministic rescaling of measurement noises w i , since it only depends on the number of performed iterations. This rescaling preserves independence as well as symmetry of noise distribution which are the only hypothesis necessary to make ( 5) and ( 6) valid. The conclusion is that the confidence region obtained in the presence of truncation is a region that contains the true parameter value with the same exact level of confidence as in the not truncated case, the only difference being the region shape. This will be further investigated in the numerical results and is a consequence of the fact that the asymptotic and the truncated amounts of collected information are not the same.

IV. NUMERICAL RESULTS

In order to show the effects of truncation on the computation of confidence regions and compare the proposed distributed schemes, simulations were performed using Matlab along with the Intlab package [START_REF] Rump | INTLAB -INTerval LABoratory[END_REF].

We considered a network of N = 100 nodes randomly deployed with uniform distribution in a square area of unit width, as shown in Fig. 1. Each node takes a single scalar measurement according to [START_REF] Verdone | Wireless Sensor and Actuator Networks: technologies, analysis and design[END_REF] distance are considered to be communicating in a error-free way.

Figs. 2 and3 show the 90% confidence regions (q = 1, m = 10), computed at node #1, after 4 and 30 consensus iterations, respectively. For efficient computation of the regions, the method of contractors described in [START_REF] Kieffer | Guaranteed characterization of exact nonasymptotic confidence regions as defined by {LSCR} and {SPS}[END_REF] was employed. The volume reduction, when a higher number of iterations is allowed, is quite evident. One can note that the true parameter is contained in the regions. The amount of data values that node #1 receives before being able to compute the confidence regions is, respectively, 3960 and 29700. The flooding algorithm is investigated in Fig. 4 where the 90% confidence region at node #1 is shown. In this case no truncation was applied and the required amount of data was 49540.

The effect of information diffusion truncation, when a constraint on the maximum amount of received data at each node is present, is investigated in the following results. As indicator of the resulting confidence region extension, we considered the diameter of the single box outer approximation of the confidence region, defined as the longest side of the outer parallelepiped that approximates the confidence region. In Fig. 5 consensus and flooding algorithms are compared, again at node #1. As can be observed, the achieved confidence region with flooding is smaller than that obtainable with consensus when the same amount of received data is considered.

As far as the mixed flooding and consensus approach is concerned, Figs. 6 and7 show the diameter of the single box outer approximation as a function of the received number of data at node #1. Fig. 6 was obtained for the topology shown in Fig. 1, while Fig. 7 was obtained with a smaller network of 50 nodes with d = 0.34, not shown here due to lack of space. As can be observed, in both cases there is a range of amount of received data wherein the mixed approach slightly outperforms flooding, which is, however, the best solution when no limitation on data exchanges is present (see the asymptotic behavior in Fig. 6 and7).

V. CONCLUSIONS

In this paper, the problem of efficient distributed evaluation of confidence regions for parameter estimations has been considered. Three different information diffusion strategies, namely flooding, average consensus and a mixed approach, have been compared in terms of confidence region outer approximation shape as a function of the amount of data required for the computation. Effects of truncation for the presented algorithms have been also addressed. Numerical results have been provided both to evaluate the truncation effect and to measure energetic efficiency of the presented algorithms. The main conclusions are that truncation does affect the region shape, but not the confidence level, and that, when the maximum amount of data exchanges is limited, an approach mixing flooding and average consensus is well suited.

  flood = n p + 1 values in the first iteration and D (last) flood = N (n p + 1) values in the last iteration.

Fig. 1 .

 1 Fig. 1. Network of 100 nodes considered in the simulations.

1 Fig. 2 .

 12 Fig.2. Projection on the (p 1 , p 2 ) and (p 2 , p 3 )-planes of the 90% confidence region at node #1 after 4 consensus iterations. Coordinates in the parameter space are denoted as p 1 , p 2 , p 3 .

Fig. 4 .

 4 Fig. 4. Projection on the (p 1 , p 2 ) and (p 2 , p 3 )-planes of the 90% confidence region at node #1 after a complete flooding run.

Fig. 5 .

 5 Fig. 5. Diameter of the single box outer approximation as a function of the amount of received data at node #1. Consensus curve is in blue.

Fig. 6 .Fig. 7 .

 67 Fig.6. Diameter of the single box outer approximation as a function of the amount of received data at node #1. The mixed approach curve is in blue.

  T . White Gaussian measurement noise was considered, with variance σ 2 = 115. The regressors are composed of random equiprobable and independent elements with values in {-1, 1}. No unit of measurement is specified because we do not want to restrict p * to any specific domain. The communication link distance is fixed to d = 0.18 to ensure network connectivity. Any two nodes under this Fig.3. Projection on the (p 1 , p 2 ) and (p 2 , p 3 )-planes of the 90% confidence region at node #1 after 30 consensus iterations.

	0.34	0.46
	0.3
	p 2	p 3	0.4
	0.26
		0.34
	with the true parameter p * = 3 ] T = [0.2, 0.3, 0.4] 0.25 2 , p * 1 , p * [p * 0.15 0.2 0.22 p 1 0.15 0.2 0.25 0.22 0.3 0.34 0.25 0.34 0.4 0.46 0.26 p 2 p 3 p 1	0.3 0.3 p 2 p 2	0.35 0.35