Pierre Gabriel

Edwige Godlewski

After a roller coaster of almost four years, here I am, writing these acknowledgements, which will be a mixture of different languages and a lot of emotions. D'abord, je tiens à remercier mes deux directeurs de thèse, je n'aurais pas pu demander mieux. Un grand merci à tous les deux pour m'avoir acceptée en thèse à mi-parcours. Régis, merci parce que même à distance tu as été très présent et parce qu'à chaque visite à Sophia ma thèse a énormément avancé. Christophe, merci d'avoir toujours trouvé du temps pour moi, même quand je me présentais à ton bureau sans rendez-vous. Ça a été un plaisir de travailler avec vous deux.

IwanttothankJeff Borggaard for his careful reading of my manuscript and for all his suggestions, english-and maths-wise. Je veux aussi remercier Frédéric Lagoutière pour son rapport très détaillé et l'intérêt porté à mon travail. Merci à Anca Belme, Carole Delenne,

Context and positioning

This PhD thesis deals with sensitivity analysis (SA) for nonlinear hyperbolic partial differential equations (PDEs). These two topics separately are very well known and have been studied thoroughly: however, when considered together, many problems arise and the literature on this subject is far from being comprehensive.

Sensitivity analysis

SA is the study of how changes in the inputs of a model affect the outputs. The sensitivity itself is defined as the derivative of the state, i.e. the solution to the PDE model considered, with respect to a parameter of interest. SA is obviously a valuable tool for engineering applications, since it allows the quantification of changes in the physical response of a system to any change of parameter values: one of its straightforward application in this direction is uncertainty quantification (UQ). SA is an efficient and deterministic method to estimate expectation and standard deviation of the state variables, as an uncertainty propagation technique. Another application of SA is optimization: sensitivities can be used to compute the gradient of a cost functional. Finally, SA methods can likewise be employed to monitor and explore interactively neighbouring solutions for a negligible computational expense and provide an answer to the question "what if...".

There are two main classes of methods to compute the sensitivities: the discretisethen-differentiate approach and the differentiate-then-discretise one. The first approach consists of, as the name says, first discretising the state PDE system and then differentiating it to get a numerical approximation of the sensitivities. In the second approach, one differentiates the state system obtaining in this way the sensitivity system, which can then be discretised. The approximated sensitivities obtained with these approaches are different, because the differentiation step and the discretisation step do not commute in general. A popular strategy that falls into the discretise-then-differentiate category is automatic differentiation [START_REF] Hovland | Automatic differentiation and Navier-Stokes computations[END_REF]: a first advantage of this is the simplification of the code that has to be written for the sensitivity; a second one is the fact that the resulting sensitivity is consistent with the discrete state solution, even if a coarser grid is used for the resolution. However, an important drawback of all the discretise-then-differentiate approaches is the need to differentiate all the computational facilitators, such as, for instance, slope limiters in a finite volume framework or even MPI communication in parallel computing. In [START_REF] Martinelli | On the use of second-order derivative and metamodel-based monte-carlo for uncertainty estimation in aerodynamics[END_REF] for instance, they quantify the error due to automatic differentiation applied to classical slope limiters for the 3D Euler system. On the other hand, the differentiate-then-discretise approach provides a sensitivity system, for which computational facilitators can be used, if needed. Both strategies are valid and they are suitable for different applications. A detailed comparison between the two for optimization problems is performed in [START_REF] Gunzburger | Perspectives in flow control and optimization[END_REF]. In this work, we will mainly focus on the differentiate-then-discretise approach.

An alternative to SA which is worth mentioning, especially for optimization, 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. Note that for the adjoint method, too, one needs to choose whether to use ad i fferentiate-then-discretise or a discretise-then-differentiate approach. The adjoint equation is independent of the input parameters, thus this approach is very efficient for optimization problems involving a large number of design parameters, as opposed to the SA approach, which requires the solution of a different sensitivity system for each parameter. However, if the PDEs considered are time-dependent, the adjoint equation should be solved backwards in time, which could lead to practical difficulties. Moreover, SA allows the computation of the derivative of the whole state and it is not restricted to functionals. Finally, we remark that the sensitivity systems are all independent of each other and therefore can be solved in parallel.

Therefore, in this work we focused our attention on the continuous sensitivity equation (CSE) method [BB97, DPB06, DP06, HEPB04], which allows to compute the derivative of the PDE solution itself, at any location and time, with respect to a single input parameter: this is done by formally differentiating the state system with respect to the parameter of interest, and then exchanging the derivatives with respect to the parameter with those in space and time, obtaining in such a way a new system of PDEs, the sensitivity system. Of course, this can be done for as many parameters as needed. This approach is a differentiate-then-discretise approach and it relies on a forward time integration. However, a certain regularity of the state solution is required, a condition that is not always met in the hyperbolic framework.

Hyperbolic equations

A one-dimensional conservation law is an equation of the following form:

∂ t U + ∂ x F(U)=0,
where ∂ t and ∂ x indicate the partial derivatives with respect to time and space, respectively, U is the conserved variable, or state, and F is the flux. It expresses the fact that, on an arbitrary domain, the time variation of the state variable U are equal to the flux of F through the boundary of the domain considered. These kinds of equations often result from the modelling of continuum physics in cases where the dissipation effects are negligible. A conservation law is said to be strictly hyperbolic if the Jacobian matrix of the flux A(U)= ∂F ∂U is R-diagonalisable, whilst if the matrix A(U) has real eigenvalues but it is not R-diagonalisable the system is weakly hyperbolic, and this distinction will be important in the next chapters.

Hyperbolic equations have been studied for a very long time, both from theoretical and numerical points of view and many books have been written on the subject [Tor09, 1.2. Thesis synopsis GR96]. For the scalar conservation laws, a technique named the method of characteristics can be used to compute an analytical solution to the equation, although it provides often an implicit solution, depending on the initial data. For systems, analytical solutions are known only for certain specific initial conditions, for instance for what it is usually called the Riemann problem: an initial value problem with a piecewise constant data presenting only one discontinuity. These problems play a crucial role in the numerical approximation of hyperbolic equations with a finite volume (FV) approach, which is, as of today, the most suitable approach for these kinds of equations thanks to its conservation property. In particular, the first step of Godunov-type method is the solution of a Riemann problem: this can be done exactly, for the problems for which the analytical solution is known, or using a so-called approximate Riemann solver. Many different approximate Riemann solvers have been developed, some based on a simplification of the structure of the solution, such as for instance the HLL solver, some based on a linearisation of the system, like for example the Roe solver.

Depending on the initial and boundary conditions, hyperbolic PDEs can present discontinuities, such as shock waves or contact discontinuities. This common feature in the hyperbolic framework is the main reason why SA in the hyperbolic equations framework is not straightforward. Indeed, such discontinuities in the solution lead to specific issues regarding SA, because they correspond to the presence of Dirac delta functions in the sensitivity fields: the CSE method briefly introduced in the previous section makes sense only under certain assumptions of regularity of the state solution, which may not be verified in the hyperbolic case. 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, where a modification of the sensitivity system was proposed, to "remove" the spikes from the numerical sensitivity solution. More specifically, a modification of HLL Riemann solver used to evaluate fluxes in a FV method was proposed in the context of Saint-Venant equations. The correction is based on the Rankine-Hugoniot conditions, which govern the state across a shock.

Thesis synopsis

In this work we deal with a hierarchy of models of increasing complexity. We start in Chapter 2 from the simplest example of nonlinear hyperbolic equation, the inviscid Burgers' equation, before focusing on one of the most known and important examples of nonlinear hyperbolic PDE system: the Euler system. The Euler system models the dynamics of a compressible material, which can be a liquid or a gas, and three physical variables are considered to describe the flow: the density of the fluid, its velocity and its pressure. At first, in Chapter 3, we deal with a version of this system, known as the p-system, which is simplified under certain physical and mathematical hypothesis: the fluid is assumed to be in barotropic conditions and Lagrangian coordinates are used; then, in Chapter 4 and Chapter 5 we deal with the complete Euler system: in particular, in Chapter 4 the Sod shock tube is studied, along with the corresponding problem for the sensitivity, and in Chapter 5 a non trivial topography and stationary solutions are considered.

As mentioned in the previous section, if standard techniques of SA, such as the CSE method, are used everywhere, even in presence of discontinuities, they provide a sensitivity which contains Dirac delta functions. Since a Dirac delta function cannot be seized numerically, this leads to a spike which is spread in the neighbourhood of the shock and which deteriorates the solution. Moreover, the spikes can change with the numerical discretisation. The use of such sensitivities is impractical for many of the applications mentioned above. Therefore, along the same lines as what is done 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], we modify the sensitivity equation in order to obtain a system which is valid also in the case of discontinuous solutions. In this work, we adopt a similar, though slightly different, point of view: we suggest to add a correction in the form of a source term to balance out the spikes. The source term is non-zero only in case of shock or contact discontinuity and has an amplitude which is proportional to the jump of the state across the discontinuity considered. The computation of the amplitude of the correction term is carried out in Introduction

Contexte et positionnement

Cette thèse porte sur l'analyse de sensibilité (AS) pour les équations aux dérivées partielles (EDP) hyperboliques non-linéaires. Les deux sujets séparément sont très bien connus et ont fait l'objet de plusieurs études approfondies : cependant, lorsqu'ils sont considérés ensemble, de nombreux problèmes se posent et la littérature sur ce sujet n'est pas importante.

Analyse de sensibilité

L'AS est l'étude de la façon dont les changements dans les entrées d'un modèle affectent la sortie. La sensibilité elle-même est définie comme la dérivée de l'état, c'est-à-dire la solution du modèle d'EDP considéré, par rapport à un paramètre d'intérêt. L'AS est évidemment un outil important pour des applications d'ingénierie, car elle permet de quantifier la réponse physique d'un système aux changements de valeurs des paramètres : une application directe dans ce sens est la quantification d'incertitude (UQ). L'AS est une méthode efficace et déterministe pour estimer l'espérance et l'écart-type des variables d'état, en tant que technique de propagation d'incertitude. Une autre application de l'AS est l'optimisation : les sensibilités peuvent être utilisées pour calculer le gradient d'une fonctionnel coût. Enfin, les méthodes d'AS peuvent également être utilisées pour l'exploration interactive des solutions voisines avec un coût de calcul négligeable et fournir une réponse à la question "que se passe-t-il si ...".

Il y a deux classes principales de méthodes pour calculer les sensibilités : l'approche discrétiser puis différencier et l'approche différencier puis discrétiser. La première approche consiste, comme son nom l'indique, à discrétiser d'abord le système d'EDP d'état, puis à le différencier pour obtenir une approximation numérique des sensibilités. Dans la seconde approche, on différencie le système d'état obtenant ainsi le système de sensibilité, qui peut ensuite être discrétisé. Les sensibilités approchées obtenues avec ces deux approches sont différentes, car l'étape de différenciation et celle de discrétisation ne commutent pas en général. La différenciation automatique [START_REF] Hovland | Automatic differentiation and Navier-Stokes computations[END_REF] est une stratégie couramment employée qui entre dans la catégorie discrétiser puis différencier : un premier avantage est la simplification du code qui doit être écrit pour la sensibilité ; un second est le fait que la sensibilité obtenue est consistante avec la solution discrète de l'état, même si un maillage plus grossier est utilisé pour la résolution. Cependant, un inconvénient important de toutes les approches discrétiser puis différencier est la nécessité de différencier tous les facilitateurs de calcul, comme, par exemple, les limiteurs de pente dans un cadre volumes finis ou même la communication MPI en calcul parallèle. Dans [START_REF] Martinelli | On the use of second-order derivative and metamodel-based monte-carlo for uncertainty estimation in aerodynamics[END_REF] par exemple, ils quantifient l'erreur due à la différentiation automatique appliquée à des limiteurs de pente classiques dans le cadre des équations d'Euler 3D. D'autre part, l'approche différencier puis discrétiser fournit un système de sensibilité pour lequel on peut, si nécessaire, utiliser des facilitateurs de calcul spécifiques. Les deux stratégies sont valables et adaptées à différentes applications. Une comparaison détaillée des deux pour les problèmes d'optimisation se trouve dans [START_REF] Gunzburger | Perspectives in flow control and optimization[END_REF]. Dans ce travail, nous nous concentrerons principalement sur l'approche différencier puis discrétiser.

Une alternative à l'AS qui mérite d'être mentionnée, spécialement pour l'optimisation, est la méthode de l'équation adjointe [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], qui introduit des variables supplémentaires, dites adjointes, afin de calculer la dérivée d'une fonctionnelle par rapport à tous les paramètres d'entrée. On note que pour la méthode adjointe aussi il existe les deux approches : différencier puis discrétiser ou discrétiser puis différencier. L'équation adjointe est indépendante du nombre de paramètres d'entrée, ce qui rend cette approche très efficace pour les problèmes d'optimisation avec beaucoup de paramètres, plutôt que l'approche AS, qui nécessite de résoudre un système de sensibilité différent pour chaque paramètre. Cependant, si les EDP considérées dépendent du temps, l'équation adjointe doit être résolue à rebours en temps, ce qui pourrait poser des difficultés pratiques. De plus, l'AS permet de calculer la dérivée de l'état et n'est pas limitée aux fonctionnelles. Enfin, nous remarquons que les systèmes de sensibilité sont tous indépendants entre eux et peuvent donc être résolus en parallèle.

Par conséquent, dans ce travail nous nous sommes concentrés sur la méthode de léquation de sensibilité continue [BB97, DPB06, DP06, HEPB04], qui permet de calculer la dérivée de la solution de l'EDP considérée, à n'importe quel endroit et temps, par rapport à un paramètre d'entrée : cela se fait en différenciant formellement le système d'état par rapport au paramètre d'intérêt, puis en échangeant les dérivées par rapport au paramètre avec celles en espace et en temps, obtenant ainsi un nouveau système d'EDP, le système de sensibilité. Bien sûr, cela peut être fait pour autant de paramètres que nécessaire. Cette approche est une approche différencier puis discrétiser et est basée sur une intégration avec avance en temps. Cependant, une certaine régularité de la solution d'état est nécessaire, condition qui n'est pas toujours vérifiée dans le cadre hyperbolique.

Équations hyperboliques

Une lois de conservation 1D est une équation de la forme suivante :

∂ t U + ∂ x F(U)=0, où ∂ t et ∂ x désignent les
dérivées partielles par rapport au temps et à l'espace, respectivement, U est la variable conservée, ou état, et F est le flux. Elle exprime le fait que, sur un domaine arbitraire, la variation temporelle de la variable d'état U est égale au flux de F à travers la frontière du domaine considéré. Ces équations proviennent souvent de la modélisation de la physique des milieux continus dans les cas où les effets dissipatifs sont négligeables. Une loi de conservation se dit strictement hyperbolique si la matrice jacobienne du flux A(U)= ∂F ∂U est R-diagonalisable, alors que si la matrice A(U) a des valeurs propres réelles mais qu'elle n'est pas R-diagonalisable, le système est faiblement hyperbolique, et cette distinction sera importante dans la suite. Les équations hyperboliques ont été étudiées depuis très longtemps, à la fois d'un point de vue théorique et numérique et de nombreux livres ont été écrits sur le sujet [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF][START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]. Pour les lois de conservation scalaires, une technique appelée méthode des caractéristiques peut être utilisée pour calculer une solution analytique de l'équation, bien qu'elle ne fournisse souvent qu'une solution implicite, selon les données initiales. Pour les systèmes, les solutions analytiques ne sont connues que pour certaines conditions initiales spécifiques, comme par exemple pour ce que l'on appelle un problème de Riemann : un problème aux valeurs initiales avec une donnée constante par morceaux qui présente une seule discontinuité. Ces problèmes jouent un rôle fondamental dans l'approximation numérique des équations hyperboliques avec une approche aux volumes finis, qui reste, à ce jour, l'approche la plus appropriée pour ces équations grâce à ses propriétés de conservation. En particulier, la première étape de la méthode de Godunov est la résolution d'un problème de Riemann : cela peut être fait exactement, pour les problèmes pour lesquels la solution analytique est connue, ou en utilisant un solveur de Riemann approché. De nombreux solveurs de Riemann approchés ont été développés, certains basés sur une simplification de la structure de la solution, comme par exemple le solveur HLL, d'autres basés sur une linéarisation du système, comme par exemple le solveur de Roe.

Selon les conditions initiales et aux bords, les EDP hyperboliques peuvent présenter des discontinuités, comme des chocs ou des discontinuités de contact. Cette caractéristique, très commune dans le cadre hyperbolique, est la raison principale pour laquelle l'AS dans le cadre des équations hyperboliques n'est pas simple. En effet, des discontinuités dans l'état causent des problématiques spécifiques en AS, car elles correspondent à la présence de distributions de Dirac dans les sensibilités : la méthode de léquation de sensibilité introduite brièvement ci-dessus ne fonctionne que sous certaines hypothèses de régularité de l'état, qui peuvent ne pas être vérifiées dans le cadre hyperbolique. Cette question a été explorée dans [BP02, MP01] d'un point de vue théorique, et plus récemment dans [Gui09, GDC09] d'un point de vue numérique, où une modification du système de sensibilité a été proposée, pour "supprimer" les pics dans la sensibilité numérique. Plus spécifiquement, une modification du solveur de Riemann HLL utilisé pour évaluer les flux dans une méthode volumes finis a été proposée dans le contexte des équations de Saint-Venant. La correction est basée sur les conditions de Rankine-Hugoniot, qui gouvernent l'état à travers un choc.

Présentation des travaux de thèse

Dans cette thèse, nous étudions une hiérarchie de modèles de complexité croissante. Nous commençons dans le chapitre 2 par l'exemple le plus simple d'équation hyperbolique non linéaire, l'équation de Burgers, avant de nous intéresser à l'un des exemples les plus connus et les plus importants des systèmes d'EDP hyperboliques non linéaires : le système d'Euler. Le système d'Euler modélise la dynamique d'un matériau compressible, qui peut être un liquide ou un gaz. Pour décrire le flux, trois variables physiques sont considérées : la densité du fluide, sa vitesse et sa pression. Dans un premier temps, dans le chapitre 3, nous traitons une version de ce système, connue sous le nom de p-système, qui est simplifiée sous certaines hypothèses physiques et mathématiques : le fluide est supposé être dans des conditions barotropes et il est décrit en coordonnées lagrangiennes ; puis, dans le chapitre 4 et le chapitre 5 nous traitons le système d'Euler complet : en particulier, dans le chapitre 4 nous étudions le tube à choc de Sod et le problème correspondant pour la sensibilité, et dans le chapitre 5 nous considérons une topographie non triviale et des solutions stationnaires.

Comme mentionné dans la section précédente, si des techniques standard d'AS, telles que la méthode de léquation de sensibilité continue, sont utilisées tout le temps, même en présence de discontinuités, elles fournissent une sensibilité qui contient des Dirac. Comme un Dirac ne peut pas être évalué numériquement, cela conduit à un pic qui s'étale dans le voisinage du choc et qui donc détériore la solution. De plus, les pics peuvent changer avec la discrétisation numérique. Ces sensibilités sont inadaptées à plusieurs des applications mentionnées ci-dessus. Par conséquent, dans le même ordre d'idées de ce qui est fait dans [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], nous modifions l'équation de sensibilité afin d'obtenir un système valable même dans le cas de solutions discontinues. Dans ce travail, nous adoptons un point de vue similaire, quoique légèrement différent : nous suggérons d'ajouter une correction sous la forme d'un terme source pour équilibrer les pics. Le terme source est non nul uniquement en cas de choc ou de discontinuité de contact et a une amplitude proportionnelle au saut de l'état à travers la discontinuité considérée. Le calcul de l'amplitude de cette correction est effectué en détail dans les chapitres suivants et il est basé sur une intégration des équations sur un volume de contrôle contenant une discontinuité. En ce qui concerne la discrétisation de ce terme source, nous remarquons qu'un détecteur de choc est nécessaire pour que le terme soit cohérent et pour éviter une surcorrection numérique.

En ce qui concerne la discrétisation de la loi de conservation, le problème d'état est bien connu et de nombreux de schémas numériques ont été développés au cours des années, mais ce n'est pas le cas pour les équations de sensibilité. De plus, considérer le système dans son ensemble n'apporte pas d'avantage puisque le système composé par l'état et la sensibilité n'est que faiblement hyperbolique en général. Pour ces raisons, notre stratégie d'un point de vue numérique est la suivante : d'abord, nous résolvons l'état en utilisant un schéma volumes finis classique ; ensuite, nous calculons le terme source, qui dépend de l'état ; enfin, un schéma numérique volumes finis adapté est appliqué au problème de sensibilité. Cette approche fonctionne bien si l'équation d'état est une équation scalaire, mais les résultats numériques montrent que pour les systèmes certaines précautions supplémentaires sont nécessaires : en particulier, la diffusion numérique joue un rôle fondamental et dégrade la convergence vers la solution exacte de la sensibilité. Par conséquent, dans ce travail nous adaptons un schéma numérique de type Godunov anti-diffusif (AD), introduit dans [START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF], et nous proposons des versions AD de différents schémas volumes finis, d'ordre différent en temps et en espace.

Enfin, nous explorons deux applications différentes de l'AS : l'optimisation et la quantification d'incertitude. Dans les deux applications nous considérons plusieurs paramètres d'entrée. Dans le chapitre 4, nous effectuons une analyse de quantification d'incertitude pour un problème de Riemann sur le système Euler complet : différentes approches d'AS sont comparées à la méthode de Monte Carlo plus coûteuse. Dans le chapitre 5 nous nous attaquons à un problème d'optimisation : en particulier, nous traitons un problème inverse pour la pression, où les paramètres d'optimisation considérés modifient directement la topographie. Le but de ces applications est de comprendre l'importance de la diffusion numérique et du terme de correction dans des situations plus réalistes et de définir des directives pratiques pour l'utilisation de l'AS.

Au cours de ma première année de thèse j'ai travaillé sur la modélisation et l'optimisation de stratégies de course, sous la direction d'Amandine Aftalion. Ce travail a abouti à un papier publié [START_REF] Fiorini | Optimization of running strategies according to the physiological parameters for a two-runners model[END_REF], qui est reproduit dans l'annexe A.

| Scalar case

This chapter deals with the nonlinear scalar case. In this simple framework we are able to introduce and tackle some of the main problems arising in the coupling of sensitivity analysis and hyperbolic equations: in particular, the definition of the source term and of a shock detector.

Analytical solution of the state equation

A scalar one-dimensional conservation law can be written in the following form:

( ∂ t u(x, t)+∂ x f (u(x, t)) = 0 x 2 R,t>0 u(x, 0) = g(x) x 2 R, (2.1)
where u : R ⇥ R + ! R is the conserved variable, f : R ! R is the flux function and g : R ! R is the initial condition. The method of characteristics can be used to solve analytically (2.1): the characteristics are curves in the plane (x, t) along which the PDE becomes an ordinary differential equation (ODE). In this case, we look for a set of curves along which the solution u is constant. Let x c (t) be the parametrisation of the curves, then one has:

0= d dt u(x c (t),t)=∂ x u(x c (t),t) dx c dt + ∂ t u(x c (t),t). (2.2) Considering that ∂ x f (u(x, t)) = f 0 (u(x, t))∂ x u(x, t)
, where f 0 = df du , and comparing the equation in (2.1) with (2.2), one can write:

dx c dt = f 0 (u(x c (t),t)) = f 0 (u(x c (0), 0)) ) dx c dt = f 0 (g(x 0 )), (2.3)
where x 0 = x c (0). Therefore, the characteristics are straight lines and their slope is f 0 (g(x 0 )). We can write the solution implicitly as follows:

u(x, t)=g(x -tf 0 (u(x, t))). (2.4)
This method is valid as long as there is no intersection among the characteristics: in case of intersection a shock is generated and the solution u(x, t) becomes discontinuous. Let us now examine whether or not it is possible for two characteristic to intersect. First, we deal with the case where g is a regular function. For this purpose, we consider two characteristics, x c,1 (t) and x c,2 (t),s u c ht h a tx c,i (0) = x i and x 1 <x 2 .I ft h et w o intersect, then one has:

x 1 + f 0 (g(x 1 ))t = x 2 + f 0 (g(x 2 ))t, therefore: t = - x 2 -x 1 f 0 (g(x 2 )) -f 0 (g(x 1 )) ,
and if f 0 • g =: h 2 C 1 for the mean value theorem 9 ξ 2 (x 1 ,x 2 ) such that:

t = - 1 h 0 (ξ)
.

Let us observe that if h is increasing the characteristics do not intersect for any time t>0, therefore there is no shock. Otherwise, we can define the breaking time t s as the smallest t for which the characteristics intersect:

t s = - 1 min h 0 (x) = - 1 min d dx f 0 (g(x))
.

(2.5)

The point x s,0 from which the shock originates is known too:

x s,0 =x + f 0 (g(x))t s , where x := argmin h 0 (x).

If g is discontinuous in a point x d , two scenarios are possible: if h(x + d ) <h(x - d ),t h e initial discontinuity is transported for all t>0 (t s =0); otherwise the initial discontinuity is smoothed out and a rarefaction wave is generated.

Once the shock is generated, it moves along a curve x s (t) such that dxs dt = σ(t), which separates the plane (x, t) into two parts, and in both of them the method of characteristics is valid. We denote with the superscript + (respectively -) the quantities in the right (respectively left) part of the plane, i.e. x>x s (t) (respectively x<x s (t)). The speed of the shock σ(t) can be computed using the Rankine-Hugoniot conditions:

σ(t)= f (u + (x s (t),t)) -f (u -(x s (t),t)) u + (x s (t),t) -u -(x s (t),t) . (2.6)
The position of the shock x s (t) can be then determined by solving the following ODE: 8 < :

dx s dt = σ(t),
x s (t s )=x s,0 .

(

This, along with some numerical methods to solve explicitly (2.4) and (2.7), is all we need to compute a reference solution. Finally, the solution of (2.1) can be written in the following compact form:

u(x, t)=u + (x, t)H(x -x s (t)) + u -(x, t)H(x s (t) -x), (2.8)
where u + and u -are obtained separately from (2.4) and H is the Heaviside function.

Sensitivity equation

Derivation of the sensitivity equation

Sensitivity analysis is the study of how variations in the output of a model can be attributed to different sources of perturbation in the model input. In this case, the model

Sensitivity equation

is a PDE such as (2.1), the input is the vector parameter a on which the initial condition g and the flux function f depend, and the output is the state u. In the following, we will consider a scalar parameter a: this is to simplify the notation, without any loss of generality from a theoretical point of view. Assuming that u is continuous and differentiable on its domain, we define its sensitivity with respect to the parameter of interest a as the derivative of u with respect to a and we use the notation u a :

u a = ∂u ∂a .
(2.9)

We now apply the CSE method and we differentiate the system (2.1) with respect to the parameter a,o b t a i n i n g :

( ∂ a (∂ t u(x, t)) + ∂ a (∂ x f (u(x, t))) = 0 x 2 R,t>0 ∂ a u(x, 0) = ∂ a g(x) x 2 R,
and exchanging the derivatives in space and time with the ones with respect to a one obtains the following equation and initial condition for the sensitivity:

( ∂ t u a + ∂ x (f a (u, u a )) = 0 x 2 R,t>0 u a (x, 0) = g a (x) x 2 R, (2.10)
where g a := ∂ a g, f a (u, u a ): =f 0 (u)u a + ∂ a f (u) and we dropped the time and space dependence in the equation for simplicity. We remark that the exchange of the derivatives can be done only under hypothesis of regularity of U, more precisely the equality of mixed partials is guaranteed only if the second partial derivatives of U are continuous.

Analytical solution of the sensitivity equation

The analytical solution of (2.10) can be computed in the regular zones starting from (2.4), where the definition (2.9) is valid. By differentiating (2.4) with respect to a,o n e finds:

u a (x, t)=g 0 (x -tf 0 (u(x, t)))(-tf 00 (u(x, t))u a )+g a (x -tf 0 (u(x, t))),
therefore one can obtain an explicit expression for the sensitivity:

u a (x, t)= g a (x -tf 0 (u(x, t))) 1+tf 00 (u(x, t))g 0 (x -tf 0 (u(x, t))) = = g a (x -tf 0 (u(x, t))) 1+th 0 (x -tf 0 (u(x, t)))
.

(2.11)

Let us observe that the denominator 1+th 0 (x -tf 0 (u(x, t))) is zero if and only if there is an intersection between two characteristics, therefore only along the shock.

If the state u is discontinuous, (2.11) is still valid on both sides of the shock. Along the shock the state u is not differentiable in the classical sense; however, it admits a Dirac distribution as weak derivative. Therefore, differentiating the compact expression (2.8) one obtains:

u a (x, t)=u + a (x, t)H(x -x s (t)) + u - a (x, t)H(x s (t) -x) +(u --u + )∂ a x s (t)δ(x s (t) -x), (2.12)
where δ is the Dirac delta function.

Correction term

In this section, we aim at proposing a new sensitivity system, whose solution does not exhibit spikes. To do that, we add a source term which should compensate the final term of (2.12), obtaining this new equation:

∂ t u a + ∂ x f a (u, u a )=S(u) x 2 R,t>0. (2.13)
Since we do not want to change the equation if the state u is regular, the source term has the following form:

S(u)=α(t)δ(x s (t) -x),
where α is the amplitude of the correction and it is computed by integrating the equation (2.13) over a control volume (x 1 ,x 2 ) ⇥ (t 1 ,t 2 ), containing a discontinuity which moves at speed σ .W eo b t a i n :

σ(t 2 -t 1 )u - a -σ(t 2 -t 1 )u + a +(t 2 -t 1 )(f a (u + ,u + a ) -f a (u -,u - a )) = Z t 2 t 1 α(t)dt,
and dividing this by (t 2 -t 1 ) and as the control volume goes to zero one has:

α(t)=σ(u - a -u + a )+f a (u + ,u + a ) -f a (u -,u - a ). (2.14)
If we differentiate the Rankine-Hugoniot conditions (2.6) with respect to a we obtain:

σ(u + a -u - a )+∂ a σ(u + -u -)+σ k (∂ x u + -∂ x u -)∂ a x s (t)= = f a (u + ,u + a ) -f a (u -,u - a )+(f 0 (u + )∂ x u + -f 0 (u -)∂ x u -)∂ a x s (t),
where the terms with ∂ a x s (t) come from the fact that the Rankine-Hugoniot are valid only if evaluated in the position of the shock. This terms are very difficult to estimate, however they are all zero in a first order finite volume framework, where the solution u is a piecewise constant, and therefore ∂ x u + =0and ∂ x u -=0 . Neglecting this terms leads to a much simpler formulation:

σ(u + a -u - a )+∂ a σ(u + -u -)=f a (u + ,u + a ) -f a (u -,u - a ).
Comparing the latter to (2.14) we have the definition of the amplitude:

α(t)=∂ a σ(u + -u -).
(2.15)

Global system

It is possible to write (2.1) and (2.10) as a system, by defining the following vectors:

U =  u u a , F(U)=  f (u) f 0 (u)u a + ∂ a f (u) , G =  g g a .
Therefore, the global system in the conservative form is the following:

( ∂ t U + ∂ x F(U)=0 x 2 R,t>0 U(x, 0) = G(x) x 2 R. (2.16) 0 ... i -1 i i +1 ... N -1 i -3 2 i -1 2 i + 1 2 i + 3 2 N -3 2 N -1 2 b c ∆x Figure 2.1 -Spatial discretisation.
We can also write (2.16) in the non conservative form:

( ∂ t U + A(U)∂ x U =0 x 2 R,t>0 U(x, 0) = G(x) x 2 R, (2.17)
where the matrix A is the Jacobian of the system and can be written:

A(U)= ∂F ∂U =  f 0 (u)0 f 00 (u)u a + ∂ a f 0 (u) f 0 (u) . (2.18)
Let us observe that A has two repeated eigenvalues, therefore it is not diagonalisable, unless f 00 (u)u a + ∂ a f 0 (u) ⌘ 0. This means that the global system (2.16) is only weakly hyperbolic. Weakly hyperbolic systems are known to develop Dirac distributions in their solution. We remark that the fact that A is non-diagonalisable is not straightforward to prove in the case of systems.

Numerical methods

The aim of this section is to present the numerical schemes that can be used to solve (2.16). From now on, we will consider the equations on an interval (b, c), which is divided into N cells, all of the same length ∆x. The subscript i indicates the i-th cell, while i ± 1 2 refers to the interfaces (see Figure 2.1). Concerning the time discretisation, we will use a variable time step ∆t, chosen such that the CFL-condition is respected. The superscript n indicates the n-th time step.

As already mentioned, we treat the state and the sensitivity system separately. For the state there is no particular issue: all the classical numerical schemes can be used. However, for the sensitivity a modified scheme is used. The suggested strategy for solving (2.16) consists of the following steps:

(i) solution of the state using numerical methods for hyperbolic equations;

(ii) computation of the correction term;

(iii) solution of the sensitivity with a corrected numerical scheme.

Numerical schemes for the state

The method chosen for the state equation is the approximate Riemann solver of Roe: the reason of this choice is that it provides the exact solution in the case of an isolated shock. For the sake of simplicity we consider the case f 0 (u) > 0. The scheme in the internal cells (i.e. i 6 =0) is the following:

x t σ = f (u + )-f (u -) u + -u - u a,R = @u R @a u a,L = @u L @a
u n+1 i = u n i + ∆t ∆x (f (u n i-1 ) -f (u n i )), (2.19)
while for the first cell we impose u n 0 = g(b) 8n.

Numerical schemes for the sensitivity

For the sensitivity equation (2.10), i.e. the one without the source term, the same scheme used for the state can be applied:

u n+1 a,i = u n a,i + ∆t ∆x (f a (u n i-1 ,u n a,i-1 ) -f a (u n i ,u n a,i )),
however, this scheme provides sensitivity with spikes. We here propose a second scheme, which comes from a quite natural discretisation of the source term:

u n+1 a,i = u n a,i + ∆t ∆x (f a (u n i-1 ,u n a,i-1 ) -f a (u n i ,u n a,i )) + σ n a,i (u n i -u n i-1 )δ n i , (2.20)
where δ n i =1if there is a shock in the i-th cell at the n-th time step, δ n i =0otherwise. One of the main difficulties of this problem is to find an efficient shock detector.

We remark that the same scheme (2.20) can be found if we consider that the solution of equation (2.13) is the derivative of the state in the regular zones without any spike along the shock. Therefore it is easy to define a Roe Riemann solver for the sensitivity (cf. Figure 2.

2).

Once the solver is defined, one can compute the average on each cell using a general Godunov method (cf. Figure 2.3), which results in the following scheme: which can be rewritten as follows:

u n+1 a,i = 1 ∆x  ∆t f (u n i ) -f (u n i-1 ) u n i -u n i-1 u n a,i-1 + ✓ ∆x -∆t f (u n i ) -f (u n i-1 ) u n i -u n i-1 ◆ u n a,i = =u n a,i + ∆t ∆x f (u n i ) -f (u n i-1 ) u n i -u n i-1 (u n a,i-1 -u n a,i ), n n +1 i i -1 σ∆t ∆x -σ∆t
u n+1 a,i = u n a,i + ∆t ∆x (f a (u n i-1 ,u n a,i-1 ) -f a (u n i ,u n a,i )) + σ n a,i (u n i -u n i-1 ), (2.21)
and one can find the correction term by difference:

σ n a,i (u n i -u n i-1 )=(f a (u n i ,u n a,i ) -f a (u n i-1 ,u n a,i-1 )) - f (u n i ) -f (u n i-1 ) u n i -u n i-1 (u n a,i -u n a,i-1 ),
and using once again the Rankine-Hugoniot conditions we have the identity:

σ n a,i (u n i -u n i-1 )=f a (u n i ,u n a,i ) -f a (u n i-1 ,u n a,i-1 ) -σ n i (u n a,i -u n a,i-1 ).

Numerical results

In this section, we show the numerical results obtained with the schemes presented in section 2.5 and we compare them with the reference solution computed as in sections 2.1 and 2.2. For this purpose, we consider the Burgers' equation, therefore the flux function is the following:

f (u)= u 2 2 ,
with different initial data g(x). We remark that f 0 (u)=u, it follows that h = g.F o r the sensitivity flux, we have:

f a (u, u a )=uu a .

Riemann problem

We start with the Riemann problem, i.e.:

g(x)= ( u L x<x c , u R x>x c , (2.22)
x s (t) 0. where x c 2 (b, c). Let us observe that if g is increasing (i.e. u R >u L ), the initial discontinuity is smoothed out, otherwise it is transported at a speed which is known explicitly in this simple case:

σ = u R + u L 2 .
(2.23)

Being σ constant, the shock propagates along a straight line of equation:

x s (t)=σt + x c . (2.24)
Since we are interested in the case of shocks, we consider u L >u R , and in particular u L =1and u R =0.1. The domain is the interval (0, 1),a n dx c =0.5. The parameter of interest a is u L , therefore the initial condition for the sensitivity g u L is the following:

g u L (x)= ( 1 x<0.5, 0 x>0.5.
(2.25)

The characteristics and the shock curve for this problem are plotted in Figure 2.4. The analytical solution for the state is the following:

u(x, t)= ( u L x<σt + x c , u R x>σt + x c . (2.26)
Differentiating (2.26) with respect to u L one obtains the analytical solution for the sensitivity, which is:

u u L (x, t)= ( 1 x<σt + x c , 0 x>σt + x c .
(2.27)

In Figure 2.6 we show the numerical solution for the state obtained with ∆x = 10 -3 compared to the analytical one, for different times: the difference between the two is due to numerical diffusion, and it gets smaller as ∆x is reduced. Knowing the analytical solution, which will be denoted by the subscript ex in the following, one can compute the error as follows: where u ex (x i ,t n ) is the exact solution in the center of the i-th cell at the n-th time step. Figure 2.8 shows the sensitivity solution without the correction term (i.e. δ i =08i): as one can see, where the state is discontinuous the numerical solution of the sensitivity has a peak, which gets bigger with time, as expected, since the coefficient of the Dirac delta function in (2.12) can be easily computed in this case and it is linear in t.

err(t n )=ku ex (x i ,t n ) -u n i k L 1 (0,1) , (2.28)
In order to apply the corrected scheme (2.20), we need to define a shock detector. In this simple case, it is easy to define a good one:

δ n i = ( 0 if u n i-1 = u n i , 1 otherwise.
(2.29) Figure 2.10 shows in blue the cells in which a shock is detected using (2.29) and in red the real position of the shock: the region is quite large due to the numerical diffusion.

The results obtained with the correction are shown in Figure 2.7. Finally, in Figure 2.9 we show the L 1 norm of the error at final time T =0.2 for the state and for the sensitivity. For reference, ku(x, T )k L 1 (0,1) =0.649 and ku u L (x, T )k L 1 (0,1) = 0.61. We remark that for this purpose we compared the numerical sensitivity obtained with the correction to the analytical sensitivity in (2.12) but without the Dirac term.

Continuous initial condition

Now we present the case of a continuous initial condition, which is plotted in Figure 2.5:

g(x)= ( A sin 2 ( ⇡ L (x -x c )+ ⇡ 2 ) if x c -L 2 <x<x c + L 2 0 otherwise. (2.30)
Let us observe that, for x>x c , g is decreasing, therefore a shock will occur. The breaking time can be computed using (2.5) and using the fact that for the Burger's equation h = g:

t s = - 1 min g 0 (x) = - 1 min{-A⇡ L sin( 2⇡ L (x -x c ))} = L Aπ .
The point from which the shock originates is:

x s,0 = x c + L 4 ✓ 1+ 2 π ◆ .
In this case, the curve x s (t) it is not known analytically because the values u + and u - are not (we recall that the solution for the state is known only in the implicit form (2.4)). However, it is possible to compute numerically x s (t) and to plot the exact solution u(x, t) as follows. One can define the following matrix:

X[n, i]=x i + t n g(x i ).
To compute a numerical approximation of x s (t n ), the strategy is the following:

(a) let n be the closest time step to the breaking time t s ;

(b) find k such that X[n, k]  x s,0 and X[n, k + 1] ≥ x s,0 ; (c) u + = g(xk +1 ) and u -= g(xk) ) σ = u + +u - 2 ; (d) x s (t n+1 )=x s (t n )+σ(t n+1 -t n );
(e) erase all the characteristics that cross the shock between time t n and t n+1 In Figure 2.11 we show the characteristics of the Burger's equation with initial condition (2.30) and the shock curve computed as explained above.

It is now possible to plot u for all time steps t n in the points X[n, i] (u(X[n, i],t n )= g(x i )), which is exact and is a good reference solution if the initial mesh is fine enough.

The same can be done for the sensitivity:

u a (X[n, i],t n )= g a (X[n, i] -t n g(x i )) 1+t n g 0 (X[n, i] -t n g(x i ))
.

In Figure 2.12 we show the numerical and the reference solution for the state, obtained with ∆x = 10 -3 .

In this case a more complex shock detector, based on the second derivative, has been used. Let D be the maximum of the second derivative of the initial condition, then δ i is defined as follows:

δ n i = 8 < : 1 if |u n i-1 -2u n i + u n i+1 | ∆x 2 >kD, 0 otherwise,
where k is a proper constant. In Figure 2.17 we show the results of this shock detector with k = 30. In Figures 2.16 and 2.15 we show the solution of the sensitivity with respect to the parameter x c , with and without the correction. The figures in the middle correspond to t = t s : as one can see, before the breaking time the solutions are the same, whilst after the breaking time the solution given by the corrected scheme does not present any peak.

Finally, in Figures 2.14 and 2.13 we show the solution of the sensitivity with respect to the parameter A, with and without the correction: in this case, the over-correction due to the fact that a shock is detected even before the breaking time, is slightly more evident.

| The p-system

In this chapter, which is an adaptation of two works [START_REF] Chalons | Sensitivity analysis for the Euler equations in Lagrangian coordinates[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], we deal with the Euler system in barotropic conditions and in Lagrangian coordinates.

Problem description

As already mentioned in the previous chapters, 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 2 R,t > 0, U(x, 0) = U 0 (x).
In fact, it is well known systems of this type, as well as the scalar conservation laws introduced in Chapter 2, can have discontinuous solutions, regardless of the regularity of the initial condition U 0 .I ft h es t a t eU is discontinuous, the sensitivity U a = ∂ a U will exhibit Dirac delta functions.

In this Chapter, we consider the barotropic Euler equations in Lagrangian coordinates, i.e. the p-system; however, everything can be extended to any hyperbolic system: for instance, the complete Euler system and the quasi 1D Euler system will be the subject of the following chapters. The choice of starting with 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 is written:

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

M(τ,u)=  0 -1 p 0 (τ )0 (3.2)
and its eigenvalues are real and distinct λ ± = ±c, where c = p -p 0 (τ ) is the Lagrangian sound speed. Therefore M is R-diagonalisable, and (3.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 (3.1), we can apply the CSE method, differentiate (3.1) with respect to a and obtain the following sensitivity equations:

( ∂ t τ a -∂ x u a =0, ∂ t u a + ∂ x (p 0 (τ )τ a )=0. (3.3)
We remark that if the pressure law depends directly on the parameter a (i.e. p = p(τ ; a)), there is an additional term in the second equation: ∂ x p a (τ ; a), where ∂ a p(τ ; a)=p a (τ ; a).

In the following, we will not consider this case. 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 0 (τ )τ a ,
and rewrite the systems (3.1) and (3.3) in a vectorial form:

( ∂ t U + ∂ x F(U)=0, ∂ t U a + ∂ x F a (U, U a )=0. (3.4)
The Jacobian matrix of the global system (3.4) is calculated as:

A(V)= ∂G ∂V = 2 6 6 4 0 -10 0 p 0 (τ )0 0 0 00 0 -1 p 00 (τ )τ a 0 p 0 (τ )0 3 7 7 5 , 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 (3.4) has two repeated eigenvalues, the same λ ± as the original system (3.1), and that the matrix A(V) is not R-diagonalisable as soon as τ a 6 =0 . This means that the global system (3.4) is only weakly hyperbolic. Therefore, the system (3.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: [BJL03, CKM12, FLF92, Jos93, LeF90, YZ12]. However, sensitivities with Dirac delta functions are unusable for many applications. 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.

Source term

In this section, we aim at proposing a new version of (3.4) which is also valid for discontinuous solutions of the state variable U. Recall indeed that (3.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).

and we suggest to differentiate them with respect to the parameter a.A sw ed ot h a t , 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 + )σ k,a + σ k (rU + -rU -)∂ a x k,s (t)= = F - a -F + a + ✓ ∂F(U + ) ∂U rU + - ∂F(U -) ∂U rU - ◆ ∂ 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 , (3.9)
with σ k,a = ∂ a σ k . Comparing the latter with (3.8), one is thus led to set

ρ ρ ρ k (t)=σ k,a (U + -U -). (3.10)
Our choice is of course valid for each k-discontinuity of the state solution, leading us to definition (3.5) 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 + )a n dσ 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 (3.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 (3.1), which is wellknown but necessary to solve (3.4). Then we differentiate 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 -p 0 (τ ) and λ 2 (U)= p -p 0 (τ ) ,
and the eigenvectors:

r 1 (U)= 2 6 6 6 4
2 p -p 0 (τ ) p 00 (τ ) - 2p 0 (τ ) p 00 (τ ) 3 7 7 7 5 and r 2 (U)= 2 6 6 6 4 - 2 p -p 0 (τ ) p 00 (τ ) - 2p 0 (τ ) p 00 (τ ) 3 7 7 7 5 ,
which are chosen in such a way that rλ i • r i =1 . Since the pairs (λ 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 presented in Figure 3.4 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 3.4 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 ). (3.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 (τ ) < 0, therefore:

u = u L ± p -(p(τ ) -p(τ L ))(τ -τ L ) . (3.12)
In order to chose the sign in (3.12), we use the Lax conditions:

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

- p -p 0 (τ ) < - p -p 0 (τ L ) ) τ < τ L ,
where we used the hypothesis p 00 (τ ) > 0. Both σ and (τ -τ L ) are negative, therefore their product is positive and from the first equation of (3.11) we can say that u L >u . We can conclude that the sign in (3.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(τ ) -p(τ L ))(τ -τ L ) . (3.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(τ ) -p(τ R ))(τ -τ R ) . (3.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 ξ 7 ! e U 1 (ξ) in a 1-rarefaction wave starting from U L respects the following equation:

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

8 > > > < > > > : τ 0 1 (ξ)= 2 p -p 0 (τ 1 ) p 00 (τ 1 ) , τ1 (ξ 0 )=τ L , ũ0 1 (ξ)=- 2p 0 (τ 1 ) p 00 (τ 1 ) , ũ1 (ξ 0 )=u L . (3.16)
The first equation can be rewritten as

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

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

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

ũ1 = u L + Z ⌧1 ⌧ L p -p 0 (τ ) dτ. If p(τ )=τ -γ ,o n ew o u l dfi n d : τ1 = ✓ γt 2 (x -x c ) 2 ◆ 1 γ+1 , ũ1 = u L + 2 p γ 1 -γ ✓ τ 1-γ 2 1 -τ 1-γ 2 L ◆ . (3.17)
Repeating exactly the same procedure, one can find the points reachable through a 2-rarefaction starting from U R :

ũ2 = u R + Z ⌧ R ⌧2 p -p 0 (τ ) dτ,
hence for our choice of p: Finally, in order to compute U ⇤ we define the two following 1-and 2-wave curves:

ũ2 = u R + 2 p γ 1 -γ ✓ τ 1-γ 2 R - τ 1-γ 2 2 ◆ , τ2 = ✓ γt 2 (x -x c ) 2
g 1 (τ ; U L )= 8 < : u L - q -(τ -γ -τ -γ L )(τ -τ L ) if τ  τ L , u L + 2 p γ 1-γ (τ 1-γ 2 -τ 1-γ 2 L ) if τ > τ L , (3.19)
g 2 (τ ; U R )= 8 < : u R + q -(τ -γ -τ -γ R )(τ -τ R ) if τ  τ R , u R + 2 p γ 1-γ (τ 1-γ 2 R -τ 1-γ 2 ) if τ > τ R , (3.20)
which are smooth functions whose derivatives with respect to τ will be denoted g 0 i . The intermediate state τ ⇤ is defined as the intersection between g 1 and g 2 ,a n do n eh a s :

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

u L -u R < 2 p γ 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: 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 e U and this concludes the computation of the analytical sensitivity. To compute τ ⇤ a and u ⇤ a ,w ed i fferentiate, with respect to a, the following equality:

U a,L = ∂U L ∂a U a,R = ∂U R ∂a .
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 0 1 (τ ⇤ ; U L )τ ⇤ a + ∂g 1 ∂τ L (τ ⇤ ; U L )τ a,L + ∂g 1 ∂u L (τ ⇤ ; U L )u a,L = = g 0 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 = @g 2 @⌧ R (τ ⇤ ; U R )τ a,R + @g 2 @u R (τ ⇤ ; U R )u a,R -@g 1 @⌧ L (τ ⇤ ; U L )τ a,L -@g 1 @u L (τ ⇤ ; U L )u a,L g 0 1 (τ ⇤ ; U L ) -g 0 2 (τ ⇤ ; U R )
.

Finally, we differentiate the state solution in the rarefaction e U given by (3.17)-(3.18).

Examples

In the numerical section of this work, we will consider as a 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 (3.17) and (3.18) one can observe that τi does not depend on those parameters, therefore: Concerning the sensitivity of the velocity, one obtains:

τa,i = ∂ τi ∂a =0 i =1, 2. x t x c U L U ⇤ U R (a) 1-shock-2-shock. x t x c U L U ⇤ U R e U (b) 1-shock-2-rarefaction x t x c U L U ⇤ U R e U (c) 1-rarefaction-2-shock. x t x c U L U ⇤ U R e U 1 e U 2 (d) 1-rarefaction-2-rarefaction.
ũa,1 = ∂ ũ1 ∂a = u a,L - p γ τ -1+γ 2 L τ a,L , ũa,2 = ∂ ũ2 ∂a = u a,R + p γ τ -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.5. This simplification is due to the fact that we are considering a reduced Euler system, under barotropic conditions (cf.

[GLC07]). 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 3.6).

(ii) if the parameter of interest is u L (or u R )w eh a v eũ 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 3.7).

Classical numerical schemes

The aim of this section is to design relevant numerical schemes for (3.6). 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 will we get a perfect agreement between exact state. We obtain the following formulas (cf. Figure 3.5 for the notations):

κ 1,j-1/2 = 8 < : u * j-1/2 -u j-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 = 8 < : u * j-1/2 -u j ⌧ 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 becomes:

U n+1 a,j = U n a,j + ∆t ∆x ⇣ κ 2,j-1/2 (U ⇤ a,j-1/2 -e U R a,j-1/2 )+c 2,j-1/2 ( e U R a,j-1/2 -U n a,j ) -κ 1,j+1/2 (U ⇤ a,j+1/2 -e U L a,j+1/2 ) -c 1,j-1/2 ( e U L a,j+1/2 -U n a,j ) ⌘ , (3.22)
where the intermediate states

U ⇤ a,j-1/2 , U ⇤ a,j+1/2 , e U R a,j-1/2
and e U L a,j-1/2 are known analytically, from section 3.3.2. We remark that the source term is encompassed in (3.22), since (3.22) comes from the exact Riemann solver of (3.6).

AR o e -t y p em e t h o d

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 = - s - p(τ n j-1 ) -p(τ n j ) τ n j-1 -τ n j , λ ROE R,j-1/2 = s - p(τ n j-1 ) -p(τ n j ) τ n j-1 -τ n j if τ n j-1 6 = τ n j and ⌥ q -p 0 (τ 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 are given as:

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 . (3.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 )). (3.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 ◆ , (3.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 (3.5). In this work, we use a very simple shock detector: in section 3.3.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
⌘ .

(3.26)

We remark that the discretisation of the source term (3.25) is such that U ⇤ a,j-1/2 = @ a U ⇤ j-1/2 , in fact differentiating (3.23) with respect to a,o n efi n d s :

@ 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 (3.23) 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 (3.26), once the shock detectors are added. Furthermore, the definition (3.26) encompasses the source term, which means that we can use the update formula (3.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 ), (3.27)
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 twostep 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 (3.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. [Bou04], 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 3.8 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 is written as:

( 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
).

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 3.8 -Second order discretisation. In red, the corresponding first order discretisation.

Numerical results

We present some numerical results obtained with the schemes described in the previous section, on a spatial domain is (0, 1), with final time 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 3.9 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 3.9-3.10. 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 3.11 we plot the convergence curves of the all the schemes for each variable: all the methods converge as expected for the state variable; however, for the sensitivity the error seems to be convergent for coarser meshes, but 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 presented here 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 (3.20). 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 0 2 (⌧ L ; U R ) ◆ ' ✓ 1 -9.35 ◆ , U a,R = ✓ 0 0 ◆ .
Figure 3.12 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 averaging operation. For this reason, in the next section we introduce a scheme without numerical diffusion in the shock.

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 that does not have any numerical diffusion in the shock. The scheme was first introduced in [CG08] 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:

x n j-1/2 = x j-1/2 + σ n j-1/2 ∆t n ,
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 at time t n+1 , U n+1 j , is chosen randomly among

U n+1 j-1 , U n+1 j
,a n dU n+1 j+1 , in agreement with their rate of presence in the cell. More 3.6. An anti-diffusive Roe-type numerical scheme precisely, given a random sequence (↵ n ) varying in (0, 1), the choice is the following:

t x x j-1/2 x j+1/2 x j-1/2 U j-1 U ⇤ j-1/2 U j e U 1 x j+1/2 U ⇤ j+1/2 U j+1 e U 1 e U 2
U n+1 j = 8 > > > < > > > : U n+1 j-1 if ↵ n+1 2 ⇣ 0, ∆t ∆x max(σ n j-1/2 , 0) ⌘ , U n+1 j if ↵ n+1 2 h ∆t ∆x max(σ n j-1/2 , 0), 1+ ∆t ∆x min(σ n j+1/2 , 0) ⌘ , U n+1 j+1 if ↵ n+1 2 h 1+ ∆t ∆x min(σ n j+1/2 , 0), 1
⌘ .

(3.28)

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]):

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

Our choice for σ n j+1/2 is the following:

σ n j+1/2 = 8 > < > : λ n j+1/2
u j >u j+1 and ⌧ j < ⌧ j+1 , -λ n j+1/2 u j >u j+1 and ⌧ j > ⌧ j+1 , 0 otherwise.

(3.29) If u is increasing, which means that a rarefaction is expected, the mesh is not modified, whilst for the case of an 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 3.13).

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

( @ t U + @ x F(U) -σ@ x U =0, @ t U + σ@ x U =0. (3.30)
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 (3.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:

e R I j (U n )=R I j (U n ) -σ n j-1/2 U ⇤ j-1/2 + σ n j+1/2 U ⇤ j+1/2
, where in the definition of σ j-1/2 we use λ j-1/2 = λ ROE j-1/2 . Then, the scheme is written:

( U n+1 j = ∆x ∆x j U n j + ∆t ∆x j e R I j (U n ), +(3.28).
(3.31)

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:

e R II j (U n )=R II j (U n ) -σ n j-1/2 U ⇤ j-1/2 + σ n j+1/2 U ⇤ j+1/2 .
Then the scheme writes:

8 > > < > > : U n+1/2 j = U n j + ∆x-∆x j ∆x j U n j + ∆t 2∆x j e R II j (U n ), U n+1 j = U n j + ∆x-∆x j ∆x j U n+1/2 j + ∆t ∆x j e R II j (U n+1/2 ), +(3.28). (3.32)
From the two-step problem (3.30) point of view, the discretisation (3.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 are shown in Figures 3.14-3.15. As one can see, removing the numerical diffusion in the shock for the state variables allows us to be more precise in the definition of the source term which, in turns, provides us with better solution for the sensitivity: the plateau in the star zone is correct. Furthermore, we show in Figure 3.16 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. For reference, k⌧ (x,

T )k L 1 (0,1) ' 0.4502, ku(x, T )k L 1 (0,1) ' 0.2357, k⌧ u L (x, T )k L 1 (0,1) ' 0.0242,a n dku u L (x, T )k L 1 (0,1) ' 0.5158.
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 3.4 and, since a = ⌧ R , the rarefaction wave splits into two discontinuities for the sensitivity as shown in Figure 3 (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 3.17 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.

| The Euler system

Introduction

In this chapter we deal with the complete Euler system. In the last part of the chapter an uncertainty quantification problem is defined.

The state system

The Euler system is written: 8 < :

@ t ⇢ + @ x (⇢u)=0, @ t (⇢u)+@ x (⇢u 2 + p)=0, @ t (⇢E)+@ x (u(⇢E + p)) = 0, (4.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 ◆ , (4.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 = q (γ -1)(H -1 2 u 2 ) . We can rewrite the system (4.1) in the vectorial form:

@ t U + @ x F(U)=0, (4.3)
where

U = 2 4 ⇢ ⇢u ⇢E 3 5 = 2 4 w 1 w 2 w 3 3 5 , F(U)= 2 4 ⇢u ⇢u 2 + p u(⇢E + p) 3 5 = 2 6 6 4 w 2 w 2 2 w 1 +(γ -1) ⇣ w 3 -1 2 w 2 2 w 1 ⌘ γ w 2 w 3 w 1 -(γ-1) 2 w 3 2 w 2 1 3 7 7 5 .
One can also write (4.1) in the nonconservative form:

@ t U + A(U)@ x U =0, (4.4)
where the Jacobian matrix A is:

A(U)= @F @U = 2 6 4 01 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 3 7 5 ,
its eigenvalues are λ 1 = u -c, λ 2 = u,a n dλ 3 = u + c and its eigenvectors are:

r 1 = 2 4 1 u -c H -uc 3 5 , r 2 = 2 4 1 u u 2 2 3 5 , r 3 = 2 4 1 u + c H + uc 3 5 .
Therefore A is R-diagonalisable and the system (4.1) is strictly hyperbolic. At last, (4.3) will be supplemented with a given initial data U(x, t = 0) = U 0 (x), 8x 2 R.

The sensitivity system

Considering only smooth solutions of (4.1), one can apply the Continuous Sensitivity Equation (CSE) [HEPB04, BB97, DPB06] method which consists in differentiating (4.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:

8 < : @ t ⇢ a + @ x (⇢u) a =0, @ t (⇢u) a + @ x (⇢ a u 2 +2⇢uu a + p a )=0, @ t (⇢E) a + @ x (u a (⇢E + p)+u((⇢E) a + p a )) = 0, (4.5)
which can be written in vectorial form as

@ t U a + @ x F a (U, U a )=0, (4.6)
where we used the following notation:

U a = @ a U = 2 4 ⇢ a (⇢u) a (⇢E) a 3 5 , F a (U, U a )=@ a F(U)= 2 4
(⇢u) a ⇢ a u 2 +2⇢uu a + p a u a (⇢E + p)+u((⇢E) a + p a ) 3 5 .

Note that differentiating (4.2) one has:

p a =(γ -1)((⇢E) a - 1 2 ⇢ a u 2 -⇢ a uu a )
which acts as a closure relation for (4.5). The initial data for the sensitivity is nothing but U a (x, t = 0) = @ a U 0 (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 U a = 2
G(V)=  F(U) F a (U, U a ) = 2
w 1 +(γ -1) ⇣ w 3 -1 2 w 2 2 w 1 ⌘ γ w 2 w 3 w 1 -(γ-1) 2 w 3 2 w 2 1
) .

w 5 γ-3 2
w 2 2 w 4 w 2 1 -(γ -3) w 2 w 5 w 1 +(γ -1)w 6 γ 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)
Therefore, the complete system writes:

( @ t V + @ x G(V)=0, V(x, 0) = V 0 (x), (4.7)
with V 0 (x)=( U 0 (x), @ a U 0 (x)) t . The Jacobian matrix of the complete system has the following form:

@G(V) @V = M(V)=  A0 BA
where A is the Jacobian matrix of the state system and B writes:

B = 2 4 00 0 (γ -3)uu a (3 -γ)u a 0 (?)( •) γu a 3 5 , (4.8)
with

(?)=- c 2 γ -1 p a p u + 3 2 (γ -2)u 2 u a + c 2 γ -1 u⇢ a ⇢ - c 2 γ -1 u a + γ u 3 ⇢ a ⇢ , and
(•)= γ 2 u 2 ⇢ a - c 2 γ -1 ⇢ a + 6 -5γ 2
u 2 ⇢ a ⇢ +(3-2γ)uu a + 3(γ -1) u⇢ a ⇢ 2 + c 2 γ -1 p a 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 (4.7) is weakly hyperbolic.

Proof. A system of the form (4.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:

p M (x)=(x -λ 1 ) 2 (x -λ 2 ) 2 (x -λ 3 ) 2 , (4.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 -λ 1 I 6 )(M -λ 2 I 6 )(M -λ 3 I 6 )=0. (4.10)
Let us write (4.10) by blocks:

 A -λ 1 I 3 0 BA -λ 1 I 3  A -λ 2 I 3 0 BA -λ 2 I 3  A -λ 3 I 3 0 BA -λ 3 I 3 =0 (4.11)
Developing the left-hand side products one obtains the following matrix:

 (A -λ 1 I 3 )(A -λ 2 I 3 )(A -λ 3 I 3 )0 (⌅)( A -λ 1 I 3 )(A -λ 2 I 3 )(A -λ 3 I 3 ) , (4.12)
where

(⌅)=B(A -λ 2 I 3 )(A -λ 3 I 3 )+(A -λ 1 I 3 )B(A -λ 3 I 3 )+(A -λ 1 I 3 )(A -λ 2 I 3 )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] 2 4 0 (3 -γ)u 2 u a -(γ-3) 2 2 u 2 u a } 3 5 + +[c -u, 1, 0] 2 4 (γ -3)uu a (γ -2)(3 -γ)u 2 u a +(γ -1)(?) M 3 5 = = - 3 2 (γ -1)u 2 u a +(γ -3)cuu a -c 2 u p a p + c 2 u ⇢ a ⇢ -c 2 u a + γ(γ -1)u 3 ⇢ a ⇢ ,
where there is no need to specify } and M. There is no reason why the quantity above 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 (4.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]. Therefore, we add a correction term in the Rankine-Hugoniot conditions of (4.5), as done in the previous chapters in order to make this system valid also in the present framework of hyperbolic systems with possibly discontinuous solutions:

S = Ns X k=1 δ k ⇢ ⇢ ⇢ k , (4.13)
where N s is the number of discontinuities, which can be either shocks or contact discontinuities, δ k = δ(x -x k,s ) is the Dirac delta function with x k,s position of the k-th shock and ⇢ ⇢ ⇢ k is the amplitude of the k-th correction which is

⇢ ⇢ ⇢ k (t)=σ k,a (U + -U -). (4.14) x t x c U L U ⇤ L U ⇤ R U R b U Figure 4
.1 -Structure of the solution for the Riemann problem for the state. We consider as parameter of interest a = p L , therefore the initial data for the sensitivity is:

x t x c U a,L U ⇤ a,L U ⇤ a,R U a,R b U a
⇢ a,L = ⇢ a,R = u a,L = u a,R = p a,R =0,p a,L =1.
This choice of initial data leads to the structure in Figures 4.1 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:

p L = p ⇤ R 0 B B @ 1 - (γ -1) c R c L ( p * R p R -1) r 2γ ⇣ 2γ +(γ + 1)( p * R p R -1) ⌘ 1 C C A -2γ γ-1 , (4.17)
where c `= r γp ⇢`, with `= L, R. In the star regions, we have:

p ⇤ L = p ⇤ R = p ⇤ , u ⇤ L = u ⇤ R = u ⇤ = c R ✓ p ⇤ p R -1 ◆ s 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:

⇢ ⇤ R = ⇢ R p ⇤ p R 1+ γ-1 γ+1 p R p * 1+ γ-1 γ+1 p * p R ! , ⇢ ⇤ L = ⇢ L ✓ p ⇤ p L ◆ 1 γ .
In the rarefaction wave, we have:

û(x, t)= 2(u ⇤ -u L ) (γ + 1)u ⇤ ✓ x -x c t ◆ +2 c L u ⇤ -u L ⇣ c L -γ+1 2 u ⇤ ⌘ (γ + 1)u ⇤ , ⇢(x, t)=⇢ L ✓ 1 -(γ -1) û(x, t) 2c L ◆ 2 γ-1 , p(x, t)=p L ✓ 1 -(γ -1) û(x, t) 2c L ◆ 2γ γ-1 .
Finally, the solution is given as:

U(x, t)= 8 > > > > > > > > < > > > > > > > > : U L x -x c < -c L t, b U -c L t<x-x c < ⇣ γ+1 2 u ⇤ -c L ⌘ t, U ⇤ L ⇣ γ+1 2 u ⇤ -c L ⌘ t<x-x c <u ⇤ t, U ⇤ R u ⇤ t<x-x c <c R q γ-1 2γ + γ+1 2γ p * p R t, U R x -x c >c R q γ-1 2γ + γ+1 2γ p * p R t. (4.18)
Sensitivity solution:b yd i fferentiating (4.17) 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 ⇤
, where:

Θ =1- (γ -1)c R ⇣ p * p R -1 ⌘ c L r 4γ 2 +2γ(γ -1) ⇣ p * p R -1 ⌘ , Ξ = c R ⇣ p * p R -1 ⌘ c a,R p 2γ c 2 L r 2γ +(γ + 1) ⇣ p * p R -1 ⌘ , Λ = p 2γ c R c L p R r 2γ +(γ + 1) ⇣ p * p R -1 ⌘ , Ψ = γ(γ + 1)c R ⇣ p * p R -1 ⌘ c L p R p 2γ ⇣ 2γ +(γ + 1) ⇣ p * p R -1 ⌘⌘ 3 2 .
In the star regions, by differentiating the corresponding state, one finds:

u ⇤ a = 2c a,L γ -1 1 - ✓ p ⇤ p L ◆ γ-1 2γ ! - c L γ ✓ p ⇤ p L ◆ -γ-1 2γ ✓ p L p ⇤ a -p ⇤ p 2 L ◆ , ⇢ ⇤ a,R = ⇢ R p ⇤ a p R ⇣ 1+ γ-1 γ+1 p R p * ⌘ ⇣ 1+ γ-1 γ+1 p * p R ⌘ + ⇢ R p ⇤ p R γ -1 γ +1 0 @ -p R p * a p * 2 ⇣ 1+ γ-1 γ+1 p * p R ⌘ -p * a p R ⇣ 1+ γ-1 γ+1 p R p * ⌘ (1 + γ-1 γ+1 p * p R ) 2 1 A , ⇢ ⇤ a,L = ⇢ L γ p L p ⇤ a -p ⇤ p 2 L ✓ p ⇤ p L ◆ 1-γ γ .
Finally, in the rarefaction:

ûa (x, t)= 2u L u ⇤ (γ + 1)u ⇤2 x -x c t +2 c a,L u ⇤2 -c a,L u L u ⇤ + c L u L u ⇤ a (γ + 1)u ⇤2 , ⇢a (x, t)=-⇢ L ✓ ûa (x, t)c L -û(x, t)c a,L c 2 L ◆✓ 1 - (γ -1)û(x, t) 2c L ◆ 3-γ γ-1 , pa (x, t)= ✓ 1 - (γ -1)û(x, t) 2c L ◆ 2γ γ-1 -p L γ ✓ ûa (x, t)c L -û(x, t)c a,L c 2 L ◆✓ 1 - (γ -1)û(x, t) 2c L ◆ γ+1 γ-1
.

The sensitivity has the same structure as the state, therefore:

U a (x, t)= 8 > > > > > > > > < > > > > > > > > : U a,L x -x c < -c L t, b U a x-xc t -c L t<x-x c < ⇣ γ+1 2 u ⇤ -c L ⌘ t, U ⇤ a,L ⇣ γ+1 2 u ⇤ -c L ⌘ t<x-x c <u ⇤ t, U ⇤ a,R u ⇤ t<x-x c <c R q γ-1 2γ + γ+1 2γ p * p R t, U a,R x -x c >c R q γ-1 2γ + γ+1 2γ p * p R t. (4.19)
We remark that if one writes the Rankine-Hugoniot conditions across the shock one finds:

-c R r γ -1 2γ + γ +1 2γ p ⇤ p R (U a,R -U ⇤ a,R )+F a (U R , U a,R ) -F a (U ⇤ R , U ⇤ a,R )=S.

Numerical methods

In this section we consider the numerical approximation of (4.16). 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, x j is the center of the j-th cell C j , whose extrema are x j-1/2 and x j+1/2 (cf. Figure 4.3). We use an adaptive time step ∆t n ,c h osen according to a CFL condition, and the intermediate times are t n+1 = t n + ∆t n .W e 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 C j at time t n .

x j-1/2

x j+1/2 ∆x x j C j 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 ,o b t a i n i n gi nt h i sw a yas o l u t i o na t time t n+1 , v(x, t n+1 )=(u(x, t n+1 ), u a (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: different choices for the solution of the Riemann problem lead to different numerical schemes.

Projection step

The projection step is usually performed by averaging the solution v(x, t n+1 ) on the cell:

V n+1 j = 1 ∆x Z x j+1/2 x j-1/2 v(x, t n+1 )dx. (4.20)
However, this projection method introduces numerical diffusion. As shown in Chapter 3, 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 [CG08] 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. In this case, we choose it equal to zero if no shock is expected at the interface j -1/2, whilst it is equal to the speed of the shock, if there is one. More details on the choice of σ n j-1/2 will be provided in the next subsections. 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 Z xj+1/2 xj-1/2 v(x, t n+1 )dx, (4.21)
where ∆x n j =x j+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 = 8 > > > < > > > : V n+1 j-1 if ↵ n+1 2 ⇣ 0, ∆t ∆x max(σ n j-1/2 , 0) ⌘ , V n+1 j if ↵ n+1 2 h ∆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 2 h 1+ ∆t ∆x min(σ n j+1/2 , 0), 1 ⌘ . (4.22)
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 for both the state and the sensitivity (cf. [START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF]):

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

In the next subsections we present some Riemann solvers for the state and for the sensitivity, and their two variations: the diffusive standard version with the projection step (4.20), and their anti-diffusive version (4.21)-(4.22) and we compare their results. The strategy is to solve the state system separately from the sensitivity system, since the global system (4.7) is only weakly hyperbolic and some of the standard approaches require strict hyperbolicity. Moreover, notice that the state system (4.3) evolves in an independent way with respect to the sensitivity system (4.6).

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. It turns out that this property plays an important role in the anti-diffusive approach (4.21)- (4.22). Indeed, it makes the method equivalent to the random choice Glimm scheme in the specific case of an isolated shock wave, and therefore is convergent to the correct solution (still in this very specific case). In addition, we remark that it would not be possible to use a solver with only one intermediate star state, such as HLL, because of the definition of the source term (4.14): two intermediate states are necessary in order to be able to compute the correction term across the contact discontinuity.

The main idea of the Roe scheme is to replace the Jacobian matrix A(U) in (4.4) with a constant matrix A(U L , U R ), 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 = 0 @ 1 ũ - c H -ũc 1 A , r2 = 0 @ 1 ũ ũ2 2 1 A , r3 = 0 @ 1 ũ +c H +ũc 1 A .
The quantities denoted with a tilde are Roe averaged quantities defined as follows:

ũ = p ⇢ L u L + p ⇢ R u R p ⇢ L + p ⇢ R , H = p ⇢ L H L + p ⇢ R H R p ⇢ L + p ⇢ R , c = s (γ -1) ✓ H - 1 2 ũ2 ◆ .
Therefore, the Roe solver consists of four constant states (U L , U ⇤ L , U ⇤ R ,a n dU R , cf. Figure 4.4) connected by three discontinuities travelling at speeds λ ROE i . To compute

λ ROE 1 λ ROE 2 λ ROE 3 x t x c U L U ⇤ L U ⇤ R U R
∆U = U R -U L = 3 X i=1 ↵ i ri . (4.23)
The relation (4.23) is used to compute the coefficients ↵ i , then one has:

U ⇤ L = U L + ↵ 1 r1 = U R -↵ 2 r2 -↵ 3 r3 , U ⇤ R = U R -↵ 3 r3 = U L + ↵ 1 r1 + ↵ 2 r2 . (4.24)
Once all the quantities U ⇤ L , U ⇤ R ,andλ ROE `are known at each interface x j-1/2 , w(x, t n+1 ) is defined and the integrals (4.20) or (4.21) can easily be computed, since w(x, t n+1 ) is a piecewise constant.

It is well known that, in case of transonic rarefaction, the Roe solver provides a nonentropic 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].

Definition of σ σ σ n j-1/2

As already said, σ n j-1/2 is 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.5). The definition of σ n j-1/2 is the following:

σ n j-1/2 = 8 > < > : λ 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 shock detectors, d `,j-1/2 =1if 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 u j <u j-1 , 0 otherwise, d 3,j-1/2 = ( 1 if ⇢ j < ⇢ j-1 and u j <u j-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.

λ ROE 1 λ ROE 3 x t x c U a,L U ⇤ a U a,R

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. 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 (U a,L , U ⇤ a ,a n dU a,R ) connected by two discontinuities travelling at speeds λ ROE 1 and λ ROE 3 (cf. Figure 4.5). 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 [HLL83] applied to a 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 -F a (U j , U a,j )+F a (U j-1 , U a,j-1 )+S j-1/2 ⌘ , (4.25)
where the source term is naturally discretised as follows:

S j-1/2 = @ a λ ROE 1,j-1/2 (U ⇤ L,j-1/2 -U j-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 (U j -U ⇤ R,j-1/2 )d 3,j-1/2 .

HLLC-type scheme

Another possible approach for the sensitivity is to keep the same structure as for the state (cf. Figure 4.4), 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 .A g a i n , 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, created from the Rankine-Hugoniot jump relations:

8 > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > : -λ 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 +F a,R | 2 -F a,L | 2 = ∆xS| 2 , (λ 2 -λ 1 )(⇢E) ⇤ a,L +(λ 3 -λ 2 )(⇢E) ⇤ a,R + λ 1 (⇢E) a,L -λ 3 (⇢E) a,R +F a,R | 3 -F a,L | 3 = ∆xS| 3 , (4.26)
where λ 1 =ũ -c, λ 2 =ũ,a n dλ 3 =ũ +c. The first three equations are the Rankine-Hugoniot condition on ⇢ across the three waves, differentiated with respect to a.N o t e 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;t h el a s tt w o 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 = 0 B B B B B B @ b 1 b 2 b 3 b 4 b 5 b 6 1 C C C C C C A = 0 B B B B B B @ @ 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 -F a,R | 2 + F a,L | 2 ∆xS| 3 -λ 1 (⇢E) a,L + λ 3 (⇢E) a,R -F a,R | 3 + F a,L | 3 1 C C C C C C A
the system can be rewritten as:

Ax = b,
where A is the following matrix:

A = 0 B B B B B B @ -λ 1 01 00 0 λ 2 -λ 2 -11 00 0 λ 3 0 -10 0 -(γ-3) 2 ũ2 (γ-3) 2 ũ2 -(2 -γ)ũ (2 -γ)ũ -(γ -1) (γ -1) 00 c c 00 00 0 0 c c 1 C C C C C C A
and we have det(A)=4c 4 (γ -1) 6 =0. The solution of the system has the following form:

x = 0 B B B B B B B B B B B B B B B B @ (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 1 C C C C C C C C C C C C C C C C A .
An alternative strategy to compute U ⇤ a,L and U ⇤ a,R is to differentiate with respect to a the following relations:

U ⇤ L = U L + ↵ 1 r 1 , U ⇤ R = U R -↵ 3 r 3 , (4.27)
obtaining

U ⇤ a,L = U a,L + ↵ a,1 r 1 + ↵ 1 r a,1 , U ⇤ a,R = U a,R -↵ a,3 r 3 -↵ 3 r a,3 , (4.28)
with

r 1 = 0 @ 1 ũ - c H -ũc 1 A , r a,1 = 0 @ 0 ũa -ca Ha -ũa c -ũc a 1 A , r 2 = 0 @ 1 ũ ũ2 2 1 A , r a,2 = 0 @ 0 ũa ũũ a 1 A , r 3 = 0 @ 1 ũ +c H +ũc 1 A , r a,3 = 0 @ 0 ũa +c a Ha +ũ a c +ũc a 1 A , 8 > > > < > > > : ↵ 2 = γ -1 c2 h (⇢ R -⇢ L )( H -ũ2 )+ũ ⇣ (⇢u) R -(⇢u) L ⌘ - ⇣ (⇢E) R -(⇢E) L ⌘i , ↵ 1 = 1 c h (⇢ R -⇢ L )(ũ +c) - ⇣ (⇢u) R -(⇢u) L ⌘ -c↵ 2 i , ↵ 3 =(⇢ R -⇢ L ) -(↵ 1 + ↵ 2 ), 8 > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > : ↵ a,2 = - 2c a (γ -1) c3 h (⇢ R -⇢ L )( H -ũ2 )+ũ ⇣ (⇢u) R -(⇢u) L ⌘ - ⇣ (⇢E) R -(⇢E) L ⌘i + γ -1 c2 h (⇢ 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 ⌘i , ↵ a,1 = - ca c2 h (⇢ R -⇢ L )(ũ +c) -((⇢u) R -(⇢u) L ) -c↵ 2 i + 1 2c h (⇢ a,R -⇢ a,L )(ũ +c)+(⇢ R -⇢ L )(ũ a +c a ) -((⇢u) a,R -(⇢u) a,L ) -ca ↵ 2 -c↵ a,2 i , ↵ 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 (4.28) solve the system (4.26).

Proof. We will prove that the star sensitivities defined in (4.28) satisfy the system (4.26).

1. First equation. Writing the first coefficient of (4.27) one easily finds ⇢ ⇤ L -⇢ L = ↵ 1 , and writing the first two coefficients of (4.28) 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 (4.26) 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

U R -U L = 3 X i=1 ↵ i r i , U a,R -U a,L = 3 X i=1 ↵ a,i r i + ↵ i r a,i .
Therefore, one has:

U ⇤ R -U ⇤ L = ↵ 2 r 2 , U ⇤ a,R -U ⇤ a,L = ↵ a,2 r 2 + ↵ 2 r a,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 (4.26) 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 +c a ).
By replacing them in the third equation of (4.26) 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 (4.26) one can easily check that the equation is always verified, since λ 2 =ũ.

4.5. Convergence tests for the numerical schemes

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/4

Convergence tests for the numerical schemes

We consider the Riemann problem of section 4.3. In Figures 4.7-4.8-4.9 we show the convergence of the different numerical schemes presented in Section 4.4. For reference, the L 1 norms of the states ⇢ ex , u ex and p ex are respectively, 0.5625, 0.2204,a n d0.5354. Furthermore, those for the sensitivities ⇢ a,ex , u a,ex and p a,ex are 0.0379, 0.1768,and0.462. Figure 4.7 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 4.8-4.9 we plot the error for the sensitivity, first with the HLLtype scheme (Figure 4.8) and then with the HLLC-type scheme (Figure 4.9): 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 .I nF i g u r e4 . 1 0w ep l o tt h es o l u t i o na tt h efi n a lt i m eT =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 4.11. In Figure 4.12 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 4.13 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. 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.

Uncertainty Quantification

Problem description

In this section, we show how SA can be used for uncertainty quantification (UQ), cf.

[PNTG01, TPB01, Del14]. Many UQ techniques have been developed during the last decades: these methods can be either probabilistic or deterministic. SA falls into the second category while the most well-known of these techniques, the Monte Carlo method, is in the first. Other UQ techniques are polynomial chaos [Wal03, XK03, KM06] and the random space partition [START_REF] Abgrall | A semi-intrusive deterministic approach to uncertainty quantification in non-linear fluid flow problems[END_REF]. The first one is based on a decomposition of the stochastic part of the solution using an orthogonal polynomial basis. Then using a Galerkin method a new system of equations is derived, that provides the coefficients of different statistical quantities, allowing in this way the computation of the statistical moments of the output. Concerning the second approach, the main idea behind it is to consider the random parameters as variable and to solve the system with a finite volume method in a higher dimensional space, the new dimension being the sum of the spatial dimension and the random space dimension. 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].

The main aim of UQ is to determine a confidence interval for the output of a model, in our case U, given the uncertainty and the error on the input parameters. According to the AIAA definition, error has a deterministic nature, while uncertainty a stochastic one. Furthermore, uncertainty can be categorised in aleatoric uncertainty and epistemic uncertainty.

In this work, we compare two different UQ methods: Monte Carlo and sensitivity analysis. Both methods aim to provide statistical quantities like moments (mean, variance, ...) of the output of the model. In the following, X will represent one of the variables, considered as random variables, i.e. X can either be ⇢, u or p,a n dX 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 these 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 that the variable X will fall within the interval. For instance, the choice  =2provides an ⇠ 95% confidence interval.

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 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 X i=1 X i , σ 2 X = 1 N -1 N X i=1 (µ X -X i ) 2 .
These estimates are good if N is sufficiently large: the slow convergence, and therefore the high computational cost, is the main limitation of the Monte Carlo method. However, this method is readily parallelisable. 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,a n dσ a the covariance matrix:

µ a = 2 6 4 µ a 1 . . . µ a M 3 7 5 , σ a = 2 6 6 6 4 σ 2 a 1 cov(a 1 ,a 2 ) ... cov(a 1 ,a M ) cov(a 1 ,a 2 ) σ 2 a 2 ... cov(a 2 ,a M ) . . . . . . . . . cov(a 1 ,a M ) ... σ 2 a M 3 7 7 7 5
,
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 X i=1 (a i -µ a i )X a i (µ a )+o(kak 2 ).
Then computing the average, since X(µ a ) and X a i (µ a ) are not random variables, at first order one gets:

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

σ 2 X = E[(X(a) -µ X ) 2 ]=E 2 4 M X i=1 X a i (µ a )(a i -µ a i ) ! 2 3 5 = = M X i=1 X 2 a i (µ a )E[(a i -µ a i ) 2 ]+ M X i,j=1 i6 =j X a i (µ a )X a j (µ a )E[(a i -µ a i )(a j -µ 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 X i=1 X 2 a i σ 2 a i + M X i,j=1 i6 =j
X a i X a j 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].

Numerical results

We applied the uncertainty quantification techniques described in the previous subsection to the test case already introduced in section 4.3. The uncertain parameters are the left and right values of the physical variables for the state, i.e.:

a =(⇢ L , ⇢ R ,u L ,u R ,p L ,p R ) t ,
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.I nF i g u r e4 . 1 4w e 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 larger in that zone. In Figures 4. 15-4.16 we show the results of the SA approach, with ∆x = 10 -3 and the diffusive first order scheme, when the sensitivity is computed without the correction term (4.13): the spikes in the neighbourhood of the shock are very different with respect to the ones we get with the Monte Carlo approach. On one hand, these peaks lead to non-physical 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 the zone enough to contain the majority of the samples: one can observe that four out of five samples fall outside of the predicted interval in the neighbourhood of the shock. For these reasons, we suggest that the corrected sensitivity are more appropriate in this context. The results obtained with the corrected sensitivities are shown in Figure 4.17: the confidence interval obtained correspond to the ones obtained with the Monte Carlo approach, apart for the shock zone. Of course, the SA approach does not capture the uncertainty in the neighbourhood 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 SA approach, and it is the case with and without correction. However, the SA approach is less expensive: the Monte Carlo approach requires 1000 solutions of the state, whilst the SA approach requires only one solution of the state and as many solutions of the sensitivity as the number of uncertain parameters, in this case 6. Furthermore, the solution of the different sensitivities can be done in parallel. Finally, in Figure 4.18 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:t h ea n t i -d i ffusive 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.

| Quasi 1D Euler system

Introduction

State and sensitivity system

The quasi-1D Euler system is written: 8 < : @ t (h⇢)+@ x (h⇢u)=0, @ t (h⇢u)+@ x (h⇢u 2 + p)=p@ x h, @ t (h⇢E)+@ x (hu(⇢E + p)) = 0,

(5.1)
where ⇢ is the density, u is the velocity, ⇢E the total energy per mass unit, p the pressure, and h = h(x) > 0 is a smooth function of the space x and it is known. This system describes a flow in a nozzle of height h and it allows us to investigate more realistic applications, while remaining in the simpler and computationally less expensive onedimensional framework. In this work, we considered the following height:

h(x)= 8 < : 2 -A sin 2 ✓ x -x c `⇡ - ⇡ 2 ◆ x c -2 <x<x c + 2 , 2 otherwise, (5.2)
which is plotted in Figure 5.1. It is described by three parameters: A is the maximal depth, x c the point of maximal depth and `the length. The system is closed by the following algebraic equation:

p =(γ -1) ✓ ⇢E - 1 2 ⇢u 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 = q (γ -1)(H -1 2 u 2 ) . We can rewrite the system (5.1) in the vectorial form:

@ t (hU)+@ x (hF(U)) = P@ x h, (5.3)
where Since h(x) does not depend on time, one can formally rewrite the system (5.3) as:

U = 2 4 ⇢ ⇢u ⇢E 3 5 F(U)=
@ t U + @ x F(U)=(P -F(U)) @ x h h . (5.4)
We remark that the left-hand side corresponds to the classical Euler system.

To obtain the sensitivity system, we differentiate (5.4) with respect to the parameter a. Let h a be the derivative of the height of the channel with respect to a; then, the sensitivity equations are:

@ t U a + @ x F a (U, U a )=(P a -F a (U, U a )) @ x h h +(P -F(U)) h@ x h a -h a @ x h h 2 , (5.5)
where we recall that

F a (U, U a )=@ a F(U)= 2 4
(⇢u) a ⇢ a u 2 +2⇢uu a + p a u a (⇢E + p)+u a ((⇢E) a + p a ) 3 5 and P =(0,p a , 0) t . Once again, the left-hand side is the sensitivity of the Euler system, and as done previously, a correction term needs to be added to the right-hand side in order to consider the possible discontinuities:

@ t U a + @ x F a (U, U a )=(P a -F a (U, U a )) @ x h h +(P -F(U, U a )) h@ x h a -h a @ x h h 2 + S(U).
(5.6)

The source S is defined as follows:

S(U)= Ns X k=1 ⇢ ⇢ ⇢ k δ(x -x s,k (t)), (5.7)
where ⇢ ⇢ ⇢ k is the amplitude of the correction, the computation of which is the object of the next subsection.

formed just by averaging over each cell the solution of the Riemann problem computed in the previous step. At the moment, we have not considered anti-diffusive schemes for this kind of problem. This is mainly motivated by two reasons: first, numerical results in the previous chapter show how the role of numerical diffusion is very important to have the expected order of convergence, but is less significant when it comes to applications; secondly, for this system we are interested in stationary solutions and an anti-diffusive scheme such as the one used in the previous chapters can significantly slow down the convergence in time if, for instance, a shock keeps jumping back and forth between two adjacent cells. For similar reasons we do not consider higher order schemes for the moment.

Stationary solutions

The system (5.4) admits stationary solutions. As a first step in order to validate the code, we want to be able to reproduce the results obtained in [START_REF] Telib | The effect of shocks on second order sensitivities for the quasi-one-dimensional euler equations[END_REF][START_REF] Giles | Analytic adjoint solutions for the quasi-onedimensional Euler equations[END_REF] for two test cases with and without shock and described hereafter in details (subsections 5.3.2-5.3.3). Let us first address the boundary conditions definition.

Boundary conditions

Inlet and outlet boundary conditions have to be considered to compute the fluxes at domain extremities. In [START_REF] Telib | The effect of shocks on second order sensitivities for the quasi-one-dimensional euler equations[END_REF][START_REF] Giles | Analytic adjoint solutions for the quasi-onedimensional Euler equations[END_REF], they consider the triplet (H, p, p tot ) (p tot being the total pressure, defined as p tot = p + 1 2 ⇢u 2 ), and depending on the test case they impose the value of these physical quantities at the inlet and/or at the outlet in a fictitious external cell, before computing the boundary flux. The quantities not imposed are extrapolated from the interior adjacent cell. The fluxes are based on the conservative variables (⇢, ⇢u, ⇢E), therefore the first thing to do is to obtain the triplet (⇢, ⇢u, ⇢E) from (H, p, p tot ):

8 > > > > < > > > > : H = E + p ⇢ , p =(γ -1) ✓ ⇢E - 1 2 ⇢u 2 ◆ , p tot = p + 1 2 ⇢u 2 , ) 8 > > > > < > > > > : ⇢ = ✓ p tot + p γ -1 ◆ 1 H ⇢u = p 2⇢(p tot -p) ⇢E = p γ -1 + 1 2 ⇢u 2 .
(

We remark two things about this:

(i) this result is valid only if u>0, otherwise we should take ⇢u = -p 2⇢(p tot -p) . In our test cases we will always consider a positive velocity;

(ii) the square root p 2⇢(p tot -p) should not present any problem because ⇢ > 0 and p tot >p . However, this second inequality is valid cell by cell: if neither p nor p tot are imposed at the boundary, the inequality is naturally verified; if one wants to impose both of them, the constraint needs to be verified; finally, if we are in a case in which we impose only one of them (and the other one's value is taken from the first cell inside the domain) the inequality is not guaranteed to be verified. To avoid this problem, if p tot <pfor some time step, we force ⇢u =0.

Concerning the sensitivity boundary conditions, for each variable we impose the same type of condition as for the state: for instance, if ⇢ is imposed as inlet external state, the corresponding sensitivity ⇢ a is imposed as inlet external value as well. The values for (⇢ a , (⇢u) a , (⇢E) a ) are obtained by differentiating with respect to a the second system of (5.11):

8 > > > > > < > > > > > : ⇢ a = ✓ p tot,a + p a γ -1 -⇢H a ◆ 1 H , (⇢u) a = ⇢ a (p tot -p)+⇢(p tot,a -p a ) ⇢u , (⇢E) a = p a γ -1 + u(⇢u) a - 1 2 ⇢ a u 2 .
Once again, this is valid only if u>0.

Isentropic transonic case

The first test case, taken from [TAI11, GP01], that we want to reproduce is an isentropic transonic case, i.e. with the following boundary conditions for the state:

( H L =4, p tot,L =2,
and all the other variables are extrapolated from the interior solution. The value of the parameters describing h are A =1, `=0.5,a n dx c =0.5. We do not have an analytical solution for the state, however the analytical Mach number (Ma = u p ⇢ p γp ) is known and it is plotted in Figure 5.4 for the case considered, compared to the numerical one obtained with our numerical scheme: the two perfectly match. The Mach number is continuous in this case, as are all the physical quantities. The passage from the subsonic regime (Ma < 1) to the supersonic one (Ma > 1) happens in the neck of the nozzle, i.e. x = x c , as expected. To validate the numerical results obtained for the sensitivity, we compute the empirical sensitivity, defined as follows:

U emp a = U(a + δa) -U(a) δa ,
for a sufficiently small δa. In Figure 5.5 we show in colour the sensitivity obtained with our numerical schemes, and in black the empirical sensitivity computed with δa = 10 -3 , on a mesh of constant spatial step ∆x = 10 -3 with a first order Roe scheme. Furthermore, we consider that at the continuous level:

ERR(U):=kU(a + δa) -U(a) -δaU a (a)k L 1 , (5.12)
which is supposed to be O(δa 2 ). Figure 5.7 shows (5.12) for all the variables for different values of δa computed with the discrete solutions. We remark that for ⇢ and u we have the expected rate of convergence for the bigger δa: the slight change in the slope for smaller δa corresponds to the point in which O(∆x r ), r being the order of the scheme used, gets comparable with O(δa 2 ) and consequently the error due to the discretisation gets comparable with the one of the Taylor expansion. We recall that the spatial step is ∆x = 10 -3 and the scheme is a first order Roe scheme. Concerning ERR(p), we do not have the expected order of convergence even for larger δa. We remark however, that the case here considered presents a transonic rarefaction and it is well known that a Roe scheme is not able to capture that. To overcome this problem, an entropy fix has been implemented: however, a small jump still remains in the state at x =0.5 (cf. left side of Figure 5.6) and this generates the peak in the sensitivity noticeable in the right side of Figure 5.6, which deteriorates the convergence.

Transonic case with sho ck

The second test case we consider is a transonic case with shock, i.e. with the following boundary conditions for the state:

8 > < > : H L =4, p tot,L =2, p R =1.6,
and all the other variables are extrapolated from the interior solution. The value of the parameters describing h are the same as the previous test case, i.e. A =1, `=0.5,a n d x c =0 .5. We show in Figure 5.8 the comparison of the analytical and numerical Mach number for the case considered. As expected, the results are identical. The Mach number and all the physical variables present a discontinuity at x ' 0.6. As in the previous case, the transition from subsonic to supersonic happens at x = x c . In Figure 5.9 we plotted the sensitivity of the three variables with respect to a = H L , and we compare the scheme with and without correction term S(U): as one cas see, the two solutions are identical in the whole domain except for the spikes at x ' 0.6 which are present if the correction term is not taken into account. In this case, we do not compute ERR(U) for different δa, because it comes from a first order Taylor expansion, which of course is not valid in case of discontinuity. Computation of r a J(U, U

a new = a -rr a J 10:

if a new / 2 A then 11:
Projection of a new onto A 12:

end if

13:

State solution ! U 14:

Evaluation of J(U) ! J new 15:

while J new >J do J = J new a = a new 22: end while where p ⇤ is the target pressure.

Test case 1

For the first test case, we consider two parameters of optimization: A and `(cf. definition of h(x) (5.2) and Figure 5.1). The set of admissible parameters A is the rectangle A =[ 0 , 2) ⇥ [0, 1]. We remark that A is open (A =2is not admitted because we ask for h(x) > 0). This is not a problem because we will consider a test case with minimum inside the domain; on the other hand, `=0is to be intended as h(x)=const. 8x 2 (0, 1). The target pressure is the pressure obtained with A =1and `=0 .5:i nt h i sw a y ,t h e optimum is known and reachable and we aim at recovering this value with our algorithm. The gradient of the cost functional is the following:

r a J(U)=  (p -p ⇤ ,p A ) L 2 (p -p ⇤ ,p `)L 2 ,
where we used the notation (•, •) L 2 to indicate the L 2 dot product. The two sensitivities U A and U `can be computed by solving the system (5.6), with

h A (x)= 8 <
:

sin 2 ✓ x -x c `⇡ - ⇡ 2 ◆ x c -2 <x<x c + 2 , 0 otherwise, and h
`(x)= ( -⇡(x-xc) `2 sin 2⇡ x-xc ` x c -2 <x<x c + 2 , 0 otherwise.
Isentropic transonic flow. First we deal with the isentropic transonic case, already introduced in the previous section. In Figure 5.10 we show the cost functional: as one can see, it is quite flat in one direction, which can significantly slow down the convergence of the gradient descent method. Furthermore, one can remark a discontinuity: this is due to the fact that, when the nozzle is too deep (i.e. bigger A), a shock occurs in the state solution. In Figure 5.11 we show the optimization algorithm steps obtained with a first order Roe scheme with a constant spatial step ∆x = 10 -2 .W e c h o s e as starting parameters (A, `)=( 0 .5, 0.2). The algorithm stops after 154 iterations at (A, `)=( 1 .0003, 0.5002), with a value of the cost functional J =2 .2013 ⇥ 10 -9 .W e remark that there is no crossing of the cost functional discontinuity. In Figure 5.12 we show the value of the cost functional J and of the euclidean norm of its gradient r a J with respect to the number of iterations in a semi-logarithmic scale. One can notice that the decrease of the cost functional is faster at the beginning and then slows down significantly, as we expected, given the shape of the cost functional.

Shocked flow. The second optimization test case we present is the shocked flow case.

Here, because of the presence of a shock in the solution, we compare the results of two different numerical schemes: with and without the correction term (5.7). In Figure 5.13 we show the cost functional: it has a quite different shape with respect to the previous test caseIn this test case, first we computed the sensitivity without the correction term, i.e. by solving the equations (5.5), and then with the correction term, i.e. by solving (5.6). This leads to different gradients: in fact, if we think of the analytical solution for the sensitivities p A and p `, with and without the Dirac delta function, there is a missing term in the second case, when computing the integrals to obtain r a J(U). On the other hand, considering the discrete solution, there is an approximation error due to the Dirac delta function, which cannot be seized numerically. In Figure 5.14 we show the steps of the gradient descent method: in blue the results obtained without the correction and in red the results obtained with the correction. As one can see from the zoom (right side of the Figure ), the convergence is smoother if the gradient is computed with the corrected sensitivities. This can be seen also from Figure 5.15, in which the cost functional and its gradient are plotted in a semi-logarithmic scale with respect to the iterations: we remark that the difference in the value of the cost functional is unimportant, while the gap between the gradients is more significant. The algorithm with the corrected sensitivities converges in 169 iterations to the point (A, `)=(1.0003, 0.5002) with a value of the cost functional J =1 .2185 ⇥ 10 -8 , as opposed to the one with the uncorrected sensitivities, which stops in 196 iteration at the point (A, `)=(1.0003, 0.5001), with J =1.6821⇥10 -9 .

Test case 2

Here, we consider two other parameters of optimization: x c and `(cf. definition of h(x) (5.2) and Figure 5.1). The set of admissible parameters A is a closed triangle:

A = {(x c , `) 2 R 2 :0 x c  1, 0  ` min(2x c , 2 -2x c )},
in order for the nozzle not to exit the domain. We remark that `=0means h(x)= const. 8x 2 (0, 1). The target pressure is the pressure obtained with x c =0 .5 and `=0 .5: in this way, the optimum is again known and reachable. The gradient of the cost functional is the following:

r a J(U)=  (p -p ⇤ ,p xc ) L 2 (p -p ⇤ ,p `)L 2 .
The two sensitivities U xc and U `can be computed by solving the system (5.6), with

h xc (x)= ( -⇡ `sin 2⇡ x-xc ` x c -2 <x<x c + 2 , 0 otherwise, and
h `(x)= ( -⇡(x-xc) `2 sin 2⇡ x-xc ` x c -2 <x<x c + 2 , 0 otherwise.
Isentropic transonic flow. The first optimization test case we deal with is the isentropic transonic case, already introduced in the previous section. In Figure 5.16 we show the cost functional: as one can see from the left plot, changes in the parameter x c affect the value of J way more than changes in the parameter `. This leads to a cost functional which is quite flat in one direction and that can significantly slow down the convergence of the gradient descent method. On the right part of Figure 5.16, it is plotted the cost functional with respect to `for x c =0 .5: as one can see, the minimum is in (x c , `)=( 0 .5, 0.5). In Figure 5.17 we show the optimization algorithm steps obtained with a first order Roe scheme with a constant spatial step ∆x = 10 -2 .W e chose as starting parameters (x c , `)=(0.4, 0.4). The algorithm stops after 88 iterations at (x c , `)=(0.499998, 0.49974), with a value of the cost functional J =4.36936 ⇥ 10 -9 . In Figure 5.18 we show the value of the cost functional J and of the euclidean norm of its gradient r a J with respect to the number of iterations in a semi-logarithmic scale. One can notice that the decrease of the cost functional is faster at the beginning and then slows down significantly, as we expected, given the shape of the cost functional.

Concerning the oscillations of the gradient, they are due to the naïf choice of the step and a relaxation step would probably guarantee a faster convergence.

Shocked flow. The second optimization test case we present is the shocked flow case.

Here, because of the presence of a shock in the solution, we compare the results obtained with and without the correction term (5.7), with the diffusive Roe first order scheme. In Figure 5.19 we show the cost functional: againt it has a quite different shape with respect to the previous isentropic transonic test case. In this test case, first we computed the sensitivity without the correction term, i.e. by solving the equations (5.5), and then with the correction term, i.e. by solving (5.6). This leads to different gradients. In Figure 5.20 we show the steps of the gradient descent method: in blue the results obtained without the correction and in red the results obtained with the correction. As one can see from the zoom (right side of the Figure ), the convergence is faster if the gradient is computed with the corrected sensitivities. This can be seen also from Figure 5.21, in which the cost functional and its gradient are plotted in a semi-logarithmic scale with respect to the iteration: the algorithm with the corrected sensitivities converges in 9 iterations to the point (x c , `)=(0.5002, 0.498539) with a value of the cost functional J =4.50852 ⇥ 10 -7 , as opposed to the one with the uncorrected sensitivities, which stops in 36 iteration at the point (x c , `)=(0.502818, 0.477847),w i t hJ =1.12238 ⇥ 10 -4 .

| Conclusion and perspectives

The first goal of this thesis was to adapt some sensitivity techniques, and in particular the continuous sensitivity equation method, to the hyperbolic framework in case of discontinuous solutions. The first step was to define a sensitivity system valid also in the case of discontinuous state: this was achieved by adding a correction term to the sensitivity equations. The definition of this correction term was given for a general hyperbolic system and, in Chapter 5, we showed how this definition is valid also in the case of non-conservative systems. The scalar case of Chapter 2 allowed us to understand more in depth the Dirac delta functions exhibited by the sensitivities in the absence of the correction term, since we were able to compute analytically their coefficient.

The sensitivity system obtained with the correction term provided us with solutions without spikes approximating the Dirac delta function. This was an objective of this thesis, because Dirac delta function cannot be seized numerically and this corrupts the solution in the neighbourhood of the shocks and makes convergence impossible. Moreover, the spikes change with the numerical discretisation.

The second aim of this thesis was to design proper numerical schemes in order to discretise and approximate the solution of the sensitivity system obtained. This did not present any particular issue in the scalar case of Chapter 2, however the numerical results presented in Chapters 3-4-5 showed that, in the case of systems, the numerical diffusion plays a very important role in the discretisation of the sensitivity, 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 intermediate zones is not correct.

To overcome this problem, we proposed some anti-diffusive numerical schemes, based on sampling techniques and inspired by Glimm random choice method: with these schemes we were able to discretise more precisely the source term and to obtain a correct solution for the sensitivity.

Two applications were tackled: an uncertainty quantification and an optimization problem. Both of them underlined the importance of the correction term: on one hand, for the uncertainty quantification the results obtained without the correction term provide non-physical estimates for some of the variables, without adding any useful information on the shock displacement; on the other hand, the optimization algorithm converges more quickly and with a smoother trajectory if the corrected sensitivity are employed to compute the gradient of the cost functional. Concerning the anti-diffusive schemes, they were only tested on the uncertainty quantification problem and they did not provide better results, in spite of their higher computational cost. This is a very interesting point because anti-diffusive schemes such as Glimm method are not extendible in 2D: comparing the results obtained with diffusive and anti-diffusive schemes, we can say that the firsts are suitable for more complex applications, even though they do not provide a solution for the sensitivity as precise as one would expect with respect to the mesh.

This work leads to some new interesting problems. Firstly, it would be interesting to investigate higher order sensitivities: the first order sensitivity exhibits more discontinuities than the state, therefore the correction term for a second order sensitivity would be even more complex. Moreover, the definition of such a correction term would add constraints of the numerical schemes used to solve the first order sensitivity. Secondly, to tackle more realistic fluid dynamics application, one should extend the methods here developed in higher spatial dimension. From a theoretical point of view, it should be straightforward. The main practical difficulties that we can foresee are the definition of a good shock detector and, if needed, the 2D extension of the anti-diffusive numerical schemes. To do this, one would need a parametric description of the shock and then some shock fitting techniques [START_REF] Paciorri | A shock-fitting technique for 2D unstructured grids[END_REF] could be used in this framework. Finally, one could perform SA on conservation laws presenting non-classical shocks [START_REF] Lefloch | Hyperbolic Systems of Conservation Laws: The theory of classical and nonclassical shock waves[END_REF] and study the sensitivity of the solutions with respect to the underlying kinetic function, which allows to select the admissible solutions.

It follows a more detailed conclusion, chapter by chapter.

Chapter 2. In this chapter we considered a scalar nonlinear hyperbolic PDE. Starting from the scalar case allowed us to introduce some of the main problems in a simpler framework. In the first part of the chapter, the analytical solution for the state equation was computed using the method of characteristics, which provided us with an implicit solution, unless the initial data was particularly simple, like in the case of a Riemann problem. Secondly, the sensitivity system was defined and its analytical solution was obtained starting from the state solution: interestingly, for the sensitivity we were able to give an explicit solution, although depending on the state, for all initial data. Then, the correction term was defined and the hyperbolicity of the global system was briefly analysed. We remarked that the solution to the sensitivity equation without the correction term presented a Dirac delta function where the state is discontinuous: in this simple scalar case, we were able to compute analytically the coefficient associated to the Dirac and one could notice that the spikes were larger if the dependence of the position of the shock on the parameter of interest was stronger. Finally, a numerical scheme was presented, along with the results obtained for the inviscid Burgers' equation with different initial data: a Riemann problem and a continuous initial conditions. We observed that, for a scalar equation, numerical diffusion does not play a fundamental role, because the problem is too simple: for instance, for the Riemann problem there is no intermediate state to be computed. On the other hand, the test case of a continuous initial data was complicated enough to tackle the problem of the definition of a shock detector: in this case it is way more difficult to numerically detect a discontinuity than in the case of a piecewise constant function (which is the solution to a scalar Riemann problem).

Chapter 3. In this chapter, which is an adaptation of [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], we dealt with the Euler equations in barotropic conditions and in Lagrangian coordinates. It is a system of two equations with two unknowns. It describes the dynamics of a compressible material in barotropic conditions (i.e. the density is a function only of the pressure), using two physical variables: the co-volume ⌧ , which is the reciprocal of the density, and the Lagrangian speed. The choice of the Lagrangian coordinates was motivated by the fact that in this way the sign of the eigenvalues of the system is known, which simplified the design of numerical schemes. In the first sections, we presented the state equations and derived the sensitivity equations in the regular case. Next, a section was devoted to the definition of the correction term in order to have a sensitivity system which was valid also in the case of discontinuous state. The correction term was defined by integrating the sensitivity equation over a control volume and then using the Rankine-Hugoniot conditions, which govern the state across a shock. We then detailed the exact resolution of the Riemann problem for the state and for the corrected sensitivity system: the state presents a classical structure, composed by two waves which can be either shocks or rarefaction waves; contrarily, the sensitivity has a quite particular structure, with two discontinuities associated with the beginning and the end of the rarefaction wave and a constant plateau inside the rarefaction itself. Then, some classical numerical methods for the state and their adaptation for the sensitivity were examined. In particular an exact Godunov scheme was implemented for both state and sensitivity, and a Roe approximate Riemann solver was implemented for the state and an adaptation of it for the sensitivity. We extended the Roe-type schemes to the second order, with a two-step Runge-Kutta method in time and a MUSCL-type approach in space. However, the second order extension in space is not classical: a precaution was necessary in order to be able to discretise the source term in the second order framework. Some numerical tests were conducted, which exhibited mesh-convergence issues, regardless of the order of the scheme. In particular, the test case of an isolated shock led us to conclude that the convergence problem was due to numerical diffusion. For this reason, an anti-diffusive version of the Roe-type schemes mentioned above was introduced, inspired to Glimm's random choice method, both for the first and second order. This allowed to overcome the mesh-convergence problems, as shown in some numerical results in the last section of the chapter.

Chapter 4. In this chapter we studied the complete Euler system. It is a system of three PDEs with four unknowns, closed by an algebraic equation. One can choose different sets of variables to describe the flow. A common choice is, for instance, to work with the conserved variables: the density of the fluid ⇢, the momentum ⇢u and the total energy per unit mass ⇢E. In this work, we presented all the results with the physical, or primitive, variables: the density ⇢, the speed u and the pressure p.I nt h efi r s tp a r t of the chapter the state system was described and the sensitivity system was derived in the regular case. It followed a detailed analysis on the hyperbolicity of the global system of state and sensitivity as a whole and we proved that, in the general case, the system is only weakly hyperbolic. Then, the case of discontinuous state was considered and the correction term for the sensitivity equations was computed, in a similar way to the previous chapter. The solution of a specific Riemann problem, known as the Sod shock tube, was detailed for the state. The choice of this case was motivated by the presence of three different waves: a rarefaction, a contact discontinuity and a shock. Starting from the analytical solution for the state, the analytical sensitivity was derived: as we observed in the previous chapter, the sensitivity exhibits two discontinuities, one at the beginning and the other at the end of the rarefaction wave; however in this case there is a proper rarefaction wave for the sensitivity, too. In the following part of the Chapter, some numerical schemes were presented: we illustrated some constraints regarding the schemes that could be used for the state. In particular, what we categorised as HLLtype Riemann solvers (i.e. with only one intermediate star state) are not suitable for our problem: if used, they would make impossible the definition and the computation of the correction term across the contact discontinuity. For this reason, the approximate Riemann solver of Roe was used. Concerning the schemes for the sensitivity, first we introduced a simpler HLL scheme, we then illustrated two approaches to obtain an HLLCtype scheme (i.e. with two different intermediate zones) and proved that they both lead to the same scheme. Diffusive and anti-diffusive versions of all of the schemes mentioned were implemented at first and second order, as done in the previous chapter, and some mesh-convergence tests were performed. The results obtained were consistent with the ones obtained in the previous chapter: the numerical diffusion affects the plateau values and this corrupts the convergence. Concerning the sensitivity, the computationally most expensive HLLC-type scheme did not provide a more accurate solution.

In the final part of the Chapter, the first application was tackled. A UQ analysis was performed on the Riemann problem: confidence intervals were calculated for all the physical variables, using the sensitivities. We computed the sensitivities in three different ways: with and without the correction term, and with the diffusive and anti-diffusive numerical schemes. We showed the confidence intervals obtained with the well-known Monte Carlo method and with the SA methods for each physical variable. In the regular zones the two approaches gave the same results, with a significantly lower computational cost of the SA method. In the neighbourhood of the discontinuity there was a loss of precision of the SA methods, due to the fact that the influence of the uncertain parameters on the position of the shock was neglected. However, the corrected sensitivities gave more reasonable results: the sensitivities obtained without the correction term, due to the presence of the spikes, provided non-physical confidence intervals for the density and the pressure. An important observation is that the anti-diffusive schemes did not provide a significantly different solution: this is interesting in the perspective of more realistic, two-dimensional applications because the random choice methods, on which our anti-diffusive methods are based, are not easily extendible in 2D.

Chapter 5. In this chapter, we investigated the quasi 1D Euler system: this is an extension of the Euler system considered in the previous chapter, but it models a flow happening in a channel whose section is not constant. To describe that, a smooth function of the space, h(x; a), depending on some parameters a, was introduced and a term depending on h(x; a) and on its derivative in space was added to the right-hand side of the Euler system. Once again, the sensitivity system was defined at first in the regular case and then a correction term was added to account for possible discontinuities in the state. Since the state system considered in this chapter was not conservative, the computation of the correction term was slightly different, because it could not rely on the Rankine-Hugoniot conditions. Nonetheless the same conclusion as in the previous cases was reached. As the conservative part of the system coincided with the Euler system of the previous chapter, the same numerical schemes were used. However, in light of the results obtained, in particular the UQ ones, we chose to focus our attention on the first order diffusive Roe scheme. A brief description of the discretisation of the non-conservative part was done. This system admits stationary solutions. We considered in particular two test cases: an isentropic transonic case and a case presenting a shock. The boundary conditions for these cases were discussed in details. We compared the numerical results obtained for the state to the analytical solution given in [START_REF] Giles | Analytic adjoint solutions for the quasi-onedimensional Euler equations[END_REF] for the Mach number. Concerning the sensitivity, we do not have an analytical solution for this case, therefore we compared the numerical sensitivities obtained with our schemes with the empirical sensitivities. The empirical sensitivities were computed by finite differences, starting from two different realisations of the state for two close enough values of the parameter of interest.

In the final part of the Chapter, a second application was tackled: we defined an optimization problem. The optimization method used was a projected gradient algorithm which used sensitivities to compute the gradient of the cost functional. The problem considered is an inverse problem, and in particular a pressure matching one: we set a target pressure p ⇤ , which is a pressure obtained with a certain vector of parameters a ⇤ , and we aimed at recovering this value with the optimization algorithm. We performed the same test in the isentropic transonic case and in the case with shock, for two different choices of optimization parameters. In the case with shock, we compared the results obtained if the sensitivities, and therefore the gradient, were computed with or without the correction term. The results showed that the correction term speeds up the convergence of the optimization algorithm, and therefore confirmed the importance of the correction.

Conclusion et perspectives

Le premier objectif de cette thèse était d'adapter certaines techniques d'AS, et en particulier la méthode de l'équation de sensibilité continue, au cadre hyperbolique en cas de solutions discontinues. La première étape consistait à définir un système de sensibilité valable également en cas d'état discontinu : ceci a été réalisé en ajoutant un terme de correction aux équations de sensibilité. La définition de ce terme de correction a été donnée pour un système hyperbolique général et, dans le chapitre 5, nous avons montré que cette définition est valable aussi dans le cas de systèmes non conservatifs. Le cas scalaire du chapitre 2 nous a permis de comprendre plus en détail les Dirac présents dans les sensibilités en l'absence du terme de correction, puisque nous avons pu calculer analytiquement leur amplitude.

Le système de sensibilité obtenu avec le terme de correction nous a fourni des solutions sans pics approchant les Dirac. Ceci était un objectif de cette thèse, car une distribution de Dirac ne peut pas être évaluée numériquement et cela dégrade la solution au voisinage des chocs et rend la convergence vers la solution analytique impossible. De plus, les pics changent avec la discrétisation numérique.

Le deuxième objectif de cette thèse était de définir des schémas numériques appropriés afin de discrétiser et d'approximer la solution du système de sensibilité obtenu. Cela n'a pas posé de problèmes particuliers dans le cas scalaire traité dans le chapitre 2, cependant les résultats numériques présentés dans les chapitres 3-4-5 ont montré que, dans le cas des systèmes, la diffusion numérique joue un rôle très important dans la discrétisation de la sensibilité, à tel point que les schémas volumes finis classiques ne convergent pas vers la solution analytique : en particulier, la valeur du plateau dans les zones intermédiaires n'est pas correcte.

Pour surmonter ce problème, nous avons proposé des schémas numériques antidiffusifs, basés sur des techniques d'échantillonnage et inspirés de la méthode de choix aléatoire de Glimm : avec ces schémas, nous avons pu discrétiser plus précisément le terme source et obtenir la bonne solution pour la sensibilité.

Nous nous sommes attaqués à deux applications : un problème de quantification d'incertitude et un d'optimisation. Les deux ont souligné l'importance du terme de correction : d'une part, pour la quantification d'incertitude, les résultats obtenus sans le terme de correction fournissent des estimations non physiques pour certaines variables, sans cependant ajouter d'informations utiles sur le déplacement du choc ; d'autre part, l'algorithme d'optimisation converge plus rapidement et avec une trajectoire plus lisse si la sensibilité corrigée est utilisée pour calculer le gradient de la fonctionnel coût. Concernant les schémas anti-diffusifs, nous ne les avons testés que sur le problème de quantification d'incertitude et ils n'ont pas donné de meilleurs résultats, malgré un coût de calcul plus élevé. C'est un point très intéressant car les schémas anti-diffusifs comme la méthode de Glimm ne sont pas facilement extensibles en deux dimensions : en comparant les résultats obtenus avec des schémas diffusifs et anti-diffusifs, on peut dire que les premiers sont adaptés à des applications plus complexes, même si ils ne fournissent pas une solution de sensibilité au niveau de précision que l'on s'attend par rapport au maillage.

Ce travail mène à des nouveaux problèmes intéressants. D'abord, il serait intéressant d'étudier des sensibilités d'ordre supérieur : la sensibilité de premier ordre présente plus de discontinuités que l'état, donc le terme de correction pour une sensibilité de second ordre serait encore plus complexe. De plus, la définition d'un tel terme de correction ajouterait des contraintes aux schémas numériques utilisés pour résoudre la sensibilité de premier ordre. Deuxièmement, pour aborder une application plus réaliste de la dynamique des fluides, on devrait étendre les méthodes développées ici en dimension spatiale supérieure. D'un point de vue théorique, cela devrait être simple. Les principales difficultés pratiques que nous pouvons prévoir sont la définition d'un bon détecteur de choc et, si nécessaire, l'extension 2D des schémas numériques anti-diffusifs. Pour cela, il faudrait avoir une description paramétrique du choc et les techniques de shock fitting [PB09] pourraient être adaptées à ce contexte. Enfin, on pourrait appliquer l'AS à des lois de conservation présentant des chocs non classiques [START_REF] Lefloch | Hyperbolic Systems of Conservation Laws: The theory of classical and nonclassical shock waves[END_REF] et étudier ainsi la sensibilité des solutions par rapport à la fonction cinétique sous-jacente qui permet de sélectionner les solutions admissibles.

Ci-dessous une conclusion plus détaillée, chapitre par chapitre.

Chapitre 2. Dans ce chapitre, nous avons considéré une EDP hyperbolique non linéaire scalaire. Partir du cas scalaire nous a permis d'affronter certains des problèmes principaux dans un cadre plus simple. Dans la première partie du chapitre, la solution analytique pour l'équation d'état a été calculée en utilisant la méthode des caractéristiques, qui nous a fourni une solution implicite, sauf dans des cas de donnée initiale particulièrement simple, comme un problème de Riemann. Deuxièmement, le système de sensibilité a été défini et sa solution analytique a été obtenue à partir de la solution d'état : il est intéressant de noter que, pour la sensibilité, nous avons pu donner une solution explicite, en fonction de l'état, pour toute donnée initiale. Ensuite, le terme de correction a été défini et l'hyperbolicité du système global a été brièvement analysée. Nous avons observé que la solution de l'équation de sensibilité sans le terme de correction présentait un Dirac là où l'état était discontinu. Dans ce cas scalaire simple, nous avons pu calculer analytiquement le coefficient associé au Dirac et nous avons constaté que les pics étaient plus petits si l'influence du paramètre d'intérêt sur la position du choc était faible. Enfin, un schéma numérique a été présenté, ainsi que les résultats obtenus pour l'équation de Burgers non visqueux avec différentes données initiales : un problème de Riemann et un problème avec conditions initiales continues. Nous avons observé que, pour une équation scalaire, la diffusion numérique ne joue pas un rôle fondamental, car le problème est trop simple : par exemple, pour le problème de Riemann, il n'y a pas d'état intermédiaire à calculer. D'autre part, le cas d'une donnée initiale continue était suffisamment compliqué pour aborder le problème de la définition d'un détecteur de choc : dans ce cas, il est beaucoup plus difficile de détecter numériquement une discontinuité que dans le cas d'une fonction constante par morceaux, qui est la solution du problème scalaire de Riemann. Chapitre 3. Dans ce chapitre, qui est une adaptation de [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], nous avons traité les équations d'Euler dans des conditions barotropes et en coordonnées lagrangiennes. Il s'agit d'un système de deux équations à deux inconnues. Il décrit la dynamique d'un matériau compressible dans des conditions barotropes (c'est-à-dire la densité est une fonction seulement de la pression), utilisant deux variables physiques : le co-volume ⌧ , qui est l'inverse de la densité, et la vitesse lagrangienne. Le choix des coordonnées lagrangiennes a été motivé par le fait que l'on connaît le signe des valeurs propres du système, ce qui a simplifié la conception des schémas numériques. Dans les premières sections, nous avons présenté les équations d'état et nous avons dérivé les équations de sensibilité dans le cas régulier. Ensuite, une section a été consacrée à la définition du terme de correction, pour avoir un système de sensibilité valable également dans le cas d'état discontinu. Le terme de correction a été défini en intégrant l'équation de sensibilité sur un volume de contrôle, puis en utilisant les conditions de Rankine-Hugoniot, qui gouvernent l'état à travers un choc. Nous avons ensuite détaillé la résolution analytique du problème de Riemann pour l'état et pour le système de sensibilité corrigé : l'état présente une structure classique, composée par deux ondes qui peuvent être soit des chocs, soit des ondes de détente ; au contraire, la sensibilité a une structure assez particulière, avec deux discontinuités associées au début et à la fin de l'onde de détente et un plateau constant à l'intérieur de la détente elle-même. Ensuite, quelques méthodes numériques classiques pour l'état et leur adaptation pour la sensibilité ont été examinées. En particulier un schéma de Godunov exact a été implémenté pour l'état et la sensibilité, et un solveur de Riemann approché, le solveur de Roe, a été implémenté pour l'état et une adaptation de celui-ci pour la sensibilité. Nous avons étendu les schémas de type Roe au second ordre, avec une méthode Runge-Kutta en temps et une approche de type MUSCL en espace. Cependant, l'extension au second ordre en espace n'est pas classique : une précaution était nécessaire pour pouvoir discrétiser le terme source au second ordre. Des tests numériques ont été faits, qui présentaient des problèmes de convergence en maillage, indépendamment de l'ordre du schéma. En particulier, le cas test d'un choc isolé nous a mené à conclure que le problème de convergence était dû à la diffusion numérique. Pour cette raison, une version anti-diffusive des schémas de type Roe mentionnés ci-dessus a été introduite, inspirée de la méthode de choix aléatoire de Glimm, à la fois pour le premier et le second ordre. Cela a permis de surmonter les problèmes de convergence en maillage, comme le montrent certains résultats numériques dans la dernière section du chapitre. Chapitre 4. Dans ce chapitre, nous avons étudié le système d'Euler complet. C'est un système de trois EDP à quatre inconnues, fermé par une équation algébrique. On peut choisir différentes variables pour décrire le flux. Un choix commun est, par exemple, de travailler avec les variables conservatives : la densité du fluide ⇢, la quantité de mouvement ⇢u et l'énergie totale par unité de volume ⇢E. Dans ce travail, nous avons présenté tous les résultats avec les variables physiques, ou primitives : la densité ⇢, la vitesse u et la pression p. Dans la première partie du chapitre, nous avons présenté le système d'état et nous avons dérivé le système de sensibilité dans le cas régulier. Il a suivi une analyse détaillée sur l'hyperbolicité du système global formé par l'état et la sensibilité et nous avons prouvé que, dans le cas général, le système n'est que faiblement hyperbolique. Ensuite, nous avons considéré le cas d'état discontinu et nous avons calculé le terme de correction pour les équations de sensibilité, de la même façon que dans le chapitre précédent. La solution d'un problème de Riemann spécifique, connu sous le nom de tube à choc de Sod, a été détaillée pour l'état. Le choix de ce cas est motivé par la présence de trois ondes différentes : une détente, une discontinuité de contact et un choc. À partir de la solution analytique de l'état, nous avons dérivé celle de la sensibilité : comme nous l'avons observé dans le chapitre précédent, la sensibilité présente deux discontinuités, une au début et une à la fin de l'onde de détente ; cependant, dans ce cas, il y a une vraie onde de détente aussi pour la sensibilité. Dans la partie suivante du chapitre, des schémas numériques ont été présentés : nous avons illustré certaines contraintes concernant les schémas qui peuvent être utilisés pour l'état. En particulier, ce que nous avons classé comme solveurs de Riemann de type HLL (c'est-à-dire avec un seul état étoile intermédiaire) ne conviennent pas à notre problème : si utilisés, ils rendraient impossible la définition et le calcul du terme de correction à travers la discontinuité de contact. Pour cette raison, le solveur de Riemann approché de Roe a été utilisé. Concernant les schémas de sensibilité, nous avons d'abord introduit un schéma HLL classique, puis nous avons illustré deux approches pour obtenir un schéma de type HLLC (c'est-à-dire avec deux zones intermédiaires différentes) et montré qu'elles conduisaient au même schéma. Des versions diffusives et anti-diffusives de tous les schémas mentionnés ont été implémentées au premier et au second ordre, comme dans le chapitre précédent, et des tests de convergence en maillage ont été également réalisés. Les résultats obtenus sont en accord avec ceux obtenus dans le chapitre précédent : la diffusion numérique affecte les valeurs du plateau et cela dégrade la convergence. En ce qui concerne la sensibilité, le schéma de type HLLC, plus coûteux en termes de calcul, n'a pas fourni une solution plus précise.

Dans la dernière partie du chapitre, nous nous sommes attaqués à la première application. Une analyse de quantification d'incertitude a été effectuée sur le problème de Riemann : nous avons calculé les intervalles de confiance pour toutes les variables physiques, en utilisant les sensibilités. Nous avons calculé les sensibilités de trois manières différentes : avec et sans le terme de correction, et avec les schémas numériques diffusifs et anti-diffusifs. Nous avons montré les intervalles de confiance obtenus avec la méthode Monte Carlo et avec les méthodes d'AS pour chaque variable physique. Dans les zones régulières, les deux approches ont donné les mêmes résultats, avec un coût de calcul significativement plus bas pour la méthode d'AS. Au voisinage de la discontinuité, il y avait une perte de précision des méthodes d'AS, due au fait que l'influence des paramètres incertains sur la position du choc a été négligée. Cependant, les sensibilités corrigées ont donné des résultats plus raisonnables : les sensibilités obtenues sans le terme de correction, en raison des pics, ont fourni des intervalles de confiance non physiques pour la densité et la pression. Une observation importante est que les schémas anti-diffusifs ne donnent pas une solution significativement différente : ceci est intéressant dans la perspective d'applications bidimensionnelles plus réalistes car les méthodes de choix aléatoires, sur lesquelles sont basées nos méthodes anti-diffusives, ne sont pas facilement extensibles en 2D. Chapitre 5. Dans ce chapitre, nous avons étudié le système d'Euler quasi-1D : il s'agit d'une extension du système d'Euler considéré dans le chapitre précédent, mais il modélise un écoulement dans un canal dont la section n'est pas constante. Pour décrire cela, une fonction lisse de l'espace, h(x; a), dépendante d'un vecteur de paramètres a,a été introduite et un terme dépendant de h(x; a) et de sa dérivée en espace ont été ajoutés à la partie droite du système d'Euler. Encore une fois, le système de sensibilité a été défini d'abord dans le cas régulier, puis un terme de correction a été ajouté pour tenir compte des éventuelles discontinuités dans l'état. Comme le système d'état considéré dans ce chapitre n'était pas conservatif, le calcul du terme de correction était légèrement différent, car il ne pouvait pas s'appuyer sur les conditions de Rankine-Hugoniot. Néanmoins, la même conclusion que dans les cas précédents a été obtenue. Comme la partie conservative du système coïncidait avec le système d'Euler du chapitre précédent, les mêmes schémas numériques ont été utilisés. Cependant, vus les résultats obtenus, en particulier ceux de la quantification d'incertitude, nous avons choisi de concentrer notre attention sur le schéma de Roe diffusif au premier ordre. Une brève description de la discrétisation de la partie non conservative a été faite. Ce système admet des solutions stationnaires. Nous avons considéré en particulier deux cas tests : un cas isentropique transsonique et un cas transsonique avec un choc. Les conditions aux limites de ces deux cas ont été discutées en détail. Nous avons comparé les résultats numériques obtenus pour l'état à la solution analytique donnée dans [START_REF] Giles | Analytic adjoint solutions for the quasi-onedimensional Euler equations[END_REF] pour le nombre de Mach. Concernant la sensibilité, nous n'avons pas de solution analytique pour ce cas, donc nous avons comparé les sensibilités numériques obtenues avec nos schémas avec les sensibilités empiriques. Les sensibilités empiriques ont été calculées par des différences finies, à partir de deux réalisations différentes de l'état pour deux valeurs proches du paramètre d'intérêt.

Dans la dernière partie du chapitre, une deuxième application a été abordée : nous avons défini un problème d'optimisation. La méthode d'optimisation utilisée est un algorithme de gradient projeté qui utilise des sensibilités pour calculer le gradient de la fonctionnelle coût. Le problème considéré est un problème inverse pour la pression : nous établissons une pression cible p ⇤ , qui est une pression obtenue avec un certain vecteur de paramètres a ⇤ , et nous visons à récupérer cette valeur avec l'algorithme d'optimisation. Nous avons effectué le même test dans le cas isentropique transsonique et dans le cas transsonique avec choc, pour deux choix différents de paramètres d'optimisation. Dans le cas du choc, nous avons comparé les résultats obtenus si les sensibilités, et donc le gradient, étaient calculés avec ou sans le terme de correction. Les résultats ont montré que le terme de correction accélère la convergence de l'algorithme d'optimisation, confirmant ainsi l'importance de la correction.

force. The aim is to solve (A.1)-(A.2)-(A.3) in such a way that, given a distance D,t h e final time T is minimal. From a mathematical point of view, it is a problem of optimal control: the propulsive force f is the control variable and the time T is the cost functional to be minimised, which depends on f through the states variables v and e. Therefore, the problem can be written as follows: where F is the set of admissible controls:

F = {f :0  f (t)  f M 8t ≥ 0}.
In his work, Keller [START_REF] Keller | Optimal velocity in a race[END_REF] claims that his energy balance takes into account only the aerobic energy (i.e. energy provided by oxygen consumption). However, Aftalion and Bonnans [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF] remark that what he encompasses in the balance is the accumulated oxygen deficit: e 0 -e(t). Therefore, e(t) in equation (A.1) is in fact the anaerobic energy (i.e. energy provided by glycogen and lactate). In order to reproduce the results of [START_REF] Hanon | Oxygen uptake and blood metabolic responses to a 400-m run[END_REF][START_REF] Hanon | Pacing strategy and VO2 kinetics during a 1500-m race[END_REF][START_REF] Hanon | Effects of optimal pacing strategies for 400-, 800-, and 1500-m races on the VO2 response[END_REF], the oxygen uptake σ introduced in [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF] is piecewise defined: in the most part of the race, σ is constant equal to its maximal value σ max , but it is increasing at the beginning of the race, and decreasing at the end (see The oxygen uptake σ as defined in (A.5) is continuous but not C 1 . In the numerical simulations, it has been smoothed, since from the physiology, it is clear that the passage from one zone to the other occurs smoothly. The sigma used is shown in Figure A.1. Let us observe that, consistently with Hanon experimental results [START_REF] Hanon | Oxygen uptake and blood metabolic responses to a 400-m run[END_REF], the functional form of σ does not change with the athlete, whilst the values of σ max and σ f do. Moreover, in [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF], a second modification is introduced to the energy equation: for sufficiently long races (longer than 1000m), it has been observed that slowing down recreates some of the anaerobic energy. Therefore, the energy equation results in: ė = σ(e)+η(v) -fv, where η depends on the acceleration v and has the following form:

η(v)= ( 0 if v>0 c ⌘ | v| 2 if v  0, (A.6)
where c ⌘ is a constant to be tuned. This leads to oscillations in the velocity profile (cf.

[1], Figure 2.4). In this paper, we will not focus on the term η and on the causes of the oscillations: we will briefly present some numerical results obtained with the term η in one simple case. For further details, see [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF].

A.2.2 Equations for two runners

In real races, the competition between runners has a fundamental impact on the strategy.

Starting from the works of Keller [START_REF] Keller | Optimal velocity in a race[END_REF], Quinn [START_REF] Quinn | The effects of wind and altitude in the 400-m sprint[END_REF] and Kyle [START_REF] Kyle | Reduction of wind resistance and power output of racing cyclists and runners travelling in groups[END_REF], Pitcher introduces a two-runners model in [START_REF] Pitcher | Optimal strategies for a two-runner model of middle-distance running[END_REF] based on the slipstream. The observation is that running behind someone can save 1 or 2 seconds per lap in middle distance races. Therefore, in the equations, the friction gets reduced when runners are just behind their competitor. Some of the quantities in this model have a subscript i, which refers to the runner i: therefore x i , v i and e i are, respectively, the position, the velocity and the energy of runner-i. All the physiological parameters which depend on the runner have the subscript i,too. Finally ,thestatev ariablex D represents the distance between the runners, defined as x D := x 2 -x 1 . The energy balance for each runner is the same as in (A.1), with a constant value of σ, while the dynamics equation incorporates an aerodynamical term. Because it is the relative position which is important, instead of using x i ,t h epo s i t i o no f each runner as parameters, we use x 1 and the relative position x D . The resulting model is the following: for i =1, 2 8 > > > > > > < > > > > > > :

ẋ1 = v 1 x 1 (0) = 0 ẋD = v 2 -v 1 x D (0) = 0 v1 = f 1 -v 1 ⌧ 1 -c 1 v 2 1 (1 -γ(e -↵(x D -β) 2 )) v 1 (0) = 0 v2 = f 2 -v 2 ⌧ 2 -c 2 v 2 2 (1 -γ(e -↵(x D +β) 2 )) v 2 (0) = 0 ėi = σ i -f i v i e i (0) = e 0 i .
(A.7)

As in Keller's model, the velocities and energies equations are normalized with respect to the mass of the runner. The term -c i v 2 i is the friction with the air. It is necessary to highlight and separate the effect of the friction with the air from the other ones, because it is the only one that is reduced while running in the slipstream of someone else. This frictional term is modulated by 1 -γ(e -↵(x D ±β) 2 ), which is shown in The parameter β represents the optimal distance a runner should keep from the other in order to obtain the maximal reduction of the air friction, while γ 2 (0, 1) is the The initial strategy is the same as the previous case: runner-2 slows down in order to stay behind; this allows them to keep the same velocity as runner-1 using a smaller force and therefore compensating the smaller σ. The speed up in the final part is less evident than it was in the previous case, because the difference between the energies is smaller. The final distance covered by runner-1 in this case is 1499.72m. The final time of runner-2 is 249.66s, compared to 251.665s if running alone. Figure A.11 shows the distance between the runners during the race: as in the previous case, the overtaking occurs late in the race.

Figures A.12-A.13 show the results for two runners who have a different value for τ . The threshold values are: τ 1 =1.33s τ 2 =1.31s.

A smaller τ indicates a bigger drop in velocity due to frictional effects, therefore, a greater force is necessary to keep the same velocity. Being behind another runner is a way to compensate this weakness: as shown in Figure A.12, runner-2 manages to keep the same velocity as runner-1, using a slightly smaller force, in spite of having a smaller τ . In the final part, when the difference in the energies is sufficiently big, runner-2 overtakes runner-1. In this case, at the end of the race, runner-1 has covered a distance [START_REF] Channel | Athletics. Men's 1500m Final IAAF Diamond League Rome 2014[END_REF][START_REF] Channel | Athletics Gymnastics Aquatics. Athletics -Men's 1500M -Beijing 2008 Summer Olympic Games[END_REF][START_REF] Channel | Sport Singapore. Athletics Men's 1500m Final (Day 6) | 28th SEA Games Singapore 2015[END_REF]. We want to point out that the athlete who won the Beijing 2008 Olympics was disqualified one year later for doping and his gold medal was reassigned. Nonetheless, it is still interesting to analyse the race, knowing that doping increases the maximal value of σ but delays the time at which the peak velocity is reached: this should lead to a slower start, but it provides a capacity to keep a higher velocity for a longer part of the race. From the three videos [START_REF] Channel | Athletics. Men's 1500m Final IAAF Diamond League Rome 2014[END_REF][START_REF] Channel | Athletics Gymnastics Aquatics. Athletics -Men's 1500M -Beijing 2008 Summer Olympic Games[END_REF][START_REF] Channel | Sport Singapore. Athletics Men's 1500m Final (Day 6) | 28th SEA Games Singapore 2015[END_REF], one can observe that the winner always runs the first part of the race behind, and this is consistent with our numerical results. The overtaking occurs at 84.6% of the race in Beijing 2008, at 96.9% in Rome 2014 and at 91.8% in Singapore 2015: these values are close to our numerical results, that vary between 87% and 99% of the race depending on the difference between the athletes.

A.5 Conclusion

In this work, we have presented a new model for a two-runners problem, starting from the single runner model of Aftalion and Bonnans [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF] and from the two-runners model of Pitcher [START_REF] Pitcher | Optimal strategies for a two-runner model of middle-distance running[END_REF], changing the optimal control problem. The key of our simulations is that they quantify very precisely in terms of physiological parameters, optimal control problems and numerical simulations, phenomena which are only qualitatively understood.

In this paper, we do not take into account the curvature of the track, which is the aim of an upcoming paper, since it requires more effort in the modelling. As expected, going from one runner to two runners does not change the main characteristics of the velocity profile individuated already in [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF]. We can still clearly distinguish the different phases of the race: the fast start, with maximal propulsive force and strong acceleration until the peak velocity is reached; an intermediate phase in which the propulsive force and the velocity first smoothly decrease and then increase; a final part at maximal force again, where the runner speeds up (final sprint), followed by a very short zero energy arc, in which there is a drop in force and velocity. Running behind someone allows to keep a velocity with a smaller propulsive force than the one needed when running alone; this leads to a smaller energy consumption, therefore the zero energy level occurs later and the speed up at the end of the race is more pronounced.

Our numerical results show that if a runner has a fast start and leads the race for the most part, even if they are slightly stronger than their opponent, at the end they are overtaken: in order to lead the race and win, the physiological difference between the athletes has to be significant. Furthermore, we have shown how runners can improve their personal best performance by exploiting the advantage of running behind someone else, who can be stronger or weaker. The most significant improvements are obtained by running behind someone stronger.

An interesting development, in order to have more realistic results, would be to include a delayed reaction term which takes into account the fact that runners cannot adapt instantaneously their strategy to changes in their competitor's strategy. This could be compared to a stochastic model. Finally, one could increase the number of runners, in order to be able to model real races more accurately. These considerations are outside the scope of this paper, but they can be important for future research.

This model suggests to use special runners to set the pace for others and help improve their racing times in training. The other major application for Olympic training could be for athletes to estimate whether they should stay behind or lead, and when, given their physiology, and that of their opponents, is the best time to overtake.

Figure 2 . 2 -

 22 Figure 2.2 -Scheme of the solver of Roe for the sensitivity.

Figure 2 . 3 -

 23 Figure 2.3 -Godunov scheme: σ is the slope of the blue line, i.e. the speed at which the information travels.

5 1 0 x 1 tFigure 2 . 4 -

 124 Figure 2.4 -Characteristic curves for the Burger's equation with initial condition (2.22).

Figure 2 . 5 -

 25 Figure 2.5 -Continuous initial condition.

(

  f) repeat (b)-(e) until the final time T .

Figure 2 . 15 -

 215 Figure 2.15 -Sensitivity solution for the problem with continuous i.c. with correction, a = x c .

Figure 2 . 16 -

 216 Figure 2.16 -Sensitivity solution for the problem with continuous i.c. without correction, a = x c .

Figure 2 . 17 -

 217 Figure 2.17 -Shock detector based on the second derivative.

Figure 3 . 2 -

 32 Figure 3.2 -Set of points reachable from u L through a shock (in red) or through a rarefaction (in blue).

Figure 3 . 3 -

 33 Figure 3.3 -Set of points reachable from u R through a shock (in red) or through a rarefaction (in blue).

Figure 3 . 4 -

 34 Figure 3.4 -Configurations for the state variable U.

Figure 3 . 9 -

 39 Figure 3.9 -Classical finite volume schemes.

1 Figure 3 . 11 -Figure 3 . 12 -

 1311312 Figure 3.10 -Classical finite volume schemes for sensitivities -zoom.

Figure 3 . 13 -

 313 Figure 3.13 -Definition of the temporary staggered mesh.

Figure 3 . 14 -

 314 Figure 3.14 -Roe-type schemes without numerical diffusion.

Figure 3 . 1 Figure 3 . 16 -Figure 3 . 17 -

 31316317 Figure 3.15 -Roe-type schemes without numerical diffusion for the sensitivity -zoom.

w 3 2 w 4 w 3 1 + γ w 2 w 6

 16

Figure 4 . 2 -

 42 Figure 4.2 -Structure of the solution for the Riemann problem for the sensitivity.

- 4 . 2 ,

 42 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 Figure4.1 for the state and Figure4.2 for the sensitivity. Let us now give the exact formulas for the state and for the sensitivity.

Figure 4 . 3 -

 43 Figure 4.3 -Spatial discretisation.

Figure 4 . 4 -

 44 Figure 4.4 -Structure of the Roe solver for the state.

Figure 4 . 5 -

 45 Figure 4.5 -Structure of the HLL-type solver for the sensitivity.

Figure 4 . 6 -

 46 Figure 4.6 -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.

Figure 4 . 7 -

 47 Figure 4.7 -Convergence test for the state.

1 Figure 4 . 8 -

 148 Figure 4.8 -Convergence test for the sensitivity -HLL-type scheme.

1 Figure 4 . 9 -

 149 Figure 4.9 -Convergence test for the sensitivity -HLLC-type scheme.

  0

Figure 4 . 12 -Figure 4 . 13 -

 412413 Figure 4.12 -First and second order schemes, without numerical diffusion. HLL-type scheme for the sensitivity.

Figure 4 . 14 -Figure 4 . 15 -

 414415 Figure 4.14 -Monte Carlo approach. Average and the average plus and minus twice the standard deviation in red. Five samples in black dashed lines.

Figure 4 . 16 -Figure 4 . 17 -

 416417 Figure 4.16 -SA approach without correction. Average and the average plus and minus twice the standard deviation in red. Five samples in black dashed lines -zoom.

Figure 4 . 18 -

 418 Figure 4.18 -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.

2 Figure 5 . 1 -

 251 Figure

Figure 5 . 4 -

 54 Figure 5.4 -Numerical results Mach number comparision -isentropic transonic case.

Figure 5 . 5 -Figure 5 . 6 -Figure 5 . 7 -Figure 5 . 8 -Figure 5 . 9 - 1 : 2 : 3 : 4 :: r r 0 7 :

 555657585912347 Figure 5.5 -Numerical vs empirical sensitivity for the isentropic transonic case.

a 18 : 19 :

 1819 new = a -rr a J State solution ! U Evaluation of J(U) ! J new 20:

Figure 5 . 10 -

 510 Figure 5.10 -Cost functional (5.14) for the isentropic transonic case, a =(A, `).

Fi gure 5 . 11 -Figure 5 . 12 -Figure 5 . 13 - 14 -Figure 5 . 15 -

 51151251314515 Figure5.12 -Value of the cost functional and norm of its gradient with respect to iterations in a semilogarithmic scale -isentropic transonic case, a =(A, `).

Figure 5 . 16 -Fi gure 5 . 17 -

 516517 Figure 5.16 -Cost functional (5.14) for the isentropic transonic case, a =(x c , `).

Figure 5 . 18 -Figure 5 . 19 -Fi gure 5 . 20 -

 518519520 Figure 5.18 -Value of the cost functional and norm of its gradient with respect to iterations in a semilogarithmic scale -isentropic transonic case, a =(x c , `).

Figure 5 . 21 -

 521 Figure 5.21 -Value of the cost functional and norm of its gradient with respect to iterations in a semilogarithmic scale -shocked flow case, a =(A, `).

T

  (f ) s.t (A.1)-(A.2), (A.4)

⌘

  Figure A.1). In fact, σ depends on five parameters: the initial value at rest σ r ,themaximalv alueσ max ,thefinalv alueσ f and two parameters ϕ and e cr which denote the transition point from one zone to another, ϕ,e cr 2 (0, 1): σ(e; σ max , σ f , σ r , ϕ,e cr )if e<e 0 e cr σ max if e 0 e cr  e  e 0 ϕ σ r + (σmax-σr)(e 0 -e) e 0 (1-') if e ≥ e 0 ϕ, (A.5)

2 )Figure A. 2 -

 22 Figure A.2 -Term that modulates the friction with the air in the two-runners model.

  Figure A.2.

Figure A. 7 -

 7 Figure A.7 -Competition between two runners with the same parameters, with the recreation term η; distance between the runners.

Figure A. 10 -

 10 Figure A.10 -Competition between two runners with different σ.

Figure A. 11 -Figure A. 13 -

 1113 Figure A.11 -Competition between two runners with different σ; distance between the runners.

Figure A. 14 -

 14 Figure A.14 -Competition between two runners with different ⌧ ; stronger runs the first part of the race behind.

Figure A. 15 -

 15 Figure A.15 -Competition between two runners with different ⌧ ; stronger runs the first part of the race behind; distance between the runners.

Remerciements

Acknowledgements: I would like to thank Frédéric Bonnans for his very helpful comments and A. Aftalion for suggesting this interesting topic of research and for her many remarks. A first version of this work, of which she is co-author, can be found on arxiv (http://arxiv.org/abs/1508.00523v1).

The definition of the source term comes from the derivation with respect to a of the Rankine-Hugoniot relations associated with (4.3):

-σ(U + -U -)+F(U + ) -F(U -)=0, which gives:

where the terms with @ a x s (t) are very difficult to estimate, however they are all zero in a first order finite volume framework, therefore they can be neglected. A special treatment, which will be detailed later, is necessary for a second order discretisation.

The new system can thus be written as: ( @ t U + @ x F(U)=0 @ t U a + @ x F a (U, U a )=S. (4.16) Before going to the design of a numerical scheme to approximate the solution of (4.16), we briefly discuss in the next section the typical structure of these solutions and compare it with the ones of (4.3). We specify the solution for a given initial data of Riemann type.

Riemann problem

In this section, we write the exact solution for the system (4.16) in a specific case (cf. [START_REF] Appel | Sensitivity calculations for conservation laws with application to discontinuous fluid flows[END_REF]), which will be 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 , r 2 ) is linearly degenerate, i.e. rλ 2 • r 2 =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 4.1 we show the structure of the state in the case of a rarefaction-contactshock. Concerning the sensitivity, it has the same structure as the state (cf. Figure 4.2 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 4.2 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:

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 + λ 1 U a,L -λ 3 U a,R + F a,R -F a,L = ∆xS, which can be rewritten as:

)+F a,R -F a,L = ∆xS. Replacing the definitions (4.28) one finds:

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

Therefore, we obtain:

which is consistent with our discretisation of the source term.

Second order MUSCL-type extension

In this section, we extend to the second order the schemes presented above. In time, we use a standard two-step Runge-Kutta method, whilst in space we use a MUSCLtype approach. In a few words the main idea of a MUSCL-type scheme is to consider a replacement of a constant value V n j in each cell by a higher order polynomial V n j (x), x 2 [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 (4.13)-(4.14) is valid only if the state is piecewise constant (we refer to the previous Chapter for more details). Therefore, we suggest 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.6). In this work, we compute the edge values with a standard approach:

V n j±1/4 = V n j ± ∆V n j , and the 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 account for the neglected terms in (4.15).

Numerical schemes

In this section we briefly describe the numerical schemes implemented for the state and for the sensitivity. We consider a uniform grid in space with a constant step ∆x, x j is the center of the j-th cell C j , whose extrema are x j-1/2 and x j+1/2 (cf. Figure 5.3). 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 C j at time t n .

We remark that the left-hand side of the state equation (5.4) is the classical Euler system (4.3) and the left-hand side of the sensitivity equation (5.6) is identical to the left-hand side of the sensitivity of the Euler system (4.6). Therefore, the same numerical schemes used to solve the Euler system and its sensitivity can be used here to discretise the conservative part. For more details, cf. Chapter 4.

Concerning the right-hand side of the state equation (5.4) and its corresponding part in the sensitivity equation (5.6), we propose a very simple discretisation: the term (P -F(U)) is evaluated in each cell, whilst for the derivative of h(x) we can either use the analytical derivative evaluated in the center of each cell, or a centred finite-difference discretisation.

As we did for all the systems considered in this thesis, for the conservative part of the system we use Godunov-type schemes, which consist of two main steps: first, the solution of the Riemann problem at each interface x j-1/2 at time t n ; then, the projection step in order to obtain a piecewise constant solution on the mesh:

• Riemann solvers. For the first step, we use two different approximate Riemann solvers, one for the state and one for the sensitivity. For the state, we choose the approximate Riemann solver of Roe. This choice is motivated by the fact that to compute a correction term of the form (5.7) we need to be very accurate for the state in the shock and this solver has the property of being exact for an isolated shock. However, different numerical schemes can be used, provided they have two intermediate states, as already discussed in the previous chapter. This means that the HLL solver cannot be used for the state. Concerning the sensitivity, we analysed two different approximate Riemann solvers in the previous chapter: an HLL-type scheme and an HLLC-type scheme. Numerical results show that these two schemes provide almost the same solution (the difference is negligible), the HLL-type scheme being less expensive from a computational point of view. For this reason, in this framework we adopted the HLL-type scheme for the sensitivity.

• Projection step. In a classical Godunov method, the projection step is usually per-

Optimization

Problem description

Sensitivity analysis can be used to solve PDE constrained optimization problems such as: min a2A J(U), subject to a system of PDEs, for instance (5.3). J is the cost functional and it depends on the vector of parameters a through the state U. The set A is the set of admissible parameters. The existence of the minimum is guaranteed under some hypotheses: continuity of the cost functional J and compactness of the set A. Moreover, if the cost functional is convex with respect to the parameters, the minimum is unique.

Let us consider a problem where the cost functional can be written as a bilinear form:

Classical optimization algorithms require the gradient of cost functional, whose i-th component can be written as:

[r a J(U)] i = @ a i J(U)=b(U, U a i ).

(5.13)

One of the most used techniques to compute such gradients is the adjoint equation method [Jam88, MP01, Pir74], which introduces additional adjoint variables to compute the derivative of any functional output with respect to all input parameters. The adjoint equation is independent of the input parameters, thus this approach is 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. Here, we use the continuous sensitivity equation method: one sensitivity system has to be solved for each parameter a i . However, all the sensitivity systems are independent of each other, therefore they can be solved in parallel.

Optimization algorithm

In this section we briefly describe the optimization algorithm used to obtain the results presented in the next section. It is a projected gradient descent method and it is detailed in Algorithm 1. First, the state U is solved and the cost functional J is evaluated in the initial parameters a. Secondly, all the sensitivity equations are solved (line 7) in order to compute the gradient of the cost functional r a J (5.13): we remark that this step is performed in parallel, all the sensitivities being independent from each other. Finally, a line search is performed (lines 15-20): notice that for this step, only the solution of the state equation is required. The parameters value is updated and the loop is repeated until convergence.

Test cases

In all the test cases presented in this section, we deal with pressure matching problems, i.e. the cost functional is the following:

Appendices

A | Modelling of running strategies

A.1 Abstract

In order to describe the velocity and the anaerobic energy of two runners competing against each other for middle-distance races, we present a mathematical model relying on an optimal control problem for a system of ordinary differential equations. The model is based on energy conservation and on Newton's second law: resistive forces, propulsive forces and variations in the maximal oxygen uptake are taken into account.

The interaction between the runners provides a minimum for staying one meter behind one's competitor. We perform numerical simulations and show how a runner can win a race against someone stronger by taking advantage of staying behind, or how they can improve their personal record by running behind someone else. Our simulations show when it is the best time to overtake, depending on the difference between the athletes. Finally, we compare our numerical results with real data from the men's 1500m finals of different competitions.

A.2 Introduction

The running strategy to win an Olympic medal is quite complex. It relies on outstanding physiology, good preparation, psychological factors and the optimal way to compete with the others to beat them. Quite a few mathematical works starting with Keller's [START_REF] Keller | Optimal velocity in a race[END_REF], have analysed running strategies for a single runner [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF][START_REF] Behncke | A mathematical model for the force and energetics in competitive running[END_REF][START_REF] Mathis | The effect of fatigue on running strategies[END_REF][START_REF] Morton | A three component model of human bioenergetics[END_REF][START_REF] Woodside | The optimal strategy for running a race (a mathematical model for world records from 50 m to 275 km)[END_REF], but very few take into account the competition situation where the point is to beat the others [START_REF] Kyle | Reduction of wind resistance and power output of racing cyclists and runners travelling in groups[END_REF][START_REF] Pitcher | Optimal strategies for a two-runner model of middle-distance running[END_REF]. The point of view of Keller [START_REF] Keller | Optimal velocity in a race[END_REF] is to write the equations governing the energy and the velocity of a single runner, starting from Newton's second law and energy conservation. He considers a simple problem, in which the athlete runs alone on a straight path of length D,a n dt h ea i mi st om i n i m i s et h et i m eT when the runner reaches the final distance D. This model matches the final times of world records. However some hypotheses are not physiologically reasonable, and this leads to a non-realistic velocity profile. Some authors have tried to improve Keller's model: Woodside in [START_REF] Woodside | The optimal strategy for running a race (a mathematical model for world records from 50 m to 275 km)[END_REF] and Mathis in [START_REF] Mathis | The effect of fatigue on running strategies[END_REF] introduce a correction by adding a fatigue term for long races. Ward-Smith [START_REF] Ward-Smith | A mathematical theory of running, based on the first law of thermodynamics, and its application to the performance of world-class athletes[END_REF] and Morton [START_REF] Morton | A three component model of human bioenergetics[END_REF][START_REF] Morton | A 3-parameter critical power model[END_REF][START_REF] Morton | The critical power and related whole-body bioenergetic models[END_REF] follow a different approach. Morton has introduced a three component model to take into account the variations in the oxygen uptake ( VO2) but the full optimal control problem is not solved. Behncke [START_REF] Behncke | A mathematical model for the force and energetics in competitive running[END_REF] incorporates the hydraulic model of Morton to a biomechanical model that extends the ones of Keller and Ward-Smith: it is more detailed in terms of resistive forces and takes into account the reaction time of the athlete. Aftalion and Bonnans [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF] improve the models of Keller [START_REF] Keller | Optimal velocity in a race[END_REF], Behncke [START_REF] Behncke | A mathematical model for the force and energetics in competitive running[END_REF] and Morton [START_REF] Morton | A three component model of human bioenergetics[END_REF][START_REF] Morton | A 3-parameter critical power model[END_REF][START_REF] Morton | The critical power and related whole-body bioenergetic models[END_REF] and assume that maximal oxygen uptake is a function of the anaerobic energy of the athlete. The aim is to match experimental measurements, in particular, in [START_REF] Hanon | Pacing strategy and VO2 kinetics during a 1500-m race[END_REF][START_REF] Hanon | Effects of optimal pacing strategies for 400-, 800-, and 1500-m races on the VO2 response[END_REF][START_REF] Thomas | Oxygen uptake response to an 800-m running race[END_REF], where Hanon et al. show how the oxygen uptake varies during races of 400, 800 and 1500m. They solve the full numerical control problem using an optimal control solver Bocop.

As pointed out by Pitcher [START_REF] Pitcher | Optimal strategies for a two-runner model of middle-distance running[END_REF], in middle distance running, it is common practice to try to position oneself behind but within striking distance of the leader for most of the race and then overtake them near the finish line. Pitcher explains that the runner behind can take advantage of the slipstream of the runner in front and relies on analyses of Pugh [START_REF] Pugh | The influence of wind resistance in running and walking and the mechanical efficiency of work against horizontal or vertical forces[END_REF] and Kyle [START_REF] Kyle | Reduction of wind resistance and power output of racing cyclists and runners travelling in groups[END_REF]. We believe that it is a combination of slipstream and psychological factors which explain why it is better to stay behind, and the equations can incorporate all this. The weakness of Pitcher's paper is that she imposes a strategy for one of the runners and allows only the second runner to have a free strategy. Therefore, in this paper, based on the recent work of Aftalion-Bonnans [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF], we extend the model of Pitcher [START_REF] Pitcher | Optimal strategies for a two-runner model of middle-distance running[END_REF] to include slipstream and psychological factors in a two runners race and we set a realistic optimal control problem with each runner having a free strategy.

A.2.1 Mathematical model for a single runner

The system of Keller couples together the velocity of the runner at instant t, v(t),t h e energy of the runner at instant t, e(t), and the propulsive force of the runner at instant t, f (t). The first equation is Newton's second law: it involves the propulsive force and the friction. Here, ⌧ is a constant coefficient which gathers together all the friction effects, supposed to be linear in v. The friction term can be modified to include air resistance [START_REF] Kyle | Reduction of wind resistance and power output of racing cyclists and runners travelling in groups[END_REF] which adds a term in -cv 2 to the first equation.

The second equation is an energy balance incorporating the oxygen uptake, σ, considered constant in Keller's paper, while the second term is the work of the propulsive force f . Both equations are normalised with respect to the runner's mass:

We use the dot notation to indicate the derivative with respect to time, i.e. v = dv dt and ė = de dt . We have set x(t) to be the position so that ẋ = v. The final time T is defined as the time to reach the distance D. Moreover, e 0 is the initial energy. It is necessary to add some physiological constraints to the system (A.1):

• the energy must be positive:

• the propulsive force has an upper bound which depends on the runner's physiology, and a positive lower bound due to the fact that they are moving forwards:

Therefore, in this model, the athlete is identified by four parameters: e 0 ,t h ei n i t i a l energy, ⌧ , the friction coefficient, σ, the oxygen uptake, and f M the maximal propulsive percentage reduction at the optimal distance. The parameter ↵ is an index of the variance of the phenomenon and is chosen large enough so that the exponential term becomes negligible as we deviate by half a meter from β. We observe that the choice of a non symmetric term would probably be more accurate and realistic, but more complicated and with more parameter to estimate: for these reasons, in this work we use the term suggested by Pitcher in [START_REF] Pitcher | Optimal strategies for a two-runner model of middle-distance running[END_REF], which gives reasonable results when compared to real races. The parameters ↵, β and γ do not depend on the runner, however the parameter c is related to the drag coefficient and depends on the shape and surface properties of the athlete's body (see [START_REF] Quinn | The effects of wind and altitude in the 400-m sprint[END_REF] for further details). For simplicity, in this work, we consider c 1 = c 2 = c. If ever the athletes have different masses, then we would have

The problem has physiological constraints:

In order to have the state equations for both runners defined on the same time interval, Pitcher chooses a fixed final time T . Moreover, she fixes the running strategy of runner 1 (therefore f 1 (t) is given), as the optimal strategy they would adopt if running alone, therefore the only control is f 2 (t). Formally, the resulting optimization problem is:

x D (T ) s.

t. (A.7)-(A.8). (A.9)

It is clear that this strategy is not realistic, because all the runners adapt their strategy according to the performances of their opponents. In fact, one of the crucial point in a two-runners problem is how to model competition between them: some possible ways are game theory or multi-objective optimization. However, in this work, as explained below, we model the competition by leaving both strategies free, therefore having two controls, f 1 and f 2 , and by encompassing in the cost functional the distance between the runners at the final time.

A.3 Mathematical model

In this work, we use Pitcher's two-runners model [START_REF] Pitcher | Optimal strategies for a two-runner model of middle-distance running[END_REF] and the single runner model of Aftalion-Bonnans [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF] to build a new model for two runners, which incorporate psychological factors.

In order to simplify the notation, we substitute the expression σ(e i ; σ max,i , σ f,i , σ r,i , ϕ i ,e cr,i ) from (A.5) with σ i (e i ). We recall that the friction in Pitcher's paper is -cv 2 i (1γe -↵(x D ±β) 2 ) and is supposed to model the slipstream. Shielding behind someone has a strong impact when there is wind or when the velocity is high, however in real races this position does not come entirely from slipstream but from strategic factors, too, for which the position one meter behind is the best. Therefore, the potential 1-γe -↵(x D ±β) 2 can also be considered to model a psychological factor which consists in trying to follow one's competitor, in order to be able to overtake. Indeed, it is a potential which has a minimum at distance β behind and decreases global friction. On the other hand, when the other runner is too far, there is no benefit. Let us observe that this model does not take into account the lateral displacement, and therefore the additional propulsive force, that is necessary to overtake. One can model the fact that overtaking requires some additional energy, by possibly using a non symmetric potential well 1γe -↵(x D ±β) 2 , with a varying ↵, which is not what we have done in the simulations presented to reduce the number of parameters involved. We obtain the following equations: for i =1, 2

(A.10)

The equations in (A.10) are defined for t 2 (0,T), where T is the time at which the first of the two runners reaches the final distance D; to model this, it is necessary to add the following boundary condition to the system:

As in the previous models, the energy has a lower bound:

The choice of the cost functional, i.e. the quantity to be minimised, is a key point.

As said before, in contrast with Pitcher's choice, in this case none of the strategies is fixed, therefore there are two controls: f 1 (t) and f 2 (t). Here, we propose to minimise the following quantity, given a proper constant weight c w > 0:

The aim of this choice is to minimise the final time of the winner, and the term c w |x D (T )| models the fact that the loser has tried to win as well. Different values of c w can lead to different results, as in real races when two runners compete against each other multiple times the outcome of the race can change.

The resulting problem is:

where F i is the set of the admissible controls, and depends on the athlete. For physiological reasons, it is necessary to impose a bound to the variations of ḟ ,i na d d i t i o n to the bounds on f already introduced, related to the fact that athletes cannot vary their propulsive force too quickly (see more details in [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF]). This leads to the following definition of F i :

where K i and f M,i are constants depending on the athlete, which model the fact that every runner has a limited maximal force (f M,i ) and cannot vary it too quickly (K i ).

The rest of the paper consists in providing numerical simulations of (A.14).

A.4 Numerical results

All the results presented in this section are obtained with the free software BOCOP [START_REF] Bonnans | BOCOP -A toolbox for optimal control problems[END_REF]. The equations are solved with a finite difference scheme (implicit Euler), while the optimization problem is solved with an iterative method, using as stopping criterion the difference between successive iterates, with a tolerance of 10 -10 . The aim of these simulations is to find out if a runner can win a race against someone stronger, by running behind the first part of the race, and to quantify in term of variation of some parameters how much weaker they can be and when the best time to overtake is.

The single runner strategy has been mathematically proven and numerically computed in [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF], and found experimentally in [START_REF] Hanon | Pacing strategy and VO2 kinetics during a 1500-m race[END_REF], and it consists in three parts:

• a first part of maximal force with a strong acceleration during which the peak velocity is achieved,

• a second part in which the propulsive force first decreases smoothly and then increases again, with the corresponding decrease and increase of the velocity,

• a final part at maximal force and maximal velocity again, until zero energy level is reached, where the velocity drops.

From the two-runners model, we expect an overall similar strategy: however, the additional term in the velocities equations encourages one of the runners to start slightly slower and to position themselves at distance β from the other. This allows them to keep the same velocity as the other runner while using a smaller propulsive force, which, in turn, leads to a lower energy consumption. We expect that, at a certain point during the race, the runner who is behind, will overtake the other by using the energy they have saved throughout the race and will be able to perform a longer final sprint. Moreover, it is reasonable to think that this moment occurs sooner if the runner who runs the first part of the race behind is stronger. Let us observe that this running strategy could also lead, in some favourable situations, to an improvement of the personal record. Nonetheless, it is important to underline that the decrease in the final time is not the main goal for a runner in some occasions, such as the Olympic finals, in which the final position is definitely more important than the final time: in [START_REF] Thiel | Pacing in olympic track races: competitive tactics versus best performance strategy[END_REF], Thiel et al. study, starting from the Beijing 2008 Olympic Games data and the world records data, how the difference in the goal affects the pacing strategies: win the race versus minimising the final time.

How to estimate the parameters values starting from a race is beyond the scope of this paper. For this reason, the reference values for the parameters used in the simulations are taken from the literature and reported in Table A.1. The initial, maximal and final value of sigma (respectively σ r , σ max and σ f ) are taken from the VO 2 values reported in [START_REF] Hanon | Pacing strategy and VO2 kinetics during a 1500-m race[END_REF]: let us observe that these values are given in ml kg -1 min -1 . In order to convert them in the unity of measurement needed (i.e. m 2 s -3 = Jk g -1 s -1 ), we consider that the uptake of 1ml of oxygen is often converted into an energy expenditure estimate of 21J. It is then necessary to convert the minutes in seconds, obtaining in this way the conversion factor 21/60. The values chosen for ϕ and e cr aim at fitting the 2 profile reported in Figure 2 of [START_REF] Hanon | Pacing strategy and VO2 kinetics during a 1500-m race[END_REF]. The values for f M and τ are strongly related: in fact, a first order approximation of the maximal velocity v peak a runner can reach is τ f M . Therefore, starting from the velocity profile plotted in Figure 1 of [START_REF] Hanon | Pacing strategy and VO2 kinetics during a 1500-m race[END_REF], one can chose a couple of reasonable values. The value of the constant c is taken from [START_REF] Quinn | The effects of wind and altitude in the 400-m sprint[END_REF], c ⌘ from [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF], ↵ and β from [START_REF] Pitcher | Optimal strategies for a two-runner model of middle-distance running[END_REF] and, finally, γ from [START_REF] Pugh | The influence of wind resistance in running and walking and the mechanical efficiency of work against horizontal or vertical forces[END_REF].

First of all, let us consider the perfectly symmetric situation, in which the two runners have the same parameters. In this case it is very influential which runner is running the first part of the race behind: this can be mathematically modelled by choosing, in a proper way, the initial guess for the variable x D given to the iterative method that solves the optimization problem. In fact, giving as initial guess β (or -β) forces runner-1 (or runner-2) to start the race more slowly and to position themselves behind for the first part. If one gives as initial guess x D =0, one finds a solution in which the runners overtake each other multiple times, which is not very realistic. The non uniqueness of the solution is not surprising, especially in a perfectly symmetric situation, such as the one considered: if a certain couple of strategies {f 1 (t),f 2 (t)} provides a minimum for the cost functional, the couple {f 2 (t),f 1 (t)} provides a minimum, too. This can be considered to model the fact that if the same two runners run against each other multiple times, the outcome of the race can change. The results of this simulation are shown in Figures A. 3-A.4: one can observe that the runner who stays behind can keep the same velocity as the other runner, while using a significantly lower propulsive force and therefore having a much lower energy consumption. All the graphs in Figure A.3 have the position on the x-axis, which means that the velocity of runner-i is plotted with respect to the position of runner-i (and it is the same for the propulsive force, the energy and sigma). The choice of plotting with respect to position, and not time, has been made because it is the most common in the sports literature. Figure A.4 shows the distance between the runners: the overtaking occurs at about 94% of the race, corresponding to 1416m. This variable is plotted with respect to a normalized time (i.e. t/T ), because it does not concern only one runner, but both of them, therefore it does not make sense to plot it with respect to the position of any of them. From this first result, it is clear that running behind someone for the most part of the race allows runners to win against athletes as strong as themselves. We now want to investigate how the strategy of runners changes if they are running alone or behind someone else. In Figure A.5, two performances of the same runner are compared: the blue line represents the optimal strategy of the runner running alone (they complete the race in 249.681s), the red line represents the strategy adopted when running behind someone as strong as themselves (they completes the race in 247.822s). This difference in the final time (almost 2s of improvement) is equivalent to a difference of about 0.05m/s in the mean velocity. However, what we have simulated is not the most favourable situation to reduce the final time, and therefore to improve the mean velocity. If the aim is exclusively to improve the personal best performance and not to win the race, the best scenario possible for a runner is to run behind someone slightly stronger for the whole race. However, let us observe that the opponent must not be too much stronger: in this case, the runner would use too much energy to stay behind and would soon reach the zero energy level, which would cause a drop in the propulsive force too soon in the race. Nevertheless, from can be found also in the single-runner problem, see [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF] for further details), and the mean value around which the propulsive force of runner-1 is oscillating is slightly bigger than the one around which the force of runner-2 is oscillating; the energy curve does not change significantly, compared to the previous results. From Figure A.7, one can notice that runner-2 overtakes at about 94.85% of the race, which corresponds to 1417m: one meter later, if compared with the case without oscillations. Let us observe that in this case the final time is smaller: in fact runner-2 completes the race in 247.75s. This decrease in the final time, when adding the recreation η, is consistent with the results for the single-runner problem presented in [START_REF] Aftalion | Optimization of running strategies based on anaerobic energy and variations of velocity[END_REF].

Finally, we can say that the recreation term does not change the strategy of the race, however the results are more difficult to read, due to the oscillations. For this reason, from now on we will present only results obtained without η.

We now want to find the threshold values, i.e. how much runners can be weaker than their opponent and still win the race by running behind. Therefore, we vary one parameter at a time, making runner-2 weaker. Figures A.8-A.9 show the results for two runners who have a different initial energy. Given the reference value for e 0 1 = 1400 (as in Table A.1), the lowest initial energy runner-2 can have, while still being able to win the race, is:

At the beginning of the race, runner-2 slows down in order to stay behind: in this way runner-2 manages to keep the same velocity as runner-1 using a smaller force; this leads to a smaller energy consumption, therefore at the end of the race runner-2 has enough energy to speed up and overtake runner-1. For the boundary condition (A.11), the time stops as soon as the first runner finishes the race. Therefore, when runner-2 reaches the of 1498.82m. The final time of runner-2 is 249.536s, while their best performance running alone is 251.83s, i.e. they have an improvement of more than 2s.

Finally, let us consider a case in which the runner who starts behind is stronger. For this purpose, we use the following parameters:

All the other parameters remain unvaried (reference values from Table A.1). Figures A. 14-A.15 show the results obtained: in this case, the difference between the athletes is much bigger than in the previous ones, and this is evident from all the curves in Figure A.14. From Figure A.15, one can notice that the overtaking occurs very early in the race: at about 87.1%, i.e. 1290m. What is interesting in this case is that runner-2 completes the race in 248.726s, which is almost 1s less than their best performance running alone: therefore, in order to improve a personal record, it is not necessary to run behind someone stronger.

Let us now compare the results obtained here with Pitcher's ones. In [START_REF] Pitcher | Optimal strategies for a two-runner model of middle-distance running[END_REF], Figure 5.2, when the weaker runner is the one with the fixed strategy, the stronger runner remains only slightly ahead of their opponent until nearly the end of the race. This is in order that the weaker runner does not gain the advantage of running in the slipstream of the stronger for a long part of the race. However, it is the runner who stays behind who should adjust their position with respect to the other one, and not the opposite. The advantage of having two strategies free allow us to avoid this unrealistic result, and in this case, we get that the weaker runner stays one meter behind, and either wins the race if the difference in energy is not too big or drops following if they do not have enough energy.

Finally, we want to analyse the strategy and to see when the overtaking occurs in real races and compare them with our results. For this purpose, we have considered