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

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

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

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

Entropy stable and positivity preserving Godunov-type schemes for multidimensional hyperbolic systems on unstructured grid

Keywords: Hyperbolic system of conservation laws, Godunov-type scheme, Simple approximate Riemann solver, Entropy stability, Lagrangian representation, Eulerian representation, Gas dynamics

This paper describes a novel subface flux-based Finite Volume (FV) method for discretizing multidimensional hyperbolic systems of conservation laws of general unstructured grids. The subface flux numerical approximation relies on the notion of simple Eulerian Riemann solver introduced in the seminal work [G. Gallice; Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates; Numer. Math., 94, 2003]. The Eulerian Riemann solver is constructed from its Lagrangian counterpart by means of the Lagrange-to-Euler mapping. This systematic procedure ensures the transfer of good properties such as positivity preservation and entropy stability. In this framework, the conservativity and the entropy stability are no more locally face-based but result respectively from a node-based vectorial equation and a scalar inequation. The corresponding multi-dimensional FV scheme is characterized by an explicit time step condition ensuring positivity preservation and entropy stability. The application to gas dynamics provides an original multi-dimensional conservative and entropy-stable FV scheme wherein the numerical fluxes are computed through a nodal solver which is similar to the one designed for Lagrangian hydrodynamics. The robustness and the accuracy of this novel FV scheme are assessed through various numerical tests. We observe its insensitivity to the numerical pathologies that plague classical face-based contact discontinuity preserving FV formulations.

Introduction

The conservation laws of gas dynamics, magneto-hydrodynamic flows and other branches of classical physics are typically expressed by nonlinear hyperbolic systems of Partial differential equations (PDEs) that prepare mathematically a host of wave interactions occurrences. The solutions of nonlinear hyperbolic systems generate jump discontinuities that propagate on as shock waves for instance. The mathematical theory may be challenging since it is necessary to confront weak solutions. In consequence, numerical simulation is often employed to tackle hyperbolic problems. A Riemann Solver (RS), should it be exact or approximated, is one of the essential tool in Computational Fluid Dynamics (CFD). In all generality, it is a method for computing the evolution of a discontinuity separating two constant physical states, namely a Riemann problem. S.K. Godunov [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] proposed the first Finite-Volume (FV) scheme to simulate compressible gas dynamics equations employing a Riemann solver that makes use of the exact solution of a Riemann problem. This first-order accurate method gained popularity and is now known as the Godunov method/scheme in CFD. However, this exact Riemann solver may become prohibitively expensive when the Riemann problem becomes too complex, leading to the subsequent development of alternative approximate and less expensive Riemann Solvers following on Godunov approach -to cite but a few, Roe [START_REF] Roe | Approximate Riemann solvers, parameter vectors and difference scheme[END_REF], Harten, Lax and van Leer (HLL) [START_REF] Harten | On upstream Differencing and Godunov-Type schemes for Hyperbolic Conservation Laws[END_REF], Einfeldt (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 HLL-Riemann solver[END_REF], Engquist [START_REF] Engquist | One-sided difference approximations for nonlinear conservation laws[END_REF], Osher [START_REF] Osher | Riemann solvers, the entropy condition and difference approximations[END_REF] and others and we refer the interested readers to the book of Toro [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] and Godlewski-Raviart [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF].

The aforementioned approximate Riemann solvers are intrinsically one dimensional in their framework and in spite of this rather good success of approximate schemes, they are presumably less efficient in multidimensional problems due to the directional bias. The standard way of extending to multidimensions with the one-dimensional schemes is through flux splitting algorithms by solving one dimensional Riemann problems in a direction normal to each cell interface for all faces independently to achieve multidimensional behavior. This notion neglects flow variations that might be propagating in a transversal direction. Therefore, researchers have been on an endeavor to develop multidimensional Riemann solvers to solve Riemann problems in different directions of the interface, hoping to better capture flow characteristics. Among other attempts in contributing to multidimensional solvers is the generic multidimensional HLLE solver of Balsara [START_REF] Balsara | Self-Similar internal Structure-Part 2-Application to Hyperbolic Conservation Laws on Unstructured Meshes[END_REF]. The base idea is to use the one-dimensional Riemann solver in all directions to coupled with multidimensional correction terms to constitute a multidimensional solver for conservative hyperbolic systems. The first truly multidimensional solver by Roe [START_REF] Deconinck | A multidimensional generalization of Roe's flux difference splitter for the Euler equations[END_REF] also brought compelling contribution to this topic. The key element of this solver is to assume continuous piecewise linear space variation of the variables defined at cell vertices, allowing for a multidimensional generalization without employing 1D Riemann problems or dimensional splitting.

By the same token, the dimensionally split approximate Riemann solvers also give rise to the occurrence of various forms of numerical shock instabilities when simulating shock wave propagation, i.e., odd-even decoupling and carbuncle, refer to [START_REF] Quirk | A contribution to the great Riemann solver debate[END_REF]. Such instabilities appear as serrated disturbance of a flow featured as a nonphysical default. The earliest report on numerical instabilities by Quirk [START_REF] Quirk | A contribution to the great Riemann solver debate[END_REF] states that low-dissipative schemes causing expansive growth of acoustic waves may be the root of the carbuncle phenomenon. Thus, there is a strong need for really multidimensional formulations that resolve all characteristic fields while instituting necessary dissipation. For instance, Rodionov [START_REF] Rodionov | Artificial viscosity Godunov-type schemes to cure the carbuncle phenomenon[END_REF] constructs a multipurpose remedy by adding artificial viscosity. Alternatively, Fleischmann [START_REF] Fleischmann | A shock-stable modification of the HLLC riemann solver with reduced numerical dissipation[END_REF] also proposes a modification of the HLLC Riemann solver with a centralized reformulation of the numerical flux that reduces the acoustic dissipation, offering a shock-stable solver. In this current work, we will be bringing into play a node-based conservation condition by involving all materials surrounding a node. Upon that, transversal information are not omitted and this strategy substantiates the pre-existing assumption that multidimensional dissipation takes effect on shock instabilities problem.

One common feature of most of the present solvers is that they were developed in the Eulerian framework point of view. Contrastingly, it is possible to adopt a different approach by tackling the resolution of the Riemann problem from the Lagrangian standpoint such as in the 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]. The main feature of Lagrangian numerical methods lies in the fact that the motion of the fluid is intrinsically linked to the geometrical transformation that follows the fluid path, providing a natural framework to track interfaces of multi-dimensional flows. The simplicity of the Lagrangian formulation not only allows for a direct estimation of waves speeds in an ordered manner but it also facilitates the study of particular properties (conservation, positivity-preserving and entropy stable). Recurrently, numerical approximations may generate negative density or pressure which then leads to instability or code crash. This event is more critical in the Lagrangian framework on account of to the moving and deforming grid during computation. The positivity-preserving property can be assessed easily in the Lagrangian framework. Once a Lagrangian Riemann solver is built, the Eulerian counter-part that inherits the properties of the Lagrangian one can be obtained [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 Nonconservatifs[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][START_REF] Chan | Positivity preserving and entropy consistent approximate Riemann solvers dedicated to the high-order MOOD-based Finite Volume discretization of Lagrangian and Eulerian gas dynamics[END_REF]. Besides Gallice, only few authors have been working adopting the resolution from a Lagrangian standpoint, for instance, Cheng and Shu [START_REF] Cheng | Positivity-preserving Lagrangian scheme for multi-material compressible scheme[END_REF] developed positivity-preserving HLLC approximate Riemann solver for the Lagrangian scheme in one and two dimensional for compressible Euler equations with general equations of state. 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 authors also worked on positivity-preserving approximate Riemann solvers conducted for both ideal gas and non ideal gas equations of state and extending to high-order accuracy with appropriate limitation.

As for the multidimensional Lagrangian case, the one-dimensional solver cannot be implemented directly due to the vast number of neighboring cells sharing a node. Hence, in this framework, the numerical fluxes are evaluated by means of an approximate Riemann solver located at the grid nodes that provides the nodal velocity required to move the grid in a compatible manner. Two main approaches are classically used, namely the staggered scheme and cell-centered FV schemes. This paper focuses on the cell-centered approach. Different techniques may be employed to build the numerical fluxes and move the grid through the use of approximate Riemann solvers with respect to the Geometrical Conservation Law. In [START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF] and [START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF], truly multidimensional Lagrangian schemes that are arbitrary for unstructured mesh are proposed. In this approach, for the twodimensional case, the FV scheme is established using subface fluxes that are expressed in terms of the difference between the velocity in the cell and the nodal velocity, instigating non-conservativity of the scheme. Consequently, the conservation of the scheme is retrieved via a nodal conservation condition stating that the sum of the subface fluxes impacting a node must be equal to zero. This gives birth to a nodal solver, the cornerstone of Lagrangian schemes in order to evaluate the velocity of displacement of the nodes, thus allowing the movement of the grid. This new approach endorsed the construction of multidimensional FV schemes for Lagrangian gas dynamics that is more robust than the classical ones.

Lately, Shen and his collaborators [START_REF] Shen | A robust and contact resolving Riemann solver on unstructured mesh, Part I, Euler method[END_REF] were inspired by the aforementioned FV schemes for Lagrangian gas dynamics and they adapted it by extending to the case of Eulerian gas dynamics. They employed a HLLC-like approximate Eulerian Riemann solver using subface fluxes of which the approximation is based on [START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF] nodal solver, making it possible to ensure the conservation of the scheme. However, the construction of their scheme does not guarantee wave speeds ordering and lacks an in-depth study of important properties such as positive density and internal energy and entropy inequality. Nevertheless, the numerical results proved to be promising with this method and encouraged us to carry out an extensive work not only to correct the previous flaws, but also to propose a generalization of formulation of multidimensional FV scheme for all types of hyperbolic systems.

In the present work, we shall describe a peculiar multidimensional subface flux-based FV scheme for which the subface flux results from an approximate Riemann solver. With reasonable assumptions on this solver, we are able to exhibit a time step condition ensuring that the cell-centered solution at time t n+1 consists of a convex combination between the solution at time t n and the intermediate states of the Riemann solver at each subfaces. An entropy inequality as well as an entropy flux can be identified. The conservation and the entropy control of this FV formulation are obtained by means of node-based conservation conditions that take into account the states and geometry of the cells surrounding a node.

The construction of approximate Riemann solvers characterized by good properties regarding positivity and entropy stems from the seminal work [START_REF] Gallice | Positive and Entropy Stable Godunov-Type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF]. More precisely, we impose those properties onto the approximate Lagrangian Riemann solver and derive its Eulerian counter part using the Lagrange-to-Euler mapping. The node-based conservation and entropy conditions for a simple multidimensional Riemann solver are then put forward once again by calling upon the Lagrangian formulation. This original framework is incorporated into the system of equations for gas dynamics to put into effect the construction, bringing about a first-order multidimensional Eulerian FV scheme that is provably preserving the physical admissible states and impervious to numerical instabilities. Numerical tests will be carried out to illustrate and validate this original FV scheme.

Therefore, the aim of this current work is to formulate a generic multidimensional FV scheme for hyperbolic systems on unstructured grids based on subface fluxes that guarantees good properties. We will be recalling the Lagrange-to-Euler mapping that allows us to develop a Lagrangian Riemann solver and the Eulerian counterpart, both sharing interesting properties. On that, a firstorder Godunov-type FV scheme which is positivity preserving and entropy stable with an explicit CFL like condition is developed. Then, the Lagrangian nodal solver is convened in order to fulfill the node-based conservation and entropy condition in the original FV scheme with the intention of working as a cure for numerical instabilities. We want assess the accuracy and the robustness of the scheme on different test cases applied to the gas dynamics system.

The rest of this paper is organized as follows. After this introduction, we present a multidimensional Godunov-type scheme for solving hyperbolic systems of conservation laws on unstructured grids. A generic FV scheme characterized by a corner flux associated to an approximate Riemann solver is developed and the conservation condition and its entropy stability is studied. The third section introduces the construction of the one-dimensional approximate Riemann solver by means of the Lagrange-to-Euler mapping. More precisely, the notion of simple approximate Riemann solver in the Lagrangian framework is derived along with some properties, and the Eulerian version is proposed by making use of the Lagrangian one as a building block. The fourth section revisits the Lagrangian node-based conservation and entropy conditions in the case of simple Riemann solvers. To complete the design, the node-based conservation condition is applied on the gas dynamics system discretized onto unstructured grids, calling upon the Lagrangian nodal solver. A test campaign to assess the performances of the scheme, especially to study instability behaviors of the newly developed scheme is provided. Concluding remarks and perspectives are drawn in the last section.

Multi-dimensional Godunov-type Finite Volume scheme for solving hyperbolic systems of conservation laws on unstructured grids

Governing equations and notation

We aim at designing multi-dimensional Finite Volume schemes for solving the hyperbolic system of conservation laws

∂U ∂t + ∇ • F(U) = 0. (1)
Here, U = U(x, t), for x ∈ R d and t ≥ 0, is the vector of conservative variables which takes values in R q and F = F(U) is the flux tensor in R q ×R d . The positive integers d and q denote respectively the space dimension and the size of the foregoing hyperbolic system. Let e j , for j = 1 . . . d, be the the j-th vector of the Cartesian basis of R d then F j = Fe j ∈ R q is the j-th component of the tensor flux and its divergence writes We assume that system (1) is equipped with the entropy-entropy flux pair (Σ, Q). Namely, U � → Σ(U) is convex and the following compatibility condition holds true for the entropy flux [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF] � ∂F j ∂U

∇ • F(U) = d � j=1 ∂F j ∂x j . x c x p l pcf ω pc ω c F pcf x - p x + p n pcf f
� t ∂Σ ∂U = ∂Q j ∂U , for j = 1 . . . d,
where

Q j = Q • e j .
We are looking for the entropic solutions of (1), that is the ones satisfying the following entropy inequality

∂Σ ∂t + ∇ • Q(U) ≤ 0, (2)
which turns into an equality for smooth solutions. Let D ⊂ R q be the domain of definition of F(U), Σ(U) and Q(U). We assume that D is a convex subset of R q .

For the sake of simplicity, we shall limit the presentation of the numerical methods to the twodimensional case, i.e., d = 2, knowing that the three-dimensional extension is quite straightforward. The computational domain is a polygonal portion of R 2 and we pave it with a collection of non overlapping polygonal cells ω c where c is the generic label of the cell. Let P(c) be the set of vertices (points) of ω c . The generic label of a point is p and x p denotes its vector position. In the counterclockwise ordered list of points of ω c , p -and p + are respectively the previous and the next points with respect to p, refer to figure 1. The subcell ω pc related to cell c and point p is the quadrangle formed by joining the cell centroid, x c , to the midpoints of [x p -, x p ], [x p , x p + ] and to x p . The set of subcells ω pc for p ∈ P(c) constitutes a partition of the cell ω c , that is,

ω c = � p∈P(c) ω pc .
We also introduce the set of faces of cell ω c and denote it F(c), for instance [x p , x p + ] belongs to F(c). Each face f of cell c is decomposed into subfaces by means of the partition of c induced by the subcells pc for p ∈ P(c). This leads us to define SF(pc) the set of subfaces attached to the corner pc, which is nothing but the set of faces of subcell ω pc impinging at point p, for instance [x p , 1 2 (x p + x p + )] belongs to SF(pc). We denote respectively by l pcf and n pcf the measure and the unit outward normal of the subface f . We observe that the set of subfaces SF(pc) for p ∈ P(c) constitutes a partition of the set of faces of ω c , that is,

F(c) = � p∈P(c)
SF(pc).

Subface-based Finite Volume discretization

Integrating the system of conservation laws (1) over ω c and employing Green formula leads to

|ω c | dU c dt + � ∂ωc F(U)n ds = 0, (3)
where

U c (t) = 1 |ωc| � ωc U(x, t) dv is the cell-averaged value of U over ω c . Employing a classical
first-order explicit time integration turns (3) into

U n+1 c -U n c + Δt |ω c | � ∂ωc F(U n )n ds = 0. (4)
Here, U n c denotes the approximation of U c (t) at time t n , and, t n+1 = t n + Δt where Δt > 0 is the time step.

The design of the Finite Volume scheme (4) requires to construct an approximation of the normal flux integral. In what follows, we are going to define an original node-based approximation of this integral term relying on the partition of ω c into subcells, that is,

� ∂ωc F(U n )n ds = � p∈P(c) � ∂ωpc∩∂ωc F(U n )n ds. ( 5
)
The surface integral term at the right-hand side of ( 5) is approximated along the subfaces as follows

� ∂ωpc∩∂ωc F(U n )n ds = � f ∈SF (pc) l pcf F pcf ,
where F pcf is the subface flux related to the subface f attached to the corner pc. Substituting this subface-based approximation of the flux into (4) yields the subface-based generic Finite Volume scheme

U n+1 c -U n c + Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf F pcf = 0, (6)
which is characterized by the subface flux F pcf displayed in figure 1 by the blue rectangles. Viewed from cell ω c , this peculiar Finite Volume discretization introduces two subface fluxes per cell face. We observe that this type of Finite Volume discretization which consists in splitting the faces into subfaces has been already utilized not only in the framework of cell-centered Lagrangian hydrodynamics [START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF][START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF][START_REF] Loubère | Handbook of Numerical Methods for Hyperbolic Problems: Basic and Fundamental Issues[END_REF] but also for developing cell-centered diffusion schemes [START_REF] Maire | A nominally second-order accurate finite volume cell-centered scheme for anisotropic diffusion on two-dimensional unstructured grids[END_REF][START_REF] Jacq | A Nominally Second-Order Cell-Centered Finite Volume Scheme for Simulating Three-Dimensional Anisotropic Diffusion Equations on Unstructured Grids[END_REF]. We note in passing that this formalism encompasses the classical face-based Finite Volume discretization [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]. It remains to provide a consistent numerical approximation of the subface flux.

Subface flux approximation by means of a Riemann solver

Let f be the generic subface attached to cell c and vertex p, characterized by its unit outward normal n pcf . We assume that the subface flux F pcf attached to the subcell f depends respectively on the adjacent cell averaged values U c , U d(c,f ) , where d(c, f ) is the neighbor of cell c such that f ⊂ (ω c ∩ ω d ), on the unit normal n pcf and also on the vectorial parameter v p ∈ R 2 related to node p. This leads us to write

F pcf = F pcf (U c , U d(c,f ) , n pcf , v p ). ( 7
)
Contrary to the classical face-based Finite Volume discretization, the foregoing subface flux expression exhibits a dependency on the nodal vectorial parameter v p , which is unknown for the moment. Such dependency has been already employed for designing flux approximation dedicated to the cell-centered discretization of Lagrangian hydrodynamics [START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF] and also more recently, inspired by the latter works, in the framework of cell-centered Eulerian hydrodynamics [START_REF] Shen | A robust and contact resolving Riemann solver on unstructured mesh, Part I, Euler method[END_REF]. In both cases, the parameter vector, v p , corresponds to a nodal approximation of the velocity field. For Lagrangian hydrodynamics discretization, the nodal velocity approximation is required to move the computational grid, whereas for Eulerian hydrodynamics its need is less obvious. We shall see later the fundamental role played by this nodal vector parameter. Bearing this in mind, the subface flux numerical approximation is constructed via the introduction of an approximate Riemann solver, which is nothing but the approximate solution of the one-dimensional Riemann problem defined in the n pcf direction (RP)

           ∂U ∂t + ∂[F n pcf (U)] ∂x n pcf = 0, U(x n pcf , 0) = � U c if x n pcf < 0, U d(c,f ) if x n pcf ≥ 0.
Here, x n pcf = x • n pcf is the coordinate in the unit normal direction, F n pcf (U) = F(U)n pcf is the projection of the tensor flux onto the unit normal direction. The resulting one-dimensional approximate Riemann solver depends on the states on both sides of the interface, on the self-similar variable ξ = xn pcf t and also on the parameter v p . With these arguments, the Riemann solver writes

W pcf = W pcf (U c , U d(c,f ) , n pcf , ξ, v p ). (8)
For all U l (left state), U r (right state), n (unit normal), ξ (self-similar variable) and v (vector parameter), we assume that the Riemann solver satisfies the following classical properties

• W pcf (U l , U r , n, ξ, v) = U l for -ξ large enough; • W pcf (U l , U r , n, ξ, v) = U r for ξ large enough; • W pcf (U, U, n, ξ, v) = U.
We also make the assumption that the approximate Riemann solver is symmetric with respect to the interface, that is

W pcf (U c , U d(c,f ) , n pcf , ξ, v p ) = W pdf (U d(c,f ) , U c , n pcf , -ξ, v p ). (9)
Finally, we express the subface flux F pcf in terms of the approximate Riemann solver W pcf as follows This subface flux approximation is obtained integrating the conservation law of (RP) over the space-time domain [-Δx l , 0]×[0, Δt] where Δx l and Δt are respectively space and time increments. This formula might be found in [START_REF] Harten | On upstream Differencing and Godunov-Type schemes for Hyperbolic Conservation Laws[END_REF].

F pcf = F(U c )n pcf - � 0 -∞ � W pcf (U c , U d(c,f ) , n pcf , ξ, v p ) -U c � dξ. ( 10
) npcf xp ωc ωd lpcf f Fpdf Fpcf
Remark 1. This choice to define the subface flux from the approximate Riemann solver is original and crucial. We note in passing that the aforementioned works employing subface flux-based discretization do not rely on such a definition. Here, the subface flux (10) might be viewed as a left-sided flux and there is absolutely no reason why F pcf = F pdf , where F pdf is the subface flux attached to subface f viewed from cell d = d(c, f ), refer to figure 2. Indeed, in general we should have F pcf � F pdf this implies that the Finite Volume scheme [START_REF] Einfeldt | On Godunov-type methods near low densities[END_REF] characterized by subface flux [START_REF] Einfeldt | On Godunov-type methods for gas dynamics[END_REF] is not conservative in the classical sense. We shall present in section 2.7 a novel framework to study the conservativity of this subface-based Finite Volume method.

Preservation of the definition domain

The notion of invariant domain is classical in the context of hyperbolic systems of conservation laws, refer for instance to [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws[END_REF]. A domain is invariant if for any initial condition U 0 belonging to the domain, the solution of the hyperbolic system under consideration remains in the domain for all time t > 0. The most practical situation corresponds to the case for which the invariant domain is convex. For instance, in the case of Lagrangian gas dynamics U = (τ, v, e) t , where τ is the specific volume, v is the velocity vector and e the total energy, and the definition domain (admissible set) D = {U such that τ ≥ 0 and e -1 2 v 2 ≥ 0} is convex. In what follows, we assume that the definition domain D of F(U), Σ(U) and Q(U) is convex and we shall study under which condition our Finite Volume discretization preserves the definition domain.

The studied Finite Volume scheme [START_REF] Deconinck | A multidimensional generalization of Roe's flux difference splitter for the Euler equations[END_REF] writes under the form

U n+1 c -U n c + Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf � F pcf -F(U n c )n pcf � = 0, (11)
where we have made use of the geometric identity

� p∈P(c) � f ∈SF (pc) l pcf n pcf = 0, (12)
to make appear the fluctuations related to each subface.

Assuming that the approximate Riemann solver is D-preserving, that is, if U n c ∈ D then W pcf (ξ) ∈ D for all ξ ∈ R, we shall determine the time step condition ensuring that the foregoing Finite Volume scheme is itself D-preserving, that is, U n+1 c ∈ D. To this end, we introduce ξ min pcf ≥ 0 such that W pcf (ξ) = U n c , for ξ < -ξ min pcf . This allows us to develop the subface flux expression as follows

F pcf =F(U n c )n pcf - � 0 -ξ min pcf (W pcf (ξ) -U n c ) dξ = F(U n c )n pcf + ξ min pcf U n c - � 0 -ξ min pcf W pcf (ξ) dξ.
Substituting the subface flux into the foregoing Finite Volume scheme yields

U n+1 c = U n c - Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf � ξ min pcf U n c - � 0 -ξ min pcf W pcf (ξ) dξ � .
Now, collecting the terms in factor of U n c we arrive at

U n+1 c =   1 - Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf ξ min pcf   U n c + Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf � 0 -ξ min pcf W pcf (ξ) dξ. (13) Assuming that U n c ∈ D and W pcf (ξ) ∈ D then U n+1 c ∈ D provided that the time step satisfies the condition Δt ≤ |ω c | � p∈P(c) � f ∈SF (pc) l pcf ξ min pcf . ( 14
)
In this case, we observe that U n+1 c is nothing but a convex combination of U n c and the intermediate states of the subface-based approximate Riemann solvers. Introducing ) satisfies

Δt c = |ω c | � p∈P(c) � f ∈SF (pc) l pcf ξ min
Σ n+1 c ≤   1 - Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf ξ min pcf   Σ n c (16) + Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf ξ min pcf Σ � 1 ξ min pcf � 0 -ξ min pcf W pcf (ξ) dξ � .
By virtue of Jensen inequality

Σ � 1 ξ min pcf � 0 -ξ min pcf W pcf (ξ) dξ � ≤ 1 ξ min pcf � 0 -ξ min pcf Σ(W pcf (ξ)) dξ.
Substituting the foregoing result into [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Nonconservatifs[END_REF], we arrive at

Σ n+1 c -Σ n c ≤ Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf � 0 -∞ (Σ(W pcf (ξ)) -Σ n c ) dξ. ( 17
)
Here, we have used the fact that if ξ ≤ ξ min pcf then W pcf (ξ) = U n c . Finally, utilizing the geometric identity [START_REF] Engquist | One-sided difference approximations for nonlinear conservation laws[END_REF] we get

� p∈P(c) � f ∈SF (pc) l pcf Q(U n c ) • n pcf = 0.
Introducing the foregoing expression into (17) leads to

Σ n+1 c -Σ n c + Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf � Q(U n c ) • n pcf - � 0 -∞ (Σ(W pcf (ξ)) -Σ(U n c )) dξ � ≤ 0. ( 18
)
This is formally the discrete counterpart of the continuous entropy inequality [START_REF] Balsara | Self-Similar internal Structure-Part 2-Application to Hyperbolic Conservation Laws on Unstructured Meshes[END_REF]. The comparison between the foregoing inequality and (2) incites us to define the subface entropy flux

Q pcf = Q(U n c ) • n pcf - � 0 -∞ (Σ(W pcf (ξ)) -Σ(U n c )) dξ. (19)
With this notation, we rewrite inequality [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and postprocessing facilities[END_REF] under the compact form

Σ n+1 c -Σ n c + Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf Q pcf ≤ 0. ( 20
)
Remark 2. At this stage nothing can be said regarding the entropy stability of the Finite Volume scheme under consideration. This crucial point shall be investigated in section 2.8.

Summary about the subface-based Finite Volume scheme

We have designed a subface-based Finite Volume scheme for which the subface flux F pcf is defined from the approximate Riemann solver W pcf . Under an explicit time step condition this Finite Volume scheme is D-preserving and satisfies a formal entropy inequality. We summarize hereafter the main characteristics of this Finite Volume scheme.

• The generic Finite Volume scheme reads

U n+1 c -U n c + Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf F pcf = 0,
where the subface flux expression in terms of the approximate Riemann solver W pcf writes

F pcf = F(U n c )n pcf - � 0 -∞ (W pcf (ξ) -U n c ) dξ.
• The time step condition to ensure that the FV scheme is D-preserving reads

Δt ≤ Δt c = |ω c | � p∈P(c) � f ∈SF (pc) l pcf ξ min pcf .
• The formal entropy inequality attached to the FV scheme under the foregoing time step condition reads

Σ n+1 c -Σ n c + Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf Q pcf ≤ 0,
where the corner entropy flux writes

Q pcf = Q(U n c ) • n pcf - � 0 -∞ (Σ(W pcf (ξ)) -Σ(U n c )) dξ.
Now, it remains to investigate not only the conservation property of the subface-based Finite Volume scheme but also its entropy stability.

Conservation property

This section aims at determining under which conditions the studied subface-based Finite Volume scheme is conservative. Assuming that the computational domain is the whole space R 2 , the subface-based Finite Volume scheme,

|ω c |(U n+1 c -U n c ) + Δt � p∈P(c) � f ∈SF (pc) l pcf F pcf = 0, is conservative if and only if � c |ω c |U n+1 c = � c |ω c |U n c ⇐⇒ � c � p∈P(c) � f ∈SF (pc) l pcf F pcf = 0.
Now, exchanging the summation over the cells with the summation over the nodes, the right-sided condition turns into

� p � c∈C(p) � f ∈SF (pc) l pcf F pcf = 0,
where C(p) is the set of cells sharing the point p.

We claim that a sufficient condition to ensure the conservativity of the subface-based Finite Volume scheme writes

� c∈C(p) � f ∈SF (pc) l pcf F pcf = 0. ( 21
)
This means that the summation over the cells c sharing point p of the fluxes attached to the subfaces impinging at p is equal to zero, refer to figure 3a, where the subface fluxes have been displayed by means of blue patches on both sides of each subface emanating from point p. Noticing that the sum over the cells c sharing p of the fluxes attached to the subfaces impinging at p is rigorously equal to the sum over the left-sided and the right-sided fluxes attached to the subfaces impinging at p leads to reformulate sufficient condition [START_REF] Godunov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF] into Here, SF (p) denotes the set of subfaces impinging at point p. For f ∈ SF(p), l pf is the length of subface f and n pf is its unit normal pointing towards the right state. In the foregoing equation,

� f ∈SF (p) l pf (F l pf + F r pf ) = 0. ( 22
)
n pcf x p l pcf f F pcf ω c (a) Cell-based notation. n pf x p f F l pf l pf F r pf (b) Face-based notation.
F l pf (resp. F r pf
) denotes respectively the left-sided (resp. right-sided) flux attached to the subface f , refer to figure 3b. By virtue of [START_REF] Després | Lagrangian systems of conservation laws[END_REF], with an obvious notation adaptation, the left and right-sided fluxes expressions in terms of the left state U lf , right state U rf , unit normal n pf and vector nodal parameter v p read

F l pf = F l pf (U lf , U rf , n pf , v p ), and F r pf = F r pf (U rf , U lf , -n pf , v p ).
Substituting this into the conservation condition [START_REF] Harten | A Random Choice Finite Difference Scheme for Hyperbolic Conservation Laws[END_REF] leads to

� f ∈SF (p) l pf � F l pf (U lf , U rf , n pf , v p ) + F r pf (U rf , U lf , -n pf , v p ) � = 0. (23)
It is worth pointing out that the subface flux on both sides of the subface f depends on the nodal vector parameter v p and thus a priori

F l pf (U lf , U rf , n pf , v p ) + F r pf (U rf , U lf , -n pf , v p ) � 0.
This prevents the node-based conservation condition [START_REF] Harten | On upstream Differencing and Godunov-Type schemes for Hyperbolic Conservation Laws[END_REF] from boiling down to the classical facebased conservation condition [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]. The node-based conservation condition consists of a system of scalar equations at node p involving not only the length of the edges impinging at this node but also the states adjacent to these edges. This system shall allow us to determine the parameter vector v p and we note that the dimension of v p is the rank of the aforementioned system. By analogy with the seminal works undertaken for cell-centered Lagrangian hydrodynamics [START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF][START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF], system [START_REF] Harten | On upstream Differencing and Godunov-Type schemes for Hyperbolic Conservation Laws[END_REF] shall be named the nodal solver for v p .

Thanks to [START_REF] Einfeldt | On Godunov-type methods for gas dynamics[END_REF] and with an obvious notation adaptation, we express the subface fluxes in (23) in terms of the approximate Riemann solver

F l pf (U lf , U rf , n pf , v p ) = F(U lf )n pf - � 0 -∞ (W l pf (U lf , U rf , ξ, v p ) -U lf ) dξ, F r pf (U rf , U lf , -n pf , v p ) = F(U rf )(-n pf ) - � 0 -∞ (W r pf (U rf , U lf , ξ, v p ) -U rf ) dξ.
Substituting these expressions of the corner fluxes into the nodal conservation condition [START_REF] Harten | On upstream Differencing and Godunov-Type schemes for Hyperbolic Conservation Laws[END_REF] and by virtue of the symmetry assumption of the approximate Riemann solver (9 24)

), i.e., W l pf (U lf , U rf , ξ, v p ) = W r pf (U rf , U lf , -ξ, v p ) we arrive at � f ∈SF (p) l pf � � 0 -∞ (W l pf (U lf , U rf , ξ, v p ) -U lf ) dξ + � ∞ 0 (W l pf (U lf , U rf , ξ, v p ) -U rf ) dξ (
+ (F(U rf ) -F(U lf )) n pf � = 0.
The Finite Volume scheme under consideration is conservative provided that the foregoing nodebased condition is fulfilled for each node.

Entropy stability

We aim at exhibiting conditions which ensure that the total entropy over the whole space R 2 is non increasing that is

� c |ω c | � Σ n+1 c -Σ n c � ≤ 0. ( 25
)
Knowing that the local in-cell inequality (20) holds true, i.e.,

Σ n+1 c -Σ n c + Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf Q pcf ≤ 0,
implies that the global entropy inequality [START_REF] Kamm | On efficient generation of numerically robust Sedov solutions[END_REF] holds true provided that

� c � p∈P(c) � f ∈SF (p) l pcf Q pcf ≥ 0,
where Q pcf denotes the subface entropy flux attached to the subface f of corner pc. Similarly to the study of the Finite Volume scheme conservation in section 2.7, we exchange the summation over the cells with the summation over the nodes in the foregoing inequality to arrive at

� p � c∈C(p) � f ∈SF (p) l pcf Q pcf ≥ 0,
where C(p) is the set of cells sharing the point p. Therefore, a sufficient condition to ensure that the Finite Volume scheme satisfies the global entropy inequality [START_REF] Kamm | On efficient generation of numerically robust Sedov solutions[END_REF] writes

� c∈C(p) � f ∈SF (p) l pcf Q pcf ≥ 0. ( 26
)
Once more, observing that the sum over the cells c sharing p of the entropy fluxes attached to the subfaces impinging at p is rigorously equal to the sum over the left-sided and the right-sided entropy fluxes attached to the subfaces impinging at p leads to reformulate sufficient condition [START_REF] Kulikowskii | Mathematical Aspects of Numerical Solution of Hyperbolic Systems[END_REF] into

� f ∈SF (p) l pf � Q l pf + Q r pf � ≥ 0, (27)
where SF(p) is the set of subfaces impinging at point p. Developing the expression of the left-sided and the right-sided subface entropy fluxes into the sufficient condition [START_REF] Loubère | Handbook of Numerical Methods for Hyperbolic Problems: Basic and Fundamental Issues[END_REF] yields

� f ∈SF (p) l pf � Q l pf (U lf , U rf , n pf , v p ) + Q r pf (U rf , U lf , -n pf , v p ) � ≥ 0. ( 28
)
Substituting the expression of the subface entropy flux [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF] in terms of the approximate Riemann solver into the nodal entropy condition [START_REF] Loubère | ReALE: A reconnection-based arbitrary-Lagrangian-Eulerian method[END_REF] and by virtue of the symmetry assumption of the approximate Riemann solver (9

), i.e., W l pf (U lf , U rf , ξ, v p ) = W r pf (U rf , U lf , -ξ, v p ) we arrive at � f ∈SF (p) l pf � � 0 -∞ (Σ(W l pf (ξ, v p )) -Σ lf ) dξ + � ∞ 0 (Σ(W l pf (ξ, v p )) -Σ rf ) dξ (29)
+ (Q(U rf ) -Q(U lf )) • n pf � ≤ 0.

Construction of approximate one-dimensional Riemann solvers by means of the Lagrange-to-Euler mapping

This section provides a general and systematic framework for constructing simple approximate Riemann solvers for one-dimensional hyperbolic systems of conservation laws written under Eulerian form. The underlying methodology stems from the transformation relating the Lagrangian and the Eulerian representations of conservation laws which is further applied to the Lagrangian Riemann solvers to deduce their Eulerian counterparts. This manner of proceeding ensures the direct transfer of the properties (conservation, positivity and entropy control) satisfied by the Lagrangian solver to its Eulerian counterpart. This transformation has been initially introduced in [START_REF] Gallice | Schémas de type Godunov entropiques et positifs préservant les discontinuités de contact[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] for designing conservative, positive and entropic simple approximate Riemann solvers, and, recently reused for the gas dynamics equations in [START_REF] Chan | Positivity preserving and entropy consistent approximate Riemann solvers dedicated to the high-order MOOD-based Finite Volume discretization of Lagrangian and Eulerian gas dynamics[END_REF]. It is worth noticing that this approach has been extended not only to the magnetohydrodynamics equations written under Powell's form but also to hyperbolic systems with source terms leading to well-balanced numerical discretizations, refer to [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Nonconservatifs[END_REF].

One-dimensional Lagrange-to-Euler mapping

This section studies the one-dimensional Riemann problem located at the subface interface which is required to construct the numerical flux approximation by means of the generic expression [START_REF] Einfeldt | On Godunov-type methods for gas dynamics[END_REF]. Let us consider the subface characterized by the unit outward n, the corresponding Riemann problem reads

(RP E )          ∂U ∂t + ∂F n (U) ∂x n = 0, where F n (U) = F(U)n, U(x n , 0) = � U l if x n < 0, U r if x n ≥ 0.
Here, x n = x • n denotes the space variable in the direction normal to the interface. This Riemann problem is also equipped with the the entropy inequality

∂Σ ∂t + ∂Q n ∂x n ≤ 0, ( 30
)
where

Q n = Q • n denotes the entropy flux.
From now on, we focus on systems of conservation laws describing physical phenomena in the domain of continuum mechanics, for instance gas dynamics, shallow water equations, Magneto-HydroDynamics (MHD), hyperelasticity... The interested reader might refer for instance to [START_REF] Kulikowskii | Mathematical Aspects of Numerical Solution of Hyperbolic Systems[END_REF] wherein several systems of conservation laws originated from continuum mechanics are described and studied. In this framework, we assume that the first component of the general system of conservation laws (1) is the mass conservation equation

∂ρ ∂t + ∇ • (ρv) = 0, (31)
where ρ > 0 is the mass density and v the material velocity. Thus, the first component of the Riemann problem (RP E ) reads

∂ρ ∂t + ∇ • (ρv n ) = 0, ( 32
)
where v n = v • n is the projection of the material velocity onto the unit normal n.

We construct the Lagrangian counterpart of (RP E ) introducing the Lagrange-to-Euler mapping

m � -→ x n (m, t) such that dx n = 1 ρ dm + v n dt,
is an exact differential and m denotes the Lagrangian mass coordinate. By construction

∂x n ∂m (m, t) = 1 ρ , and ∂x n ∂t (m, t) = v n ,
and since dx n is an exact differential the following compatibility condition holds true

∂τ ∂t - ∂v n ∂m = 0, (33)
where τ = 1 ρ is the specific volume. This is nothing but the Lagrangian mass/volume equation. Let us point out that we utilize the same notation for the Lagrangian and the Eulerian time. Moreover, the same physical quantity might be indifferently expressed either in terms of the Lagrangian coordinates (m, t) or in terms of the Eulerian ones (x n , t) knowing that x n = x n (m, t). This amounts to write formally

U(m, t) = U(x n (m, t), t).
Taking the time derivative of the foregoing identitity holding m fixed, i.e., the Lagrangian time derivative, and applying the chain rule leads to

∂U ∂t (m, t) � � � m = ∂U ∂t (m, t) � � � xn + v n ∂U ∂x n (x n , t) � � � t .
Therefore, the following identity holds true

ρ ∂ ∂t (τ U)(m, t) � � � m = ∂U ∂t (m, t) � � � xn + ∂(v n U) ∂x n (x n , t) � � � t .
Substituting the Eulerian time derivative thanks to (RP E ) in the foregoing equation leads to

∂(τ U) ∂t + ∂ ∂m (F n -v n U) = 0. ( 34
)
This system is the Lagrangian counterpart of (RP E ). However, one notices that its first component is trivial, we thus replace it by [START_REF] Noh | Errors for calculations of strong shocks using an artificial viscosity and an artificial heat flux[END_REF]. Finally, to complete the definition of the Lagrange-to-Euler mapping let us introduce the expressions of the Lagrangian variable and flux in terms of their Eulerian counterparts

V = τ (U -ρe 1 ) + τ e 1
, where e 1 = (1, 0, . . . , 0) t ,

G n = F n -v n U -v n e 1 .
The foregoing formulas define the Lagrange-to-Euler transformation which allows us to deduce the Lagrangian Riemann problem

(RP L )          ∂V ∂t + ∂G n (V) ∂m = 0, V(x n , 0) = � V l if m < 0, V r if m ≥ 0,
which is the Lagrangian counterpart of (RP E ). Here, V = V(m, t) and G n are respectively the vector of conservative variables and the flux vector written under Lagrangian representation. The foregoing Lagrangian Riemann problem is also equipped with the Lagrangian entropy inequality

∂σ ∂t + ∂q n ∂m ≤ 0, ( 35
)
where the Lagrangian entropy and entropy flux pair are written in terms of their Eulerian counterparts σ = τ Σ and

q n = Q n -v n Σ. Let W E (U l , U r , ξ E ),
where ξ E = xn t , denotes the Eulerian Riemann solver which consists of an approximate solution to (RP E ). Similarly, we define W L (V l , V r , ξ L ), where ξ L = m t , the Lagrangian Riemann solver which also consists of an approximate solution to (RP L ). We assume that the Eulerian approximate Riemann solver fulfills the basic requirements

• W E (U l , U r , ξ E ) = U l for -ξ E large enough, • W E (U l , U r , ξ E ) = U r for ξ E large enough, • W E (U, U, ξ E ) = U.
We prescribe similar assumptions for the Lagrangian Riemann solver.

Consistency of approximate Riemann solvers with the integral forms of the conservation law

and the entropy inequality Here, we recall the fundamental notion of consistency that has been initially introduced in the seminal works [START_REF] Harten | A Random Choice Finite Difference Scheme for Hyperbolic Conservation Laws[END_REF][START_REF] Harten | On upstream Differencing and Godunov-Type schemes for Hyperbolic Conservation Laws[END_REF]. Integrating the conservation law (RP) E over [-Δx l , 0] × [0, Δt] and replacing U(x n , t) by its approximation

W E (U l , U r , ξ E ) leads to � 0 -Δx l � W E (U l , U r , x n Δt ) -U l � dx n + � Δt 0 � F n (U(0 -, t)) -F n (U l ) � dt = 0.
This incites us to define the left-sided flux as follows

F - n = F n (U l ) - 1 Δt � 0 -Δx l � W E (U l , U r , x n Δt ) -U l � dx n .

Now, making the change of variable ξ E = xn

Δt in the foregoing integral and noticing that for -ξ E large enough W E (U l , U r , ξ E ) = U l , we rewrite the left-sided flux employing the compact formula

F - n = F n (U l ) - � 0 -∞ (W E (U l , U r , ξ E ) -U l ) dξ E . ( 36
)
Similarly, the right-sided flux writes

F + n = F n (U r ) + � +∞ 0 (W E (U l , U r , ξ E ) -U r ) dξ E . ( 37
)
Subtracting [START_REF] Pandolfi | Numerical Instabilities in Upwind Methods: Analysis and Cures for the Carbuncle Phenomenon[END_REF] to [START_REF] Peery | Blunt-body flow simulations[END_REF] yields

F + n -F - n = (38) � 0 -∞ (W E (U l , U r , ξ E ) -U l ) dξ E + � +∞ 0 (W E (U l , U r , ξ E ) -U r ) dξ E + F n (U r ) -F n (U l ).
The approximate Riemann solver, W E , is consistent with the integral form of the Riemann problem, (RP E ), if and only if the left and right-sided fluxes are equal that is

� 0 -∞ (W E (U l , U r , ξ E ) -U l ) dξ E + � +∞ 0 (W E (U l , U r , ξ E ) -U r ) dξ E (39) +F n (U r ) -F n (U l ) = 0.
Then, the approximate Riemann solver induces a Finite Volume Godunov-type scheme which is conservative by construction.

Similarly to what has been done for the flux, we also introduce the left-sided and the right-sided entropy fluxes

Q - n =Q n (U l ) - � 0 -∞ [Σ(W E (U l , U r , ξ E )) -Σ(U l )] dξ E , ( 40a
)
Q + n =Q n (U r ) + � +∞ 0 [Σ(W E (U l , U r , ξ E )) -Σ(U r )] dξ E . ( 40b
)
Then, the approximate Riemann solver, W E , is consistent with the integral form of the entropy inequality [START_REF] Maire | A nominally second-order accurate finite volume cell-centered scheme for anisotropic diffusion on two-dimensional unstructured grids[END_REF] if and only if

Q + n -Q - n ≤ 0.
Namely, substituting (40a) and (40b) into the foregoing difference yields the inequality

� 0 -∞ [Σ(W E (U l , U r , ξ E )) -Σ(U l )] dξ E + � +∞ 0 [Σ(W E (U l , U r , ξ E )) -Σ(U r )] dξ E (41) + Q n (U r ) -Q n (U l ) ≤ 0.
If this inequality holds true, the Riemann solver induces an entropic Finite Volume scheme. Obviously, similar definitions might be introduced for the Lagrangian Riemann solver W L (V l , V r , ξ L ) and we omit it for the sake of conciseness.

Lagrangian simple Riemann solvers and construction of their Eulerian counterparts via the

Lagrange-to-Euler mapping

Definition of the simple Lagrangian Riemann solver

In this section, we focus on a particular class of approximate Riemann solvers initially introduced in [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Nonconservatifs[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] and named simple Riemann solvers. The Lagrangian Riemann solver W L (V l , V r , ξ L ), where ξ L = m t , represents an approximate solution of the Lagrangian Riemann problem (RP L ). It is a simple Riemann solver if and only if it consists of m + 1 constant states V k , k = 1 . . . m + 1 separated by m discontinuities of slopes λ k , k = 1, . . . , m in the (m, t) plane. More precisely,

W L (V l , V r , m t ) =      V 1 = V l if m t < λ 1 , V k if λ k-1 ≤ m t < λ k , k = 2, . . . , m, V m+1 = V r if λ m ≤ m t .
Here, the λ k for k = 1, . . . , m are the Lagrangian wave speeds in ascending order in the (m, t) plane and thus homogeneous to m t . As its name suggests, the simple solver represents the simplest form of approximate Riemann solver. Among others, Roe [START_REF] Roe | Approximate Riemann solvers, parameter vectors and difference scheme[END_REF], HLL [START_REF] Harten | On upstream Differencing and Godunov-Type schemes for Hyperbolic Conservation Laws[END_REF] and HLLC [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] solvers are famous examples of simple approximate Riemann solvers. On the other hand, Godunov exact Riemann solver [START_REF] Godunov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF] and Osher solver [START_REF] Osher | Upwind schemes and boundary conditions with applications to Euler equations in general geometries[END_REF] are not.

Consistency of the simple Lagrangian solver with its underlying conservation law and entropy inequality

Similarly to what has been presented in section 3.2 the general expressions of the left and right-sided Lagrangian fluxes write

G - n =G n (V l ) - � 0 -∞ (W L (V l , V r , ξ L ) -V l ) dξ L , ( 42a
)
G + n =G n (V r ) + � +∞ 0 (W L (V l , V r , ξ L ) -V r ) dξ L . ( 42b
)
Replacing W L by its expression in the foregoing formulas leads to the explicit expressions of the left and right-sided fluxes in terms of the intermediate states and the wave speeds

G - n =G n (V l ) - m � k=1 λ (-) k (V k+1 -V k ), (43a)
G + n =G n (V r ) - m � k=1 λ (+) k (V k+1 -V k ), (43b)
where for any real, x, we denote by x (+) = 1 2 (|x| + x) and x (-) = 1 2 (|x|x) respectively its positive and negative part. We claim that the simple Lagrangian Riemann solver W L is consistent with the integral form of the conservation law (RP L ) if and only if

G - n = G + n . Indeed subtracting (43a) to (43b) this amounts to write - m � k=1 λ k (V k+1 -V k ) + G n (V r ) -G n (V l ) = 0. (44)
Therefore, the numerical flux at the interface writes

G n = 1 2 [G n (V l ) + G n (V r )] - 1 2 m � k=1 |λ k |(V k+1 -V k ).
As such, we have recovered the classical expression of the numerical flux decomposed into a centered part plus a viscous part, refer for instance to [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. It remains to investigate the consistency of the simple Lagrangian Riemann solver with the integral form of the Lagrangian entropy inequality [START_REF] Osher | Upwind schemes and boundary conditions with applications to Euler equations in general geometries[END_REF]. First, we compute the expression of the left and the right-sided Lagrangian entropy fluxes defined respectively by

q - n =q n (V l ) - � 0 -∞ [σ(W L (V l , V r , ξ L )) -σ(V l )] dξ L , q + n =q n (V r ) + � +∞ 0 [σ(W L (V l , V r , ξ L )) -σ(V r )] dξ L .
Replacing W L by its expression in terms of the intermediate states and the waves speeds yields

q - n =q n (V l ) - m � k=1 λ (-) k (σ k+1 -σ k ), q + n =q n (V r ) - m � k=1 λ (+) k (σ k+1 -σ k ),
where σ k = σ(V k ). We recall that the Riemann solver W L is consistent with the entropy inequality [START_REF] Osher | Upwind schemes and boundary conditions with applications to Euler equations in general geometries[END_REF] if and only if q + nq - n ≤ 0. Substituting the foregoing explicit expressions of the left and rightsided entropy fluxes in the aforementioned inequality, we thus claim that the simple Lagrangian Riemann solver W L is consistent with the entropy inequality if and only if

- m � k=1 λ k (σ k+1 -σ k ) + q n (V r ) -q n (V l ) ≤ 0. ( 45
)
Remark 3. It is worth pointing out that very often the Lagrangian systems of conservation laws governing physical phenomena in the domain of continuum mechanics, e.g., gas dynamics, MHD, hyperelasticity, are characterized by a zero entropy flux, i.e., q n = 0. This remarkable property which has been characterized in [START_REF] Després | Lagrangian systems of conservation laws[END_REF] considerably simplifies the study of Lagrangian simple approximate Riemann solvers from the point of view of entropy stability.

Construction of the simple Eulerian Riemann solver from its Lagrangian counterpart

Now, we construct the simple Eulerian Riemann solver W E from the simple Lagrangian Riemann one W L employing the Lagrange-to-Euler transformation introduced in section 3.1. This methodology has been initially introduced in [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Nonconservatifs[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] to deduce the Eulerian Riemann solver from its Lagrangian counterpart. Following [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Nonconservatifs[END_REF] we suppose that the Riemann solver W L satisfies the assumptions

• (H 1 ) λ k (τ k+1 -τ k ) + v n,k+1 -v n,k = 0, for k = 1, . . . , m. • (H 2 ) τ k ≥ 0, for k = 1, . . . , m.
Hypothesis (H 1 ) is nothing but the weak form of the volume/mass conservation equation [START_REF] Noh | Errors for calculations of strong shocks using an artificial viscosity and an artificial heat flux[END_REF] written across each discontinuity of speed λ k for k = 1, . . . , m. By virtue of (H 1 ) for any k = 1, . . . , m

v n,k + λ k τ k = v n,k+1 + λ k τ k+1 .
This in turn allows to define the Eulerian wave speeds

Λ k = v n,k + λ k τ k = v n,k+1 + λ k τ k+1 , for k = 1 . . . m. ( 46
)
Now, observing that

Λ k+1 -Λ k = v n,k+1 + λ k+1 τ k+1 -v n,k+1 -λ k τ k+1 = τ k+1 (λ k+1 -λ k )
and by virtue of (H 2 ), we deduce that the ordering of the Eulerian wave speeds is similar to that of the Lagrangian ones. Bearing this in mind, the Eulerian simple approximate Riemann solver

W E (U l , U r , ξ E ), where ξ E = xn t , is deduced from its Lagrangian counterpart W L (V l , V r , ξ L )
, where ξ L = m t , employing the Euler-to-Lagrange transformation which maps the Lagrangian vector of variables V onto its Eulerian counterpart U(V) = ρ(V-τ e 1 )+ρe 1 . Applying this transformation to the intermediate states of

W L yields W E (U l , U r , ξ E ) =      U 1 = U l = U(V l ) if ξ E < Λ 1 , U k = U(V k ) if Λ k-1 ≤ ξ E < Λ k for k = 2, . . . , m, U m+1 = U r = U(V r ) if Λ m ≤ ξ E .
Here, the Eulerian wave speeds are deduced from the Lagrangian ones by means of [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]. Thanks to (H 1 ), one can easily demonstrate that the weak form of the Eulerian mass conservation [START_REF] Menikoff | Errors when shock waves interact due to numerical shock width[END_REF] holds true, that is

-Λ k (ρ k+1 -ρ k ) + ρ k+1 v n,k+1 -ρ k v n,k = 0, for k = 1 • • • m. ( 47
)
This is nothing but the Eulerian version of (H 1 ).

Fundamental property relating the Eulerian and the Lagrangian fluxes

We finish this section by presenting the fundamental property of the simple Riemann solvers which allows to show that the difference between the left and the right-sided Eulerian fluxes is rigorously equal to the difference between the left and the right-sided Lagrangian ones provided that the underlying Eulerian and Lagrangian Riemann solvers are deduced one from the other by means of the Lagrange-to-Euler transformation. Substituting the expression of the simple Eulerian Riemann solver into [START_REF] Quirk | A contribution to the great Riemann solver debate[END_REF] then the difference of the fluxes becomes

F + n -F - n = - m � k=1 Λ k (U k+1 -U k ) + F n (U r ) -F n (U l ). ( 48
)
Now, substituting the expression of the Eulerian wave speeds [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] as functions of the Lagrangian ones and invoking the definition of the Eulerian intermediate states and fluxes allows to form their Lagrangian counterparts by means of the Lagrange-to-Euler transformation

F + n -F - n = - m � k=1 λ k (τ k+1 U k+1 -τ k U k ) + G n,r -G n,l + (v n,r -v n,l )e 1 = - m � k=1 λ k (V k+1 -V k ) + G n,r -G n,l + m � k=1 � λ k (τ k+1 -τ k ) + v n,k+1 -v n,k �������������������������������������������������������������������������������� =0 thanks to (H1) � e 1 = - m � k=1 λ k (V k+1 -V k ) + G n,r -G n,l .
Finally, we arrive at the fundamental formula

- m � k=1 Λ k (U k+1 -U k ) + F n (U r ) -F n (U l ) = - m � k=1 λ k (V k+1 -V k ) + G n (V r ) -G n (V l ). ( 49
)
This amounts to write that

F + n -F - n = G + n -G - n . ( 50
)
This result is a consequence of the construction of the Eulerian simple solver from the Lagrangian one utilizing the Lagrange-to-Euler mapping and (H 1 ) assumption. We point out that formula (49) has been already introduced in [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Nonconservatifs[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] to demonstrate the equivalence of the Lagrangian and Eulerian Riemann solvers regarding the consistency properties provided that assumption (H 1 ) holds true.

Remark 4. It is worth noting that the difference between the left and the right-sided Eulerian fluxes is always equal to the difference between the left and the right-sided Lagrangian fluxes, i.e. (49) always holds true, even if the approximate Riemann solver is not consistent with the integral form of the conservation law.

Following the same methodology, we compute the difference between the Eulerian left and rightsided entropy fluxes by subtracting (40a) to (40b) for the simple Eulerian Riemann solver

W E Q + n -Q - n = - m � k=1 Λ l (Σ k+1 -Σ k ) + Q n (U r ) -Q n (U l ).
Recalling that Eulerian wave speed satisfies

Λ k = λ k τ k + v n,k = λ k τ k+1 + v n,k+1
, and, the Eulerian entropy and entropy flux are expressed in terms of their Lagrangian counterparts by Σ = ρσ and

Q n = q n + ρσv n yields 20 - m � k=1 Λ k (Σ k+1 -Σ k ) + Q n (U r ) -Q n (U l ) = - m � k=1 λ k (σ k+1 -σ k ) + q n (V r ) -q n (V l ). ( 51
)
This amounts to write

Q + n -Q - n = q + n -q - n .

Revisiting the node-based conservation and entropy conditions in the case of simple Riemann solvers

The conservation and entropy properties of the generic Finite Volume scheme (6) developed in section 2.6 rely respectively on the node-based conservation condition [START_REF] Jacq | A Nominally Second-Order Cell-Centered Finite Volume Scheme for Simulating Three-Dimensional Anisotropic Diffusion Equations on Unstructured Grids[END_REF] and the node-based entropy condition [START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF]. We have shown that these node-based conditions are sufficient to ensure that the aforementioned Finite Volume scheme is conservative and satisfies an entropy inequality under the time step condition [START_REF] Gaitonde | High-Speed Viscous Flows Past Blunt Bodies and Compression Corner with Flux-Split Methods[END_REF]. Let us recall that for a generic node p these node-based conditions are written in terms of the left and right-sided subface fluxes, F

Expression of the node-based conditions for a simple Riemann solver

Here, we shall further develop these node-based conservation and entropy conditions in the particular cases for which W l pf is a simple approximate Riemann solver. Consequently, W l pf consists of m + 1 constant states U k for k = 1 . . . m + 1 separated by m waves of speeds

Λ k for k = 1 . . . m W l pf (U lf , U rf , n pf , ξ, v p ) =      U 1 = U lf if ξ < Λ 1 , U k if Λ k-1 < ξ ≤ Λ k , k = 2, . . . , m, U m+1 = U rf if Λ m ≤ ξ.
Here, the approximate Riemann solver is attached to the subface f whose unit normal n pf points towards the cell characterized by the right state U rf . The self-similar variable reads ξ =

x•n pf t .
In what follows, for the sake of simplicity, we shorten the notation of the Riemann solver into W l pf (U lf , U rf , ξ). Before proceeding any further, let us introduce the following identity which shall be useful and holds true for any real valued function f

� 0 -∞ � f (W l pf (U lf , U rf , ξ)) -f (U lf ) � dξ+ � ∞ 0 � f (W l pf (U lf , U rf , ξ)) -f (U rf ) � dξ = (52) - m � k=1 Λ k [f (U k+1 ) -f (U k )] .
Its straightforward application to the node-based conservation condition [START_REF] Jacq | A Nominally Second-Order Cell-Centered Finite Volume Scheme for Simulating Three-Dimensional Anisotropic Diffusion Equations on Unstructured Grids[END_REF] turns it into

� f ∈SF (p) l pf    - � m � k=1 Λ k (U k+1 -U k ) � l,r + [F(U rf ) -F(U lf )] n pf    = 0, (53)
Similarly, by virtue of identity (52) the node-based entropy condition (29) becomes

� f ∈SF (p) l pf    - � m � k=1 Λ k (Σ(U k+1 ) -Σ(U k )) � l,r + [Q(U rf ) -Q(U lf )] • n pf    ≤ 0. ( 54
)

Lagrangian equivalence of the node-based conditions

So far we have written the expressions of the node-based conservation and entropy conditions obtained for a simple Eulerian Riemann solver. These expressions might be further simplified employing the fundamental property attached to simple Eulerian and Lagrangian Riemann solvers put in relation through the Lagrange-to-Euler mapping and the (H 1 ) hypothesis, refer to section 3.3. To this end, let us introduce W l pf,L (V lf , V rf , ξ L ), the Lagrangian counterpart of the Eulerian simple solver W l pf (U lf , U rf , ξ). The Lagrangian simple Riemann solver writes

W l pf,L (V lf , V rf , ξ L ) =      V 1 = V lf if ξ L < λ 1 , V k if Λ k-1 ≤ ξ L < λ k , k = 2, . . . , m, V m+1 = V rf if λ m ≤ ξ L .
Here, the Lagrangian self-similar variable is defined by ξ L = m t where m is the mass Lagrangian coordinate defined from the Eulerian coordinate x n pf . The Lagrangian intermediate states V k are connected to the Eulerian intermediate states through the Lagrange-to-Euler mapping and the Lagrangian wave speeds λ k are related to their Eulerian counterpart thanks to (H 1 ) hypothesis. Now, with an obvious notation adaptation, the fundamental relations (49) for the flux and the entropy flux (51) turn into

- m � k=1 Λ k (U k+1 -U k ) + F n,r -F n,l = - m � k=1 λ k (V k+1 -V k ) + G n,r -G n,l , - m � k=1 Λ k (Σ k+1 -Σ k ) + Q n,r -Q n,l = - m � k=1 λ k (σ k+1 -σ k ) + q n,r -q n,l .
Finally, utilizing these identities, the node-based Eulerian conservation and entropy conditions for the Finite Volume scheme boil down to a Lagrangian expression. More precisely, the Eulerian node-based conservation condition (53) is equivalent to the Lagrangian one

� f ∈SF (p) l pf    - � m � k=1 λ k (V k+1 -V k ) � l,r + G n pf ,r -G n pf ,l    = 0. (55)
Similarly the Eulerian node-based entropy condition (54) is equivalent to the Lagrangian one

� f ∈SF (p) l pf    - � m � k=1 λ k (σ(V k+1 ) -σ(V k )) � l,r + [q(V r ) -q(V l )] • n pf    ≤ 0. ( 56
)
The interest of the above formulations is their simplicity as they are expressed only in terms of Lagrangian fluxes and unknowns. Particularly, for fluid dynamics systems having a zero entropy flux Lagrangian formulation, refer to remark 3, the entropy condition boils down to the node-based condition

� f ∈SF (p) l pf � m � k=1 λ k (σ(V k+1 ) -σ(V k )) � l,r ≥ 0. ( 57
)
These Lagrangian formulations of the node-based conservation and entropy conditions shall be investigated further in the next section presenting the application of this original framework to the gas dynamics system.

Application to the system of gas dynamics

In this section, we aim at describing the main building blocks of an original Finite Volume scheme for discretizing the gas dynamics system onto unstructured general grids.

Governing equations

The gas dynamics system of conservation laws expresses the conservation of mass, momentum and total energy and is written under Eulerian representation as follows

∂U ∂t + ∇ • F(U) = 0.
Here, U = U(x, t), where x ∈ R d , is the vector of conservative variables. This vector writes U = (ρ, ρv, ρe) t ∈ R d+2 where ρ is the mass density, v the velocity vector and e the specific total energy. The physical flux is described by the (d + 2) × d tensor

F(U) =   ρv t ρv ⊗ v + pI d ρev t + pv t   ,
where p denotes the thermodynamic pressure. The specific internal energy is given by ε = e -1 2 v 2 . Let η be the specific physical entropy, and τ = 1 ρ the specific volume, we make the fundamental assumption that (τ, η) � → ε(τ, η) is strictly convex which is equivalent to assume that (τ, ε) � → η(τ, ε) is strictly concave, refer to [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]. We work with the particular entropy, entropy flux pair (Σ, Q) = (-ρη, -ρvη) and thus the gas dynamics system of conservation laws is equipped with the entropy inequality ∂ρη ∂t

+ ∇ • (ρηv) ≥ 0. ( 58
)
The thermodynamic closure of this system of conservation laws is ensured by means of the complete equation of state

p(τ, η) = - ∂ε ∂τ , θ(τ, η) = ∂ε ∂η . ( 59
)
We make the classical assumption that the absolute temperature is strictly positive: θ > 0. By virtue of (59), writing the differential of ε(τ, η) leads to the fundamental Gibbs relation

θ dη = p dτ + dε. ( 60
)
The convexity of the specific internal energy with respect to the specific volume allows us to define the isentropic sound speed

a 2 τ 2 = - ∂p ∂τ = ∂ 2 ε ∂τ 2 . ( 61
)

Description of the one-dimensional Eulerian and Lagrangian systems of conservation laws

Let n be the unit normal to a generic interface along which we shall define Eulerian and Lagrangian one-dimensional problems in the normal direction and their associated Riemann solvers and finally construct the numerical approximation of the subfluxes of our Finite Volume scheme. In what follows, we shall restrict our developments to the bidimensional space, i.e., d = 2. In this framework, t is the unit vector such that (t, n) is a direct orthonormal basis attached to the generic interface. The normal and tangential components of the velocity write respectively v n = v • n and v t = v•t, and, obviously, v = v n n+v t t. The vector of conservative variables and the flux projected onto the normal direction n write

U =     ρ ρv n ρv t ρe     , F n = Fn =     ρv n ρv 2 n + p ρv n v t ρv n e + pv n     .
Therefore, the one-dimensional Eulerian system associated to the gas dynamics system in the n direction, where

x n = x • n, reads ∂U ∂t + ∂F n (U) ∂x n = 0.
This system is hyperbolic and admits the four following eigenvalues v na, v n with multiplicity 2 and v n + a, and is equipped with the entropy inequality

∂ρη ∂t + ∂ ∂x n (ρηv n ) ≥ 0. ( 62
)
Employing the Lagrange-to-Euler mapping introduced in section 3.1 we can derive the corresponding one-dimensional Lagrangian system

∂V ∂t + ∂G n (V) ∂m = 0, ( 63
)
where m is the Lagrangian mass coordinate related to the Eulerian coordinate x n . The Lagrangian vector of conservative variables and the Lagrangian flux write

V =     τ v n v t e     , G n =     -v n p 0 pv n     .
The Lagrangian system is also hyperbolic and admits the four eigenvaluesa τ , 0 with multiplicity 2 and a τ . The selection of physically admissible weak solutions is ensured supplementing this system of conservation laws by the entropy inequality ∂η ∂t ≥ 0.

(64)

The Lagrangian simple Riemann solver and its main properties 5.3.1. Construction of the Lagrangian Riemann solver

We consider the Lagrangian simple approximate Riemann solver corresponding to the foregoing one-dimensional Lagrangian gas dynamics system. This approximate Riemann solver structure mimics the continuous structure of the one-dimensional Lagrangian system of conservation laws. It is naturally composed of four states V l , V � l , V � r and V r separated respectively by discontinuities of speeds -λ l , 0 and λ r in the (m, t) plane. Here, λ l and λ r are positive real parameters which shall be constrained to ensure the positivity and entropy stability properties of the Riemann solver following the methodology introduced initially in [START_REF] Gallice | Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Nonconservatifs[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] and revised recently in [START_REF] Chan | Positivity preserving and entropy consistent approximate Riemann solvers dedicated to the high-order MOOD-based Finite Volume discretization of Lagrangian and Eulerian gas dynamics[END_REF]. Bearing this in mind the Lagrangian Riemann solver writes

W L � V l , V r , m t � =          V l if m t ≤ -λ l , V � l if -λ l < m t ≤ 0, V � r if 0 < m t ≤ λ r , V r if λ r < m t .
The states components write V s = (τ s , v n,s , v t,s , e s ) t and the intermediate states components

V � s = (τ � s , v � n,s , v � t,s , e � s ) t for s = l, r.
Assuming (H 1 ) hypothesis is satisfied yields

λ l (τ � l -τ l ) -(v � n,l -v n,l ) = 0, 0(τ � r -τ � l ) -(v � n,r -v � n,l ) = 0, -λ r (τ r -τ � r ) -(v n,r -v � n,r ) = 0.
This implies v � n,r = v � n,l and we denote v � n the common value of the velocity, i.e., v

� n = v � n,r = v � n,l
. This shows that (H 1 ) ensures the continuity of kinematic velocity through the contact waves which is a rather satisfying physical behavior. Finally, the foregoing system boils down to

λ l (τ � l -τ l ) -(v � n -v n,l ) = 0, (65a) λ r (τ � r -τ r ) + v � n -v n,r = 0. ( 65b
λ l (V � l -V l ) + G n,l -G n,l = 0, (66a)
-λ r (V r -V � r ) + G n,r -G n,r = 0, (66b)
where G n,s = G n (V s ) for s = l, r. Combining the first components of (66a), (66b) with (65a), (65b) leads to v n,l = v � n = v n,r . On the other hand, we make the structural assumption that the intermediate total energy flux writes under the form

(pv n ) s = p s v � n , for s = l, r. ( 67
)
Gathering these results, we arrive at the following expressions of the intermediate states and fluxes

V � s =     τ � s v � n v � t,s e � s     , and G n,s =     -v � n p s 0 p s v � n     , for s = l, r.
(H 1 ) hypothesis and the structural assumption (67) for the total energy flux allow us to reduce the number of scalar unknowns to 9, which must satisfy 8 scalar equations corresponding to the vectorial equations (66a) and (66b), that is

(S l )          λ l (τ � l -τ l ) -(v � n -v n,l ) = 0, λ l (v � n -v n,l ) + p l -p l = 0, λ l (v � t,l -v t,l ) = 0, λ l (e � l -e l ) + p l v � n -p l v n,l = 0, (S r )          λ r (τ � r -τ r ) + v � n -v n,r = 0, λ r (v � n -v n,r ) -(p r -p r ) = 0, λ r (v � t,r -v t,r ) = 0, λ r (e � r -e r ) -(p r v � n -p r v n,r ) = 0.
We observe that the tangential velocity is conserved through the left and right discontinuities, i.e., v � t,s = v t,s for s = l, r since λ s > 0. It is worth noticing that, the intermediate normal velocity v � n might be viewed as a parameter in terms of which the 8 remaining unknowns might be expressed.

Positivity preserving and entropy control of the Lagrangian Riemann solver

We briefly describe how to ensure that the foregoing simple Lagrangian Riemann solver ensures not only the positivity of the intermediate specific volume, i.e., τ � s ≥ 0, and internal energy, i.e., ε � s ≥ 0, but also the intermediate entropy control, i.e., η � sη s ≥ 0, for s = l, r. For a more detailed presentation of this topic the interested reader might refer to [START_REF] Chan | Positivity preserving and entropy consistent approximate Riemann solvers dedicated to the high-order MOOD-based Finite Volume discretization of Lagrangian and Eulerian gas dynamics[END_REF]. First, eliminating the normal velocity increment v � nv n,s between the mass/volume equation and the momentum equation of (S s ) yields for s = l, r

p s -p s = -λ 2 s (τ � s -τ s ). (68)
Further, dot-multiplying the momentum equation of (S s ) by 1 2 (v � n + v n,s ), and subtracting it to the total energy equation provides us the internal energy equation for s = l, r

ε � s -ε s + p s + p s 2 (τ � s -τ s ) = 0, ( 69
) since ε � s = e � s -1 2 (v � n,s ) 2 -1 2 (v � t,s ) 2
for s = l, r and v � t,s = v t,s . Finally, substituting (68) into (69) leads to the expression of the post-discontinuity internal energy for s = l, r

ε � s = ε s -p s (τ � s -τ s ) + λ s 2 (τ � s -τ s ) 2 . ( 70
)
This equation is fundamental since it allows a straightforward derivation of a positivity condition for those specific internal energies. It also facilitates the study of the entropy production related to the simple approximate Riemann solver. Equation (70) shows that the specific internal energy ε � s is a convex quadratic function with respect to τ � sτ s , and thus it is always greater than its minimum value

ε � s ≥ ε s - p 2 s 2λ 2 s .
Therefore, the specific internal energy ε � s is positive provided that the wave speed λ s satisfies the condition

λ s ≥ p s √ 2ε s .
This condition has been already proposed 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, noticing that for a convex equation of state, i.e., τ � → p(τ, η) strictly convex, there holds a 2 τ 2 ≥ p 2 2ε [START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF], the foregoing positivity condition turns into

λ s ≥ a s τ s , for s = l, r. ( 71
)
The positivity conditions of the specific volumes are readily obtained written under the form

λ l ≥ - v � n -v n,l τ l , and λ r ≥ v � n -v n,r τ r .
It appears that these conditions are parametrized by the normal velocity v � n . We are then able to gather the positivity conditions of specific internal energy and specific volume into the global conditions

λ l ≥ max � a l τ l , - v � n -v n,l τ l � , and λ r ≥ max � a r τ r , v � n -v n,r τ r � . ( 72
)
Finally, we address the entropy control of the Lagrangian approximate Riemann solver. Decomposing the specific internal energy variation ε � sε s across the discontinuities into an isentropic process followed by an isochoric one, we are able to express the entropy production in terms of the specific volume variation τ � sτ s . The study of this entropy production term [START_REF] Chan | Positivity preserving and entropy consistent approximate Riemann solvers dedicated to the high-order MOOD-based Finite Volume discretization of Lagrangian and Eulerian gas dynamics[END_REF] shows that specific entropy increases across the discontinuity, i.e., η � sη s ≥ 0 provided that the wave speed satisfies the condition

λ 2 s ≥ a 2 (τ s , η s ) τ 2 s for all τ s ∈ (τ s , τ � s ), for s = l, r. ( 73
)
We note that this condition has been also derived in the framework of relaxation scheme utilizing a relatively cumbersome approach, refer to [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws[END_REF].

The corresponding simple Eulerian Riemann solver

Relying on hypothesis (H 1 ) we can deduce the Eulerian wave speeds Λ l , Λ 0 and Λ r from their Lagrangian counterparts setting

Λ l = v n,l -λ l τ l = v � n -λ l τ � l , Λ 0 = v � n , Λ r = v � n + λ r τ � r = v n,r + λ r τ r . ( 74
)
It is clear that the Eulerian wave speeds are ordered, i.e. Λ l ≤ Λ 0 ≤ Λ r , provided that the Lagrangian approximate Riemann solver is positivity preserving, i.e., τ � s ≥ 0. This holds true granted that the Lagrangian wave speeds satisfy (72). Bearing this in mind, we are able to deduce the Eulerian approximate Riemann solver from its Lagrangian counterpart as follows

W E � U l , U r , x n t � =          U l if xn t ≤ Λ l , U � l = U(V � l ) if Λ l < xn t ≤ Λ 0 , U � r = U(V � r ) if Λ 0 < xn t ≤ Λ r , U r if Λ r < xn t .
Here, V � → U(V) is the Lagrange-to-Euler mapping introduced in section 3 which allows us to define straightforwardly the intermediate states of the simple Eulerian solver from its Lagrangian counterpart. Consequently, the Eulerian intermediate states read

U � s = (ρ � s , ρ � s v � n , ρ � s v � t,s , ρ � s e � s ) t knowing that ρ � s = (τ � s ) -1 for s = l, r.
We observe that the Eulerian approximate Riemann solver is also parametrized by the normal star-velocity. More importantly, by construction, the Eulerian approximate Riemann solver has the same properties than its Lagrangian counterpart, i.e., it preserves the positivity of mass density, specific energy and ensures entropy increase under specific conditions on the Lagrangian wave speeds, refer to (73).

Consistency of the Lagrangian Riemann solver with its underlying conservation law

Let us investigate the consistency of the simple Lagrangian Riemann solver with the onedimensional conservation law (63), refer to section 3.3.2. To this end, we sum relations (66a) and (66b) characterizing the intermediate fluxes G n,l , G n,r and we get

G n,r -G n,l = λ l (V � l -V l ) -λ r (V r -V � r ) + G n,r -G n,l . (75)
On the other hand, utilizing the expression of the components of G n,r and G n,l we arrive at

G n,r -G n,l = (p r -p l )     0 1 0 v � n     . ( 76
)
Now, we address the consistency of the simple Lagrangian Riemann solver with its underlying conservation law by studying the solutions of G n,r -G n,l = 0. To this end, we compute the difference p rp l summing the second equations of (S l ) and (S r )

p r -p l = (λ l + λ r ) � v � n - � λ l v n,l + λ r v n,r λ l + λ r - (p r -p l ) λ r + λ l �� . ( 77
)
This equation incites us to introduce the normal velocity

v n = λ l v n,l + λ r v n,r λ l + λ r - (p r -p l ) λ r + λ l , ( 78
)
which is nothing but the normal velocity of the classical Godunov acoustic solver, refer for instance to [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. Bearing this in mind, we claim that

• If v � n = v n ,
then the simple Lagrangian Riemann solver is consistent with its underlying conservation law. Thus, the simple Lagrangian Riemann solver induces a classical conservative Godunov-type Finite Volume scheme.

• If v � n � v n ,
then the simple Lagrangian Riemann solver is not consistent with its underlying conservation law and thus not does not induce a conservative Godunov-type Finite Volume scheme.

In what follows, we investigate further the second case for which in general v � n � v n and we shall demonstrate how to retrieve a global conservation property for the Finite Volume scheme by means of the node-based conservation conditions introduced in section 4.

The nodal solver 5.4.1. Expression of the node-based conservation condition

We have seen in section 4, that the global conservation of the Finite Volume scheme is ensured provided that the Lagrangian node-based condition (55) is satisfied. This condition reads

� f ∈SF (p) l pf    - � m � k=1 λ k (V k+1 -V k ) � l,r + G n pf ,r -G n pf ,l    = 0.
We are going to develop it in the case of gas dynamics using the properties of the simple Lagrangian Riemann solver constructed previously. First, we express the generic term between curly brackets at the left-hand side of the foregoing conservation condition with an obvious notation adaptation

- � m � k=1 λ k (V k+1 -V k ) � l,r + G n pf ,r -G n pf ,l = λ l (V � l -V l ) -λ r (V r -V � r ) + G n,r -G n,l , (79)
The right-hand side of the foregoing equation coincides precisely with the right-hand side of (75) which has been derived for studying the consistency of the simple Lagrangian solver with its underlying conservation laws. Thus, by virtue of ( 75) and (76), equation ( 79) turns successively into

- � m � k=1 λ k (V k+1 -V k ) � l,r + G n pf ,r -G n pf ,l =G n,r -G n,l =(p r -p l )     0 1 0 v � n     = (p rf -p lf )     0 1 0 v � n pf     .
We observe that the last line has been written using once more an obvious notation adaptation. Namely, p rf and p lf denote the left and the right-sided interfacial pressures attached to subface f . Finally, we arrive at the conclusion that our Finite Volume scheme is conservative provided that the following node-based condition is satisfied

� f ∈SF (p) l pf (p rf -p lf )     0 1 0 v � n pf     = 0. (80)
Realizing that the second and third components of the foregoing vector correspond to n pf , then condition (80) implies

� f ∈SF (p) l pf (p rf -p lf )n pf = 0. (81)
We observe that v � n pf is still an unknown parameter attached to each subface impinging at node p. Thus, the number of unknown parameters at node p is much greater than the number of equations given by the conservation condition (80). Therefore, to close this system of equations, we assume that the parameter v � n pf is the projection of the unknown nodal vector v p onto the unit normal n pf , that is

v � n pf = v p • n pf . ( 82
)
This fundamental assumption drastically reduces the number of unknowns to the vectorial unknown v p , which can be interpreted as an approximation of the nodal velocity. With this assumption the conservation condition (81) is equivalent to the conservation condition (80). Thanks to (77) and (78) we are able to express the interface pressures difference into the node-based conservation which becomes

� f ∈SF (p) l pf (λ lf + λ rf )(v p • n pf -v n pf )n pf = 0,
where v n,pf is obtained from (78) with obvious notation adaptation

v n pf = λ lf v n pf ,l + λ rf v n pf ,r λ lf + λ rf - p rf -p lf λ lf + λ rf .
Finally, the node-based conservation condition (81) boils down to the system

� f ∈SF (p) l pf (λ lf + λ rf )(n pf ⊗ n pf )v p = � f ∈SF (p) l pf (λ lf + λ lf )v n pcf n pf . ( 83
)
This system always admits a unique solution which provides an approximation of the nodal velocity v p . It is thus called a nodal solver. We point out that the foregoing system has been already obtained when constructing a cell-centered Finite Volume discretization of multidimensional Lagrangian hydrodynamics [START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF]. It has been also retrieved in [START_REF] Shen | A robust and contact resolving Riemann solver on unstructured mesh, Part I, Euler method[END_REF] for designing a Finite Volume scheme for Eulerian gas dynamics.

Expression of the node-based entropy condition

Gathering the foregoing results, we are now in position to express what has become the nodalbased entropy condition (57). With the present simple Lagrangian Riemann solver the left-hand side of (57) writes

� f ∈SF (p) l pf � m � k=1 λ k (σ(V k+1 ) -σ(V k )) � l,r = � f ∈SF (p) l pf [λ l (η � l -η l ) + λ r (η � r -η r )] l,r .
Here, we have use σ = -η, where η is the physical concave entropy. Consequently the node-based entropy condition for our Finite Volume scheme turns into

� f ∈SF (p) l pf [λ l (η � l -η l ) + λ r (η � r -η r )] l,r ≥ 0, (84)
which is satisfied provided that the simple Lagrangian Riemann solver is entropy stable, refer to condition (73). This inequality might be satisfied relatively simply by utilizing a classically tuning of the wave speeds of the Riemann solver similarly to what has been undertaken in the one-dimensional framework, refer to [START_REF] Chan | Positivity preserving and entropy consistent approximate Riemann solvers dedicated to the high-order MOOD-based Finite Volume discretization of Lagrangian and Eulerian gas dynamics[END_REF].

Summary of the Eulerian multi-dimensional Finite Volume scheme

We recall that the multi-dimensional Finite Volume scheme writes

U n+1 c -U n c + Δt |ω c | � p∈P(c) � f ∈SF (pc) l pcf F pcf = 0,
where F pcf is the left-sided flux with respect to the subface f and the unit outward normal n pcf . The left-sided flux (36) for a simple Eulerian Riemann solver is obtained by

F pcf ≡ F - n pcf = F n pcf (U c ) - � m � k=1 Λ (-) k (U k+1 -U k ) � c,d . ( 85
)
Here, the subscripts c and d correspond to the cells located respectively on the left and the right sides of the subface f with respect to the unit normal n pcf , refer to figure 2. On the other hand, the right-sided flux writes

F + n pcf = F n pcf (U d ) - � m � k=1 Λ (+) k (U k+1 -U k ) � c,d . ( 86
)
Taking the arithmetic average of the left and the right-sided fluxes allows us to define the averaged flux on the subface f

F n pcf = 1 2 � F n pcf (U c ) + F n pcf (U d ) � - 1 2 � m � k=1 |Λ k |(U k+1 -U k ) � c,d . ( 87
)
This is the numerical flux corresponding to the classical face-based Finite Volume method. Since, in this case, the numerical flux depends uniquely on the two states adjacent to the subface, we call it two-point flux.

F + n pcf -F - n pcf =G + n pcf -G - n pcf , thanks to (50) =(p pcf,r -p pcf,l )     0 1 0 v p • n pcf     , thanks to (76),
where p pcf,r and p pcf,l are the left and right-sided interface pressures attached to the subface f . Finally, combining the foregoing equation with (87) we arrive at

F - n pcf = 1 2 � F n pcf (U c ) + F n pcf (U d ) � - 1 2
� m � k=1 |Λ k |(U k+1 -U k ) � c,d (88)
- 1 2 (λ pcf,l + λ pcf,r ) � v p • n pcf -v n pcf �     0 1 0 v p • n pcf     .
This explicit expression of the numerical flux shows that in the particular case where v p • n pcf = v n pcf we retrieve the classical conservative face-based Finite Volume scheme. In the general case, v p • n pcf � v n pcf and the numerical flux depends not only on the two states adjacent to the subface but also on the states surrounding node p through the expression of the nodal velocity, refer to (83), hence the name multi-point flux. We have displayed in figure 4 the stencils of the Finite Volume scheme over the fragment of a polygonal grid for the two-point flux and the multi-point flux approximations. We observe that the stencil corresponding to multi-point flux approximation consists of the neighboring cells that share a node with the target cell.

With an obvious notation adaptation, the Eulerian wave speeds (74) related to the subface characterized by the unit normal n pcf read

Λ pcf,l =v n c • n pcf - λ pcf,l ρ n c = v p • n pcf - λ pcf,l ρ � pcf,l , Λ pcf,0 =v p • n pcf , Λ pcf,r =v p • n pcf + λ pcf,r ρ � pcf,r = v n d • n pcf + λ pcf,r ρ n d .
Here, -λ pcf,l ≤ 0 and λ pcf,r ≥ 0 are respectively the left and right-sided wave speeds of the underlying Lagrangian Riemann solver. The parameters λ pcf,l and λ pcf,r have been monitored such that the intermediate densities ρ � pcf,l and ρ � pcf,r are positive, refer to (72), and thus the Eulerian wave speeds are ordered as follows

Λ pcf,l ≤ Λ pcf,0 ≤ Λ pcf,r .
Gathering the foregoing results, the time step condition [START_REF] Gaitonde | High-Speed Viscous Flows Past Blunt Bodies and Compression Corner with Flux-Split Methods[END_REF] to ensure that the Finite Volume scheme is D-preserving becomes

Δt ≤ |ω c | � p∈P(c) � f ∈SF (pc) l pcf Λ (-) pcf,l , where Λ (-) pcf,l = 1 2 �� � �v n c • n pcf - λ pcf,l ρ n c � � � -v n c • n pcf + λ pcf,l ρ n c
� . This time step condition ensures that

U n+1 c
is a convex combination of U n c and all the intermediate states of the approximate Riemann solvers attached to each subface of cell ω c . This convex combination property implies the preservation of positivity for the mass density and the specific internal energy provided that the approximate Riemann solver is positivity preserving. This is indeed the case when the Lagrangian wave speeds λ pcf,l and λ pcf,r satisfies (72). Finally, observing that Λ

(-) pcf,l ≤ |v n c • n pcf | + λ pcf,l ρ n c leads to the practical time step condition Δt ≤ |ω c | � p∈P(c) � f ∈SF (pc) l pcf � |v n c • n pcf | + λ pcf,l ρ n c � . ( 89
)

Numerical results

From now on, we name 'two-point' scheme the classical conservative face-based Finite Volume scheme, while the proposed method based on node-based conservation is referred to as the 'multipoint' scheme. In this section we run an extensive numerical test case campaign to demonstrate the performance of both the first-order explicit two-point and multi-point schemes. These test cases are ran on various types of grids that can be made of triangular, quadrangular and polygonal cells. The CFL number is set to 0.5 for all cases. For all the test cases the triangular meshes are obtained employing the open source mesh generator Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and postprocessing facilities[END_REF], the structured quadrangular meshes are directly computed by the Finite Volume code and the polygonal grids result from of an in-house mesher based on Voronoi tessellation [START_REF] Loubère | ReALE: A reconnection-based arbitrary-Lagrangian-Eulerian method[END_REF].

Odd-even decoupling

The Odd-Even decoupling problem is a difficult simulation for certain numerical schemes. Quirk in [START_REF] Quirk | A contribution to the great Riemann solver debate[END_REF] proposed this so-called odd-even decoupling test, where one simulates the propagation of a planar shock but on a perturbed Cartesian mesh for which only the center line is slightly shifted. The computational domain is defined by Ω = [0, 800] × [0, 20] with a Cartesian uniform structured mesh with Δx = Δy = 1, with I = 800 cells in x-direction, J = 20 in y-direction. A cell is then labelled with i, j indices as usual for Cartesian mesh. The center-line, at y = 10, coincides with the j 0 labeled line of mesh with j 0 = 10. The perturbation of amplitude 10 -6 is alternating as follows for line j 0 :

x i,j0 ←x i , y i,j0 ←y j0 + (-1) i 10 -6 . The shock wave is traveling from left to right with a Mach number of Ma = 6, the initial domain is filled with a diatomic gas with (ρ 0 , u 0 , v 0 , p 0 , γ) = (1, 0, 0, 1, 1.4). The left side inflow values are determined using the Rankine-Hugoniot relations:

u s = Ma √ γ, ρ ∞ = (γ + 1)Ma 2 (γ -1)Ma 2 + 2 , u ∞ = u s 2(Ma 2 -1) (γ + 1)Ma 2 , p ∞ = 2γMa 2 -(γ -1) (γ + 1) ,
The left Boundary Condition (BC) corresponds to inflow boundary, top and bottom are wall-type BCs while the right one is left open. The final time is t final = 50. Obviously the exact solution is a vertical Ma = 6 shock wave traveling at speed u s . The instability reveals itself through symmetry perturbation of the shock wave as it travels along the duct. A sensitive measure of this perturbation, see [START_REF] Rodionov | Artificial viscosity Godunov-type schemes to cure the carbuncle phenomenon[END_REF], consists in computing the deviation ε 0 of the numerical solution:

ε 0 = max i,j (|ρ i,j -ρ i |), with ρ i = 1 J J � j=1 ρ i,j . (90)
The density color maps at the stopping time are displayed in figure 5a respectively for the twopoint (top) and the multi-point (bottom) schemes. The density contours have been plotted using 20 iso-lines in the range [1.4, 7.3]. The development of the odd-even instability is clear for the twopoint scheme whereas it seems to be absent for the multi-point scheme. In figure 5b one displays the plot of ε 0 as a function of the distance traveled by the shock wave, X s (t) = u s t, for the twopoint scheme (red curve) and the multi-point one (blue curve). The presence of the instability can be quantified with the plot of ε 0 , refer to figure 5b. We observe an exponential increase up to ε 0 > 1 = Δx for the two-point scheme, while the multi-point scheme value remains consistently at the order 10 -4 � Δx.

A more demanding version of the odd-even decoupling problem has been proposed by Rodionov in [START_REF] Rodionov | Artificial viscosity Godunov-type schemes to cure the carbuncle phenomenon[END_REF] by changing the seed of the instability. Here, we simulate the propagation of the same shock wave. However a small perturbation of the form The instability is much more pronounced with higher frequency for the two-point scheme. Contrarily the multi-point scheme does not present such an amplification of the initial perturbations, see figure 6a. This is also observed on the ε 0 curves, refer to figure 6b.

Hypersonic flow over half cylinder

The hypersonic flow over a half cylinder test case is a well-documented test case to challenge numerical methods. In particular for this flow some schemes may develop the infamous carbuncle phenomena when classical face-based Finite Volume upwind schemes are employed. Instead of having a smooth bow shock profile upstream of the half cylinder, the carbuncle failing features a pair of oblique shock ahead of the stagnation region, compromising the overall flow predictions around the cylinder, refer for instance to [START_REF] Peery | Blunt-body flow simulations[END_REF][START_REF] Gaitonde | High-Speed Viscous Flows Past Blunt Bodies and Compression Corner with Flux-Split Methods[END_REF][START_REF] Pandolfi | Numerical Instabilities in Upwind Methods: Analysis and Cures for the Carbuncle Phenomenon[END_REF] just to cite few references. Here, following [START_REF] Rodionov | Artificial viscosity Godunov-type schemes to cure the carbuncle phenomenon[END_REF], we simulate an inviscid flow at Mach Ma = 20 around a half cylinder blunt body subject to a incoming hypersonic flow characterized by (ρ 0 , u 0 , v 0 , p 0 , γ) = (1, Ma √ γ, 0, 1, 1.4). The steadystate resulting flow is simulated by means of an explicit time marching procedure, that is, the simulation ends when the residual is 6 orders of magnitude smaller than its initial value. The computational domain covers a large enough domain which contains half of a cylinder centered at the origin with a radius R = 1, and a left incoming hypersonic flow. At the cylinder surface a walltype boundary condition is considered, while bottom/upper boundary conditions are free outflow and inflow condition at the left boundary. Three types grid are tested to assess the robustness of our unstructured Finite Volume multi-point scheme: a triangular grid is composed of 5671 unstructured triangles, a quadrangular one made of 5000 structured quadrangles, and, a polygonal grid with 5632 unstructured polygons. Firstly, this test case is ran with the two-point scheme over the triangular grid, and, as can be observed in figure 7a the carbuncle instability clearly develops. On the contrary, such instability does not develop with the multi-point scheme over the triangular, quadrangular and polygonal grids, refer respectively to figures 7b, 8a and 8b. A quantitative comparison of these numerical results can be achieved plotting the pressure coefficient, C p = p-p0 1 2 ρ0u 2 0 , at the wall with respect to the angular position. The numerical C p is then compared to the approximate analytical value coming from the Newtonian theory, see [START_REF] Anderson | Hypersonic and high-temperature gas dynamics[END_REF]. The pressure coefficient, C p , computed by the multi-point scheme matches the Newtonian theory for angles below 40 • , while the two-point scheme results are polluted by the carbuncle effect, refer to figure 9a. A rather good convergence is observed for the multi-point scheme in figure 9b for successively refined structured grids made of 25 × 50, 50 × 100 and 100 × 200 quadrangular cells.

We pursue our numerical study comparing both two-point and multi-point schemes over the same polygonal grid for the half cylinder test case. It is worth mentioning that the polygonal grid results from a Voronoi tessellation. For this particular grid, a generic polygonal cell has exactly the same number of face-based or vertex-based neighbors. This implies, that the two-point Finite Volume scheme and the multi-point Finite Volume scheme share the very same stencil. The density map obtained with both schemes (top/bottom panels: two/multi-point schemes) over the polygonal grid under consideration is displayed in figure 10. 20 isolines ranging from ρ = 1 to ρ = 6.2 are plotted. Although the two-point scheme does not exhibit large spurious spike as the one observed in figure 7a, the isolines exhibit emerging instabilities behind the shock, see also the 3D view in figures 10b. This test case demonstrates that the multi-point solver does not generate any carbuncle-like instability on any type of grid. As such one has to be careful when praising the low-dissipation property of a solver such as the two-point one.

Sedov blast wave problem

The Sedov test case is a point explosion problem for which an exact solution can be derived, refer to [START_REF] Sedov | Similarity and Dimensional Methods in Mechanics[END_REF]. In the case of two-dimensional planar geometry, the problem consists of a cylindrical explosion generating a diverging shock wave. The computational domain Ω

= [-1.2, 1.2] × [-1.2, 1.2]
is initially filled with a perfect gas at rest characterized by the initial conditions (ρ 0 , u 0 , v 0 , p 0 , γ) = (1, 0, 0, 10 -6 , 7 5 ). The point explosion is initiated by an energy deposition in the vicinity of the origin, i.e., the pressure in the cells in contact with the origin is set to the value 0.397056. Reflective boundary conditions are applied on all boundaries. With this particular setup taken from [START_REF] Kamm | On efficient generation of numerically robust Sedov solutions[END_REF], the cylindrical shock wave radius is r shock = � x 2 + y 2 = 1 at the final time t final = 1 with a peak density ρ max = γ+1 γ-1 = 6. The exact density as a function of radius is plotted using a continuous black line in figure 12.

Similarly to what has been done for the previous test cases, we compared the numerical solutions obtained by the two-point and multi-point schemes on 400 × 400 uniform quadrangles in figure 11a and in figure 11b respectively. Then, the scattered plot of density with respect to the radius of each cell center is displayed in figures 12a and 12b for both schemes versus the analytical solution. These results exhibit the occurrence of a spurious effect along the x and y axes for the two-point scheme. Notice that the full 2π Sedov problem is actually run so that the boundary conditions cannot be blamed. The scattered plots also confirm the loss of cylindrical symmetry for the two-point scheme. Contrarily, the multi-point scheme does not present such parasitical effect, and, as such preserves the cylindrical symmetry, while being slightly more dissipative. These observations have already been done in [START_REF] Rodionov | Artificial viscosity Godunov-type schemes to cure the carbuncle phenomenon[END_REF] for instance. In order to confirm such a good behavior, the multi-point scheme is used on the following successively refined quadrangular grids: 100 × 100, 200 × 200 and quarter of a circle of radius of r = 1. A perfect gas with γ = 5/3 is initially assigned with a uniform density ρ 0 = 1 and a unit radial inward pointing velocity, hence the velocity components are initialized with v 0 (x, y) = (-x/r, -y/r) where r = � x 2 + y 2 denotes the radius and �v 0 � = 1. The initial pressure is p 0 = 10 -6 everywhere. A diverging cylindrical shock wave is generated at the origin. The state behind the shock wave is constant, whereas the state in front is not anymore due to the cylindrical converging flow. Thus, the exact solution at time t > 0 writes (ρ, v, p) ex (t) =

     � ρ 0 � γ+1 γ-1 � 2 , 0, 1 2 ρ 0 (γ + 1) 2 γ -1 � , if r < r s , � ρ 0 (1 -(t/r)) 2 , v 0 , 10 -6 � , if r ≥ r s . ( 91
)
Here, the radius of the shock wave is r s (t) = v s t with v s = 1 2 (γ -1)�v 0 �. With this setup and the final time t final = 0.6, we arrive at r s = 0.2 and a post-shock state characterized by the values ρ ex = 16, v ex = 0 and p ex = 16 3 . Symmetry BCs are prescribed on the axis x = 0 and y = 0, whereas at the outer radius, the exact space/time dependent velocity is imposed employing the analytical solution (91). A radial/polar grid made of N × N quadrangles is constructed such that Δr = 1/N and Δθ = π/(4N ) are constants. Let us notice that the cells in contact with the origin are triangles which are considered as degenerated quadrangles. We present the numerical density (colors and isolines) for a N = 200 × 200 polar grid at the final time t final = 0.6 for the two-point and the multi-point schemes in figures 15a and 15b respectively. The two-point scheme produces strong post-shock instabilities. We have verified that those instabilities are not generated by inappropriate boundary conditions. Contrarily the multi-point scheme presents a smooth symmetrical solution. Consistently with the results of Sedov problem, we present the scattered plot of the density with respect to the radius of the cell center for both schemes in figure 16a. Again we observe that the multi-point scheme handles such cylindrical flow without any loss of symmetry. Finally, a grid convergence analysis for the multi-point scheme is performed using the sequence of grids characterized by 50 2 , 100 2 and 200 2 polar cells. The results for the density are plotted in figure 16b. We point out that the multi-point scheme indeed converges towards the exact solution. Let us notice that the undershoot close to the origin corresponds to the so called wall-heating effect, refer to [START_REF] Menikoff | Errors when shock waves interact due to numerical shock width[END_REF][START_REF] Rider | Revisiting wall heating[END_REF].

Forward-facing step

Next we run the forward facing step problem proposed by Woodward and Colella in [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF]. This test case is a Mach 3 wind tunnel with a step. The computational domain is given by Ω = [0, 3] × [0, 1] [0.6, 3] × [0, 0.2] with the following initial conditions ρ = γ, v = (3, 0) t , p = 1, γ = 7/5. The final time is t final = 4 and reflective boundary conditions are applied on the upper and lower boundaries of the domain, whereas inflow and outflow boundary conditions are applied at the left entry and the right exit respectively. The solution presents several shock waves further interacting with the wall boundaries. Once again we use this test case to compare both the two-point and multi-point schemes on meshes constituted of quadrangles. In figures 17a and 17b, (resp. figures 18a and 18b) we present the results for the two-point schemes (resp. multi-point schemes) employing a structured grid made of N = 334043 quadrilateral cells. The density gradient (numerical Schlieren) and density contours (colors and iso-contours) are both illustrated. The general shape and position of the multiple shocks appear correctly captured by the two schemes. However, we also observe in details that the two-point scheme develops some spurious phenomena along the flat step up to the reflection Y-shaped region, and, along the primary bow shock, especially close to the y = 0 axe. The numerical Schlieren plot enhances these behaviors.

On the contrary the multi-point scheme in figures 18a and 18b does not present this pathological behavior.

If at first glance the two-point scheme may seem less dissipative and simpler in its design. However it reveals to develop unacceptable spurious phenomena which drastically or viciously pollute the obtained numerical solution. Contrarily the multi-point scheme behaves far better and may seem an appropriate robust, entropic and positivity preserving first-order Eulerian scheme on general grids.

Conclusion and perspectives

We have developed a novel subface flux-based multi-dimensional Finite Volume scheme for solving hyperbolic systems of conservation laws over general unstructured meshes. The subface flux approximation relies on simple approximate Riemann solvers. This original FV method is characterized by an explicit time step condition which ensures the preservation of the definition domain. Contrarily to the classical face-based Finite Volume schemes, the notions of conservativity and entropy stability are defined by means of node-based conditions.

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 derive approximate Riemann solvers in the normal direction to the subface by means of the fundamental Lagrange-to-Euler mapping. As such, the Eulerian solvers are constructed directely from their Lagrangian counter-parts. Doing so, we exhibit a fundamental link between Eulerian and Lagrangian numerical fluxes.

The node-based conservation condition implies that the stencil onto which the scheme operates covers all cells in contact with the current one leading to a genuinely multidimensional scheme, hence the name 'multi-point' scheme.

The application to the system of gas dynamics has been thoroughly depicted. First, Lagrangian and Eulerian simple Riemann solvers linked together through the Euler-to-Lagrange mapping have been exhibited. The associated Eulerian scheme is multi-point, conservative, entropic, and positive by construction under a well-defined time step condition. Moreover, the wave-speeds of the Riemann solver are appropriately ordered by construction. This Eulerian numerical scheme has been implemented in 2D in its first-order explicit version on unstructured grids. A list of test cases are presented and the numerical results of the multi-point scheme have been compared to the ones of the classical two-point scheme. We have observed that the multi-point scheme is robust and seems insensitive to spurious numerical phenomena like the carbuncle instability. As such the multi-point scheme would be an appropriate first-order robust candidate scheme to build upon high-order accurate extensions.

In the future, we will investigate the extension to higher-orders of accuracy, namely secondand third-orders. Also, we plan to design and implement similar concepts in a three-dimensional geometry on unstructured grids. For steady aerodynamics applications, we also intend to develop a time marching algorithm relying on an implicit version of the scheme. At last, we shall study the application of this novel methodology more complex physics such as the shallow water equations with source terms and the magneto-hydrodynamics.

Figure 1 :

 1 Figure 1: Geometrical entities attached to the polygonal cell ωc.

Figure 2 :

 2 Figure 2: Interface between cell ωc and cell ω d .

W

  pcf , and assuming that the Riemann solver preserves the domain of definition, we claim that the Finite Volume scheme is D-preserving under the global time-step condition Δt ≤ min c Δt c . (15) 2.5. Entropy inequality Here, we derive the entropy flux approximation attached to our Finite Volume scheme. Assuming that U n c ∈ D, W pcf (ξ) ∈ D and Δt satisfies the time step condition (14), then by virtue of (13), pcf (ξ) dξ, which thus belongs to D. Under the foregoing assumptions and thanks to the convexity of the entropy, Σ n+1 c = Σ(U n+1 c

Figure 3 :

 3 Figure 3: Fragment of the computational grid in the vicinity of point p.

  conservation condition and also in terms of the left and right-sided subface entropy fluxes Q l pf and Q r pf for the entropy condition, refer respectively to sections 2.7 and 2.8. The aforementioned left and right sided subface fluxes and subface entropy fluxes are attached to subface f impinging at p, refer to figure 3b, and defined through the approximate Riemann solver, W l pf , refer to (10) and (19).

  ) Now, following [4] we complete the Lagrangian Riemann solver characterization introducing the intermediate fluxes for s = l, r G n,s = (-v n,s , p s , 0, (pv n ) s ) t . These intermediate fluxes are nothing but the left and right-sided fluxes defined respectively by (43a) and (43b). Therefore G n,l ≡ G n and G n,r ≡ G + n satisfy the system

(

  a) Two-point flux approximation. (b) Multi-point flux approximation.

Figure 4 :

 4 Figure 4: Stencil representation of the Finite Volume scheme in terms of the subface flux approximation over the fragment of a polygonal grid.

  (a) Density contours at t final = 50 of two-point (top) and multi-point (bottom) schemes.

  solver (b) Comparison of the solution deviations ε0.

Figure 5 :

 5 Figure 5: Classical odd-even decoupling problem.

  (a) Density contours at t final = 50 of two-point (top) and multi-point (bottom) schemes.

  solver (b) Comparison of the solution deviations ε0.

Figure 6 :

 6 Figure 6: Modified odd-even decoupling problem.

  (a) Triangular grid, two-point scheme. (b) Triangular grid, multi-point scheme.

Figure 7 :

 7 Figure 7: Hypersonic flow over half cylinder. Density maps over triangular grids for both two-point and multi-point schemes with 20 isolines over the interval [1, 6.2].

  (a) Quadrangular grid, multi-point scheme.(b) Polygonal grid, multi-point scheme.

Figure 8 :

 8 Figure 8: Hypersonic flow over half cylinder. Density maps over quadrangular and polygonal grids for both multipoint scheme with 20 isolines over the interval [1, 6.2].

  Multi-point solver -1250 quads Multi-point solver -5000 quads Multi-point solver -20000 quads (b) Grid convergence.

Figure 9 :

 9 Figure 9: Hypersonic flow over half cylinder. Various plots of the pressure coefficient Cp at the wall with respect to the angular position employing the foregoing grids.

( a )

 a Left half : Two-point scheme, Right half : Multipoint scheme. (b) Left: Two-point, Right: Multi-point.

Figure 10 :

 10 Figure 10: Hypersonic flow over half cylinder. Comparison between density map for the two-point and multi-point schemes over the same polygonal grid and the 3D extruded views for the density field obtained with both schemes.

Figure 17 :

 17 Figure 17: Forward-facing step on a structured grid composed of N = 334043 quadrilateral cells: Numerical results employing the two-point scheme.

  (a) Schlieren plot. (b) Density map.

Figure 18 :

 18 Figure 18: Forward-facing step on a structured grid composed N = 334043 quadrilateral cells: Numerical results employing the multi-point scheme.

xi0,j = x i0,j + 10 -4 (2ζ j -1)

400×400. The resulting scattered plots are displayed in figure 12c, where we observe a convergence of the solution without spurious effects. To assess the capability of the numerical methods in maintaining symmetry for shock wave propagation on irregular grids, we run the Sedov problem on the setup introduced in [START_REF] Dobrev | Curvilinear Finite elements for Lagrangian hydrodynamics[END_REF]. A Cartesian grid is divided into four quadrants with the following grid spacing:

. We once again present the numerical density map for both schemes in figure 13 and the scattered plots for both schemes in figure 14. We observe that the radial nature is preserved with the multi-point scheme even though the shock wave is crossing regions of with different aspect ratios and resolutions. Similar instabilities are observed for the two-point scheme.

Noh implosion problem

The problem designed by Noh [START_REF] Noh | Errors for calculations of strong shocks using an artificial viscosity and an artificial heat flux[END_REF] is a well known test case used to validate numerical schemes in the regime of infinite strength shock wave. The initial computational domain is defined by a 38

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐ The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: