Agnes Chan
email: agnes.chan@math.u-bordeaux.fr

Gérard Gallice
email: gerard.gallice@gmail.com

Raphaël Loubère
email: raphael.loubere@math.u-bordeaux.fr

Pierre-Henri Maire
email: pierre-henri.maire@cea.fr

Positivity preserving and entropy consistent approximate Riemann solvers dedicated to the high-order MOOD-based Finite Volume discretization of Lagrangian and Eulerian gas dynamics

Keywords: simple approximate Riemann solver, Positivity preserving and entropy consistent Riemann solver, Lagrangian gas dynamics, Eulerian gas dynamics, High-order Finite Volume discretization

In this paper we propose to revisit the notion of simple Riemann solver both in Lagrangian and Eulerian coordinates following the seminal work of Gallice "Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates" in Numer. Math., 94, 2003. We provide in this work the relation between the Eulerian and Lagrangian forms of systems of conservation laws in 1D. Then an approximate (simple) Lagrangian Riemann solver for gas dynamics is derived based on the notions of positivity preservation and entropy control. Its Eulerian counterpart is further deduced. Next we build the associated 1D first-order accurate cell-centered Lagrangian Godunov-type Finite Volume scheme and show numerically its behaviors on classical test cases. Then using the Lagrangian-Eulerian relationships, we derive and test the Eulerian Godunov-type Finite Volume scheme, which inherits by construction the properties of the Lagrangian solver in terms of positivity preservation and well-defined CFL condition. At last we extend this Eulerian scheme to arbitrary orders of accuracy using a Runge-Kutta time discretization, polynomial reconstruction and an a posteriori MOOD limiting strategy. Numerical tests are carried out to assess the robustness, accuracy, and essentially non-oscillatory properties of the numerical methods.

Introduction

A Riemann Solver (RS), should it be exact or approximated, is one of the classical tool in Computational Fluid Dynamics (CFD). In all generality it is a method for computing the numerical flux at the interface separating two constant physical states by a discontinuity. S.K. Godunov in an enlightening and inspiring article [START_REF] Godunov | Finite difference methods for the computation of discontinuous solutions of the equations of fluid dynamics[END_REF] has developed a fruitful line of reasoning adopting the solution of the Riemann's initial value problem as a key tool to develop his finite volume (FV) numerical method to simulate compressible gas-dynamics equations. This method has become the famous and world-wide celebrated Godunov method/scheme in CFD [START_REF] Godunov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF][START_REF] Godunov | Reminiscences about difference schemes[END_REF]. Later on, some approximate Riemann solvers have been designed or developed by famous researchers following on Godunov's approach, to cite but a few: Roe [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF], Harten, Lax and van Leer (HLL [START_REF] Harten | On upstream Differencing and Godunov-Type schemes for Hyperbolic Conservation Laws[END_REF]), Einfeld (HLLE [START_REF] Einfeldt | On Godunov-Type Methods for Gas Dynamics[END_REF]), Munz (HLLEM [START_REF] Einfeldt | On Godunov-type methods near low densities[END_REF]), Toro (HLLC [START_REF] Toro | Restoration of the contact surface in the Harten-Lax-van Leer Riemann solver[END_REF]), Engquist [START_REF] Engquist | One-sided difference approximations for nonlinear conservation laws[END_REF], Osher ([48]) and others. We refer the interested readers to the book of Toro [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] and the references therein.

Of all the aforementioned solvers, the most well-known is probably the one due to Roe. It is based on a linearisation of the compressible gas dynamics equations. By construction, it unfortunately has the drawback to create nonphysical results, i.e entropy violating discontinuous waves, or to generate the infamous carbuncle phenomenon. Note that the effect of this drawback can be lessened by implementing an appropriate entropy fix. Over the years, the HLL [START_REF] Harten | On upstream Differencing and Godunov-Type schemes for Hyperbolic Conservation Laws[END_REF] and HLLE [START_REF] Einfeldt | On Godunov-Type Methods for Gas Dynamics[END_REF] Riemann solvers also gained considerable popularity thanks to their simplicity. However, they are considered an incomplete Riemann solvers due to the structure of the solution which neglects some intermediate characteristic fields. This also has the downfall of causing excessive numerical dissipation. The modified version of HLL, namely the HLLC [START_REF] Toro | Restoration of the contact surface in the Harten-Lax-van Leer Riemann solver[END_REF] Riemann solver, takes into account the presence of those intermediate waves, i.e contact discontinuities and shear waves, thus making it a complete Riemann solver. However, in its current form it does not include an exact computation of the waves speeds that are simply estimated using uppermost and lowermost bounds, refer for instance to [START_REF] Batten | On the choice of wavespeeds for the HLLC Riemann solver[END_REF].

Nonetheless, in recent years, there have not been many major breakthroughs in one-dimensional Riemann solver design. Many researchers rather worked on developing multi-dimensional ones. Among others, Balsara [START_REF] Balsara | Multidimensional riemann problem with self-similar internal structure part 1 application to hyperbolic conservation laws on structured meshes[END_REF] introduced the Multi-dimensional, Self-similar, strongly-Interacting, Consistent Riemann solver (MuSIC). The base idea is that it uses the one-dimensional Riemann Solver in all directions to constitute a multi-dimensional solver for conservative hyperbolic systems. Other researchers also use the same strategy, relying on some sort of flux splitting algorithm, using a one-dimensional Riemann Solver in all directions independently to achieve multi-dimensional behaviour.

This shows that conducting research on how to improve and add properties to one-dimensional Riemann solvers, like the present work has some value and can be ground work for further multi-dimensional extensions as well. One common feature of all the aforementioned solvers is that they were all developed adopting the Eulerian framework point of view. It is nonetheless possible to adopt a different approach like in the works of Gallice [START_REF] Gallice | Schémas de type Godunov entropiques et positifs préservant les discontinuités de contact[END_REF][START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Non-conservatifs[END_REF][START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF], who tackles the resolution of the Riemann problem from the Lagrangian standpoint. Starting therefrom, a direct estimation of the waves speeds can be derived naturally. Once built in the Lagrangian framework, the Riemann Solver in the Eulerian framework is obtained easily and inherits the properties of the Lagrangian one. A particular property that will be studied is referred to as positivity preservation of state variables, especially the density and internal energy, which is fundamental in establishing a robust numerical scheme. Recurrently, numerical approximations may generate negative states which then leads to instability or code crash. This event is more critical in the Lagrangian framework due to the moving and deforming grid during computation. Besides Gallice, only few authors have been working adopting the Lagrangian standpoint, for instance, [START_REF] Munz | On Godunov-type schemes for Lagrangian gas dynamics[END_REF], the positivity preserving property of the HLL [START_REF] Harten | On upstream Differencing and Godunov-Type schemes for Hyperbolic Conservation Laws[END_REF] and HLLE [START_REF] Einfeldt | On Godunov-Type Methods for Gas Dynamics[END_REF] solvers is assessed in the Lagrangian framework. Cheng and Shu [START_REF] Cheng | Positivity-preserving Lagrangian scheme for multi-material compressible flow[END_REF] also developed positivity-preserving Lagrangian scheme in a direct ALE (Arbitrary-Lagrangian-Eulerian) point of view. In [START_REF] Vilar | Positivity-preserving cell-centered Lagrangian schemes for multi-material compressible flows: From first-order to high-orders. Part I: The one-dimensional case[END_REF], the counterpart of the HLLC solver in term of the Lagrangian mass coordinate is demonstrated to be positivity-preserving under particular definitions of the left and right wave speeds.

To improve the precision of a numerical scheme, a well-known series of papers has established a solid background to get a high-order space and time discretization. The MUSCL scheme elaborated by Van Leer [START_REF] Van Leer | Towards the ultimate conservative difference scheme V: A second order sequel to Godunov's method[END_REF] is a monotone scheme with second-order accuracy in space which achieved popularity quickly thanks to its elegant simplicity. Another study was done by Woodward and Colella [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF] on the piecewise parabolic method (PPM) that gives a third-order accuracy in space. The idea of both methods is to replace the piecewise constant approximation of the Godunov scheme by high-order reconstructed cell states used as building blocks for any Riemann solver. Since then, numerous authors worked on high-order reconstruction strategies. To cite but a few, the ENO [START_REF] Harten | Uniformly high-order accurate non-oscillatory schemes I[END_REF][START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic Conservation Laws[END_REF][START_REF] Abgrall | On essentially non-oscillatory schemes on unstructured meshes: analysis and implementation[END_REF][START_REF] Casper | A finite-volume high-order ENO scheme for two-dimensional hyperbolic systems[END_REF] and WENO [START_REF] Cockburn | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF][START_REF] Titarev | Finite-volume WENO schemes for three-dimensional conservation laws[END_REF][START_REF] Dumbser | High-order unstructured Lagrangian one-step WENO finite volume schemes for non-conservative hyperbolic systems: Applications to compressible multi-phase flows[END_REF] interpolators are widely used thanks to their good properties in the presence of strong discontinuities, as well as the ADER approach [START_REF] Titarev | ADER schemes for three-dimensional nonlinear hyperbolic systems[END_REF][START_REF] Hidalgo | ADER schemes for nonlinear systems of stiff advection-diffusion-reaction equations[END_REF][START_REF] Balsara | Efficient implementation of ADER schemes for euler and magnetohydrodynamical flows on structured meshes speed comparisons with rungekutta methods[END_REF] discontinuous Galerkin [START_REF] Reed | Triangular mesh methods for neutron transport equation[END_REF][START_REF] Cockburn | The Runge-Kutta local projection P1-Discontinuous Galerkin finite element method for scalar conservation laws[END_REF][START_REF] Jiang | On a cell entropy inequality for discontinuous Galerkin methods[END_REF] and the residual distributive scheme [START_REF] Roe | Fluctuations and signals -a framework for numerical evolution problems[END_REF][START_REF] Abgrall | Residual distribution schemes for conservation laws via adaptive quadrature[END_REF][START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF].

In the present paper, following the seminal work of Gallice [START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF], we shall derive a simple approximate Lagrangian Riemann solver which preserves contact discontinuities. Monitoring the wave speeds of the approximate Riemann solver, we shall provide refined positivity and entropy stablity conditions compared to those initially described in [START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF]. Then, following the general methodology proposed in [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Non-conservatifs[END_REF][START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF], we deduce the simple approximate Eulerian Riemann solver from the Lagrangian one. It naturally inherits the positivity and entropy stability properties from its Lagrangian counterpart. The resulting entropic Godunovtype Finite Volume Eulerian scheme is nothing but the one described in [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Non-conservatifs[END_REF][START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF], however it has been obtained employing a more intuitive approach. The Finite Volume schemes in both Lagrangian and Eulerian frameworks are first order accurate in both space and time to begin with. Then, we pursue introducing their high-order accurate extensions. Polynomial reconstruction using a least-squares minimization algorithm is employed to gain high accuracy in space, further completed with Runge-Kutta time integrators [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF]. Robustness of the high-order computations and its ability to inherit the expected properties from the first order scheme are achieved thanks to the MOOD paradigm [START_REF] Clain | A high-order finite volume method for systems of conservation laws -multi-dimensional optimal order detection (MOOD)[END_REF] based on a posteriori limiting and local order decrementing. In this MOOD paradigm, of paramount importance, is the so-called 'parachute' scheme. This is a robust and dissipative scheme which is provably preserving the physical admissible states. The first-order Eulerian scheme based upon the Lagrangian Riemann solver developed in this paper will play this crucial role. Numerical tests will be carried on to illustrate and validate this family of high-order Eulerian numerical schemes. Therefore the aims of this paper are

• recalling the Euler-Lagrangian mapping which allows us to develop a Lagrangian Riemann solver and its Eulerian counterpart, both sharing interesting properties such as positivity preservation;

• developping upon it a first order Godunov type FV scheme which is provably 1st-order accurate, positivity preserving and entropy stable with an explicit and reasonnable CFL like condition;

• use this scheme as the parachute scheme of a high order MOOD extension;

• validate the Lagrangian, Eulerian and high-order Eulerian schemes on demanding test cases.

The rest of this paper is organized as follows. After this introduction, in the second section we present the relation between the Eulerian and Lagrangian systems of conservation laws and their mapping. The third section deals with the Godonuv-type FV scheme in Lagrangian frame. More precisely the notion of simple approximate Riemann solver for Lagrangian gas dynamics is derived, along with some properties and numerical validations. The Eulerian frame version is then proposed by employing the previously designed Lagrangian Riemann solver as a constitutive building brick in the fourth section. Accordingly numerical validation is also provided. The extension of this Eulerian scheme to high-order of accuracy both in space and time is detailed in section 6 relying on classical Runge-Kutta time discretization, polynomial reconstruction and a posteriori MOOD limiting strategy. Again, a test campaign to assess the performances of this new family of high accurate Eulerian schemes is provided. Conclusions and perspectives are drawn in section 7.

Relation between the Eulerian and Lagrangian forms of systems of conservation laws in 1D

This section consists in recalling formally the links between the smooth and discontinuous solutions of the Lagrangian and Eulerian representations of one-dimensional gas dynamics equations. In the sequel, these notions shall be useful to derive approximate Riemann solvers in both representations. The reader interested in the rigorous mathematical treatment of the links between Eulerian and Lagrangian weak solutions of gas dynamics might refer to [START_REF] Wagner | Equivalence of the Euler and Lagrangian Equations of Gas Dynamics for Weak Solutions[END_REF][START_REF] Peng | Euler-Lagrange change of variables in conservation laws[END_REF]. We also quote [START_REF] Gallice | Matrices de Roe pour des lois de conservation générales sous forme eulérienne ou lagrangienne : applications à la dynamique des gaz et à la magnétohydrodynamique[END_REF][START_REF] Gallice | Système d'Euler-Poisson, magnétohydrodynamique et schémas de Roe[END_REF] wherein a general framework relating Eulerian and Lagrangian representations of general systems of one-dimensional conservations law is introduced to derive Roe matrices in a systematic manner.

One-dimensional gas dynamics written under Eulerian representation

In this paragraph, we briefly recall the one-dimensional gas dynamics equation written under Eulerian representation. Namely, the conservation equations of mass, momentum and total energy are viewed in terms of the current (spatial) position, x, of the material particle and x is referred to as the Eulerian coordinate. Denoting the time, t > 0, in this framework the gas dynamics equations are written under the form of the system of conservation laws

∂u ∂t + ∂ ∂x f (u) = 0. (1)
Here, u = u(x, t) and f = f (u(x, t)) are respectively the vector of conservative variables and the flux vector.

More precisely

u = (ρ, ρu, ρe) t , f = (ρu, ρu 2 + p, ρue + pu) t ,
where ρ denotes the mass density, u is the velocity, e is the specific total energy and p the pressure. The thermodynamic closure of this system of conservation laws shall be precised in the sequel. For the moment, we shall focus on the form of this system under Lagrangian representation assuming that we are considering smooth solutions of (1).

Lagrange-Euler mapping

Let X denotes the initial position of the material particle located at x at time t > 0. We represent the motion of this material particle with respect to time introducing the Lagrange-Euler mapping (LEM)

Φ : X → x = Φ(X, t). (2)
By definition, this mapping satisfies Φ(X, 0) = X and its Jacobian, J, reads J(X, t) = ∂Φ ∂X .

Noticing that J(X, 0) = 1, a continuity argument leads us to assume that J(X, t) > 0 for t > 0. This ensures that the LEM is a one-to-one mapping. The partial derivative of Φ with respect to t holding X fixed is nothing but the kinematic velocity

u(X, t) = ∂Φ ∂t | X .
Time differentiating the Jacobian of the LEM and utilizing the definition of the kinematic velocity yields the Geometric Conservation Law (GCL)

∂J ∂t - ∂u ∂X = 0. ( 3
)
Let us point out that the kinematic velocity and the GCL have been derived by means of the LEM and independently of the system of conservation laws [START_REF] Abgrall | On essentially non-oscillatory schemes on unstructured meshes: analysis and implementation[END_REF]. Obviously, the kinematic velocity coincides with the fluid velocity initially introduced in (1). Physical variables can be represented regardless in the Lagrangian form or in the Eulerian one. By an abuse of notation, we shall set

u(X, t) = u(Φ(X, t), t) = u(x, t),
in the sequel of the paper, knowing that x = Φ(X, t). Time differentiating the foregoing identity, holding X fixed, we arrive at

∂u(X, t) ∂t | X = ∂u(x, t) ∂t | x + u ∂u(x, t) ∂x . (4)
This equation expresses nothing but the fact that the lagrangian time derivative (holding X fixed) coincides with the so-called material derivative. Finally, combining identity (4) and the GCL, we arrive at

J -1 ∂(Ju) ∂t | X = ∂u ∂t | x + ∂(uu) ∂x . (5)
This last identity will be of great help to transform the Eulerian system of conservation laws (1) into its Lagrangian counterpart.

One-dimensional gas dynamics written under Lagrangian form

Substituting identity ( 5) into (1) leads formally to

J -1 ∂(Ju) ∂t | X + ∂ ∂x ( f (u) -uu) = 0.
Recalling that J dX = dx, the above equation turns into

∂(Ju) ∂t | X + ∂ ∂X ( f (u) -uu) = 0. (6)
From now on, we omit to specify that the time derivative is taken holding X fixed. Let us point out that the first component of this system, i.e., the mass conservation equation, boils down to the trivial equation ∂(ρJ) ∂t = 0, which after time integration yields

(ρJ)(X, t) = ρ 0 (X), (7)
where ρ 0 (X) > 0 denotes the initial mass density distribution. Thanks to mass conservation, the Jacobian rewrites J = ρ 0 τ, where τ = 1 ρ is the specific volume assuming that ρ > 0. Substituting this expression of the Jacobian into the GCL leads to

ρ 0 ∂τ ∂t - ∂u ∂X = 0. ( 8
)
Finally, gathering the foregoing results we are now in position to write the Lagrangian form of the onedimensional gas dynamics

ρ 0 ∂U ∂t + ∂F ∂X = 0, (9)
where U = (τ, u, e) t and F = (-u, p, pu) t are respectively the vector of conservative variables and the flux vector in the Lagrangian frame. We observe that the second and third components of the flux and the conservative variables written in Lagrangian and Eulerian forms are related by

U = τu, F = f -uu.
Before proceeding any further, let us provide the thermodynamic closure of the gas dynamics equations. First, we recall that the specific total energy is the sum of the specific internal energy, ε, plus the specific kinetic energy, 1 2 u 2 , that is e = ε + 1 2 u 2 . We assume that the specific internal energy is a strictly convex function with respect to the specific volume, τ, and to the specific entropy, η. This means that the Hessian of ε(τ, η) is positive definite

∂ 2 ε ∂τ 2 > 0, ∂ 2 ε ∂η 2 > 0, ∂ 2 ε ∂τ 2 ∂ 2 ε ∂η 2 - ∂ 2 ε ∂τ∂η 2 > 0.
This amounts to say that (τ, ε) → η(τ, ε) is strictly concave, refer to [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]. The thermodynamic closure is completed by expressing the pressure, p, and the temperature, θ, in terms of the partial derivatives of the specific internal energy

p(τ, η) = - ∂ε ∂τ , θ(τ, η) = ∂ε ∂η . ( 10
)
This corresponds to the complete equation of state, refer to [START_REF] Menikoff | The Riemann problem for fluid flows of real materials[END_REF]. We also assume that the temperature is strictly positive, i.e., θ > 0. The convexity assumption allows us to define the isentropic sound speed

a τ = - ∂p ∂τ 1 2 = ∂ 2 ε ∂τ 2 1 2
.

With this thermodynamic closure, one can easily demonstrate that the one-dimensional Lagrangian gas dynamics equations consist of an hyperbolic system characterized by the real eigenvalues (-

a ρ 0 τ , 0, a ρ 0 τ ).
The interesting reader might find the detailed demonstration, for instance, in [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]. Moreover, the complete equation of state implies the fundamental Gibbs relation

θ dη = p dτ + dε. ( 11
)
Observing that dε = deu du, the substitution of the Lagrangian gas dynamics equations into the Gibbs relation leads to the supplementary conservation law

ρ 0 ∂η ∂t = 0. ( 12
)
This latter demonstrates that the entropy is conserved for smooth flows. More precisely, time integrating the above equation yields η(X, t) = η(X, 0). Namely, the entropy of a material particle at time t > 0 is equal to its initial value. On the other hand, for non smooth flows, the physically admissible piecewise continuous solutions have to satisfy the entropy inequality

ρ 0 ∂η ∂t ≥ 0. ( 13
)
The Eulerian counterpart of the Lagrangian entropy conservation equation ( 12) is readily obtained recalling that ρJ = ρ 0 and thanks to identity ( 5)

∂(ρη) ∂t | x + ∂(ρηu) ∂x = 0. ( 14
)
This equation expresses the conservation of entropy for smooth flows in the Eulerian frame, or equivalently, that the entropy is conserved along the trajectory of material particles. For non smooth flows, the Eulerian counterpart of the Lagrangian entropy inequality [START_REF] Cheng | Positivity-preserving Lagrangian scheme for multi-material compressible flow[END_REF] reads

∂(ρη) ∂t | x + ∂(ρηu) ∂x ≥ 0. ( 15
)

Discontinuous solutions in Lagrangian and Eulerian frames

The aim of this paragraph is to recall briefly the Rankine-Hugoniot (RH) condition related to the system of conservation laws under consideration for both Lagrangian and Eulerian representations. The RH condition expresses the admissibility of piecewise discontinuous solutions with respect to the system of conservation laws written in weak (integral) form, refer for instance to [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]. In addition, we are going to exhibit the relationship between the Lagrangian and the Eulerian discontinuity speeds. This relationship will be of great interest, in what follows, for deriving the wavespeeds of the Eulerian approximate Riemann solver in terms of its Lagrangian counterpart.

Lagrangian Rankine-Hugoniot condition

Let us consider the discontinuity curve defined in the (X, t) plane by the equation X = Ξ(t) and the related piecewise discontinuous function

U(X, t) =        U l if X < Ξ(t), U r if X > Ξ(t).
This is a weak solution of the Lagrangian system of conservation laws [START_REF] Bouchut | Semi-discrete entropy satisfying approximate Riemann solvers. The case of the Suliciu relaxation approximation[END_REF] if and only if the Lagrangian RH condition holds true

- dΞ dt ρ 0 U + F = 0.
Here, U = U r -U l is the jump of U across the discontinuity. Thanks to [START_REF] Boscheri | A second-order cell-centered lagrangian ADER-MOOD finite volume scheme on multidimensional unstructured meshes for hydrodynamics[END_REF], Lagrangian mass conservation RH condition reads -dΞ dt ρ 0 = 0. Furthermore, recalling the identity ab = a b + b a , where a = 1 2 (a l + a r ), leads to rewrite the Lagrangian RH condition under the form

-ρ 0 dΞ dt U + F = 0. (16)
Let Σ = ρ 0 dΞ dt be the mass flux swept by the Lagrangian discontinuity. For Σ = 0 it is a contact discontinuity, whereas for Σ > 0 (resp. Σ < 0), it is a rightgoing (resp. leftgoing) shock wave. With this notation, the Lagrangian RH condition reads

-Σ U + F = 0, where U = (τ, u, e) t and F = (-u, p, pu) t . (17)
It is well known that the weak solutions satisfying the RH condition are not unique. Nevertheless the physically admissible solutions are selected supplementing the RH condition [START_REF] Cockburn | The Runge-Kutta local projection P1-Discontinuous Galerkin finite element method for scalar conservation laws[END_REF] by the RH entropy inequality

-Σ η ≥ 0. ( 18
)
This RH entropy inequality simply expresses that the entropy of material particles is increasing through the shock wave consistently with the second law of thermodynamics since a shock wave is an irreversible thermodynamic process.

Eulerian Rankine-Hugoniot condition

We consider the Eulerian representation of the discontinuity curve defined in the (x, t) plane by the equation x = ξ(t) and the piecewise discontinuous function

u(x, t) =        u l if x < ξ(t), u r if x > ξ(t).
This piecewise discontinuous function is a weak solution of the Eulerian system of conservation laws (1) if and only if the Eulerian RH condition holds true

- dξ dt u + f = 0.
Let σ = dξ dt denotes the Eulerian discontinuity speed, then the Eulerian RH condition turns into σ u + f = 0, where u = (ρ, ρu, ρe) t and f = (ρu, ρu 2 + p, ρue + pu) t .

The selection of the physically admissible weak solutions is enforced through the use of the Eulerian RH entropy inequality

-σ ρη + ρηu ≥ 0. ( 20
)

Fundamental relation between Lagrangian and Eulerian discontinuity speeds

A weak solution of the gas dynamics equations is characterized by the RH condition, which can be written under Lagrangian representation [START_REF] Cockburn | The Runge-Kutta local projection P1-Discontinuous Galerkin finite element method for scalar conservation laws[END_REF] as well as under Eulerian representation [START_REF] Dumbser | High-order unstructured Lagrangian one-step WENO finite volume schemes for non-conservative hyperbolic systems: Applications to compressible multi-phase flows[END_REF]. We relate the discontinuity speeds Σ and σ by comparing the first component of ( 17) and ( 19)

-Σ τ -u = 0, (21a) -σ ρ + ρu = 0. ( 21b
)
Before proceeding any further, we introduce some useful notations. Let α and φ be real numbers. We denote respectively (α l , φ l ) and (α r , φ r ) the values taken by (α, φ) on the left and the right side of the discontinuity.

Assuming that the α-weights sum to one, i.e., α l + α r = 1, the weighted averages of φ write

φ = α l φ l + α r φ r , φ = α r φ l + α l φ r .
This allows us to express the jump of a product thanks to the identity

φψ = φ ψ + ψ φ . ( 22
)
Applying this identity to (τ, ρ) and observing that τρ = 1 yields ρ τ = -τ ρ . Substituting this result into (21a) we get ρ u = Στ ρ .

On the other hand, applying identity [START_REF] Engquist | One-sided difference approximations for nonlinear conservation laws[END_REF] to (21b) leads to

-σ ρ + ρ u + u ρ = 0.
Substituting the expression of the velocity jump [START_REF] Fambri | Space-time adaptive ADER-DG schemes for dissipative flows: Compressible Navier-Stokes and resistive MHD equations[END_REF] into the foregoing equation we finally arrive at

(-σ + u + τΣ) ρ = 0.
This equation holds true regardless the mass density jump. Therefore, the term between parentheses is equal to zero, which provides us the fundamental relation between the Eulerian and Lagrangian discontinuity speeds σ = τΣ + u.

Utilizing the definition of the underline average and thanks to (21a) one can easily show that σ satisfies also

σ = τ l Σ + u l = τ r Σ + u r . (25)
Multiplying respectively the first equality by ρ l and the second one by ρ r leads to express the mass flux swept by the Lagrangian discontinuity, Σ, in terms of its Eulerian velocity

Σ = -ρ l (u l -σ) = -ρ r (u r -σ). (26)
Let us point out that the mass flux swept by the Lagrangian discontinuity has exactly the opposite sign to the mass flux crossing the Eulerian discontinuity.

Approximate Riemann solvers for gas dynamics

Here, inspired by the seminal works of Gallice [START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF], we start by designing a family of approximate Riemann solvers for gas dynamics equations written under Lagrangian representation, which depends on two parameters. Then, we describe how to monitor these parameters to ensure good theoretical properties such as positivity of specific volume and internal energy and also an entropy inequality. Finally, we arrive at the description of a family of Eulerian approximate Riemann solvers simply deduced from their Lagrangian counterparts by mimicking the Lagrange Euler mapping at the discrete level. In this framework, the good properties of the Lagrangian approximate Riemann solver are directly transferred to the Eulerian one.

X dX dt = -z l ρ 0 l dX dt = z r ρ 0 r dX dt = 0 t U ⋆ r U ⋆ l 0 -∆X l ∆X r U l U r ∆t

Simple approximate Riemann solver for Lagrangian gas dynamics

The Riemann problem related to Lagrangian gas dynamics reads

ρ 0 ∂U ∂t + ∂F(U) ∂X = 0, U(X, 0) =        U l if X < 0, U r if X ≥ 0. (27)
Let W U l , U r , X t be an approximate solution to the Riemann problem [START_REF] Gallice | Schémas de type Godunov entropiques et positifs préservant les discontinuités de contact[END_REF] defined by

W U l , U r , X t =                                    U l if X t < - z l ρ 0 l , U l if - z l ρ 0 l < X t < 0, U r if 0 < X t < z r ρ 0 r , U r if z r ρ 0 r < X t . (28)
Following Gallice [START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF], W is named simple Riemann solver since it consists of four constant states (U l , U l , U r , U r ) separated by the three discontinuities characterized by the speeds (-

z l ρ 0 l , 0, z r ρ 0 r
). Here, z l and z r are given strictly positive parameters. It is worth noticing that the three discontinuity speeds of W are prescribed to mimic the three eigenvalues of the continuous problem [START_REF] Bouchut | Semi-discrete entropy satisfying approximate Riemann solvers. The case of the Suliciu relaxation approximation[END_REF]. The structure of the simple Riemann solver for Lagrangian gas dynamics is displayed in figure 1. The intermediate states write U s = (τ s , u s , e s ) t and e s = ε s + 1 2 (u s ) 2 for s = l, r. For now on, for any variable or equation with subscript s we assume that it holds for s = l, r.

The consistency of the approximate Riemann solver W with the integral form of the Riemann problem ( 27) is ensured provided that the intermediate states satisfy the equation

F r -F l = -z l (U l -U l ) + z r (U r -U r ), (29)
where F s = F(U s ).This consistency condition results from the integration of (27) over the space-time interval [-∆X l , ∆X r ] × [0, ∆t] displayed in figure 1 and the substitution of the exact solution for the approximate Riemann solver, refer for instance to [START_REF] Harten | On upstream Differencing and Godunov-Type schemes for Hyperbolic Conservation Laws[END_REF].

The consistency of the approximate Riemann solver W with the integral form of the entropy inequality ( 13) is ensured provided that the intermediate states satisfy the inequality

z l (η l -η l ) + z r (η r -η r ) ≥ 0, ( 30
)
where η s = η(U s ) and η s = η(U s ) for s = l, r. This inequality results from the integration of ( 13) over [-∆X l , ∆X r ] × [0, ∆t] and the substitution of the exact solution for the approximate Riemann solver. Let us emphasize that the consistency conditions ( 29) and ( 30) are crucial to ensure that the approximate Riemann solver induces an entropic Godunov-type Finite Volume scheme, refer to [START_REF] Trangenstein | Numerical Solution of Hyperbolic Partial Differential Equations[END_REF]. Now, we are in position to determine the intermediate states U s for s = l, r satisfying the foregoing consistency conditions. There are many solutions since we have 6 unknowns for 3 equations and 1 inequality. We choose to construct the intermediate states by introducing the supplementary unknowns F s named the intermediate fluxes. The intermediate states, U s , and the intermediate fluxes, F s , are linked by the 3 following jump relations written across each discontinuity of W from left to right, refer to figure 1

z l (U l -U l ) + F l -F l = 0, (31a) 0(U r -U l ) + F r -F l = 0, (31b) -z r (U r -U r ) + F r -F r = 0. ( 31c
)
Summing the above equations leads straightforwardly to the consistency condition [START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF]. This is one of the advantages of the introduction of the intermediate fluxes. Combining the 3 equations of system (31) yields the explicit expression of the intermediate flux at the interface

X t = 0 F = 1 2 (F l + F r ) = 1 2 (F l + F r ) - z l 2 (U l -U l ) - z r 2 (U r -U r ). (32)
This shows that the numerical flux is expressed in terms of the states of Riemann solver. Moreover, the numerical viscosity is directly governed by the parameters (z l , z r ). However, the problem is still underdetermined since we have only 9 equations for 12 unknowns, assuming that z l and z r are given parameters. At this point, mimicking the continuous flux, we make the structural assumption that the intermediate fluxes components writes F s = (-u s , p s , (pu) s ) t recalling the intermediate state components U s = (τ s , u s , e s ) t for s = l, r. We arrive at 10 unknowns pointing out that we have used the same velocity unknown for the second component of the intermediate state as well as for the first component of the intermediate flux.

Then, (31b) implies that u l = u r , p l = p r and (pu) l = (pu) r , which leads us to write the intermediate flux F = F l = F r = (-u , p , (pu) ) t . The state vector turns into U s = (τ s , u , e s ) t for s = l, r. Thus, we end up with 7 unknowns for 6 equations given by (31a) and (31c). Finally, assuming that (pu) = p u , we arrive at the following closed systems respectively for the left intermediate states

z l (τ l -τ l ) -(u -u l ) = 0, (33a)
z l (u -u l ) + p -p l = 0, (33b)
z l (e le l ) + p up l u l = 0, (33c) and the right ones

z r (τ r -τ r ) + u -u r = 0, (34a) z r (u -u r ) -(p -p r ) = 0, (34b) z r (e r -e r ) -(p u -p r u r ) = 0. ( 34c
)
We end up with 6 equations for the 6 unknowns (τ l , τ r , u , e l , e r , p ). Next,let us derive expressions of these unknows in terms of the parameters z l and z r . These explicit expressions shall be usefull to exibit explicit conditions on the parameters to enforce positivity of specific volumes, specific internal energies and entropy control on the intermediate states of our approximate Riemann solver. We start by summing (33b) and (34b) and we readily get the explicit expression of u in terms of the parameters z l and z r u = z l u l + z r u r z l + z r -p z l + z r .

Then, substituting this expression into (33a) and (34a) yields the explicit expression of τ l and τ r in terms of parameters z l and z r as

τ l = z l z l + z r       τ l - p z 2 l       + z r z l + z r τ l + u z l , (36a)
τ r = z l z l + z r τ r + u z r + z r z l + z r τ r + p z 2 r . ( 36b
)
It is also interesting to derive the balance of internal energy for the left and right intermediate states. To this end, we compute the balance of kinetic energy dot-multiplying the balance of momentum (33b) (resp.

(34b)) by the centered velocity 1 2 (u + u l ) (resp. 1 2 (u + u r )). Then, subtracting the kinetic energy balance on the left (resp. right) side to the total energy balance (33c) (resp. (34c)) we arrive at

ε s -ε s + p + p s 2 (τ s -τ s ) = 0. (37)
Here, we have used the fact that z l and z r are strictly positive parameters and e s = ε s + 1 2 u .This equation is similar to the well known Hugoniot equation that characterizes the locus of the states resulting from a shock wave. However, this is not exactly the Hugoniot equation since p does not correspond to a thermodynamic pressure, i.e., p p(τ s , η s ). In fact, p is directly expressed by solving the foregoing systems of equation. For instance, substituting (33a) (resp. (34a)) into (33b) (resp. (34b)) leads to the expression of pp s in terms of τ s -τ s and z s p -

p s = -z 2 s (τ s -τ s ). ( 38
)
This equation shows that p is a decreasing function with respect to τ s ,which is a physical admissible behavior. We are now in position to derive explicit conditions on the parameters z l and z r to enforce positivity properties.

Positivity properties and inequality entropy of the Lagrangian simple Riemann solver 3.2.1. Positivity of the specific internal energy

Let us derive an explicit condition on z s to ensure the positivity of ε s .We express ε s in terms of ε s and τ s -τ s substituting p expression [START_REF] Hidalgo | ADER schemes for nonlinear systems of stiff advection-diffusion-reaction equations[END_REF] into the energy balance ( 37)

ε s = ε s -p s (τ s -τ s ) + z s 2 (τ s -τ s ) 2 . ( 39
)
Hence, the specific internal energy is a quadratic convex function with respect to the specific volume jump. It is thus always greater than or equal to its minimum value

ε s ≥ ε s - p 2 s 2z 2 s .
Therefore, the positivity of specific internal energy ε s is ensured provided that parameter z s satisfies the condition

z s ≥ p s √ 2ε s . ( 40
)
This condition is not overly restrictive. Indeed, for a perfect gas equation of state, p = (γ -1)ρε, where γ > 1 is the polytropic index, the foregoing condition boils down to z s ≥ γ -1 2γ ρ s a s . This shows that the parameter z s only needs to be greater than or equal to a fraction of the acoustic impendance to ensure the positivity of specific internal energy. Condition ( 40) is relatively well known and has already been given in [START_REF] Vilar | Positivity-preserving cell-centered Lagrangian schemes for multi-material compressible flows: From first-order to high-orders. Part I: The one-dimensional case[END_REF]. Moreover, if the equation of state is such that the isentrope curves are convex, i.e., the function τ → p(τ, η) is strictly convex with respect to τ, then the following inequality holds true [45]

a 2 τ 2 ≥ p 2 2ε .
Combining this result with condition (40) leads to claim that for a convex equation of state, the approximate Riemann solver preserves the positivity of internal energy provided parameter z s fulfills

z s ≥ ρ s a s . (41)
This last condition simply expresses that the parameters should be greater than or equal to the acoustic impedance. We note in passing that the classical Godunov acoustic solver relying on the particular choice z s = ρ s a s should preserve the positivity of the specific internal energy for any convex equation of state.

Positivity of the specific volume

In this paragraph, we derive explicit conditions on z l , z r to ensure the strict positivity of τ l and τ r . First, by virtue of (36a) and (36b), it is obvious that τ l and τ r result in a convex combination of the terms between parantheses. Thus, the positivity of the specific volume holds true provided that these terms between parentheses are positive. This leads us to claim that the positivity of τ l and τ r is taken for granted if z l and z r are defined by

z l = max       ρ l a l , p + τ l , - u τ l       , and z r = max       ρ r a r , -p + τ r , - u τ r       , (42)
where for x ∈ R, (x) + denotes the positive part of x, i.e., (x) + = 1 2 (x + |x|). Let us point out that the above condition provides us an explicit definition of the parameters z l and z r in terms of the initial left and right states. In addition, these parameters are always greater than or equal to the left and right acoustic impendances. We observe that this condition has been already proposed in [START_REF] Bouchut | Semi-discrete entropy satisfying approximate Riemann solvers. The case of the Suliciu relaxation approximation[END_REF] where approximate Riemann solvers are constructed by means of relaxation schemes. In the sequel of the paper, this approach to monitor the parameters (z l , z r ) will be referred to as method 1.

An alternative approach to enforce the positivity of the specific volumes, (τ l , τ r ), has been proposed in [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Non-conservatifs[END_REF]. It consists in assuming that the parameters, (z l , z r ), are linked by the ratio r = z r z l . Substituting z r = rz l (resp. z l = 1 r z r ) into (36a) (resp. (36b)) leads to the following expressions for τ l (resp. τ r )

τ l = (1 + r)τ l z 2 l + r u z l -p (1 + r)z 2 l , τ r = (1 + r)τ r z 2 r + u z r + r p (1 + r)z 2 r .
We observe that the positivity of τ l (resp. τ r ) is strictly equivalent to the positivity of the quadratic function with respect to z l (resp. z r ) present at the numerator of the above rational expressions. Introducing the discriminants

∆ l = r 2 u 2 + 4(1 + r)τ l p , ∆ r = u 2 -4r(1 + r)τ r p ,
we arrive at the following conditions to ensure the positivity of the specific volumes.

• The positivity of τ l holds true if

∆ l ≤ 0 or if ∆ l > 0 and z l > z + l = 1 2(1 + r)τ l (-r u + ∆ l ), z r > rz + l . (43)
• The positivity of τ r holds true if

∆ r ≤ 0 or if ∆ r > 0 and z r > z + r = 1 2(1 + r)τ r (-u + ∆ r ), z l > 1 r z + r . (44)
For practical applications, the parameter r is defined as the ratio of the acoustic impendances, that is, r = ρ r a r ρ l a l . This method of monitoring the parameters (z l , z r ) to enforce the positivity of the specific volumes will be referred to as method 2 in the sequel of the text. Let us note that this approach might be also employed to enforce the positivity of the specific energies ε l and ε r , the interested reader should refer to [START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF] for more details.

The last approach to enforce the positivity of the specific volumes, (τ l , τ r ), consists in analyzing the signs of the right-hand side of (36a) and (36b) as functions of unknowns z l and z r with parameters u , p , τ l and τ r . The subsequent developments are detailed in Appendix 8. This last method to monitor (z l , z r ) will be referred to as method 3.

Entropy control

We aim at deriving a condition on z s parameters to ensure the entropy inequality, η s ≥ η s , across the discontinuities of the approximate Riemann solver. First, to obtain an expression of the entropy jump η s -η s , we decompose the thermodynamic process (τ s , η s ) → (τ s , η s ) into

• the isentropic process: (τ s , η s ) → (τ s , η s ), followed by

• the isochoric process: (τ s , η s ) → (τ s , η s ).
By virtue of this decomposition, the specific internal energy jump writes

ε(τ s , η s ) -ε(τ s , η s ) = isochoric process ε(τ s , η s ) -ε(τ s , η s ) + isentropic process ε(τ s , η s ) -ε(τ s , η s ) = η s η s ∂ε ∂η (τ s , η) dη + τ s τ s ∂ε ∂τ (τ, η s ) dτ, = η s η s θ(τ s , η) dη - τ s τ s p(τ, η s ) dτ.
The second line at the right-hand side of the above equation results from the definition of the complete equation of state [START_REF] Casper | A finite-volume high-order ENO scheme for two-dimensional hyperbolic systems[END_REF] recalling that the absolute temperature, θ, is strictly positive. Substituting the expression of the internal energy jump (39) into the foregoing equation yields

η s η s θ(τ s , η) dη = τ s τ s p(τ, η s ) dτ -p s (τ s -τ s ) + z s 2 (τ s -τ s ) 2 ,
which is rearranged into Knowing that θ > 0, the entropy inequality, η s -η s ≥ 0 holds true if and only if z s satisfies the condition

η s η s θ(τ s , η) dη = τ s τ s (p(τ, η s ) -p(τ s , η s )) dτ + z s 2 (τ s -τ s ) 2 . ( 45
)
z 2 s ≥ - 2 (τ s -τ s ) 2 τ s τ s (p(τ, η s ) -p(τ s , η s )) dτ, for s = l, r. (46)
The right-hand side of the above inequality is always non negative since ∂p ∂τ < 0 by virtue of the convexity assumption made on the equation of state and thus the above condition is always well defined. Thanks to the mean value theorem, it is also possible to reformulate condition (46) into

z 2 s ≥ a 2 (τ s , η s ) τ2 s
, where τs ∈ (τ s , τ s ) for s = l, r.

Although being simpler, the above condition is not very useful for practical application since it remains implicit by construction. We acknowledge that a similar condition has been derived in [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF] in the framework of relaxation schemes. We conclude this paragraph performing an asymptotic analysis of the entropy production (45) when the specific volume jump is small, i.e., τ s -τ s 1. In this case, the Taylor expansion of the right-hand side of ( 45) up to third-order leads to 3 , where τs ∈ (τ s , τ s ).

η s η s θ(τ s , η) dη = 1 2 z 2 s - a 2 s τ 2 s (τ s -τ s ) 2 + 1 6 ∂ 2 p ∂τ 2 (τ s , η s )(τ s -τ s )
This asymptotic analysis shows that the acoustic Godunov solver characterized by z s = a s τ s has an entropy production of third-order with respect to the specific volume jump. In addition, for a convex equation of state, the second partial derivative of pressure with respect to specific volume is strictly positive and the entropy production term is negative across compressive discontinuities which is not consistent with the second law of thermodynamics. For a perfect gas equation of state characterized by the polytropic index γ, the isentrope curve writes p(τ s , η s ) = p s τ s τ s γ for s = l, r. Substituting this function into (45) yields the expression of the entropy production term

η s η s θ(τ s , η) dη = a 2 s f τ s τ s , where f (x) = 1 γ(γ -1) (1 -x 1-γ ) + 1 γ (1 -x) + 1 2 z s ρ s a s 2 (1 -x) 2 .
We have displayed in figure 2 the entropy production with respect to the ratio τ s τ s for the three following values of the parameter: z s = ρ s a s , z s = 2ρ s a s and z s = 4ρ s a s . It is clear that the Godunov acoustic solver characterized by z s = ρ s a s exhibits a negative entropy production in the compressive zone, i.e., for τ s ≤ τ s , whereas for z s sufficiently greater than the acoustic impedance, ρ s a s , the entropy production remains positive almost everywhere in the compressive zone.

This concludes the design of the Riemann solver in the Lagrangian framework. Next we present its Eulerian counter-part.

Simple approximate Riemann solver for Eulerian gas dynamics

The Riemann problem for gas dynamics written under Eulerian form reads

∂u ∂t + ∂ f (u) ∂x = 0, u(x, 0) =        u l if x < 0, u r if x ≥ 0. ( 48
)
For the sake of completeness, let us recall that u = (ρ, ρu, ρe) t is the vector of conservative variables and f (u) = (ρu, ρu 2 + p, ρue + pu) t is the flux vector. The aim of this section is to derive a simple Riemann solver, w u l , u r , x t , approximating the solution of the Riemann problem [START_REF] Osher | Riemann solvers, the entropy condition and difference approximations[END_REF]. We acknowledge that there exists many approximate Riemann solvers for Eulerian gas dynamics, refer for instance to [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. Here, we present an Eulerian approximate Riemann solver which is deduced from its Lagrangian counterpart, W, constructed in paragraph 3.1. This Eulerian approximate Riemann solver is constructed following the general methodology introduced in [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Non-conservatifs[END_REF][START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF], which allows not only to derive Eulerian approximate Riemann solvers from their Lagrangian counterparts but also to transfer the positivity and entropy stability properties. Therefore, the present Eulerian approximate solver coincides with the one introduced in [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Non-conservatifs[END_REF][START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF] in a more general and theoretical framework. Let us point out that its has been derived employing a different and more intuitive approach. This Eulerian simple solver, w, consists of 4 constant states (u l , u l , u r , u r ) separated by the 3 discontinuities characterized by the speeds (S l , S , S r )

w u l , u r , x t =                          u l if x t < S l , u l if S l < x t < S , u r if S < x t < S r , u r if S r < x t . ( 49
)
The structure of the simple Riemann solver, w, is displayed in figure 3. Our starting point is the Lagrangian approximate Riemann solver, W, composed of the four states (U l , U l , U r , U r ), separated by the three discontinuities of speeds

      - z l ρ 0 l , 0, z r ρ 0 r       .
Here, we assume that the parameters z l , z r are such that the positivity of specific volumes, internal energies and entropy inequalities are taken for granted.

The construction of w is deduced directly from W. We have seen in section 2.4 that the Lagrangian (Σ) and the Eulerian (σ) speeds of an admissible discontinuity are related by where the superscripts l and r refer to the left and right sides of the discontinuity under consideration.

σ = τ l Σ + u l = τ r Σ + u r , x t 0 dx dt = S ⋆ dx dt = S r dx dt = S l u l u r u ⋆ r u ⋆ l -∆x l ∆x r
Mimicking the foregoing formula, we shall express the speeds of the Eulerian discontinuities, (S l , S , S r ), in terms of the speeds of the Lagrangian discontinuities (-

z l ρ 0 l , 0, z r ρ 0 r
). Replacing the mass flux swept by the Lagrangian discontinuity, Σ, respectively by -z l , 0 and z r we obtain the expression of the Eulerian speeds

S l = -z l τ l + u l = -z l τ l + u , (50a)
S = u , (50b)
S r = z r τ r + u = z r τ r + u r . ( 50c
)
We observe that the wave speeds S l , S and S r are fully determined by means of the knowledge of the intermediate states of the Lagrangian simple solver since by construction z l (τ l -τ l ) -(uu l ) = 0 and z r (τ r -τ r ) + uu r = 0 holds true, refer to (33a) and (34a). Finally, comparing (50a), (50b) and (50c), it is clear that the the wave speeds ordering

S l ≤ S ≤ S r , (51)
holds true provided that the parameters z l , z r and the specific volumes τ l , τ r are strictly positive. This is achieved when z l and z r are determined for instance by means of condition [START_REF] Kemm | A simple diffuse interface approach for compressible flows around moving solids of arbitrary shape based on a reduced Baer-Nunziato model[END_REF]. Assuming that this condition holds true, we define the intermediate mass density

ρ s = 1 τ s . Substituting it into (50a), (50b)
and (50c) allows us to express the Lagrangian mass fluxes (z l , 0, z r ) in terms of the Eulerian wave speeds as follows

z l = ρ l (u l -S l ) = ρ l (u -S l ), ( 52a
) 0 = ρ l (u -S ) = ρ r (u -S ), (52b)
z r = -ρ r (u -S r ) = -ρ r (u r -S r ). ( 52c
)
We conclude that the mass flux swept by the Lagrangian discontinuities has exactly the opposite sign to the mass flux crossing the Eulerian discontinuities. This result has been already observed in section 2.4 wherein we have derived the Eulerian Rankine-Hugoniot conditions from the Lagrangian ones, refer to [START_REF] Gallice | Système d'Euler-Poisson, magnétohydrodynamique et schémas de Roe[END_REF].

We are now in position to derive the intermediate states and fluxes of the Eulerian simple solver from the knowledge of the Lagrangian simple solver utilizing the expressions of the Lagrangian mass fluxes in terms of the Eulerian wave speeds. First, let us recall that the Lagrangian solver intermediate states and fluxes are determined by means of the jump relations [START_REF] Godunov | Finite difference methods for the computation of discontinuous solutions of the equations of fluid dynamics[END_REF]

z l (U l -U l ) + F l -F l = 0, 0(U r -U l ) + F r -F l = 0, -z r (U r -U r ) + F r -F r = 0.
Here, U s = (τ s , u s , e s ) t , U s = (τ s , u , e s ) t , F s = (-u s , p s , p s u s ) t and F s = (-u , p , p u ) t . The star quantities are fully determined by the systems of equations ( 33) and ( 34) once the parameters (z l , z r ) are prescribed. Now, substituting the expressions of z l , 0 and z r given respectively by (52a), (52b) and (52c) into the foregoing equations leads to

-S l (ρ l U l -ρ l U l ) + ρ l u U l + F l -(ρ l u l U l + F l ) = 0, (53a) -S (ρ r U r -ρ l U l ) + ρ r u U r + F r -(ρ l u U l + F l ) = 0, (53b) -S r (ρ r U r -ρ r U r ) + ρ r u r U r + F r -(ρ r u U r + F r ) = 0. ( 53c
)
These vectorial equations might be interpreted as jump equations across the Eulerian discontinuities of speeds (S l , S , S r ). However, one readily observes that the first component of the 3 foregoing equations is trivial. Thus, this first component cannot contribute to the definition of the intermediate states and fluxes of the approximate Riemann solver. On the other hand, rearranging (52a), (52b) and (52c) we arrive at

-S l (ρ l -ρ l ) + ρ l u -ρ l u l = 0, -S (ρ r -ρ l ) + ρ r u -ρ l u = 0, -S r (ρ r -ρ r ) + ρ r u r -ρ r u = 0.
This is nothing but the jump relations across the discontinuities of speeds (S l , S , S r ) related to the Eulerian mass conservation equation. Finally, gathering the second and the third components of system (53) with the foregoing system yields

-S l (u l -u l ) + f l -f l = 0, (54a) -S (u r -u l ) + f r -f l = 0, (54b) -S r (u r -u r ) + f r -f r = 0, (54c)
where on the one hand

u s = (ρ s , ρ s u s , ρ s e s ) t , f s = (ρ s u s , ρ s u 2 s + p s , ρ s u s e s + p s u s ) t ,
and on the other hand

u s = (ρ s , ρ s u , ρ s e s ) t , f s = (ρ s u , ρ s (u ) 2 + p , ρ s u e s + p u ) t , for s = l, r. (55)
The Eulerian intermediate states and fluxes are fully determined from the Lagrangian Riemann solver. Moreover, the second and third components of the Eulerian intermediate states and fluxes might be formally related to their Lagrangian counterparts by u s = ρ s U s and f s = ρ s u U s + F s for s = l, r. Summing equations (54a), (54b) (54c) leads to

f r -f l = S l (u l -u l ) + S (u r -u l ) + S r (u r -u r ), (56)
which shows the consistency of the approximate Riemann solver w with the integral form of the Riemann problem [START_REF] Osher | Riemann solvers, the entropy condition and difference approximations[END_REF]. Indeed, equation ( 56) is nothing but the integration of (48) over the space-time interval [-∆x l , ∆x r ] × [0, ∆t] displayed in figure 3. The interface flux, f , at x t = 0 is expressed in terms of the wave speeds as follows

f =                f l if 0 < S l < S < S r , f l if S l < 0 < S < S r , f r if S l < S < 0 < S r , f r if S l < S < S r < 0. ( 57
)
This expression can be rewritten under the compact form

f = 1 2 ( f l + f r ) - | S l | 2 (u l -u l ) - | S | 2 (u r -u l ) - | S r | 2 (u r -u r ). ( 58
)
The numerical dissipation is governed by the wave velocities S l , S and S r .

Remark 1 (Similarity with HLLC approximate Riemann solver). The present Eulerian approximate Riemann solver is very similar to the famous HLLC approximate Riemann solver introduced by Toro and his co-authors in the seminal work [START_REF] Toro | Restoration of the contact surface in the Harten-Lax-van Leer Riemann solver[END_REF] and also presented in details in [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. First, the structure of the approximate Riemann solver composed of 4 states separated by 3 waves is exactly the same. Second, the jump relations (54) satisfied by the intermediate states and fluxes across each wave are also identical to those satisfied by the states and fluxes of the HLLC approximate Riemann solver. Therefore, the structure of the resulting interface flux [START_REF] Toro | Anomalies of conservative methods: analysis, numerical evidence and possible cures[END_REF] is the same than the one of the HLLC solver. On the other hand, our Eulerian approximate Riemann solver has been entirely deduced from its Lagrangian counterpart utilizing the relationship (50) between Lagrangian and Eulerian wave speeds. This approach provides us a self-consistent computation of the Eulerian wave speeds, which are ordered by construction. This also allows the transfer of the properties (positivity preserving and entropy stability) of the Lagrangian approximate Riemann solver to its Eulerian counterpart. Let us point out that the wave speeds estimate in the context of the HLLC approximate Riemann solver has been and still remains an active topic of research, refer for instance to [START_REF] Batten | On the choice of wavespeeds for the HLLC Riemann solver[END_REF][START_REF] Toro | Bounds for Wave Speeds in the Riemann Problem: Direct Theoretical Estimates[END_REF].

It remains to assess the consistency of the Eulerian approximate Riemann solver, w, with the integral form of the Eulerian entropy inequality [START_REF] Clain | The mood method for the non-conservative shallow-water system[END_REF]. First, we assume that the Lagrangian approximate Riemann solver, W, is consistent with the integral form of the Lagrangian entropy inequality, namely [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF] holds true. This might be obtained prescribing parameters z l , z r to satisfy [START_REF] Nogueira | High accurate SPH method with multidimensional optimal order detection[END_REF]. Let us compute the entropy production across each wave by transforming its Lagrangian expression into its Eulerian counterpart simply by substituting the expressions of z l and z r in terms of S l and S r given by (52a) and (52c). For the leftgoing wave we arrive at

z l (η l -η l ) = ρ l (u -S l )η l -ρ l (u l -S l )η l = -S l (ρ l η l -ρ l η l ) + ρ l η l u -ρ l η l u l , (59)
whereas for the right-going wave we obtain

z r (η r -η r ) = -ρ r (u -S r )η r + ρ r (u r -S r )η r = -S r (ρ r η r -ρ r η r ) + ρ r η r u r -ρ r η r u . ( 60
)
Summing these two equations leads to

z l (η l -η l ) + z r (η r -η r ) = -S l (ρ l η l -ρ l η l ) -S (ρ r η r -ρ l η l ) -S r (ρ r η r -ρ r η r ) + ρ r η r u r -ρ l η l u l , (61)
since u = S thanks to (52b). Finally, we claim that the Eulerian approximate Riemann solver is consistent with the integral form of the Eulerian entropy inequality if and only if the Lagrangian approximate Riemann solver is consistent with the integral form of the Lagrangian entropy inequality, and, there holds

-S l (ρ l η l -ρ l η l ) -S (ρ r η r -ρ l η l ) -S r (ρ r η r -ρ r η r ) + ρ r η r u r -ρ l η l u l ≥ 0. ( 62
)
This concludes the design of a Riemann solver in Eulerian framework derived from its Lagrangian counter-part. Let us now design their associated Godunov numerical schemes.

t X i-1 2 X i+ 1 2 t n t n+1 U ⋆ l,i-1 2 U ⋆ r,i-1 2 U n i+1 U n i -z l,i-1 2 z r,i-1 2 -z l,i+ 1 2 z r,i+ 1 2 U ⋆ r,i+ 1 2 U n i-1 U ⋆ l,i+ 1 2 X i

Godunov-type Finite Volume scheme for one-dimensional Lagrangian gas dynamics

The aim of this section is to provide a first-order space and time Finite Volume discretization of the one-dimensional gas dynamics equations written under Lagrangian formalism. The resulting numerical method is a Godunov-type scheme since the updated cell-averaged value are computed by combining the cell-interface approximate Riemann solvers that have been constructed in section 3.1.

Governing equations and notation

The one-dimensional Lagrangian gas dynamics equations read

ρ 0 ∂U ∂t + ∂F ∂X = 0,
where U = (τ, u, e) t , F(U) = (-u, p, pu) t , e = ε + 1 2 u 2 and ρ 0 (X) > 0 is the initial mass density. The set of physically admissible states for this system of conservation laws writes

A = U = (τ, u, e) t , τ > 0, ε > 0 . ( 63
)
The computational domain

Ω = [X min , X max ] is partitioned into N c non overlapping cells [X i-1 2 , X i+ 1 2 ]
, where X i+ 1 2 denotes the position of a generic node. The cell mass is constant and defined by

m i = X i+ 1 2 X i-1 2 ρ 0 (X) dX.
At time t n , we assume that the solution of the foregoing system of conservation laws is piecewise constant over each cell and defined by the mass-averaged value

U n i = 1 m i X i+ 1 2 X i-1 2 ρ 0 (X)U(X, t n ) dX. ( 64
)

Godunov-type scheme

Let us compute the discrete solution of the system of conservation laws at time t n+1 = t n + ∆t in terms of the piece-wise solution at time t n , where ∆t > 0 denotes the time step. This will be done in combining the approximate Riemann solvers located at X i-1 2 and X i+ 1 2 , refer to figure 4. The time step, ∆t, has been chosen sufficiently small to ensure that the rightgoing wave emanating from X i-1 2 does not interact with the leftgoing wave emanating from X i+ 1 2 . The updated mass-averaged value writes

m i U n+1 i = X i X i-1 2 ρ 0 W       (U n i-1 , U n i , X -X i-1 2 ∆t       dX + X i+ 1 2 X i ρ 0 W       U n i , U n i+1 , X -X i+ 1 2 ∆t       dX, ( 65
)
where W is the approximate Riemann solver defined by [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Non-conservatifs[END_REF] and

X i = 1 2 (X i-1 2 + X i+ 1 2 )
. Performing the computation of the integrals, we arrive at the Godunov-type Finite Volume scheme for Lagrangrian gas dynamics

m i U n+1 i = z r,i-1 2 ∆tU r,i-1 2 + [m i -(z r,i-1 2 + z l,i+ 1 2 )∆t]U n i + z l,i+ 1 2 ∆tU l,i+ 1 2 . ( 66
)
It is worth noting that from the above equation one can recover the classical Finite Volume written in flux form. To this end, let us rearrange (66) as follows

m i U n+1 i = m i U n i -z r,i-1 2 ∆t U n i -U r,i-1 2 + z l,i+ 1 2 ∆t U l,i+ 1 2 -U n i . (67)
Now, applying the jump relations (31c) and (31a) respectively at X i-1 2 and X i+ 1 2 leads to

-z r,i-1 2 U n i -U r,i-1 2 = F i-1 2 -F n i , z l,i+ 1 2 U l,i+ 1 2 -U n i = -F i+ 1 2 + F n i ,
which allows us to introduce the numerical fluxes at cell interfaces

F i-1 2 and F i+ 1 2 .
Finally gathering the foregoing results we arrive at the flux form Finite Volume scheme

U n+1 i = U n i - ∆t m i F i+ 1 2 -F i-1 2 . ( 68
)
The cell interface flux is uniquely defined by

F i+ 1 2 =               -u i+ 1 2 p i+ 1 2 p i+ 1 2 u i+ 1 2               . (69)
Here, u i+ 1 2 and p i+ 1 2 are obtained by solving [START_REF] Godunov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF] and [START_REF] Gottlieb | Strong stability-preserving high-order time discretization methods[END_REF] in terms of the parameters (z l , z r ) at X i+ 1 2 .

Time step monitoring

Dividing (66) by m i yields

U n+1 i = z r,i-1 2 ∆t m i U r,i-1 2 +       1 - (z r,i-1 2 + z l,i+ 1 2 )∆t m i       U n i + z l,i+ 1 2 ∆t m i U l,i+ 1 2 . ( 70
)
It is clear that

U n+1 i is a convex combination of (U r,i-1 2 , U n i , U l,i+ 1 2
) provided that the coefficients of the combination are non negative. This implies the following time step monitoring to ensure the convex combination property ∆t ≤ m i z r,i-

1 2 + z l,i+ 1 2 ≡ ∆t i . (71)
In the case of the Godunov acoustic solver the parameters (z l , z r ) are equal to the cell acoustic impedance, i.e., z r,i-1 2 = z l,i+ 1 2 = ρ n i a n i . Using mass conservation the cell mass writes m i = ρ n i ∆x n i , where

∆x n i = x n i+ 1 2 -x n i-1 2 20
denotes the Eulerian cell width, and the above time step condition boils down to the well known CFL-like condition ∆t ≤

∆x n i 2a n i . (72)
Here, the Eulerian position of the cell interfaces are determined through the use of the discrete trajectory equation, i.e., x n+1

i+ 1 2 = x n i+ 1 2 + ∆tu i+ 1 2 .
In practice we compute the time-step as the minimum over all cells of the values ∆t i supplemented with a security factor 0 < CFL ≤ 1: are respectively convex combinations of (τ

∆t = CFL min i        m i z r,i-1 2 + z l,i+ 1 2        . (73
r,i-1 2 , τ n i , τ l,i+ 1 2
), (u

r,i-1 2 , u n i , u l,i+ 1 2
) and (e

r,i-1 2 , e n i , e l,i+ 1 2
).

First, let us consider the case of the specific volume, assuming that τ n i > 0. This means that the positivity of τ n+1 i holds true provided that τ r,i-1 2 > 0 and τ l,i+ 1 2 > 0. These latter conditions are satisfied when the parameters (z l , z r ) are determined according to the explicit formulas derived in section 3.2.2. Now, let us investigate the positivity of

ε n+1 i = e n+1 i - 1
2 (u n+1 i ) 2 . Scheme (70) provides

u n+1 i =α r,i-1 2 u i-1 2 + α i u n i + α l,i+ 1 2 u i+ 1 2 , e n+1 i =α r,i-1 2 e r,i-1 2 + α i e n i + α l,i+ 1 2 e l,i+ 1 2 ,
where α r,i-

1 2 = z r,i-1 2 ∆t m i , α l,i+ 1 2 = z l,i+ 1 2 ∆t m i and α i = 1 -(α r,i-1 2 + α l,i+ 1 2 ). Convexity of function x -→ x 2 implies 1 2 u n+1 i 2 ≤ 1 2 α r,i-1 2 u i-1 2 2 + 1 2 α i u n i 2 + 1 2 α l,i+ 1 2 u i+ 1 2 2
.

This shows the loss of kinetic energy induced by the averaging process of the cell velocity. This is probably the main source of numerical dissipation inherent to the Godunov scheme. Subtracting the kinetic energy to the total energy at time t n+1 leads to the updated internal energy which therefore satisfies

ε n+1 i ≥ α r,i-1 2 ε r,i-1 2 + α i ε n i + α l,i+ 1 2 ε l,i+ 1 2 .
This latter equation shows that the kinetic energy is converted into internal energy via the averaging procedure of the Godunov scheme. Finally, if the time step condition (71) holds true and if we assume that

ε n i > 0 then ε n+1 i > 0 provided that ε r,i-1 2 > 0 and ε l,i+ 1 2
> 0. These latter conditions are satisfied when the parameters (z l , z r ) are determined according to the explicit formulas derived in section 3.2.1.

Entropy inequality

The cell-averaged entropy [X i-1 2 , X i+ 1 2 ] at time t is given by

η i (t) = 1 m i X i+ 1 2 X i-1 2 ρ 0 (X)η (U(X, t)) dX. (74)
Knowing that U(X, t) is piecewise constant with respect to X, leads to write

η n i = η(U n i ) and η n+1 i = η U n+1 i
, where U n+1 i is computed from the Godunov-type scheme (70). Then, under the time step condition (71), the updated cell-averaged entropy satisfies

η n+1 i = η α r,i-1 2 U r,i-1 2 + α i U n i + α l,i+ 1 2 U l,i+ 1 2 , ≥ α r,i-1 2 η U r,i-1 2 + α i η U n i + α l,i+ 1 2 η U l,i+ 1 2 , ←-thanks to η concavity ≥ α r,i-1 2 η r,i-1 2 + α i η n i + α l,i+ 1 2 η l,i+ 1 2 .
Now, observing that α r,i-1 2 + α i + α l,i+ 1 2 = 1 we arrive at the inequality

η n+1 i -η n i ≥ α r,i-1 2 (η r,i-1 2 -η n i ) + α l,i+ 1 2 (η l,i+ 1 2 -η n i ). ( 75
)
This shows that the cell entropy increase with respect to time is controlled by the entropy increase across the right-sided and left-sided waves. Finally, the Godunov-type scheme satisfies the entropy inequality

η n+1 i -η n i ≥ 0, (76)
provided that the Riemann approximate solver is entropic, that is η

r,i-1 2 ≥ η n i and η l,i+ 1 2
≥ η n i . These latter conditions are taken for granted provided that the parameters (z l , z r ) satisfy the constraint (47).

1D Lagrangian numerical validation

Here, we present several classical test cases to assess the properties of the 1D Lagrangian scheme based on the approximate Riemann solver developed in section 3.1. The main purpose is to observe that the numerical method is stable and preserves the positivity of the specific volume and energy. The CFLis set to 0.9 otherwise noticed and used to run the classical Sod shock tube problem, an extreme double rarefaction problem and the Le Blanc shock tube. The initial conditions for these test cases are gathered in table 1. First of all, we perform a comparison between the three available methods for the determination of (z l , z r ).

Name

Left First, we run a test case with simple waves -the Sod shock tube and then, a test case with interacting waves -the Woodward-Collela blast wave test.

Sod shock tube. We run the planar Sod shock tube problem to assess the ability of the methods to capture one-dimensional simple waves. The exact solution for this one-dimensional Riemann problem can be derived using for instance [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. On the left of figure 5 we present the results for the density variable on a mesh refinement situation (N = 100 to 400 cells) and on the right the comparison of numerical density with different algorithms to determine the impedances are employed. The numerical solutions are compared to the exact solution (black line). Extremely few differences are observed for all three methods. Woodward-Collela. The Woodward-Collela blastwave test [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF] is a double shock tube case that simulates the interaction of simple waves. Two shock waves and two contact discontinuities develop and propagate towards the wall boundary conditions and reflect from them. These initial simple waves further interact creating a more complex flow pattern. The reference solution is obtained by a Lagrangian numerical scheme for 5000 cells (black line). In figure 6 we present the results for the density variable on a mesh refinement situation (N = 100 to 400 cells) using the first method to compute z l and z r . We then present the comparison of the numerical density obtained by different algorithms to determine the impedances in figure 7 on N = 200 cells. Figure 8 (a) to (c) shows a grayscale map of the normalized right impedance z r /ρ r a r for method 1, 2 and 3 respectively, illustrating the evolution of the right impedance of each cell at each time step. We observe that for all methods, the right impedance appears to match the trajectory of the right moving shock and its subsequent interactions. Method 3 in particular show more patterns and this is due to the construction of this method that is more sensitive. Because the differences in the numerical solution remain small, we select Method 1 for the rest of this work.

Le Blanc shock tube. This test is an extreme version of a shock tube (density jump is 10 3 , pressure jump is 10 9 ) generating violent waves, which, however are still simple waves that can be exactly computed [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. The initial conditions are recalled in table 1. The left and middle images in figure 9 present the results for the density and internal energy variables respectively, on a mesh refinement situation (N = 900 to 3600 cells) versus the exact solution (black line). The mesh refined solution seems to converge towards the exact solution with classical low convergence speed as is expected for this several Riemann problem. On the right is a grayscale map of the normalized right impedance z r /ρ r a r , illustrating the evolution of the right impedance for each cell at each time step. Once again we observe that the right impedance matches the trajectory of the right moving shock as expected for this problem involving only simple waves and only one shock wave.

Extreme double rarefaction. This problem is inspired from the 123 problem [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. The latter is one of benchmark tests presenting near vacuum state as it involves two rarefaction fans moving in opposite directions therefore emptying the central zone where a trivial steady contact discontinuity remains. Here we consider an extreme version where the initial condition on a computational domain Ω = [0, 2] is in table 1. Outflow boundary conditions are considered. In this extreme configuration with γ = 7/5, the resolution of the Riemann leads to values of density ρ * 1.124 × 10 -6 and pressure p * 1.875 × 10 -9 and ε * 1.192 × 10 -3 , see figure 10 at final time. In figure 10 we plot the results for successively refined grids. The numerical solutions are compared to the exact solution (black line). The spurious peak in the internal energy profile is a classical flaw of any numerical schemes due to entropy dissipation. Nevertheless we do not observe any lack of positivity for the numerical scheme, although the smallest numerical densities reach 8.65 × 10 -3 for 800 cells. Moreover the numerical solution seems to converge towards the exact one.

x i-1 2 x i+ 1 2 t n t n+1 u n i+1 S r,i+ 1 2 u n i-1 x i u n i S r,i-1 2 u ⋆ r,i+ 1 2 u ⋆ l,i+ 1 2 u ⋆ l,i-1 2 u ⋆ r,i-1 2 S l,i-1 2 S ⋆ i+ 1 2 S ⋆ i-1 2 S l,i+ 1 2

Godunov-type Finite Volume scheme for one-dimensional Eulerian gas dynamics

The aim of this section is to provide a first-order space and time Finite Volume discretization of the one-dimensional gas dynamics equations written under Eulerian formalism as a derived version from the Lagrangian one. The resulting numerical method is a Godunov-type scheme since the updated cell-averaged value are computed by combining the cell-interface approximate Riemann solvers that have been constructed in section 3.3.

Governing equation and notation

The one-dimensional Eulerian gas dynamics equations read

∂u ∂t + ∂ f ∂x = 0,
where u = (ρ, ρu, ρe) t , f (u) = (ρu, ρu 2 + p, ρeu + pu) t , where e = ε + 1 2 u 2 . The set of physically admissible states for this system of conservation laws writes

A E = u = (ρ, ρu, ρe) t , ρ > 0, ε > 0 . ( 77
)
The computational domain ω = [x min , x max ] is partitioned into n c non overlapping cells [x i-1 2 , x i+ 1 2 ], where x i+ 1 2 denotes the fixed position of a generic node. The cell volume is

∆x i = x i+ 1 2 -x i-1 2 .
At time t n , we assume that the solution of the foregoing system of conservation laws is piecewise constant over each cell and defined by the cell-averaged value

u n i = 1 ∆x i x i+ 1 2 x i-1 2 u(x, t n ) dx. ( 78
)

Godunov-type scheme

We define the discrete solution of the system of conservation laws at time t n+1 = t n + ∆t in terms of the piecewise solution at time t n combining the approximate Riemann solvers located at x i-1 2 and x i+ 1 2 , refer to figure 11. The time step, ∆t, has been chosen sufficiently small to ensure that the rightgoing wave emanating from x i- 1 2 does not interact with the leftgoing one emanating from x i+ 1 2 . Then, the updated cell-averaged value writes

∆x i u n+1 i = x i x i-1 2 w u n i-1 , u n i , x -x i-1 2 ∆t dx + x i+ 1 2 x i w u n i , u n i+1 , x -x i+ 1 2 ∆t dx. ( 79
)
The first term at the right-hand side corresponds to the integral of the approximate Riemann solver at

x i-1 2 over [x i-1 2 , x i ] where x i = 1 2 (x i-1 2 + x i+ 1
2 ), whereas the second term corresponds to the integral of the approximate Riemann solver at x i+ 1 2 over [x i , x i+ 1 2 ]. After some algebra we arrive at the explicit expressions of these integrals

x i x i-1 2 w u n i-1 , u n i , x -x i-1 2 ∆t dx = S + l,i-1 2 ∆t u n i-1 + (S ,+ i-1 2 -S + l,i-1 2 )∆t u l,i-1 2 + (S + r,i-1 2 -S ,+ i-1 2 )∆t u r,i-1 2 + (x i -x i-1 2 -S + r,i-1 2 ∆t)u n i , x i+ 1 2 x i w u n i , u n i+1 , x -x i+ 1 2 ∆t dx = (x i+ 1 2 -x i -S - l,i+ 1 2 ∆t)u n i + (S - l,i+ 1 2 -S ,- i+ 1 2 )∆t u l,i+ 1 2 + (S ,- i+ 1 2 -S - r,i+ 1 2 )∆t u r,i+ 1 2 + S - r,i+ 1 2 ∆t u n i+1 .
Here, S l,i+ 1 2 , S i+ 1 2 and S r,i+ 1 2 are the waves speeds of the Eulerian approximate Riemann solver at x i+ 1 2 . These Eulerian speeds are computed from the underlying Lagrangian approximate Riemann solver according to respectively (50a), (50b) and (50c). In addition, x + = 1 2 (| x | +x) denotes the positive part of the real x, whereas x -= 1 2 (| x | -x) denotes its negative part. These are non negative functions of x respectively increasing and decreasing. We observe that the positive part (resp. negative part) of the waves emanating from x i-1 2 (resp. x i+ 1 2 ) contributes to update the solution at time t n+1 . Summing the foregoing results and collecting the terms in factor of the intermediate states, we finally obtain

u n+1 i = ∆t ∆x i S + l,i-1 2 u n i-1 + ∆t ∆x i (S ,+ i-1 2 -S + l,i-1 2 )u l,i-1 2 + ∆t ∆x i (S + r,i-1 2 -S ,+ i-1 2 )u r,i-1 2 (80) + 1 - ∆t ∆x i (S + r,i-1 2 + S - l,i+ 1 2
)

u n i + ∆t ∆x i (S - l,i+ 1 2 -S ,- i+ 1 2 )u l,i+ 1 2 + ∆t ∆x i (S ,- i+ 1 2 -S - r,i+ 1 2 )u r,i+ 1 2 + ∆t ∆x i S - r,i+ 1 2 u n i+1 .
The underlying approximate Riemann solver, w, being consistent with the integral form of the system of conservation laws, the foregoing numerical scheme is of Godunov-type and can be equivalently written under the flux form

∆x i (u n+1 i -u n i ) + ∆t ( f i+ 1 2 -f i-1 2 ) = 0, (81)
where f

i+ 1 2
is the numerical flux at interface x i+ 1 2 resulting from the Eulerian approximate Riemann solver. According to [START_REF] Toro | Anomalies of conservative methods: analysis, numerical evidence and possible cures[END_REF] the expression of the interface flux, with an obvious notation adaptation, reads

f i+ 1 2 = 1 2 ( f n i + f n i+1 ) - | S l,i+ 1 2 | 2 (u l,i+ 1 2 -u n i ) - | S i+ 1 2 | 2 (u r,i+ 1 2 -u l,i+ 1 2 ) - | S r,i+ 1 2 | 2 (u n i+1 -u r,i+ 1 2
).

(82)

Time step monitoring

We shall compute the time step, ∆t, to ensure that the updated cell-averaged value,

u n+1 i , is a convex combination of u n i-1 , u l,i-1 2 , u r,i-1 2 , u n i , u l,i+ 1 2 , u r,i+ 1 2
and u n i+1 . This amounts to investigate the positivity of the coefficients of the linear combination present at the right-hand side of (80). Firstly, the positive part and the negative part of the wave speeds are non negative. Secondly, the wave speeds of the approximate Riemann solvers at x i-1 2 and x i+ 1 2 are ordered according to S l,i±

1 2 ≤ S i± 1 2 ≤ S r,i± 1 2 .
Bearing this in mind and recalling that x → x + (resp. x → x -) is an increasing (resp. decreasing) function yields

S + l,i± 1 2 ≤ S ,+ i± 1 2 ≤ S + r,i± 1 2 , and S - l,i± 1 2 ≥ S ,- i± 1 2 ≥ S - r,i± 1 2 .
Finally, all the coefficients at the right-hand side of (80) are unconditionally positive except the one in factor of u n i . This latter coefficient is positive if the time step ∆t satisfies

∆t ≤ ∆x i S + r,i-1 2 + S - l,i+ 1 2 . ( 83
)
Therefore u n+1 i is a convex combination of u n i and the intermediate states of the approximate Riemann solvers located at x i± 1 2 provided that ∆t satisfies condition (83), which can be viewed as a CFL-like condition. Developing the expression of the Eulerian wave speeds in terms of the underlying Lagrangian mass fluxes, i.e., S r,i-

1 2 = u n i + z r,i-1 2 ρ n i and S l,i+ 1 2 = u n i - z l,i+ 1 2 ρ n i , we get S + r,i-1 2 + S - l,i+ 1 2 = 1 2 | u n i - z l,i+ 1 2 ρ n i | + | u n i + z r,i-1 2 ρ n i | + z r,i-1 2 ρ n i + z l,i+ 1 2 ρ n i ≤| u n i | + z r,i-1 2 + z l,i+ 1 2 ρ n i .
Employing this lower bound leads to replace (83) by a more explicit condition

∆t ≤ ∆x i | u n i | + z r,i-1 2 + z l,i+ 1 2 ρ n i . (84)
Likewise for the Lagrangian case the time-step is practically computed as

∆t = CFL ∆x i max i | u n i | + z r,i-1 2 + z l,i+ 1 2 ρ n i . ( 85
)
In the following we focus on the positivity and entropy issues.

Positivity of mass density and specific internal energy

The Finite Volume scheme (81) for Eulerian gas dynamics has been constructed through the use of the Eulerian approximate Riemann solver described in section 3.3, which inherits all the properties from its Lagrangian counterpart. More precisely, provided that the parameters (z l , z r ) of the Lagrangian Riemann solver satisfies condition [START_REF] Kemm | A simple diffuse interface approach for compressible flows around moving solids of arbitrary shape based on a reduced Baer-Nunziato model[END_REF] the positivity of specific volume and specific internal energy is ensured. Therefore, the waves speeds of the Riemann solvers are ordered according to [START_REF] Reed | Triangular mesh methods for neutron transport equation[END_REF] and the positivity of mass density and specific energy for the intermediate states is also ensured. Under the time step condition (83), u n+1 i is a convex combination of u n i and the intermediate states of the approximate Riemann solvers located at the cell interfaces, refer to (80). Therefore, if ρ n i > 0 and ε n i > 0 then ρ n+1 i > 0 and ε n+1 i > 0. This shows that Finite Volume scheme (81) preserves the positivity of mass density and specific internal energy under the time step condition (83).

Entropy inequality

In this paragraph, we aim at constructing the discrete entropy inequality satisfied by the Godunov-type Finite Volume scheme (81) following [START_REF] Chalons | Nonlinear Hyperbolic Systems of Conservation Laws[END_REF]. We shall start by recalling the expression of the updated cell-

averaged value u n+1 i u n+1 i = ∆t ∆x i S + l,i-1 2 u n i-1 + ∆t ∆x i (S ,+ i-1 2 -S + l,i-1 2 )u l,i-1 2 + ∆t ∆x i (S + r,i-1 2 -S ,+ i-1 2 )u r,i-1 2 + 1 - ∆t ∆x i (S + r,i-1 2 + S - l,i+ 1 2
)

u n i + ∆t ∆x i (S - l,i+ 1 2 -S ,- i+ 1 2 )u l,i+ 1 2 + ∆t ∆x i (S ,- i+ 1 2 -S - r,i+ 1 2 )u r,i+ 1 2 + ∆t ∆x i S - r,i+ 1 2 u n i+1 .
Let us recall that under time step condition (83),

u n+1 i is a convex combination of u n i-1 , u l,i-1 2 , u r,i-1 2 , u n i , u l,i+ 1 2 , u r,i+ 1 2
and u n i+1 . The strict concavity of u → (ρη)(u) with respect to u, refer to [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF], and Jensen's inequality yield the following inequality for the cell-averaged entropy time-increment

(ρη)(u n+1 i ) -(ρη)(u n i ) ≥ - ∆t ∆x i S + l,i-1 2 (ρη)(u l,i-1 2 ) -(ρη)(u n i-1 ) (86)
-

∆t ∆x i S ,+ i-1 2 (ρη)(u r,i-1 2 ) -(ρη)(u l,i- 1 2
)

+ S + r,i-1 2 (ρη)(u n i ) -(ρη)(u r,i- 1 2
)

+ ∆t ∆x i S - l,i+ 1 2 (ρη)(u l,i+ 1 2
)

-(ρη)(u n i ) + ∆t ∆x i S ,- i+ 1 2 (ρη)(u r,i+ 1 2
)

-(ρη)(u l,i+ 1 2 ) + S - r,i+ 1 2 (ρη)(u n i+1 ) -(ρη)(u r,i+ 1 2 )
.

It remains to exhibit a lower bound to the right-hand side of the foregoing inequality. This is achieved invoking the consistency of the underlying Eulerian approximate Riemann solver with the integral form of the Eulerian entropy inequality [START_REF] Clain | The mood method for the non-conservative shallow-water system[END_REF]. Namely, (62) holds true

-S l (ρ l η l -ρ l η l ) -S (ρ r η r -ρ l η l ) -S r (ρ r η r -ρ r η r ) + ρ r η r u r -ρ l η l u l ≥ 0.
Substituting the foregoing inequality at the right-hand side of (86) for the terms located at x i-1 2 (resp. x i+ 1 2 ) and developing the expressions of the negative and the positive parts of the wave speeds, we finally obtain the cell entropy inequality

(ρη)(u n+1 i ) -(ρη)(u n i ) + ∆t ∆x i (ρηu) i+ 1 2 -(ρηu) i-1 2 ≥ 0, (87)
where the entropy flux, ρηu, at

x i-1 2 writes (ρηu) i+ 1 2 = 1 2 (ρηu)(u n i ) + (ρηu)(u n i+1 ) - | S l,i+ 1 2 | 2 (ρη)(u l,i+ 1 2
)

-(ρη)(u n i ) (88) - | S i+ 1 2 | 2 (ρη)(u r,i+ 1 2
)

-(ρη)(u l,i+ 1 2 ) - | S r,i+ 1 2 | 2 (ρη)(u n i+1 ) -(ρη)(u r,i+ 1 2
) .

This concludes the design of the Eulerian scheme. The next section focuses on its numerical validation.

1D Eulerian numerical validation

In this section, we test the first-order Finite Volume scheme which has been previously described for Eulerian gas dynamics. Here, our main purpose is to ensure that this Eulerian numerical method preserves the theoretical properties and behaves accordingly. We fix the CFL to 0.9 and the initial conditions of the test cases are all gathered in table 1. Modified Sod shock tube. The modified Sod test simulates the impact of a dense fluid onto a light one at rest. The simulation is ran on the interval of [0, 1] with the interface located at x = 0.3. The initial set up reads ρ l = 1, u l = 0.75, p l = 1, ρ r = 0.125, u r = 0, p r = 0.1. A piston boundary condition is imposed on the left with u = u l = 0.75, while u = 0 is imposed on the right. This test case ensures that the sonic glitch in the rarefaction wave is absent. The numerical results at final time t final = 0.2 with N = 100 to 400 cells are presented in figure 12. The numerical solution seems to converge towards the exact one. More important we do not observe any sonic glitch in the rarefaction wave, see [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF].

Collela-Woodward blastwave. In figure 13, we present the numerical density obtained for N = 800 to 6400 cells, with a zoom of the central area. The quality of the solution improves with the mesh refinement. Then we present on the middle panel the time-step evolution which clearly adapts to the numerical solution. On the right panel, the plot of the normalized right impedance as a function of time allows us again to track the main shock waves and their interactions. Le Blanc shock tube. We simulate the Le Blanc shock tube with the Eulerian scheme. Figure 14 (left panel) shows the numerical density plotted in the log scale for the Le Blanc shock tube problem for N = 900, 1600 and 3200 cells. The panel in the middle presents the time-step evolution showing that the main adaptation is related to the separation of the simple waves Then the evolution of the normalized right impedance that traces the path of the right moving shock is presented on the right panel. Again beause only one single shock is present in this problem, then only one curves can be observed. Extreme double rarefaction. We present the results for the extreme double rarefaction problem in figure 15. No positivity issue is encountered and the minimal density is about ρ i = 1.48 × 10 -3 for N = 800 cells (ρ i = 2.43 × 10 -3 for 400 cells and ρ i = 9.09 × 10 -3 for 100 cells). The mesh convergence results for the density, velocity and specific internal energy show that the numerical scheme improves the quality of the numerical solution when the mesh size decreases. Test case suite of Toro . In his book [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] Toro has gathered some now-classical test cases. All tests are simulated on domain Ω = [0, 1] with γ = 7/5, a CFL set to 0.9 and the number of cells is 100. The initial configurations are recalled in table 2 where the discontinuity is located at x 0 and the final time is set to t final .

We refer the readers to [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] for the phenomenological description of the test cases. Remark that the test labeled 1 corresponds to modified Sod test case that we have already presented. The results are gathered in figure 16 for tests 2 to 5 for the current Eulerian scheme. We present the density, velocity and internal energy variables against the exact solution. The present numerical results are perfectly consistent with what can be obtained by a first-order accurate Eulerian Finite Volume scheme. Next in figure 17 we present the numerical results obtained by the current scheme and the classical HLL one on tests 6 and 7, which correspond respectively to a stationnary contact wave and a moving contact wave. The current numerical scheme is able to exactly resolve a stationary contact by construction (left panel), and can capture the moving contact discontinuity but with some numerical dissipation. Notice that the HLL scheme has no embedded mechanism to capture such stationary contact, while the family of HLLC schemes has, likewise for our approach.

Label

Left We would like to emphasize that the Lagrangian scheme naturally adapts its mesh to the flow: some smaller and compressed cells are encountered after a shock wave and expanded ones after a rarefaction wave. Moreover moving contact discontinuities are exactly preserved. This however comes with the drawback of entropy deposition even in rarefaction waves. On the contrary, an Eulerian scheme considers a fixed mesh, hence ensuring a better accuracy in rarefaction waves, but less accurate shocks and diffused moving contact discontinuities. However, the extension to multi-dimensions and high-orders of accuracy is somewhat simpler with Eulerian schemes since in this formalism the mesh is not moving. This set of test cases have validated the first order Eulerian FV scheme and its theoretical behaviors are numerically reproduced or observed. The next step of this work consists in using this scheme as a so-called 'parachute' scheme in an a posteriori MOOD limited high-order scheme.

Extension to high order of accuracy

In this section we present an extension of the previous first-order Eulerian scheme into a high-order accurate space/time scheme. Time discretization proceeds via classical Strong Stability preserving Runge-Kutta (SSPRK) scheme [START_REF] Gottlieb | Strong stability-preserving high-order time discretization methods[END_REF][START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF]. The space discretization relies on classical polynomial reconstruction with not-so-classical a posteriori MOOD loop limiting procedure, see [START_REF] Clain | A high-order finite volume method for systems of conservation laws -multi-dimensional optimal order detection (MOOD)[END_REF]. The MOOD paradigm for FV schemes is built upon a try-and-fail concept. Starting from data, (u n i ) 1≤i≤N at t n a high-order scheme produces a 'candidate' solution (u * ,n+1 i ) 1≤i≤N at t n+1 . This candidate solution is tested against some 'detection' criteria which determine if the computed cell value at t n+1 is physically and numerically valid or not. In the case of validity the cell is accepted, otherwise the cell is flagged as troubled/bad and sent back to t n for re-computation. The re-computation is performed using a lower order scheme from a 'cascade' of schemes ordered from the more accurate to the more robust one. The latest scheme of this cascade is called the 'parachute' and must always produce a valid solution according to the detection criteria. This parachute scheme must be extremely robust and provably adapted to the detection criteria. In the previous sections we have designed such a first-order Eulerian scheme with good properties: conservation, robustness, positivity preservation under explicit CFL, entropy inequality, etc. The entire iterative procedure is called a MOOD loop and always converges to an acceptable numerical solution for which some cells have been updated with a high order scheme, while others with the parachute. Figure 18 shows an illustration to build a high-order numerical scheme with an a posteriori MOOD loop.

Several works using an a posteriori MOOD loop in different numerical contexts can be found for La- grangian, ALE, Discontinuous Galerkin limiter, Smoother-Particle-Hydrodynamics, Adaptative Mesh Refinement, WENO schemes, slope limiter, etc. [START_REF] Nogueira | High accurate SPH method with multidimensional optimal order detection[END_REF][START_REF] Vilar | A posteriori correction of high-order discontinuous galerkin scheme through subcell finite volume formulation and flux reconstruction[END_REF][START_REF] Boscheri | A second-order cell-centered lagrangian ADER-MOOD finite volume scheme on multidimensional unstructured meshes for hydrodynamics[END_REF][START_REF] Boscheri | Direct Arbitrary-Lagrangian-Eulerian ADER-MOOD finite volume schemes for multidimensional hyperbolic conservation laws[END_REF][START_REF] Zanotti | Space-time adaptive ADER Discontinuous Galerkin finite element schemes with a posteriori sub-cell finite volume limiting[END_REF][START_REF] Desveaux | An entropic MOOD scheme for the euler equations[END_REF][START_REF] Loubère | A new family of high order unstructured MOOD and ADER finite volume schemes for multidimensional systems of hyperbolic conservation laws[END_REF][START_REF] Farmakis | WENO schemes on unstructured meshes using a relaxed a posteriori MOOD limiting approach[END_REF][START_REF] Kitamura | Simple a posteriori slope limiter (post limiter) for high resolution and efficient flow computations[END_REF], and, physical contexts like Euler or Navier-Stokes equations, turbulence simulation, shallow-water equations, astrophysics, MHD, Baer-Nunziato multi-phase model, etc. [START_REF] Jiang | Efficient methods with higher order interpolation and MOOD strategy for compressible turbulence simulations[END_REF][START_REF] Clain | The mood method for the non-conservative shallow-water system[END_REF][START_REF] Zanotti | Space-time adaptive ADER Discontinuous Galerkin finite element schemes with a posteriori sub-cell finite volume limiting[END_REF][START_REF] Fambri | Space-time adaptive ADER-DG schemes for dissipative flows: Compressible Navier-Stokes and resistive MHD equations[END_REF][START_REF] Kemm | A simple diffuse interface approach for compressible flows around moving solids of arbitrary shape based on a reduced Baer-Nunziato model[END_REF] .

High-order space discretization

Let us denote by o = d + 1 > 0 the target order of accuracy fixed by the user. At discrete time t n the polynomial reconstruction operator R in cell Ω i is the operator which given a large enough set of 2K i neighbor mean values on stencil S d i = {u i±k , k = 1, . . . , K i } produces the polynomial u i ∈ P d (Ω). We skip the time exponent in this section as no confusion can occur. This polynomial has a degree d ≥ 0 and is such that its mean value exactly matches u i and best fits the neighbor mean values u i±k , that is

R : x ∈ Ω i , u i , d, S d i -→ u i (x). ( 89
)
Indeed it gathers several useful properties: mathematical firmly based grounds, Lagrangian↔Eulerian mapping, demonstrated ordered wave speeds, positivity under explicit CFL condition, entropy inequality and validation of its numerical behaviors on classical test cases (carbuncle-free, robustness, efficiency);

(iii) a cascade of schemes ranked from the more accurate and prone to spurious instability to the more robust one, that is the parachute scheme. This cascade may be parametrized by the polynomial degree d, for instance we can consider the following cascade:

d = 4 → 3 → 2 → 1 → 0;
Detection criteria. The physical admissibility criteria P are firmly based on the system of PDEs solved, they are obviously the positivity of the density and internal energy for the Eulerian scheme, hence

P : ρ * i > 0, and ε * i > 0 (92)
Next, the numerical criteria ensuring the ENO behavior are more subjective and somewhat attached to the developer preferences. Here we rely on a relaxed discrete maximum principle (RDMP) but applied to the output data of the nodal solver. Indeed the nodal solver is (one of) the constitutive brick, the locality of which is restricted to the neighbor cells. Moreover the outcome of the nodal solver are the nodal velocity and sub-pressures, which further serve to determine the impedances, and, ultimately the star states and the numerical fluxes to update the state vector. A candidate variable q * p fulfills the RDMP if

N : -δq p + min k∈C(p) q n k ≤ q * p ≤ max k∈C(p) q n k + δq p , (93)
where C(p) is the neighborhood associated to point p, δq p = max 10 -4 , 10 -3 (M pm p ) , and, M p and m p are the maximal and minimal local values of q m p = min k∈C(p) q n k and M p = max k∈C(p) q n k . If a nodal solver result does not fulfill (93) then the surrounded cells are marked as 'bad' ones. At last we also test if the numerical method has generated any un-representable value such as NaN. If a cell has been declared as 'invalid', then only this cell and its neighbors are sent back to t n for a recomputation with local lower order reconstructions. Possibly some cells are re-computed several times, and, may ultimately be updated with the parachute scheme. The troubled/bad cells are the cells which do not fulfill at least one of these detection criteria

D :         
P : physically admissible, i.e positivity under explicit ∆t with respect to (92) C : representable in a computer, i.e no Nan. N : numerically valid, i.e essentially non oscillating with respect to (93).

The two neighbors of a bad cell are also recomputed, and, as such flagged as bad ones without polynomial degree decrementing. This type of limiting is referred to an a posteriori procedure because the candidate solution at t n+1 is always confronted to the detection criteria and must pass them. Fortunately in general only few cells need to be re-computed which renders this a posteriori MOOD loop relatively efficient in comparison to other limiting approaches.

MOOD loop. The MOOD loop always converges in a finite number of iterations. Indeed there is only a few and fixed number of schemes in the cascade and a finite number of cells. Moreover a cell is always considered as valid if its reconstruction degree has dropped down to d i = 0. Consequently the number of MOOD loop iterations is necessarily finite. The final valid solution is then constituted of pieces of high accurate and low accurate approximations fulfilling the detection criteria. As a consequence the detection criteria do play a paramount role in the overall quality of the solution.

High-order time discretization

For the time discretization we rely on successive explicit sub-steps, which altogether allows to reach a nominal dth accuracy. Depending on the order of accuracy d, we employ the family of classical Runge-Kutta (RK) explicit schemes [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF]. For d = 0 a forward Euler scheme is used, while for d = 1 an explicit midpoint method (predictor-corrector or equivalently RK2 scheme) is employed. When d = 2 the classical SSPRK3 scheme is chosen. Beyond d ≥ 3 there exists no SSPRK method anymore and we then rely either on RK3 with appropriate time-step limitation or on a RK4 scheme, see table 4 for the associated Butcher tableaus. We denote by K > 1 the number of RK iterates of the time discretization scheme, each indexed by κ and 1 ≤ κ ≤ K. The restriction on the time step ∆t is the one derived in the first sections which ensures the positivity with a safety coefficient CFL ≤ 1. The time step ∆t ≡ ∆t κ=1 is determined at the end of the first RK iterate as a function of the minimal/maximal impedances z κ=1 l,i+1/2 and z κ=1 r,i+1/2 . It may happen that, during the physical evolution the time step required at RK iterate κ is such that: ∆t κ < ∆t. In this case the cell is flagged as a 'bad' one and will be limited (in space and time) by the MOOD loop, possibly dropping its time accuracy to a first order scheme for which ∆t ≡ ∆t κ=1 always ensures the positivity.

c 1 a 11 a 12 • • • a 1s c 2 a 21 a 22 • • • a 2s . . . . . . . . . . . . c s a s1 a s2 • • • a ss b 1 b 2 • • • b s 0 0 0 1 1 0 1/2 1/2 0 0 0 0 1 1 0 0 1/2 1/4 1/4 0 1/6 1/6 2/3 0 0 0 0 0 1/2 1/2 0 0 0 1/2 0 1/2 0 0 1 0 0 1 0 1/6 1/3 1/3 1/6 Generic RK RK2 SSPRK3 RK4

1D numerical validation of high order extension.

In the previous numerical sections we have tested the first order Lagrangian and Eulerian schemes. Here we solely focus on test cases with more complex interacting wave problems, such as the Collela-Woodward blastwave and the Shu-Osher test. Above all, we perform a convergence analysis of our numerical scheme using a smooth solution of the gas dynamics equations.

The following schemes are under scrutiny

• Eul-0: 1st order Eulerian scheme;

• Eul-X: (X+1)th order (unlimited) Eulerian scheme;

• Eul-XMOOD: (X+1)th order Eulerian scheme with MOOD limiting, with Eul-0 as parachute scheme.

which will be compared to some classical ones such as: Lag-0, the 1st order Lagrangian scheme from section 4, and, Eul-1Lim, the nominally 2nd order Eulerian scheme i.e with piece-wise linear reconstruction on primitive variables supplemented with classical van Leer slope limiter.

Density sine wave. This paragraph shows that the high order reconstructions and RK schemes allow respectively for 1 st , 2 nd , 3 rd , 4 th and 5 th order convergence in space for a smooth solution of the 1D compressible Euler equations. We consider the advection of a sine wave on a [0, 1] mesh with periodic boundaries. The sine wave is characterized by: The results of the convergence analysis performed with polynomial reconstruction going from d = 0 to d = 4 are displayed in table 5, and the L 2 error is shown in figure 19. The expected rate of convergence is reached and confirms the high-order accuracy for our numerical scheme.

ρ 0 = 1 + 0.1 sin(2πX) , u 0 = 1 , p 0 = 1/γ , (
Next we make use of this test case to compare the efficiency of the high-order schemes, i.e. the accuracy as a function of the computation time. In figure 20, we plot the relative L 2 error versus the CPU time for the Eul-0, Eul-1, Eul-2, Eul-3 and Eul-4 schemes. Some efficiency plots have been extrapolated to ease the visual comparison. In table 6, the efficiencies of the different numerical schemes are compared for a fixed error (ε 0 = 10 -8 ), and, then for a fixed execution time (10ms). For large enough CPU time or small error the high-order schemes out-perform low-order ones.

Eul-0

Eul- Modified smooth solution problem. The second test case is a smooth test case derived in the isentropic case, with a polytropic index of γ = 3 for the perfect gas EOS. In this situation, the characteristic curves of the Euler equations become straight lines. The governing equations reduce to two Burgers equations, allowing us to solve this problem analytically. For the numerical simulation, the initial condition is modified to yield a more challenging example, as

X ∈ [0, 1] , ρ 0 = 1 + 0.1 sin(2πX) , u 0 = 0 , p 0 = ρ γ 0 , (95)
with periodic boundary conditions. The initial density and pressure being very close to zero would confirm the positivity-preserving property of the numerical scheme. Figure 21 display the results of the third-order scheme with 100 cells on the [0, 1] domain with a final time of t final = 0.8 and CFL = 0.1. Making use of the analytical solution, the global truncation error corresponding to the first-to fourth-order scheme are shown in table 7 and in figure 22. These results confirm the high accuracy of this scheme. This problem involves small scales after the shock has interacted with the sine wave. These small scale features can be captured either with a fine enough mesh or with high order accurate method. Figure 23 shows the comparison of numerical solution using different polynomial reconstruction from Eul-0 to Eul-4MOOD and a reference solution with 200, 400 and 800 cells. As expected the accuracy of the numerical solution improves with the order of polynomial reconstruction. The reference solution taken here is the numerical solution of the Eul-0 scheme with N = 10000 cells (the CPU time is about 6100µs). Figure 24 exhibits the associated cell reconstruction type with each time step and as expected, the shock wave is captured in its motion. Then in table 8 we present the percentages of cells updated with each polynomial degree from 0 to 4 and for three different meshes. We observe that about 90% of cells are updated with the highest order scheme, the 10% left are updated with generally the parachute scheme. In the very same table we gather the CPU time needed to get the final solution as a function of the nominal accuracy of the scheme employed. d = 0 means that the first order scheme is run, while d = 4 means the 5th order MOOD scheme is used. As is observed the cost of second order scheme is about 3 times more expensive than a first order one (respectively 6, 7 and 9 times for third, fourth and fifth order ones).

Woodward-Colella blastwave. The complex flow pattern in the Woodward-Collela blastwave problem is a good test case to evaluate the high-order extension of our numerical scheme. Similar to the Shu-Osher problem, we present the numerical density with a zoom at the central area obtained by the Eul-0, Eul-1MOOD, Eul-2MOOD and Eul-3MOOD schemes compared with a reference solution with 200, and 400 cells in figure 25. The reference solution taken is the numerical solution computed with the Lagrangian solver. Once

Conclusion and Perspectives

Following the methodology introduced in the pioneering work [START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF], we have constructed a Lagrangian simple approximate Riemann solver which preserves contact discontinuities. The monitoring of wave speeds allows to derive explicit conditions for ensuring positivity and entropy stability. The resulting Godunov-type Finite Volume scheme preserves the positivity of specific volume and internal energy and also satisfies an entropic inequality provided a time step condition is satisfied.

Employing the general formalism described in [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Non-conservatifs[END_REF][START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF] an Eulerian simple approximate Riemann solver is deduced from its Lagrangian counterpart. In this framework, the Eulerian wave speeds are deduced from the Lagrangian ones and are naturally ordered provided that the Lagrangian approximate Riemann solver preserves the positivity of specific volume. Moreover, this framework provides also the transfer of the good properties (positivity and entropy stability) of the Lagrangian approximate Riemann solver to its Eulerian counterpart. We observe that the proposed approximate Riemann solver is nothing but the one initially introduced in [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Non-conservatifs[END_REF][START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF] in a more general and theoretical context. Here, we have proposed a less formal derivation based on more intuitive arguments. We also want to acknowledge that the numerical flux induced by present Eulerian approximate solver has the same structure than the one induced by the famous HLLC solver. However, the main difference relies on the fact that here the wave speeds are consistently derived and ordered. This positive and entropic Eulerian simple approximate Riemann solver is the cornerstone upon which we build a positive and entropic Godunov-type Finite Volume scheme provided an adhoc time step condition is fulfilled.

The robustness and the accuracy of our basic first-order Finite Volume schemes are assessed against various classical and demanding numerical tests: Sod and modified shock tube, extreme double rarefaction waves, Woodward-Collela blastwave and Le Blanc shock tube. We witness pertinent results, i.e., the positivity-preserving and entropy-stability properties are correctly enforced, and mesh convergence was also achieved. Furthermore, we note that the positivity-preserving property is mainly enforced during compression, thus, the trajectory of shock waves can be tracked along with the evolution of the waves speeds.

Due to its good properties, the first-order Eulerian scheme is then used in a last part as a parachute scheme for a high-order extension under MOOD paradigm. We have implemented high-order cell polynomial reconstruction in space, the Runge-Kutta method in time and the MOOD a posteriori limiting. In this part, we demonstrate that the high-order extension has similar positivity-preserving property as the firstorder numerical scheme while improving the numerical accuracy. The nominal accuracy is extended up to 5th order and may be decremented to the 1st order when needed. Two test cases with smooth solutions are first presented as a proof of rate of convergence and to validate the cell polynomial reconstruction. On the advanced test cases presenting complex wave interactions, the numerical results improve with the order of the reconstruction and the MOOD limiting ensures the robustness of the scheme.

In the future, we plan to implement a similar concept, meaning developing a cell-centered Eulerian positivity-preserving Riemann solver derived from the Lagrangian framework, in a multi-dimensional geometry on structured and unstructured mesh. We will study the ability of the numerical scheme to cure numerical shock instabilities such as the odd-even decoupling and carbuncle phenomenon [START_REF] Quirk | A contribution to the great riemann solver debate[END_REF], and others [START_REF] Toro | Anomalies of conservative methods: analysis, numerical evidence and possible cures[END_REF]. We also consider achieving multi-dimensional high-order extension with a MOOD limiting which may be an interesting challenge.

Figure 1 :

 1 Figure 1: Simple Riemann solver for Lagrangian gas dynamics.

Figure 2 : 5 )

 25 Figure 2: Plots of entropy production for a perfect gas equation of state (γ = 7 5) for z s = ρ s a s , z s = 2ρ s a s and z s = 4ρ s a s .

Figure 3 :

 3 Figure 3: Simple Riemann solver for Eulerian gas dynamics.

Figure 4 :

 4 Figure 4: Lagrangian approximate Riemann solvers at cell interfaces.

) 4 . 4 .

 44 Positivity of specific volume and internal energy If the time step condition (71) is taken for granted then through the Finite Volume scheme (70) τ n+1 i , u n+1 i and e n+1 i

Figure 5 :

 5 Figure 5: Sod shock tube problem -Numerical density of the 1st order Lagrangian scheme -Left : Mesh convergence from N = 100 to 400 cells of the density variable over the full domain -Right : Comparison of different methods to determine the impedances z l , z r : Method 1 (magenta diamond), Method 2 (green square) and Method 3 (blue circle).

Figure 6 :

 6 Figure 6: Woodward-Collela blastwave problem -1st order Lagrangian scheme using the first method to determine the impedances -Left : Mesh convergence from N = 100 to 400 cells of the numerical density over the full domain -Right : zoom on the central area

Figure 7 :

 7 Figure 7: Woodward-Collela blastwave problem -1st order Lagrangian scheme -Mesh N = 200 -Left : Comparison of different methods to determine the impedances z over the full domain such as Method 1 represented in magenta diamond, Method 2 represented in green square and Method 3 represented in blue circle -Right : zoom on the central area.

Figure 8 :

 8 Figure 8: Woodward-Collela blastwave problem -1st order Lagrangian scheme -Evolution of the normalized right impedance z r with each time step: (a) Method 1, (b) Method 2 and (c) Method 3.

  mass flux log(zr/(ρrar))

Figure 9 :

 9 Figure 9: Le Blanc shock tube problem -1st order Lagrangian numerical scheme with Method 1 (to determine the impedance z) for N = 900, 1800 and 3600 cells -Left: Mesh convergence of numerical density -Middle: Mesh convergence for the numerical specific internal energy -Right: Evolution of the normalized right impedances z r for each time step.

Figure 10 :

 10 Figure 10: Extreme double rarefaction problem -1st order Lagrangian scheme -Density (left), velocity (middle) and specific internal energy (right) -Mesh convergence from N = 100 to 400 cells.

  t

Figure 11 :

 11 Figure 11: Eulerian approximate Riemann solvers at cell interfaces.

Figure 12 :

 12 Figure12: Modified Sod problem -Eulerian numerical scheme -Density, velocity, and specific internal energy for the exact solution (black line) vs the numerical results (symbols) with mesh convergence such as N = 100 cells (red cross), 200 (blue circle) and 400 (green square).

Figure 13 :

 13 Figure 13: Woodward-Collela blastwave problem -Eulerian numerical scheme -Density for the reference solution (black line) vs the numerical results -Left: Mesh convergence of the density variable on N = 800 cells (red cross), 1600 (blue circle), 3200 (green square) and N = 6400 cells (magenta diamond) -Middle: Time step evolution -Right: Evolution of the normalized right impedance z r with each time step with N = 1600.

  mass flux log(zr/(Ωrar))

Figure 14 :

 14 Figure 14: Le Blanc shock tube problem -Eulerian numerical scheme -Left : Density plotted on the log scale with mesh convergence on N = 900 (red cross), N = 1800 (blue circle) and N = 3600 (green squares) cells -Middle: Time step evolution -Right: Evolution of the normalized right impedance z r with each time step with N = 1800.

Figure 15 :

 15 Figure 15: Extreme 123 problem -Eulerian numerical scheme -Density, velocity, and specific internal energy for the exact solution (black line) vs the numerical results (symbols) -Top panels: numerical results for N = 100 cells -Bottom panels: mesh convergence, N = 100 cells (red cross), 200 (blue circle) and 400 (green square).

Figure 16 :

 16 Figure 16: Test 2 to 5 Toro (top to bottom) -Eulerian numerical scheme N = 100 cells, CFL = 0.9 -Density, velocity, and specific internal energy (left to right) for the exact solution (black line) vs the numerical results (symbols).

Figure 17 :

 17 Figure17: Test 6 and 7 of Toro (left to right) -Eulerian numerical schemes N = 100 cells, CFL = 0.9 -Density for the current scheme (red crosses) against the classical HLL scheme (blue bullets) and the exact solution (black line).

Figure 18 :

 18 Figure18: Sketch of the high order numerical scheme based on polynomial reconstructions (blue) and a posteriori MOOD loop (green) embracing the general scheme (black).

Figure 19 :

 19 Figure 19: Density sine wave test case -Eulerian numerical scheme -L 2 error for the 1 st , 2 nd , 3 rd , 4 th and 5 th order schemes Numerical results for N = 50, N = 100, N = 200 and N = 400 cells.

Figure 20 :

 20 Figure 20: Density sine wave test case -Efficiency of the Eulerian numerical schemes -L 2 error versus CPU time for 1 st , 2 nd , 3 rd , 4 th and 5 th order schemes.

Figure 21 :

 21 Figure 21: Smooth solution test case -Third order (d = 3) Eulerian numerical scheme -Density, velocity, and specific internal energy for the exact solution (black line) vs the numerical results (symbols) -Numerical results for N = 100 cells (red cross) -The curves are almost superimposed.

Shu-

  Osher oscillatory test case. This test[START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF] is a 1D hydrodynamic shock tube. The downstream flow has a sinusoidal density fluctuation ρ = 1 -ε sin(λπx) with a wave length λ = 5 and an amplitude ε = 0.2. A Mach 3 shock front is initially located at x = -4 on domain X ∈ [-5; 5] and a CFL of 0.1. The left and the right states are as follows: ρ l = 3.857143 , u l = 2.629369 , p l = 10.33333 , ρ r = 1 + 0.2 sin(5πx) , u r = 0 , p r = 1 (96)

Figure 22 :Table 8 :

 228 Figure22: Modified smooth solution test case -Eulerian numerical scheme -L 2 error for the 1 st , 2 nd , 3 rd and 4 th order accurate schemes -Numerical results for N = 50, 100, 200, 400 and 800 cells.Percentage of cell reconstruction type(%)Comparison of CPU time (µs)d = 0 d = 1 d = 2 d = 3 d = 4 d = 0 d = 1 d = 2 d = 3 d = 4 N = 200 3.0 0.25 3.0 0.25 92.50 129 288 423 869 1254 N = 400 10.25 0 0.5 0.25 89.75 248 980 1716 1791 2292 N = 800 10.5 0 0.15 0.1 89.75 500 2130 3454 4079 4388 Table 8: Shu-Osher oscillatory test case -High-order extension of the Eulerian numerical scheme with 200, 400 and 800 cells respectively -Percentage of cells with 1 st , 2 nd , 3 rd and 4 th order cell reconstruction (Left) -Comparison of the CPU time with 1 st , 2 nd , 3 rd and 4 th order cell reconstruction for 200, 400 and 800 cells (Right).

Figure 23 :

 23 Figure 23: Shu-Osher oscillatory test case -High-order extension of the Eulerian numerical scheme with a mesh of N = 200 (top), N = 400 (middle) and N = 800 (bottom) -Evolution of the polynomial reconstruction degree d of each cell with each time step -Black cells are cells decremented to d = 0 , red cells are cells decremented to d = 1, blue cells are cells decremented to d = 2, green cells are cells decremented to d = 3 and white cells are cells that stayed at d = 4.

Table 1 :

 1 Initial conditions for the 1D test cases, recall ρ = 1 τ .

		state	Right state	Domain Discont. Final time EOS
		(ρ, u, p) l		(ρ, u, p) r	Ω	x 0	t final	γ
	Sod	(1, 0, 1)	(0.125, 0, 0.1)	[0, 1]	0.5	0.2	7/5
	Extreme Rarefact.	(1, -3.5, 0.4)	(1, 3.5, 0.4)	[0, 2]	1	0.15	7/5
	Le Blanc	(1, 0, 2/3 × 10 -1 ) (10 -3 , 0, 2/3 × 10 -10 )	[0, 9]	3	6.0	5/3
	Name	Left state	Middle state Right state Domain Discont. Final time EOS
		(ρ, u, p) l	(ρ, u, p) m	(ρ, u, p) r	Ω	x 0 , x 1	t final	γ
	Woodward-Collela (1, 0, 1000)	(1, 0, 0.01)	(1, 0, 100)	[0, 1]	x 0 = 0.5	0.038	7/5
						x 1 = 0.9

Table 2 :

 2 Initial conditions for the 1D test cases of Toro from[START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF].

		state	Right state	Discont. Final time
		(ρ, u, p) l	(ρ, u, p) r	x 0	t final
	2	(1, -2.0, 0.4)	(1, 2, 0.4)	0.5	0.15
	3	(1, 0, 1000)	(1, 0, 0.01)	0.5	0.012
	4	(5.99924, 19.5975, 460.894) (5.99242, -6.19633, 46.0950)	0.4	0.035
	5	(1, -19.59745, 1000)	(1, -19.59745, 0.01)	0.8	0.012
	6	(1.4, 0, 1)	(1, 0, 1)	0.5	2.0
	7	(1.4, 0.1, 1)	(1, 0.1, 1)	0.5	2.0

Table 4 :

 4 Butcher tableau for the Runge-Kutta explicit time discretization schemes and three examples used in this work.

Table 5 :

 5 L 1 , L 2 and L ∞ norm errors on density ρ between the numerical solution and the exact solution of a density since wave problem until t final = 1.

Table 6 :

 6 Density sine wave test case -Efficiency of the high-order Eulerian numerical schemes -Comparison of efficiency for a fixed error (horizontal line) and at a fixed CPU time (vertical line).

Table 7 :

 7 L 1 -, L 2and L ∞ -norm errors on density ρ between the numerical solution and the exact solution of the modified smooth solution problem until t final = 0.8.

Cell indexes

Stencil size Linear system size Degree d i -3 i -2 i -

Table 3: Centered stencils S d i used for the polynomial reconstructions of degree d, their size and the linear system size associated to the reconstruction operator.

u i is expressed under a Taylor basis as u i (x) = u i + d m=1 r i,m ψ i,m (x) where the unknowns are the polynomial coefficients r i,m and the polynomial basis functions {ψ i,m } m=1,...,d are such that for all m:

is the cell centroid. The polynomial coefficients are then optimized in the least-squares sense by the minimization of the cost function

which amounts to solve an over-determined linear system if the size of the stencil S d i is such that K i > d, and usually we take K i 1.5d. The size of the linear system is (2K i ) × d. The reconstruction operator is accurate at order (d + 1) for smooth enough solution. For the sake of simplicity the stencils S d i are always the centered ones and adapt their width as d changes, see table 3. In this work the primitive variables are reconstructed. Notice that conservative or characteristics variables could also be considered leading to slightly more or less oscillatory results. The reconstructed values are used in substitution to the mean values within the Riemann solvers as

Also remark that the reconstruction operates at time t n , and, once a maximal polynomial degree is set at d ≥ 0, then at maximum (d + 1) reconstructions are possible for degrees 0 ≤ k ≤ d. Obviously for d = 0 the reconstruction is trivially equal to the cell mean value.

a posteriori high-order limiting: detection, parachute scheme and MOOD loop

As is well known, the polynomial reconstruction operator suffers from the Runge phenomena in the case of non-smooth enough solution, leading to oscillatory polynomials u i as soon as d ≥ 1. In classical second-order FV schemes (d = 1) a so-called slope (or flux) limiting procedure is supplemented to avoid those spurious oscillations. Unfortunately, there exists no agreement on how to limit higher order polynomial reconstructions. In this work we rely on the a posteriori Multi-dimensional Optimal Order Detection (MOOD) algorithm [START_REF] Clain | A high-order finite volume method for systems of conservation laws -multi-dimensional optimal order detection (MOOD)[END_REF]. For this paradigm three entities are needed:

(i) a list of detection criteria D to assess the validity of a numerical solution. These criteria must ensure the physical admissibility of the numerical solution (i.e the positivity). Once the physical admissibility is assured then some numerical acceptability properties may be required. Often an essentially-non-oscillatory (ENO) behavior, and, possibly the computer representation validity (i.e no NaN) are enforced. Bad/troubled cells are gathered into set B.

(ii) a parachute scheme which must be a genuinely robust scheme producing a numerical solution fulfilling the detection criteria D. Here we consider the first order Eulerian scheme developed in this work. more, we confirm the mesh convergence and the accuracy convergence as the order of the polynomial reconstruction increases. At last in figure 26 we present the polynomial degree used in each cell at each time-step for N = 200 (left) and 800 cells (right). It can be observed that the main waves are followed by the lower polynomial degree (black and red) while the maximal accuracy (white) is maintained away from them as expected.

Percentage of cell reconstruction type (%) Comparison of CPU time (ms)

Then in table 9 we present the percentages of cells updated with each polynomial degree from 0 to 3 and for three different meshes. We observe that about 85 -93% of cells are updated with the highest order scheme, the remaining 7 -15% left are updated with the other schemes, mainly the parachute one. At last the CPU times show that the cost of second order scheme is about 5 times more expensive than a first order one (respectively 11 times and 16 times for third and fourth order ones).

8. Appendices: Positivity of specific volumes (τ l , τ r )

Let us study under which conditions on z l and z r , the specific volumes τ l and τ r remain positive assuming that τ l and τ r are strictly positive, recalling that

These explicit expressions of the specific volumes τ l and τ r have been obtained substituting the expression of u (35) into the left and right jump relations, respectively (33a) and (34a). To ease the study let us rename the unknowns as x = z l > 0, y = z r > 0 and the parameters are for s = l, r

Notice that the signs of a and b are a priori unspecified but τ l > 0 and τ r > 0, and, a l = a r , b l = b r if and only if τ l = τ r . Then (97) recasts into

Hence depending on the signs of a and b we may have different situations to consider, each L and R functions having different positive regions. Because L and R must be positive simultaneously, the intersection of these regions will be carried on as a last step. Notice that the terms x(x + y) ≥ 0 and y(x + y) ≥ 0 because (x, y) ∈ R + × R + , therefore the negative parts of L(x, y) and R(x, y) are necessarily resulting from the terms a s yb s . Let us first study L(x, y) as a function of the sign of a l and b l .

Case ++: a l > 0, b l > 0. In this case we have only b l which contributes for the negativity:

The hyperbole L intersects the x-axis at x = √ b l and the y-axis at y = b l a l , and to ensure the positivity it is sufficient that (x, y) belongs to the following positivity region

Case +-: a l > 0, b l < 0. L is always positive in this case, hence C +- L = (R 2 ) + .

Case -+: a l < 0, b l > 0. In this case L presents a vertical asymptote at location

. The positivity region is then

Case --: a l < 0, b l < 0. This last case corresponds to a branch of hyperbole intersecting the y-axis and with a vertical asymptote at x = -a l . Hence

The same study can be made for function R, the situation being equivalent inverting the role of x for y, and, b l for -b r , therefore the set C ++ R , C +- R , C -+ R and C -- R can be derived accordingly. Moreover the constraints on the acoustic impedances state that x ≥ ρ l a l and y ≥ ρ r a r . 50

At last gathering all constraints in the case of compression we can define the regions for (x, y) = (z l , z r ) ensuring the positivity of the specific volumes where β l = max( √ b l , -a l ) and β r = max( √ -b r , -a r ). As final remarks, let us point that these sufficient constrains could be more precisely derived but at expensive cost. More importantly, recall that these constrains only apply in a compression situation.