C Chalons

R Duvigneau

C Fiorini

Sensitivity analysis and numerical diffusion effects for hyperbolic PDE systems with discontinuous solutions. The case of barotropic Euler equations in Lagrangian coordinates

Sensitivity analysis (SA) is the study of how the output of a mathematical model is affected by changes in the inputs. SA is widely studied, due to its many applications: uncertainty quantification, quick evaluation of close solutions, and optimization, to name but a few. In this work we show that the classic SA techniques, in particular the continuous sensitivity equation (CSE) method, cannot be used if the mathematical model is a system of hyperbolic partial differential equations (PDEs) with discontinuous solutions. The problem arises from the fact that the CSE method requires the differentiation of the state variable: if the latter is discontinuous, this in turn generates Dirac delta functions in the sensitivity. The focus of the first part of this work is to define a system of sensitivity equations valid also in case of discontinuous state: in order to do that, we add a correction term based on the Rankine-Hugoniot conditions. In the second part of this work we illustrate with some numerical tests how some classical finite volume schemes (an exact Godunov method and a Roe-type method) do not converge to the analytical solution for the sensitivity, due to the numerical diffusion: for this reason, we present an anti-diffusive numerical scheme, which provides the correct results for the sensitivity. In this work we carry out the computation in detail for the barotropic Euler equations in Lagrangian coordinates, but the approach is general and can be applied to any hyperbolic system with discontinuous solutions.

Introduction

Sensitivity analysis (SA) concerns the quantification of changes in Partial Differential Equations (PDEs) solution due to perturbations in the model input. It is obviously a valuable tool for engineering applications, which allows to quantify the physical response of a system to any change of parameter values (geometry, boundary conditions, etc). The underlying concepts have been used for a long time in optimal design methods, which replace now the traditional "trials and errors" approach, to determine rigorously the optimal parameters of a system. Beside optimization, SA methods can also be carried out to measure the performance loss due to an unexpected perturbation of the operational conditions, in the framework of uncertainty quantification approaches. Finally, SA methods can likewise be employed to monitor and explore interactively neighboring solutions for a negligible computational expense.

Two strategies can be implemented to compute such sensitivities. The first and most popular one is the adjoint equation method [START_REF] Jameson | Aerodynamic design via control theory[END_REF][START_REF] Mohammadi | Applied Optimal Shape Design for Fluids[END_REF][START_REF] Pironneau | On optimum design in fluid mechanics[END_REF], which introduces additional adjoint variables to compute the derivative of any functional output with respect to all input parameters. The adjoint equation is independent from the input parameters, yielding this approach very efficient for optimization problems involving a large number of design parameters. However, the adjoint equation should be solved backwards in time, which could lead to practical difficulties for unsteady problems. The second strategy is the sensitivity equation method [START_REF] Borggaard | A PDE sensitivity equation method for optimal aerodynamic design[END_REF][START_REF] Duvigneau | An improved continuous sensitivity equation method for optimal shape design in mixed convection[END_REF][START_REF] Duvigneau | A sensitivity equation method for fast evaluation of nearby flows and uncertainty analysis for shape parameters[END_REF][START_REF] Hristova | A continuous sensitivity equation method for time-dependent incompressible laminar flows[END_REF], which allows to compute the derivative of the PDE solution itself, at any location and time, with respect to a single input parameter. This approach is therefore better suited to the exploration of neighboring solutions. Note that it relies on a forward time integration and, thus, seems to be well adapted to time-dependent hyperbolic PDE systems, which constitute the context of the present work.

However, a major difficulty arises when one considers the application of the sensitivity equation method to hyperbolic PDE systems: in such a context, discontinuities in the solution can appear, even if the initial solution is regular. This corresponds to the generation of shock waves in compressible flows for instance. Such discontinuities in the solution lead to specific issues regarding sensitivity analysis, because they correspond to the presence of Dirac functions in the sensitivity fields. This question has been explored in [START_REF] Bardos | A formalism for the differentiation of conservation laws[END_REF][START_REF] Mohammadi | Applied Optimal Shape Design for Fluids[END_REF] with a theoretical viewpoint, and more recently in [START_REF] Guinot | Upwind finite volume solution of sensitivity equations for hyperbolic systems of conservation laws with discontinuous solutions[END_REF][START_REF] Guinot | An approximate riemann solver for sensitivity equations with discontinuous solutions[END_REF] with a numerical viewpoint. Indeed, the capture of the Dirac peak in the sensitivity solution by numerical schemes is intractable in practice. Therefore, a modification of the sensitivity system was proposed in the later references, to remove the Dirac functions from the numerical sensitivity solution. More specifically, a modification of HLL Riemann solver used to evaluate fluxes in a finite-volume method was proposed, in the context of Saint-Venant equations.

In this work, we adopt a similar point of view: we propose to define and approximate numerically a modified system of sensitivity equations which is valid also when the state is discontinuous. To correct some shortcomings appearing in the sensitivity solution, reported in [START_REF] Guinot | Upwind finite volume solution of sensitivity equations for hyperbolic systems of conservation laws with discontinuous solutions[END_REF][START_REF] Guinot | An approximate riemann solver for sensitivity equations with discontinuous solutions[END_REF], some alternative formulations to define the sensitivity solution update are studied, ranging from Godunov method, first-and second-order Roe-type solvers, to an "anti-diffusive" scheme. These methods are tested numerically in the context of the barotropic Euler equations in Lagrangian coordinates, but the proposed formulations are general and could be used for other hyperbolic systems of conservation laws.

This article is organized as follows: in the first sections, we present the state equations, derive the sensitivity equations and the modification of the sensitivity equations to account for the Dirac functions. Then, we detail the exact resolution of the Riemann problem for the state and sensitivities. In section 5 the implementation of two classical methods (Godunov and Roe-type solvers) is examined and some numerical tests are conducted, which exhibit grid-convergence issues. For this reason, an "anti-diffusive" scheme is introduced in section 6 to alleviate these difficulties and the accuracy of the proposed approach is demonstrated for some cases. Finally, we present an application to uncertainty quantification and we compare the SA results with the Monte Carlo method.

Problem description

As already mentioned in the previous section, standard SA methods can be used only if the solution U is regular enough [START_REF] Bardos | A formalism for the differentiation of conservation laws[END_REF], which is usually not the case for hyperbolic systems of the general form

∂ t U + ∂ x F(U) = 0, x ∈ R, t > 0, U(x, 0) = U 0 (x).
In fact, it is well known this kind of systems can have discontinuous solutions, regardless of the regularity of the initial condition U 0 . If the state U is discontinuous, the sensitivity U a = ∂ a U will exhibit Dirac delta functions. Here and throughout this work, a denotes the parameter of interest which may vary and induce a non trivial sensitivity. We remark that we consider a scalar parameter a only to simplify the notation: everything we do in this work can be done for a vector of parameters a as well, with no added complexity from a theoretical point of view.

In this work, we consider the barotropic Euler equations in Lagrangian coordinates, i.e. the psystem; however, everything can be extended to any hyperbolic system. The choice of the p-system is motivated by the fact that, although quite simple, it presents all the main features of hyperbolic systems: this allows us to solve the state problem easily and to focus on the sensitivity problem.

The system writes:

∂ t τ -∂ x u = 0, ∂ t u + ∂ x p(τ ) = 0, (1)
where τ > 0 is the co-volume (i.e. τ = 1 ρ , and ρ is the density of the fluid), u is the Lagrangian velocity and the pressure p(τ ) is a function only of τ . We assume p (τ ) < 0 and p (τ ) > 0. The Jacobian matrix of the system is the following:

M(τ, u) = 0 -1 p (τ ) 0 (2)
and its eigenvalues are real and distinct λ ± = ±c, where c = -p (τ ) is the Lagrangian sound speed. Therefore M is R-diagonalisable, and ( 1) is strictly hyperbolic. In this work we will consider p(τ ) = τ -γ , where γ = 1.4 is the heat capacity ratio.

If we consider smooth solutions of (1), we can apply the CSE method, differentiate (1) with respect to a and obtain the following sensitivity equations:

∂ t τ a -∂ x u a = 0, ∂ t u a + ∂ x (p (τ )τ a ) = 0. (3)
In order to introduce a more compact notation, we define the state and sensitivity vectors and their fluxes:

U = τ u , F(U) = -u p(τ ) , U a = τ a u a , F a (U, U a ) = -u a p (τ )τ a ,
and rewrite the systems (1) and (3) in a vectorial form:

∂ t U + ∂ x F(U) = 0, ∂ t U a + ∂ x F a (U, U a ) = 0. (4)
The Jacobian matrix of the global system (4) writes:

A(V) = ∂G ∂V =     0 -1 0 0 p (τ ) 0 0 0 0 0 0 -1 p (τ )τ a 0 p (τ ) 0     , with V = U U a , G = F F a .
One can remark that A(V) is a lower triangular block matrix whose diagonal blocks are identical to each other and to the state system's Jacobian matrix. We observe that the global system (4) has two repeated eigenvalues, the same λ ± as the original system (1), and that the matrix A(V) is not R-diagonalisable as soon as τ a = 0. This means that the global system (4) is only weakly hyperbolic. Therefore, the system (4) as it is will provide us, in case of discontinuous state U, with a sensitivity U a presenting Dirac delta functions, in addition to the usual discontinuity, so that these solutions have to be interpreted in the sense of measures. We refer for instance the reader to the following papers, and the references therein: [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF][START_REF] Chalons | Beyond pressureless gas dynamics: quadrature-based moment models[END_REF][START_REF] Forestier | Multivalued solutions to some non-linear and non-strictly hyperbolic systems[END_REF][START_REF] Joseph | A riemann problem whose viscosity solutions contain δ-measures[END_REF][START_REF] Lefloch | An existence and uniqueness result for two nonstrictly hyperbolic systems, ima volumes in math. and its appl. 27[END_REF][START_REF] Yang | New developments of delta shock waves and its applications in systems of conservation laws[END_REF]. However, sensitivities with Dirac delta functions are unusable for many applications: since a Dirac delta function cannot be seized numerically, the spike in the numerical solution mimicking it is spread over multiple cells, leading to a corrupt solution in the neighbourhood of the shock. For this reason, we add a correction term to the sensitivity equations, as done in [START_REF] Guinot | Upwind finite volume solution of sensitivity equations for hyperbolic systems of conservation laws with discontinuous solutions[END_REF]. The definition of a proper correction term is the subject of the next section.

x 1

x 2 x k,s (t 1 )

t 1 t 2 σ k

Source term

In this section, we aim at proposing a new version of (4) which is also valid for discontinuous solutions of the state variable U. Recall indeed that (4) has been derived assuming formally that the solution is smooth whilst hyperbolic equations are well known to develop discontinuities in finite time even for smooth initial data U(x, t = 0) = U 0 (x). In order to compensate the Dirac delta functions that appear in the solutions U a of (4) when U is discontinuous, we add to (3) a source term S of the following form:

S = Ns k=1 δ k ρ ρ ρ k , (5)
where N s is the number of discontinuities in the state solution U, ρ ρ ρ k is the amplitude of the k-th correction (to be computed), and δ k is the Dirac delta function δ k = δ(x -x s,k ), where x s,k is the position of the k-th discontinuity. The new version of (4) we are going to consider thus writes:

∂ t U + ∂ x F(U) = 0, ∂ t U a + ∂ x F a (U, U a ) = S. (6)
Let us motivate our choice and define ρ ρ ρ k by considering a control volume (x 1 , x 2 ) × (t 1 , t 2 ) as in Figure 1, which contains only the k-th discontinuity, propagating at speed σ k . We integrate the second equation of (6) over the control volume:

x2 x1 t2 t1 ∂ t U a + ∂ x F a (U, U a )dxdt = x2 x1 t2 t1 δ(x -x s,k )ρ ρ ρ k (t)dxdt,
and we obtain

x2 x1 U a (x, t 2 ) -U a (x, t 1 )dx + t2 t1 F a (x 2 , t) -F a (x 1 , t)dt = t2 t1 ρ ρ ρ k (t)dt, (7)
where we used the simplified notation F a (x, t) = F a (U(x, t), U a (x, t)). We divide ( 7) by (t 2 -t 1 ) and as the size of the control volume tends to zero we have:

ρ ρ ρ k (t) = (U - a -U + a )σ k + F + a -F - a , (8)
where the plus (respectively minus) stands for the value of the variables to the right (respectively left) of the discontinuity. The relation [START_REF] Duvigneau | A sensitivity equation method for fast evaluation of nearby flows and uncertainty analysis for shape parameters[END_REF] gives a natural meaning of ρ ρ ρ k in terms of a defect measure of the Rankine-Hugoniot conditions for [START_REF] Bouchut | Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources[END_REF]. Now, we want to define ρ ρ ρ k so that the new model with the source term is valid also in the case of discontinuous state. We start from the Rankine-Hugoniot conditions for the state variable U across a discontinuity:

(U --U + )σ k = F --F + ,
and we suggest to differentiate them with respect to the parameter a. As we do that, we should consider the fact that the Rankine-Hugoniot conditions are valid only at the discontinuity location x k,s (t), which depends on the parameter a. Therefore, we obtain

(U - a -U + a )σ k + (U --U + )∂ a σ k + σ k (∇U + -∇U -)∂ a x k,s (t) = = F - a -F + a + ∂F(U + ) ∂U ∇U + - ∂F(U -) ∂U ∇U -∂ a x k,s (t).
The terms depending on ∂ a x k,s (t) are very difficult to estimate. However we remark that, thanks to the presence of the gradients, they are zero if we consider that the solution U is constant in the left and right neighbourhoods of the shock. This is verified in a standard first order finite volume approach. We obtain therefore a simpler formula:

(U - a -U + a )σ k + (U --U + )σ k,a = F - a -F + a , (9)
with σ k,a = ∂ a σ k . Comparing the latter with [START_REF] Duvigneau | A sensitivity equation method for fast evaluation of nearby flows and uncertainty analysis for shape parameters[END_REF], one is thus led to set

ρ ρ ρ k (t) = σ k,a (U + -U -). ( 10
)
Our choice is of course valid for each k-discontinuity of the state solution, leading us to definition [START_REF] Chalons | Sensitivity analysis for the Euler equations in Lagrangian coordinates[END_REF] where the sum is taken over the number of discontinuities. Note also that by construction, if a triple (U -, U + , σ) is associated with an admissible discontinuity with a left (respectively right) state U -(respectively U + ) and σ is the speed of propagation, then the triple

(U - a , U + a , σ) with U - a = ∂ a U -and U + a = ∂ a U + is
also admissible in the sense that it satisfies the generalised Rankine-Hugoniot relations imposed by S. In other words, the sensitivity solution of ( 6) is obtained by differentiating the state solution with respect to a when the solution is smooth or discontinuous with constant left and right states. As far as the initial condition is concerned, we have:

U(x, t = 0) U a (x, t = 0) = U 0 (x) ∂ a U 0 (x) .

Exact solution of the Riemann problem

In this section, we present the exact resolution of the Riemann problem for the state and the sensitivity, associated with the initial data:

U(x, 0) = U L x < x c , U R x > x c , U a (x, 0) = U a,L x < x c , U a,R x > x c ,
for a given x c . First, we compute the solution of the state system (1), which is well-known but necessary to solve (4). Then we derive it with respect to the parameter of interest a to obtain the sensitivity. As we will see, the sensitivity exhibits interesting and non trivial behaviours.

The state variable

We recall that the eigenvalues of the Jacobian matrix of the state system are:

λ 1 (U) = --p (τ ) and λ 2 (U) = -p (τ ) ,
and the eigenvectors:

r 1 (U) =      2 -p (τ ) p (τ ) - 2p (τ ) p (τ )      and r 2 (U) =      - 2 -p (τ ) p (τ ) - 2p (τ ) p (τ )     
, which are chosen in such a way that ∇λ i • r i = 1. Since the couples (λ i , r i ) are both genuinely non linear, the waves associated can be either shocks or rarefaction waves. The structure of the analytical solution of the state is resumed in Figure 2 and it consists of two waves, whose speeds can be computed exactly.

In order to give more details on this structure, which will be necessary to explain the structure of the sensitivity, let us consider the plane (τ, u) and the points U L and U R : starting from U L we need to reach U R passing from an intermediate state U * using shocks and rarefaction waves, see Figure 2 for the notations. First, we compute which points U are reachable through a shock of speed σ from U L . Across a shock, the Rankine-Hugoniot conditions are valid, therefore:

-u + u L = σ(τ -τ L ) p(τ ) -p(τ L ) = σ(u -u L ). ( 11
)
Finding σ from the first equation and replacing it in the second one, one has:

(u -u L ) 2 = -(p(τ ) -p(τ L ))(τ -τ L ),
and we observe that the right-hand side is always positive because p (τ ) < 0, therefore:

u = u L ± -(p(τ ) -p(τ L ))(τ -τ L ) . (12)
In order to chose the sign in [START_REF] Guinot | An approximate riemann solver for sensitivity equations with discontinuous solutions[END_REF], we use the Lax conditions:

λ 1 (U) < σ < λ 1 (U L ),
which implies that σ is negative (because λ 1 is negative), and:

--p (τ ) < --p (τ L ) ⇒ τ < τ L ,
where we used the hypothesis p (τ ) > 0. Both σ and (τ -τ L ) are negative, therefore their product is positive and from the first equation of [START_REF] Guinot | Upwind finite volume solution of sensitivity equations for hyperbolic systems of conservation laws with discontinuous solutions[END_REF] we can say that u L > u. We can conclude that the sign in ( 12) is a minus and therefore the points reachable from U L through a shock are on the curve of equation:

u = u L --(p(τ ) -p(τ L ))(τ -τ L ) . ( 13
)
Repeating everything for the 2-wave one finds that the states U reachable from U R through a shock are those on the curve of equation:

u = u R + -(p(τ ) -p(τ R ))(τ -τ R ) . (14)
We now repeat the same procedure to compute the set of points reachable through a rarefaction wave from U L and from U R as we did for the shock.

The self-similar solution ξ → U 1 (ξ) in a 1-rarefaction wave starting from U L respects the following equation:

U 1 (ξ) = r 1 ( U 1 (ξ)), U 1 (ξ 0 ) = U L , (15)
with ξ = λ 1 (U(ξ)) for all ξ, which gives the following system:

         τ 1 (ξ) = 2 -p (τ 1 ) p (τ 1 ) , τ1 (ξ 0 ) = τ L , ũ 1 (ξ) = - 2p (τ 1 ) p (τ 1 ) , ũ1 (ξ 0 ) = u L . (16)
The first equation can be rewritten as

(p (τ 1 )) = 2 -p (τ 1 ),
and its solution writes p (τ 1 ) = -ξ 2 , therefore:

τ1 = (p ) -1 (-ξ 2 ).
After a change of variable, the second equation writes:

dũ 1 dτ 1 = -p (τ 1 ) ,
hence, one has the curve of points reachable through a 1-rarefaction starting from U L :

ũ1 = u L + τ1 τ L -p (τ ) dτ.
If p(τ ) = τ -γ , one would find:

τ1 = γt 2 (x -x c ) 2 1 γ+1 , ũ1 = u L + 2 √ γ 1 -γ τ 1-γ 2 1 -τ 1-γ 2 L . ( 17
)
Repeating exactly the same, one can find the points reachable through a 2-rarefaction starting from U R :

ũ2 = u R + τ R τ2 -p (τ ) dτ,
hence for our choice of p:

ũ2 = u R + 2 √ γ 1 -γ τ 1-γ 2 R - τ 1-γ 2 2 , τ2 = γt 2 (x -x c ) 2 1 γ+1 . (18)
Finally, in order to compute U * we define the two following 1-and 2-wave curves:

g 1 (τ ; U L ) =    u L --(τ -γ -τ -γ L )(τ -τ L ) if τ ≤ τ L , u L + 2 √ γ 1-γ (τ 1-γ 2 -τ 1-γ 2 L ) if τ > τ L , (19)
g 2 (τ ; U R ) =    u R + -(τ -γ -τ -γ R )(τ -τ R ) if τ ≤ τ R , u R + 2 √ γ 1-γ (τ 1-γ 2 R -τ 1-γ 2 ) if τ > τ R , (20)
which are smooth functions whose derivatives with respect to τ will be denoted g i . The intermediate state τ * is defined as the intersection between g 1 and g 2 , and one has:

u * = g 1 (τ * ; U L ) = g 2 (τ * ; U R ).
Newton's method can be used to compute τ * . We remark that there is no intersection between g 1 and g 2 under the following condition:

u L -u R < 2 √ γ 1 -γ τ 1-γ 2 L + τ 1-γ 2 R .

The sensitivity variable

As already explained, to compute the sensitivity we differentiate with respect to a the state solution, this means that the initial data for the state and the sensitivity are linked by the following relation:

U a,L = ∂U L ∂a U a,R = ∂U R ∂a .
Furthermore, the sensitivity has the same two-wave structure as the state and the waves travel at the same speed as for the state. Therefore, we need to compute the derivative of U * and U and this concludes the computation of the analytical sensitivity. To compute τ * a and u * a , we differentiate, with respect to a the following equality:

g 1 (τ * ; U L ) = g 2 (τ * ; U R )
and we obtain (recall that

U L = (τ L , u L ) T and U R = (τ R , u R ) T ): u * a = g 1 (τ * ; U L )τ * a + ∂g 1 ∂τ L (τ * ; U L )τ a,L + ∂g 1 ∂u L (τ * ; U L )u a,L = = g 2 (τ * ; U R )τ * a + ∂g 2 ∂τ R (τ * ; U R )τ a,R + ∂g 2 ∂u R (τ * ; U R )u a,R ,
which gives an explicit solution (although dependent on U * ) for τ * a :

τ * a = ∂g2 ∂τ R (τ * ; U R )τ a,R + ∂g2 ∂u R (τ * ; U R )u a,R -∂g1 ∂τ L (τ * ; U L )τ a,L -∂g1 ∂u L (τ * ; U L )u a,L g 1 (τ * ; U L ) -g 2 (τ * ; U R )
.

Finally, we differentiate the state solution in the rarefaction U given by ( 17)-( 18).

Examples

In the numerical section of this work, we will consider as parameter of interest the initial data, which means that a can either be τ L , u L , τ R , u R or a combination of them, and from ( 17) and ( 18) one can observe that τi does not depend on those parameters, therefore:

τa,i = ∂ τi ∂a = 0 i = 1, 2.
Concerning the sensitivity of the velocity, one obtains:

ũa,1 = ∂ ũ1 ∂a = u a,L - √ γ τ -1+γ 2 L τ a,L , ũa,2 = ∂ ũ2 ∂a = u a,R + √ γ τ -1+γ 2 R τ a,R .
Interestingly, we remark that the sensitivity is constant in the rarefaction zone of the state variable, which means that for the sensitivity this zone corresponds to at most two discontinuities propagating with velocities given by the extreme left and right velocities of the rarefaction in the state variable, see Figure 3. This simplification is due to the fact that we are considering a reduced Euler system, under barotropic conditions (cf. [START_REF] Guinot | Sensitivity equations for hyperbolic conservation law-based flow models[END_REF]). In particular, there are two cases:

(i) if the state presents a 1-rarefaction (respectively a 2-rarefaction) and the parameter of interest a is τ L (respectively, τ R ), the wave associated with the rarefaction in the sensitivity splits in two discontinuities, as explained above (cf. Figure 4).

(ii) if the parameter of interest is u L (or u R ) we have ũa,1 = u a,L and ũa,2 = u a,R , therefore the wave associated with the rarefaction becomes a single discontinuity for the sensitivity, travelling at the more internal velocity of the state rarefaction wave (cf. Figure 5).

Classical numerical schemes

The aim of this section is to design relevant numerical schemes for [START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF]. As we will see, this task is not easy and requires a nice discretisation of S in order to avoid Dirac delta functions and it is necessary to control numerical diffusion across the shocks where this term is active. Only under these conditions we will get a perfect agreement between exact and numerical solutions. Let us first introduce our notation, although quite classical: we use a constant space step ∆x and a varying time step ∆t n . The mesh interfaces are denoted x j+1/2 = j∆x, the cells C j = [x j-1/2 , x j+1/2 ], the cell centres x j and the intermediate times t n+1 = t n + ∆t n , where ∆t n is chosen according to the usual CFL condition. In the following subsections, we will briefly introduce two classical schemes for the state, Godunov and Roe method, and we will adapt them to the sensitivity.

x t x c U L U * U R (a) 1-shock-2-shock. x t x c U L U * U R U (b) 1-shock-2-rarefaction x t x c U L U * U R U (c) 1-rarefaction-2-shock. x t x c U L U * U R U 1 U 2 (d) 1-rarefaction-2-rarefaction.
Figure 2: Configurations for the state variable U.

κ 1 κ 2 x t x c U a,L U * a U a,R (a) 1-shock-2-shock. x t x c U a,L U * a U a,R U R a κ 1 κ 2 c 2 (b) 1-shock-2-rarefaction. x t x c U a,L U L a U * a U a,R κ 2 κ 1 c 1 (c) 1-rarefaction-2-shock. x t x c U a,L U * a U R a U L a U a,R κ 1 c 1 κ 2 c 2 (d) 1-rarefaction-2-rarefaction.
Figure 3: Corresponding configurations for the sensitivity U a .

x t x c U a,L U * a U a U a,R (a) Corresponding configuration to state case (b)-Figure 2 if a = τ R . x t x c U a,L U * a U a U a,R (b) Corresponding configuration to state case (c)-Figure 2 if a = τ L . x t x c U a,L U a U * a U a,R (c) Corresponding configuration to state case (d)-Figure 2 if a = τ L . x t x c U a,L U * a U a U a,R (d) Corresponding configuration to state case (d)-Figure 2 if a = τ R .
Figure 4: Corresponding configurations for the sensitivity U a , example (i).

x t x c U a,L U * a U a,R (a) Corresponding configuration to state case (a)-Figure 2 if a = u L or a = u R . x t x c U a,L U * a U a,R (b) Corresponding configuration to state case (b)-Figure 2 if a = u L or a = u R . x t x c U a,L U * a U a,R (c) Corresponding configuration to state case (c)-Figure 2 if a = u L or a = u R . x t x c U a,L U * a U a,R (

The Godunov method

In this section, we present an exact Godunov-type method. Being the state equations (1) conservative, the classic update formula can be used:

U n+1 j = U n j - ∆t ∆x (F(U * j+1/2 ) -F(U * j-1/2 )), (21)
where U * j-1/2 is the exact intermediate state known implicitly from ( 19)-( 20), with U L = U j-1 and U R = U j [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF].

The update formula ( 21) cannot be applied to the sensitivity variable, because of the source term. However, as explained in section 4.2, the structure of the sensitivity is made of discontinuities only. Therefore, we can directly compute the average on each cell, if the slopes of the red lines and the solid blue lines in Figure 3 are known at each interface j -1/2. The slopes of the red lines are computed from the Rankine-Hugoniot conditions, while the ones of the blue lines are the eigenvalues evaluated in the correct state. We obtain the following formulas (cf. Figure 3 for the notations):

κ 1,j-1/2 =    u * j-1/2 -uj-1 τj-1-τ * j-1/2 if 1-shock at interface j -1/2, λ 1 (U * j-1/2 ) if 1-rarefaction at interface j -1/2, κ 2,j-1/2 =    u * j-1/2 -uj τj -τ * j-1/2 if 2-shock at interface j -1/2, λ 2 (U * j-1/2 ) if 2-rarefaction at interface j -1/2, c 1,j-1/2 = κ 1,j-1/2 if 1-shock at interface j -1/2, λ 1 (U j-1 ) if 2-rarefaction at interface j -1/2, c 2,j-1/2 = κ 2,j-1/2 if 2-shock at interface j -1/2, λ 2 (U j ) if 2-rarefaction at interface j -1/2.
Then, the update formula for the sensitivity writes:

U n+1 a,j = U n a,j + ∆t ∆x κ 2,j-1/2 (U * a,j-1/2 -U R a,j-1/2 ) + c 2,j-1/2 ( U R a,j-1/2 -U n a,j ) -κ 1,j+1/2 (U * a,j+1/2 -U L a,j+1/2 ) -c 1,j-1/2 ( U L a,j+1/2 -U n a,j ) , (22)
where the intermediate states U * a,j-1/2 , U * a,j+1/2 , U R a,j-1/2 and U L a,j-1/2 are known analytically, from section 4.2. We remark that the source term is encompassed in [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF], since [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF] comes from the exact Riemann solver of (6).

A Roe-type method

First order

In this section we illustrate a Roe-type Riemann solver, consisting of three constant states (which we denote U L , U * and U R for the state, and U a,L , U * a and U a,R for the sensitivity), connected by two discontinuities travelling at velocities

λ ROE L,j-1/2 = -- p(τ n j-1 ) -p(τ n j ) τ n j-1 -τ n j , λ ROE R,j-1/2 = - p(τ n j-1 ) -p(τ n j ) τ n j-1 -τ n j if τ n j-1 = τ n j and ∓ -p (τ n j ) otherwise.
In the following, we will use the notation λ ROE j-1/2 = λ ROE R,j-1/2 = -λ ROE L,j-1/2 . The Harten, Lax and van Leer consistency relations [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF] for the state at the interface j -1/2 write:

U * j-1/2 = 1 2 (U n j-1 + U n j ) - F(U n j ) -F(U n j-1 ) 2λ ROE j-1/2 . ( 23
)
Since U * j-1/2 and λ ROE j-1/2 are known at each interface, we can write the following update formula for the state:

U n+1 j = U n j + ∆t ∆x (λ ROE j-1/2 (U * j-1/2 -U n j ) + λ ROE j+1/2 (U * j+1/2 -U n j )). (24)
Writing the integral conditions for the sensitivity with the source term, one obtains:

F a (U n j , U n a,j ) -F a (U n j-1 , U n a,j-1 ) -∆xS n j-1/2 = λ ROE j-1/2 (U n a,j + U n a,j-1 ) -2λ ROE j-1/2 U * a,j-1/2 ,
from which we have the following form for U * a,j-1/2 :

U * a,j-1/2 = U n a,j-1 + U n a,j 2 - F a (U n j , U n a,j ) -F a (U n j-1 , U n a,j ) 2λ ROE j-1/2 + ∆xS n j-1/2 2λ ROE j-1/2 .
The source term is discretised as follows:

S n j-1/2 = λ ROE a,j-1/2 (U n j-1 -U * j-1/2 ) d 1,j-1 ∆x + (U n j -U * j-1/2 ) d 2,j ∆x , ( 25
)
where d ,j is a shock detector which is equal to 1 if there is an -shock in the j-th cell, it is zero elsewhere and d ,j /∆x approximates numerically the Dirac δ k in the definition of the source term [START_REF] Chalons | Sensitivity analysis for the Euler equations in Lagrangian coordinates[END_REF]. In this work, we use a very simple shock detector: in section 4.1 we showed that the velocity u is decreasing across a shock, whilst the co-volume τ decreases across a 1-shock, and it increases across a 2-shock. Based on this, we set:

d 1,j = 1 if u j > u * j+1/2 and τ j > τ * j+1/2 , 0 otherwise, d 2,j = 1 if u j < u * j-1/2 and τ j > τ * j-1/2 , 0 otherwise.
Finally, U * a,j-1/2 is computed as follows:

U * a,j-1/2 = 1 2 (U n a,j-1 + U n a,j ) - F a (U n j , U n a,j ) -F a (U n j-1 , U n a,j ) 2λ ROE j-1/2 + λ ROE a,j-1/2 2λ ROE j-1/2 (U n j-1 -U * j-1/2 )d 1,j-1 + (U n j -U * j-1/2 )d 2,j . (26)
We remark that the discretisation of the source term ( 25) is such that U * a,j-1/2 = ∂ a U * j-1/2 , in fact differentiating [START_REF] Yang | New developments of delta shock waves and its applications in systems of conservation laws[END_REF] with respect to a, one finds:

∂ a U * j-1/2 = U n a,j-1 + U n a,j 2 - F a (U n j , U n a,j ) -F a (U n j-1 , U n a,j ) 2λ ROE j-1/2 + F(U n j ) -F(U n j-1 ) 2λ ROE j-1/2 λ ROE a,j-1/2 λ ROE j-1/2 .
Using again [START_REF] Yang | New developments of delta shock waves and its applications in systems of conservation laws[END_REF] one has:

U * a = U n a,j-1 + U n a,j 2 - F a (U n j , U n a,j ) -F a (U n j-1 , U n a,j ) 2λ ROE j-1/2 + λ ROE a,j-1/2 2λ ROE j-1/2 (U n j-1 + U n j -2U * j-1/2 ),
which is equal to (26), once the shock detectors are added. Furthermore, the definition (26) encompasses the source term, which means that we can use the update formula (24) for the sensitivity, too. Finally, in order to prepare the second order extension in space, we define a residual as follows:

R I j (U n ) = λ ROE j-1/2 (U * j-1/2 -U n j ) + λ ROE j+1/2 (U * j+1/2 -U n j ), (27)
x x j-3/2 x j-1/2 x j+1/2 x j+3/2 U j-5/4 U j-3/4 U j-1/4 U j+1/4 U j+3/4 U j+5/4
Figure 6: Second order discretisation. In red, the corresponding first order discretisation.

where

R I j (U n ) is a more compact notation for R I (U j-1 , U j , U j+1
). This allows us to write the update formulas in the following way:

U n+1 j = U n j + ∆t ∆x R I j (U n ), U n+1 a,j = U n a,j + ∆t ∆x R I j (U n a ).
Furthermore, it will be useful for the numerical schemes introduced hereafter.

Second order

We extend this scheme to the second order: for the time discretisation we use a two-step Runge-Kutta method, whilst in space we propose a MUSCL-type scheme with some minor modifications in order to have a second order discretisation of the source term. In particular, we remark that ( 10) is valid only if the solution U is locally constant to the left and to the right of the shock, which is true for a first order approximation but not for a second order, in which, classically, the numerical solution is a piecewise affine function. To overcome this problem, we suggest to consider the numerical solution to be a piecewise constant function on half of every cell (cf. [START_REF] Bouchut | Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources[END_REF], section 2.8): the value in the left half (respectively right half) of the j-th cell is denoted U j-1/4 (respectively U j+1/4 ), as shown in Figure 6 and they are computed as in a classical MUSCL approach:

U n j±1/4 = U j ± ∆U n j ,
and a usual choice for ∆U n j is given by a slope-limiter procedure. In this work we use the so-called minmod limiter:

∆U n = 1 2 minmod(U n j+1 -U n j , U n j -U n j-1 ),
where the function minmod is defined as follows:

minmod(a, b) = sgn(a) min(|a|, |b|) if ab > 0, 0 otherwise.
This interpretation of the second order allows us to define the source term as we did for the first order, however we need to consider an additional Riemann problem for each cell. This leads to the following definition of the residual:

R II j (U n ) = (λ ROE j-1/2 (U * j-1/2 -U n j-1/4 )+λ ROE j+1/2 (U * j+1/2 -U n j+1/4 ))+λ ROE j (2U * j -U j-1/4 +U j+1/4 ),
where all the λ ROE and the U * are computed from the extrapolated values U n j±1/4 . Finally, the second order scheme writes:

U n+1/2 j = U n j + ∆t 2∆x R II (U n ), U n+1 j = U n j + ∆t ∆x R II (U n+1/2 ), U n+1/2 a,j = U n a,j + ∆t 2∆x R II (U n a ), U n+1 a,j = U n a,j + ∆t ∆x R II (U n+1/2 a
).

Numerical results

We present some numerical results obtained with the schemes described in the previous section. The spatial domain is (0, 1) and final time is T = 0.03. We consider Riemann problems with x c = 0.5. First, we consider a 1-shock-2-rarefaction case, with the following initial conditions for the state:

U L = 0.7 0 , U R = 0.2 0 .
The parameter of interest is a = u L , so that the initial conditions for the sensitivity are:

U a,L = 0 1 , U a,R = 0 0 .
Figure 7 shows the state variables u and τ and their sensitivities u a and τ a at the final time T . Since the state is a quite classical problem, it is not surprising that all the methods provide very similar solutions one to another. As for the sensitivity, we remark that the modified formulation is able to remove the peak which approximate the Dirac delta function, located at x ≈ 0.4 and evident in the scheme without correction term, whose label is "S = 0" in Figures 78. However, even with the addition of the source term, the sensitivity solution have two issues: first, the discontinuity associated with the state rarefaction is not well captured; secondly, the value of the plateau in the star zone is not the analytical one. Out of these two problems, the first is the less important one, for two reasons: the fact that the state rarefaction splits into two discontinuity for the sensitivity is typical to the PDEs system considered, it does not happen, for instance, in the case of the complete Euler system; furthermore, the numerical solution converges to the analytical one as ∆x goes to 0, meaning that this issue can be solved by using a finer mesh or a higher-order scheme.

The second problem is more critical and we believe that numerical diffusion is the cause of it. In Figure 9 we plot the convergence curves for the all the schemes and for each variable: as one can see, all the methods converge as expected for the state variable; however, for the sensitivity the error seems to be convergent only for coarser meshes, and it reaches a plateau for finer ones. This can be explained if we split the error into two parts: the part concentrated in the rarefaction zone, which is the bigger one in the coarse meshes, converges; however when this part reaches the same order of magnitude as the error in the star zone, which is constant, the plateau is reached. The second test case here presented is an isolated 2-shock for the state as well as for the sensitivity. In order to have an isolated shock we choose the following initial data:

U L = 0.2 g 2 (τ L ; U R ) 0.2 -1.56 , U R = 0.5 -3 ,
where g 2 is the 2-wave curve defined in [START_REF] Mohammadi | Applied Optimal Shape Design for Fluids[END_REF]. As parameter of interest a we choose the arc length of the curve g 2 , which yields the following initial data for the sensitivity:

U a,L = 1 g 2 (τ L ; U R ) 1 -9.35 , U a,R = 0 0 .
Figure 10 shows the results for the state and the sensitivity obtained with a mesh ∆x = 10 -3 : one can notice a spurious wave in the state which does not affect the value in the star zone. However, in the sensitivity this spurious wave is amplified; moreover, the value in the star zone is not correct. Considering the fact that the approximate Riemann solver of Roe is exact in the case of an isolated shock (as well as the exact Godunov solver), the error is necessarily introduced in the average step of the numerical methods and therefore it is due to the numerical diffusion which comes along with the average operation. For this reason, in the next section we introduce a scheme without numerical diffusion in the shock. x j-1/2

U j-1 U * j-1/2 U j U 1 x j+1/2 U * j+1/2 U j+1 U 1 U 2
Figure 11: Definition of the temporary staggered mesh.

the sensitivity: the plateau in the star zone is correct. Furthermore, we show in Figure 14 the convergence results of the classical Roe-type schemes with diffusion compared to the same schemes without diffusion: the latter show a good convergence rate even for the sensitivity variables. We now present another test case with initial data:

U L = 0.7 0 , U R = 0.2 0 , U a,L = 0 0 , U a,R = 1 0 ,
therefore the parameter of interest a is in this case τ R . The initial data for the state is the same as in the previous test case, meaning that we are in configuration (b) of Figure 2 and, since a = τ R , the rarefaction wave splits into two discontinuities for the sensitivity as shown in Figure 5-(a). For this test case we chose a bigger final time (T = 0.07) so that the two extremes of the rarefaction wave could be well separated, in order to attenuate the effect of the numerical diffusion in the middle. We also changed the starting point of the discontinuity (x c = 0.3) in order to have the second discontinuity associated with the rarefaction still in the domain at the final time. The results shown in Figure 15 are obtained with a mesh ∆x = 10 -4 : even in this particular case, with three discontinuities, we are able to approximate well the sensitivity provided that the mesh is fine enough.

Uncertainty quantification

In this section, we show how SA can be used for uncertainty quantification (UQ). The main aim of UQ is to determine a confidence interval for the output of a model, in our case U, given the uncertainty on the input parameters. In this work, we compare two different UQ methods: Monte Carlo and sensitivity analysis. Both methods aim at providing statistical quantities like moments (mean, variance, ...) of the output of the model. In the following, X will stand for one of the variables, considered as random variables, i.e. X can either be τ or u, and X a the corresponding sensitivity. We use the notation µ X to indicate the expected value of the variable X and σ 2 X for its variance. Once this two quantities are known, one can build a confidence interval for the variable X as: CI X = [µ X -κσ X , µ X + κσ X ]. The coefficient κ regulates the amplitude of the interval and it is related to the probability for the variable X to actually fall in the interval. For instance, the choice κ = 2 provides a 95% confidence interval.

Monte Carlo method. Here we briefly introduce the Monte Carlo method. The Monte Carlo method is a probabilistic technique: to obtain an estimate of the average and of the standard deviation one needs to perform multiple simulations. Let a be the vector of uncertain parameters, with a known distribution. Then, N random samples a i are drawn from this distribution, and for each a i the corresponding solution X i is computed. Then, the unbiased average and variance estimators are used:

µ X = 1 N N i=1 X i , σ 2 X = 1 N -1 N i=1 (µ X -X i ) 2 .
These estimates are satisfactory if N is sufficiently large: the slow convergence, and therefore the high computational cost, is probably the main limitation of the Monte Carlo method. Sensitivity analysis method. SA is a deterministic approach to estimate the average µ X and the variance σ 2 X of the output X. Let µ a be the average of the uncertain vector a, and σ a the covariance matrix:

µ a =    µ a1 . . . µ a M    , σ a =      σ 2 a1 cov(a 1 , a 2 ) . . . cov(a 1 , a M ) cov(a 1 , a 2 ) σ 2 a2 . . . cov(a 2 , a M ) . . . . . . . . . cov(a 1 , a M ) . . . σ 2 a M     
, where M is the number of uncertain parameters, µ ai the average of the i-th parameter, σ 2 ai its variance and cov(•, •) the covariance. One can obtain the following first order estimates of the average and the variance of the variable X:

µ X = X(µ a ), σ 2 X = M i=1 X 2 ai σ 2 ai + M i,j=1 i =j
X ai X aj cov(a i , a j ).

Higher order estimates require higher order sensitivities [START_REF] Martinelli | On the use of second-order derivative and metamodel-based monte-carlo for uncertainty estimation in aerodynamics[END_REF].

Test cases. We applied these techniques to two different test cases. First, we considered a Riemann problem as the one shown in Figures 1213, with uncertain left and right values, i.e. the vector of uncertain parameters is a = (τ R , τ L , u R , u L ) t , with average µ a = (0.2, 0.7, 0, 0) t . We consider that all the parameters are independent of each other, i.e. the covariance matrix is diagonal, and that each parameter has a variance of 0.1% of its average, except for the velocity, whose variance is chosen as 0.0001 for both the left and right values.

In Figure 16 we show the results of the Monte Carlo approach: the average and the average plus and minus twice the standard deviation (i.e. κ = 2) are plotted in red, five samples are plotted in black. These results are obtained with N = 500 samples, on a mesh with ∆x = 2 × 10 -3 using a Roe second order anti-diffusive scheme. As one can see, the average process smudges the shock and the standard deviation is bigger in that zone. In Figure 17 we show the results of the sensitivity approach, on the same mesh and with the same numerical scheme, when the sensitivity is computed with the correction term.

The second test case is a problem with the following initial data for the state:

τ (x, 0) =          τ L x < x c -2 , (m -τ L ) sin 2 x-xc π + pi 2 + τ L x c -2 ≤ x ≤ x c , (m -τ R ) sin 2 x-xc π + pi 2 + τ R x c ≤ x ≤ x c + 2 , τ R x ≥ x c + 2 , u(x, 0) = 0.
The initial data for τ is plotted in Figure 18. The vector of uncertain parameters is a = (τ R , τ L , x c , m, ) t , with average µ a = (0.4, 0.7, 0.5, 0.05, 0.25) t . Each parameter has a variance of 0.1% and they are all independent of each other, leading to a diagonal covariance matrix σ a .

In Figure 19 we show the results of the Monte Carlo approach: the average and the average plus and minus twice the standard deviation (i.e. κ = 2) are plotted in red, five samples are plotted in black. These results are obtained with N = 500 samples, on a mesh with ∆x = 2 × 10 -3 using a Roe second order anti-diffusive scheme. As one can see, the average process smudges the shock and the standard deviation is bigger in that zone. In Figure 20 we show the results of the sensitivity approach, on the same mesh and with the same numerical scheme, when the sensitivity is computed with the correction term.

One can notice that in the regular zones the two methods provide exactly the same results for the Riemann problem, and similar results for the second test case considered, for which the intervals computed with the sensitivity method are slightly larger in the middle zone certainly due to nonlinear effects which are not taken into account by this first-order estimate of the moments. In the discontinuity, the change in the shock speed is neglected with the sensitivity approach, because the average is approximated with only one realization of the state, and this is why the samples do not fall in the confidence interval in that zone for both test cases. However, the sensitivity approach is less expensive: the Monte Carlo approach requires 500 solutions of the state, whilst the SA approach necessitates only one solution of the state and as many solutions of the sensitivity as the number of uncertain parameters, i.e. 4 for the first test case and 5 for the second. Figure 20: SA approach with correction, second order anti-diffusive scheme. Average and the average plus and minus twice the standard deviation in red. Five samples in black dashed lines

Conclusion and discussion

The first goal of this work was to define a sensitivity system providing a solution which is suitable for applications, i.e. without peaks approximating the Dirac delta function: this was achieved by adding a properly defined source term. However, the numerical results presented in section 5.3 show that the numerical diffusion plays a very important role in the discretisation of the sensitivity system, to such an extent that classical finite volume schemes do not converge to the analytical solution: in particular, the value of the plateau in the star zone is not correct. To overcome this problem, we propose a numerical scheme based on sampling techniques, which does not have any numerical diffusion: with this scheme we are able to discretise more precisely the source term and to obtain a correct solution for the sensitivity, too. Currently, we are extending this to the complete Euler system. This will allow us to tackle more realistic problems and to use the sensitivity analysis for the applications described in the introduction.

Figure 1 :

 1 Figure 1: Control volume in light blue, x k,s (t) is the position of the k-th discontinuity at time t.

  d) Corresponding configuration to state case (d)-Figure 2 if a = u L or a = u R .

Figure 5 :

 5 Figure 5: Corresponding configurations for the sensitivity U a , example (ii).

Figure 7 :Figure 8 :

 78 Figure 7: Classical finite volume schemes.

Figure 16 :Figure 17 :Figure 18 :Figure 19 :

 16171819 Figure16: Monte Carlo approach for the Riemann problem, second order anti-diffusive scheme. Average and the average plus and minus twice the standard deviation in red. Five samples in black dashed lines

6 An anti-diffusive Roe-type numerical scheme

Since we believe that the failure of the previous schemes is caused by the numerical diffusion in the shock, we present a scheme which does not have any numerical diffusion in the shock. The scheme was first introduced in [START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF] and here we adapt it to the sensitivity problem. It is a modified Godunov method and it can be coupled with any Riemann solver, in this work we couple it with the Roe-type method proposed in the previous section. In fact, the first step is to solve a Riemann problem at each interface, as for a standard Godunov method. The difference between the two methods is in the average step: instead of averaging on the cells [x j-1/2 , x j+1/2 ], a new temporary mesh is defined, whose j-th cell is denoted [x n j-1/2 , x n j+1/2 ], and the average is performed on this mesh. The new mesh is non uniform and it is defined as follows:

where σ n j-1/2 is a proper speed and it depends on the problem. The average operation on the modified mesh provides us with a piecewise constant solution on the new mesh, which we denote U n+1 j

. The final step of this method is to go back to the initial mesh, i.e. compute U n+1 j starting from U n+1 j , and this is done using a sampling technique: the value of the solution on the j-th cell

, and U n+1 j+1 , in agreement with their rate of presence in the cell. More precisely, given a random sequence (α n ) varying in (0, 1), the choice is the following:

The sampling technique mimics the classical averaging if (α n ) is a well distributed random sequence, for instance α n ∼ U(0, 1), or if it is a deterministic low discrepancy sequence, such as the van der Corput sequence (cf. [START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF]):

where i k = 0, 1 is the binary expansion of the integers. Our choice for σ n j+1/2 is the following:

If u is increasing, which means that a rarefaction is expected, the mesh is not modified, whilst in case of expected shock the mesh follows it: in this way one never performs the average across a shock and therefore there is no numerical diffusion (cf. Figure 11).

Remark. Considering only the initial (non moving) mesh, we remark that this method can also be understood as solving the following two-step problem:

The first step is equivalent to solving the Riemann problems at each interface and performing the average on the initial uniform mesh, whilst the second step is equivalent to the sampling (28).

First order formulas

Here, as already said, we couple this anti-diffusive approach with the same Roe-type approximate Riemann solver we presented in the previous section, so we define the following residual:

where in the definition of σ j-1/2 we use λ j-1/2 = λ ROE j-1/2 . Then, the scheme writes:

Second order formulas

One can also couple this approach with the second order Roe-type scheme from the previous section. As for the first order, we define the new residual as follows:

Then the scheme writes:

From the two-step problem (30) point of view, the discretisation (32) is a second order discretisation of the first step followed by the second step, i.e. the sampling technique, which remains unvaried.

Numerical results of the anti-diffusive method

The results of the anti-diffusive method for the test case are shown in Figures