Christophe Chalons
email: christophe.chalons@uvsq.fr

Path-conservative In-cell Discontinuous Reconstruction schemes for non conservative hyperbolic systems

We are interested in the numerical approximation of discontinuous solutions in non conservative hyperbolic systems. We introduce the basics of a new strategy based on in-cell discontinuous reconstructions to deal with this challenging topic, and apply it to a 2x2 non conservative toy model, and a 3x3 gas dynamics system in Lagrangian coordinates. The strategy allows in particular to compute exactly isolated shocks. Numerical evidences are proposed.

Introduction

In this paper, we are interested in the numerical approximation of non conservative hyperbolic systems of the form

{ ∂ t u + A(u)∂ x u = 0, x ∈ R, t ∈ R +,⋆ , u(x, 0) = u 0 (x), (1.1)
where u(x, t) ∈ R p is the unknown and u 0 the initial data, supplemented with an initial condition u(x, 0) = u 0 (x), x ∈ R,

and the validity of an entropy inequality

∂ t U(u) + ∂ x F(u) ≤ 0. (1.3)
Here (U, F) is an entropy-entropy flux pair, that is to say T ∇UA = T ∇F with U strictly convex. By non conservative, we mean that A is not a Jacobian matrix. This does not prevent some of the equations of (1.1) from being in conservation form, but we assume that they are not all conservation laws. It turns out that the theoretical and numerical study of such systems is a very difficult task as briefly recalled now. Theoretical aspects. Let us first review the main theoretical aspects of non conservative hyperbolic systems. Generally speaking, hyperbolic systems develop discontinuous solutions for large times (see [START_REF] Lax | Hyperbolic systems of conservation laws and the mathematical theory of shock waves[END_REF]), so that solutions in a weak sense are considered. When the model is made of conservation laws, solutions are usually defined in the sense of distributions and, under the validity of an entropy inequality, existence and uniqueness results are proved for initial data close to a constant state (see for instance Liu [START_REF] Liu | The riemann problem for general systems of conservation laws[END_REF], [START_REF] Liu | The deterministic version of the glimm scheme[END_REF], Glimm [START_REF] Glimm | Solutions in the large time for nonlinear hyperbolic systems of equations[END_REF], Lax [START_REF] Lax | Hyperbolic systems of conservation laws[END_REF], [START_REF] Lax | Hyperbolic systems of conservation laws and the mathematical theory of shock waves[END_REF], or LeFloch [START_REF] Lefloch | Hyperbolic systems of conservation laws : the theory of classical and nonclassical shock waves[END_REF] for a review and extensions). In the case of a non conservative system made of one or several non conservation laws, the distribution theory does not apply anymore. Dal Maso, LeFloch and Murat proposed in [START_REF] Maso | Definition and weak stability of a non conservative product[END_REF] a definition of the non conservative product A(u)∂ x u which extends the notion of weak solution of conservation laws. More precisely, they introduce the paths theory to define A(u)∂ x u thanks to a family of paths ϕ : [0, 1] × Ω × Ω → Ω satisfying the consistency property ϕ(0, u, u) = u, ϕ(1, u, u) = u, for all (u, u) ∈ Ω × Ω.

Under specific assumptions given in [START_REF] Maso | Definition and weak stability of a non conservative product[END_REF], the non conservative product A(u)∂ x u at a given point x 0 separating two constant states u 0 = (u 0 , v 0 ) and u 1 = (u 1 , v 1 ) is defined by

[A(u)∂ x u] ϕ = ∫ 1 0
A(ϕ(s, u 0 , u 1 )) ∂ϕ ∂s (s, u 0 , u 1 )ds δ x0 .

(1.4) Said differently, the shock is admissible provided that the generalized Rankine-Hugoniot relation

-σ(u 1 -u 0 ) + ∫ 1 0
A(ϕ(s, u 0 , u 1 )) ∂ϕ ∂s (s, u 0 , u 1 )ds = 0 (1.5) holds true, where σ denotes the speed of propagation of the shock. Similarly to the conservative setting, this definition leads to existence and uniqueness results of weak solutions to (1)-( 2) but a first difficutly is clearly to define the relevant path according to the physics of the model under consideration.

Numerical aspects (in brief ). The numerical approximation of discontinuous solutions in non conservative systems is a very difficult task. The main reasons are the deep sensitiveness of the standard methods with respect to the choice of the path and the usual discretisation parameters, see for instance [START_REF] Hou | Why nonconservative schemes converge to wrong solutions: Error analysis[END_REF], [START_REF] Castro | Why many theories of shock waves are necessary: Convergence error in formally path-consistent schemes[END_REF], [START_REF] Berthon | Why many theories of shock waves are necessary: kinetic relations for non-conservative systems[END_REF], [START_REF] Lefloch | Numerical methods with controlled dissipation for small-scale dependent shocks[END_REF] and the references therein, as well as the lack of a Lax-Wendroff type convergence result. In particular, it is not guaranteed that the converged solution satisfies the path theoretical requirement (1.5). The literature is large on the topic but the proposed schemes are often not satisfying in the sense that either they work only for some very particular systems or small amplitude shocks, or they involve some random sampling techniques which are difficult to extend in several space dimensions. Without any attempt to be exhaustive, we refer for instance the reader to [START_REF] Berthon | Nonlinear projection methods for multi-entropies navier-stokes systems[END_REF], [START_REF] Berthon | Schéma nonlinéaire pour l'aproximation numérique d'un système hyperbolique non conservatif[END_REF], [START_REF] Chalons | Navier-stokes equations with several independent pressure laws and explicit predictor-corrector schemes[END_REF], [START_REF] Audebert | Hybrid godunov-glimm method for a nonconservative hyperbolic system with kinetic relations[END_REF], [START_REF] Castro | Entropy conservative and entropy stable schemes for nonconservative hyperbolic systems[END_REF], [START_REF] Fjordholm | Accurate numerical discretizations of non-conservative hyperbolic systems[END_REF], [START_REF] Parès | Numerical methods for non-conservative hyperbolic systems: a theoretical framework[END_REF], [START_REF] Chalons | A new comment on the computation of non conservative products using roe-type path conservative schemes[END_REF] and the references therein where different models and numerical approaches have been considered. Among these methods, the most recent and complete theory is probably the so-called path-conservative schemes theory, developed by C. Pares [START_REF] Parès | Numerical methods for non-conservative hyperbolic systems: a theoretical framework[END_REF] and collaborators. However, it was proved in [START_REF] Abgrall | A comment on the computation of non-conservative products[END_REF], [START_REF] Castro | Why many theories of shock waves are necessary: Convergence error in formally path-consistent schemes[END_REF] that the consistency defnition provided by the pathconservative formalism is not always enough to ensure the convergence to the expected solution. This is especially true in the case of small-scale dependent solutions of interest in the present paper, again because of a lack of control of the numerical diffusion. Nevertheless, we will see that when combined with a suitable in-cell discontinuous reconstruction strategy, the path-conservative formalism allows to control the numerical diffusion in numerical shocks.

General context. The present contribution follows a series of recent works on this topic, and more precisely the two comments on the computation of non conservative products recently given in [START_REF] Abgrall | A comment on the computation of non-conservative products[END_REF] and [START_REF] Chalons | A new comment on the computation of non conservative products using roe-type path conservative schemes[END_REF]. In a few words, the authors consider in [START_REF] Abgrall | A comment on the computation of non-conservative products[END_REF] the gas dynamics equations in Lagrangian coordinates and show numerically that path-conservative schemes are not convergent to the correct solution when applied to a non conservative version of these equations. This fact was explained theoretically in [START_REF] Castro | Why many theories of shock waves are necessary: Convergence error in formally path-consistent schemes[END_REF]. In [START_REF] Chalons | A new comment on the computation of non conservative products using roe-type path conservative schemes[END_REF], the authors consider the same set of equations and show how to slightly modify the usual path-conservative schemes to compute correctly the solutions of this non conservative formulation. The proposed modification is based on a new averaging procedure of the path-conservative schemes and relies on both the introduction of modified averaging cells and a random sampling at each time step. The numerical results are really convincing and a convergence result is proved for isolated shocks. This shows that if the averaging procedure is dealt with care, then the path-conservative approximate Riemann solvers can be a powerful tool for the purpose of computing non conservative shocks. This was actually the main message of [START_REF] Chalons | A new comment on the computation of non conservative products using roe-type path conservative schemes[END_REF]. However and as already stated above, the averaging procedure proposed in [START_REF] Chalons | A new comment on the computation of non conservative products using roe-type path conservative schemes[END_REF] relies on a random sampling and it is well-known from the work by Collela [18] that the computation a shocks with Glimm's random choice type methods is difficult to extend in several space dimensions. Therefore it could be a strong limitation for future works.

Objective of the paper. The aim of the present contribution is to propose a new averaging strategy based on in-cell discontinuous reconstruction in order to get rid of random sampling and modified cells. As we will see, it allows to follow isolated shocks exactly, and provides convergent results to the correct solution for more general initial data. In-cell discontinuous reconstruction techniques used in the present paper were developed by F. Lagoutière and B. Després to reduce the numerical diffusion in the transport of discontinuous solutions of linear and non linear equations, see for instance [START_REF] Després | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF], [START_REF] Lagoutière | A non-dissipative entropic scheme for convex scalar equations via discontinuous cell-reconstruction[END_REF], [START_REF] Lagoutière | Stability of reconstruction schemes for scalar hyperbolic conservations laws[END_REF] before being extended to different setting. In particular, in [START_REF] Boutin | Convergent and conservative schemes for nonclassical solutions based on kinetic relations[END_REF] and [START_REF] Aguillon | Capturing nonclassical shocks in nonlinear elastodynamic with a conservative finite volume scheme[END_REF], the authors define a conservative scheme which is based on in-cell discontinuous reconstructions of non classical shocks for approximating the solutions of non convex scalar conservation laws and non genuinely non linear systems of conservation laws. Again, the striking feature of the strategy is to allow for a perfect control of the numerical diffusion associated with the non classical discontinuities. More precisely, it allows for the exact computation of such isolated simple waves. In [START_REF] Chalons | General constrained conservation laws. Application to pedestrian flow modeling[END_REF], [START_REF] Chalons | A conservative scheme for non-classical solutions to a strongly coupled PDE-ODE problem[END_REF] and [START_REF] Villa | Moving bottlenecks for the Aw-Rascle-Zhang traffic flow model[END_REF], the authors succeeded in extending this approach based on in-cell reconstructions to constrained (scalar or systems of) conservation laws in traffic modeling. In the present contribution, we aim at considering a first step towards the extension of in-cell discontinuous reconstructions towards non conservative systems. Despite the present contribution is only the very beginning in the development of this strategy applied to non conservative systems, we believe that it might be considered as a relevant alternative to numerical methods involving random sampling, which are so far the only ones for which convergence results can be proved.

Outline of the paper. The outline of the paper is as follows. In section 2, we consider a non conservative toy model and show how the in-cell discontinuous reconstruction strategy can be used to define a relevant projection onto the set of piecewise constant solutions at each time step and therefore to properly compute the shock discontinuities. Note that for this toy model, the exact Riemann solver will be used to define the in-cell reconstructions. At last, Section 3 considers the non conservative gas dynamics equations in Lagrangian coordinates and show how the discontinuous reconstruction strategy can be combined with the use of an approximate Riemann solver while keeping the same accuracy in the shocks computations. The last section gives the main conclusions and perspectives of this work.

Application to a non conservative toy model

In this section, we are interested in the numerical approximation of the weak solutions of the following non conservative system of two partial differential equations :

     ∂ t u + ∂ x u 2 2 + u∂ x v = 0, ∂ t v + ∂ x v 2 2 + v∂ x u = 0, (x, t) ∈ R × R + , ( 2.1)
where u = (u, v) t belongs to the state space Ω = {u ∈ R 2 , u + v > 0}. This system can be given the condensed form (1.1) where the non Jacobian matrix A is defined by

A(u) = ( u u v v ) . (2.2)
Such a model, which consists of two coupled Burgers equations, is probably the simplest example of non conservative model. It has already been studied for instance in [START_REF] Berthon | Schéma nonlinéaire pour l'aproximation numérique d'un système hyperbolique non conservatif[END_REF] and can be understood as a simplified two-fluid model where u and v denote the velocity of each fluid.

Let us state useful properties of the model (see again [START_REF] Berthon | Schéma nonlinéaire pour l'aproximation numérique d'un système hyperbolique non conservatif[END_REF], or [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] for the basic definitions), the proof of which is left to the reader.

Lemma 2.1. System (1.1) is strictly hyperbolic over Ω with eigenvalues

λ 1 (u) = 0 < λ 2 (u) = u + v,
and eigenvectors

r 1 (u) = (1, -1) t , r 2 (u) = (u, v) t .
The first characteristic field is linearly degenerate and the second characteristic field is genuinely nonlinear. Moreover, the Riemann invariants are respectively given by

I 1 (u) = u + v, I 2 (u) = u/v.
Remark. We have implicitly assumed v ̸ = 0 in the definition of I 2 . In the case v = 0 and for u ∈ Ω, the Riemann invariant is given by I 2 (u) = v/u. Lemma 2.2. Smooth solutions of (1.1) obey the following additional conservation laws

∂ t (u + v) + ∂ x (u + v) 2 2 = 0, ∂ t ( v u + v ) = 0. (2.3)
More generally, for any convex function f from R to R, smooth solutions of (1.1) satisfy

∂ t f (u + v) + ∂ x ( ∫ u+v sf ′ (s)ds ) = 0. (2.4)
In other words, the mapping (u, v) → f (u + v) is an entropy of (1.1).

In the forthcoming developments, the initial-value problem (1.1)-(1.2) is supplemented with the validity of the entropy inequality

∂ t f (u + v) + ∂ x ( ∫ u+v sf ′ (s)ds ) ≤ 0 (2.5)
in the usual distributional sense and for any convex function f from R to R. As discussed in the introduction, such an entropy inequality is sufficient to prove existence and uniqueness of solutions close to a constant states when the system is conservative. Here, we are clearly in a non conservative setting and according to the path theory of [START_REF] Maso | Definition and weak stability of a non conservative product[END_REF], an additional information encompassed in the so-called paths is needed for the problem to be well-posed.

A first example of family of paths. Following Volpert [START_REF] Volpert | The space bv and quasilinear equations[END_REF], one can choose for ϕ the straight lines family given by ϕ(s, u, u

) = u + s(u -u), ∀ u, v ∈ Ω, ∀ s ∈ [0, 1].
In this case, (1.5) writes

     -σ(u 1 -u 0 ) + (u 1 + u 0 ) 2
( (u 1 + v 1 ) -(u 0 + v 0 ) ) = 0, -σ(v 1 -v 0 ) + (v 1 + v 0 ) 2
( (u 1 + v 1 ) -(u 0 + v 0 ) ) = 0.
(2.6)

Let us assume that u 0 ̸ = u 1 . Then, if u 0 +v 0 ̸ = u 1 +v 1 , it is easy to check that (2.6) can be equivalently written

{ v 1 u 0 = v 0 u 1 , σ = 1 2 ( (u 0 + v 0 ) + (u 1 + v 1 ) ) , ( 2.7)
while in the case u 0 +v 0 = u 1 +v 1 , meaning that the value of the first Riemann invariant I 1 is the same, we get σ = 0 and the discontinuity is a contact discontinuity associated with the first characteristic field.

Note that the condition u 0 + v 0 > 0 implies existence and uniqueness of u 1 = (u 1 , v 1 ) satisfying (2.7) for any given u 0 = (u 0 , v 0 ) and σ. Conversely, condition u 1 + v 1 > 0 implies existence and uniqueness of u 0 = (u 0 , v 0 ) satisfying (2.7) for any given u 1 = (u 1 , v 1 ) and σ. In the following, we will use the notation

{ u 0 = φ(u 1 , σ), σ = 1 2 ( (u 0 + v 0 ) + (u 1 + v 1 ) ) ,
where φ will be called a kinetic function.

A second example of family of paths. Following LeFloch [START_REF] Lefloch | Shock waves for nonlinear hyperbolic systems in nonconservative form[END_REF] and Sainsaulieu [START_REF] Sainsaulieu | Travelling waves solutions of convection-diffusion systems whose convection terms are weakly nonconservative[END_REF], ϕ can also be implicitly defined by adding a second-order diffusion tensor to (2.1) :

     ∂ t u + ∂ x u 2 2 + u∂ x v = ε 1 ∂ xx (u + v), ε 1 > 0, ∂ t v + ∂ x v 2 2 + v∂ x u = ε 2 ∂ xx (u + v), ε 2 > 0.
(2.8)

In this case, a shock discontinuity (σ, u 0 , u 1 ) is said to be admissible if there exists a travelling wave solution of (2.8) such that :

u(x, t) = u(ξ), ξ = x -σt, lim ξ→-∞ u(ξ) = u 0 , lim ξ→+∞ u(ξ) = u 1 .
(2.9)

It is shown in LeFloch [START_REF] Lefloch | Shock waves for nonlinear hyperbolic systems in nonconservative form[END_REF] how to derive a family of paths consistant with this definition. Berthon [START_REF] Berthon | Schéma nonlinéaire pour l'aproximation numérique d'un système hyperbolique non conservatif[END_REF] used this definition and showed for system (2.1) that for any u 0 in Ω and σ in ](u 0 + v 0 )/2, (u 0 + v 0 )[, there exists a unique state u 1 ̸ = u 0 in Ω and a unique travelling wave solution (up to a translation) satisfying (2.9) and such that the generalized Rankine-Hugoniot conditions (1.5) write

     v 1 = ε 2 ε 1 + ε 2 ( 2σ -(u 0 + v 0 ) ) + ε 1 v 0 -ε 2 u 0 ε 1 + ε 2 e ( 2-2(u0+v0)/σ ) , σ = 1 2
( (u 0 + v 0 ) + (u 1 + v 1 )
) ,

(2.10)

or equivalently      v 0 = ε 2 ε 1 + ε 2 ( 2σ -(u 1 + v 1 ) ) + ε 1 v 1 -ε 2 u 1 ε 1 + ε 2 e ( 2-2(u1+v1)/σ ) , σ = 1 2
( (u 0 + v 0 ) + (u 1 + v 1 )
) .

(2.11)

We note in particular that the exit state u 1 actually depends on the shape of the diffusion tensor, and more precisely on the ration ε 2 /ε 1 . This is a characteristic of non conservative systems as illustrated in various contributions on this subject, see for instance Raviart and Sainsaulieu [START_REF] Raviart | A nonconservative hyperbolic system modeling spray dynamics. part 1 : solution of the riemann problem[END_REF], Sainsaulieu [START_REF] Sainsaulieu | Travelling waves solutions of convection-diffusion systems whose convection terms are weakly nonconservative[END_REF], Berthon and Coquel [START_REF] Berthon | Shock layers for turbulence models[END_REF], [START_REF] Berthon | Multiple solutions for compressible turbulent flow models[END_REF], Chalons and Coquel [START_REF] Chalons | The riemann problem for the multi-pressure euler system[END_REF], [START_REF] Chalons | Navier-stokes equations with several independent pressure laws and explicit predictor-corrector schemes[END_REF]... We alsor refer to Berthon, Coquel and LeFloch [START_REF] Berthon | Why many theories of shock waves are necessary: kinetic relations for non-conservative systems[END_REF]. Here again, we will use the notation

{ u 0 = φ(u 1 , σ), σ = 1 2 ( (u 0 + v 0 ) + (u 1 + v 1 )
) .

A path-conservative in-cell discontinuous numerical scheme

Let us now turn to the numerical approximation of the solutions of our toy model. We first introduce some notations and briefly recall the usual Godunov scheme. As we will see, this scheme fails in approximating the correct shock solutions defined by a family of paths ϕ, but it will be useful for approximating the smooth parts of the solutions, in particular the rarefaction waves. We thus motivate and describe the proposed in-cell discontinuous reconstruction strategy which allows in particular to compute exactly any isolated admissible shock. This property is the key property to explain the success of the approach for general initial data. We introduce a constant space step ∆x and constant time step ∆t and we set ν = ∆t/∆x. The mesh interfaces are defined by x j+1/2 = j∆x for j ∈ Z and the intermediate times by

t n = n∆t for n ∈ N.
As usual in the finite volume framework, we seek at each time t n for an approximation u n j of the solution in the interval [x j-1/2 , x j+1/2 ), j ∈ Z. Therefore, a piecewise constant approximate solution x → u ν (x, t n ) of the solution u is given by

u ν (x, t n ) = u n j for all x ∈ C j = [x j-1/2 ; x j+1/2 ), j ∈ Z, n ∈ N.
When n = 0, we set

u 0 j = 1 ∆x ∫ x j+1/2
x j-1/2 u 0 (x)dx, for all j ∈ Z.

Failure of the classical Godunov scheme

The classical Godunov scheme is composed of two steps : a first step in which the solution evolves in time according to the PDE model under consideration, and a second step of projection onto piecewise constant functions.

Step 1 : Evolution in time In this first step, one solves the following Cauchy problem

{ ∂ t u(x, t) + A(u(x, t))∂ x u(x, t) = 0, x ∈ R, u(x, 0) = u ν (x, t n ), (2.12)
with the given family of paths for times t

∈ [0, ∆t]. Recall that x → u ν (x, t n ) is piecewise constant. Then, under the usual CFL restriction ∆t ∆x max{|λ i (u)|, i = 1, 2} ≤ 1 2 , ( 2.13)
for all the u under consideration, the solution of (2.12) is known by gluing together the solutions of the Riemann problems set at each interface :

u(x, t) = u r ( x -x j+1/2 t ; u n j , u n j+1 ) for all (x, t) ∈ [x j , x j+1 ] × [0, ∆t], (2.14)
where

(x, t) → u r ( x t ; u L , u R ) denotes the self-similar solution of the Riemann problem    ∂ t u(x, t) + A(u(x, t))∂ x u(x, t) = 0, x ∈ R, t ∈ R +,⋆ u(x, 0) = { u L if x < 0, u R if x > 0,
given appendix A whatever u L and u R are in the phase space Ω. Recall that this solution actually depends on the family of paths under consideration.

Step 2 : Projection In order to get a piecewise constant approximate solution on each cell C j at time t n+1 , the solution x → u(x, ∆t) given by (2.14) is simply averaged on C j , as expressed by the following update formula:

u n+1 j = 1 ∆x ∫ x j+1/2 x j-1/2 u(x, ∆t)dt, j ∈ Z.
(2.15)

In the following, it will be useful to write (2.15) equivalently as

u n+1 j = 1 2 ( u n+1 j,L + u n+1 j,R ) , j ∈ Z, (2.16)
with

u n+1 j,L = 2 ∆x ∫ xj x j-1/2 u r ( x -x j-1/2 ∆t ; u n j-1 , u n j )dx (2.17)
and

u n+1 j,R = 2 ∆x ∫ x j+1/2 xj u r ( x -x j+1/2 ∆t ; u n j , u n j+1 )dx. ( 2

.18)

As illustrated on Figure 1 obtained with initial data

u 0 (x) = (u, v) 0 (x) = { (6, 5) if x < 0.5, (0.7, 0.3) if x > 0.5, (2.19)
and the second family of path with ϵ 2 = ϵ 1 , the numerical results provided by this scheme are not satisfactory when a shock is present in the solution in the sense that the intermediate state is different from the exact one. On the contrary, if we consider for instance ϵ 2 = 10ϵ 1 and As clearly explained for instance in [START_REF] Berthon | Shock layers for turbulence models[END_REF], [START_REF] Berthon | Multiple solutions for compressible turbulent flow models[END_REF], [START_REF] Berthon | Schéma nonlinéaire pour l'aproximation numérique d'un système hyperbolique non conservatif[END_REF], [START_REF] Chalons | The riemann problem for the multi-pressure euler system[END_REF], [START_REF] Chalons | Navier-stokes equations with several independent pressure laws and explicit predictor-corrector schemes[END_REF],the main reason of this failure is the excessive numerical diffusion of the Godunov scheme across the shocks, which disagrees with the underlying regularization operator at the discrete level. In other words, the numerical diffusion plays a crucial role and must be controlled to make the approximate and exact solutions coincide. If the numerical diffusion does not exactly mimic the action of the regularization operator, the numerical solutions disagree with the exact solutions. This is observed with the usual Godunov scheme but also with any standard finite volume scheme. The sensitiveness with respect to the numerical diffusion is typical of non conservative systems, but also appear in conservative systems (when the matrix A is the Jacobian matrix of a flux function), when the system is hyperbolic but has at least one characteristic field that is neither genuinely nonlinear, nor linearly degenerate, or when it is not hyperbolic but mixed hyperbolic-elliptic. Such systems need also to be closed by a kinetic relation, which is similar to the previous notion of path, and can give rise to the so-called non classical shock waves, see for instance [START_REF] Lefloch | Hyperbolic systems of conservation laws : the theory of classical and nonclassical shock waves[END_REF]. From a numerical point of view, similar issues to those already discussed come out and approximating non classical shocks is challenging because of the dependence on the underlying diffusion mechanisms. Again, standard techniques are useless and a deeper analysis shows that the failure can be related to the (un)control of the numerical diffusion.

u 0 (x) = (u, v) 0 (x) = { (1, 2) if x < 0.5, (5, 1) if x > 0.5, (2.20)
In order to overcome this difficulty, a new numerical approach was first proposed in [START_REF] Boutin | Convergent and conservative schemes for nonclassical solutions based on kinetic relations[END_REF] (see also [START_REF] Aguillon | Capturing nonclassical shocks in nonlinear elastodynamic with a conservative finite volume scheme[END_REF]) compute non classical solutions to scalar conservation laws. The proposed scheme is fully conservative on fixed meshes and has the property of exactly capturing isolated non classical shocks. For such isolated discontinuities, the underlying numerical diffusion thus reduces at the minimum, namely at one point, unlike standard finite difference schemes. The method is based on an in-cell discontinuous reconstruction technique performed in each computational cell that may contain a non classical shock. The next section proposes to extend this approach to the present setting of a non conservative system in order to properly compute the underlying (in some sense non classical) shocks on a fixed mesh and satisfy at the same time the family of path.

In-cell discontinuous reconstruction

Overview of the strategy. In the previous sections, it was shown that the Godunov scheme is not a good candidate when shocks are present in the solution, but it works correctly when the solution is smooth. Therefore, we will first propose to keep on using the Godunov scheme "far away" from shock discontinuities. On the contrary, in the vicinity of shock discontinuities, we will follow the same approach as in [START_REF] Boutin | Convergent and conservative schemes for nonclassical solutions based on kinetic relations[END_REF] which consists in adding details in the piecewise constant representation of the approximate solution on each cell C j . More precisely, we will reconstruct discontinuities in the relevant cells C j and use them to define u n+1 j instead of simply using the constant values u n j-1 , u n j and u n j+1 like in the Godunov scheme. As we will see hereafter, such an approach will allow to exactly compute isolated shock discontinuities in the sense that for such solutions u n j will equal the average of the exact solution on the cell C j . The corresponding numerical discontinuity will then be diffused on one cell at most. Such a sharp control of the numerical diffusion is at the core of the success of the strategy.

The reconstruction procedure. It is now a matter of defining which cells are to be concerned with the reconstruction procedure as well as the reconstructed discontinuities themselves, but also the strategy to evaluate u n+1 j using the new details provided by the discontinuous reconstructions. Let us consider the cell C j and proceed as follows. Assume that at time t n ,

(u + v) n j-1 > (u + v) n j+1 . (2.21)
According to the Riemann solver, we consider that a shock discontinuity is expected to appear locally around the cell C j and to develop at the next times t > t n . Indeed, such a shock is present in the Riemann solution associated with the inital states u n j-1 and u n j+1 . Hence and with clear notations, we are tempted to introduce in the cell C j the left and right states u j,l = u * (u n j-1 , u n j+1 ) and u n j,r = u n j+1

x n j-1

x n j+1 x n j d n,u u n j-1 u n j+1 u n j u * (u n j-1 , u n j+1 )
Figure 3: Reconstruction of a shock in cell C j . Example of the u component, assuming that the cell j starts at 0 and has length 1; otherwise, just replace d n,u by x j-1/2 + d n,u ∆x. A similar drawing could be done for the v component.

of the shock which is expected to be present in the Riemann solution associated with u n j-1 and u n j+1 . Since we are considering the cell C j , we require that the reconstructed discontinuity between u j,l and u j,r is located inside C j at a position

xu j = x j-1/2 + d n,u j ∆x, (2.22)
for the u component, and

xv j = x j-1/2 + d n,v j ∆x, (2.23)
for the v component, for some d n,u j and d n,v j in [0, 1]. Note indeed that in general, we will consider that the positions of the discontinuities may be different for both components u and v, see Figure 3. Regarding the position of the discontinuities in the cell, it is natural to impose that the reconstruction procedure has to be conservative, which writes

d n,u j u n j,l + (1 -d n,u j )u n j,r = u n j , (2.24)
or equivalently,

d n,u j = u n j,r -u n j u n j,r -u n j,l , (2.25)
for the u component, and

d n,v j v n j,l + (1 -d n,v j )v n j,r = v n j , (2.26)
or equivalently,

d n,v j = v n j,r -v n j v n j,r -v n j,l , (2.27)
for the v component. Clearly, it is possible to reconstruct such discontinuities inside the cell C j provided that

0 ≤ d n,u j = u n j,r -u n j u n j,r -u n j,l ≤ 1, (2.28)
and

0 ≤ d n,v j = v n j,r -v n j v n j,r -v n j,l ≤ 1, (2.29)
which gives two additional conditions for the in-cell reconstruction procedure to make sense.

To conclude the definition of the reconstruction strategy, let us mention that still according to the Riemann solver, it is natural to consider that the speed of propagation σ j,l,r of the reconstructed discontinuity equals σ(u n j-1 , u n j+1 ) for both components u and v, where of course σ(u n j-1 , u n j+1 ) denotes the exact value of the speed of propagation of the shock in the Riemann solution associated with u n j-1 and u n j+1 .

Update formulas. At this stage, the reconstructed discontinuity is completely defined, as well as the reconstruction criterion (2.21), (2.28) and (2.29) for this reconstruction to take place. It thus remains to define the update formula for u n+1 j , as well as the influence of the reconstruction on the update formulas of u n+1 j-1 and u n+1 j+1 . Since σ(u n j-1 , u n j+1 ) > 0, note from now on that for the sake of simplicity and in order to avoid dealing with the interaction of two reconstructed discontinuities in adjacent cells, no reconstruction will be considered in the cell C j+1 .

The cell C j . In this cell, we consider that the system under consideration is completely solved by the reconstructed discontinuity, and thus writes

∂ t u = -[A(u)∂ x u] u ϕ δ x-x u j =σ(u n j-1 ,u n j+1 )t
for the u component, and

∂ t v = -[A(u)∂ x u] v ϕ δ x-x v j =σ(u n j-1 ,u n j+1 )t
for the v component, where with clear notations [A(u)∂ x u] u,v ϕ is given by

-[A(u)∂ x u] u ϕ = -σ(u n j-1 , u n j+1 )(u n j,r -u n j,l )
for the u component, and

-[A(u)∂ x u] v ϕ = -σ(u n j-1 , u n j+1 )(v n j,r -v n j,l )
for the v component. Integrating in space and time, we get for the u component

u n+1 j = u n j - 1 ∆x ∫ x j+1/2 x j-1/2 ∫ t n +∆t t n [A(u)∂ x u] u ϕ δ x-x u j =σ(u n j-1 ,u n j+1 )t , namely u n+1 j = u n j - σ(u n j-1 , u n j+1 )(u n j,r -u n j,l ) ∆x × min(∆t, ∆t u ) (2.30)
where ∆t u is the time needed by the reconstructed discontinuity in u to reach the interface x j+1/2 , that is to say

∆t u = 1 -d n,u j σ(u n j-1 , u n j+1 ) ∆x.
For the v component

v n+1 j = v n j - 1 ∆x ∫ x j+1/2 x j-1/2 ∫ t n +∆t t n [A(u)∂ x u] v ϕ δ x-x v j =σ(u n j-1 ,u n j+1 )t , namely v n+1 j = v n j - σ(u n j-1 , u n j+1 )(v n j,r -v n j,l ) ∆x × min(∆t, ∆t v ) (2.31)
where ∆t v is the time needed by the reconstructed discontinuity in v to reach the interface x j+1/2 ,

∆t v = 1 -d n,v j σ(u n j-1 , u n j+1 ) ∆x.
Formulas (2.30) and (2.31) are equivalent to set u n+1 j = u n j,l if ∆t is greater than the times needed by the reconstructed discontinuities on u and v to reach the interface x j+1/2 . If not, they are equivalent to average the reconstructed discontinuities at their new position in the cell C j after moving at velocity σ j,l,r for a time of length ∆t.

The cell C j+1 . If ∆t is greater than the times needed by the reconstructed discontinuities on u and v to reach the interface x j+1/2 , it is clear that the reconstructed discontinuities are expected to influence the update formulas on the cell C j+1 . However, under the CFL condition (2.13), the reconstructed discontinuities in u and v in the cell C j cannot reach the middle point x j+1 of the cell C j+1 and may thus influence the half interval [x j+1/2 , x j+1 ) only. Since no reconstruction is considered in the cell C j+1 , we consider the usual update formula

u n+1 j+1 = 1 2 ( u n+1 j+1,L + u n+1 j+1,R
) but with (component by component)

u n+1 j+1,L = u n j+1 - 2σ(u n j-1 , u n j+1 )(u n j,r -u n j,l ) ∆x × ( ∆t -min(∆t, ∆t u ) )
in order to take into account the propagation of the reconstructed discontinuities inside the half interval [x j+1/2 , x j+1 ). Note that compared to the usual Godunov scheme, the value of u n+1 j+1,R is unchanged.

To conclude the proposed numerical scheme, let us underline that when no reconstruction takes place in the cells C j-1 and C j , we use the classical Godunov scheme, namely

u n+1 j = 1 2
( u n+1 j,L + u n+1 j,R
) with the definitions (2.17) and (2.18).

Summary. To sum up, the update value of a given cell C j is kept unchanged with respect to the Godunov scheme if no reconstruction takes place in the cells C j-1 and C j , the update value of a given cell C j is completely changed if a reconstruction takes place in the cell C j , and the update value of a given cell C j is partially changed if no reconstruction takes place in the cell C j but a reconstruction takes place in the cell C j-1 . In this case, u n+1 j,L is changed but not u n+1 j,R .

At last, recall that a reconstruction is considered in the cell C j if and only if the criterion (2.21), (2.28) and (2.29) are satisfied and the criterion (2.21), (2.28) and (2.29) adapted to the cell C j-1 are not satisfied.

Following [START_REF] Boutin | Convergent and conservative schemes for nonclassical solutions based on kinetic relations[END_REF], let us prove an important property of the proposed scheme, which explains the very good results obtained in the next section. The result states that isolated shock discontinuities are exactly captured by the scheme and contain no spurious numerical diffusion.

Theorem. Assume that u 0 j = u L if j ≤ 0, u 0 j = u R if j ≥ 1 and that u L and u R are two constant states in the phase space Ω such that they can be joined by an admissible shock discontinuity. In other words, the Riemann solution associated with these left and right states is given by u(x, t) = u L if x < σt and u(x, t) = u R if x > σt where σ is the speed of propagation given by σ = 1 2

( (u L + v L ) + (u R + v R ) )
according to the exact Riemann solver. Then the proposed scheme provides an exact numerical solution on each cell C j in the sense that

u n j = 1 ∆x ∫ x j+1/2 x j-1/2 u(x, t n )dx, j ∈ Z, n ∈ N. (2.32)
In particular, the numerical discontinuity is diffused on one cell at most.

Proof. Let us first notice that there is no relevant reconstruction in the first time iteration since the only cells which could be affected are j = 0 and j = 1 but by conservation we necessarily have

d 0,u 0 = d 0,v 0 = 1 and d 0,u 1 = d 0,v 1 = 0.
In other words, considering a reconstructed discontinuity in these cells gives back the original averaged value. The Godunov scheme is then used during the first step and as an immediate consequence, equality (3.13) is proved for the first iterate by definition of the Godunov scheme. Note that we have in particular

u 1 1 = u R -σ ∆t ∆x (u R -u L )
component by component.

Let us now see what happens in the next time iteration. It is first clear from above that only C 1 is to be concerned with a reconstructed discontinuity between u L and u R . Interestingly, by conservativity the reconstructed discontinuities in u and v are necessarily located at the exact position of the solution, namely at the position x = x j-1/2 + σ∆t. In other words, we reconstruct the exact solution at time t = ∆t. To get the required identity (2.32) for the second iterate, it is sufficient to focus on the two cells C 1 and C 2 (the other ones are trivial) and for instance on the u variable (the v variable can be dealt with in a similar way). Let us first assume that ∆t u ≤ ∆t. The numerical schemes gives

u 2 1 = u 1 1 - σ(u R -u L ) ∆x × ∆t,
that is to say

u 2 1 = u R -σ ∆t ∆x (u R -u L ) - σ(u R -u L ) ∆x × ∆t = u R -σ 2∆t ∆x (u R -u L )
and u 2 2 = u 1 2 = u R , which clearly coincides with the average of the exact solution after two time steps ∆t on the cells C 1 and C 2 . Let us now assume that ∆t u ≥ ∆t so that the exact shock will pass through the interface x 1+1/2 and be located at position

x = x 1+1/2 + σ(∆t -∆t u )
in the cell C 2 . On the other hand, the numerical schemes gives

u 2 1 = u 1 1 - σ(u R -u L ) ∆x × ∆t u ,
that is to say

u 2 1 = u R -σ ∆t ∆x (u R -u L ) - σ(u R -u L ) ∆x × ∆x -σ∆t σ = u L and u 2 2 = 1 2 ( u R - 2σ(u R -u L ) ∆x × (∆t -∆t u ) + u R ) ,
or equivalently

u 2 2 = u R -σ (∆t -∆t u ) ∆x (u R -u L )
which again clearly coincides with the average of the exact solution on the cells C 1 and C 2 after two time steps. And the process is going on in a similar way for the next time iterations, which proves the result.

Numerical experiments

In this section, we illustrate the behavior of the proposed scheme based on in-cell discontinuous reconstructions.

Test 1. In this first test case, we consider an isolated shock associated with the second family of path with ϵ 1 = ϵ 2 and associated with the left and right states of the following initial data,

u 0 (x) = (u, v) 0 (x) = { (u L , v L ) if x < 0.5, (u R , v R ) if x > 0.5, =
{ (0, 1) if x < 0.5, (-0.00670855951629595, 0.50670855951629590) if x > 0.5.

(2.33) The speed of propagation is σ = 3/4. As we can see on Figure 4, and in agreement with our theorem, the numerical solution is exact and contains only one point of numerical diffusion. On Figure 5, we plot the numerical entropy dissipation D n with respect to the time t n and defined by

D n = 1 2 ( ∑ j ∆x(u n+1 j + v n+1 j ) ) 2 - ∑ j ∆x (u n j + v n j ) 2 2 + ∆t ( (u R + v R ) 3 3 - (u L + v L ) 3 3
) (2.34)
where the sum is taken over the mesh cells. We observe that it is nonpositive as expected.

-

u 0 (x) = (u, v) 0 (x) = { (u L , v L ) if x < 0.4 or x > 0.6, (u R , v R ) if otherwise, (2.35)
with (u L , v L ) = (0, 1) and (u R , v R ) = (5, 2), and again the second family of path with ϵ 1 = ϵ 2 . On Figure 7 we compare the numerical quantities u + v obtained with our scheme and the classical Godunov scheme. Recall that this quantity is conserved so that both methods are expected to give the same solution. Again, we clearly see that the new scheme is less diffusive than the Godunov one at the point of discontinuity of the N -wave profile.

Application to the gas dynamics equations in Lagrangian coordinates

In this section, we apply the in-cell reconstruction technique to the following gas dynamics equations in Lagrangian coordinates :

   ∂ t τ -∂ x u = 0, ∂ t u + ∂ x p = 0, ∂ t E + ∂ x pu = 0, (3.1)
where τ > 0 represents the inverse of a density, u is the velocity and p = p(τ, e) > 0 is the pressure.

Here e > 0 denotes the internal energy and satisfies E = e + u 2 /2. For the sake of simplicity, we consider a perfect gas equation of state p(τ, e) = (γ -1)e/τ where γ > 1. Recall that (3.1) is strictly hyperbolic with eigenvalues λ 0 = 0 and λ ± = ±c, c = √ γp/τ , and that the characteristic field associated with λ 0 is linearly degenerate and the ones associated with λ ± are genuinely nonlinear [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]. On the other hand, the admissible solutions of (3.1) are selected by the Lax entropy inequalities, which here are equivalent to σ(τ + -τ -) > 0 where τ + and τ -are the left and right states of the underlying discontinuity, and σ its speed of propagation.

At this stage, (3.1) is written in a classic conservative form which does not raise any difficulty from a numerical point of view since usual Godunov-type schemes can be used, see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] again. However, the following non conservative formulation of (3.1) can be easily obtained

   ∂ t τ -∂ x u = 0, ∂ t u + ∂ x p = 0, ∂ t e + p∂ x u = 0, (3.2)
where only the last equation on the total energy has been replaced with a last equation on the internal energy. Setting u = (τ, u, e), the matrix A(u) is given by In order to define the admissible solutions of (3.2), we consider again the path theory of Dal Maso, LeFloch and Murat. Here, a very simple choice of path is defined for all u 0 and u 1 such that σ(τ 1 -τ 0 ) > 0 in a linear ay with respect to τ , u and p, namely

A(u) =   0 -1 0 ∂ τ p(τ,
   τ (s) = τ 0 + s(τ 1 -τ 0 ), u(s) = u 0 + s(u 1 -u 0 ), p(s) = p 0 + s(p 1 -p 0 ),
for all s ∈ [0, 1]. Actually, it turns out that easy calculations show that the generalized jump relations (1.5) of the path theory boil down to the classic Rankine-Hugoniot relations applied to (3.1), namely

   σ(τ 1 -τ 0 ) + (u 1 -u 0 ) = 0, -σ(u 1 -u 0 ) + (p 1 -p 0 ) = 0, -σ(E 1 -E 0 ) + (p 1 u 1 -p 0 u 0 ) = 0, (3.3) or equivalently      σ(τ 1 -τ 0 ) + (u 1 -u 0 ) = 0, -σ(u 1 -u 0 ) + (p 1 -p 0 ) = 0, -σ(e 1 -e 0 ) + 1 2 (p 1 + p 0 )(u 1 -u 0 ) = 0. (3.4)
In other words and with such a choice of path, both conservative and non conservative formulations (3.1) and (3.2) select the same solutions.

A Roe-type path-conservative approximate Riemann solver

We begin with the definition of a Roe-type path-conservative approximate Riemann solver associated with (3.2) and a given path ϕ. According to [START_REF] Toumi | A weak formulation of roe's approximate riemann solver[END_REF] and [START_REF] Parès | Numerical methods for non-conservative hyperbolic systems: a theoretical framework[END_REF], it is based on a Roe linearization A ϕ such that 1. for all u L and u R , A ϕ (u L , u R ) has 3 distinct eigenvalues, 2. for all u, A ϕ (u, u) = A(u), 3. for all u L and u R ,

A ϕ (u L , u R )(u R -u L ) = ∫ 1 0 A(ϕ(s, u L , u R ))∂ s ϕ(s, u L , u R )ds.
The three properties are satisfied if we set

A ϕ (u L , u R ) = A(u), u = u(u L , u R ) = (τ , u, e) with τ = τ L + τ R 2 , u = u L + u R 2 , e = p τ γ -1 and p = p L + p R 2 ,
see [START_REF] Munz | On godunov-type schemes for lagrangian gas dynamics[END_REF]. The approximate Riemann solution constructed from the Roe linearization is the solution of

   ∂ t u(x, t) + A ϕ (u L , u R )∂ x u(x, t) = 0, u(x, t = 0) = { u L if x < 0, u R if x > 0,
given by

u(x/t; u L , u R ) =        u L if x/t < -σ(u L , u R ), u * L if -σ(u L , u R ) < x/t < 0, u * R if 0 < x/t < σ(u L , u R ), u R if x/t > σ(u L , u R ), (3.5)
where the left and right intermediate states are easily obtained from the left and right eigenvectors l k and r k , k = 1, 2, 3 of A ϕ (u L , u R ), respectively, namely

u * L = (u R , l 1 )r 1 + 3 ∑ k=2 (u L , l k )r k , u * R = 2 ∑ k=1 (u R , l k )r k + (u L , l 3 )r 3 , and σ(u L , u R ) = c(u(u L , u R )) = √ γp/τ (u L , u R ).
For the sake of clarity, it will be useful to have in mind the wave pattern of this solution, which is recalled on the next figure.

-σ(u L , u R ) σ(u L , u R ) 0 u * L u * R u L u R x = 0
Note that once the solution is defined, one can denote by x → ũ(x, t) the piecewise constant approximate solution obtained by glueing together the Roe-type approximate solutions at each interface, that is to say ũ(x, t) = u((x -x j+1/2 )/t; u n j , u n j+1 ) for all (x, t) ∈ [x j , x j+1 ) × [0, ∆t), j ∈ Z, n ∈ N. One can also define a Roe-type path-conservative scheme according to [START_REF] Parès | Numerical methods for non-conservative hyperbolic systems: a theoretical framework[END_REF] as any Godunov-type scheme by averaging the solution on each cell [x j-1/2 , x j+1/2 ), namely

u n+1 j = 1 ∆x ∫ x j+1/2 x j-1/2 ũ(x, ∆t)dx = 1 2 (u n+1 j,L + u n+1 j,R ), (3.6)
with

u n+1 j,L = 2 ∆x ∫ xj x j-1/2 ũ(x, ∆t)dx = 2 ∆x ( σ(u n j-1 , u n j )∆tu * R (u n j-1 , u n j ) + ( ∆x 2 -σ(u n j-1 , u n j )∆t)u n j )
and

u n+1 j,R = 2 ∆x ∫ x j+1/2 xj ũ(x, ∆t)dx = 2 ∆x ( σ(u n j , u n j+1 )∆tu * L (u n j , u n j+1 ) + ( ∆x 2 -σ(u n j , u n j+1 )∆t)u n j ) , under the CFL restriction ∆t max j∈Z |σ(u n j , u n j+1 )| ≤ ∆x 2 . ( 3.7)
In the sequel, we will also use the notation σ n j+1/2 = σ(u n j , u n j+1 ). However, such a Roe-type path conservative scheme fails in computing correctly the discontinuous solutions of our system, see [START_REF] Abgrall | A comment on the computation of non-conservative products[END_REF], and we now aim at applying the in-cell discontinuous reconstruction method instead. As we will see, such a strategy allows to obtain a perfect agreement between the exact and numerical solutions, and even the exact capture of isolated discontinuities.

A path-conservative in-cell discontinuous numerical scheme

The design principle is the same as for the toy model in section 2.1. The main differences here are the following. First, the exact Riemann solver will be replaced with a Roe-type approximate Riemann solver and second, the local Riemann solutions at each interface may contain two discontinuities propagating with velocities having opposite signs, unlike the toy model for which only one discontinuity propagating with a positive speed could occur. Apart from this, the idea is actually the same, namely to keep on using the classic Godunov type scheme given above "far away" from shock discontinuities, and to add details in the piecewise constant representation of the approximate solution in the vicinity of shock discontinuities.

Reconstruction procedure. Let us first define which cells j are to be concerned with the reconstruction procedure. We consider the cell C j and proceed as follows. Assume that at time t n ,

u n j-1 > u n j+1 . (3.8)
According to the entropy condition σ(τ + -τ -) > 0 and the Rankine-Hugoniot relation σ(τ + -τ -) = -(u + -u -) across a shock discontinuity, we consider that a shock discontinuity is expected to appear locally around the cell C j when (3.8) holds true. This is quite natural since such a shock is actually present in the Riemann solution associated with the inital states u n j-1 and u n j+1 and which will develop at the next times t > t n . Hence, we are tempted to introduce in the cell C j a discontinuity given by the Roe-type path-conservative approximate Riemann solver proposed in the previous section. More precisely and with clear notations the left and right states u j,l and u j,r of the reconstructed solution are defined by

u j,l = u n j-1 and u n j,r = u * L (u n j-1 , u n j+1 ) if (τ n j+1 -τ n j-1 ) < 0, and u j,l = u * R (u n j-1 , u n j+1 ) and u n j,r = u n j+1 if (τ n j+1 -τ n j-1 ) > 0.
The speed of propagation σ j,l,r of the reconstructed discontinuity on the cell j is naturally defined by -σ(u n j-1 , u n j+1 ) if (τ n j+1 -τ n j-1 ) < 0 and σ(u n j-1 , u n j+1 ) if (τ n j+1 -τ n j-1 ) > 0, where of course ±σ(u n j-1 , u n j+1 ) refer to the speeds of propagation of the discontinuities in the Roe-type approximate Riemann solver associated with the initial states u n j-1 and u n j+1 . Since we are considering the cell C j , we also require that the reconstructed discontinuity associated with those left and right states is located inside C j at a position

xα j = x j-1/2 + d n,α j ∆x, (3.9)
with α = τ, u, e and for some d n,α j in [0, 1] which may vary with α. In order to define d n,α j , it is natural to impose that the reconstruction procedure is conservative, namely

d n,α j α n j,l + (1 -d n,α j )α n j,r = α n j , (3.10)
or equivalently,

d n,α j = α n j,r -α n j α n j,r -α n j,l , (3.11)
for the α = τ, u, e component. Clearly, it is possible to reconstruct the discontinuities provided that

0 ≤ d n,α j ≤ 1, (3.12)
which gives three additional conditions for the in-cell reconstruction procedure to make sense.

At last, for the sake of simplicity and in order to avoid dealing with the interaction of two reconstructed discontinuities in adjacent cells, no reconstruction will be considered in the cell C j if (3.8) and (3.12) adapted to the cell C j+1 hold true and σ j+1,l,r < 0, while no reconstruction will be considered in the cell C j-1 if (3.8) and (3.12) adapted to the cell C j-1 hold true and σ j-1,l,r > 0.

Remark. In practice, we also impose to the reconstructed states to be admissible in the sense τ > 0 and e > 0, which is not guaranteed by the Roe approximate solver.

Update formulas.

Let us now give the update formulas for u n+1 j , as well as the influence of the in-cell reconstruction on the update formulas for u n+1 j-1 and u n+1 j+1 since σ j,l,r may be positive or negative. We follow exactly the same approach as for the toy model, which leads to the reconstructed discontinuity propagates with a positive speed in the cell C j , and that no reconstruction will be considered in the cell C j-1 if the reconstructed discontinuity propagates with a negative speed in the cell C j .

The case σ j,l,r > 0 and the cell C j . We set

α n+1 j = α n j - σ j,l,r (α n j,r -α n j,l ) ∆x × min(∆t, ∆t α )
where ∆t α is the time needed by the reconstructed discontinuity in α = τ, u, e to reach the interface x j+1/2 , namely

∆t α = 1 -d n,α j σ j,l,r ∆x.
The case σ j,l,r > 0 and the cell C j+1 . Under the CFL condition (2.13), the reconstructed discontinuities in the cell C j cannot reach the middle point x j+1 of the cell C j+1 and may thus influence only the half interval [x j+1/2 , x j+1 ). Since no reconstruction is considered in the cell C j+1 , we consider the usual update formula

u n+1 j+1 = 1 2 ( u n+1 j+1,L + u n+1 j+1,R
) but with (component by component)

α n+1 j+1,L = α n j+1 - 2σ j,l,r (α n j,r -α n j,l ) ∆x × ( ∆t -min(∆t, ∆t α )
) .

Note that compared to the usual Godunov scheme, the value of α n+1 j+1,R will be changed if and only if an in-cell reconstruction takes place in cell C j+2 .

The case σ j,l,r > 0 and the cell C j-1 (and no reconstruction in this cell). Under the CFL condition (2.13), the reconstructed discontinuities in the cell C j cannot influence the cell C j-1 farther than x j-1 . Since we consider the case where no reconstruction is considered in the cell C j-1 , we consider the usual update formula

u n+1 j-1 = 1 2 ( u n+1 j-1,L + u n+1 j-1,R
) but with (component by component)

α n+1 j-1,R = α n j-1 - 2σ j,l,r (α n j-1 -α n j,l ) ∆x ∆t.
The case σ j,l,r < 0 and the cell C j . We follow the same lines as above which leads to the same update formulas for τ , u and e, namely The case σ j,l,r < 0 and the cell C j-1 . Under the CFL condition (2.13) and as before, the reconstructed discontinuities in the cell C j cannot reach the middle point x j-1 of the cell C j-1 and may thus influence the half interval [x j-1 , x j-1/2 ) only. Since no reconstruction is considered in the cell C j-1 , we consider the usual update formula

α n+1 j = α n j - σ j,l,r (α n j,r -α n j,l ) ∆x × min(∆t,
u n+1 j-1 = 1 2 ( u n+1 j-1,L + u n+1 j-1,R
) but with (component by component)

α n+1 j-1,R = α n j-1 - 2σ j,l,r (α n j,r -α n j,l ) ∆x × ( ∆t -min(∆t, ∆t α ) ) .
Note that compared to the usual Godunov scheme, the value of α n+1 j-1,L will be changed if and only if an in-cell reconstruction with positive speed of propagation takes place in cell C j-2 .

The case σ j,l,r < 0 and the cell C j+1 (and no reconstruction in this cell). Under the CFL condition (2.13), the reconstructed discontinuities in the cell C j cannot influence the cell C j+1 farther than x j+1 . Since we consider the case where no reconstruction is considered in the cell C j+1 , we consider the usual update formula

u n+1 j-1 = 1 2 ( u n+1 j-1,L + u n+1 j-1,R
) but with (component by component)

α n+1 j+1,L = α n j,r - 2σ j,l,r (α n j+1 -α n j,r ) ∆x ∆t.
The case with no reconstruction in the cells C j and C j±1 . At last and to conclude the proposed numerical scheme, let us mention that when no reconstruction takes place in the cells C j-1 , C j and C j+1 , we use the classical Godunov-type scheme, namely

u n+1 j = 1 2 ( u n+1 j,L + u n+1 j,R
) .

Similarly to the toy model, one can easily prove that the scheme satisfies by construction the following theorem.

Theorem. Assume that u 0 j = u L if j ≤ 0, u 0 j = u R if j ≥ 1 and that u L and u R are two constant states in the phase space Ω such that they can be joined by an admissible (entropic) shock discontinuity. In other words, the Riemann solution associated with these left and right states is given by u(x, t) = u L if x < σt and u(x, t) = u R if x > σt where σ is the speed of propagation given by

σ = ± √ - p R -p L τ R -τ L .
Then the proposed scheme provides an exact numerical solution on each cell C j in the sense that

u n j = 1 ∆x ∫ x j+1/2 x j-1/2 u(x, t n )dx, j ∈ Z, n ∈ N. (3.13)
In particular, the numerical discontinuity is diffused on one cell at most.

Numerical experiments

We now propose several test cases to illustrate the behavior of the scheme. The adiabatic coefficient is set to γ = 1.4. We compare the solutions with the ones given by the original path-conservative scheme applied to (3.1) or by a classical conservative scheme applied to (3.2). The domain is [0, 1] and the CFL restriction is 0.45. The first two cases are such that exact solutions are either an isolated discontinuity or two shock discontinuities starting from the same right state. The last test case is inspired from the first test case of [START_REF] Abgrall | Numerical approximation of a compressible multiphase system[END_REF] and has a large pressure jump.

Test 1. The first test case is an isolated shock associated with the initial data (τ, u, p) 0 (x) = (2.09836065573770281, 2.3046638387921279, 1.0) if x < 0.5, (8.0, 0.0, 0.1) otherwise.

The speed of propagation is 0.3905124837953326544238 and the final time of the simulation is t = 0.5. We clearly see on Figure 8 that the original path-conservative scheme fails while Figure 9 shows a perfect agreement between our scheme and the exact solution. Recall that our scheme is exact in this case and therefore captures the discontinuity with only one point of numerical diffusion.

Test 2. The second test case is a Riemann problem leading to three waves, namely two shocks and one contact discontinuity, and corresponds to the following initial data, (τ, u, p) 0 (x) = (5.0, 3.323013993227, 0.481481481481) if x < 0.5, (8.0, 0.0, 0.1) otherwise.

The 3-shock is the same as in the previous test case. The density of the first shock goes from 5.0 to 3.0 and its speed of propagation is 0.509175077217. The final time is 0.5. Again, we observe on Figures 10 and 11 that the original path-conservative scheme fails and the new one succeeds and properly computes the shocks without numerical diffusion.

D n = ∑ j ∆x E n+1 j - ∑ j ∆x E n j + ∆t (p R u R -p L u L ) (3.14)
where the sum is taken over the mesh cells. This quantity is clearly zero for a conservative scheme. It is also expected to converge to zero with the mesh size for a convergent non conservative scheme since our choice of path is equivalent to the classical Rankine-Hugoniot relations applied to the conservative system. It is actually the case for our scheme based on in-cell reconstructions. Interestingly, we observe that the energy dissipation oscillates around zero for a given mesh (the amplitude goes to zero with the mesh size). At last, we conclude this section with a more difficult test case taken from [START_REF] Abgrall | Numerical approximation of a compressible multiphase system[END_REF] and with a large pressure jump in the initial data given by (τ, u, p) 0 (x) = (1/1185, 0, 2.0e11) if x < 0.5, (1/1185, 0, 1.0e5) otherwise, where the density, velocity and pressure units are respectively kg/m 3 , m/s, and P a. The final time of simulation is 2e -8 and the mesh is made of 1000 points. We compare on Figure 13 the solution given by our scheme with the one given by the classical path-conservative scheme but applied directly to the conservative variable τ , u and E, so that it approximates correctly the solution in this case since both the system and the scheme are conservative. Here again, we see that our (non conservative) scheme gives similar results and thus is also able to properly approximate the exact solution.

Conclusion and perspectives

We have introduced the basics of the so-called path-conservative in-cell discontinuous reconstruction schemes for the numerical approximation of shock solutions in non conservative systems. By basics, we mean that it has been applied to quite simple systems, namely a toy model and the non conservative gas dynamics equations in Lagrangian coordinates. The first (respectively second) system has one (resp. two) characteristic fields leading to shocks, but in both cases the sign of the corresponding characteristic speed is known a priori. The next steps are to consider systems for which the sign of the characteristic speeds depends on the state value, like for instance the gas dynamics equations in Eulerian coordinates, but also to consider non conservative systems which do not admit an equivalent conservative formulation like the one considered in the present paper. One can think for instance of systems arising in turbulence modeling or geophysical flows. It is also important to note that a key property to make the in-cell discontinuous reconstruction approach successful, and in particular to have the validity of the theorem stating that it gives the exact solution in the case of isolated shocks, lies in the fact that the underlying Riemann solver (exact or of path-conservative Roe type) is able to provide an exact solution in such a case of an isolated shock discontinuity. This therefore emphasizes the need for the development of approximate and entropy-satisfying Riemann solvers which are able to exactly reproduce isolated shocks, which is to be proposed in a forthcoming contribution too. At last, the scheme is first-order accurate in its present form, although it is ∞-accurate for isolated shocks. The extension to higher order of accuracy is also a current investigation. In particular, using the general high-order path-conservative formalism provides a nice opportunity to extend the present approach and design new high-order finite volume solvers that do not introduce any numerical viscosity on the propagation of isolated shocks, and to explore the extension to multidimensional problems.

Acknowledgments. The author is thankful to Frédéric Coquel and Pascal Jaisson for useful discussions on this work. The author is also very grateful to Manuel J. Castro Díaz and Tomás Morales de Luna for pointing out strong connections between the path-conservative methods and the numerical methods proposed here, and suggesting extensions to high-order accuracy and multidimensional setting in that direction.

A Appendix : exact Riemann solver for the non conservative toy model

In this appendix, we briefly give the solution to the Riemann problem (1.1)-(2.5) with initial condition given by u(x, 0) = u 0 (x) =

{ u L if x < 0, u R if x > 0, (A.1)
for two constant states u L and u R in Ω. By Lemma 2.1, this solution is expected to be made of two simple waves, namely a stationary contact discontinuity associated with λ 1 from u L to an intermediate state u ⋆ and a nonlinear wave associated with λ 2 from u ⋆ to u R . The latter is either a shock discontinuity satisfying the generalized Rankine-Hugoniot relations (1.5) and the entropy inequality (2.5) in the sense of distributions, or a rarefaction wave. Let us go further into details.

Contact discontinuities. As is customary, the set C 1 (u L ) of admissible states u ⋆ that can be joined to u L on the right by a contact discontinuity associated with λ 1 is defined thanks to the Riemann invariants. Here we get

C 1 (u L ) = {u ⋆ = (u ⋆ , v ⋆ ) t ∈ Ω, I 1 (u ⋆ ) = I 1 (u L )} or equivalently C 1 (u L ) = {u ⋆ = (u ⋆ , v ⋆ ) t ∈ Ω, u ⋆ + v ⋆ = u L + v L }.
Given u ⋆ in C 1 (u L ), the stationary contact discontinuity solution of (1.1) is then defined by

u(x, t) = { u L if x < 0, u ⋆ if x > 0.
Rarefaction waves. The set R 1 (u R ) of admissible states u ⋆ that can be joined to u R on the left by a rarefaction wave associated with λ 2 is also defined thanks to the Riemann invariants, together with the compatibility condition λ 2 (u ⋆ ) ≤ λ 2 (u R ). More precisely, we have

R 1 (u R ) = {u ⋆ = (u ⋆ , v ⋆ ) t ∈ Ω, I 2 (u ⋆ ) = I 2 (u R ), λ 2 (u ⋆ ) ≤ λ 2 (u R )} or equivalently R 1 (u R ) = {u ⋆ = (u ⋆ , v ⋆ ) t ∈ Ω, v R u ⋆ = u R v ⋆ , u ⋆ + v ⋆ ≤ u R + v R }.

Figure 1 :

 1 Figure 1: u (left) and v (right) -Contact discontinuity followed by a shock wave -Final time t = 0.05 -1000-point mesh

Figure 2 :

 2 Figure 2: u (left) and v (right) -Contact discontinuity followed by a rarefaction wave -Final time t = 0.05 -1000-point mesh

Figure 4 :Figure 5 :

 45 Figure 4: u (left) and v (right) -Isolated shock -Final time t = 0.2 -100-point mesh

  0

Figure 6 :

 6 Figure 6: u (left) and v (right) -Contact discontinuity followed by a shock -Final time t = 0.05 -100-point mesh

Figure 7 :

 7 Figure 7: u + v -Periodic simulation -Final time t = 1 -600-point mesh

  ∆t α ) component by component, where ∆t α is now the time needed by the reconstructed discontinuity in α to reach the interface x j-1/2 , namely ∆t α = d n,α j |σ j,l,r | ∆x.

Figure 8 :

 8 Figure 8: τ -Test 1 -Classical path-conservative scheme -Final time t = 0.5 -300-point mesh

Figure 9 :Figure 10 :

 910 Figure 9: τ (top left), u (top right) and p (bottom) -Test 1 -Our scheme -Final time t = 0.5 -300-point mesh

Figure 11 :

 11 Figure 11: τ (top left), u (top right) and p (bottom) -Test 2 -Our scheme -Final time t = 0.5 -300-point mesh

3 Figure 12 :

 312 Figure 12: Energy dissipation -Test 2 -Our scheme -Final time t = 0.5 -300-point mesh

Figure 13 :

 13 Figure 13: τ (top left), u (top right) and p (bottom) -Test 3 -Final time t = 2e-8 -1000-point mesh

Given u ⋆ in R 1 (u R ), the rarefaction fan solution of (1.1) is then defined by

where we have set ξ = x/t for t > 0 and where u ⋆ (ξ) is defined by

We refer for instance to [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] for more details.

Shock discontinuities.

As motivated above, the set S 2 (u R ) of admissible states u ⋆ in Ω that can be joined to u R on the left by a shock discontinuity propagating at velocity σ is admissible provided that both the generalized Rankine-Hugoniot relations (1.5) and the entropy inequality (2.5) in the distributional sense hold true. More precisely, u ⋆ has to satisfy

where ϕ and f respectively denote a family of paths and any convex function. Adding the two components of the generalized Rankine-Hugoniot relations in (A.2) gives

which in passing does not depend on the family of paths ϕ anymore, and then

Note that we have implicitly assumed that u ⋆ + v ⋆ ̸ = u R + v R in order to deal with a true shock discontinuity and not a contact discontinuity. Then, since

the definition of σ above and the mean value theorem give

for some s in between (u ⋆ + v ⋆ ) and (u R + v R ), that is to say

By convexity of f , it is thus clear that the entropy inequality in (A.2) is equivalent to

The set S 2 (u R ) is then defined by

for a given family of paths. Given u ⋆ in S 2 (u R ), the shock solution of (1.1) is then defined by

The Riemann solution. Glueing together the simple waves associated with λ 1 and λ 2 and for a given family of paths ϕ, we get that the Riemann solution to (2.1)-(2.5)-(A.1) is given as follows :

-

with ξ = x/t and where u ⋆ and u ⋆ (ξ) are respectively defined by

which gives in particular In particular, we still have σ = (u⋆+v⋆)+(u R +v R ) 2 .


Pierre Schneider

Univ. Artois, UR 4027, Centre de Recherche et d'Études Histoire et Sociétés (CREHS), F-62000 Arras, France

Trade exchanges between the Mediterranean and the Indian Ocean (3 rd cent. BC -3 rd cent. AD) The history of maritime trade begins under the rule of Ptolemy II , who initiated Greek seafaring on the Red Sea . Indeed, after Alexander the Great, war elephants had become part of Hellenistic armies. Ptolemy II, who could not obtain Indian war elephants, decided to have them hunted in East Africa and brought back to Egypt. This policy was taken up by his successors. As a consequence, stopovers and hunting bases were founded along the western coast of the Red Sea and the coast of present-day Somaliland, as far as Cape Gardafui, allowing the Greeks to improve their knowledge of East Africa and the Red Sea navigation. From the beginning of the 2nd century BC, however, war elephants were not used any longer. Yet, this did not put an end to Greek presence in the Red Sea, since merchant (Greek : emporoi) took the place from hunters. However, whereas they could sail as far as the Bab el-Mandeb Strait and Aden, they were prevented from going further by the Arabs, until around 125 BC, when a shipwrecked Indian captain was rescued by the coastguard of King Ptolemy VIII. He promised in return to show how to get to India by direct route through the southwest monsoon and how to sail back when the monsoon had reversed. In fact, the Greeks were now able to reach Indian markets. This was a decisive step in the Indian Ocean trade, which increased in volume and changed in nature, with new 'exotic' products being introduced into the Mediterranean. According to Strabo, there were fewer than twenty ships crossing the Bab el-Mandeb Strait before the battle of Actium (30 BC). But, he adds, once the Roman administration was established in Egypt, one hundred and twenty ships crossed it each year. In fact, the annexation of Egypt was a turning point, given that the Romans, who had been demanding Indian Ocean goods for more than a century, now controlled Alexandria and the Red Sea ports. This caused imports to rocket and led to the significant increase in maritime traffic referred to by Strabo. However, it should be borne in mind that the Romans were not the only ones be involved in this trade system : they had to share benefits with other players (Indians, Persians, Palmyrenians, Aksumites, etc.). The sea routes are described in a document dating to the mid-first century AD (entitled Periplus of the Erythraean Sea), whose author describes the main ports of call in the western Indian Ocean, from the Egyptian ports of the Red Sea as far as Cape Delgado (Mozambique) and the Ganges delta. He lists the goods that can be brought there and the cargo that could be loaded . He provides a description of three main maritime routes -but one must keep in mind the Indian Ocean trade also included land and river routes. Some emporoi were content to stay in the Red Sea, loading goods at Adoulis (Zulah, Eritrea) while others, more interested in aromatics, perfumes and spices, headed for Mouza (al-Mukha, Yemen). Other merchants, on the other hand, crossed the straits to make the East African circuit. Leaving Berenice in July, they sailed a coastal trade, selling and buying as they went, until they reached the furthest port of call (Rhapta, Mozambique). Finally, there were emporoi that headed for North-West India (Barygaza (Bharuch)) or South India (e.g., Mouziris, in the Kochi region), with large ships. The outbound journey started in July and India could be reached in October. The ships would leave back around December-January and arrive around March, at which time Mediterranean shipping would resume. This route was the the privilege of those who could benefit from the financial support of the senatorial elite.