Jérôme Charton

Stefka Gueorguieva

Pascal Desbarats

SINEFITTING: ROBUST CURVATURE ESTIMATION ON SURFACE TRIANGULATION

Keywords: triangular meshes, generalized curvature estimation, surface fitting, curvature measures

A novel algorithm for robust curvature estimation based on sinusoidal curve fitting is proposed. The evaluation of this algorithm is presented on analytical surface triangulations by comparing it with other recognized fitting methods according to three criteria: convergence, precision and robustness. By experimenting on various data, we show that the Sinefitting algorithm is less affected by errors in vertex normal estimation.

RELATED WORK

Knowledge about the geometric shape of an object is based on the understanding of the differential structure of the object boundary surface: the principal curvatures and directions, the Gauss and the Mean curvatures of the boundary surface. We are interested in estimating this differential structure of the underlying smooth surface from a given triangulation. Since the pioneer works of [START_REF] Besl | Invariant surface characteristics for 3d object recognition in range images[END_REF][START_REF] Sander | Generic curvature features from 3-d images[END_REF][START_REF] Flynn | On reliable curvature estimation[END_REF], curvature estimation is a central issue in a great number of research. For surveys we refer the reader to [START_REF] Maltret | Discrete curvatures and applications : a survey[END_REF][START_REF] Petitjean | A survey of methods for recovering quadrics in triangle meshes[END_REF][START_REF] Surazhsky | A comparison of gaussian and mean curvatures estimation methods on triangular meshes[END_REF]. Generalized curvatures, convergence and measure stability are covered by [START_REF] Guoliang | Convergence analysis of a discretization scheme for gaussian curvature over triangular surfaces[END_REF][START_REF] Morvan | Generalized Curvatures[END_REF][START_REF] Chazal | Stability of curvature measures[END_REF]. The approaches considered in this paper are based on the fitting paradigm: first, a quadratic form is constructed at each vertex of the surface triangulation and then the local differential structure is derived. Chen et al [START_REF] Chen | Intrinsic surface properties from surface triangulation[END_REF] use the normal curvature as stated by the Meusnier and Euler's theorem and locally fit a set of circles through the neighbour vertices. McIvor et al [START_REF] Mcivor | A comparison of local surface geometry estimation methods[END_REF] use a quadratic surface fitting for determining the principal frame and the rotated principal quadric. Different variants of this method are applied depending on the type of the quadratic form (simple SQFA, extended or full). More generally, Cazals et al [START_REF] Cazals | Estimating differential quantities using polynomial fitting of osculating jets[END_REF] make use of osculating jets defined as truncated Taylor expansions. For Taubin [START_REF] Taubin | Estimating the tensor of curvature of a surface from a polyhedral approximation[END_REF] the quadratic form is expressed as an integral representation that is used to obtain the curvature tensor. Langer et al [START_REF] Langer | Exact and interpolatory quadratures for curvature tensor estimation[END_REF] use integral representations of the Gauss and the Mean curvatures. Cohen-Steiner et al [START_REF] Cohen | Restricted delaunay triangulations and normal cycle[END_REF] use integrals of specific differential forms on their normal cycles. Fitting enables precise estimation of curvatures but it is very sensitive to the surface discretization and to the distribution of the edge directions in the vertex neighbourhood. Indeed, these methods make use of the surface normal approximation which strongly depends on the regularity of the triangles in this neighbourhood.

In the current paper, a novel method for curvature estimation, called SineFitting, is presented. At each vertex of the surface we construct a sinusoidal curve to fit the directional angles from the target vertex to its neighbours. Therefore, in contrast to Chen's method, there is no need to restrict the choice of neighbours to specific configurations (pairs of geometrically opposite vertices). Moreover, computations in under or oversampled neighbourhoods remain robust while maintaining high precision. The elaborated method has a linear convergence rate and is not acutely affected by errors in estimation of the normal.

OVERVIEW OF THE SINEFITTING METHOD

We refer the reader to [START_REF] Spivak | A Comprehensive Introduction to Differential Geometry[END_REF] for detailed discussion on surface differential geometry. Let S be a surface, P a target vertex from S and N the unit normal vector to S in P. Let τ be the tangent plane of S in P, see Fig. 1. Given a unit vector T in τ, the osculating plane ι through N and T intersects S in a curve c. Let n be the unit normal vector of c in P and  the angle between n and N , )

( max T , T = θ  ) ( ) ( θ •sin² k + θ •cos² k = k min max T (2)
The proposed algorithm, called SineFitting algorithm, is based on (1) and (2) in order to evaluate the principal directions min T and , T max the principal curvatures k min and k max , the Gauss and the Mean curvatures, k G and k H , in P. Our approach is based on a two steps procedure: First, the normal curvature k T is evaluated according to (1) for curves locally fitting the normal sections of the surface. Second, a sinusoidal curve is constructed to approximate the computed values of k T following (2). Principal curvatures and directions are calculated for specified values of the sinusoidal amplitude and frequency. Let T S be a triangulation of S, N(P) a neighbourhood of P, N={P i , PP i  T S , i=0,…,m-1} and k i the normal curvature of S along PP i . Let P i * be the projection of P i on τ, M i the middle point of PP i and M i the median of PP i , shown in Fig. 1. Let us construct a circle σ passing through P and P i , and centred at

O i , O i = N l ∩ m i , where N l is the straight line supporting N . Let n be the unit normal of σ, i i σ PO PO = n
, and k σ be the normal curvature of σ,

i σ PO = k 1 .
Let us now apply (1) to k i and k σ :

)) ( .cos( n , N k = k σ i  (3)
The unit vectors N and n are aligned with opposite directions and thus

1 )) ( cos(   = n , N
. The curvature k i can be expressed then as:

i σ i PO = k = k 1   (4)
Let ψ be the angle  (M i , P, O i ). From the right-angled triangle Δ(M i , P, O i ) it follows that

i i i i PO PP • PO PM = ψ 2 ) cos(  (5)
The value of cos(ψ) can be estimated from the scalar product of N  and i PP and by the substitution in (4) we finally compute the estimation of k i :

² 2 i i i PP PP • N • = k (6)
It should be noticed that the estimation ( 6) is similar to the estimation in [START_REF] Taubin | Estimating the tensor of curvature of a surface from a polyhedral approximation[END_REF]. The above construction gives a geometrical insight of the normal curvature estimation.

The second step of our algorithm includes the construction of a sinusoidal curve that approximates the normal curvature in N(P). Our goal is to rise the constraint of using min T and max T in (2), and fit a linear sinusoidal expression. According to (2) we have:

    i min i max i θ •sin² k + θ •cos² k = k (7)
where

) ( max i i T , T = θ  . We introduce the angle ϕ ϕ =  ( 0 T , max T ) (8)
where 0 T is any vector from N * (P), N * (P)

= {PP i * , i=0, …, m-1}. In (7), i θ is substituted by α i , α i = θ i -ϕ.
Then ( 7) is rewritten as:

) ( ) (   + α •sin² k + + α •cos² k = k i min i max i (9)
Next we substitute ϕ, k min and k max as:

0 2 arctan > a , a b =         (10) ²) ² ( b a c = k min   (11)   b² + a² + c = k max ( 12
)
The equation ( 9) is rewritten as:

c + α b• + α a• = k i i i ) sin(2 ) cos(2 (13)
Now we are ready to compute a sinusoidal curve to approximate in a least square way the estimated from (6) normal curvatures. The intuition of this approach is to first estimate a sine wave shape of arbitrary phase. The peaks of its deviation will correspond to the values of k min and k max . The principal directions min T and max T will correspond to the phases of the leftmost peaks in relation to the origin.

The SineFitting algorithm is given in Algo. 1 [START_REF] Goldman | Curvature formulas for implicit curves and surfaces[END_REF], k max [START_REF] Guoliang | Convergence analysis of a discretization scheme for gaussian curvature over triangular surfaces[END_REF] 14: Calculate max T (8) & ( 10)

                                    * * 2 * 1 m m 2 2 1 1 ... 1 ) 2 sin( ) 2 cos( 1 ... ... 1 ) 2 sin( ) 2 cos( 1 ) 2 sin( ) 2 cos( m k k k = c b a • α α α α α α 13: Calculate k min
15: Calculate c = k + k = k , •k k = k max min H max min G 2

EXPERIMENTAL SETUP

The discussion below is detailed on three particular examples, a sphere, S sphere , a right circular cylinder, S cylinder and a trigonometric bivariate function, S trigonometric defined as follows:

S sphere : x² + y² + z² = 2² (14) S cylinder : x² + y² = 2², -2 ≤ z ≤ 2 (15) S trigonometric : 0.1(cos(xπ) + cos(xπ)) (16)
The full experimentation data set covering all geometric configurations according to the surface classification given in [START_REF] Flynn | On reliable curvature estimation[END_REF], is accessible at http://dept-info.labri.fr/~charton/curvature_analysis/. The surface triangulation in use is the square split sampling defined in [START_REF] Hamann | Curvature approximation for triangulated surfaces[END_REF].

The proposed comparative analysis includes methods acknowledged as representative in both discrete and continuous approaches: the discrete approach formulates a closed form for differential geometry operators that works directly on the discrete representation as the method proposed by Desbrun, Meyer et al [START_REF] Desbrun | Discrete differential-geometry operators for triangulated 2-manifolds[END_REF][START_REF] Meyer | Discrete differential operators for computer graphics[END_REF] and abbreviated as SDA. The continuous approach includes a two stage procedure. First an entity fitting is processed: fitting of surfaces, as the simple quadratic fitting method of McIvor et al [START_REF] Mcivor | A comparison of local surface geometry estimation methods[END_REF] abbreviated SQFA, or fitting of curves, as [START_REF] Chen | Intrinsic surface properties from surface triangulation[END_REF], or fitting of the curvature tensor, as [START_REF] Taubin | Estimating the tensor of curvature of a surface from a polyhedral approximation[END_REF] and [START_REF] Langer | Exact and interpolatory quadratures for curvature tensor estimation[END_REF]. Then fitted entities are "interrogated" in order to evaluate the principal curvatures and directions, the Gauss and the Mean curvatures. We study the performances of the SineFitting method with respect to three criteria: the pointwise convergence, the precision and the robustness.

Pointwise convergence

Let P be a target vertex on the surface S, B(P, r) a ball centred at P with radius r and c B(P, r) an intersection curve, c B(P, r) = S ∩ B(P, r). Let us consider N(P) with central vertex at P and neighbour vertices 5 0 = i i } {P on c B(P, r) . The pointwise convergence tests for the target vertex P consist in checking if the estimated values of the curvatures converge to the exact values when r→0.

Precision

The estimated curvatures are compared with the exact ones computed from the curvature formulas for implicit surfaces given in see in [START_REF] Goldman | Curvature formulas for implicit curves and surfaces[END_REF]. For S sphere curvatures are constant, k min = k max = k H = 0.5, and k G = 0.25. For S cylinder any point not belonging to the bottom and the up circle sides has k min = 0, k max = 0.5, k H = 0.25 and k G = 0. For S trigonometric the symbolic computations are performed using Maple17.

Robustness

In order to test robustness, we investigate four strategies for the sampling of

5 0 = i i } {P
on c B(P, r) inspired from the experiments of [START_REF] Garimella | Curvature estimation for unstructured triangulations of surfaces[END_REF], [START_REF] Gatzke | Curvature maps for local shape comparison[END_REF], [START_REF] Gatzke | Estimating curvature on triangular meshes[END_REF] and [START_REF] Langer | Exact and interpolatory quadratures for curvature tensor estimation[END_REF]:

The regular neighbourhood N reg (P) corresponds to a regular sampling around P and enables to study convergence when no degenerate triangles occur in the vicinity of the target vertex.

In the irregular neighbourhood N irreg (P), pairs of vertices P i and P i+2 , i=0,1,2, are aligned. This geometric configuration corresponds to the "regular vertex" following [START_REF] Langer | Exact and interpolatory quadratures for curvature tensor estimation[END_REF]. Being less constrained than N reg (P), it focuses on the direction distribution and downplays the distances to the neighbours. The regular neighbourhood with angle perturbation N regδAngle (P) is constructed from N reg (P) by displacement of a single vertex P j toward or away from its neighbours on the surface. The regular neighbourhood with distance perturbation N regδDist (P) is constructed from N reg (P) by displacement of a single vertex P j towards or away from P on the surface.

RESULTS

Pointwise convergence

Examples of the pointwise convergence, in P(2, 0, 0) on S sphere are provided in Fig. 3. The vertical axis indicates the variation of the mean curvature k H . The center of the k H variation interval is the exact value of k H . The horizontal axis corresponds to the radius r of B(P, r). According to our experimentation all methods converge to the exact curvature values when the theoretical normal is used no matter the perturbations on the neighbourhood N(P). The rates of convergence are linear except for the SQFA method.

When the normal is approximated with an area weighted normal and the neighbourhood vertices are displaced for N regδDist (P), see Fig. 3(c), only SDA converges to the exact value.

The pointwise convergence in P(2, 0, 0) on S cylinder is shown in Fig. 4 and Fig. 5. For this example, all methods converge to the exact values with the theoretical normal except for the case of N reg (P), see Fig. 4(a). The methods SQFA, Langer's and the SineFitting converge to the exact values when the theoretical normal is used for all types of N ( P ) . The best convergence rate is achieved by the SineFitting method.

The results on S cylinder when the area weighted normal is used do not converge to the exact values when the neighbourhood is perturbed. Moreover, the type of perturbation, depending on the distance between the target vertex and the neighbours, or on the angles adjacent to the target vertex, does not have the same impact on the algorithms. As for example, Chen's and Langer's methods are precise for k H with angle perturbation, N irreg (P), while approximating H k with N regδAngle (P) and N regδDist (P). The extremely erroneous value for k H , k H = 0.277778, is produced by the Taubin's method on N irreg (P).

Precision

The precision is illustrated in Fig. 7 with respect to the Gauss curvature evaluation on S trigonometric . As one can see in Fig. 7(b,c,f), Taubin's, Chen's and Langer's methods lack of precision in some points of inflection. The SineFitting method surpasses in precision the other methods: the error deviation is minimal for the different neighbourhood perturbations.

Robustness

The sensitivity of curvature estimation to the precision of the normal is illustrated in Fig. 6 and the area-weighted ones. SQFA and Langer's methods are the more sensitive to the normal precision as seen in Fig. 6(k)(l). The proposed SineFitting algorithm handles the irregularities in the valence and the shape of the surface triangulation when the theoretical normal is available and remains stable when the normal is approximated.

CONCLUSION AND PERSPECTIVES

In the present article we propose a method for the evaluation of the principal directions and curvatures, the Gauss and the Mean curvatures of surface triangulations. The elaborated method is based on a curvature Sinefitting algorithm and allies the advantages of the curve and surface fitting methods in processing the irregular data sampling. Three fundamental treatments are identified during the curvature evaluation at a target vertex: extraction of the neighbourhood, estimation of the normal and principal directions and curvature computation. Key attention is kept to the role of the normal estimation. A comparative analysis with other widely used methods is provided. Experimental results on chosen data set enhance the SineFitting performances with respect to the pointwise convergence, the precision and the robustness of the computation. One future extension of our work is to improve the normal estimation on meshes. In addition, we intend to extend our method to compute curvature on a wider scale by examining multi-ring neighbourhood. Our goal is to revisit geometric saliency and develop segmentation techniques based on precise and robust curvature estimation.

Figure 1 :

 1 Figure 1: Local surface geometry around point P According to Meusnier & Euler's theorem, the normal curvature T k of S in P along T could be defined as    cos . k = k T (1) where k denotes the curvature of c in P and N • n = ) cos( . While T is rotated around N , an infinite number of normal curvatures k T could be defined. The extreme values of k T , k min and k max , are achieved along the principal directions, min T and . max T Euler's theorem gives the relation between k T , k min , k max and the angle θ ,

Figure 2 : 1 12 :

 2112 Figure 2: SineFitting curvature estimation over one ring neighbourhood

Figure 3 :

 3 Figure 3: Sphere pointwise convergence test for mean curvature computation with area-weighted normal.

Figure 4 :

 4 From a quantitative point of view, by counting the total number of successful tests, SQFA and the SineFitting perform better than the others. The convergence rate makes the difference: the quadratic convergence rate of the SQFA makes it less attractive, favouring the SineFitting method. Cylinder pointwise convergence test for mean curvature computation with theoretical normal.

Figure 5 :

 5 . The regions where the estimated values are equal to the exact ones are coloured in green. With blue and red colours the regions where curvatures are respectively underestimated or overestimated are shown. All examples in Fig.6(a)-(f) use the theoretical normal while those in Fig.6(g)-(l) exploit the area-weighted normal approximation. As one can see the errors in the curvature estimation are structurally related to the quality of the underlying region triangulation, vertex valence for the poles of the sphere, and triangle quality all over the triangulation. The SDA outperforms all algorithms having similar results for both normals, the theoretical .Cylinder pointwise convergence test for mean curvature computation with area-weighted normal.

Figure 6 :Figure 7 :

 67 Mean curvature evaluation based on theoretical and area-weighted normal (Gauss curvature evaluation based on theoretical normal


A C Bayeul-Lainé

P Dupont

L Miccoli

G Cavazzini
email: giovanna.cavazzini@unipd.it

A Dazin

G Pavesi

G Bois

Keywords: Pump, diffuser, PIV measurements, numerical calculations, probes

Introduction

The paper presents the results of numerical analysis on the local and global internal flow behaviour at the inlet of the vaned diffuser of a radial flow pump model, taking into account the effect of fluid leakages for various flow rates and a given rotation of speed. For each flow rate, numerical calculations were performed both with two different boundary conditions:

-Without any leakage effects.

-With calculated leakage effects The numerical simulations were realized with the two commercial codes: i-Star CCM+ 8.02.011 (at LML), ii-CFX 10.0 (at University of Padova). RANS unsteady calculations, with a k- RNG model were performed with Star CCM+. Fully unsteady calculations of the whole pump were done with CFX with DES turbulentmethod (Cavazzini [1]). For each flow rate, different angular positions of the impeller are considered. First part of the paper shows global comparisons between numerical and experimental results already presented in ref [START_REF] Pavesi | Analysis of Rotor-stator interactions effects within the vaned diffuser of a radial flow pump[END_REF][START_REF] Cavazzini | Analysis of the rotor-stator interaction in a radial flow pump[END_REF], for which the effects of fluid leakage due to the gap between the rotating and fixed part of the pump model were found to be interesting to be analysed. Second part is devoted to the local comparisons of flow structures at the inlet section of the diffuser, without and with leakages only with Star CCM+.

Leakage modeling

Two sources of fluid leakage occur (Fig. 1):

-The first one between the inlet and the impeller (leakage1) -The second one between the impeller and the diffuser (leakage 2)

Fig. 1 Numerical model

The two fluid leakages are taken account with the two codes: -In CFX (Fig. 2) , the seal systems both at the impeller inlet and outlet were modelled (At the labyrinth of the seal system close to the impeller inlet the mass flow rates determined from the experimental data was prescribed. At the impeller outlet the leakage mass flow rate was controlled by the known experimental pressure in the large plenums upstream the labyrinth).

FLUID LEAKAGE EFFECT ON ANALYSIS OF A VANED DIFFUSER OF SHF PUMP

Inlet zone

Impeller zone

Diffuser zone zone

Leakage 1

Leakage 2 -In Star CCM+, the sketches of the gaps are really taken account (Fig. 3, [START_REF] Ac Bayeul | Numerical and experimental investigations in a vaned diffuser of SHF impeller : fluid leakage effect,21 ème Congrès Français de Mécanique[END_REF]), a complete meshing was realized including the real geometry of the gaps and external domain far upstream the labyrinth. In this respect, boundary conditions include the relative casing of the model.

Global results

Global results of total theoretical head pump are in good agreement between Euler one dimensional approach and frozen rotor calculation. Unsteady calculation results, performed for the two nondimensional flow rates, give the same level, as it can be seen in figures 4 in which: These tendencies were also observed in previous study [START_REF] Ac Bayeul | Numerical and experimental investigations in a vaned diffuser of SHF impeller : fluid leakage effect,21 ème Congrès Français de Mécanique[END_REF].

Q*=Q/Q i (1)
 is =CQ U 2 2 /2 (2)  t =dp t /(U 2 2 /2) (3)

Local velocity and pressure fields

Among all the numerical calculations that have been performed, only results for two specific flow rates are shown in this paper. These flow rates correspond to a positive and a negative leakage flow between impeller outlet and diffuser inlet. Other constant values such as non-dimensional static pressure and absolute angle are also shown respectively on figures 10 and 11. All the figures represent the diffuser blade to blade inlet section from hub to shroud. For each figure, three different impeller positions, relative to the diffuser leading edge section, are given corresponding to positions P3, P5 to P7 as shown in figure 7. It has been found, from calculation results, that leakage effects do not strongly affect tangential components of absolute velocity. So, only radial and axial velocity components have been chosen instead of tangential ones, in order to compare more clearly both results without and with leakage.

Leakage is found to have an important effect mainly near the diffuser pressure side-shroud corner section. This can be seen on figures 8a and 8b for low mass flow rate corresponding to positive leakage. For high mass flow rate (figures 8c and 8d), the differences are more pronounced near the suction mid span region. Axial velocity components show more local differences as shown in figures 9a to 9d. Effect of impeller blade location on axial velocity differences is more visible near the suction side shroud corner since this location corresponds to low static pressure values especially for low mass flow rates for which incidence angle on the diffuser leading edge is important. For high mass flow rate, the incidence angle becomes negative and consequently, local static pressure is higher near the suction side as shown in figures 10. As a consequence, the values of the absolute angle (figures 11) mainly differ in the suction side region and extend on the whole passage section at low mass flow rate. At high mass flow rate, slight modification occurs at the hub to shroud mid plane near the suction side region.

Conclusion

Fig. 2 Fig. 3

 23 Fig.2 Fluid leakage models in CFX

Fig. 4

 4 Fig.4 Pump non-dimensional head curves Fluid leakages are calculated in each code for three different locations (inlet, hub side and shroud side) for the two codes. In star CCM+ frozen rotor and unsteady calculations are performed. Analytical solutions, based on singularities, are also calculated. All results are presented in Figure 5, fluid leakage are considered positive if they income inside the domain. The values of non-dimensional flow rates of fluid leakages are defined by equation (4).

  can be observed, at the inlet gap, all methods lead to quite same results. In case of Star CCM+, fluid leakages are smaller at high flow rate. Concerning, leakages at hub and shroud sides, the amount of flow rate depends on the flow rate in the pump. At low values of Q*, a positive leakage between impeller and diffuser was observed. At high values the tendency inverses as shown in figure6.

Fig. 5 Fig. 6

 56 Fig.5 non-dimensional fluid leakages

Fig. 7

 7 Fig.7 Impeller blade positions Non-dimensional values of velocity are calculated with equations (5)and (6):

Fig 8a P7Fig 8c P7Fig 8b P7Fig 8d P7Fig 9a P7Fig 9c P7Fig 9b P7Fig 9d P7Fig 10a P7Fig 10c P7Fig 10b P7Fig 10d P7FigP7FigP7FigP7Fig

 8a8c8b8d9a9c9b9d10a10c10b10d Fig 8a Non-dimensional radial Velocity Vr* -Q*=0.766case without leakage


Nina Aguillon

A Non-dissipative Reconstruction Scheme for the Compressible Euler Equations

Keywords: Shock capturing scheme, Euler equations, slowly moving shocks, wall heating phenomenon 2010 Mathematics Subject Classification: 35L65, 35L67, 65M06

We present a finite volume scheme, first on the Burgers equations, then on the Euler equations, based on a conservative reconstruction of shocks inside each cells of the mesh. Its main features are the following:

• the scheme is exact whenever the initial datum is a pure shock, in the sense that the approximate solution is the exact solution averaged over the cells of the mesh;

• in general, the scheme has a very low numerical diffusion and the shocks have a width of one or two cells;

• no spurious oscillations in the momentum appear behind slowly moving shocks, which is not the case in most of the scheme developed so far;

We also present prospective result on the full Euler equations with energy. The wall heating phenomenon, which is an artificial elevation of the temperature when a shock reflects on a wall, is also drastically diminished.

Introduction

When approaching the solution of hyperbolic system, each finite volume scheme introduces a small amount of dissipation. The main effect of this numerical viscosity is to smear out the discontinuities over several cells. But in some cases, even the simplest finite volume scheme, as Lax-Friedrichs's or Godunov's, fails to approach the solution in a satisfactory way. In [START_REF] Quirk | A contribution to the great Riemann solver debate[END_REF], Quirk gives a list of such failings. We recall two of them in the context of gas dynamics. The first one is the apparition of an artificial spike in the momentum, followed by spurious oscillations, when computing a slowly moving shock (i.e. a shock that travels slowly compared to the speed of sound). The second one is the wall heating phenomenon, which is a nonphysical spike on the internal energy when a shock reflects on a solid wall, or when two shocks interact. It is worth noticing that for the slowly moving shocks as for the shock reflections, higher order schemes can lead to less satisfactory results than the first order schemes. In both those cases, the numerical viscosity is held responsible (see [START_REF] Shi | The effects of numerical viscosities. I. Slowly moving shocks[END_REF][START_REF] Arora | On postshock oscillations due to shock capturing schemes in unsteady flows[END_REF] and [START_REF] Quirk | A contribution to the great Riemann solver debate[END_REF][START_REF] Noh | Errors for calculations of strong shocks using an artificial viscosity and an artificial heat flux[END_REF]). Several fixes have been proposed, based on the idea of adding a sufficient amount of dissipation to the scheme, either "by hand" or by interlacing two different schemes, see for example [START_REF] Donat | Capturing shock reflections: an improved flux formula[END_REF].

In this paper, we propose a finite volume scheme that is exact when the initial datum is a pure shock (or contact discontinuity). As a consequence, it approaches exactly slowly moving shocks, and no momentum spikes appear when computing a Riemann problem containing a slowly moving shock. The wall heating phenomenon seems to be eliminated as well. Another nice feature of our scheme is that it reconstructs perfectly isentropic compression. The numerical simulations show that our scheme is of order 2 on Riemann problems containing only shocks, which is, up to our knowledge, a novelty.

The main idea of our scheme is the following. Suppose that the initial datum corresponds to a pure shock located at an interface of the mesh. After one iteration in time, at time ∆t, this shock has moved inside a cell. As finite volume scheme are designed to be conservative, the numerical solution will be average in that cell, creating an intermediate value that does not correspond to a pointwise value of the exact solution. Roughly speaking, our scheme consists in rebuilding the initial shock inside that cell. Then, the reconstructed shocks are advected during a time ∆t and the corresponding fluxes are computed. This type of scheme has been first introduced in [START_REF] Després | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF] on scalar conservation law with an increasing convex flux , and then used to compute non classical shocks in [START_REF] Boutin | Convergent and conservative schemes for nonclassical solutions based on kinetic relations[END_REF]. A general framework for framework and a stability analysis can be found in [START_REF] Lagoutière | Non-dissipative entropy satisfying discontinuous reconstruction schemes for hyperbolic conservation laws[END_REF] and [START_REF] Lagoutière | Stability of reconstruction schemes for scalar hyperbolic conservation laws[END_REF]. Other types of discontinuous reconstruction has been applied to systems, and particularly to gas dynamics, in [START_REF] Alouges | The reservoir technique: a way to make Godunov-type schemes zero or very low diffuse. Application to Colella-Glaz solver[END_REF] and [START_REF] Li | Numerically neither dissipative nor compressive scheme for linear advection equation and its application to the Euler system[END_REF].

Our aim is to generalized this idea to a wider class of hyperbolic equations. In comparison with [START_REF] Lagoutière | Stability of reconstruction schemes for scalar hyperbolic conservation laws[END_REF], we would first like to suppress the assumption of monotonicity on the flux, and allow waves to travel in every direction. We must carefully avoid waves interaction near an interface. This can be done in a very simple way with a moving mesh, where the interfaces are enough bend, so that every wave travel slower than the mesh. The Rusanov scheme and its higher order extensions can be interpreted this way. Second, we attempt to apply this method to system of conservation laws. The solution of the Riemann problem consists in the succession of several waves (shocks, contact discontinuities or rarefactions) and it is a priori unclear how to reconstruct the solution. In every cell, we choose to reconstruct the solution as a single shock or contact discontinuity.

The paper is organized as follows. In the first section, we present the reconstruction scheme on a single conservation law with convex flux, typically the Burgers equation. We describe the three steps of the scheme: reconstruction, advection and march in time, and insist on the necessity to use a moving mesh to avoid wave interaction. In the second section, we extend this scheme to the barotropic Euler equations. The key point is to choose the wave that has to be reconstructed. We present a simple criterion that we believe is relevant and robust. The reconstruction scheme degenerates to Lax-Friedrichs scheme in smooth regions. We explain how to coupled it with higher order schemes to regain accuracy in those regions, and present various numerical tests. Remarkably, the spurious oscillations behind a slowly moving shock are completely eliminated. The third and last section is devoted to prospective results on the Euler equations with energy. We present some simulations of shock reflections, where the wall heating phenomenon is eliminated. This results are encouraging and show that the reconstruction scheme do extend to the full Euler equation. However this simulations are done with a rather abrupt detection of the contact discontinuities, and could probably be improved with a smarter criterion to decide which wave should be reconstructed. In a forthcoming version of this paper, we will discuss the two dimensional extension on the shallow water equation.

Acknowledgments

The author heartily thanks Frédéric Lagoutière for his enlightening advice and guidance during the achievement of this work.

Presentation of the scheme on scalar conservation law

We consider the Cauchy problem for a scalar conservation law

∂ t u + ∂ x f (u) = 0 for t ∈ R + and x ∈ R u |t=0 = u 0 . ( 1
)
where the flux f is C 1 (R) and the initial data u 0 are in L ∞ (R). We are interested in entropy solution, i.e. weak solutions of (1) that verify the additional entropy inequality

∂ t S(u) + ∂ x G(u) ≤ 0 for t ∈ R + and x ∈ R. (2)
The functions (S, G) are an entropy-entropy flux pair: both S and G are C 1 , S is convex and G ′ = S ′ f ′ . If f and S are strictly convex, then this solution is unique, as proven for example in [START_REF] Krushkov | First-order quasilinear equations in several independent variables[END_REF][START_REF] Godlewski | Hyperbolic systems of conservation laws, volume 3/4 of Mathématiques & Applications[END_REF].

Let us now consider a finite volume scheme for (1). We denote by t 0 = 0 < t 1 < • • • < t n < • • • the successive times where the solution is computed, and by ∆t k = t k+1 -t k the corresponding time steps. We subdivide the space R in cells having all the same size ∆x. We denote by x j = j∆x the centers of the cells and by x j+1/2 = x j + ∆x 2 their extremities. Finally, we denote by u n j an approximation of the exact solution at time t n , averaged on cell j:

u n j ≈ 1 ∆x x j+1/2 x j-1/2 u(t n , x)dx.
Integrating (1) on the slab [t n , t n+1 ] × [x j-1/2 , x j+1/2 ], we obtain the family of finite volume schemes

   u n+1 j = u n j -∆t n ∆x f n j+1/2 -f n j-1/2 , n ≥ 0, u 0 j = 1 ∆x x j+1/2 x j-1/2 u 0 (x)dx, (3)
where the numerical flux f n j+1/2 is an approximation of

t n+1 t n f (u(t, x j+1/2
))dt that characterized the scheme. It remains to explain how the fluxes are chosen in our scheme. There is not much difference with the case of linear advection presented, in term of discontinuous reconstruction scheme, in [START_REF] Boutin | Convergent and conservative schemes for nonclassical solutions based on kinetic relations[END_REF]. The idea is depicted on Figure 1. We consider that the value u n j in cell j is actually coming from the averaging of a shock between u n j-1 and u n j+1 located somewhere inside the cell. Starting from u n j we can reconstruct this discontinuity. By conservation of mass, it should lie at a distance

d n j = ∆x u n j+1 -u n j u n j+1 -u n j-1
of the left interface. It falls inside the cell if and only if u n j takes value between its neighbors u n j+1 and u n j-1 . We do not want to reconstruct a shock that is not entropy satisfying. Therefore, as f is convex, we accept the reconstruction only if u n j-1 ≥ u n j+1 . We now let the reconstructed discontinuities evolve during a time ∆t and compute the fluxes at each interface. Suppose for a minute that min x∈R f ′ (u 0 (x)) > 0. Then every reconstructed shock travels to the right and, if two such shocks interact, the resulting waves will also travel to the right. As a consequence, if we perform a reconstruction in cell j, we know that the reconstructed shock will cross the right interface at time

u j u j+1 u j-1 u j,R u j,L j j + 1 j -1 d j ∆x -d j
∆t j+1/2 = ∆x -d n j σ n j where σ n j = f (u n j-1 )-f (u n j+1 ) u n j-1 -u n j+1
is the speed of the reconstructed shock. As everything is coming from the left, the flux passing through the j + 1/2 interface is f (u n j+1 ) before the crossing time and f (u n j-1 ) afterwards (if the crossing time is smaller than the time step ∆t). Suppose now that it exists a sonic point, i.e. that f ′ (u 0 ) changes sign. If u 0 is decreasing, we can reconstruct in every cell with entropy satisfying shocks. Then by convexity of f , σ n j is bigger than σ n j+1 and it is possible that the two waves interact near their common interface in x j+1/2 , as depicted on Figure 2. Moreover, the resulting wave can travel in any direction, and it becomes impossible to compute the flux through this interface without solving the new Riemann problem created by the interaction of the waves.

j j j + 1 j + 1 σ n j σ n j+1
Figure 2: Interaction of waves near an interface. On the left, the bold lines represent the reconstruction and the arrows the speed of the shocks. On the left, the auto similar waves created by the reconstructed shocks are drew in the (x, t) plane.

To overcome this difficulty we use a moving mesh. The property of propagation at finite speed states that the speed of waves are smaller than V waves = max j∈Z |f ′ (u j )|. We introduce a mesh speed V mesh ≥ V waves and bend the interface according to it. It follows that the shock reconstructed in cell C j will cross the left interface. Moreover, even though two adjacent waves interact, the resulting waves will also travel slower than V mesh . Therefore, they will not catch up with the interface, as depicted on Figure 3. At the next time iteration, the mesh is bent in the other direction so that the space grid is unchanged after two time steps.

Figure 3: A mesh with interfaces traveling faster than the maximum of wave speeds. Two interactions between adjacent waves occur, but the resulting perturbations of the solution do not cross the interfaces. This use of a moving mesh makes possible the computation of the flux without resolving any wave interaction. The very same idea is the key stone of the central schemes. Indeed, if |V mesh | > V waves , the evaluation of the solution of Riemann problems along the x = V mesh t becomes a trivial task. This is used in the Rusanov scheme, where V mesh is locally computed, and in its higher order extension initiated in [START_REF] Nessyahu | Nonoscillatory central differencing for hyperbolic conservation laws[END_REF] .

Reconstruction scheme

Summarizing the previous section, we split the reconstruction scheme into 3 steps.

Step 1: Reconstruction We first try to reconstruct, for every j ∈ Z, the value u n j as the average of an entropy satisfying shock between u n j-1 and u n j+1 . If it is not possible, nothing is done. In other words, we compute the distance to the left interface of the cell by conservation of mass

d n j = ∆x u n j+1 -u n j u n j+1 -u n j-1 , (4)
and we accept the reconstruction if and only if it corresponds to an entropy satisfying shock lying inside the cell. We introduce the reconstructed states

u n j,L = u n j-1 and u n j,R = u n j+1 if 0 < d n j < ∆x and u n j-1 > u n j+1 , u n j,L = u n j and u n j,R = u n j , otherwise , (5)
and the speed associated to the reconstructed shock

σ n j = f (u j,L )-f (u j,R ) u j,L -u j,R if 0 < d n j < ∆x and u n j-1 > u n j+1 , f ′ (u j ) otherwise. ( 6
)
Step 2: Computation of the fluxes We suppose at time t n , the mesh speed V n mesh is given such that

|V n mesh | ≥ V n waves = max j∈Z |f ′ (u n j )|. (7)
In practice, V n mesh will change of sign at each iteration. The finite volume scheme is now obtained by integrating (1) on the slab

{(t, x) : t n ≤ t ≤ t n+1 , x n j-1/2 + tV n mesh ≤ x ≤ x n j+1/2 + tV n mesh },
where the extremities of the space cells are now defined from time to time by x n+1 j+1/2 = x n j+1/2 + V n mesh ∆t. Therefore f n j+1/2 is now an approximation of

t n+1 t n f (u(t, x n j+1/2 + V n mesh t)) -V n mesh u(t, x n j+1/2 + V n mesh t)dt.
Our construction of the numerical fluxes is as follow. If V n mesh < 0, the reconstructed shocks cross the right interface x = x n j+1/2 + V n mesh t. The crossing time ∆t n+1/2 j+1/2 verifies

x n j-1/2 + d j + σ n j ∆t n+1/2 j+1/2 = x n j+1/2 + V n mesh ∆t n+1/2 j+1/2 ,
and is therefore equals to

∆t n+1/2 j+1/2 = ∆x -d n j σ n j -V n mesh .
Consequently, we set

∆tf n j+1/2 = min(∆t, ∆t n+1/2 j+1/2 ) f (u n j,R ) -V n mesh u n j,R + (∆t -min(∆t, ∆t n+1/2 j+1/2 )) f (u n j,L ) -V n mesh u n j,L . (8)
If V n mesh > 0, then the reconstructed shocks cross the left interface and we denote by

∆t n+1/2 j-1/2 = d j V n mesh -σ n j
the crossing time. The flux is given by

∆tf n j-1/2 = min(∆t, ∆t n+1/2 j-1/2 ) f (u n j,L ) -V n mesh u n j,L + (∆t -min(∆t, ∆t n+1/2 j-1/2 )) f (u n j,R ) -V n mesh u n j,R . (9)
Step 3: March in time We update compute u at the next time by applying

u n+1 j = u n j - ∆t n ∆x f n j+1/2 -f n j-1/2 , n ≥ 0.
Remark 1.1. When no reconstruction is performed, either because u n j-1 and u n j+1 are linked by a rarefaction wave, or because d j does not belong to [0, ∆x], the value of d j does not matter. The flux degenerates toward a simple left or right decentering

f n j+1/2 = f (u n j ) -V n mesh u n j if V n mesh < 0 f (u n j+1 ) -V n mesh u n j+1
if V n mesh > 0 which also corresponds to the fluxes of the staggered Lax-Friedrichs scheme, which is convergent under Hypothesis (7).

Proposition 1.2. The reconstruction scheme is consistent: if

u n j = u n j+1 , then f n j+1/2 = f (u n j ) -V n mesh u n j . Proof. Indeed, if u n j = u n j+1 , no reconstruction is performed.
The stability and convergence of reconstruction scheme for scalar conservation law has recently been studied in [START_REF] Lagoutière | Non-dissipative entropy satisfying discontinuous reconstruction schemes for hyperbolic conservation laws[END_REF], where two extensions of the scheme inside rarefaction waves are also proposed.

Remark 1.3. The use of a staggered mesh imposes the CFL condition

∆t n ≤ ∆x |V n mesh | + V waves .
Proposition 1.4. When the initial data correspond to a pure shock, i.e. has the form

u |t=0 (x) = u l 1 x<0 + u r 1 x>0 ,
with u l > u r , then the reconstruction scheme is exact, in the sense that

u n j = 1 ∆x x j+1/2 x j-1/2 u exa (n∆t, x)dx.
The function

u exa (t, x) = u l 1 x< f (u l )-f (ur ) u l -ur t + u r 1 x> f (u l )-f (ur ) u l -ur t
is the exact solution.

Proof. The scheme is exactly build to that purpose. For the sake of simplicity we suppose that u l > u r > 0. Therefore for every u ∈ [u r , u l ], every wave generated by a Riemann problem between u l an u or between u and u r has a positive speed and it is not necessary to use a moving mesh. We set V n mesh = 0 for all n. Let us proceed by induction. Suppose that for some n ≥ 0, for every j ∈ Z,

u n j = 1 ∆x x j+1/2 x j-1/2 u exa (n∆t, x)dx.
It is true for n = 0 as soon as the initial sampling is done correctly. We denote by j 0 the cell where the shock is located and by δ ∈ [0, ∆x] its distance to x j 0 -1/2 . We have

     u n j = u l if j < j 0 , u n j 0 = δ ∆x u l + ∆x-δ ∆x u r if j = j 0 , u n j = u r if j > j 0 .
As a consequence, the scheme reconstructs the shock in cell j 0 with d j 0 = δ, and it has the correct speed

σ j 0 = f (u r ) -f (u l ) u r -u l .
Moreover, as u n j 0 -2 = u n j 0 -1 and u n j 0 +1 = u n j 0 +2 , no reconstruction is performed in cells j 0 -1 and j 0 + 1, and

f n j 0 -1/2 = f (u l ) = u 2 l 2 , while f n j 0 +3/2 = f (u r ) = u 2 r 2 .
We now compute the flux f n j 0 +1/2 . The reconstructed discontinuity crosses the right interface at time

∆t n+1/2 j 0 +1/2 = (∆x -δ) σ j 0 .
If ∆t j 0 +1/2 is smaller than ∆t, then Formula (8) rewrites

∆tf n j 0 +1/2 = ∆t n+1/2 j 0 +1/2 f (u r ) + (∆t -∆t n+1/2 j 0 +1/2 )f (u l ) = (∆x -δ) σ j 0 [f (u r ) -f (u l )] + ∆tf (u l ) = (∆x -δ)(u r -u l ) f (u r ) -f (u l ) [f (u r ) -f (u l )] + ∆tf (u l ).
We obtain

f n j 0 +1/2 = ∆x -δ ∆t (u r -u l ) + f (u l ),
and therefore we have

u n+1 j 0 = u n+1 j 0 - ∆t ∆x (f n j 0 +1/2 -f n j 0 -1/2 ) = δ ∆x u l + ∆x -δ ∆x u r - ∆t ∆x ∆x -δ ∆t (u r -u l ) + f (u l ) -f (u l ) = u l
and, if we denote by δ ′ the quantity δ + ∆tσ j 0 -∆x,

u n+1 j 0 +1 = u n+1 j 0 +1 - ∆t ∆x (f n j 0 +3/2 -f n j 0 +1/2 ) = u r - ∆t ∆x f (u r ) - ∆x -δ ∆t (u r -u l ) + f (u l ) = δ ′ ∆x u l + ∆x -δ ′ ∆x u r .
The crossing time ∆t n+1/2 j 0 +1/2 is smaller than ∆t if and only if δ + σ j 0 ∆t is larger than ∆x. Therefore, in that case at time n + 1, the shock has moved inside the cell j 0 + 1 and its distance to the left interface is exactly δ ′ . On the other hand, when ∆t n+1/2 j 0 +1/2 ≥ 0, the shock is still in cell j 0 at time n + 1. Its distance to the left interface is δ ′ = δ + ∆tσ j 0 . In that case, the flux f n j 0 +1/2 is simply f (u r ), and we easily obtain u n+1 j 0 +1 = u r and

u n+1 j 0 = u n+1 j 0 - ∆t ∆x (f n j 0 +1/2 -f n j 0 -1/2 ) = δ ∆x u l + ∆x -δ ∆x u r - ∆t ∆x (f (u r ) -f (u l )) = δ ′ ∆x u l + ∆x -δ ′ ∆x u r ,
which concludes the proof.

Numerical test: isentropic compression

We perform a numerical test that exhibit the non diffusive behavior of the reconstruction scheme. We consider the Burgers equation f (u) = u 2 2 with the initial data

u 0 (x) =      3 if x ≤ -3; 3 -(x + 3) if -3 ≤ x ≤ -1; 1 if 1 ≤ x.
The exact solution is given by

u(t, x) = 3 * 1 x≤-3+3t + 3 -x-(3+3t) (1-t) 1 -3+3t<x<-1+t + 1 -1+t≤x if t < 1 3 * 1 x≤2(t-1) + 1 2(t-1)<x if t ≥ 1 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 1 1.5 2 2.5 3 Reconstruction Godunov Exact T=0 T=0.5 T=1 T=1.5
Figure 4: Simulation of an isentropic compression with the Godunov scheme and the reconstruction scheme

We compare the Godunov scheme and the reconstruction scheme on the interval [-4, 2] discretized with 100 cells. The CFL number is 0.4. As depicted on Figure 4, the shock is perfectly reconstructed at time T = 1 with the reconstruction scheme, while it is smeared out on 10 cells by the Godunov's scheme. After time T = 1, the shock is exactly advected by the reconstruction scheme, while the Godunov scheme's is still diffusing it. The apparition of steps inside the isentropic compression is very similar to those which appears in the smooth regions in [START_REF] Després | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF] and [START_REF] Lagoutière | Non-dissipative entropy satisfying discontinuous reconstruction schemes for hyperbolic conservation laws[END_REF]. It does not affect the order of convergence of the scheme. Entropy fixes have been proposed in [START_REF] Lagoutière | Non-dissipative entropy satisfying discontinuous reconstruction schemes for hyperbolic conservation laws[END_REF], and no such steps appear in [START_REF] Li | Numerically neither dissipative nor compressive scheme for linear advection equation and its application to the Euler system[END_REF] 2 Extension to the isothermal Euler equation

Adaptation of the reconstruction scheme

This section is devoted to the generalization of the reconstruction scheme to a particular system of conservation laws, the isothermal Euler equations. It describes the evolution of an inviscid compressible gas, having density ρ and velocity u, when the pressure law is p = p(ρ) = c 2 ρ. Here, c is the speed of sound. This pressure law to avoids us to deal with vacuum, but the scheme extend straightforwardly to pressure laws p(ρ) = aρ γ with a > 0 and γ > 1 when there is no vacuum. The Cauchy problem for this system writes

     ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + c 2 ρ) = 0, ρ |t=0 = ρ 0 and u |t=0 = u 0 . ( 10
)
This system is strictly hyperbolic, the eigenvalues of the jacobian matrix being u -c and u + c. The corresponding fields are genuinely nonlinear. This system has been widely studied. We are here interesting in the structure of the Riemann problem that can be found, for example, in [START_REF] Evans | Partial differential equations[END_REF]. We recall that a Riemann problem is (10), with initial data

ρ 0 (x) = ρ L 1 x<0 + ρ R 1 x>0 , u 0 (x) = u L 1 x<0 + u R 1 x>0 .
The solution consists of the succession of two waves, rarefactions or shocks, separated by an intermediate state. As a consequence, unless we start with very specific initial data, there is no reason for states (ρ n j-1 , u n j-1 ) and (ρ n j+1 , u n j+1 ) to be linked by an entropy satisfying shock, and the Riemann problem between these two states is more likely to develop the full pattern of two waves. How to decide for a reconstruction procedure in this framework? We decided to try to reconstruct the state (ρ n j , u n j ) as the average of one of the waves appearing in the Riemann problem between its neighbor states (ρ n j-1 , u n j-1 ) and (ρ n j+1 , u n j+1 ). This wave must be an entropy satisfying shock and must somehow be the prevailing wave in the Riemann problem. We use the following lemma to decide wether we try to reconstruct and which wave is chosen.

Lemma 2.1. If u L > u R and ρ L < ρ R , then the Riemann problem between (ρ L , u L ) and (ρ R , u R ) contains a 1-shock . If u L > u R and ρ L > ρ R , it contains a 2-shock.
Proof. Let us denote by (ρ I , u I ) the intermediate state of this Riemann problem. As the velocity increases through rarefaction waves, the Riemann problem contains at least one shock. Moreover, the density decreases through 1-rarefaction waves and 2-shocks, and increases through 1-shocks and 2-rarefaction waves. Hence, if ρ L < ρ R , the Riemann problem contains a 1-shock, while if ρ L > ρ R , it contains a 2-shock.

The scheme is built as follow.

Step 1: Reconstruction We decide wether we should try to reconstruct in cell j using Lemma 2.1. We denote by (ρ * j , u * j ) the intermediate state in the Riemann problem between (ρ n j-1 , u n j-1 ) and (ρ n j+1 , u n j+1 ). The candidates for the reconstructed states are:

(ρ n j,L , ūn j,L ) =      (ρ n j-1 , u n j-1 ) if u n j-1 > u n j+1 and ρ n j-1 < ρ n j+1 ; (ρ * j , u * j ) if u n j-1 > u n j+1 and ρ n j-1 > ρ n j+1 ; (ρ n j , u n j ) otherwise; (11)
and

(ρ n j,R , ūn j,R ) =      (ρ * j , u * j ) if u n j-1 > u n j+1 and ρ n j-1 > ρ n j+1 ; (ρ n j+1 , u n j+1 ) if u n j-1 > u n j+1 and ρ n j-1 < ρ n j+1 ; (ρ n j , u n j ) otherwise. ( 12
)
It corresponds respectively to a 1-shock reconstruction, a 2-shock reconstruction and no reconstruction.

Remark 2.2. If by extraordinary (ρ n j-1 , u n j-1 ) and (ρ n j+1 , u n j+1 ) were linked by a single shock, we will have

(ρ n j,L , ūn j,L ) = (ρ n j-1 , u n j-1 ) and (ρ n j,R , ūn j,R ) = (ρ n j+1 , u n j+1 ).
This is crucial to approach exactly pure shocks.

We denote by σ j the associated speed of the shock, set arbitrarily to 0 when no reconstruction is performed:

σ j =            u n j-1 -c ρ * j ρ n j-1 if u n j-1 > u n j+1 and ρ n j-1 > ρ n j+1 ; u n j+1 + c ρ * j ρ n j+1 if u n j-1 > u n j+1 and ρ n j-1 < ρ n j+1 ; 0 otherwise; (13)
We then compute the distances to the left interface. We now have two distances, one for the conservation law on the density ρ and one for the conservation law on the momentum q = ρu:

d n,ρ j = ∆x ρn j,R -ρn j ρn j,R -ρn j,L and d n,q j = ∆x qn j,R -qn j qn j,R -qn j,L . (14)
We insists on the fact that even though we reconstruct a different discontinuity in ρ and q, it has the same speed σ n j in the two variables. We must now decide whether the reconstruction is accepted or not. The most natural choice is to accept the reconstruction when 0 < d n,ρ j < ∆x and 0 < d n,q j < ∆x, i.e. when we are able to place the discontinuity inside the cell, such that both ρ and q are conserved in the cell. This choice will be referred to as the fully conservative reconstruction. However, we also consider the case where the reconstruction is accepted whenever 0 < d n,ρ j < ∆x. No condition is required on d q j,n , which means that the reconstruction on the momentum can be nonconservative. This choice will be referred to as the half conservative reconstruction. The numerical tests below will show the interest of this less severe choice. Definition 2.3. For the isothermal Euler equation (10), the left and right reconstructed states are:

• for the fully conservative reconstruction,

(ρ n j,L , u n j,L , ρ n j,R , u n j,R ) = (ρ n j,L , ūn j,L , ρn j,R , ūn j,R ) if 0 < d n,ρ j < ∆x and 0 < d n,q j < ∆x; (ρ n j , u n j , ρ n j , u n j ) otherwise. (15)
• for the half conservative reconstruction,

(ρ n j,L , u n j,L , ρ n j,R , u n j,R ) = (ρ n j,L , ūn j,L , ρn j,R , ūn j,R ) if 0 < d n,ρ j < ∆x; (ρ n j , u n j , ρ n j , u n j ) otherwise. ( 16
)
where ρn j,L , ūn j,L , ρn j,R and ūn j,R are defined in (11,12), and d n,ρ j , d n,q j are defined in (14).

Step 2: Computation of the fluxes

The fluxes are computed exactly as in the previous section, Equations ( 8) and (9). The only difference is that in general, d n,ρ j = d n,q j and we now have two crossing times

     ∆t ρ j+1/2 = ∆x-d n,ρ j σ j -V mesh ; ∆t q j+1/2 = ∆x-d n,q j σ j -V mesh ; or      ∆t ρ j-1/2 = d n,ρ j V mesh -σ j ; ∆t q j-1/2 = d n,q j V mesh -σ j ,
for V mesh < 0 and V mesh > 0 respectively. The fluxes are given by

             ∆tf n,ρ j+1/2 = min(∆t, ∆t ρ j+1/2 ) f ρ (ρ n j,R , q n j,R ) -V mesh ρ n j,R
+ (∆t -min(∆t, ∆t ρ j+1/2 )) f ρ (ρ n j,L , q n j,L ) -V mesh ρ n j,L ; ∆tf n,q j+1/2 = min(∆t, ∆t q j+1/2 ) f q (ρ n j,R , q n j,R ) -V mesh q n j,R + (∆t -min(∆t, ∆t q j+1/2 )) f q (ρ n j,L , q n j,L ) -V mesh q n j,L ,

when V mesh is positive, and

             ∆tf n,ρ j-1/2 = min(∆t, ∆t ρ j-1/2 ) f ρ (ρ n j,L , q n j,L ) -V mesh ρ n j,L + (∆t -min(∆t, ∆t ρ j-1/2 )) f ρ (ρ n j,R , q n j,R ) -V mesh ρ n j,R ; ∆tf n,q j-1/2 = min(∆t, ∆t q j-1/2 ) f q (ρ n j,L , q n j,L ) -V mesh q n j,L + (∆t -min(∆t, ∆t q j-1/2 )) f q (ρ n j,R , q n j,R ) -V mesh q n j,R ; (18)
when V mesh is negative. Here, f ρ and f q denote the two components of the flux: f ρ (ρ, q) = q and f q (ρ, q) = q 2 ρ + c 2 ρ. Step 3: March in time Eventually, the conservative variables are updated to the next time step:

ρ n+1 j = ρ n j -∆t ∆x (f n,ρ j+1/2 -f n,ρ j-1/2 ), q n+1 j = q n j -∆t ∆x (f n,q j+1/2 -f n,q j-1/2 ). ( 19
)
Proposition 2.4. The reconstruction scheme is consistant, and exact on pure shocks.

Proof. We proceed once again by induction, and focus on the case of a 1-shock. We denote by (ρ L , q L ) (resp. (ρ R , q R )) the left (resp. the right) density and momentum. Suppose that at the n-th iteration, the scheme gave the exact average of the solution. Denote by j 0 the cell where the shock lies, and by δ its distance to the right interface of the j 0 -th cell. In other words, we have

     ρ n j = ρ L if j < j 0 , ρ n j 0 = δ ∆x ρ L + ∆x-δ ∆x ρ R if j = j 0 , ρ n j = ρ R if j > j 0 , and
     q n j = q L if j < j 0 , q n j 0 = δ ∆x q L + ∆x-δ ∆x q R if j = j 0 , q n j = q R if j > j 0 .
Let us first check that u n j 0 -1 ≥ u n j 0 ≥ u n j 0 +1 . This insures that a 1-shock is detected in cell j 0 with Lemma 2.1. We denote by β the quantity α ∆x . We have

u j 0 -u L = βq L + (1 -β)q R βρ L + (1 -β)ρ R - q L ρ L = ρ L (βq L + (1 -β)q R ) -q L (βρ L + (1 -β)ρ R ) ρ L (βρ L + (1 -β)ρ R ) = (1 -β)(ρ L q R -ρ R q L ) ρ L (βρ L + (1 -β)ρ R ) = (1 -β)ρ L ρ R (u R -u L ) ρ L (βρ L + (1 -β)ρ R ) ≤ 0,
and we obtain similarly that u j 0 is bigger than u R . Therefore the 1-shock is detected, and with no change compared to the scalar case, its position and speed are correctly reconstructed. It remains to prove that no reconstruction is performed on cells j 0 -1 and j 0 + 1. In both those cells, Lemma 2.1 detects a 1-shock. Therefore in cell j 0 -1, the left reconstructed density is ρ L and, by conservation of ρ, no reconstruction can be performed in this cell. In cell j 0 + 1, no reconstruction is performed if the Riemann problem between (ρ j 0 , q j 0 ) and (ρ R , q R ) contains a 1-shock and a 2-rarefaction. Indeed in that case, the density is increasing, and both the left and the right reconstructed densities are smaller than ρ j 0 +1 = ρ R , which blocks the reconstruction. Let us prove that this Riemann problem cannot contains a 1-shock and a 2-shock. If it was possible, the middle state (ρ * , q * ) would belong to the light gray region in Figure 5. Indeed, we would have u * > u R and ρ * > ρ R . On the other hand, (ρ * , q * ) is also on the 1-wave curve of (ρ j 0 , q j 0 ). This state lies on the chord joining (ρ L , q L ) and (ρ R , q R ), and hence the 1-wave curve of (ρ j 0 , q j 0 ) (dark gray curve in the Figure 5) is below the 1-wave curve of (ρ L , q L ). This last curve does not pass through the light gray area, because it is concave and its slope in ρ R is smaller than u R , and we obtain a contradiction.

(ρ L , q L ) (ρ j 0 , q j 0 ) (ρ R , q R ) (ρ * , q * ) ρ q Figure 5:
The Riemann problem between cells j 0 and j 0 + 2 cannot consist of a 1-shock and a 2-shock Once we notice that at each time step, the two reconstructed distances d n,ρ j and d n,q j are equal, the remaining of the proof is identical to that of the scalar case. The data on top of Figure 6 are prepared to correspond to a 1-shock moving at speed 0.1:

u L = 0.1 + c * ρ R
ρ L and u R = 0.1 + c * ρ L ρ R , which yields q L ≈ 2.4361 and q R ≈ 3.2361. It crosses a cell in approximately 6.25 iterations. Inthe bottom of Figure 6, the momentum is modified: q L ≈ 2.6361 and q R ≈ 1.2361. When the initial data consist of a Riemann problem only containing shocks, even if one of them is a slowly moving one, the reconstruction scheme is of order 2. This is illustrated on the left of Figure 7. On the right of this figure, we can see that the error is entirely due to the first iterations in time. Once the two shocks are separated, the shocks are better approximated, and no diffusion appears. The data are the same than in the bottom of Figure 6. An explanation of the superiority of the half conservative reconstruction scheme over the other schemes (including the fully conservative reconstruction scheme) might be linked to the study of slowly moving shocks by Jin and Liu in [START_REF] Shi | The effects of numerical viscosities. I. Slowly moving shocks[END_REF]. In this paper, they proved with a traveling wave analysis that a momentum spike appears in the viscous shock profile for the Euler equations with a linear viscosity

∂ t ρ + ∂ x (ρu) = ǫ∂ xx ρ, ∂ t (ρu) + ∂ x (ρu 2 + c 2 ρ) = ǫ∂ xx q,
when the density has a monotonous profile. When the viscosity is Navier-Stokes-like, the momentum remains monotonous . We believe that the half conservative reconstruction scheme is the only one that has a numerical viscosity looking like the physical Navier-Stokes' viscosity. Indeed, as the reconstruction in the density is performed whenever it is possible (while there is an additional constraint on the momentum in the fully conservative reconstruction), the numerical viscosity is likely to be zero on the mass conservation law as often as possible. Another argument is that if the momentum spike appears in the fully conservative reconstruction scheme, it blocks the reconstruction and the scheme will later behave like the Lax-Friedrichs scheme, while the half conservative reconstruction scheme is more flexible, and will continue to reconstruct near the shock.

Use of other schemes in the rarefaction waves

When no shocks are detected, or when the reconstruction is not accepted, the flux degenerates toward the Lax-Friedrichs flux. As a consequence, it is very diffusive inside the rarefaction waves. A simple cure is to the replace, whenever it is chosen, the Lax-Friedrichs flux by a more accurate one. The fluxes write (for the half conservative scheme):

f n j+1/2 = f n,REC j+1/2 if (0 < d n j < dx and V n mesh < 0) or (0 < d n j+1 < dx and V n mesh > 0), f n,N T j+1/2
otherwise, where f n,REC j+1/2 is defined by ( 17) and (18), while f n,N T j+1/2 is another flux on a staggered grid. In our implementation we use the simplest version of the Nessyahu and Tadmor scheme [START_REF] Nessyahu | Nonoscillatory central differencing for hyperbolic conservation laws[END_REF]. There is no Riemann problems to solve, so it almost does not impact the computation time. As expected, this coupled scheme behaves like the Nessyahu-Tadmor scheme in the rarefaction wave, and like the reconstruction scheme on the shock. This is illustrated by the Figure 8 below.

Prospective results on the full Euler equations

Let us now consider the full gas dynamics equation for an ideal gas:

     ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p) = 0, ∂ t E + ∂ x (u(E + p)) = 0,
where ρ is the density, u is the fluid's velocity and E is the total energy per unit volume. The pressure p writes p = p(ρ, e) = e(γ -1)ρ, where e = E ρ -1 2 u 2 is the specific internal energy and γ is the ratio of specific heat. The complete solution of the Riemann problem can be found, for example, in [START_REF] Eleuterio | Riemann solvers and numerical methods for fluid dynamics[END_REF]. The novelty is that the 2-wave is a contact discontinuity. We denote by U * ,L/R j = (ρ * ,L j , ρ * ,R j , u * j , p * j ) the left and right density, velocity and pressure appearing around the contact discontinuity in the Riemann problem between U j-1 = (ρ n j-1 , u n j-1 , p n j-1 ) and U j+1 = (ρ n j+1 , u n j+1 , p n j+1 ). Let us denote by C CF L the CFL number. The wave to be reconstruct is chosen as follow.

• If u n j-1 ≥ u n j+1 , ρ n j-1 ≤ ρ n j+1 , p n j-1 ≤ p n j+1 and if |ρ j-1 -ρ * ,L j | > C CF L max(|ρ * ,L j -ρ * ,R j |, |ρ * ,R j -ρ j+1 |),
with is a criterion on the force of the wave. Note that with this criterion, some shocks are not detected, while without it, some contact discontinuities are detected as shocks. We observe numerically that using the CFL number tunes the scheme pretty well. Instead of this criterion, we also tried a "two shots reconstruction", in which we first try to reconstruct everywhere with contact discontinuities. If this reconstruction fails (i.e. if we do not have 0 < d ρ j < ∆x and 0 < d E j < ∆x), then we attempt to reconstruct 1-and 3-shocks using Lemma 2.1. The numerical results are very similar to those presented below.This very same difficulty to detect contact discontinuities led us to add more constraints to accept a reconstruction. The first one is similar to the previous 2 × 2 system. We do not impose a conservation constraint on the momentum to mimic the Navier Stockes viscosity. The density and velocity are monotonous along the profile of a viscous shock, which justify the second constraint. Eventually, the last constraint forces the pressure to remain positive.

We perform various numerical tests to compare the reconstruction scheme with other classical schemes. On all these figures, Rec indicates the half conservative reconstruction scheme and Rec+NT indicates the half conservation reconstruction scheme coupled with the Nessyahu-Tadmor scheme. Theses schemes are compared with the Godunov and Rusanov schemes, and with the simplest version of the Nessyahu-Tadmor scheme [START_REF] Nessyahu | Nonoscillatory central differencing for hyperbolic conservation laws[END_REF] and the simplest version of the MUSCL scheme [vL97] Example 3.1. On Figure 9, we compare numerous schemes on a Riemann problem with ρ L = 5.99924 u L = 19.5975 and p L = 460.894, ρ R = 5.99242 u R = -6.19633 and p R = 46.0950.

The solution consists of three discontinuities moving to the right (cf the book of Toro [START_REF] Eleuterio | Riemann solvers and numerical methods for fluid dynamics[END_REF]). We took a CFL number of 0.4 and discretized the interval in 400 cells. The shocks are very well computed, an we observe an improvement on the contact discontinuity. This improvement is often much better when the reconstruction scheme is coupled with a higher order scheme. Example 3.3. The Figure 12 presents the results at time 1.8 of an entropy satisfying shock interacting with a sine wave, presented for example in [START_REF] Liska | Comparison of several difference schemes on 1D and 2D test problems for the Euler equations[END_REF]. The mesh has 400 cells, the CFL number is 0.45 and the initial datum is:

     ρ 0 (x) = 3.897143 * 1 x<-4 + (1 + 0.2 sin(5x))1 x≥-4 ; u 0 (x) = 2.629369 * 1 x<-4 ; p 0 (x) = 10.33333 * 1 x<-4 + 1 x≥-4 .
The reference solution is, once again, the result given by the Nessyahu-Tadmor scheme with 30 000 cells and a CFL number of 0.48. This test case is a slight perturbation of the pure slowly moving shock introduced by [START_REF] Quirk | A contribution to the great Riemann solver debate[END_REF], where ρ R = p R = 1. Once again, no spurious oscillation appears. Example 3.5. Another nice feature of the reconstruction scheme is that is seems to drastically diminish the wall heating phenomenon. It occurs when a shock reflects on a solid wall, and takes the form of a hollow in a density, or a spike in the temperature, near the wall. We tested our scheme on two cases considered by Donat and Marquinat in [START_REF] Donat | Capturing shock reflections: an improved flux formula[END_REF], who proposed a cure by interlacing two schemes. For those two tests only, γ = 5/3. The first case is a Riemann problem developing two symmetric shocks. The initial datum is

ρ L = 1, u L = 4, p L = 1 and ρ R = 1, u R = -4, p R = 1 .
The simulation is running with 200 cells and a CFL number of 0.4. The results, shown on Figure 14, show that the wall heating phenomenon is drastically diminish with the reconstruction scheme. The second test is the reflection of a gas of density 1, pressure 0.001 and velocity 1 on a solid wall on its right. On Figure 15 is a zoom around the wall, and we can clearly see the wall heating phenomenon and the resulting spurious oscillations for the Godunov's and MUSCL's schemes, and the good behavior of the reconstruction scheme. We took a CFL number of 0.45 and 1 000 cells.

Conclusions and perspectives

In this paper, we presented a finite volume scheme in which the mean value inside each cell is replaced, whenever it is possible, by a single shock. The fluxes are computed by letting these discontinuities evolved during the time step. For convex scalar conservation laws and for the barotropic Euler equation, we proved that the scheme is exact on pure shocks. In the latest case, we use a robust criterion to detect shocks and decide which wave should be reconstruct. The lack of such a criterion for the full gas dynamics explains why we added extra conditions in that case. However, numerical results are encouraging, especially on the problematic test cases of slowly moving shocks and shock reflections. Indeed, the spurious oscillations in the momentum and the hollow in the density are eliminated. Our main perspectives are first to find a more suitable criterion on the full gas dynamics, and second to use an approximate Riemann solver.

Figure 1 :

 1 Figure 1: Conservative reconstruction of a shock.

Figure 6 :

 6 Figure 6: Comparison of the half conservative (HCRec) and fully conservative (FCRec) reconstruction schemes with classical schemes on a unique slowly moving shock (top) and on a Riemann problem consisting of a slowly moving 1-shock followed by a 2-shock (bottom).

Figure 7 :

 7 Figure 7: Left: order of the half conservative reconstruction on the Riemann problem only containing shocks. Right: L 1 error through time for different CFL conditions, for a mesh of 500 cells.

Figure 8 :

 8 Figure8: Comparison of the reconstruction scheme (Rec), the Nessyahu-Tadmor scheme and the coupled scheme (Rec+NT) on a Riemann problem consisting of a 1-shock and a 2-rarefaction wave. The mesh has 100 cells and the CFL number is 0.1.

Figure 9 :

 9 Figure 9: Density at time 0.035 for the first test case.

Figure 11 :

 11 Figure 11: Density at time 0.038 in the blast wave test case, with a mesh of 400 cells (top), 1 200 cells (middle) and 2 000 cells (bottom).

Figure 12 :

 12 Figure 12: Density at time 1.8 for the shock entropy wave interaction

Figure 13 :

 13 Figure 13: Momentum at time 0.3 for a Riemann problem with a slowly moving shock.

Figure 14 :Figure 15 :

 1415 Figure 14: Density and internal energy at time 0.1 for a Riemann problem with two symmetric shocks.

  Example 3.2. We compare those same schemes in the Colella and Woodward blast wave

			4 5 6 7 6		Rusanov Rec Reference Godunov Rusanov Rec Godunov Reference									6 4 5 7 6		Rec+NT Rec+NT NessyahuTadmor NessyahuTadmor Reference MUSCL MUSCL Reference
		density	3 5												density	3 5
			2 4													2 4
		density	1												density	1
			3													3
			0	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1		0	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
			2						X							2						X
		14 1												14 1
		12												12
		10												10
			8 0	0	0.1	0.2	0.3	0.4	0.5 X	0.6	0.7	0.8	0.9	1		0 8 0	0.1	0.2	0.3	0.4	0.5 X	0.6	0.7	0.8	0.9	1
	velocity	2 4 6 7											velocity	4 6 2 7
		-6 -4 -2 0 6		Rusanov Rec Reference Godunov Rusanov Rec Godunov Reference								-6 -4 -2 0 6		Rec+NT Rec+NT NT NessyahuTadmor MUSCL Reference MUSCL Reference
			5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1		5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
									X													X
		1200 4 1400 density								Rusanov Rec		1200 1400 4 density							Rec+NT NessyahuTadmor
			3								Godunov			3							MUSCL
		1000								Reference		1000							Reference
	internal energy	600 800 1 2											internal energy	600 800 1 2
		400												400
		200 0.5 0	0.55	0.6		0.65	0.7	0.75	0.8		0.85	0.9	200 0.5 0	0.55	0.6		0.65	0.7	0.75	0.8	0.85	0.9
									X													X
			0													0
				0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1			0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
									X													X
			7													7
	Rusanov Figure 10: From top to bottom, density, velocity and internal energy at time 0.026 in the Rec+NT Rec NT
	6 blast wave test case, with a mesh of 400 cells. Godunov Reference		6		MUSCL Reference
			5													5
			4													4
		density													density
			3													3
			2													2
			1													1
			0.5 0	0.55	0.6		0.65	0.7	0.75	0.8		0.85	0.9		0.5 0	0.55	0.6		0.65	0.7	0.75	0.8	0.85	0.9
									X													X
														22

Slowly moving shocks

Many finite volume schemes fails to approach correctly shocks that move slowly compared to the speed of sound. Typically, a spike in the momentum appears in the first iterations in time, and is then, by conservation of the momentum, counterbalance by a hollow. The spike oscillates through time in an almost periodic manner, where the period corresponds to the time that the shock needs to cross an entire cell. Therefore, even though the oscillations are diffused by the scheme, it is a constant source of error that is blamed for slow convergence to the steady state (see [START_REF] Noh | Errors for calculations of strong shocks using an artificial viscosity and an artificial heat flux[END_REF] and [START_REF] Menikoff | Errors when shock waves interact due to numerical shock width[END_REF]). Moreover, higher order schemes tend to better preserve the spurious oscillations (that should not be here) better than first order scheme. This problem was first report by Colella and Woodward in [START_REF] Woodward | The numerical simulation of twodimensional fluid flow with strong shocks[END_REF]. In [START_REF] Roberts | The behavior of flux difference splitting schemes near slowly moving shock waves[END_REF] and [START_REF] Arora | On postshock oscillations due to shock capturing schemes in unsteady flows[END_REF], numerous numerical tests and comparison between schemes are performed. Unlike all schemes tested in those papers, the half and fully conservative reconstruction schemes are exact not only on steady shocks, but on all shocks (see Proposition 2.4). This can be checked on the top of Figure 6, where the spurious oscillations created by all the other schemes can be seen. This good property is inherited when computing a Riemann problem containing a slowly moving shock with the half conservative reconstruction scheme, but lost with the fully conservative reconstruction scheme, as depicted on the bottom of Figure 6. Those simulations are run on a mesh with 200 cells, with a CFL number of 0.45. In both simulations, ρ L = 1 and ρ R = 20 and the speed of sound is c = 0.5.

we set the desired reconstructed states to be Ūj,L = U j-1 and Ūj,R = U * ,L j .

In other words, we will try to reconstruct a 1-shock.

we set the desired reconstructed states to be

In other words, we will try to reconstruct a 3-shock.

• Otherwise, and if

we set Ūj,L = U * ,L j and Ūj,R = U * ,R j . and try to reconstruct as a 2-contact discontinuity.

Then, we compute the distances d ρ j , d q j and d E j by conservation of mass, momentum and total energy inside the j-th cell. The reconstruction is accepted when the following conditions are fulfilled:

• Both d ρ j and d E j are between 0 and ∆x;

• Both triplets (ρ n j-1 , ρ n j , ρ n j+1 ) and (u n j-1 , 1 ∆x

x n j-1/2 u rec (x)dx, u n j+1 ) are monotonous;

• Eventually, e rec remains positive on the cell.

Here, u rec and e rec are the piecewise constant reconstructed velocities and internal energy on the j-cell. The reconstruction scheme appears to be much more complicated in this case. This is due to the lack of criterion to detect a dominant contact discontinuity. Indeed, Lemma 2.1 (even with an additional test on the pressure) is still valid for this system as p and u remains constant through contact discontinuities. But it is not entirely satisfactory to decide which wave prevails. Indeed, plenty Riemann problems with a 1-shock and a 2-contact discontinuity verify the inequalities

even though the shock is small and the rarefaction is strong. This is why we added additional constraints like

ρ 0 (x) = 1; u 0 (x) = 0; p 0 (x) = 1 000 * 1 x<0.1 + 0.01 * 1 0.1<x<0.9 + 100 * 1 x≥0.9 .

The solution is computed on the interval [0, 1], with reflective boundary conditions at the two extremities of the interval. The reference solution is obtained by running the Nessyahu Tadmor scheme on 30 000 cells. A more accurate reference solution can be found in [START_REF] Woodward | The numerical simulation of twodimensional fluid flow with strong shocks[END_REF].

On Figures 10, we used a mesh containing 400 cells, a CFL number set to 0.45 and the final time is T = 0.026. We plot the density, velocity and internal energy. The non-dissipative feature of the reconstruction scheme is particularly obvious on this last plot. On Figure 11, the finite time is T = 0.038. We took a CFL number of 0.48 and 400, 1 200 and 2 000 cells.

Observe that even with few points, the discontinuities are sharply captured, even though no reconstruction is performed in the middle area, where we recover the behavior of Lax Friedriechs' or Nessayahu-Tadmor's schemes.