M L Rosinberg

G Tarjus

Comment on "Thermodynamic uncertainty relation for time-delayed Langevin systems"

. Here we show that the derivation is erroneous.

An important recent development in the field of stochastic thermodynamics has been the discovery of the so-called thermodynamic uncertainty relations (TURs) that provide general lower bounds on the fluctuations of time-integrated currents in nonequilibrium systems (see e.g. [2] and references therein). Such relations have been so far established for Markov processes only but in a recent work [1] Vu and Hasegawa have presented an extension to time-delayed Langevin systems in a steady state. If correct, this would be an interesting result since delays are ubiquitous in real-world processes, for instance in biology. Unfortunately, the arguments in [1] are incorrect, as we show in the present comment.

The steady-state TUR for a general Markovian dynamics is expressed as

2 (T ) ≡ Θ 2 -Θ 2 Θ 2 ≥ 2 Σ , (1)
where Θ is an arbitrary current integrated over some observation time T and Σ is the total entropy production accumulated by T (in units where Boltzmann's constant is set to k B = 1). According to [1], this relation remains valid for a time-delayed Langevin dynamics provided Σ is replaced by a "generalized" dissipation Σ g (defined by Eq. ( 13) in [1] or Eq. ( 6) below). This is an intriguing result, but we here show that it follows from an incorrect treatment of the non-Markovian nature of the dynamics. Specifically, the original Langevin equation for the N -dimensional random variable x(t) (cf. Eq. ( 2) in [1]),

ẋ = F(x, x τ ) + √ 2Dξ , ( 2
)
where τ is the delay, x τ ≡ x(t -τ ), F(x, x τ ) is a drift force, and ξ is a Gaussian white noise, has been mistakenly replaced by

ẋ = F(x) + √ 2Dξ , (3)
where F(x) is the (instantaneous) effective force defined by F(x)P ss (x) = F(x, x τ )P ss (x, t; x τ , t -τ ) dx τ (here, P ss (x) and P ss (x, t; x τ , t -τ ) are the steady-state one time and two-time probability distributions, respectively). This replacement allows Vu and Hasegawa to * Electronic address: mlr@lptmc.jussieu.fr express the probability density of a stochastic trajectory as

P(Γ) ∝ exp -(1/4D) T 0 ẋt -F(x t ) 2 dt
and to obtain a lower bound on 2 (T ) by repeating the derivation performed in [3] for a Markovian Langevin dynamics. The point we want to stress is that P(Γ) is not the probability of observing a trajectory generated by the non-Markovian dynamics described by Eq. ( 2) in a steady state, despite the fact that Eq. ( 3) leads to the same probability distribution P ss (x) as Eq. ( 2). The same mistake was made in Ref. [4] and signaled in [5] where P(Γ) as given above was shown to differ from the exact path probability computed for a linear time-delayed Langevin equation (in the case T ≤ τ ). In other words, we argue that the inequality derived in [1] applies to an effective stochastic dynamics that is not the true one. For the same reason, and contrary to the claim in Ref. [4], which is repeated in [1], the quantity ∆S tot g (cf. Eq. ( 9) in [1]) does not satisfy an integral fluctuation theorem (IFT) with the actual dynamics described by (2). In fact, as shown in [5], there is another candidate for the entropy production in time-delayed systems, which is obtained from time inversion and satisfies a proper IFT.

To illustrate our point, we explicitly show that 2/ Σ g is not a lower bound on the squared relative uncertainty 2 (T ). To this aim, we consider a two-dimensional version of Eq. (2) with

F(x, x τ ) = -a 11 x 1 -a 12 x 2,τ -a 21 x 1,τ -a 22 x 2 , (4)
and we choose Θ = -

T 0 {[a 11 x 1 (t) + a 12 x 2 (t)] • ẋ1 (t) + [a 12 (t)x 1 (t) + a 22 x 2 (t)] • ẋ2 (t)
}dt as the current, where • denotes the Stratonovich product. The model studied in section IV.C of [1] corresponds to the symmetric case a 11 = a 22 and a 12 = -a 21 . We here focus on the model recently studied in [6] in which there is no feedback from 1 to 2. Specifically, we take a 11 = a, a 22 = b, a 12 = -c, and a 21 = 0. Note that these models are exactly solvable in a steady state due to the linearity of the force F(x, x τ ) and the Gaussian character of the white noise, which makes all probability distributions Gaussian. Therefore, there is no need to restrict the study to the small-τ limit, as done in [1].

In particular, using the same method as [7], one can easily compute the steady-state correlation functions φ ij (t) ≡ x i (0)x j (t) for 0 ≤ t ≤ τ . For instance, we arXiv:1810.12467v1 [cond-mat.stat-mech] 30 Oct 2018

find φ 21 (t) = Dc/[b(a + b)]e b(t-τ ) , from which we get 1 T Θ = c φ21 (0) = D c 2 a + b e -bτ . (5)
The calculation of Σ g , defined in Ref. [1] as

1 T Σ g ≡ 1 D F(x t ) • ẋt , (6)
is also quite easy because the effective force is linear, i.e.,

F 1 = -K 11 x 1 -K 12 x 2 , F 2 = F 2 = -ax 2 ,
P ss (x) ∝ exp[-(1/2)x T σ -1
x] and σ is the covariance matrix with elements σ ij ≡ φ ij (0). K 11 and K 12 are then expressed in terms of the σ ij 's. This eventually yields

K 11 = ab[(a + b) 2 e bτ + c 2 e -bτ ]/[(ab + b 2 + c 2 )(a + b)e bτ -ac 2 e -bτ ], K 12 = -bc[(a + b) 2 + c 2 ]/[(ab + b 2 + c 2 )(a + b)e bτ -ac 2 e -bτ ]
, and in turn

Σ g /T = bc 2 [(a + b) 2 + c 2 ]e -bτ (a + b)[(ab + b 2 + c 2 )(a + b)e bτ -ac 2 e -bτ ] . (7)
In the more general case of the force defined by Eq. ( 4), solving the Fokker-Planck equation does not fully determine F(x), but one can then use the expression of the transition probability of Gaussian stationary processes in terms of the correlation functions (see Eq. (A1) in [6]). (In passing, we also note that F(x) at the order τ is not obtained by simply taking the τ = 0 limit of the transition probability, as defined in Eq. (30) in [1]. For instance, in the model considered in section IV.C of [1], the exact calculation shows that the coefficient of x 1 in F 1 , and of

x 2 in F 2 , is -a + b 2 τ + O(τ 2
). Accordingly, one should have A = a -b 2 τ in the expression (43) of P ss (x), implying that the variance of x 1 and x 2 increases with τ instead of decreasing. This error suggests that the small-τ limit is also incorrect in the two other examples considered in [1]. However, this may be undetectable at the scale of the figures displayed in [1].) Finally, we compute the variance of Θ, and for simplicity we focus on the long-time limit.

Then lim

T →∞ T -1 [ Θ 2 -Θ 2 ] = χ Θ (0)
where χ Θ (k) is the scaled cumulant generating function defined by χ Θ (k) = lim T →∞ T -1 ln e kΘ . A standard calculation using discrete Fourier series (see e.g. [8]) yields

χ Θ (k) = -1/(2π) ∞ 0 dω ln[1 -F k (ω)] with F k (ω) = 4kDc 2 ω[a sin(ωτ )+ω cos(ωτ )+kDω]/[(a 2 +ω 2 )(b 2 +ω 2 )]. This leads to χ Θ (0) = D 2 c 2 b(a + b) 3 (a + b)(2ab + 2b 2 + c 2 ) + c 2 [b(1 + 2bτ ) -a(1 -2bτ )]e -2bτ . (8)
An example of the behavior of the quantity R Θ ≡ lim T →∞ T [ 2 (T ) -2/ Σ g ] as a function of τ is shown in Fig. 1. We observe that R Θ becomes negative for large values of τ , thus invalidating the TUR derived in [1] (more generally, the parabolic lower bound (25) on χ Θ (k) is invalid). On the other hand, as expected, R Θ is always positive if 2 (T ) is calculated with the effective stochastic dynamics defined by Eq. ( 3). We have confirmed these results by performing numerical simulations of the two dynamics. In conclusion, the extension of the TUR to timedelayed Langevin systems is still an open problem. Whether or not the connection between TUR and Fisher information recently discussed in [9][10][11] offers a possible solution remains to be seen.

FIG. 1 :

 1 FIG. 1: (Color on line) RΘ ≡ limT →∞ T [ 2 (T ) -2/ Σg ] as a function of τ calculated with the original time-delayed dynamics [Eq. (2)] (solid black line) and the effective Markovian dynamics [Eq. (3)](dashed red line). The model parameters are a = 1, b = 0.2, c = 0.5.

  and the unknown coefficients K 11 and K 12 can be readily obtained by solving the steady-state Fokker-Planck equation i=1,2 ∂ xi [-F i (x)P ss (x) + D∂ xi P ss (x)] = 0 where


Camilla Fiorini

Christophe Chalons

Régis Duvigneau

A modified sensitivity equation method for the Euler equations in presence of shocks

The Continuous Sensitivity Equation (CSE) method allows to quantify how changes in the input of a Partial Differential Equation (PDE) model affect the outputs, by solving additional PDEs obtained by differentiating the model. However, this method cannot be used directly in the framework of hyperbolic PDE systems with discontinuous solution, because it yields Dirac delta functions in the sensitivity solution at the location of state discontinuities. This difficulty is well known from theoretical viewpoint, but only a few works can be found in the literature regarding the possible numerical treatment. Therefore, we investigate in this study how classical numerical schemes for compressible Euler equations can be modified to account for shocks when computing the sensitivity solution. In particular, we propose the introduction of a source term, that allows to remove the spikes associated to the Dirac delta functions in the numerical solution. Numerical studies exhibit a strong impact of the numerical diffusion on the accuracy of this strategy. Therefore, we propose an anti-diffusive numerical scheme coupled with the approximate Riemann solver of Roe for the state problem. For the sensitivity problem, two different numerical schemes are implemented and compared: one which takes into account the contact wave and another that neglects it. The effects of the numerical diffusion on the convergence of the schemes with respect to the grid are discussed. Finally, an application to uncertainty propagation is investigated and the different numerical schemes are compared.

Introduction

The study of how changes in the inputs of a model affect the outputs is critical for several engineering processes, such as design optimization or uncertainty quantification. This task is usually referred as sensitivity analysis (SA) and can be done in many ways, depending on the nature of the model, the amplitude of the perturbations considered, their deterministic or stochastic nature, etc [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF]. In the present work, we consider only systems governed be Partial Differential Equations (PDEs) and we focus on the estimation of the derivative of the PDE solution with respect to an input parameter. This approach is intrinsically local and only make sense for perturbations of small amplitude, especially for highly non-linear models. The estimation of the derivative is achieved by solving a set of additional PDEs obtained by differentiating the original PDE model with respect to a single input parameter of interest. This approach is referred as the sensitivity equation method for PDE models [START_REF] Borggaard | A PDE sensitivity equation method for optimal aerodynamic design[END_REF], and is closely related to the linear perturbation method [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF].

In the specific case of PDE models, there are two main classes of methods to compute the sensitivities: the discretise-then-differentiate approach and the differentiate-then-discretise one. Both strategies have advantages and disadvantages and both are valid and are suitable for different applications. The differentiate-then-discretise approach is usually considered as more flexible, because it does not require the knowledge of how the original PDE model is solved, and is qualified as non-intrusive. On the contrary, the discretise-then-differentiate approach necessitates the knowledge of the discretized equations, but yields a set of consistent derivatives. A detailed comparison between the two for optimization problems is done in [START_REF] Gunzburger | Perspectives in flow control and optimization[END_REF]. In this work, we focus on the differentiate-then-discretise approach, referred in this context as the Continuous Sensitivity Equation (CSE) 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]. Therefore, the sensitivity equations are obtained by formally differentiating the PDE model with respect to the parameter of interest, and then by exchanging the derivatives with respect to the parameter with the ones in space and time, yielding a new system of PDEs that should be discretized and solved numerically.

However, this method works only under certain assumptions of regularity of the state solution, which may not be verified in the hyperbolic framework. In fact, if this technique is directly applied to hyperbolic equations in case of discontinuous solutions, Dirac delta functions will appear in the sensitivity. 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][START_REF] Chalons | Sensitivity analysis and numerical diffusion effects for hyperbolic PDE systems with discontinuous solutions. The case of barotropic Euler equations in Lagrangian coordinates[END_REF] with a numerical viewpoint. While some authors have adapted their numerical strategies to handle the Dirac delta functions in the solution [START_REF] Clarisse | A Godunov-type method in Lagrangian coordinates for computing linearly-perturbed planar-symmetric flows of gas dynamics[END_REF][START_REF] Després | Robust Uncertainty Propagation in Systems of Conservation Laws with the Entropy Closure Method[END_REF][START_REF] Godlewski | The linearized stability of solutions of nonlinear hyperbolic systems of conservation laws: A general numerical approach[END_REF], others have proposed a modification of the sensitivity system to "remove" the spikes from the numerical sensitivity solution [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], while maintaining the original solution in the regular regions. This is mainly motivated by the observation that the spikes can difficultly be seized numerically, even if they are physical, and do not interact well with classical numerical schemes [START_REF] Appel | Sensitivity calculations for conservation laws with application to discontinuous fluid flows[END_REF][START_REF] Guinot | Sensitivity equations for hyperbolic conservation law-based flow models[END_REF][START_REF] Mohammadi | Applied Optimal Shape Design for Fluids[END_REF]. We already contributed to these investigations, in particular in the context of the barotropic Euler system in Lagrangian coordinates, i.e. the p-system [START_REF] Chalons | Sensitivity analysis and numerical diffusion effects for hyperbolic PDE systems with discontinuous solutions. The case of barotropic Euler equations in Lagrangian coordinates[END_REF]. In this paper, we extend the proposed methodology to the complete compressible Euler system. Firstly, to remove the barotropic condition the additional energy equation must be considered, and this leads to the presence of a third wave, which is a contact discontinuity. Secondly, in Lagrangian coordinates the sign of the speed of the waves is known, which is not the case in Eulerian coordinates. These facts lead to a slightly more complicated design of the numerical schemes. Numerical results show that the numerical diffusion plays an important role in this framework, so particular attention is given to the design of anti-diffusive numerical schemes. Finally, another objective of this paper is to investigate, for a simple problem of uncertainty propagation, the impact of removing the spikes in the sensitivity solution and compare the different schemes designed in this context.

The paper is organised as follows: in the first sections, we introduce the state equations and derive the sensitivity equations. Then, the modification of the sensitivity equations to account for the Dirac delta functions is presented and the new sensitivity system is introduced. Next, we detail the exact resolution of the Riemann problem for the state and sensitivity in a specific case, known as the Sod shock tube problem. Some diffusive and anti-diffusive numerical schemes are illustrated: in particular, for the state a Roe Riemann solver is proposed, and two different schemes are designed for the sensitivity. Some numerical convergence tests are conducted, which exhibit grid-convergence issues of the diffusive schemes and a faster convergence for the anti-diffusive ones. Finally, an uncertainty quantification problem is defined and the results of the diffusive and anti-diffusive, with and without correction term, are compared to the results of the Monte Carlo method.

The state system

The Euler system writes:

   ∂tρ + ∂x(ρu) = 0, ∂t(ρu) + ∂x(ρu 2 + p) = 0, ∂t(ρE) + ∂x(u(ρE + p)) = 0, ( 1
)
where ρ is the density, u is the velocity, ρE the total energy per volume unit, and p the pressure.

The system is closed by the following algebraic equation:

p = (γ -1) ρE - 1 2 ρu 2 , (2)
where γ = 1.4 is the heat capacity ratio. We introduce two other quantities which will be useful in the following: the total enthalpy H = E + p ρ and the speed of sound c = (γ -1)(H -1 2 u 2 ) . We can rewrite the system (1) in the vectorial form:

∂tU + ∂xF(U) = 0, (3)
where

U =   ρ ρu ρE   =   w1 w2 w3   , F(U) =   ρu ρu 2 + p u(ρE + p)   =     w2 w 2 2 w 1 + (γ -1) w3 -1 2 w 2 2 w 1 γ w 2 w 3 w 1 -(γ-1) 2 w 3 2 w 2 1     .
One can also write (1) in the nonconservative form:

∂tU + A(U)∂xU = 0, (4)
where the Jacobian matrix A writes:

A(U) = ∂F ∂U =   0 1 0 γ-3 2 u 2 (3 -γ)u γ -1 γ-2 2 u 3 -c 2 u γ-1 3-2γ 2 u 2 + c 2 γ-1 γu   ,
its eigenvalues are λ1 = u -c, λ2 = u, and λ3 = u + c and its eigenvectors are:

r1 =   1 u -c H -uc   , r2 =   1 u u 2 2   , r3 =   1 u + c H + uc   .
Therefore A is R-diagonalisable and the system (1) is strictly hyperbolic. At last, (3) will be supplemented with a given initial data U(x, t = 0) = U0(x), ∀x ∈ R.

The sensitivity system

Considering only smooth solutions of (1), one can apply the Continuous Sensitivity Equation (CSE) [START_REF] Hristova | A continuous sensitivity equation method for time-dependent incompressible laminar flows[END_REF][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] method which consists in differentiating (1) with respect to the parameter of interest a. One can then formally exchange the derivatives in time and space with the ones with respect to a (see [START_REF] Bardos | A formalism for the differentiation of conservation laws[END_REF] for the theoretical aspects) and obtain the following sensitivity system:

   ∂tρa + ∂x(ρu)a = 0, ∂t(ρu)a + ∂x(ρau 2 + 2ρuua + pa) = 0, ∂t(ρE)a + ∂x(ua(ρE + p) + u((ρE)a + pa)) = 0, (5)
which can be written in vectorial form as ∂tUa + ∂xFa(U, Ua) = 0, [START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF] where we used the following notation:

Ua = ∂aU =   ρa (ρu)a (ρE)a   , Fa(U, Ua) = ∂aF(U) =   (ρu)a ρau 2 + 2ρuua + pa ua(ρE + p) + u((ρE)a + pa)   .
Note that differentiating (2) one has:

pa = (γ -1)((ρE)a - 1 2 ρau 2 -ρauua)
which acts as a closure relation for [START_REF] Chalons | Sensitivity analysis and numerical diffusion effects for hyperbolic PDE systems with discontinuous solutions. The case of barotropic Euler equations in Lagrangian coordinates[END_REF]. The initial data for the sensitivity is Ua(x, t = 0) = ∂aU0(x).

The global system

In order to write the global system, i.e. the state and sensitivity system, in a more compact way, we introduce the following vectors:

V = U Ua =         w1 w2 w3 w4 w5 w6         , G(V) = F(U) Fa(U, Ua) =             w2 w 2 2 w 1 + (γ -1) w3 -1 2 w 2 2 w 1 γ w 2 w 3 w 1 -(γ-1) 2 w 3 2 w 2 1 ) w5 γ-3 2
w 2 2 w 4 w 2 1 -(γ -3) w 2 w 5 w 1 + (γ -1)w6 γ w 3 w 5 w 1 -γ w 2 w 3 w 4 w 2 1 -3 2 (γ -1)
w 2 2 w 5 w 1 + (γ -1)
w 3 2 w 4 w 3 1 + γ w 2 w 6 w 1             .
Therefore, the complete system writes:

∂tV + ∂xG(V) = 0, V(x, 0) = V0(x), (7)
with V0(x) = (U0(x), ∂aU0(x)) t . The Jacobian matrix of the complete system has the following form:

∂G(V) ∂V = M(V) = A 0
B A where A is the Jacobian matrix of the state system and B writes:

B =   0 0 0 (γ -3)uua (3 -γ)ua 0 ( ) (•) γua   , (8)
with

( ) = - c 2 γ -1 pa p u + 3 2 (γ -2)u 2 ua + c 2 γ -1 uρa ρ - c 2 γ -1 ua + γ u 3 ρa ρ , and
(•) = γ 2 u 2 ρa - c 2 γ -1 ρa + 6 -5γ 2
u 2 ρa ρ + (3 -2γ)uua + 3(γ -1) uρa ρ 2 + c 2 γ -1 pa p .
The matrix M has three repeated eigenvalues, which are the eigenvalues of the matrix A. More precisely, one can prove the following result.

Proposition. 1. The global system ( 7) is weakly hyperbolic.

Proof. A system of the form ( 7) is weakly hyperbolic if its Jacobian matrix has real eigenvalues and it is not R-diagonalisable. We want to investigate whether or not the matrix M is R-diagonalisable.

A matrix is diagonalisable if and only if its minimal polynomial splits in distinct roots. Since the characteristic polynomial of the matrix M is the following:

pM (x) = (x -λ1) 2 (x -λ2) 2 (x -λ3) 2 , (9)
the minimal polynomial, in order to have distinct roots, can be at most of degree 3. Therefore, if M is diagonalisable, it must be:

(M -λ1I6)(M -λ2I6)(M -λ3I6) = 0. ( 10
)
Let us write (10) by blocks:

A -λ1I3 0 B A -λ1I3 A -λ2I3 0 B A -λ2I3 A -λ3I3 0 B A -λ3I3 = 0 (11)
Developing the left-hand side products one obtains the following matrix:

(A -λ1I3)(A -λ2I3)(A -λ3I3) 0 ( ) (A -λ1I3)(A -λ2I3)(A -λ3I3) , (12)
where

( ) = B(A -λ2I3)(A -λ3I3) + (A -λ1I3)B(A -λ3I3) + (A -λ1I3)(A -λ2I3)B.
The top-left and bottom-right coefficients are equal to the characteristic polynomial of A evaluated in A, thus they are zero. Therefore, the matrix M is diagonalisable if and only if ( ) = 0. Let us compute the coefficient (1, 1) of ( ):

( ) (1,1) = 0 + [c -u, 1, 0]   0 (3 -γ)u 2 ua -(γ-3) 2 2 u 2 ua ♦   + +[c -u, 1, 0]   (γ -3)uua (γ -2)(3 -γ)u 2 ua + (γ -1)( )   = = - 3 2 (γ -1)u 2 ua + (γ -3)cuua -c 2 u pa p + c 2 u ρa ρ -c 2 ua + γ(γ -1)u 3 ρa ρ ,
where there is no need to specify ♦ and . There is no reason why the quantity should be always be zero. Therefore, the matrix is not diagonalisable and the complete system is not hyperbolic in general. However, as the eigenvalues are real, the system is weakly hyperbolic.

Source term

The sensitivity system (5) was derived assuming that the state solution U is regular. However, this is not generally true for hyperbolic systems such as the one considered [START_REF] Bardos | A formalism for the differentiation of conservation laws[END_REF]: if the state is discontinuous, the sensitivity exhibits Dirac delta functions. As said earlier, different choices are possible: some authors have tried to adapt their numerical schemes in order to deal with the Dirac functions [START_REF] Clarisse | A Godunov-type method in Lagrangian coordinates for computing linearly-perturbed planar-symmetric flows of gas dynamics[END_REF][START_REF] Després | Robust Uncertainty Propagation in Systems of Conservation Laws with the Entropy Closure Method[END_REF][START_REF] Godlewski | The linearized stability of solutions of nonlinear hyperbolic systems of conservation laws: A general numerical approach[END_REF], some others added a correction term to the sensitivity equations [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]. We decide to adopt the second strategy, as done in [START_REF] Chalons | Sensitivity analysis and numerical diffusion effects for hyperbolic PDE systems with discontinuous solutions. The case of barotropic Euler equations in Lagrangian coordinates[END_REF], that leads to an accurate sensitivity almost everywhere in the domain, except for the discontinuity points. The correction term that we add to the sensitivity equations has the following form:

S = Ns k=1 δ k ρ ρ ρ k , ( 13
)
where Ns is the number of discontinuities, which can be either shocks or contact discontinuities,

δ k = δ(x -x k,s (t))
is the Dirac delta function with x k,s (t) position of the k-th shock and ρ ρ ρ k is the amplitude of the k-th correction. To compute the amplitude ρ ρ ρ k (t), we start by integrating the sensitivity equations with the source term on a control volume which contains a single discontinuity travelling at speed σ k . As the control volume goes to zero, one has:

ρ ρ ρ k = (U - a,k -U + a,k )σ k + Fa(U + k , U + a,k ) -Fa(U - k , U - a,k ), (14)
where

U + k,a (respectively U - k,a
) is the value of the sensitivity to the right (respectively left) of the k-th discontinuity. Then, one writes the Rankine-Hugoniot relations associated with (3)

-σ k (U + k -U - k ) + F(U + k ) -F(U - k ) = 0, where U + k (respectively U - k )
is the value of the state to the right (respectively left) of the k-th discontinuity. If we differentiate these conditions with respect to a, we obtain:

(U - k,a -U + k,a )σ k + (U - k -U + k )σ k,a + σ k (∂xU + -k ∂xU - k )∂ax k,s (t) = = Fa(U - k , U - a,k ) -Fa(U + k , U + a,k ) + ∂F(U + k ) ∂U ∂xU + k - ∂F(U - k ) ∂U ∂xU - k ∂ax k,s (t), (15)
where σ k,a := ∂aσ k . Replacing (15) into [START_REF] Glimm | Solutions in the large for nonlinear hyperbolic systems of equations[END_REF], one obtains the following definition of ρ ρ ρ k , which does not depend on the sensitivity itself:

ρ ρ ρ k = (U + -U -)σ k,a + σ k (∂xU + -∂xU -)∂ax k,s (t) - ∂F(U + ) ∂U ∂xU + - ∂F(U -) ∂U ∂xU -∂ax k,s (t).
(16) The terms depending on ∂ax k,s (t) are very difficult to estimate. However, one can remark that all these terms contain ∂xU ± k : therefore, in the following sections, when we design first order finite volume schemes for the sensitivity we will consider the simpler expression

ρ ρ ρ k (t) = σ k,a (U + k -U - k ), (17)
since the solution, in a first order finite volume framework, is piecewise constant on the cells and therefore ∂xU ± k = 0. Replacing (17) into [START_REF] Duvigneau | An improved continuous sensitivity equation method for optimal shape design in mixed convection[END_REF] gives the following definition of the source term (valid only for piecewise constant functions):

S = Ns k=1 σ k,a (U + k -U - k )δ k . (18)
A special treatment, which will be detailed later, is necessary for a second or higher order discretisation, where the discrete solution is not constant within each cell and therefore ∂xU ± k = 0. The new system thus writes:

∂tU + ∂xF(U) = 0 ∂tUa + ∂xFa(U, Ua) = S. (19)
In the next section, we design a numerical scheme to approximate the solution of [START_REF] Gunzburger | Perspectives in flow control and optimization[END_REF]. The analytical solution for a given initial data of Riemann type is detailed in appendix A.

x j-1/2

x j+1/2 ∆x

x j C j

Numerical methods

In this section we consider the numerical approximation of [START_REF] Gunzburger | Perspectives in flow control and optimization[END_REF]. We derive first and second order Roe-type numerical schemes and we pay particular attention to the numerical diffusion effects induced by these approaches. Indeed and as we will see it may prevent the numerical solution from converging to the correct solution. We consider a uniform grid in space with a constant step ∆x, xj is the center of the j-th cell Cj, whose extrema are x j-1/2 and x j+1/2 (cf. Figure 1). We use an adaptive time step ∆t n , chosen according to a CFL condition, and the intermediate times are t n+1 = t n + ∆t n . We indicate with V n j = (U n j , U n a,j ) t the average value of the state and the sensitivity in the cell Cj at time t n .

We use Godunov-type schemes, which consist of two main steps: first, one solves the Riemann problem at each interface x j-1/2 at time t n , obtaining in this way a solution at time t n+1 , v(x, t n+1 ) = (u(x, t n+1 ), ua(x, t n+1 )) t ; the second step is to project v(x, t n+1 ) in order to obtain a piecewise constant solution on the mesh. How to compute v(x, t) is the topic of the next subsections: first we describe an approximate solver for the state and then two for the sensitivity. Different choices for the solution of the Riemann problem lead to different numerical schemes. Then, we explain two different projection techniques: the classical one and an anti-diffusive one. Finally, we explain how to extend the schemes to higher order in space, focusing in particular on the second order.

Riemann solver for the state

First, we consider the state system, for which the classical numerical schemes can be used: in this work we used the approximate Riemann solver of Roe, because it has the property of being exact for an isolated shock and we want to be as precise as possible in the shocks. In addition, we remark that it would not be possible to use a solver with only one intermediate star state, such as HLL (Harten, Lax and van Leer [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]), because of the definition of the source term [START_REF] Guinot | An approximate riemann solver for sensitivity equations with discontinuous solutions[END_REF]: two intermediate states are necessary in order to be able to compute the correction term across the contact discontinuity (cf. Figures 23for the structure of different solvers).

The main idea of the Roe scheme is to replace the Jacobian matrix A(U) in ( 4) with a constant matrix A(UL, UR), obtaining in this way a linearised system, whose solution to the Riemann problem can be computed exactly. For the Euler system, a proper linearisation is provided by Roe in the original paper [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF]. Furthermore, there is no need to assemble the matrix, it is sufficient to know its eigenvalues and eigenvectors, which are the following:

λ ROE 1 = ũ -c, λ ROE 2 = ũ, λ ROE 3 = ũ + c, r1 =   1 ũ - c H -ũc   , r2 =   1 ũ ũ2 2   , r3 =   1 ũ + c H + ũc   .
The quantities denoted with a tilde are Roe averaged quantities defined as follows:

ũ = √ ρL uL + √ ρR uR √ ρL + √ ρR , H = √ ρL HL + √ ρR HR √ ρL + √ ρR , c = (γ -1) H - 1 2 ũ2 .
Therefore, the Roe solver consists of four constant states (UL, U * L , U * R , and UR, cf. The relation ( 20) is used to compute the coefficients αi, then one has:

∆U = UR -UL = 3 i=1 αiri. ( 20
) λ ROE 1 λ ROE 2 λ ROE 3 x t x c U L U * L U * R U R
λ ROE 1 λ ROE 3 x t x c U a,L U * a U a,R
U * L = UL + α1r1 = UR -α2r2 -α3r3, U * R = UR -α3r3 = UL + α1r1 + α2r2. (21)
Once all the quantities U * L , U * R , and λ ROE are known at each interface x j-1/2 , u(x, t n+1 ) can be built by juxtaposition of the solutions of each Riemann problem.

It is well known that, in case of transonic rarefaction, the Roe solver provides a non-entropic solution. To overcome this problem, we implemented the entropic fix proposed in [START_REF] Harten | Self adjusting grid methods for one-dimensional hyperbolic conservation laws[END_REF].

Riemann solvers for the sensitivity

For the sensitivity we propose two different strategies. Indeed and as explained in the previous section, for the state it is necessary to use a Riemann solver with two different star states, in order to be able to compute the source term across the contact discontinuity. However, for the sensitivity an HLL-type approach can be used, which gives a first strategy. Another possible strategy is to keep for the sensitivity the same structure as for the state, and therefore to have an HLLC-type scheme (Harten, Lax and van Leer Contact [START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF]). A third possibility which we will not analyse here, explored in detail in [START_REF] Appel | Sensitivity calculations for conservation laws with application to discontinuous fluid flows[END_REF], is to rewrite the sensitivity flux in such a way that the same Roe Riemann solver used for the state can be applied for the sensitivity. Let us now describe the two possibilities considered in detail.

HLL-type scheme

The first Riemann solver proposed for the sensitivity has a simpler structure than the state solver: we neglect the contact discontinuity, therefore the solver consists only of three constant states (Ua,L, U * a , and Ua,R) connected by two discontinuities travelling at speeds λ ROE 1 and λ ROE 3 (cf. Figure 3). The star value of the sensitivity U * a at the interface j -1/2 can be computed directly from the Harten, Lax and van Leer conditions [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] applied to system of conservation laws with source terms. We get:

U * a,j-1/2 = 1 λ ROE 3 -λ ROE 1 λ ROE 3 U n a,j -λ ROE 1 U n a,j-1 -Fa(Uj, Ua,j) + Fa(Uj-1, Ua,j-1) + S j-1/2 , (22)
where the source term is discretised as follows:

S j-1/2 = ∂aλ ROE 1,j-1/2 (U * L,j-1/2 -Uj-1)d 1,j-1/2 + ∂aλ ROE 2,j-1/2 (U * R,j-1/2 -U * L,j-1/2 ) +∂aλ ROE 3,j-1/2 (Uj -U * R,j-1/2 )d 3,j-1/2 ,
where d ,j-1/2 are shock detectors, d ,j-1/2 = 1 if there is an -shock at the interface j -1/2, it is zero otherwise. They are based on the fact that the velocity u is always decreasing across a shock, whilst the density ρ is increasing across a 1-shock and it is decreasing across a 3-shock:

d 1,j-1/2 = 1 if ρj > ρj-1 and uj < uj-1, 0 otherwise, d 3,j-1/2 = 1 if ρj < ρj-1 and uj < uj-1, 0 otherwise.
Furthermore, we remark that there is no need for a contact detector because it is known that the middle wave is always a contact discontinuity. Such a discretisation of the source term comes directly from [START_REF] Guinot | Sensitivity equations for hyperbolic conservation law-based flow models[END_REF], if one considers the fact that a Riemann problem can have at most three discontinuities.

HLLC-type scheme

Another possible approach for the sensitivity is to keep the same structure as for the state (cf. Figure 2), with the same speeds of propagation for the three discontinuities. We need to compute the two intermediate constant states U * a,L and U * a,R . Again, a possible strategy to compute U * a,L and U * a,R is to follow the Harten, Lax and van Leer formalism with source term and to impose the following linear system, made of Rankine-Hugoniot jump relations:

                                   -λ1(ρ * a,L -ρa,L) + (ρu) * a,L -(ρu)a,L = ∂aλ1(ρ * L -ρL), -λ2(ρ * a,R -ρ * a,L ) + (ρu) * a,R -(ρu) * a,L = ∂aλ2(ρ * R -ρ * L ), -λ3(ρa,R -ρ * a,R ) + (ρu)a,R -(ρu) * a,R = ∂aλ3(ρR -ρ * R ), (γ-3) 2 ũ2 (ρ * a,R -ρ * a,L ) + (2 -γ)ũ((ρu) * a,R -(ρu) * a,L ) +(γ -1)((ρE) * a,R -(ρE) * a,L ) = ∂aλ2((ρu) * R -(ρu) * L ), (λ2 -λ1)(ρu) * a,L + (λ3 -λ2)(ρu) * a,R + λ1(ρu)a,L -λ3(ρu)a,R +Fa,R|2 -Fa,L|2 = ∆xS|2, (λ2 -λ1)(ρE) * a,L + (λ3 -λ2)(ρE) * a,R + λ1(ρE)a,L -λ3(ρE)a,R +Fa,R|3 -Fa,L|3 = ∆xS|3, ( 23
)
where λ1 = ũ -c, λ2 = ũ, and λ3 = ũ + c. The first three equations are the Rankine-Hugoniot condition on ρ across the three waves, differentiated with respect to a. Note that summing up these equations gives the integral condition of the Harten, Lax and van Leer formalism of the density variable. The fourth equation is the Rankine-Hugoniot condition on ρu for the linearised system differentiated with respect to a; the last two equations are the integral conditions on the sensitivities (ρu)a and (ρE)a. If we define the following vectors

x = (ρ * a,L , ρ * a,R , (ρu) * a,L , (ρu) * a,R , (ρE) * a,L , (ρE) * a,R ) t b =         b1 b2 b3 b4 b5 b6         =         ∂aλ1(ρ * L -ρL) + (ρu)a,L -λ1ρa,L ∂aλ2(ρ * R -ρ * L ) ∂aλ3(ρR -ρ * R ) -(ρu)a,R + λ3ρa,R ∂aλ2((ρu) * R -(ρu) * L ) ∆xS|2 -λ1(ρu)a,L + λ3(ρu)a,R -Fa,R|2 + Fa,L|2 ∆xS|3 -λ1(ρE)a,L + λ3(ρE)a,R -Fa,R|3 + Fa,L|3        
the system can be rewritten as:

Ax = b,
where A is the following matrix:

A =         -λ1 0 1 0 0 0 λ2 -λ2 -1 1 0 0 0 λ3 0 -1 0 0 -(γ-3) 2 ũ2 (γ-3) 2 ũ2 -(2 -γ)ũ (2 -γ)ũ -(γ -1) (γ -1) 0 0 c c 0 0 0 0 0 0 c c        
and we have det(A) = 4c 4 (γ -1) = 0. The solution of the system has the following form:

x =                    (2c+ũ)b 1 +(c+ũ)b 2 +ũb 3 -b 5 2c 2 -ũb 1 +(c-ũ)b 2 +(2c-ũ)b 3 +b 5 2c 2 (ũ 2 +cũ)b 1 +(ũ 2 -c 2 )b 2 +(ũ 2 -c)ũb 3 +(c-ũ)b 5 2c 2 -(ũ 2 +cũ)b 1 +(c 2 -ũ 2 )b 2 +(cũ-ũ 2 )b 3 +(c+ũ)b 5 2c 2 (γ-1)(ũ 3 +cũ 2 )b 1 +((γ-1)ũ 3 +2(2-γ)c 2 ũ)b 2 +(γ-1)(ũ 3 -cũ 2 )b 3 -2c 2 b 4 -(γ-1)ũ 2 b 5 +2(γ-1)cb 6 4(γ-1)c 2 -(γ-1)(ũ 3 +cũ 2 )b 1 -((γ-1)ũ 3 +2(2-γ)c 2 ũ)b 2 +(γ-1)(cũ 2 -ũ 3 )b 3 +2c 2 b 4 +(γ-1)ũ 2 b 5 +2(γ-1)cb 6 4(γ-1)c 2                   
An alternative strategy to compute U * a,L and U * a,R is to differentiate with respect to a the following relations:

U * L = UL + α1r1, U * R = UR -α3r3, (24)
obtaining

U * a,L = Ua,L + αa,1r1 + α1ra,1, U * a,R = Ua,R -αa,3r3 -α3ra,3, (25)
with

r1 =   1 ũ - c H -ũc   , ra,1 =   0 ũa -ca Ha -ũac -ũca   , r2 =   1 ũ ũ2 2   , ra,2 =   0 ũa ũũa   , r3 =   1 ũ + c H + ũc   , ra,3 =   0 ũa + ca Ha + ũac + ũca   ,          α2 = γ -1 c2 (ρR -ρL)( H -ũ2 ) + ũ (ρu)R -(ρu)L -(ρE)R -(ρE)L , α1 = 1 c (ρR -ρL)(ũ + c) -(ρu)R -(ρu)L -cα2 , α3 = (ρR -ρL) -(α1 + α2),                                        αa,2 = - 2ca(γ -1) c3 (ρR -ρL)( H -ũ2 ) + ũ (ρu)R -(ρu)L -(ρE)R -(ρE)L + γ -1 c2 (ρa,R -ρa,L)( H -ũ2 ) + (ρR -ρL)( Ha -2ũũa) +ũa (ρu)R -(ρu)L -(ρE)R -(ρE)L + ũ (ρu)a,R -(ρu)a,L -(ρE)a,R -(ρE)a,L , αa,1 = - ca c2 (ρR -ρL)(ũ + c) -((ρu)R -(ρu)L) -cα2 + 1 2c (ρa,R -ρa,L)(ũ + c) + (ρR -ρL)(ũa + ca) -((ρu)a,R -(ρu)a,L) -caα2 -cαa,2 , αa,3 = (ρa,R -ρa,L) -(αa,1 + αa,2).
The next proposition states that the two strategies to define U * a,L and U * a,R are equivalent. Proposition. 2. The star sensitivities (25) solve the system [START_REF] Iooss | A Review on Global Sensitivity Analysis Methods[END_REF].

Proof. We will prove that the star sensitivities defined in (25) satisfy the system (23).

1. First equation. Writing the first coefficient of (24) one easily finds ρ * L -ρL = α1, and writing the first two coefficient of (25) one finds:

ρ * a,L -ρa,L = αa,1, (ρu) * a,L -(ρu)a,L = αa,1(ũ -c) + α1(ũa -ca).
We now replace these three expressions in the first equation of ( 23) and we obtain:

-λ1αa,1 + αa,1(ũ -c) + α1(ũa -ca) = ∂aλ1α1,
which is always verified, since λ1 = ũ -c.

2. Second equation. We recall that

UR -UL = 3 i=1 αiri, Ua,R -Ua,L = 3 i=1 αa,iri + αira,i.
Therefore, one has:

U * R -U * L = α2r2, U * a,R -U * a,L = αa,2r2 + α2ra,2,
which gives us the following relations:

ρ * R -ρ * L = α2, ρ * a,R -ρ * a,L = αa,2, (ρu) * a,R -(ρu) * a,L = αa,2 ũ + α2 ũa.
We now replace them in the second equation of ( 23) and we obtain:

-λ2αa,2 + αa,2 ũ + α2 ũa = ∂aλ2α2,
which is always verified, since λ2 = ũ.

3. Third equation. As we did for the first two equations, one can find the three following expressions:

ρR -ρ * R = α3, ρa,R -ρ * a,R = αa,3, (ρu)a,R -(ρu) * a,R = αa,3(ũ + c) + α2(ũa + ca).
By replacing them in the third equation of ( 23) one can easily check that the equation is always verified, since λ3 = ũ + c.

4. Fourth equation. As we did for the previous equations, one can find the three following expressions:

(ρu) * R -(ρu) * L = α2 ũ, ρ * a,R -ρ * a,L = αa,2, (ρu) * a,R -(ρu) * a,L = αa,2 ũ + α2 ũa, (ρE) * a,R -(ρE) * a,L = αa,2 ũ2 2 + α2 ũũa.
By replacing them in the fourth equation of ( 23) one can easily check that the equation is always verified, since λ2 = ũ.

5. Fifth and sixth equations. The last two equations are the last two components of the following vectorial equation:

(λ2 -λ1)U * a,L + (λ3 -λ2)U * a,R + λ1Ua,L -λ3Ua,R + Fa,R -Fa,L = ∆xS,
which can be rewritten as:

λ1(Ua,L -U * a,L ) + λ2(U * a,L -U * a,R ) + λ3(U * a,R -Ua,R) + Fa,R -Fa,L = ∆xS.
Replacing the definitions (25) one finds:

-λ1(αa,1r1 + α1ra,1) -λ2(αa,2r2 + α2ra,2) -λ3(αa,3r3 + α3ra,3) + Fa,R -Fa,L = ∆xS.

We recall that by definition of Roe fluxes, one has:

FR -FL = 3 i=1 αiλiri ⇒ Fa,R -Fa,L = 3 i=1 αa,iλiri + αiλa,iri + αiλira,i.
Therefore, we obtain:

∆xS = 3 i=1 αiλa,iri = λa,1(U * L -UL) + λa,2(U * R -U * L ) + λa,3(UR -U * R ),
which is consistent with our discretisation of the source term.

Finally, one can obtain ua(x, t n+1 ) by juxtaposition, as we did for the state.

Projection step

The projection step is usually performed by averaging on the cell the solution v(x, t n+1 ), whose components u(x, t n+1 ) and ua(x, t n+1 ) can be computed as described in the previous sections.

V n+1 j = 1 ∆x x j+1/2 x j-1/2 v(x, t n+1 )dx. ( 26
)
We remark that the integral ( 26) is easy to compute, v being piecewise constant. However, this projection method introduces numerical diffusion. As shown in [START_REF] Chalons | Sensitivity analysis and numerical diffusion effects for hyperbolic PDE systems with discontinuous solutions. The case of barotropic Euler equations in Lagrangian coordinates[END_REF], numerical diffusion plays a fundamental role in the discretisation of the sensitivity, especially across shocks. For this reason, we propose another projection method, introduced in [START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF] and inspired by Glimm's method [START_REF] Glimm | Solutions in the large for nonlinear hyperbolic systems of equations[END_REF][START_REF] Chorin | Random choice solution of hyperbolic systems[END_REF]. First, we define a staggered mesh, whose cells will be denoted C n j , as follows:

C n j = (x n j-1/2 , xn j+1/2 ), xn j-1/2 = x j-1/2 + σ n j-1/2 ∆t n ,
where σ n j-1/2 is a proper speed, defined in order to avoid averaging across a shock. Numerical results show that there is no need to move the mesh for the contact discontinuity (cf. section 4). The definition of σ n j-1/2 is the following:

σ n j-1/2 =      λ ROE 1,j-1/2 if d 1,j-1/2 = 1, λ ROE 3,j-1/2 if d 3,j-1/2 = 1, 0 otherwise,
where d ,j-1/2 are the shock detectors defined earlier.

The second step is to perform the average on the staggered mesh, obtaining in this way an intermediate solution V

n+1 j : V n+1 j = 1 ∆x n j xj+1/2 xj-1/2 v(x, t n+1 )dx, (27)
where ∆x n j = xj+1/2 -xj-1/2 . Finally, the last step is a sampling step, in order to go back to the initial uniform grid. Let (βn) be a random sequence varying in (0, 1), for instance βn ∼ U([0, 1]); then:

V n+1 j =      V n+1 j-1 if βn+1 ∈ 0, ∆t ∆x max(σ n j-1/2 , 0) , V n+1 j if βn+1 ∈ ∆t ∆x max(σ n j-1/2 , 0), 1 + ∆t ∆x min(σ n j+1/2 , 0) , V n+1 j+1 if βn+1 ∈ 1 + ∆t ∆x min(σ n j+1/2 , 0), 1 . (28)
We remark that one βn is drawn at each time step and it is the same for all the cells. The method is proven to be convergent even if a low discrepancy deterministic sequence is used. In this work, we use the van der Corput sequence (cf. [START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF]):

βn = m k=0 i k 2 -(k+1) , n = m k=0 i k 2 k ,
where i k = 0, 1 is the binary expansion of the integers.

Second order extension

In this section, we extend to the second order the schemes presented above. In time, we use a standard two steps Runge-Kutta method, whilst in space we use a MUSCL-type (Monotonic Upstream-centered Scheme for Conservation Laws, [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF]) approach, inspired from [START_REF] Chalons | Sensitivity analysis and numerical diffusion effects for hyperbolic PDE systems with discontinuous solutions. The case of barotropic Euler equations in Lagrangian coordinates[END_REF]. In a few words (we refer to [START_REF] Chalons | Sensitivity analysis and numerical diffusion effects for hyperbolic PDE systems with discontinuous solutions. The case of barotropic Euler equations in Lagrangian coordinates[END_REF] for more details) the main idea of a MUSCL-type scheme is to consider in replacement of a constant value V n j in each cell, a higher order polynomial V n j (x), x ∈ [x j-1/2 , x j+1/2 ]. The edge values V n j (x j+1/2 ), V n j+1 (x x+1/2 ) are used as left and right values for the Riemann problem at the interface j + 1/2; the Riemann problem is then solved as explained in the previous section. However, the definition of the source term ( 13)-( 17) is valid only if the state is piecewise constant (cf. [START_REF] Chalons | Sensitivity analysis and numerical diffusion effects for hyperbolic PDE systems with discontinuous solutions. The case of barotropic Euler equations in Lagrangian coordinates[END_REF]). Therefore, we suggest to consider a piecewise constant state on half of each cell: these two constant values will be denoted V n j±1/4 and correspond to the edge values V n j (x j±1/2 ) (see Figure 4). In this work, we compute the edge values with a standard approach: and usual choice for ∆V n j is to use a slope-limiter procedure, for instance: We remark that this approach leads to an additional Riemann problem in the middle of the cell: in this way we are able to extend the scheme accuracy to second-order, while keeping piecewise constant representations in the cells, which is necessary to make the terms containing ∂xU ± vanishing in (15).

V n j±1/4 = V n j ± ∆V n j , x x j-3/2 x j-1/2 x j+1/2 x j+3/2 V j-5/4 V j-3/4 V j-1/4 V j+1/4 V j+3/4 V j+5
∆V n j = 1 2 minmod(V n j+1 -V n j , V n j -V n j-

Summary

Here, we briefly sum up all the ingredients introduced in this section and we clarify how they can be combined to obtain different numerical schemes:

• Order of the scheme. In this paper we focused on first and second order schemes, but higher order can be used bearing in mind that the state needs to be piecewise constant for the definition of the source term to be valid.

• Riemann solver for the state. In this paper we proposed the Roe Riemann solver for the state, but the only constraint is to use a solver with two intermediate states U * L and U * R (for instance, HLL could not be used for the state).

• Riemann solver for the sensitivity. We designed two different numerical schemes for the sensitivity: an HLL and HLLC-type scheme.

• Type of projection. Either the classical projection or the anti-diffusive projection can be used. Numerical results in the next section will show that, for this problem, diffusive scheme does not converge to the analytical solution.

Finally, we remark that these four choices are independent of each other. In table [START_REF] Abgrall | A semi-intrusive deterministic approach to uncertainty quantification in non-linear fluid flow problems[END_REF] we summarize all the combinations used in this work, with the labels used in the next sections.

Convergence tests for the numerical schemes

We consider the Riemann problem described in appendix A. The initial data for the state on the physical variables is the following: ρL = 1, uL = 0, pL = 1, ρR = 0.125, uR = 0, pR = 0.1.

We consider as parameter of interest a = pL, therefore the initial data for the sensitivity is:

ρa,L = ρa,R = ua,L = ua,R = pa,R = 0, pa,L = 1.

In Figures 5-6-7 we show the convergence of the different numerical schemes presented in Section 3.

Figure 5 shows the convergence for the state: the rate of convergence is the expected one; one can remark that the antidiffusive schemes are slightly less precise than the diffusive ones. In Figures 67we plot the error for the sensitivity, first with the HLL-type scheme (Figure 6) and then with the HLLC-type scheme (Figure 7): considering two different star regions for the sensitivity does not seem to make much difference; however one can remark the same effect shown in [START_REF] Chalons | Sensitivity analysis and numerical diffusion effects for hyperbolic PDE systems with discontinuous solutions. The case of barotropic Euler equations in Lagrangian coordinates[END_REF] for a simpler system: the diffusive schemes do not converge for the sensitivity, this is especially evident for the variable ρa.

In Figure 8 we plot the solution at the final time T = 0.1, obtained with a mesh ∆x = 10 -3 with the first order schemes, both diffusive and antidiffusive (for the sensitivity, the HLL-type scheme has been used): one can notice that the plateau in the right-star zone is not properly captured by the diffusive scheme. This does not change as one refines the mesh, nor with a higher order scheme, as one can see from Figure 9. In Figure 10 we compare the antidiffusive schemes, first and second order: for the state, the difference is noticeable mainly in the contact discontinuity (therefore only for ρ), whilst for the sensitivity the difference is significant in the neighbourhood of the discontinuities before and after the rarefaction. Finally, in Figure 11 we compare the HLL and the HLLC-type schemes for the sensitivity: as anticipated by the error plots, the two schemes are almost equivalent in terms of results; the difference between the solutions provided by the two second order schemes in L ∞ -norm is 0.0096, 0.0139, and 0.0062 respectively for ρa, ua, and pa and this is why they are almost indistinguishable in Figure 11. For this reason, the use of HLL-type scheme is preferable, being less expensive from a computational point of view and less complicated to implement. 5 Uncertainty Quantification

Problem description

In this section, we show how the estimated sensitivities can be used for uncertainty propagation, i.e. estimate statistical moments of the solution accounting for some uncertain parameters [START_REF] Putko | Approach for uncertainty propagation and robust design in cfd using sensitivity derivatives[END_REF][START_REF] Turgeon | Sensitivity and uncertainty analysis for variable property flows[END_REF][START_REF] Delenne | Propagation de la sensibilité dans les modèles hydrodynamiques[END_REF]. We want to provide a demonstration of one of the many possible applications of the CSE method, in order to underline the potentials of the proposed approach and its limitations, too. We also intend to quantify the impact of removing the spikes from the sensitivity solution. Many techniques have been developed during the last decades to propagate uncertainty through PDE models: these methods can be either probabilistic or deterministic. The proposed method based on derivatives estimation falls into the second category, while the most well-known of these techniques, the Monte Carlo method, is in the first. Other techniques are for instance polynomial chaos [START_REF] Walters | Towards stochastic fluid mechanics via polynomial chaos[END_REF][START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF][START_REF] Knio | Uncertainty propagation in CFD using polynomial chaos decomposition[END_REF][START_REF] Després | Robust Uncertainty Propagation in Systems of Conservation Laws with the Entropy Closure Method[END_REF], or the random space partition [START_REF] Abgrall | A semi-intrusive deterministic approach to uncertainty quantification in non-linear fluid flow problems[END_REF]. A very good review and comparison of many techniques with applications to fluid dynamics can be found in [START_REF] Walters | Uncertainty analysis for fluid mechanics with applications[END_REF]. A typical objective of uncertainty propagation is to determine a confidence interval for the output of a model, in our case U, given the uncertainty on the input parameters. This estimation is part of the broader domain of Uncertainty Quantification (UQ), which also includes the identification of the most critical uncertain parameters, their ranking, and the analysis of the variability of the output.

In this work, we compare the Monte Carlo approach and a sensitivity-based estimation. In the following, X will stand for one of the variables, i.e. X can either be ρ, u or p, and Xa the corresponding sensitivity. We use the notation µX to indicate the average 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:

CIX = [µX -κσX , µX + κσX ]. (29)
We remark that ( 29) is valid only for gaussian data. 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 κ = 1.96 provides a 95% confidence interval, while κ = 2.58 a 99% one.

Monte Carlo method.

Here we briefly introduce the Monte Carlo method, for more details see for instance [START_REF] Christian | Monte Carlo statistical methods[END_REF]. 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 random simulations. Let a be the vector of uncertain parameters, with a known distribution. Then, N random samples ai are drawn from this distribution, and for each ai the corresponding solution Xi is computed. Then, the unbiased average and variance estimators are used:

µX = 1 N N i=1 Xi, σ 2 X = 1 N -1 N i=1 (µX -Xi) 2 .
These estimates are good 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-based method. Once the sensitivities of the solution with respect to the input parameters are known, a deterministic estimation of the average µX and of the variance σ 2 X of the output X can be easily obtained. Let µa be the average of the uncertain vector a and σa the covariance matrix:

µa =    µa 1 . . . µa M    , σa =      σ 2 a 1 cov(a1, a2) . . . cov(a1, aM ) cov(a1, a2) σ 2 a 2 . . . cov(a2, aM ) . . . . . . . . . cov(a1, aM ) . . . σ 2 a M     
, where M is the number of uncertain parameters, µa i the average of the i-th parameter, σ 2 a i its variance and cov(•, •) the covariance. Let us consider the first order Taylor expansion for the variable X with respect to the vector of parameters a:

X(a) = X(µa) + M i=1 (ai -µa i )Xa i (µa) + o( a 2 ).
Then computing the average, since X(µa) and Xa i (µa) are not random variables, at first order one gets:

µX = E[X(a)] = X(µa) + M i=1 Xa i (µa)E[ai -µa i ] = X(µa),
because E[(ai -µa i )] = 0. In the same way, one can compute the variance:

σ 2 X = E[(X(a) -µX ) 2 ] = E   M i=1 Xa i (µa)(ai -µa i ) 2   = = M i=1 X 2 a i (µa)E[(ai -µa i ) 2 ] + M i,j=1 i =j Xa i (µa)Xa j (µa)E[(ai -µa i )(aj -µa j )].
Therefore, we 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 a i σ 2 a i + M i,j=1 i =j
Xa i Xa j cov(ai, aj).

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].

Numerical results

We applied the uncertainty propagation techniques described in the previous subsection to the test case described in appendix A. The uncertain parameters are the left and right values of the physical variables for the state, i.e.: a = (ρL, ρR, uL, uR, pL, pR) t , and have a Gaussian distribution with the following average and covariance matrix: µa = (1, 0.125, 0, 0, 1, 0.1) t , σa = diag(0.001, 0.000125, 0.0001, 0.0001, 0.001, 0.0001).

This choice means that all the parameters are uncorrelated and we chose as their variance the 0.1% of their average, except for the velocity, whose average is 0. We remark that the fact that the parameters follow a Gaussian distribution does not say anything about the distribution of the output of the model. Therefore, before using the confidence interval [START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF], one should check that the output is Gaussian, too. In Figures-12-13-14, three histograms are shown for each physical variable: for x = 0.35 (i.e. in the middle of the rarefaction wave), for x = 0.6 (i.e. in the middle plateau) and for x = 0.85 (i.e. close to the shock position). The histograms are obtained by computing the analytical solution for 5000 different values of the parameter vector a, which are sampled from its distribution. The histograms are then normalised with respect to the probability density function, using the MATLAB option 'normalization', 'pdf'. As one can see, the output is Gaussian when far from the shock; close to the shock, two distinct groups of points can be identified. The distribution predicted with the SA is drawn in red: the prediction is wrong in the neighbourhood of the shock, as as; it fits perfectly in the plateau and it is slightly shifted for the rarefaction. This shift is due to the prediction for the average, which is based entirely on the state and caused by the fact that it is only a first order approximation. However, one can remark that the variance is correctly estimated using the sensitivities. In the following, we use the expression [START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF], expecting however a loss of precision in the neighbourhood of the shock.

In Figure 15 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 = 1000 samples, on a mesh with ∆x = 10 -3 using a Roe first order diffusive scheme. As one can see, the average process smudges the shock and the standard deviation is bigger in that zone. In fact, this area of large variance around the shock location is related to the delta Dirac function in sensitivity solution, but this peak is "smooth" due to the fact that Monte-Carlo approach accounts for flow non-linearities and is not limited to the first-order estimate of the perturbation. In Figures 1617we show the results of the sensitivity-based approach, with ∆x = 10 -3 and the diffusive first order scheme, when the sensitivity is computed without the correction term (13): the spikes in the neighbourhood of the shock are very different with respect to the ones we get with the Monte Carlo approach. As explained above, these spikes are due to the default of high-order terms in the Taylor expansion. On one hand, these peaks lead to nonphysical values for the solution (in particular, the confidence intervals contains negative values for the pressure and for the density); on the other hand, they do not enlarge sufficiently the zone to contain the majority of the samples: one can observe that four out of five samples fall out of the predicted interval in the neighbourhood of the shock. Nevertheless, we observe that the non-linear effects are located in a small region around the shock and elsewhere confidence interval is well estimated by the sensitivity-based method. The results obtained with the corrected sensitivities are shown in Figure 18: the confidence interval obtained correspond to the ones obtained with the Monte Carlo approach, apart for the shock zone. Of course, the sensitivity-based approach does not capture the uncertainty in the neighborhood of the shock, because it neglects the dependence of the speed of the shock on the parameters. This is why most of the samples fall out of the zone predicted with the sensitivity-based approach, and it is the case with and without correction. However, the correction avoids non-physical values in the confidence interval. Although the sensitivity-based approach is not able to account for non-linear effects, contrary to the Monte-Carlo method, it is far less expensive: the Monte Carlo approach requires 1000 solutions of the state, whilst the sensitivity-based approach only one solution of the state and as many solution of the sensitivity as the number of uncertain parameters, in this case 6. Therefore, this approach can still be interesting for computationally-demanding problems, for which the use of the Monte-Carlo method cannot be envisaged. Finally, in Figure 19 we show the results obtained with the anti-diffusive scheme: the difference with respect to the diffusive scheme is not significant. This is a good news for possible future developments in 2D: the anti-diffusive scheme is very difficult to adapt in higher dimensional spaces; in fact the Glimm method has been proven not to work in a two-dimensional space. With these results, we underline how the numerical diffusion plays an important role in the convergence of the scheme, but it is not so significant for the final application.

The numerical results of the SA method applied to an uncertainty propagation problem show the potential and the limits of the method: it is really affordable from a computational point of view, with the trade off of being less precise than, for instance, a Monte Carlo method in the discontinuous zones. If one accepts this loss of precision, the proposed approach remains accurate in the regular zone and is highly competitive thanks to its very low computational cost.

Conclusion

In this work, we extended to the complete Euler system the method proposed in [START_REF] Chalons | Sensitivity analysis and numerical diffusion effects for hyperbolic PDE systems with discontinuous solutions. The case of barotropic Euler equations in Lagrangian coordinates[END_REF] for the p-system. The definition of the source term does not differ significantly and the same shock detectors can be used in this case. We remark that a contact discontinuity detector is not necessary, since the middle wave is always a contact discontinuity. However, in this more complex case, the form of the proposed source term precludes the application of some well-known and widely used numerical schemes such as the HLL scheme, and all the approximate Riemann solvers with only one middle state. The numerical results show that the numerical diffusion in the shocks plays an important role and corrupts the convergence to the correct solution even for the complete system. However, we remark that the expected convergence rate can be achieved without removing the numerical diffusion in the contact discontinuity, which simplifies the definition of the staggered mesh. Currently, we are extending this Figure 19: SA approach with correction, anti-diffusive scheme. Average and the average plus and minus twice the standard deviation in red. Five samples in black dashed lines to the quasi 1D Euler system and we are dealing with some applications, such as optimization and uncertainty quantification: the results obtained in those frameworks show the importance of the correction term. We consider as parameter of interest a = pL, therefore the initial data for the sensitivity is:

ρa,L = ρa,R = ua,L = ua,R = pa,R = 0, pa,L = 1.

This choice of initial data leads to the structure in Figures 20-21, for the state as well as for the sensitivity: the 1-wave is a rarefaction and the 3-wave is a shock. For the notation, please refer to Figure 20 for the state and Figure 21 for the sensitivity. Let us now give the exact formulas for the state and for the sensitivity.

State solution: the exact solution for the physical variables is given in [START_REF] Appel | Sensitivity calculations for conservation laws with application to discontinuous fluid flows[END_REF]. Every variable is given as a function of the pressure in the right-star zone p * R , which is computed numerically from the following implicit relation:

pL = p * R     1 - (γ -1) c R c L ( p * R p R -1) 2γ 2γ + (γ + 1)( p * R p R -1)     -2γ γ-1 , (30)
where c = γp ρ , with = L, R. In the star regions, we have:

p * L = p * R = p * , u * L = u * R = u * = cR p * pR -1 2 γ(γ + 1) p * p R + γ(γ -1)
, because the velocity u and the pressure p are Riemann invariants across the 2-wave; as for the density ρ, we have: .

Finally, the solution writes:

U(x, t) =                  UL x -xc < -cLt, U -cLt < x -xc < γ+1 2 u * -cL t, U * L γ+1 2 u * -cL t < x -xc < u * t, U * R u * t < x -xc < cR γ-1 2γ + γ+1 2γ p * p R t, UR x -xc > cR γ-1 2γ + γ+1 2γ p * p R t. (31)
Sensitivity solution: here we are solving the second part of system 19, that is the one with the source term. The source term was designed in such a way that the solution for the sensitivity is the derivative of the state solution in the regular zones and there are no Dirac delta function where the state is discontinuous. Therefore, by differentiating [START_REF] Turgeon | Sensitivity and uncertainty analysis for variable property flows[END_REF] with respect to a, one obtains the following explicit formula for p * a,R : .

p * a,R = p * a,L = p * a = 1 + Θ 1-3γ γ-1 Ξp * Θ -2γ γ-1 + Θ 1-3γ γ-1 (Λ -Ψ)p *
The sensitivity has the same structure as the state, therefore:

Ua(x, t) =                  Ua,L x -xc < -cLt, Ua x-xc t -cLt < x -xc < γ+1 2 u * -cL t, U * a,L γ+1 2 u * -cL t < x -xc < u * t, U * a,R u * t < x -xc < cR γ-1 2γ + γ+1 2γ p * p R t, Ua,R x -xc > cR γ-1 2γ + γ+1 2γ p * p R t. (32)
We remark that if one writes the Rankine-Hugoniot conditions across the shock one finds:

-cR γ -1 2γ + γ + 1 2γ

p * pR (Ua,R -U * a,R ) + Fa(UR, Ua,R) -Fa(U * R , U * a,R ) = S.

Figure 1 :

 1 Figure 1: Spatial discretisation.

Figure 2 )

 2 connected by three discontinuities travelling at speeds λ ROE i . To compute the star states U * L and U * R , first we decompose the jump UR -UL along the eigenvectors of the Jacobian matrix A:

Figure 2 :

 2 Figure 2: Structure of the Roe solver for the state.

Figure 3 :

 3 Figure 3: Structure of the HLL-type solver for the sensitivity.

  1 ), where minmod(a, b) = sgn(a) min(|a|, |b|) if ab > 0, 0 otherwise.

1 Figure 5 : 1 Figure 6 : 1 Figure 7 :Figure 8 :

 1516178 Figure 5: Convergence test for the state. The Roman numerals I and II stand for the order of the scheme. AD stands for anti-diffusive.

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure9: Second order schemes, with and without numerical diffusion. HLL-type scheme for the sensitivity.

Figure 12 :Figure 13 :

 1213 Figure 12: Distribution of ρ(x, T ; a) for three different values of x at final time T = 0.2

Figure 14 :Figure 15 :

 1415 Figure 14: Distribution of p(x, T ; a) for three different values of x at final time T = 0.2

Figure 16 :Figure 17 :Figure 18 :

 161718 Figure 16: SA approach without correction. Average and the average plus and minus twice the standard deviation in red. Five samples in black dashed lines

Figure 20 :Figure 21 :

 2021 Figure 20: Structure of the solution for the Riemann problem for the state.

.

  In the rarefaction wave, we have:û(x, t) = 2(u * -uL) (γ + 1)u * x -xc t + 2 cLu * -uL cL -γ+1 2 u * (γ + 1)u * , ρ(x, t) = ρL 1 -(γ -, t) = pL 1 -(γ -1) û(x, t) 2cL 2γ γ-1

- 1 )cR p * p R - 1 cL 4γ 2 - 1 ,- 1 3 2 .

 112112 + 2γ(γ -1) p * p R + 1)cR p * p R -1 cLpR √ 2γ 2γ + (γ + 1) p * p RIn the star regions, by differentiating the corresponding state, one finds:

Table 1 :

 1 Figure 4: MUSCL discretisation. Dashed red line: first order discretisation. Dotted blue line: classical second order discretisation. Solid black line: second order discretisation used in this work. Numerical schemes summary.

	/4

A Solution of the Riemann problem

In this appendix, we write the exact solution for the system [START_REF] Gunzburger | Perspectives in flow control and optimization[END_REF] in a specific case (cf. [START_REF] Appel | Sensitivity calculations for conservation laws with application to discontinuous fluid flows[END_REF]), which was used as a test case to check the convergence of the numerical schemes proposed. We consider a Riemann problem, i.e.:

The general solution for this kind of problem is quite complicated, especially for the sensitivity (the last three components of V). First, we study the state (the first three components of V): the pair (λ2, r2) is linearly degenerated, i.e. ∇λ2 • r2 = 0, therefore the middle wave is always a contact discontinuity; concerning the 1-wave and the 3-wave, they are genuinely nonlinear therefore they can either be shocks or rarefaction waves. In Figure 20 we show the structure of the state in the case rarefaction-contact-shock. Concerning the sensitivity, it has the same structure as the state (cf. Figure 21 in the case rarefaction-contact-shock): the middle wave is always a contact wave, and the 1-and 2-wave are of the same type as for the state. The only difference is that the sensitivity presents discontinuities in the two extrema of the rarefaction fan (and this is why in Figure 21 the external lines of the rarefaction fan are thicker).

In the following, we illustrate this analysis of the wave structure by giving the detailed solution for the state and for the sensitivity in a specific case. The initial data for the state on the physical variables is the following: ρL = 1, uL = 0, pL = 1, ρR = 0.125, uR = 0, pR = 0.1.