Dans cette thèse, l'apport le plus important a consisté en l'implémentation d'algorithmes modernes adaptés aux architectures massivement parallèles, dans un logiciel industriel dédié aux études de sûreté nucléaire, le code Cathare. Ce logiciel est dédié à la simulation des écoulements diphasiques au sein d'un réacteur nucléaire en conditions nominales ou accidentelles. L'implémentation de ces nouvelles techniques représentent en soi une contribution importante dans la physique des réacteurs car il permettra de déterminer, avec un temps de calcul réduit et de façon précise, l'état du coeur au cours d'accidents graves. Un effort particulier a été mené pour paralléliser de manière efficace la variable temporelle par l'algorithme pararéel. Pour cela, nous avons proposé une méthode pararéelle qui intègre de façon plus optimisée la présence de schémas en temps multi-pas. En effet, cette famille de schémas permet d'obtenir une approximation d'ordre supérieur à celui d'un schéma en temps à un pas. Cependant l'initialisation de la propagation en temps en chaque fenêtre doit être choisie avec soin. L'idée principale de ce nouveau schéma est de définir une approximation consistante des solutions permettant l'initialisation des propagations en temps, permettant ainsi à l'algorithme de converger vers la solution avec la précision voulue. Cette méthode a ensuite été appliquée sur deux cas tests représentatifs des défis numériques rencontrés dans la simulation des écoulements diphasiques dans le cadre des études de sûreté nucléaire.

La seconde partie de cette thèse est consacrée au développement de méthodes numériques permettant de traiter les difficultés numériques spécifiques aux modèles diphasiques avec un temps de calcul réduit. Dans cette partie, on développe un cadre d'analyse rigoureux pour l'étude des schémas volumes finis sur grille décalée comme celui utilisé dans le code Cathare. Les schémas décalés sont en pratique plus précis pour les fluides quasi incompressibles et sont couramment utilisés dans la communauté thermohydraulique. Cependant, pour les fluides compressibles, les études de stabilité ont été historiquement menées par une approche heuristique et par le réglage de paramètres numériques. Cette question est abordée par l'analyse des opérateurs de diffusion numérique qui permettent de porter un nouveau regard sur les schémas décalés. Cela nous permet de montrer que les schémas décalés classiques sont linéairement stables L 2 uniquement lorsque les vitesses sont de signe constant. On propose une classe de schémas décalés linéairement stables L 2 ainsi qu'une classe de schémas décalés entropiques. Ces nouvelles classes sont construites à l'aide d'un opérateur de diffusion numérique particulier et sont mieux adaptées aux modèles diphasiques pour lesquels les vitesses phasiques changent fréquemment de signe. Ces méthodes ont été appliquées au système d'Euler isentropique sur des cas tests analytiques et nous pensons que les développements actuels permettront à l'avenir son utilisation dans des cas plus réalistes et complexes, comme la simulation des écoulements diphasiques au sein d'une installation nucléaire.

A second phase of our work has been to explore numerical methods that could handle better the numerical difficulties that are specific to two-phase flows with a lower computational cost. This part of the thesis has been devoted to the understanding of the theoretical properties of finite volume schemes on staggered grids such as the one used in the Cathare code. Staggered schemes are known to be more precise for almost incompressible flows in practice and are very popular in the thermal hydraulics community. However, in the context of compressible flows, their stability analysis has historically been performed with a heuristic approach and the tuning of numerical parameters. This question has been addressed by analysing their numerical diffusion operator that gives new insight into these schemes. For classical staggered schemes, the stability is obtained only in the case of constant sign velocities. We propose a class of linearly L 2 -stable staggered schemes and a class of entropic staggered schemes. These new classes are based on a carefully chosen numerical diffusion operator and are more adapted to two-phase flows where phasic velocities frequently change signs. These methods have been successfully applied in analytical cases (involving Euler equations) and we expect that the present developments will allow its use in more realistic and complex cases in the future, like the one of the simulation of two-phase flows within a nuclear reactor during an accidental scenario.

List of Figures

Introduction (Version française) Contexte industriel

On appelle "accident grave" ou "accident de fusion du coeur" d'un réacteur nucléaire à eau pressurisé un accident au cours duquel le combustible du réacteur est significativement dégradé avec fusion plus ou moins étendue du coeur du réacteur. La fusion résulterait d'une absence prolongée de refroidissement du coeur par le fluide caloporteur et consécutivement à une augmentation importante de la température des crayons combustibles dénoyés. C'est un type d'accident qui, en raison des mesures de prévention mises en place, ne peut survenir qu'à la suite d'une accumulation de dysfonctionnements (défaillances multiples, humaines ou matérielles).

Si la dégradation du coeur ne peut pas être arrêtée dans la cuve du réacteur par refroidissement du coeur dégradé (renoyage dans la cuve par le fluide caloporteur), l'accident peut à terme conduire à une perte de l'intégrité du confinement et à des relâchements importants de substances radioactives dans l'environnement. En raison des répercussions importantes qu'aurait un tel rejet, des efforts importants sont consacrés à l'étude de ce type de scénario pour pouvoir en limiter les conséquences (approche déterministe). L'étude des accidents de fusion du coeur passe en premier lieu par l'identification des principaux scénarios pouvant conduire à ce type d'accident. En complément des études déterministes, des études probabilistes de sûreté sont également menées. La méthode consiste à analyser de façon exhaustive tous les scénarios accidentels envisageables, d'estimer, souvent en les regroupant par famille, leur probabilité d'occurrence et les conséquences associées à l'intérieur de l'installation (fusion du coeur) ou à l'extérieur (rejets radioactifs dans l'environnement).

Dans le domaine des accidents graves, les phénomènes physiques mis en jeu sont extrêmement complexes. Les objectifs de la recherche sont donc de parvenir à comprendre au mieux ces phénomènes physiques et de développer des modèles applicables aux réacteurs. Ces modèles, regroupés au sein de codes de calcul informatiques, doivent permettre de prévoir le déroulement d'un accident grave. Comme il est impossible d'effectuer, dans ce domaine, des essais à taille réelle et de reproduire toutes les situations envisageables, il est nécessaire de réaliser des essais élémentaires, permettant d'étudier séparément chaque phénomène physique. Le tout doit se faire à des échelles compatibles avec les capacités techniques et économiques des installations, tout en restant représentatives pour l'extrapolation à l'échelle du réacteur. Les réacteurs expérimentaux constituent notamment des équipements privilégiés pour étudier le comportement des combustibles nucléaires en régime accidentel.

Dans ce contexte, le logiciel Cathare est un code thermohydraulique décrivant le réacteur nucléaire à l'échelle système, développé par le CEA depuis 1979. Ce code est dédié aux études de sûreté pour les Réacteurs à Eau sous Pression et la validation des prodédures post-accidentelles. L'objectif de cette introduction n'est pas d'effectuer une liste exhaustive de tous les scénarios accidentels faisant l'objet d'études de sûreté mais de décrire un type de scénario appelé APRP (pour Accidents de Perte de Réfrigérant Primaire) afin d'illustrer de manière synthétique la démarche alliant campagnes expérimentales et codes de simulation numérique.

Accidents de perte de réfrigérant primaire (APRP): L'événement initiateur de ces accidents est une brèche dans la paroi du circuit primaire, [START_REF]Etudes et approches de la gestion des accidents graves pour les réacteurs à eau sous pression du parc français[END_REF]. La brèche provoque une fuite de réfrigérant primaire et une dépressurisation du circuit primaire. Plusieurs scénarios sont à distinguer selon l'état initial du réacteur, l'emplacement et la taille de la brèche. En cas d'APRP, la dépressurisation du circuit primaire entraîne l'arrêt automatique du réacteur, puis le démarrage automatique de l'injection de sécurité.

Les fonctions à assurer par les systèmes de protection et de sauvegarde pour limiter les conséquences de l'accident sont les suivantes: la maîtrise de la réactivité le maintien de l'inventaire en eau dans la cuve du réacteur l'évacuation de la puissance résiduelle dégagée par le combustible La maîtrise de la réactivité est assurée par l'arrêt automatique du réacteur et l'injection d'eau borée dans le coeur. Le maintien de l'inventaire en eau dans la cuve du réacteur est assuré par le système d'injection de sécurité. L'évacuation de la puissance résiduelle dégagée par le combustible est assurée par le refroidissement de l'eau circulant dans la cuve. Les scenarii accidentels menant à une fusion du coeur supposent la défaillance de l'un ou de plusieurs des systèmes de sauvegarde et sont toujours associés à une défaillance du maintien d'eau en quantité suffisante dans le circuit primaire pour refroidir le coeur.

Le réacteur de recherche PH ÉBUS [START_REF]Etudier les situations accidentelles[END_REF] est un réacteur expérimental construit en 1977 sur le Centre d'études de Cadarache. Il a été conçu pour étudier le comportement des combustibles des centrales nucléaires dans des situations accidentelles du type perte de réfrigérant primaire pouvant aller jusqu'à la fusion du combustible. La principale problématique associée à l'étude des accidents de perte de réfrigérant est celle de la dégradation du combustible et de ses conséquences: à partir de quelle température, au bout de combien de temps a-t-on rupture de la gaine du combustible ou pire, fusion du coeur? Quel est le relâchement de produits de fission associé à ces deux phénomènes? Le réacteur PH ÉBUS entre dans la catégorie des réacteurs d'essais en sûreté. L'objectif de ce programme était l'étude du comportement du combustible des réacteurs à eau sous pression (REP) dans des situations de perte de réfrigérant primaire correspondant à une situation accidentelle faisant suite à un fonctionnement en conditions nominales. Cet accident était suivi de la mise en oeuvre du refroidissement de secours. La phénoménologie étudiée était liée à l'accident de référence des REP, qui ne va pas jusqu'à la fusion du coeur. Deux objectifs étaient recherchés: évaluer les marges relatives aux deux principaux critères retenus, dans le cadre du dimensionnement du réacteur: la température maximale et l'oxydation maximale des gaines valider les codes de comportement du combustible utilisés par l'analyse de sûreté et, en particulier, le module combustible du code de calcul Cathare.

De même, le logiciel Cathare a été validé sur 20 scénarii accidentels basés sur 10 installations expérimentales.

Par ailleurs, le code Cathare est également utilisé au sein de simulateurs de réacteurs. Le simulateur est notamment l'outil de formation au fonctionnement normal et accidentel des réacteurs à eau pressurisée pour des ingénieurs de l'IRSN ainsi que pour des opérateurs. Le rôle de l'opérateur au sein d'une centrale est d'activer les procédures de conduite post-accidentelles et de maintenir le réacteur en conditions de fonctionnement nominales. Parmi les moyens disponibles pour effectuer des études de sûreté des réacteurs, SOFIA -Simulateur d'Observation du Fonctionnement Incidentel et Accidentel [START_REF]un simulateur pour améliorer la sûreté des réacteurs à eau sous pression[END_REF] est un système informatique permettant le calcul et le suivi en temps réel de l'évolution des paramètres physiques d'un réacteur nucléaire. Il permet de simuler les défaillances de matériel et les actions des opérateurs, d'arrêter le calcul pour examiner l'état de l'installation à un instant donné et de revenir en arrière, afin de modifier le scénario. Le simulateur SOFIA utilise le code de calcul Cathare en temps réel.

Objectifs de la thèse

Le code Cathare est un code à l'échelle système et décrit l'ensemble d'un réacteur nucléaire à l'aide d'un assemblage de conduites, de cuves et de pompes dont la taille est d'environ 10 mètres (taille de la cuve: 13m × 5m × 5m). Pour simuler l'ensemble de l'installation nucléaire, le code Cathare possèdent des restrictions au niveau de la taille des cellules du maillage: en prenant en compte la taille importante des circuits, la taille d'une maille peut être relativement grande, et peut aller de quelque centimètres à un mètre. Après discrétisation des variables, les simulations numériques font intervenir entre 10 2 et 10 3 inconnues et jusqu'à un million de pas de temps. Une méthode de décomposition de domaine en espace est implémentée dans le code Cathare et les performances atteignent actuellement un plateau malgré qu'il y ait des ressources informatiques supplémentaires. Les performances de cette méthode de décomposition de domaine en espace sont limitées car les maillages utilisés dans les simulations du code Cathare sont peu raffinés afin de représenter l'ensemble du réacteur nucléaire. Le premier objectif de cette thèse est de proposer une nouvelle stratégie de parallélisation, complémentaire à la méthode de décomposition de domaine en espace. Pour cette raison, on propose d'élaborer une méthode de décomposition de domaine en temps. De plus, on souhaite appliquer cette stratégie de manière non intrusive et utiliser le code Cathare en boîte noire. Pour cela, on se base sur l'algorithme pararéel. Les résultats à ce sujet sont présentés dans les chapitres 2 et 3.

Le restant des chapitres de cette thèse est motivé par le besoin en méthodes numériques nouvelles pour mieux traiter les difficultés numériques spécifiques aux modèles diphasiques avec un coût en temps de calcul raisonnable. Cette seconde partie est dédiée à l'analyse des aspects théoriques des schémas volumes finis sur grille décalée, comme celui utilisé dans le code Cathare. L'objectif est de développer une méthode d'analyse de stabilité rigoureuse des schémas décalés classiques et de proposer une nouvelle classe de schémas décalés stable L 2 . De même, on souhaite développer une nouvelle classe de schémas décalés entropiques. La méthodologie développée dans le manuscrit est très générale et pourrait s'appliquer aux modèles diphasiques.

Dans les paragraphes suivants, nous présentons un bref résumé de chaque chapitre de ce manuscrit.

Résumé des résultats par chapitres

Partie I: Chapitre 1

Ce premier chapitre a essentiellement pour but de résumer les connaissances actuelles que l'on peut trouver dans la bibliographie au sujet du modèle diphasique utilisé dans le code Cathare: le modèle bifluide à 6 équations. La plupart de ce qui est donc présenté n'est pas nouveau, mais il nous a semblé intéressant de présenter cette compilation d'informations pour positionner le modèle et les méthodes numériques du code Cathare dans la littérature sur la simulation des écoulements diphasiques.

Après la présentation du modèle bifluide à 6 équations, nous rappellerons les principaux résultats théoriques concernant l'hyperbolicité du modèle. Nous présenterons ensuite la structure complexe de la solution du modèle qui est spécifique aux modèles diphasiques.

Les techniques de discrétisation des variables du modèle seront ensuite présentées en insistant particulièrement sur le traitement du terme de convection par les schémas décalés les plus répandus. Bien que les modèles diphasiques héritent des connaissances actuelles sur la modélisation et les méthodes numériques pour les fluides monophasiques, ils possèdent néanmoins plusieurs spécificités. Nous présenterons certaines d'entre elles avec notamment le traitement des produits non conservatifs, la configuration des phases évanescentes et le traitement des termes sources discontinus.

La dernière partie du chapitre est consacrée aux techniques d'accélération actuellement déployées dans le code Cathare. On présentera également les facteurs qui limitent les performances de ces méthodes malgré la disponibilité de ressources informatiques supplémentaires. On présente ensuite la stratégie de parallélisation en temps choisie pour le code Cathare. La parallélisation de la variable temporelle est particulièrement délicate étant donné que le temps est séquentiel par nature. Malgré cela, plusieurs stratégies ont été proposées à ce sujet-là dans la littérature (see [START_REF] Chartier | A parallel shooting technique for solving dissipative ODE's[END_REF], [START_REF] Gander | 50 years of time parallel time integration[END_REF]). Nous nous sommes concentrés sur la méthode pararéelle car c'est celle qui donne les meilleures performances sur des applications complexes (voir notamment [START_REF] Baffico | Parallel-in-time molecular-dynamics simulations[END_REF], [START_REF] Fischer | A parareal in time semi-implicit approximation of the Navier-Stokes equations[END_REF], [START_REF] Samaddar | Parallelization in time of numerical simulations of fully-developed plasma turbulence using the parareal algorithm[END_REF]).

Partie I: Chapitre 2

Le deuxième chapitre résume les moyens mis en oeuvre pour implémenter l'algorithme pararéel au code Cathare de manière non intrusive. Nous présenterons les deux outils que nous avons développé durant la thèse pour appliquer l'algorithme pararéel au code Cathare: d'une part, une maquette du code Cathare restreinte à un cas test et d'autre part une librairie qui utilise le code Cathare en boîte noire de manière parallèle. Les deux cas tests que nous avons étudiés sont représentatifs des défis numériques rencontrés dans la simulation des écoulements diphasiques dans le cadre des études de sûreté. Ces défis numériques comprennent notamment les phases évanescentes pour lesquelles une des phases liquide ou vapeur disparaît dans une partie du domaine ou encore la simulation d'une brèche dans le circuit primaire d'une réacteur causant ainsi une dépressurisation rapide dans le système.

La contribution principale de ce travail est l'adaptation de l'algorithme pararéel à l'architecture logicielle du code Cathare et sa discrétisation en temps de manière non intrusive, sans modifier les fichiers sources du code, dans le but de réduire le temps de calcul et de se rapprocher d'une réponse en temps réel du code.

Finalement, dans les exemples numériques que nous avons traité, l'utilisation de l'algorithme pararéel peut accélérer les calculs d'environ un facteur 3 avec 25 processeurs. Du point de vue de l'efficacité, ces résultats ne sont pas aussi compétitifs que les méthodes de décomposition de domaine en espace fournies, mais comme il sera expliqué au chapitre 2, il existe des raisons théoriques qui expliquent la relativement basse efficacité de la méthode pararéelle. Pour cette raison, cette méthode devient intéressante pour atteindre des performances additionnelles dans un contexte où les autres techniques de parallélisation dont on peut disposer atteignent saturation.

Partie I: Chapitre 3

Comme il sera expliqué aux chapitres 1 et 2, la discrétisation en temps dans le code Cathare est basée sur un schéma en temps à deux pas. Un schéma en temps multi-pas permet d'obtenir une approximation d'ordre supérieur à celui d'un schéma en temps à un pas cependant l'initialisation de la propagation en temps en chaque fenêtre doit être définie avec rigueur. Lorsque le solveur fin et/ou le solveur grossier est un schéma en temps multi-pas, il est nécessaire de définir une approximation consistante des solutions intervenant dans l'initialisation du schéma fin pour chaque fenêtre en temps. Autrement, l'erreur commise à l'initialisation sera propagée sur l'intervalle de temps et empêchera l'algorithme pararéel de converger vers la solution avec la précision souhaitée. Dans l'objectif d'aborder ce problème, nous présentons dans ce chapitre une nouvelle variante de l'algorithme pararéel adaptée à ce type de discrétisation et qui permet de converger vers la solution cible avec un taux de convergence similaire à celui de l'algorithme pararéel classique. Un effort particulier a été réalisé afin de construire un algorithme adapté aux schémas multi-pas de manière non intrusive dans les solveurs grossier et fin. Cela permet au code Cathare d'être utilisé en boîte noire, assurant ainsi la portabilité de ce nouvel algorithme. Concernant la méthode d'initialisation, l'algorithme pararéel multi-pas est plus consistant avec le schéma sous-jacent. Nous montrons à l'aide de résultats théoriques et numériques que les propriétés de précision et de convergence de l'algorithme pararaéel multi-pas sont compétitives lorsque l'on initialise rigoureusement chaque fenêtre en temps.

Part II: Chapitre 4

La seconde partie de cette thèse est consacrée au développement d'un cadre d'analyse rigoureux pour l'étude des schémas volumes finis sur grille décalée comme celui utilisé dans le code Cathare. En particulier, la méthode présentée pourrait être appliquée dans le futur au code Cathare pour traiter les spécificités numériques propres aux modèles diphasiques. Pour développer un tel outil, il a été nécessaire tout d'abord d'étudier préalablement certains aspects théoriques et c'est ce qui est présenté dans la deuxième partie de ce manuscrit. Plusieurs exemples numériques simples seront aussi présentés, notamment sur la résolution du système d'Euler isentropique, dans le but d'illustrer la technique proposée ainsi que ses performances.

Les schémas décalés sont en pratique plus précis pour les fluides quasi incompressibles et sont couramment utilisés dans la communauté thermohydraulique ( [START_REF]A numerical fluid dynamics calculation method for all flow speeds[END_REF][START_REF] Prosperetti | Computational methods for multiphase flow[END_REF][START_REF]Discrétisation des équations du module 3d de cathare 2 & comparaison avec le module 1d -etat des lieux pour la version v2.5 2 et perspectives[END_REF]). Cependant, dans le contexte des fluides compressibles, les études de stabilité ont été historiquement menées par une approche heuristique et par le réglage de paramètres numériques ( [START_REF] Hirt | Heuristic stability theory for finite difference equations[END_REF]). Dans [START_REF]On some implicit and semi-implicit staggered schemes for the shallow water and euler equations[END_REF][START_REF] Herbin | Staggered schemes for all speed flows[END_REF][START_REF] Herbin | A kinetic energy preserving convection operator for the mac discretization of compressible navier-stokes equations[END_REF], les auteurs construisent des schémas décalés conservatifs avec des preuves rigoureuse de stabilité: le caractère entropique et l'inégalité discrète de l'énergie cinétique. Néanmoins, le caractère bornée de l'entropie n'implique pas nécessairement que la solution reste bornée et c'est particulièrement le cas pour le système d'Euler complet. Par conséquent, on étudie dans ce chapitre la stabilité L 2 linéaire des schémas décalés.

Cette question est abordée par l'analyse des opérateurs de diffusion numérique qui permet de porter un nouveau regard sur les schémas décalés. On développe d'abord la forme de l'opérateur de diffusion numérique pour les schémas décalés classiques, ce qui permet de montrer que ces schémas sont linéairement stables L 2 uniquement lorsque les vitesses sont de signe constant. On propose ensuite une classe de schémas décalés linéairement stables L 2 . Cette nouvelle classe est construite à l'aide d'un opérateur de diffusion numérique particulier et est mieux adaptée aux modèles diphasiques pour lesquels les vitesses phasiques changent fréquemment de signe.

Un schéma numérique appartenant à cette nouvelle classe de schémas décalés a été implémenté avec succès pour la simulation d'un problème de Riemann et nous pensons que les développements actuels permettront à l'avenir son utilisation dans des cas plus réalistes et complexes.

Introduction (English version) Industrial context

We call "serious accident" or "meltdown accident" of a pressurised water reactor an accident where the fuel of the reactor is significantly degraded with a more or less extensive meltdown of the reactor core. This meltdown would be the consequence of a prolonged absence of the coolant in the core causing an increase of the dewater fuel rods temperature. Due to the actual prevention measures, this type of accidents may only occur after an accumulation of malfunctions (multiple failures, human or material).

If the damages can not be stopped in the vessel by the cooling of the damaged core (reflooding of the vessel by the coolant) then the accident could lead to a loss of containment and to significant releases of radioactive substances in the environment. Due to the major consequences of such a sequence of events, several studies of these accidents are conducted to limit their consequences (deterministic approach). These studies firstly start by the identification of the main scenarii that can lead to the reactor core meltdown. Complementing the deterministic studies, probabilistic safety studies are also conducted. They consist in analysing exhaustively every possible accidental scenarii, in estimating its probability of occurence and their consequences inside the nuclear power plant (core meltdown) or outside (release of radioactive substances in the environment).

The physical phenomena encountered in nuclear safety are extremely complex. This research field is dedicated to the understanding of these phenomena and to develop physical and mathematical models applicable to reactors. These models are grouped within softwares for numerical simulation and allow to predict the course of a serious accident. Since it is impossible in nuclear safety to reproduce all possible accidental situations on real size experimental installations, it is necessary to perform several elementary experiments allowing to study seperately each physical phenomenon. These steps must be at scales that are compatible with the technical and economical constraints while remaining close to the behavior of the nuclear power plant. The experimental reactors constitute in particular a privileged equipment to study the behavior of nuclear fuel in accidental situation and represent a considerable economic and human investment.

In this context, the Cathare code (Code for Analysis of THermalhydraulics during Accident and for Reactor safety Evaluation) is a thermalhydraulic code describing a nuclear reactor at the system scale, developed by CEA since 1979 as part of an agreement between CEA, EDF, AREVA and IRSN. This software is dedicated to the safety studies of Pressurised Water Reactor and the validation of emergency procedures during an accidental scenario.

The objective of this introduction is not to give an exhaustive account of all the possible accidental scenarii studied in nuclear safety but to describe a type of accidents called LOCA (Loss Of Coolant Accident) to illustrate synthetically the approach that couples experimental campaigns and numerical simulation softwares.

Loss Of Coolant Accident (LOCA):

The initiating event of these accidents is a breach in the primary circuit, [START_REF]Etudes et approches de la gestion des accidents graves pour les réacteurs à eau sous pression du parc français[END_REF]. This breach causes a leak of the coolant and a depressurisation of the primary circuit. Many accidental scenarii are then possible according to the initial state of the reactor, the location and the size of the breach. In the case of a LOCA, the depressurisation in the primary circuit will generate an automatic shutdown of the reactor and an automatic start of backup cooling system.

To limit the consequences of these accidents, protection and backup systems will ensure the following functions: control of the reactivity maintaining the water inventory in the reactor vessel evacuate the residual heat released by the fuel rods The control of the reactivity is ensured by the automatic shutdown of the reactor and by the injection of borated water within the core. The water inventory in the vessel is maintained by the safety injection of water. The residual heat released by the fuel rods is evacuated by the cooling of the water circulating in the vessel. Accidental scenarii leading to a core meltdown are the consequences of the failure of one or several backup systems and are always associated to a failure in maintaining enough coolant in the primary circuit to cool the reactor core.

The research reactor PHEBUS [START_REF]Etudier les situations accidentelles[END_REF] is an experimental reactor built in 1977 in the CEA research Center of Cadarache. It is dedicated to the study of the fuel rods behavior in nuclear power plants under accidental conditions of LOCA type. The main issue associated to this type of accidents is the fuel rods degradation and its consequences: from which temperature, after how long do we have the break of the fuel rod cladding or worse core meltdown? How much fission products are released due to these two phenomena? The PHEBUS reactor is an experimental reactor for safety studies. The objective of this program is to study the behavior of fuel rods for pressurised water reactors under LOCA type accidents especially. There are two main goals: evaluate the relative margins for two selected criteria, in order to sizing the reactor: the maximum temperature and the maximum oxidation of the fuel rod cladding validate the softwares used in the nuclear safety studies, particularly the fuel module of the Cathare code.

Likewise, the Cathare code has been validated on 20 accidental scenarii based on 10 experimental installations.

In addition, the Cathare code is daily used on reactor simulators for the training of operators. During an accidental scenario, the role of the operator is to activate emergency procedures to keep the reactor in nominal working conditions. Reactor simulators are used by IRSN and other French and foreign organisations to train their engineers. Among the available tools to make nuclear safety studies, SOFIA -Simulateur d'Observation du Fonctionnement Incidentel et Accidentel [START_REF]un simulateur pour améliorer la sûreté des réacteurs à eau sous pression[END_REF], is a computer system allowing the real-time tracking of the evolution of many physical parameters in a nuclear reactor. It simulates the material failures and actions of the operator, it also stops the computation to analyse the installation at a given time and step back to modify the scenario. The SOFIA simulator uses the Cathare code in real-time.

Aims of this thesis

The Cathare code is a system code and describes the whole reactor as an assembly of pipes, vessels and pumps whose sizes are around 10 meters (size of the vessel: 13m × 5m × 5m). To simulate the whole system, the Cathare code has restrictions on the mesh size: taking into account the important size of the circuits, the mesh size used can be relatively large, from a few centimeters to a meter. After the discretisation of all the variables, typical cases involve up to 10 2 or 10 3 cells with 3D elements and involve up to a million of numerical time steps. A space domain decomposition method is implemented in the Cathare code and reaches its limits in its ability to use the entire computational resources. The scalability properties of the space domain decomposition method are limited by the small number of cells in the meshes of the Cathare simulations. We seek in this work to investigate a novel strategy of parallelisation to complement the actual parallelism in the Cathare code. For this reason, if we have more processors at our disposal and wish additional speed-ups, the parallelisation of other variables needs to be addressed. Our purpose is to design a strategy of time domain decomposition. We would like to use a non intrusive approach where the Cathare code is used as a black box. In this context we investigate the ability of the parareal method to match our requirements. This work constitutes the first part of this thesis and is presented in chapters 1, 2 and 3.

The remaining chapters of this thesis are motivated by the need of novel numerical methods that could handle better the numerical difficulties that are specific to two-phase flows with a lower computational cost. This second part is dedicated to the understanding of the theoretical properties of finite volume schemes on staggered grids such as the one used in the Cathare code. We seek to develop a rigorous framework for the stability analysis of classical staggered schemes and to propose a class of L 2 -stable staggered schemes. In addition, we seek to derive a class of entropic staggered schemes. The procedure derived in the thesis is very general and could be applied to two-phase flows models.

In the following section, a summary of every chapter will be provided.

Summary of the results by chapters Part I: Chapter 1

The aim of this first chapter is twofold: first, it is intended to provide a bibliographical overview of the two-phase flow model used in the Cathare code, namely the six-equation two-fluid model. Most of what is stated here is not new but it seemed interesting to us to present this compilation of information to position the model and numerical methods in the Cathare code within the existing litterature of two-phase flows simulation.

After introducing the six-equation two-fluid model, the main theoretical results regarding the hyperbolicity of the model will be presented. We will then explain the complex nature of the solution displayed by the model that is specific to two-phase flows models.

We will continue by recalling the existing discretisation techniques of the variables involved in the equations and a special emphasis will be put ont the treatment of the convection term for the most widespread staggered schemes. Even if two-fluid models inherit achievements performed in the modeling, mathematical theory and numerical methods for single-phase flows, they however display many specific difficulties. We will discuss some difficulties existing in the two-phase flow models such as the presence of non conservative products, the configuration of the vanishing phase and the handling of discontinuous source terms.

The last part of the chapter is devoted to the existing acceleration techniques that are actually available in the Cathare code. First, these methods will be presented as well as their limits in the ability to use the entire computational resources. We finally present the strategy followed in the Cathare code to parallelise the time variable. The parallelisation of the time variable is particularly involved given the sequential nature of time. Despite this, several strategies have been proposed in the literature (see [START_REF] Chartier | A parallel shooting technique for solving dissipative ODE's[END_REF], [START_REF] Gander | 50 years of time parallel time integration[END_REF]). We have focused on the parareal in time method because it is the one that seems to provide the best performances with many applications (see [START_REF] Baffico | Parallel-in-time molecular-dynamics simulations[END_REF], [START_REF] Fischer | A parareal in time semi-implicit approximation of the Navier-Stokes equations[END_REF], [START_REF] Samaddar | Parallelization in time of numerical simulations of fully-developed plasma turbulence using the parareal algorithm[END_REF] among many others).

Part I: Chapter 2

The second chapter summarises the special efforts we made to implement the parareal algorithm to the Cathare code in a non intrusive way. We will present the two computational tools we developped during the PhD in order to apply the parareal algorithm to the Cathare code: firstly a numerical clone of Cathare that is restricted to one test case and secondly through a library that uses the Cathare code as a black box in a parallel way. The two test cases we investigate are representative of the numerical challenges for the simulation of two phase flows in the context of safety studies. Numerical challenges include for instance, the vanishing phase issue where one of the two phases liquid or gas disappears in some parts of the domain or the simulation of a breach in the primary circuit that causes a fast depressurisation within the reactor core.

The main contribution of this work has been to adapt the parareal algorithm to the architecture of the software and to its time discretisation in a non intrusive way, without modification of the source files of the Cathare code, in order to reduce the computational time and get closer to a real-time response of the code.

In the numerical examples treated, the use of the parareal algorithm can speed-up the calculations by a factor of about 3 with 25 processors. From an efficiency point of view, these results are not as competitive as the high efficiency that domain decomposition methods provide, but, as it will be explained in chapter 2, there are theoretical reasons that explain the relatively low efficiency of the parareal method. Because of this fact, parareal is a useful technique to obtain additional speed-ups in the context where other more efficient parallelisation techniques have reached saturation.

Part I: Chapter 3

As will be presented in detail in chapters 1 and 2, the time discretisation of the Cathare code relies on a two-step time scheme. A multi-step time scheme can potentially bring higher approximation orders than plain one-step methods but the initialisation of the time propagation in each time window needs to be appropriately chosen. When the fine and/or coarse propagators is a multi-step time scheme, we need to choose a consistent approximation of the solutions involved in the initialisation of the fine solver at each time window. Otherwise, an initialisation error would be propagated over the whole time interval and would prevent the parareal algorithm to converge towards the solution with the desired accuracy.

In an attempt to address this issue, this chapter presents a new variant of the parareal algorithm, adapted to this type of discretisation, and that ensures to recover the target solution with a convergence rate similar to the one of the classical parareal algorithm. A special effort has been made to design an algorithm adapted to this type of discretisation without being intrusive in the coarse or fine solvers. This allows to the Cathare code to be treated as a black box, which ensures the portability of this new algorithm. With regard to the initialisation procedure, the multi-step parareal algorithm is more consistent with the underlying time scheme. We show both theoretically and numerically that the accuracy and convergence of the multi-step parareal algorithm are very competitive when we choose carefully the initialisation of each time window.

Part II: Chapter 4

The second part of this thesis is devoted to the development of a rigorous framework for the analysis of finite volume schemes on staggered grids such as the one used in the Cathare code. In particular, the family of schemes presented here could be applied in the future to the Cathare code to handle the numerical difficulties specific to two-phase flows models. The derivation of the method has required the analysis of some theoretical aspects beforehand and this is what is presented in this second part of the manuscript. Nevertheless, some analytical numerical examples will be presented on the solution of the isentropic Euler equations with the purpose of illustrating the technique and its performances.

Staggered schemes are known to be more precise for almost incompressible flows in practice and are very popular in the thermal hydraulics community ( [START_REF]A numerical fluid dynamics calculation method for all flow speeds[END_REF][START_REF] Prosperetti | Computational methods for multiphase flow[END_REF][START_REF]Discrétisation des équations du module 3d de cathare 2 & comparaison avec le module 1d -etat des lieux pour la version v2.5 2 et perspectives[END_REF]). However, in the context of compressible flows, their stability analysis has historically been performed with a heuristic approach and the tuning of numerical parameters ( [START_REF] Hirt | Heuristic stability theory for finite difference equations[END_REF]). Yet the conservative staggered schemes presented in [START_REF]On some implicit and semi-implicit staggered schemes for the shallow water and euler equations[END_REF][START_REF] Herbin | Staggered schemes for all speed flows[END_REF] are proven to be entropic and to satisfy a kinetic energy preservation [START_REF] Herbin | A kinetic energy preserving convection operator for the mac discretization of compressible navier-stokes equations[END_REF]. Unfortunately, the boundedness of the entropy does not necessarily imply the boundedness of the solution and this is particularly the case for the full Euler system. Hence, we investigate in this chapter the linear L 2 -stability of staggered schemes. This question has been addressed by analysing their numerical diffusion operator that gives a new insight into these schemes. We first derive the numerical diffusion operator for classical staggered schemes and show that the L 2 stability is obtained only in the case of constant sign velocities. We then propose a class of linearly L 2 -stable staggered schemes. This new class is based on a carefully chosen numerical diffusion operator and is more adapted to two-phase flows where phasic velocities frequently change signs.

A scheme belonging to this new class of staggered schemes has been successfully applied to the simulation of a Riemann problem and we expect that the present developments will allow its use in more realistic and complex cases in the future.

Part II: Chapter 5

In the last chapter, we address the question of the entropic character of staggered schemes. We derive a discrete entropy balance for a class of staggered schemes. On this basis, we give explicit constraints on the coefficients of the numerical diffusion matrix to ensure the dissipation of the discrete entropy. The procedure is very general and could be applied to two-phase flows models. We first investigate this strategy on the isentropic Euler system. We then implement a scheme belonging to this new class of entropic staggered schemes for the simulation of a Riemann problem that displays transonic rarefaction waves. These numerical results illustrate the ability of our method to capture the correct entropic solution in a stable way.

Part I

Numerical models for two-phase flows for safety studies This first chapter is intended to be a bibliographic summary about the six-equation two-fluid model: an overview of some theoretical results, discretisation, numerical and HPC methods to solve the model are presented. We will also show the numerical difficulties that are specific to two-phase flows and represent an obstruction to obtain a satisfactory simulation at a reasonnable computational cost. In this context, the main contribution of this work has been to explore parallel acceleration techniques to reduce this computational time (see chapters 2 and 3) and also to explore numerical methods that could handle better these numerical difficulties with a lower computational cost (see part II of this manuscript).

The six-equation two-fluid model 1.Balance equations

Considering a two-phase flow of vapour (subscript v ) and liquid (subscript l ) phases, we are interested in determining unknown physical quantities that are the vapour and liquid volume fractions

(α v , α l ∈ [0, 1]), average pressures (p v ∈ R + and p l ∈ R + ), densities (ρ v ∈ R + and ρ l ∈ R + ), temperatures (T v ∈ R + and T l ∈ R + ), velocities ( u v ∈ R 3 and u l ∈ R 3 ), enthalpies (h v ∈ R + and h l ∈ R + ) and total energies (E v ∈ R + and E l ∈ R + ).
The two-fluid model equations are obtained by averaging the balance equations for each separated phase, using space, time or ensemble averaged quantities (see [START_REF] Drew | Theory of multicomponent fluids[END_REF][START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flow[END_REF][START_REF] Ishii | Thermo-Fluid Dynamics of Two-Phase Flow[END_REF]). The mass, momentum and energy balance equations for each phase then read as follows:

                                               ∂α v ρ v ∂t + ∇ • (α v ρ v u v ) = Γ v , ∂α l ρ l ∂t + ∇ • (α l ρ l u l ) = Γ l , ∂α v ρ v u v ∂t + ∇ • (α v ρ v u v ⊗ u v ) + α v ∇p v = ∇ • (α v ( τv + t τv )) + α v ρ v g + F int v + F w v + Γ v u int v , ∂α l ρ l u l ∂t + ∇ • (α l ρ l u l ⊗ u l ) + α l ∇p l = ∇ • (α l ( τl + t τl )) + α l ρ l g + F int l + F w l + Γ l u int l , ∂α v ρ v E v ∂t + p int v ∂α v ∂t + ∇ • (α v ρ v u v (E v + p v ρ v )) = α v ρ v u v • g + F int v • u int v + σ Q v +∇ • (α v (q v + t q v )) + Γ v (h int v + ( u int v ) 2 2 ) + Q w v , ∂α l ρ l E l ∂t + p int l ∂α l ∂t + ∇ • (α l ρ l u l (E l + p l ρ l )) = α l ρ l u l • g + F int l • u int l + σ Q l +∇ • (α l (q l + t q l )) + Γ l (h int l + ( u int l ) 2 2 ) + Q w l .
(1.1) One also has to consider the closure relation:

α v + α l = 1,
which expresses the fact that the whole available volume is filled up with the fluids. The interface interaction transfers consist in the following terms: Γ k the interfacial mass transfer terms, F int k the interfacial momentum transfer terms, which are either of algebraic or differential nature. They result of the viscous stresses at the interface and pressure forces normal to the interface. They can be further divided into a non-viscous component F nv k containing for instance the interfacial pressure term, the virtual mass force and a viscous component such as the interfacial friction term. The interfacial pressure term takes the following form:

F p k = (p int k -p k ) ∇α k , u int the interfacial velocity, σ Q k the interfacial heat transfer terms.
The wall interaction modeling, denoted by the subscript w, includes wall heat transfer and wall frictional forces:

Q w k the wall heat transfer rates, F w k the wall frictional forces. τk are the Reynolds stresses and q k are the turbulent heat fluxes.

Equations of state

The number of unknown variables is in general greater than the number of equations in the system of PDEs. Therefore, it is necessary to consider supplementary constitutive equations such as equations of state. In general, such equations depend on the specific two-phase flow model and on the flow regimes. The equations of state considered are:

ρ k = ρ k (p, e k ).
The Cathare code uses tabulated equations of state based on the industrial formulation IAPWS (The International Association for the Properties of Water and Steam, [START_REF] Wagner | The IAPWS industrial formulation 1997 for the thermodynamic properties of water and steam[END_REF]). These equations of state compute all the thermodynamical variables necessary for the simulation thanks to experimental measurements of pressures and temperatures and polynomial interpolations.

Closure laws

The jump conditions at the interface are:

Mass transfer: k=g,l Γ k = 0.
Interface momentum transfer:

k=g,l F int k = 0.
Energy transfer:

k=g,l (Γ k h int k + σ Q k ) = 0.
In our work in this thesis, we have worked under the following assumptions:

We assume here that the Reynolds stresses and turbulent heat fluxes are negligible, as we are considering convection driven flow where viscosity tensor play a minor role.

There is only one bulk average pressure in the system. We assume that the pressures relaxation time is negligible, and pressure equilibrium is considered to be reached in the flow:

p v (ρ v , E v ) = p l (ρ l , E l ) = p.
We neglect surface tension phenomena in the two-fluid model:

p int v = p int l = p int .
We neglect as many authors do the pressure default (pp int ) in the energy equation.

Providing the general expressions of the momentum transfer terms is a difficult task with many ongoing works. We summarise in the sequel some expressions that are commonly employed:

Interfacial pressure term If we neglect virtual mass force, the Cathare model for the interface pressure force is given by [START_REF] Bestion | The physical closure laws in the cathare code[END_REF] ∆p = p -

p int = δ α v α l ρ v ρ l α v ρ l + α l ρ v || u v -u l || 2 .
The coefficient δ can be gauged such that the two-fluid model becomes hyperbolic ( [START_REF] Ndjinga | Influence of interfacial pressure on the hyperbolicity of the two-fluid model[END_REF]). Let us note that the Cathare model for the interfacial pressure term is slightly different in the case of stratified flows. We will give more details about this point in the section dedicated to the hyperbolicity of the six equation two-fluid model.

Virtual mass force

The dynamic drag or the transient forces that result from relative acceleration of the phases, usually also called the added mass effects is modeled in Cathare with the following formula, issued from [20]

F vm v = F vm l = -β vm α v α l ρ m ∂ u v ∂t + u v • ∇ u v - ∂ u l ∂t + u l • ∇ u l ,
where β vm is the virtual mass coefficient depending on the flow regime.

Interfacial friction term or Drag force

The drag model depends on the flow regime (dispersed, stratified, ...) and has the following general form:

F D = 1 2 C D a i ρk(α)||u v -u l || 2 ,
where C D is the drag coefficient and k(α) is a function that strongly couples the phases when one of them tends to disappear. More details will be given in section 1.2.5 dedicated to one important numerical difficulty of the two-fluid model: the vanishing phase.

Interfacial velocity

Assuming that the no slip condition is satisfied, the interfacial velocities for momentum and energy transfer are equal for both phases u int v = u int l = u int . The formulation used in Cathare is a volume fraction averaged formulation:

u int = α v u l + α l u v .

Boundary conditions and initial condition

We will assume that the domain Ω is a bounded open set of R 3 . We denote the boundary ∂Ω with n(X) the outward unit normal to ∂Ω at point X. We define the following partitions of ∂Ω:

∂Ω = ∂Ω in ∪ ∂Ω out ∪ ∂Ω wall .
We list here some of the usual boundary conditions that are associated to problem (1.1).

Inlet and outlet boundary conditions

The flow velocity profile is specified at inlet boundaries to model the incoming of liquid or vapour in the domain:

u k (t, X) • n(X) = u in k (t), ∀t ∈ [0, T ] and X ∈ ∂Ω in .
The pressure is specified at outlet boundaries to model the outgoing of liquid or vapour in the domain. Outlet boundary conditions for the pressure give informations external to the domain: depending on the test case, the outlet pressure can be equal to the atmospheric pressure (1×10 5 P a) or to the pressure within a nuclear reactor under nominal working conditions (155 × 10 5 P a). p(t, X) = p out (t), ∀t ∈ [0, T ] and X ∈ ∂Ω out .

In practice, one generally uses a ghost cell formulation at the boundary. The inlet boundary condition consists in imposing a Dirichlet condition for the velocity and a Neumann condition for the other variables (pressure, volume fractions, enthalpies). For Neumann boundary conditions, we prescribe the gradient normal to the boundary of a variable at the boundary. Usually, we take the state of the ghost cell equal to the internal state.

Wall boundary conditions

This is commonly known as no-slip boundary condition. It specifies the conditions for velocity components at the wall. The normal component is set to zero since the wall is static:

u k (t, X) • n(X) = 0, ∀t ∈ [0, T ] and X ∈ ∂Ω wall .
The tangential component is set to the velocity of the wall:

u k (t, X) • t(X) = u wall k (t), ∀t ∈ [0, T ] and X ∈ ∂Ω wall .
Heat transfer through the wall can be specified or set to zero in the case of adiabatic walls.

Periodic boundary conditions

Periodic boundary conditions are often used to simulate a large system by modeling a small part that is far from its edge. These conditions consist in enforcing a relation of the form:

u k (t, X) = u k (t, X ),
For X and X in ∂Ω. For instance, if Ω = [0, L] in a 1D case, the condition would read:

u k (t, 0) = u k (t, L), ∀t ∈ [0, T ].
The initial condition depends on the situation under consideration. In the analysis of reactor cores, what one wishes in the end is to understand the connection between a stationary state and some transient state. In practise, there are two ways of initialising a simulation:

First one can start with constant values of the unknowns. However this can set the system in a state that is very far from the stationnary state. Furthermore the transient dynamics can be very different from the one observed when one studies nuclear accident.

The second option is the one used in thermalhydraulics (in particular in the Cathare code). It consists in computing first a stationnary condition in a normal regime using the values of the boundary condition then changing the source term and the boundary condition to reflect the occurence of a nuclear accident.

Hyperbolicity of the model

When considering equal pressure law and neglecting the interfacial pressure term and the virtual mass term, the original two-fluid model (1.1) is not unconditionnally hyperbolic in the low Mach flow regime in nuclear safety studies:

u r = u v -u l
c sm , with the mixture sound speed:

c sm = α v ρ v + α l ρ l α v ρ l c v + α l ρ v c l ,
where c v , c l are the phasic sound speeds. Taking into account interfacial pressure default or the virtual mass the system becomes hyperbolic in the flow regime of interest. The reader is referred to [START_REF] Stuhmiller | The Influence of Interfacial Pressure Forces on the Character of Two-Phase Flow Model Equations[END_REF][START_REF] Bestion | The physical closure laws in the cathare code[END_REF][START_REF] Toumi | An Approximate Linearized Riemann Solver for a Two-Fluid Model[END_REF][START_REF] Amine | Modélisation et analyse numérique des écoulements diphasiques en déséquilibre[END_REF][START_REF] Toumi | Approximate Riemann solvers and flux vector splitting schemes for two-phase flow[END_REF][START_REF] Ndjinga | Influence of interfacial pressure on the hyperbolicity of the two-fluid model[END_REF] for details. In the sequel, we only study the effect of the interfacial pressure correction on the hyperbolicity of the model.

In this section, we first give the Taylor expansion of the spectrum for the compressible six-equation model in 1.2.1. The complex expression of the eigenvalues makes difficult a rigorous analysis of the compressible model, hence we turn our attention to the study of the low Mach regime in 1.2.2. In section 1.2.2, we show the effect of the interfacial pressure term in the incompressible case and give closure laws that ensure the hyperbolicity of the incompressible model. In section 1.2.3, we list some closure laws for the interfacial pressure term that ensure the hyperbolicity of the compressible model. In section 1.2.4, we summarise the analysis of the Riemann problem for two incompressible phases. This analysis also takes into account the case of a vanishing phase where one of the volume fraction α k goes to zero. In the last section 1.2.5, we show the singularity that appears in the case of a vanishing phase and the difficulty to ensure the positivity of the volume fractions α k .

Eigenstructure of the compressible six-equation model

In [START_REF]Quelques aspects de modélisation et d'analyse des systèmes issus des écou-lements diphasiques[END_REF], the author gives the spectral properties of the three-dimensional two-fluid model: six real eigenvalues are trivially computed and the remaining eigenvalues are the roots of a 4-th degree polynomial that is exactly the characteristic polynomial of the one dimensional isentropic two-fluid system:

                   ∂α v ρ v ∂t + ∂α v ρ v u v ∂x = 0, ∂α l ρ l ∂t + ∂α l ρ l u l ∂x = 0, ∂α v ρ v u v ∂t + ∂α v ρ v u 2 v ∂x + α v ∂p ∂x + ∆p ∂α v ∂x = 0, ∂α l ρ l u l ∂t + ∂α l ρ l u 2 l ∂x + α l ∂p ∂x + ∆p ∂α l ∂x = 0. (1.2)
Denoting the unknown variable U = (α v ρ v , α v ρ v u v , α l ρ l , α l ρ l u l ), the system (1.2) can be rewritten as the quasi-linear form:

∂U ∂t + A(U ) ∂U ∂x = 0
where A(U ) is the Jacobian matrix of system (1.2).

For practical purposes, one usually does not find the exact solution of the fourth degree characteristic polynomial. Instead, following the works in [START_REF] Toumi | An Approximate Linearized Riemann Solver for a Two-Fluid Model[END_REF] and [START_REF] Evje | Hybrid flux-splitting schemes for a common two-fluid model[END_REF], the authors suggest using a perturbation method to compute approximate eigenvalues. Denoting the perturbation parameter:

= u v -u l c sm = u r c sm
where c sm is an approximate mixture sound speed defined by:

c sm = α v ρ v + α l ρ l α v ρ l c v + α l ρ v c l
The analysis of the eigenvalues for the six-equation two-fluid model is approximately made around the mechanical equilibrium, i.e. the eigenvalues are computed as a perturbation of the relative velocity in comparison to the mixture sound speed, ( = u vu l c sm ). The first order approximation of the compressibel two-fluid system eigenvalues gives:

       λ 1,4 = α v ρ l u v + α l ρ v u l α v ρ l + α l ρ v ∓ c sm + O( 2 )
, pressure or acoustic waves

λ 2,3 = α v ρ l u l + α l ρ v u v α v ρ l + α l ρ v ∓ 1 α v ρ l + α l ρ v ∆p - α v α l ρ v ρ l α v ρ l + α l ρ v u 2 r + O( 2 ), void waves (1.
3) We see from (1.3) that for small relative velocities, the two eigenvalues λ 1,4 are always real and have the order of magnitude of the mixture speed of sound c sm . In contrast, if ∆p < α v α l ρ v ρ l α v ρ l + α l ρ v u 2 r , the two other eigenvalues λ 2,3 are complex. We will discuss the influence of the interfacial pressure term on the hyperbolicity of the system in the next section. Assuming that the four eigenvalues are real, we can see that the eigenvalues may easily change sign. The two-fluid model displays a complicated eigenstructure and its spectrum is more complex than the one of for single-phase flows. For single-phase flows, both mathematical theory and numerical methods of the Euler system have been studying well by numerous authors in the litterature. This system is a system of hyperbolic conservation laws with two Genuinely Non Linear fields and one Linearly Degenerate field (in one dimension). Even if the two-fluid models inherit achievements obtained in the single phase flow including of modeling, mathematical theory and numerical methods, however, the two-phase flow models possess many of specific difficulties due to existence of two phases in the same domain of interest and their interactions as well. We will see in section 1.2.4 that for the incompressible limit of system (1.2) the study of the spectrum proves that the eigenvalues are not a priori ordered and that the characteristic fields are neither GNL nor LD.

Hyperbolicity of the incompressible model

From the system (1.3), we see that there is a critical ∆p which ensures the positivity of the value under the square root. In this section, we show the influence of the interfacial pressure term on the hyperbolicity of the two-fluid model in the incompressible case. To understand the role of the interfacial pressure term to get the hyperbolicity of the system, we show the eigenstructure of a reduced system where both phases are assumed incompressible in (1.2), derived in [START_REF]Quelques aspects de modélisation et d'analyse des systèmes issus des écou-lements diphasiques[END_REF].

∂ ∂t ρ ρu + A red ∂ ∂x ρ ρu = 0 0 , where: ρ ρu = α v ρ l + α l ρ v ρ v u v -ρ l u l , (1.4)
with

A red =     α l ρ v u v + α v ρ l u l α v ρ l + α l ρ v - α v α l (ρ v -ρ l ) α v ρ l + α l ρ v ρ v ρ l (u v -u l ) 2 (ρ v -ρ l )(α v ρ l + α l ρ v ) - ∆p α v α l (ρ v -ρ l ) α l ρ v u v + α v ρ l u l α v ρ l + α l ρ v     ,
and its eigenvalues are:

λ inc = α l ρ v u v + α v ρ l u l α v ρ l + α l ρ v ± 1 α v ρ l + α l ρ v ∆p - α v α l ρ v ρ l (u v -u l ) 2 α v ρ l + α l ρ v .
The incompressible limit of the two-fluid model is hyperbolic provided:

∆p ≥ α v α l ρ v ρ l (u v -u l ) 2 α v ρ l + α l ρ v (1.5)
This type of coefficient is used in the industrial code Cathare (see [START_REF] Bestion | The physical closure laws in the cathare code[END_REF]) to make the compressible two-fluid model hyperbolic.

If we remove the interfacial pressure term by setting ∆p = 0, the matrix A red of system (1.4) becomes:

A red =     α l ρ v u v + α v ρ l u l α v ρ l + α l ρ v - α v α l (ρ v -ρ l ) α v ρ l + α l ρ v ρ v ρ l (u v -u l ) 2 (ρ v -ρ l )(α v ρ l + α l ρ v ) α l ρ v u v + α v ρ l u l α v ρ l + α l ρ v     (1.6)
In (1.6), A red is a 2 × 2 matrix with two diagonal coefficients that are equal. Such a matrix has real eigenvalues provided the extradiagonal terms have the same sign. The efficiency of the parameter ∆p in making the system hyperbolic comes from the fact that it removes the component:

ρ v ρ l (u v -u l ) 2 (ρ v -ρ l )(α v ρ l + α l ρ v )
that gives rise to complex eigenvalues, as we can see in (1.4). Thus to make the system hyperbolic, another possibility is to cancel this component by assuming an algebraic pressure disequilibrium. Without the assumption of pressure equality, the evolution equation of ρu would be:

∂ ρu ∂t + ∂ ∂x ( ρu) 2 2(ρ v -ρ l ) + ∂ ∂x p v -p l - ρ v ρ l (u v -u l ) 2 (ρ v -ρ l )(α v ρ l + α l ρ v ) = 0
An alternative to the use of an interfacial pressure correction ∆p is the assumption that there is an algebraic pressure disequilibrium having the form:

p v -p l = ρ v ρ l (u v -u l ) 2 (ρ v -ρ l )(α v ρ l + α l ρ v ) (1.7)
Hence, one can use a pressure disequilibrium of the form (1.7) or a parameter ∆p of the form (1.5) to ensure the hyperbolicity of the isentropic two-fluid model for two incompressible phases. In [START_REF] Ndjinga | A 2 × 2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF], the authors prove the existence and uniqueness of an admissible solution to the Riemann problem for the isentropic two-fluid model for two incompressible phases with a pressure disequilibrium of the form (1.7). In section 1.2.4, we detail the properties of this model and the results obtained in [START_REF] Ndjinga | A 2 × 2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF].

Hyperbolicity of the compressible model

In this section, we give some closures laws for the isentropic two-fluid model (1.2) in the compressible case that ensure the hyperbolicity of the system for a range of relative velocities u r .

For small relative velocities (u vu l ), closure laws were proposed in [START_REF] Stuhmiller | The Influence of Interfacial Pressure Forces on the Character of Two-Phase Flow Model Equations[END_REF] and [START_REF] Bestion | The physical closure laws in the cathare code[END_REF] that ensure the hyperbolicity of the system (1.2) with the following form:

∆p = δ α v α l ρ v ρ l (u v -u l ) 2 α v ρ l + α l ρ v , with: δ > 1
In the more general case of large relative velocities, the hyperbolicity of the isentropic system (1.2) has been studied in [START_REF] Amine | Modélisation et analyse numérique des écoulements diphasiques en déséquilibre[END_REF], [START_REF] Ndjinga | Influence of interfacial pressure on the hyperbolicity of the two-fluid model[END_REF], [START_REF]Quelques aspects de modélisation et d'analyse des systèmes issus des écou-lements diphasiques[END_REF], [START_REF] Vernier | Validation de Modèles Multichamps dans le Logiciel OVAP[END_REF] with:

∆p = α v α l ρ v ρ l α v ρ l + α l ρ v (u v -u l ) 2 + 1 c 2 v ρ v - α v α l ρ v ρ l α v ρ l + α l ρ v (u v -u l ) 4 .
This closure law guarantees the hyperbolicity for the relative velocities up to the sound speed of vapour phase c v . On the other hand, the authors in [START_REF] Amine | Modélisation et analyse numérique des écoulements diphasiques en déséquilibre[END_REF][START_REF]Quelques aspects de modélisation et d'analyse des systèmes issus des écou-lements diphasiques[END_REF] introduces the interfacial of the form:

∆p = ρ v (u v -u l ) 2 .
This closure law guarantees the hyperbolicity in the region where |u vu l | ≤ c v .

In the next section, we will present some recent existence results obtained on a system for two incompressible phases.

The Riemann problem for two incompressible phases

In [START_REF] Ndjinga | A 2 × 2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF], the authors study a 2 × 2 system of conservation laws that models the dynamic of two incompressible phases. They establish their model from the following system:

       ∂ t α v ρ v + ∂ x (α v ρ v ) = 0 ∂ t α l ρ l + ∂ x (α l ρ l ) = 0 ∂ t (α v ρ v u v ) + ∂ x (α v ρ v u 2 v ) + α v ∂ x p v = α v ρ v g ∂ t (α l ρ l u l ) + ∂ x (α l ρ l u 2 l ) + α l ∂ x p l = α l ρ l g (1.8)
where the non zero pressure difference p vp l takes the form (1.7).

Assuming in (1.8) that both phases are incompressible, the resulting system focuses on the study of the void waves that determine the composition of the mixture:

       ∂ t α + ∂ x α(1 -α)ω α(ρ l -ρ v ) + ρ v = 0 ∂ t ω + ∂ x ω 2 2(ρ v -ρ l ) = (ρ v -ρ l )g (1.9)
We introduce the unknown variable vector U = α ω , and the space of admissible states H = {(α, ω), α ∈ [0, 1], ω ∈ R}. We also denote the sets H + and H -such that:

H ± = {(α, ω), ω ∈ R ± \ {0} , α ∈ (0, 1)}
Many specificities arise from this incompressible model: the two eigenvalues λ 1 and λ 2 are not a priori ordered:

λ 1 = ω ρ v -ρ l 1 - ρ v ρ l (α(ρ l -ρ v ) + ρ v ) 2 , λ 2 = ω ρ v -ρ l
the characteristic fields associated to λ 1 and λ 2 are genuinely nonlinear in each domain H + and H -but are neither genuinely nonlinear nor linearly degenerate in general:

∇λ 1 • r 1 = -2ρ v ρ l ω (α(ρ l -ρ v ) + ρ v ) 3 , ∇λ 2 • r 2 = ρ v ρ l ω (ρ v -ρ l )
where r 1 , r 2 are the eigenvectors of ∇F associated to λ 1 , λ 2 .

the system is strictly hyperbolic in H ± and in general weakly hyperbolic on the domain H Considering the Riemann problem for the conservative system (1.9) in the case g = 0 with a piecewise constant initial data: 3a in [START_REF] Ndjinga | A 2 × 2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF]).

U 0 (x) = U L (α L , ω L ) if x ≤ 0 U R (α R , ω R ) if x >

Vanishing phase

In the two-fluid model, the total boiling or condensation of one phase will arise a singularity. The absent phase is called vanishing phase or ghost phase. It poses a difficulty in the two-fluid model owing to its independent velocities. The singularity arises when one computes the absent phase velocity using the conservative variables

u k = α k ρ k u k α k ρ k as α k ρ k → 0.
Studying the one dimensional two-fluid model, we are interested in the mathematical properties of the vanishing phase which is assumed to be vapour: α = α v = 0. Let us take into account the isentropic model (1.2), when α = 0, the Jacobian matrix becomes:

A α=0 =       0 0 1 0 0 0 0 1 -u 2 v 0 2u v 0 ρ l c 2 l ρ v c 2 l -u 2 l 0 2u l      
which has four real eigenvalues: u v , u v , u l + c l , u lc l but the matrix is not diagonalizable because there are only three eigenvectors:

v u l ±c l = t (0, 0, 1, u l ± c l ) , v uv = t ρ v ρ l ( (u v -u l ) 2 c 2 l -1), u v ρ v ρ l ( (u v -u l ) 2 c 2 l -1), 1, u v
The hyperbolicity of system (1.2) is then broken for α v = 0. One of the objectives of [START_REF] Ndjinga | A 2 × 2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF] was to prove that the positiveness of the volume fractions α k is inherent to the model. At the discrete level, it is not always the case that a numerical method capture positive volume fraction. In the thermalhydraulic platform Cathare, an interfacial friction term is used to ensure the positiveness of the volume fractions. We explain here the role of the interfacial friction term in ensuring that the volume fractions α k (t) ∈ [0, 1], ∀t ∈ [0, T ]. We consider in the sequel that the unknowns are smooth enough:

     ∂ t u v + u v ∂ x u v + ∂ x p ρ v + ∆p α v ρ v ∂ x α v = C D ρ v ||u v -u l ||(u v -u l ) ∂ t u l + u l ∂ x u l + ∂ x p ρ l + ∆p α l ρ l ∂ x α l = C D ρ l ||u l -u v ||(u l -u v ) (1.10)
We now write the momentum equation for the relative velocity u r = u vu l :

∂ t (u v -u l ) + u v ∂ x u v -u l ∂ x u l + 1 ρ v - 1 ρ l ∂ x p = -C D 1 ρ v + 1 ρ l ||u v -u l ||(u v -u l )
Since the variables u v , u l , p and their derivatives are assumed bounded, when the drag coefficient C D goes to infinity, the relative velocity u r verifies the ODE:

∂ t u r = -C D ||u v -u l ||(u v -u l )
Hence, the system tends to an equilibrium of the velocities with u r → 0. The numerical difficulty of the vanishing phase is thus handled by the Cathare code by imposing the equality of the two phasic velocities through an interfacial friction term that dominates the momentum equations when one of the volume fraction goes to 0. Assuming that the vapour phase is the vanishing phase and neglecting the interfacial friction coefficient C D , the corresponding momentum equation in (1.10) shows that the velocity of the vanishing phase follows a Burgers equation and does not have to be equal to the velocity of the remaining phase.

Discretisation of the two-fluid model

In general, there exists two families of numerical methods for the simulation of two-phase flows. Firstly, colocated schemes are generally used on unstructured meshes where the unknowns are located in the same place (cell-centered). In the litterature, many authors developed Riemann solvers (either Godunov-type methods or Roe-type schemes or Osher schemes or AUSM schemes) for the simulation of two-phase flows dealing with the numerical challenges encountered: vanishing phase [START_REF] Cordier | Phase Appearance or Disappearance in Two-Phase Flows[END_REF][START_REF] Coquel | A numerical method using upwind schemes for the resolution of two-phase flows[END_REF][START_REF] Munkejord | A musta scheme for a nonconservative two-fluid model[END_REF][START_REF] Shekari | Application of Osher and PRICE-C schemes to solve compressible isothermal two-fluid models of two-phase flow[END_REF][START_REF] Ambroso | A Godunov-type method for the sevenequation model of compressible two-phase flow[END_REF][START_REF] Evje | Hybrid flux-splitting schemes for a common two-fluid model[END_REF][START_REF] Ndjinga | A 2 × 2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF], non conservative products [START_REF] Munkejord | A musta scheme for a nonconservative two-fluid model[END_REF][START_REF] Morin | A roe scheme for a compressible sixequation two-fluid model[END_REF][START_REF] Evje | Hybrid flux-splitting schemes for a common two-fluid model[END_REF][START_REF] Toumi | An Approximate Linearized Riemann Solver for a Two-Fluid Model[END_REF] and stiff source terms [START_REF] Ndjinga | A 2 × 2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF][START_REF] Nguyen | On the mathematical analysis and the numerical simulation of boiling flow models in nuclear power plants thermal hydraulics[END_REF]. We can also mention the VFFC scheme ( [START_REF] Ghidaglia | Une méthode volumes finis à flux caractéristiques pour la résolution numérique des systèmes hyperboliques de lois de conservation[END_REF][START_REF] Ghidaglia | On the numerical solution to two fluid models via a cell centered finite volume method[END_REF]) that has a generic formulation by contrast with the Roe scheme ( [START_REF] Roe | Approximate riemann solvers, parameter vectors and difference schemes[END_REF]) that is applied under some algebraic conditions on the system. This category of schemes is robust but present a lack of accuracy for low Mach number/almost incompressible flows. Corrections are proposed in [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF][START_REF] Dellacherie | Analysis of godunov type schemes applied to the compressible euler system at low mach number[END_REF] to overcome this issue but generate an instability with checker-board type oscillations. We can also mention the pressure-based methods with colocated variables such as the one used in the platform Neptune-CFD (developed by EDF and CEA, [START_REF] Archambeau | Comparative study of pressure correction and Godunov-type schemes on unsteady compressible cases[END_REF]). This type of schemes does not suffer from a lack of precision for low Mach number flows but show checker-board type oscillations. A way to overcome this issue is to use Rhie and Chow type corrections to avoid spurious oscillations, [START_REF] Eymard | A class of collocated finite volume schemes for incompressible flow problems[END_REF].

On the other hand, staggered schemes are used on structured meshes with unknows located either on edges or cell centers. This category of schemes has a good behaviour for almost incompressible flows. The space discretisation in the Cathare code is based on a staggered scheme. The main drawback of this family of schemes is the handling of complex geometries since its use is limited to structured meshes. In this section, we first introduce the family of staggered schemes as initially proposed in [START_REF]A numerical fluid dynamics calculation method for all flow speeds[END_REF][START_REF] Harlow | Numerical calculation of almost incompressible flow[END_REF] (see section 1.3.1.1). We also discuss the properties of these numerical methods depending on the discretisation of the convection term for the isentropic Euler system (see section 1.3.1.2). In section 1.3.1.3, we present a family of staggered schemes derived with a new approach that discretises the conservative form of the Euler system. In section 1.3.2, we detail the discretisation of the one dimensional Cathare model to compare it with the existing litterature on staggered schemes. In the last section, we explain some numerical challenges of two-phase flows simulation by addressing:

strategies of discretisation for non conservative products occuring in two-fluid models (see section 1.

An introduction to ICE schemes

We refer here to the numerical methods designed on staggered grids to simulate compressible flows at low Mach number. The velocity unknowns are located at cell interfaces whilst the density and pressure unknowns are located at cell centers. The expression of the products ρu, ρ∂ t u and ρu 2 then raises an issue since velocity and density are located in different places that is adressed in a different way by different authors.

Historically the Marker and Cell numerical scheme (MAC) [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow with free surface[END_REF] was designed for incompressible flows. Then the Implicit Continuous-fluid Eulerian (ICE) [START_REF]A numerical fluid dynamics calculation method for all flow speeds[END_REF][START_REF] Harlow | Numerical calculation of almost incompressible flow[END_REF] method was designed for compressible flows at low Mach numbers as well as high Mach numbers. The historical ICE method [START_REF]A numerical fluid dynamics calculation method for all flow speeds[END_REF][START_REF] Harlow | Numerical calculation of almost incompressible flow[END_REF] discretises the conservative Navier-Stokes equations in a way that reduces to the MAC method for incompressible flows. In the seminal papers [START_REF]A numerical fluid dynamics calculation method for all flow speeds[END_REF][START_REF] Harlow | Numerical calculation of almost incompressible flow[END_REF], the velocity unknowns are first eliminated then the pressure unknowns are determined in an iterative process and the velocity unknowns follow. For the velocity elimination to be valid the momentum fluxes are treated in a semi-explicit way (see section 1.3.1.1 for more details).

The ICE method encountered a considerable success with numerous variants (explicit/implicit, with/without prediction correction steps) and became popular in the thermal hydraulics community [START_REF] Ninokata | Sabena: Subassembly boiling evolution numerical analysis[END_REF][START_REF] Prosperetti | Computational methods for multiphase flow[END_REF]. There the Navier-Stokes equations are discretised in non conservative form, which makes the velocity elimination easier. The mass flux is upwinded for better stability but the approach still consists in eliminating the velocity unknowns in order to first retrieve the pressure (see 1.3.1.2 for more details). However for the elimination to be rigourously valid the momentum flux should be entirely explicit, which yields a restriction on the time step.

Herbin et al ( [START_REF] Herbin | Staggered schemes for all speed flows[END_REF], [START_REF]On some implicit and semi-implicit staggered schemes for the shallow water and euler equations[END_REF]) proposed an approach that does not rely on velocity elimination and enables full implicitation. They were able to derive rigourous proofs of stability (see section 1.3.1.3 for more details).

We will consider in the sequel the 1D isentropic Euler equations in conservative form .11) and in non conservative form

     ∂ t ρ + ∂ x q = 0 ∂ t q + ∂ x q 2 ρ + ∂ x p = 0. ( 1
   ∂ t ρ + ∂ x (ρu) = 0 ρ∂ t u + 1 2 ρ∂ x u 2 + ∂ x p = 0 , (1.12)
We will give the expression of each scheme on the Euler equations which gives a first insight into the schemes. The discretisation of the momentum convective term u ∂u ∂x should verify the conservation of h = p + 1 2 ρu 2 . Some non conservative schemes discretise the form 1 2 ∂ x u 2 and satisfy the Bernoulli principle. Other non conservative schemes discretise the form u∂ x u and do not naturally recover the Bernoulli principle at the discrete level. In the sequel, we give the expression of several schemes proposed in the litterature in the single-phase case and the extension to two-phase flows can be found in references therein.

The original scheme of Harlow and Amsden

We summarise here the numerical method for all flow speeds described in the seminal article [START_REF]A numerical fluid dynamics calculation method for all flow speeds[END_REF]. In [START_REF]A numerical fluid dynamics calculation method for all flow speeds[END_REF], the stability analysis used the heuristic approach of [START_REF] Hirt | Heuristic stability theory for finite difference equations[END_REF] by estimating the numerical diffusion of the scheme and tuning the numerical viscosity to make it positive. A mass diffusion coefficient τ and artificial viscosity coefficients λ and µ are considered for example in [START_REF] Harlow | Numerical calculation of almost incompressible flow[END_REF]. For simplicity of the exposure, we neglect the various artificial viscosity terms, and present the scheme only in 1D. The discrete 1D isentropic Euler equations take the form

ρ n+1 i -ρ n i t + (ρu) n+1 i+ 1 2 -(ρu) n+1 i-1 2 x = 0 (1.13) (ρu) n+1 i+ 1 2 -(ρu) n i+ 1 2 t + ρ n+1 i+1 (u 2 ) n i+1 -ρ n+1 i (u 2 ) n i x + p n+1 i+1 -p n+1 i x = 0 (1.14)
The expression of the cell centered velocity u 2 i required in the momentum equation raises an issue of interpolation between face and cell that is adressed in different ways by different authors. There are at least four historical types of interpolation formula for u i (see [START_REF]A numerical fluid dynamics calculation method for all flow speeds[END_REF] page 207) :

Centered: u 2 i = u i-1 2 + u i+ 1 2 2 2 ZIP: u 2 i = u i-1 2 u i+ 1 2 Partial Donor: u 2 i =      u i-1 2 u i-1 2 + u i+ 1 2 2 if u i-1 2 + u i+ 1 2 > 0 u i+ 1 2 u i-1 2 + u i+ 1 2 2 if u i-1 2 + u i+ 1 2 < 0 Complete Donor: u 2 i = u 2 i-1 2 if u i-1 2 + u i+ 1 2 > 0 u 2 i+ 1 2 if u i-1 2 + u i+ 1 2 < 0
In the next section, we list some choices of discretisation for the convection term and discuss the properties of the resulting schemes.

The traditional non conservative discretisation

Other forms of the ICE method used in the thermal hydraulics community ([98] section 11.2, Cathare, Sabena [START_REF] Ninokata | Sabena: Subassembly boiling evolution numerical analysis[END_REF]) discretise the non conservative version (1.12) of the Euler equations. The schemes take the generic form

         ρ n+1 i -ρ n i t + 1 x (ρ up i+ 1 2 u n+1 i+ 1 2 -ρ up i-1 2 u n+1 i-1 2 ) = 0 ρi+ 1 2 u n+1 i+ 1 2 -u n i+ 1 2 t + 1 2 (ρ∂ x u 2 ) i+ 1 2 + 1 x (p n+1 i+1 -p n+1 i ) = 0 . (1.15) with ρi+ 1 2 = ρ i + ρ i+1 2
For stability reasons, the mass flux ρu at the cell interfaces is generally defined using an upwind density ρ up (see [START_REF] Prosperetti | Computational methods for multiphase flow[END_REF] section 11.2) defined on faces as :

ρ up i+ 1 2 =    ρ i if u i+ 1 2 > 0 ρ i+1 if u i+ 1 2 ≤ 0
Here are some examples of expressions for the momentum convection term:

semi implicit discretisation (Cathare 3D module section 11.5 page 339 of [START_REF] Lavialle | Cathare 2 v2.5 3mod3.1 code: General description (partners version)[END_REF], [START_REF] Chandesris | 3d module of cathare : improvement of the momentum convective term discretization[END_REF]) In this case the density is explicit in the mass flux ρu, in the velocity evolution terms ρ∂ t u and the convection term ρu∂ x u. The scheme takes the form

         ρ n+1 i -ρ n i t + 1 x (ρ up,n i+ 1 2 u n+1 i+ 1 2 -ρ up,n i-1 2 u n+1 i-1 2 ) = 0 ρn i+ 1 2 u n+1 i+ 1 2 -u n i+ 1 2 t + 1 2 ρ up,n i+ 1 2 (∂ x u 2 ) i+ 1 2 + 1 x (p n+1 i+1 -p n+1 i ) = 0 , (1.16) with (∂ x u 2 ) i+ 1 2 =        1 x i+ 1 2 ((u 2 ) n+1 i+ 1 2 -(u 2 ) n i-1 2 ) if u i+ 1 2 > 0 1 x i+ 1 2 ((u 2 ) n+1 i+ 3 2 -(u 2 ) n i+ 1 2 ) if u i+ 1 2 ≤ 0 (1.17)
The Bernouilli principle is recovered at the discrete level in the case of stationary incompressible non viscous flows. Note that in 2D/3D, the density is upwinded in axial momentum terms of the form ρu 2 but centered in non axial terms of the form ρuv.

The treatment in (1.17) allows for a stronger implicitation of the velocity and larger time steps whilst retaining the elimination of the velocity unknowns. However this approach introduces a time consistency error as we detail in the following.

We apply a standard truncation error analysis to the discrete momentum equation in (1.16) assuming constant densities. This corresponds to the term:

∂u ∂t + u ∂u ∂x , in the case u i+ 1 2 > 0: u(x i+1/2 , t n+1 ) -u(x i+1/2 , t n ) ∆t + 1 2 u 2 (x i+1/2 , t n+1 ) -u 2 (x i-1/2 , t n ) ∆x = ∂u ∂t (x i+1/2 , t n ) + O(∆t) + 1 2 ∂ ∂x u 2 (x i+1/2 , t n ) + O(∆x) + ∆t ∆x u(x i+1/2 , t n ) ∂u ∂t (x i+1/2 , t n ) + O(∆t)
When the time step ∆t and the mesh size ∆x goes to zero, the consistency error tends to ∂u ∂t + u ∂u ∂x with an additional term in bold that prevent the scheme to be consistent with the momentum equation.

fully implicit discretisation (Cathare 1D module section 11.3 page 320 of [START_REF] Lavialle | Cathare 2 v2.5 3mod3.1 code: General description (partners version)[END_REF], [START_REF]Discrétisation des équations du module 3d de cathare 2 & comparaison avec le module 1d -etat des lieux pour la version v2.5 2 et perspectives[END_REF]) In this case the density is treated implicitly in the discrete mass flux ρu, in the discrete evolution terms ρ∂ t u and in the discrete convection term ρu∂ x u. The scheme takes the generic form

         ρ n+1 i -ρ n i t + 1 x (ρ up,n+1 i+ 1 2 u n+1 i+ 1 2 -ρ up,n+1 i-1 2 u n+1 i-1 2 ) = 0 ρn+1 i+ 1 2 u n+1 i+ 1 2 -u n i+ 1 2 t + ρn+1 i+ 1 2 (u∂ x u) i+ 1 2 + 1 x (p n+1 i+1 -p n+1 i ) = 0 , with (u∂ x u) i+ 1 2 =        1 x i+ 1 2 u n+1 i+ 1 2 (u n+1 i+ 1 2 -u n+1 i-1 2 ) if u n+1 i+ 1 2 > 0 1 x i+ 1 2 u n+1 i+ 1 2 (u n+1 i+ 3 2 -u n+1 i+ 1 2 ) if u n+1 i+ 1 2 ≤ 0
This treatment allows for a stronger implicitation of the velocity and larger time steps, whilst retaining the elimination of the velocity unknowns. The Bernouilli principle is not recovered at the discrete level in the case of stationary incompressible non viscous flows.

explicit discretisation (Sabena [START_REF] Ninokata | Sabena: Subassembly boiling evolution numerical analysis[END_REF]) In this case the density is explicit in the mass flux ρu, in the velocity evolution terms ρ∂ t u and the convection term ρu∂ x u. The scheme takes the form

         ρ n+1 i -ρ n i t + 1 x (ρ up,n i+ 1 2 u n+1 i+ 1 2 -ρ up,n i-1 2 u n+1 i-1 2 ) = 0 ρn i+ 1 2 u n+1 i+ 1 2 -u n i+ 1 2 t + ρn i+ 1 2 ( 1 2 ∂ x u 2 ) n i+ 1 2 + 1 x (p n+1 i+1 -p n+1 i ) = 0 , (1.18) with 1 2 (ρ∂ x u 2 ) i+ 1 2 =        1 2 ρ n i 1 x i+ 1 2 ((u 2 ) n i+ 1 2 -(u 2 ) n i-1 2 ) if u i+ 1 2 > 0 1 2 ρ n i+1 1 x i+ 1 2 ((u 2 ) n i+ 3 2 -(u 2 ) n i+ 1 2 ) if u i+ 1 2 ≤ 0
The Bernouilli principle is recovered at the discrete level in the case of stationary incompressible non viscous flows.

The explicit treatment introduces constraints on the time step.

explicit convection ([98] section 11.2.1, equation 11. [START_REF] Bal | A "parareal" time discretization for non-linear pde's with application to the pricing of an american put[END_REF] In this particular case the density is explicit in the mass flux ρu, in the velocity evolution terms ρ∂ t u and the convection term ρu∂ x u. The scheme takes the form

         ρ n+1 i -ρ n i t + 1 x (ρ up,n i+ 1 2 u n+1 i+ 1 2 -ρ up,n i-1 2 u n+1 i-1 2 ) = 0 ρn i+ 1 2 u n+1 i+ 1 2 -u n i+ 1 2 t + ρn i+ 1 2 (u∂ x u) n i+ 1 2 + 1 x (p n+1 i+1 -p n+1 i ) = 0 , with (ρu∂ x u) i+ 1 2 =          1 2 ρ n i + ρ n i+1 2 1 x i+ 1 2 u n i+ 1 2 (u n i+ 1 2 -u n i-1 2 ) if u n i+ 1 2 > 0 1 2 ρ n i + ρ n i+1 2 1 x i+ 1 2 u n i+ 1 2 (u n i+ 3 2 -u n i+ 1 2 ) if u n i+ 1 2 ≤ 0
With this treatment the velocity unknowns are easily eliminated using the momentum equation. This numerical scheme differs from the previous explicit discretisation (1.18) in the treatment of the convection term. In (1.18), the form 1 2 ∂ x u 2 is used and in the present scheme the form u∂ x u is discretised. The Bernouilli principle is not recovered at the discrete level in the case of stationary incompressible non viscous flows.

In most of the methods the velocity unknowns are first eliminated, and the resulting system is solved in a way that is compatible with the incompressible regime. The drawback of this approach is that there are constraints on the discretisation of the momentum flux ρ u ⊗ u and of the viscous terms µ u for the elimination of the velocity unknowns to be possible. The elimination can take place rigorously speaking when the convective flux ρ u ⊗ u and the viscous terms µ u are explicit in time. However explicit discretisations yield time step limitations.

The recent scheme of Herbin, Latché et al

In the past decade, Herbin, Latché and their coauthors have proposed a new approach with rigorous proofs of stability. They discretise the conservative form of the Euler equations (equation 1.11) with a conservative scheme. Their approach does not rely on velocity elimination and thus explicit and implicit variants are possible. The different variants include one step ([65] section 2.1, [66] section 3.1) and prediction/correction steps ([65] section 2.2, [66] section 4.1) variants, fully implicit ([66] section 3, [65] section 2.1), semi implicit and almost explicit [START_REF]On some implicit and semi-implicit staggered schemes for the shallow water and euler equations[END_REF] (all but the pressure gradient are explicit-in-time) variants.

For simplicity we present the discrete equation of the fully implicit variant ([66] section 3, [65] section 2.1) for the 1D isentropic Euler equations in conservative form:

ρ n+1 i -ρ n i t + ρ up,n+1 i+ 1 2 u n+1 i+ 1 2 -ρ up,n+1 i-1 2 u n+1 i-1 2 x = 0 (1.19) ρn+1 i+ 1 2 u n+1 i+ 1 2 -ρn i+ 1 2 u n i+ 1 2 t + ρu n+1 i+1 u up,n+1 i+1 -ρu n+1 i u up,n+1 i x + p n+1 i+1 -p n+1 i x = 0. (1.20)
The pressure p i and the density ρ i are located at the cell centers whereas the velocity u i+ 1 2 are located at the cell interfaces. The expression of the products ρu, ρ∂ t u and ρu 2 between the velocity located at cell interfaces and the density located at cell centers thus has to be defined through interpolation formula.

The mass flux ρu at the cell interfaces is defined using an upwind density ρ up i+ 1 2 defined as :

ρ up i+ 1 2 =    ρ i if u i+ 1 2 > 0 ρ i+1 if u i+ 1 2 ≤ 0 = ρ i + ρ i+1 2 + sign(u i+ 1 2 ) ρ i -ρ i+1 2 , (1.21)
which is the sum of a centered and an upwind terms.

The expression of ρi+ 1 2 in the discrete momentum equation accounts for an average of the neighbouring densities

ρi+ 1 2 = 1 2 (ρ i + ρ i+1 ). (1.22)
The expression of ρu in the discrete momentum equation is

ρu i = 1 2 (ρ up i-1 2 u i-1 2 + ρ up i+ 1 2 u i+ 1 2
).

(1.23)

The upwind velocity u up i at cell centers is defined as :

u up i =    u i-1 2 if ρu i > 0 u i+ 1 2 if ρu i ≤ 0 = u i-1 2 + u i+ 1 2 2 + sign(ρu i ) u i-1 2 -u i+ 1 2 2 , (1.24)
which is the sum of a centered and an upwind terms. It is possible to use a centered velocity ū instead of the upwind velocity u up (see [START_REF]On some implicit and semi-implicit staggered schemes for the shallow water and euler equations[END_REF]).

We presented the conservative staggered schemes studied in [START_REF] Herbin | Staggered schemes for all speed flows[END_REF], [START_REF]On some implicit and semi-implicit staggered schemes for the shallow water and euler equations[END_REF] that is proven to be entropic and to satisfy a kinetic energy preservation. In the next section, we introduce the one dimensional Cathare scheme in the two-phase flows configuration.

The one dimensional Cathare scheme

In this section, we first give some details on the one dimensional two-fluid model before giving the discrete equations of Cathare. We presented a generic two-fluid model in the section 1.1 with the system 1.1. In the sequel, we write the one dimensional Cathare model which is a variable cross section area model.

                   ∂Aα k ρ k ∂t + ∂Aα k ρ k u k ∂x = Γ k Aα k ρ k ∂u k ∂t + u k ∂u k ∂x + Aα k ∂p ∂x = Aα k ρ k g + F int k + F w k + Aσ M k u int A ∂ ∂t α k ρ k H k + u 2 k 2 + ∂ ∂x Aα k ρ k u k H k + u 2 k 2 = Aα k ∂p ∂t + Aα k ρ k u k g + Γ k H int k + (u int k ) 2 2 + Q w k + σ Q k (1.25)
where A is the cross section area of the channel. In the sequel, we specify the discretisation of some terms of the system (1.25):

In the mass balance laws, the void fraction (α up k ) i+1/2 and density (ρ up k ) i+1/2 at the cell edge i + 1/2 for the term:

∂Aα k ρ k u k ∂x i+1/2
are determined with an upwinding approach. (α up k ) i+1/2 is upwinded according to the velocity (u k ) i+1/2 and the density (ρ up k ) i+1/2 according to the quantity (

α up k ) i+1/2 (u k ) i+1/2 , (α up k ) i+1/2 = (α k ) i -α inf k , (u k ) i+1/2 > 0 (α k ) i+1 -α inf k , (u k ) i+1/2 < 0 (ρ up k ) i+1/2 = (ρ k ) i , (α up k ) i+1/2 (u k ) i+1/2 > 0 (ρ k ) i+1 , (α up k ) i+1/2 (u k ) i+1/2 < 0 (1.26)
where α inf k are constants for the minimal values of α k . The Cathare code considers that the volume fraction of each phase admit a non-zero minimum value and therefore, also a maximum value very slightly lower than one. It is therefore accepted that a flow is never completely single-phase, the residual phase represents respectively bubbles or tiny drops in a practically single-phase liquid or vapour zone.

The mass and energy balance equations in the one dimensional Cathare model (1.25) are written in a conservative form, whereas the momentum balance equation is written in a non conservative form. To evaluate (ρ k ) i+1/2 , ( ᾱk ) i+1/2 the cell-centered values at the edge i + 1/2 for the term:

Aα k ρ k ∂u k ∂t + u k ∂u k ∂x i+1/2 ,
an average of the adjacent values is taken weighted by the cell volumes.

(ρ k ) i+1/2 = V ol i ρ i + V ol i+1 ρ i+1 V ol i + V ol i+1
The convection term:

u k ∂u k ∂x i+1/2 ,
is determined with an upwinding approach,

u k ∂u k ∂x i+1/2 = (u k ) i+1/2 ((u k ) i+1/2 -(u k ) i-1/2 ), (u k ) i+1/2 > 0 (u k ) i+1/2 ((u k ) i+3/2 -(u k ) i+1/2 ), (u k ) i+1/2 < 0
The value of the cross section area is carefully chosen in the term

Aα k ρ k u k ∂u k ∂x
to ensure the Bernoulli principle for single phase flows with varying cross section.

A Bern,i+1/2 =          Āi+1/2 2 A i+1/2 + A i-1/2 A i-1/2 , (u k ) i+1/2 > 0 Āi+1/2 2 A i+1/2 + A i+3/2 A i+3/2 , (u k ) i+1/2 < 0 with Āi+1/2 = V ol i + V ol i+1 2∆x
In the energy balance laws, the velocity (ū k ) n+1 i at the cell i for the time derivative term:

∂ ∂t α k ρ k H k + u 2 k 2 i = (α k ) n+1 i (ρ k ) n+1 i (H k ) n+1 i + 1 2 ((ū k ) n+1 i ) 2 -(α k ) n i (ρ k ) n i (H k ) n i + 1 2 ((ū k ) n i ) 2 ∆t
is computed with an interpolation using the adjacent velocity at edges and the void fractions. In this interpolation, the void fractions are explicit in time

(ū k ) n+1 i = ((α k ) n i + (α k ) n i-1 )A i-1/2 (u k ) n+1 i-1/2 + ((α k ) n i + (α k ) n i+1 )A i+1/2 (u k ) n+1 i+1/2 ((α k ) n i + (α k ) n i-1 )A i-1/2 + ((α k ) n i + (α k ) n i+1 )A i+1/2
(1.27)

In the time derivative term of the energy equations the approximation involves the terms (ū k ) n+1 i and (ū k ) n i . From the interpolation formula (1.27), we see that the evaluation of (ū k ) n+1 i is made with the adjacent volume fractions at time T n and with the volume fractions at time T n-1 for the term (ū k ) n i . Hence, the approximation of the time derivative at time T n+1 in the energy equations is a function of the unknowns at times T n+1 , T n and T n-1 . From this reason, the Cathare time scheme is a two-step time scheme. This point will raise additional difficulties in the application of the time parallel algorithm, the parareal method. This requires an adaptation of the parareal algorithm to this type of multi-step time schemes that should not be intrusive in the Cathare code. We will detail this aspect in the next chapter.

All the terms are implicit in time in the system (1.25) except for the interpolated velocity at nodes in the energy equations.

Difficulties of two-phase flow models

Even if two-fluid models inherit achievements obtained in the single-phase flow modeling, mathematical theory and numerical methods, however, they possess many specific difficulties due to the existence of two phases in the same domain and their interactions. In this section we will discuss some difficulties in general existing in the two-phase flow models such as the presence of non conservative products, the configuration of the vanishing phase and the handling of discontinuous source terms. We illustrate the mathematical challenge in discretising non conservative products on the term α k ∂p ∂x appearing in the the momentum equations of the six-equation two-fluid model (1.1) (see section 1.3.3.1). Then we comment in section 1.3.3.2 how the Cathare code handles numerically vanishing phases. Finally in section 1.3.3.3, we present the challenges coming from the discretisation of stiff source terms.

Non conservative products in the two-fluid model

A theory of hyperbolic conservation laws, studied in depth in the literature, can be found in [START_REF] Serre | System of Conservation Laws 1 : Hyperbolicity, Entropies, Shock Waves[END_REF][START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF][START_REF]Finite Volume Methods for Hyperbolic Problems[END_REF][START_REF] Lefloch | Hyperbolic Systems of Conservation Laws, the Theory of Classical and Nonclassical Shock Waves[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balances scheme for sources[END_REF][START_REF] Bressan | Hyperbolic Systems of Conservation Laws : The One-dimentional Cauchy Problem[END_REF][START_REF] Dafermos | Hyperbolic conservation laws in continuum physics[END_REF][START_REF] Desprès | Systèmes hyperboliques de lois de conservation[END_REF][START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]. Such a theory gives a fundamental understanding and main ideas for plenty of numerical methods to solve a hyperbolic system of conservation laws, i.e. find a weak solution in the sense of distributions. However, our system (1.1) possesses non conservative products and is not therefore a conservative system. A discontinuous solution would lead to the product of two distributions that is not well-defined. In order to study the weak solutions of a non conservative hyperbolic system, one may consider different approaches.

In general, the most popular approach to deal with non conservative products is the theory of non conservative hyperbolic systems studied by Dal Maso et al in [START_REF] Maso | Definition and weak stability on non conservative products[END_REF]. In the classical theory of conservation laws, a shock wave depends merely on its left state and right state, the definition of Dal Maso et al depends on the choice of a specific path which connects a left state to a right state around a shock wave for non conservative products. Defining an appropriate path requires realistic physical information which is not easy, especially in a complicated two-fluid model. Moreover, once the appropriate path is chosen, different numerical methods may converge to different solutions, see comments in [4] and references therein. Therefore, one may prefer to solve the two-fluid model by using a simpler consideration of the non conservative product. For example, in [START_REF] Toumi | An Approximate Linearized Riemann Solver for a Two-Fluid Model[END_REF], the authors rewrite the 1D two-fluid model (1.2) and choose jump conditions based on a particular case where the system has a conservative form. The resulting two-fluid model includes two conservation laws of mass of each phase, one conservation law of mixture of momentum and one equation for the liquid velocity assumed incompressible. More precisely, it consists in the following equations:

∂α v ρ v ∂t + ∂α v ρ v u v ∂x = 0 (1.28) ∂α l ρ l ∂t + ∂α l ρ l u l ∂x = 0 (1.29) ∂α v ρ v u v + α l ρ l u l ∂t + ∂α v ρ v u 2 v + α l ρ l u 2 l + p ∂x = 0 (1.30) ∂u l ∂t + ∂ ∂x u 2 l 2 + p ρ l = 0. (1.31)
In [START_REF] Toumi | An Approximate Linearized Riemann Solver for a Two-Fluid Model[END_REF], they propose a formula to locally linearize the term α k in the product α k ∂p ∂x around a discontinuity. This local linearisation denoted αk is chosen such that the original isentropic model (1.2) (neglecting ∆p∂ x α k ) and the system (1.28-1.31) have the same Rankine-Hugoniot relation. After calculations, [START_REF] Toumi | An Approximate Linearized Riemann Solver for a Two-Fluid Model[END_REF] finds:

αl = 2α L l α R l α L l + α R l , αv = 1 -αl .
This formula is then applied to the simulation of two-fluid model, [START_REF] Toumi | Approximate Riemann solvers and flux vector splitting schemes for two-phase flow[END_REF]. It is important to note that although ∆p was neglected in the identification of jump conditions for the systems (1.2) and (1.28-1.31), its contribution is fundamental in practice to obtain real characteristic waves. On the other hand the non conservative term in the energy equations is written as a spacial derivative in [START_REF] Toumi | An Approximate Linearized Riemann Solver for a Two-Fluid Model[END_REF] using the assumption of incompressible liquiq together with the liquid mass equation:

∂ t α l = -∂ x (α l u l ).
Other authors focusing on the numerical methods usually neglect the product of p∂ t α k in the energy equations, see for example [START_REF] Munkejord | A musta scheme for a nonconservative two-fluid model[END_REF][START_REF] Morin | A roe scheme for a compressible sixequation two-fluid model[END_REF].

Cathare treatment of the vanishing phase

An elementary test case: the oscillating manometer:

In this section, we illustrate the role of the interfacial friction term on a two-phase test case where one of the phases disappears. The equations used are the simplified Cathare equations, (1.25).

Transfers between phases (mass and heat) and wall friction are neglected in the system (1.25).

The only source terms we consider in this test case are the gravity force and the interfacial forces

F int k .
In [START_REF] Hewitt | [END_REF], the oscillating manometer is proposed as a numerical benchmark test for system codes to test the ability of each numerical sheme to preserve system mass and to retain the gas-liquid interface. It consists in a U-shaped tube manometer which is connected at the top, so that a closed system is formed. The system contains initially gas and liquid with the liquid forming We seek to show the influence of the interfacial friction term on the behaviour of the scheme when one of the phases disappears in some parts of the domain. Previously in section (1.2.5), we saw that the interfacial friction coefficient C D goes to infinity in the configuration of a vanishing phase, as a consequence the system tends to an equilibrium of the velocities. The expression of the interfacial friction term depends on the flow regime (bubbly, annular, dispersed,...) and on the geometry:

C D = f (α k , ρ k , σ, µ k , D h )
where: σ is the surface tension, D h the hydraulic diameter and µ k are liquid and vapour dynamical viscosities.

The equation of motion of the vapour-liquid interface in each leg of the manometer is the following:

d 2 x dt 2 + 2gx L = 0 x(t = 0) = 0, dx dt (t = 0) = u 0 .
The solution to this equation is:

x(t) = 2u 0 L 2g sin 2g L t
and the liquid velocity is :

u L (t) = u 0 cos 2g L t , (1.32)
where g designates the acceleration of gravity and L designates the water length equal to 12 m.

For these simulations, the variation versus time of the liquid velocity at the bottom of the manometer (at one edge) is plotted. In Figure 1.2, we compare the reference solution of the oscillating manometer (1.32) with several numerical solutions computed with the one dimensionnal Cathare scheme. These results show that the numerical solution converge to another solution that is damped. This can be explained by the treatment of the vanishing phase used in Cathare. When the interfacial friction term dominates the system (section 1.2.5), it behaves as if a diffusion term was present in the model that damped the solution and prevent the numerical solution to converge towards the reference solution (1.32) of the test case.

An industrial test case: the Water-Packing

The treatment of the vanishing phase is a challenge for every software that simulates two-phase flows. In the context of nuclear safety studies, a well known problem linked to vanishing phases is the Water-Packing. It occurs during the simulation of a water level rise due to condensation as well as in other more complex situations. The typical situation in which this problem was studied is a vertical tube initially filled with superheated steam, connected to a steam tank at its top, is gradually filled with very cold liquid from below. As the liquid front progresses upwards, the vapour condenses and the resulting local depression aspirates steam from the tank and the liquid front gradually heats up with the condensation. Provided you know the rate of phase change by condensation, the analytical solution to this problem is simple with regard to the pressure at the bottom of the tube. Assuming that the front rises slowly, the problem is quasi-static and the pressure at the bottom of the tube then increases regularly due to the weight of the liquid column which increases by the addition of liquid injection and condensation. On the other hand, if one wishes to simulate this flow numerically, one runs up against difficulties: the field of physical pressure is disturbed by parasitic waves of strong amplitudes. This is the phenomenon of Water-Packing. The simulation of this kind of stratified flow is of great interest. In fact this type of situation is encountered in a more complex form in accidental transient calculations (for example, the filling phase of the reactor core after a large breach). The calculation of the water level rise in a simple situation is therefore used to analyse and solve a problem occurring in much more complex situations. This problem occurs in all the codes of thermalhydraulics dedicated to the simulation of accidental transients in the circuits of a nuclear reactor. It has been the subject of numerous studies and it is used to evaluate codes in test batteries called benchmarks. The Water-Packing benchmark was formalized by V. H. Ransom in 1987, in [START_REF] Hewitt | Multiphase Science and Technology[END_REF], under the title "expulsion of steam by cold water" (see Figure 1.3). To understand the origin of this phenomenon, we recall how the Cathare code Figure 1.3: Setting of the Water-Packing benchmark in [START_REF] Hewitt | Multiphase Science and Technology[END_REF] handles vanishing phase. The principle consists in considering that the volume fraction of each phase admit a non-zero minimum value and therefore, also a maximum value very slightly lower than one. It is therefore accepted that a flow (or a flow zone) is never completely single-phase, the residual phase represents respectively bubbles or tiny drops in a practically single-phase liquid or vapour zone. So even when a single-phase flow (or flow zone) is simulated, all the equations are solved. Also the mechanical equilibrium is forced: the residual phase has the same velocity as the dominant phase, consequence of an interfacial friction made artificially very large (section 1.2.5). This choice was made in the Cathare code since its initiation and has been valid for a large number of use cases. However it raises problems in the representative test case of Water-Packing. Physically, this choice amounts to say that the residual bubbles (resp. drops) in the liquid (resp. vapour) are entrained by the dominant liquid phase (resp. vapour).

In the context of test case (1.3), the liquid front is rising thus the liquid velocity is positive below the interface. The steam condenses at the interface with cold liquid, thereby creating an intake of additional steam from the top of the tube. The vapour velocity is negative when above the interface. The residual liquid velocity in the vapour phase is then negative. Therefore, when a cell close to the interface is filled up with liquid, liquid velocity suddenly changes direction. A vector node located above the interface at a negative liquid velocity will become positive when it is joined by the interface. Artificial pressure spikes are then observed whenever the liquid level fills up a cell (see Figure 1.5, red curve).

Here is an explanatory diagram of the behaviour of the liquid velocity during the advance of the liquid front(see Figure 1.4). Currently, the Cathare code uses an Anti-Water-Packing correction, specially dedicated to the treatment of this problem. However, it comes with a large number of shortcomings with respect to the balance laws and returns sometimes erroneous values. The principle consists in locating the cell(s) in which the interface is located, then greatly reduce the rate of condensation. As a consequence, the pressure shows artificial spikes with lower amplitude (see Figure 1.5, blue curve) because by decreasing the rate of condensation, the velocity at which the vapour is aspirated decreases too and the velocity of the residual liquid phase just above the interface decreases too. The inversion of the liquid velocity still takes place but it goes from a weakly negative velocity to a stronger positive one. In summary, the inversion of the liquid velocity occurs more gradually. The upwinding approach in the Cathare numerical scheme plays an important role in solving this problem. Indeed, during the simulation of a counter-current flow, the upwinding strategy can be inconsistent with the dynamics we want to capture. Each momentum equation is upwinded according to the sign of its phasic velocity (reference section schema de Cathare). This approach is consistent when bot phasic velocities have the same sign but may fail to capture the void waves when the phasic velocities have opposite signs. In [START_REF] Ndjinga | A 2 × 2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF], the authors propose colocated schemes to capture void waves without forcing a mechanical equilibrium with a large interfacial friction. This scheme is based on upwinding according to the volume fraction wave speed that is different from the phasic velocities.

Discontinuous source terms

In the thermalhydraulics of nuclear reactors, two-fluid models display stiff source terms S(U, x). The stiffness of these source terms has different origins. First, the heat source Φ is localised on the Figure 1.5: Pressures spikes with or without Cathare anti-Water-Packing correction, from [START_REF] Mohaupt | Étude du problème numérique de Water-Packing[END_REF] core of the reactor which yields a discontinuity in space. Moreover, the dry-out of the Uranium rod when the temperature reaches a threshold (the critical heat flux) yields a discontinuity of Φ as a function of the temperature. Secondly, the boiling of the fluid is a stiff phenomena with a threshold that is the boiling temperature (or enthalpy). For these reasons the simulation of boiling of two-phase flows is challenging both from a mathematical and a numerical point of view. The source term S(U, x) being discontinuous both in U and x makes it impossible to use Cauchy-Lipschitz type of theorems for the existence of solutions even for ODEs. However, there are particular cases where a unique solution may exist (see [START_REF] Bressan | Unique solutions for a class of discontinuous differential equations[END_REF]). The source term S(U, x) being discontinuous in U makes numerical approximation more difficult in the numerical simulation because of the stiffness of the solution. Classical approaches to deal with stiff source term assume it is Lipschitz in the variable U .

In this section, we introduce non-homogeneous hyperbolic systems of conservation laws, i.e. we take into account a non zero source term S as below

∂U ∂t + ∂ x F (U ) = S(U, x), x ∈ R, t ≥ 0, (1.33)
The source term S is usually a function of the unknown vector U and spacial variable x, S(U, x).

A simple approach to solve the non homogeneous system of conservation laws is to include a source term in the right hand side as follows:

U n+1 i -U n i ∆t + F i+1/2 -F i-1/2 ∆x = S n i , (1.34)
where S n i is an approximation of

S i = 1 ∆x 1 ∆t t n+1 t n x i+1 x i Sdxdt (1.35)
and F i±1/2 is some interfacial flux function.

The numerical scheme (1.34) is a classical one solving the non homogeneous system of conservation laws. Sometimes both the flux divergence ∂ x F (U ) and the source term S(U, x) are discretised independently and in practice, S n i is simply considered as the source function at the average value U i , i.e. S n i = S(U n i , x i ). However this approach may generate instabilities in the simulation of the system (1.33), especially for stiff source terms functions S. We are interested in the capture of the stationary regime of a two-phase flow characterised by the stationary equation:

∂ x F (U ) = S. (1.36)
In some cases, for example with stiff source terms, using S n i = S(U n i , x i ) in (1.35) implies the instability of the numerical solution at the stationary state. In order to improve the numerical simulation, one suggests either upwinding the source terms ( [START_REF] Bermudez | Upwind method for hyperbolic conservation laws with source terms[END_REF]) or developping well-balanced schemes ( [START_REF]Nonlinear stability of finite volume methods for hyperbolic conservation laws and wellbalances scheme for sources[END_REF]) in the sense that it preserves the stationary state. In a two-fluid model, due to the complexity of the flux function and the lack of regularity of source terms, it seems difficult to construct a well-balanced scheme.

Acceleration techniques for the simulation of two phase flows

In the numerical resolution of partial differential equations, we generally have to solve linear or nonlinear systems arising from the discretisation. The large size of these systems and the fact that they are ill conditionned make a global resolution difficult. In this section, we firstly summarise in 1.4.1 the resolution method used by the Cathare code to solve the nonlinear system arising from the discretisation of the two-fluid model. Then we present the actual acceleration methods available in the Cathare code in 1.4.2. Finally, we introduce in 1.4.3 the time parallel algorithm that we consider in this thesis, the parareal algorithm.

Solution algorithm for the two-fluid model

After discretisation in space and time of the two-fluid model, we obtain a non linear system that is solved with a Newton method. Before describing the implementation of this method, we recall here the dependance of the different equations according to the implicited variables: the mass and energy balance equations at cell i depend on the pressure, void fraction and enthalpies defined at cell i and on the velocities defined at the edges. the momentum balance equations at edge i + 1/2 depend on the velocities defined at edges i + 1/2, i -1/2 and i + 3/2 and on the cell centered unknowns defined at the adjacent cells of the edge i + 1/2.

After applying the semi-implicit time scheme of the Cathare code, we obtain a system of the following form:

U n+1 -U n ∆t + A(U n+1 , U n ) = S(U n )
This non linear system is solved by a Newton method:

δU k+1 ∆t + J(U k , U n )δU k+1 = S(U n , U k ), where: δU k+1 = U k+1 -U k and : U k+1 = (p k+1 , α k+1 v , H k+1 l , H k+1 v , u k+1 l , u k+1 v ) (1.37)
The increments of the principal variables in the Cathare code are denoted δU k+1 and the terms of the Jacobian matrix J(U k , U n ) are computed as follows:

J i,j = ∂E i U k j
Since the Cathare time discretisation is semi-implicit, this term corresponds to the derivative of the balance equation E i with respect to the principal variable U k j that is implicited, for the kth Newton iteration. Hence, the Cathare code uses the derivatives with respect to the variables (p,

α v , H l , H v , u l , u v ) instead of the conservative variables (α v ρ v , α l ρ l , u v , u l , α v ρ v E v , α l ρ l E l ).
The resulting linear system (2.1) at each Newton iteration is solved by a Gauss elimination to obtain a system with pressure increments only. Another specificity of the Cathare code is that it uses a semi-implicit pressure solver without splitting techniques for the velocity (no prediction/correction steps), and no explicit construction of an elliptic problem on the pressure.

In the Cathare software, the junction is an element of the geometry that links two other elements of the geometry: a 1D element (a pipe, ...) with a 3D element (vessel, ...) or a 1D/3D element with a boundary condition element (see Figure 1.6). The junctions store the informations about the edge that is common to the two elements linked by this junction (two velocities and derivatives of the momentum equations associated to this edge with respect to the principal unknowns). Hence the system depends on the N i pressure increments belonging to the internal cells and on the N j pressure increments belonging to the junctions. The unknowns associated to the junctions are called external variables. The dependance on the external pressure increments appears when we write the momentum balance equations associated to the junctions. The resulting system is of the following form:

A 11 A 12 A 21 A 22 ∆p(N i ) ∆p(N j ) = S 1 S 2 (1.38)
The method adopted in the Cathare code consists in eliminating the internal pressure increments ∆p(N i ) and then solve the problem on the external pressure increments ∆p(N j ):

A 22 -A 21 A -1 11 A 12 ∆p(N j ) = S 2 -A 21 A -1
11 S 1 This linear system is then solved by a direct method with the library LAPACK BLAS.

Actual acceleration methods in Cathare

In [START_REF] Ruby | Quest for the real-time for the safety analysis code cathare 2 used in the post-accident simulator sipa[END_REF], the authors present the strategy implemented in the Cathare code in order to make its use compatible with a real-time response. This work was carried out within the context of the integration of the Cathare code in the SIPA simulator (Post-Accident Simulator) for training and engineering studies on nuclear PWR reactors under nominal conditions or accidental transients. In the sequel, we list the steps of the solution algorithm described in section 1.4.1, that were parallelised:

Each element of the Cathare geometry (see figure 1.6) assembles a block of the Jacobian matrix in parallel. Hence, for each element, the matrix block depends on the N i pressure increments belonging to the internal cells of the element and on the N j pressure increments belonging to the junctions associated to the element (step (i) in figure 1.7). An element of the Cathare geometry can be seen as a subdomain of a spacial domain decomposition method.

Each element eliminates the internal variables of system (1.38) and obtains a system depending on the junction variables ∆p(N j ) only (step (ii) in figure 1.7).

The pressure increments ∆p(N j ) for all the junctions of the Cathare geometry are computed by a Gauss elimination. It is performed by an iterative algorithm that successively eliminates the blocks of the Jacobian matrix corresponding to a junction common to two elements. This process goes on, as long as the non-eliminated junctions are common to at least one other element. The order used to eliminate the block matrices is called an elimination tree. It depends on the reactor meshing and is optimised by the Cathare code for a sequential resolution (step (iii) in figure 1.7).

The increments of the other principal variables (velocities, volume fractions and enthalpies) are computed in parallel over all the elements (step (iv) in figure 1.7).

In figure 1.7, we see the performances of the actual parallel method implemented in the Cathare code. We see in the left figure that the computational cost of the step (i) in the solution algorithm represents 65% of the global computational time of the simulation. Hence, the parallelisation of step (i) over the elements of the geometry offers good speed up performances: a speed up of 10 is obtained for the parallelisation of step (i) on 12 processors where the global speed up is 6. This actual acceleration method has been successfully applied on industrial test cases for the simulation of a nuclear reactor under an accidental scenario. However, these parallel methods are implemented with OpenMP that allows shared memory parallelism. Their use is limited to a level of parallelism of about 20 processors or of at most all the processors of a standard desk computer. Moreover, the performances of this strategy reaches its limit when the geometry of the simulation includes 1D elements with ∼ 10 2 cells and 3D elements with ∼ 10 3 . This imbalance of tasks between the processors can damage the speed up performances. We illustrate this behaviour in the left figure 1.8 where the computational time of one step of the Jacobian assembling is measured for one Newton iteration in one time step. An optimisation of the parallel method that we have previously described has been proposed to overcome the issue of tasks imbalance between processors. This new algorithm allows to assign elements of the geometry to the processors with load balancing, knowing the computational time of the previous time step. In the right figure 1.8, we see that this algorithm allows to significantly improve the load balancing between processors and thus the performances of the parallel method.

Time domain decomposition: the parareal algorithm

The scalability properties of the space domain decomposition method implemented in the Cathare code are limited by the small number of cells in the meshes of the Cathare simulations. Several works are in progress to propose novel parallel algorithms in the Cathare development group. For example, the possibility to assign a 3D element to many threads by dividing the element on subdomains is actually investigated. Since the number of elements in a Cathare simulation is usually smaller than the number of available processors, we seek in this work to investigate a novel strategy of parallelisation to complement the actual parallelism in the Cathare code. For this reason, if we have more processors at our disposal and wish additional speed-ups, the parallelisation of other variables needs to be addressed. Several approaches have been proposed over the years to decompose the time direction when solving a partial differential equation (see [START_REF] Nievergelt | Parallel methods for integrating ordinary differential equations[END_REF], [START_REF] Chartier | A parallel shooting technique for solving dissipative ODE's[END_REF], [START_REF] Erhel | Algorithme parallèle pour le calcul d'orbites -Parallélisation à travers le temps[END_REF], [START_REF] Emmet | Toward an efficient parallel in time method for partial differential equations[END_REF], [START_REF] Gander | Paraexp : A parallel integrator for linear initial-value problems[END_REF], [START_REF] Christlieb | Parallel high-order integrators[END_REF], [START_REF] Falgout | Parallel time integration with multigrid[END_REF], [START_REF] Gander | 50 years of time parallel time integration[END_REF]). Of these, the parareal algorithm, whose performances we explore in this work, was proposed two decades ago by [START_REF] Lions | Résolution par un schéma en temps "pararéel[END_REF] and has received an increasing amount of attention in the last years. The parareal method can also be cast into the category of multiple shooting type methods that were introduced in [START_REF] Nievergelt | Parallel methods for integrating ordinary differential equations[END_REF] (we refer to [START_REF] Gander | Analysis of the parareal time-parallel time-integration method[END_REF] for a detailed discussion about the several possible interpretations of the parareal method). The parareal method has been applied successfully to a number of applications (see [START_REF] Baffico | Parallel-in-time molecular-dynamics simulations[END_REF], [START_REF] Fischer | A parareal in time semi-implicit approximation of the Navier-Stokes equations[END_REF], [START_REF] Samaddar | Parallelization in time of numerical simulations of fully-developed plasma turbulence using the parareal algorithm[END_REF], [START_REF] Lunet | Time-parallel simulation of the decay of homogeneous turbulence using Parareal with spatial coarsening[END_REF], [START_REF] Samaey | A micro/macro parallel-in-time (parareal) algorithm applied to a climate model with discontinuous non-monotone coefficients and oscillatory forcing[END_REF] among many others), demonstrating its versatility. Theoretical advances on this method include stability analysis ( [START_REF] Bal | A "parareal" time discretization for non-linear pde's with application to the pricing of an american put[END_REF], [START_REF] Staff | Stability of the parareal algorithm[END_REF], [START_REF]On the convergence and the stability of the Parareal algorithm to solve partial differential equations[END_REF], [START_REF] Dai | Stable Parareal in time method for first-and second-order hyperbolic systems[END_REF]), its coupling with spatial domain decomposition methods ( [START_REF] Maday | The Parareal in Time Iterative Solver: a Further Direction to Parallel Implementation[END_REF], [START_REF] Guetat | Méthode de parallélisation en temps: Application aux méthodes de décomposition de domaine[END_REF]) and control problems ( [START_REF] Maday | Monotonic time-discretized schemes in quantum control[END_REF], [START_REF] Maday | The Parareal in Time Iterative Solver: a Further Direction to Parallel Implementation[END_REF]). To see how the method works and how it has been applied to the two-fluid model, we write the system (1.1), after the discretisation in space on N degrees of freedom:

∂U ∂t + A(t, U ) = 0, t ∈ [0, T ], U (t = 0) = U 0 (1.39) A : R × R N → R N , U ∈ R N (1.40)
where U is the unknown, in our case, U = (p, α v , h l , h v , u v , u l ). Let us assume that we have two propagators G and F to solve (1.39). For any given t ∈ [0, T ], s ∈ [0, Tt] and any function w in a Banach space, G(t, s, w) (respectively F (t, s, w)) takes w as an initial value at time t and propagates it at time t + s.

G(T n , ∆T, U n ) computes a coarse approximation of U (T n+1 ) with initial condition U (T n ) U n (low computational cost) F (T n , ∆T, U n ) computes a more accurate approximation of U (T n+1 ) with initial condition U (T n ) U n (high computational cost)
The fine propagator F can perform the propagation of the phenomenon with small time steps δt and with very accurate physics described by A. On the other hand, the coarse approximation G does not need to be as accurate as F and can be chosen much less expensive, by the use of a scheme with a much larger time step ∆T δt or by treating reduced physics. In addition to these two propagators F and G, the parareal in time algorithm is based on the division of the time interval [0, T ] into N sub-intervals [0, T ] = ∪ N -1 n=0 [T n , T n+1 ] that will each be assigned to a processor P n . The parareal algorithm applied to (1.39) is an iterative technique where, at each iteration k, the value U (T n ) is approximated by U n k with an accuracy that tends to the one achieved by the fine solver when k increases. U n k is obtained by the recurrence relation:

U n+1 k+1 = G(T n , ∆T, U n k+1 ) + F (T n , ∆T, U n k ) -G(T n , ∆T, U n k ) (1.41)
Starting from U n+1 0 = G(T n , ∆T, U n 0 ). From formula (1.41), one can see by recursion that the method is exact after enough iterations. Indeed, for any n > 0, U n n = F (T 0 , T n -T 0 , U 0 ). However, convergence of U n k to F (T 0 , T n -T 0 , U 0 ) goes much faster than this. While the main results about the convergence properties of the method were studied in depth several years ago (see [START_REF] Lions | Résolution par un schéma en temps "pararéel[END_REF], [START_REF] Baffico | Parallel-in-time molecular-dynamics simulations[END_REF], [START_REF] Bal | A "parareal" time discretization for non-linear pde's with application to the pricing of an american put[END_REF]), more recent efforts ( [START_REF] Minion | A hybrid parareal spectral deferred corrections method[END_REF], [START_REF] Aubanel | Scheduling of tasks in the parareal algorithm[END_REF], [START_REF] Emmet | Toward an efficient parallel in time method for partial differential equations[END_REF], [START_REF] Berry | Event-based parareal: A data-flow based implementation of parareal[END_REF], [START_REF] Maday | An adaptive parareal algorithm[END_REF]) focus on the algorithmics to implement it in order to improve the speed-up provided by the original algorithm.

Time domain decomposition for hyperbolic problems

An instability of the parareal algorithm may appear when it is applied to convection dominated problems. In [START_REF]On the convergence and the stability of the Parareal algorithm to solve partial differential equations[END_REF], the author points out the need of a strong regularity on the initial condition to ensure the stability of the parareal algorithm. When the problem is parabolic, the smoothing character of the problem prevents the appearance of instabilities even if the initial condition is not regular enough but when the problem is hyperbolic, an initial condition with high frequencies components may trigger the instability. Moreover in [START_REF]On the convergence and the stability of the Parareal algorithm to solve partial differential equations[END_REF], the author shows the influence of the numerical diffusion on this instability. This aspect has also been explored in [START_REF] Ruprecht | Convergence of Parareal with spatial coarsening[END_REF], [START_REF] Fischer | A parareal in time semi-implicit approximation of the Navier-Stokes equations[END_REF], [START_REF] Lunet | Time-parallel simulation of the decay of homogeneous turbulence using Parareal with spatial coarsening[END_REF], [START_REF] Ruprecht | Wave propagation characteristics of Parareal[END_REF].

Others strategies have been studied to treat the instability of the original parareal. They propose to improve the coarse approximations at every parareal iteration using the previous fine solutions to overcome the instability issues with hyperbolic equations. In [START_REF] Farhat | Time-parallel implicit integrators for the near-real-time prediction of linear structural dynamic responses[END_REF] and [START_REF] Gander | Analysis of a Krylov subspace enhanced Parareal algorithm for linear problems[END_REF], the parareal solution is projected in a Krylov subspace generated by the set of fine solutions from the previous parareal iterations. In [START_REF] Chen | On the use of reduced basis methods to accelerate and stabilize the Parareal method[END_REF], the authors use a subspace thanks to a reduced basis built from the matrix made of the previous fine solutions. In the same spirit, in [START_REF] Dai | Stable Parareal in time method for first-and second-order hyperbolic systems[END_REF], the authors propose a parareal algorithm with an additional step that project the parareal solution in a manifold to ensure the conservation of invariants (for example, the conservation of the Hamiltonian).

Other contributions propose cures to the parareal algorithm with an algebraic viewpoint. In [START_REF] Chen | An adjoint approach for stabilizing the Parareal method, Rapport technique[END_REF], the authors formulate the method with an iteration matrix and propose preconditionners to improve the behaviour of the parareal algorithm.

In [START_REF] Gander | Analysis of the Parareal algorithm applied to hyperbolic problems using characteristics[END_REF], the author shows the difficulties of the parareal algorithm to converge in a reasonable number of iterations on the advection equation and the wave equation using the method of characteristics.

In [START_REF] Nielsen | Communication-aware adaptive parareal with application to a nonlinear hyperbolic system of partial differential equations[END_REF], the authors propose an algorithm called Communication Aware Adaptive Parareal (CAAP) to speed up the non linear shallow water equation beyond what is possible using spatial domain decomposition methods alone.

We will also note that the application of other parallel in time algorithms is not straightforward on convection dominated problems. In [START_REF] Sterck | Parallel-in-time multigrid with adaptive spatial coarsening for the linear advection and inviscid Burgers equations[END_REF], the authors propose MGRIT algorithms that improve stability and scalability for the resolution of the advection equation with 1st-order numerical scheme. In [START_REF] Sterck | Convergence analysis for parallel-in-time solution of hyperbolic systems[END_REF], analysis tools are proposed to understand the source of the instability issues on MGRIT algorithms and on the parareal algorithm using the multigrid interpretation. In [START_REF] Krzysik | On selecting coarse-grid operators for Parareal and MGRIT applied to linear advection[END_REF], a convergence analysis is conducted to propose criteria for coarse-grid operators involved in the MGRIT and parareal algorithms. This strategy can ensure stability and scalability for the time parallelisation of the advection equations with high order discretisations. In this chapter, we present the two strategies we developed during the PhD in order to apply the parareal algorithm to the Cathare code: a numerical clone of Cathare that is restricted to one test case and a library that uses the Cathare code in a non intrusive way. The main contribution of this work has been to adapt the parareal algorithm to the architecture of the software and to its time discretisation in a non intrusive way, without any changes of the source files of the Cathare code, in order to reduce the computational time and get closer to a real-time response of the code. In section 2.1, we present the challenges of implementing the parareal algorithm on an industrial software, the Cathare code, in a non intrusive way. In section 2.2, we introduce the new algorithm we designed to handle multi-step time schemes such as the one used within the Cathare code. In the last section, we report the speed up performances obtained on two test cases: the oscillating manometer using the numerical clone of the Cathare code and the Omega test case using the Parareal library that uses the Cathare software. Each test case is representative of the numerical challenges for the simulation of two phase flows in the context of safety studies such as vanishing phases and discontinuous source terms.

The Parareal library for the Cathare code

The development of the Cathare code started 30 years ago for the simulation of nuclear reactors under nominal or accidental situations. The software was designed for engineers in nuclear energy and experts in the thermalhydraulics of accidental scenarii in nuclear power plants. Hence, the terminology used in Cathare is directly linked to the physics of nuclear reactors. For example, to define the geometry of the simulation, the user has to define a reactor, a primary and a secondary circuit and hydraulic elements representing pumps or steam generators. Generally, the definition of the simulation one seeks to run is composed of two main blocks: First, informations are given for the description of the circuits. Each element is defined with the reference to a hydraulic element or Cathare element (1-D, 0-D, 3-D, junction, etc...). Geometrical parameters and meshing are defined. Elements are connected to constitute elementary circuits (one primary and several secondary circuits). Heat exchangers between elementary circuits are defined.

The characteristics of the calculation are then specified. Successive directives are given corresponding to the different steps of the simulation:

operation of the initial state process, -time propagation called transient calculation, -time step control, -events occuring during the calculation (safety injection, break opening, valves, ...).

The actual parallel methods available in the Cathare code (see section 1.4 chapter 1) are implemented with OpenMP that allows shared memory parallelism. Their use is limited to a level of parallelism of about 20 processors or of at most all the processors of a standard desk computer. In this work, all the developments of the parareal algorithm applied to the Cathare code were made with MPI that allows distributed memory parallelism. Hence, the parareal algorithm allows to run a Cathare simulation in parallel on many processors located in different computers, including supercomputers.

In the sequel, we summarise the obstructions linked to the Cathare structure and the strategies we chose to handle them. These cures were of two natures: computer science and algorithmic. We think that the experience of implementing a parallel algorithm to an industrial software in a non intrusive way will be instructive for the future developments of parallel techniques in the Cathare code. This is why we list the obstructions that are only of a computer science nature. In the sequel, we list the adaptations and adjustments necessary to develop a library that applies the parareal algorithm to the Cathare code where it is used as a black box without modifying the source files of the software. The data structure in the Cathare code mainly lies on an array of values containing the informations related to the mesh, the unknowns, the coefficients of the Jacobian matrix and of the right hand side. This array stores the principal variables of pressure, volume fraction and liquid and vapour velocities and enthalpies for every cell of the mesh. The Cathare array also contains auxiliary variables that are necessary to carry on the simulation. These auxiliary values are of different natures: variables computed from the principal unknowns variables tracked during the simulation: for example, the water level in an element of the system boolean variables giving the flow regime in an area of the system: bubbly, annular, dispersed, seperated phases, etc.

The first obstruction we met is that we can only define one system (one reactor and one array of values) in a Cathare simulation while in the context of the parareal algorithm we use two systems: a coarse one and a fine one. The first difficulty is to run a simulation with coarse and fine systems in the same processor. There are two possible approaches:

When running the simulation over N processors, we assign one processor ("the coarse processor") to the coarse propagation for all time windows [T n , T n+1 ] and the remaining processors ("the fine processors") run the fine propagation only, over one time window. The coarse processor will also transmit the coarse approximation at times T n to the corresponding fine processor. The main drawback of this approach is the multiplication of communications and the use of N + 1 ressources for a parallelism over N time windows.

We use the same array of values for the coarse and the fine propagations in each processor. This is made possible by the use of the Parareal library that is independent of the Cathare code and contains the parareal algorithm. This library collects the coarse and fine propagations made by the Cathare code then apply the parareal corrections and finally send to the Cathare code the updated initial conditions for each processor [T n , T n+1 ]. Depending on whether the Cathare will make a coarse or a fine propagation, the Parareal library will send the data array associated to the coarse or the fine solver (for example, the coarse or fine time step and the suitable initial condition). For the implementation of the parareal algorithm applied to the Cathare in a non intrusive way, we chose this option.

The data structure of the Cathare code represents a challenge for memory storage with up to 10 6 variables and a size of about 8 Go. To reduce the cost of communication, we choose to exchange only the principal variables between processors. Some of the auxiliary variables are then computed from the principal unknowns. This choice will have consequences on the accuracy due to the error made on the reconstruction of the auxiliary variables. This question is addressed in section 2.3.2.

Data set

A Cathare data set can be split in different blocks of instructions: definition of the geometry and the mesh calculation of an initial condition: computation of a stationary state starting from the instructions of the user: pressures and temperatures for some specific cells, for example, or rotating speed of a pump or flow direction a loop for the time integration from T 0 to T f

The first two points are called "Initialisation" step and the last one "Time integration" step. For the application of the parareal algorithm, we split the Cathare data file into two files for the "Initialisation" step and for the "Time integration" step. Each processor P n has an initialisation file and a file for the time integration between T n and T n+1 . The initialisation file of processor P 0 is particular since it builds the geometry and mesh, calculates the initial state of the simulation and sends them to the other processors. The initial state will not be used in the remaining processors but this communication allows to initialise the size of the Cathare arrays in each processor. After that, P 0 makes a coarse propagation and sends G(T 0 , T 1 , U 0 0 ) to P 1 . The other processors receive the informations for the geometry and mesh and their local initial condition G(T n-1 , T n , U n-1 0 ) to run their own coarse propagation. We met two obstructions at this stage:

The calculation in [T 1 , T 2 ] and generally for [T n , T n+1 ], with n = 0, has to start from an initial time T n different from 0. It is unusual for the Cathare code to start a simulation from a time T n = 0. The consequence is: instead of beginning the calculation with the coarse solution sent from the previous time window as an initial condition, the processor uses the initial state sent from rank 0 processor. To treat this difficulty, the time window makes one time step and delete the computed solution. This allows to fix the initial time T n of the time window to the correct value.

The size of the arrays for every processors is computed at the moment processor P 0 sends to P i , i = 1, • • • , N -1 the geometry, mesh and initial state. If some variables are declared between the initialisation step and the time loop then the array in processor P 0 will have the correct size but the arrays in the remaining processors will still have the size of the state sent by P 0 during the initialisation step. Hence, the arrays in processors P i , i = 1, • • • , N -1 have a different size from the one in processor P 0 . This will generate conflicts during the communications between processors and the copy of arrays. Thus we need to adapt the "Initialisation" data file of processor P 0 to handle this.

Obstructions linked to the time discretisation

The Cathare time discretisation is based on a two-step time scheme (see Chapter 1, section 1.3.2 on the Cathare discretisation methods). This leads us to design a new variant of the parareal algorithm that takes more efficiently into account the presence of multi-step time schemes. We detail this aspect in section 2.2. Moreover, the Cathare time scheme uses an adaptive time step.

In a Cathare simulation, the user imposes an initial time step ∆t 0 and a maximal time step ∆t max . At each iteration, the code gets closer to ∆t max by multiplying by 2 the current time step while respecting the CFL condition integrated in the code. If in a time iteration, the Newton method does not converge after the maximal number of iterations, the current time step is reduced by multiplying it by 2 3 . Thus, we need to communicate the fine and coarse time steps between processors. For the coarse propagation, we can transfer the time step computed at the current parareal iteration because it is a sequential step. Since the fine propagation is made in parallel, the time window [T n , T n+1 ] at the k-th parareal iteration uses the fine time step from the time window [T n-1 , T n ] computed at the (k -1)-th parareal iteration.

A numerical clone of the Cathare code

In this section, we present the numerical clone of the Cathare code, called MiniCathare, implemented during the PhD. This tool allowed us to make a first trial on the efficiency of the parareal algorithm to speed up a two-phase test case using the Cathare model and numerical scheme. The test case we consider here is the oscillating manometer (see section 2.3.1 for the description of the test case and the numerical results obtained with the numerical clone). Hence, MiniCathare is restricted to one test case. We give in the sequel some details about this numerical tool: MiniCathare only makes the time propagation in a time intervall [T 0 , T f ] and does not compute an initial state like the Cathare code. The initial condition in MiniCathare is extracted from the Cathare code after the initial state process.

MiniCathare is implemented in C++ and the parareal algorithm in MPI.

MiniCathare allows to make convergence tests which is not possible with the Cathare code since there is a limitation on the mesh size in Cathare. In Chapter 1, section 1.3.3.2, we show the convergence properties of the Cathare numerical scheme on the oscillating manometer by comparing several numerical solutions, obtained by MiniCathare, to the reference solution of this test case.

The only difference between MiniCathare and the Cathare code is the solution of the non linear problem at each time step:

U n+1 -U n ∆t + A(U n+1 , U n ) = S(U n )
In both cases, this non linear system is solved by a Newton method:

δU k+1 ∆t + J(U k , U n )δU k+1 = S(U n , U k )
, where:

δU k+1 = U k+1 -U k and : U k+1 = (p k+1 , α k+1 v , H k+1 l , H k+1 v , u k+1 l , u k+1 v ) (2.1)
The increments of the principal variables in Cathare are denoted δU k+1 and the terms of the Jacobian matrix J(U k , U n ) are computed as follows:

J i,j = ∂E i U k j
Since the Cathare time discretisation is semi-implicit, this term corresponds to the derivative of the balance equation E i with respect to the principal variable U k j that is implicited for the k-th Newton iteration. The resulting linear system (2.1) is solved differently in Cathare and MiniCathare:

-Cathare, by a Gauss elimination, obtains a system with pressure increments δp k+1 only and then solves the problem in pressure with a direct linear solver (Lapack-Blas)

-MiniCathare assembles the whole Jacobian matrix and solves the complete linear system with an iterative linear solver (PETSC library)

Multi-step variant of the parareal algorithm

Several approaches have been proposed over the years to decompose the time direction when solving a partial differential equation (see [START_REF] Gander | 50 years of time parallel time integration[END_REF] for an overview). Of these, the parareal in time algorithm, which performances we explore in this work, has received an increasing amount of attention in the last twenty years with many applications (see [START_REF] Baffico | Parallel-in-time molecular-dynamics simulations[END_REF], [START_REF] Fischer | A parareal in time semi-implicit approximation of the Navier-Stokes equations[END_REF], [START_REF] Samaddar | Parallelization in time of numerical simulations of fully-developed plasma turbulence using the parareal algorithm[END_REF] among many others). In the sequel, we recall the classical parareal algorithm ( [START_REF] Lions | Résolution par un schéma en temps "pararéel[END_REF], [START_REF] Baffico | Parallel-in-time molecular-dynamics simulations[END_REF], [START_REF] Bal | A "parareal" time discretization for non-linear pde's with application to the pricing of an american put[END_REF]) and present the multi-step variant we will apply to the Cathare code.

Original parareal algorithm and notations

After the discretisation of a PDE in space, we obtain an ODE system of the form:

∂u ∂t + A(t, u) = 0, t ∈ [0, T ], u(t = 0) = u 0 (2.2)
where A : R × R N → R N , and N denotes the number of degrees of freedom.

Let G and F be two propagators such that, for any given t ∈ [0, T ], s ∈ [0, Tt] and any function w in a Banach space, G(t, s, w) (respectively F (t, s, w)) takes w as an initial value at time t and propagates it at time t + s. The full time interval is divided into N c sub-intervals [T n , T n+1 ] of size ∆T that will each be assigned to a processor. The algorithm is defined using two propagation operators:

G(T n , ∆T, u n ) computes a coarse approximation of u(T n+1 ) with initial condition u(T n ) u n (low computational cost) F (T n , ∆T, u n ) computes a more accurate approximation of u(T n+1 ) with initial condition u(T n ) u n (high computational cost)
Starting from a coarse approximation u n 0 at times T 0 , T 1 , • • • , T N c , obtained using G, the parareal algorithm performs for k = 0, 1, • • • the following iteration:

u n+1 k+1 = G(T n , ∆T, u n k+1 ) + F (T n , ∆T, u n k ) -G(T n , ∆T, u n k )
In the parareal algorithm, the value u(T n ) is approximated by u n k at each iteration k with an accuracy that tends rapidly to the one achieved by the fine solver, when k increases. The coarse approximation G can be chosen much less expensive than the fine solver F by the use of a scheme with a much larger time step (even δT = ∆T ) δT δt (time step of the fine solver) or by using a reduced model. All the fine propagations are made in parallel over the time windows and the coarse propagations are computed in a sequential way but have a low computational cost. We refer to [START_REF] Maday | An adaptive parareal algorithm[END_REF] about the parallel efficiency of parareal and a recent work offering a new formulation of the algorithm to improve the parallel efficiency of the original one. The main convergence properties were studied in [START_REF] Gander | Analysis of the parareal time-parallel time-integration method[END_REF] and stability analysis was made in [START_REF] Staff | Stability of the parareal algorithm[END_REF], [START_REF]Parallelization in time of (stochastic) ordinary differential equations[END_REF].

Adaptation to multi-step time schemes

Case of a two-step fine time scheme

In the sequel, we will consider the case where the fine solver is a two-step time scheme. Hence we will use the following notation for the fine solver F that takes two initial values:

F (t, s, x, y), t ∈ [0, T ], s ∈ [0, T -t[
and the initial values x, y are in a Banach space U.

Example 2.1. If one solves (2.2) with a multi-step time scheme as fine propagator F like the second-order BDF method:

3 2 u j+1 -2u j + 1 2 u j-1 = -δtA(u j+1 , t j+1 ), j = 1, • • • , N f , t j+1 -t j = δt
Here the fine solver reads as: u j+1 = F (t j , δt, u j-1 , u j ). Now, we apply the parareal algorithm with a coarse grid: T 0 , • • • , T N c where:

T n+1 -T n = ∆T = N f δt.
Then we can write:

u(T n + jδt) u n,j , j = 1, • • • , N f , n = 1, • • • , N c .
In order to perform the fine propagation, in a given time window [T n , T n+1 ], we only need the local initial condition u n k and a consistent approximation of u(T nδt).

In [START_REF] Audouze | Symplectic multi-time step parareal algorithms applied to molecular dynamics[END_REF], the authors propose a consistent approximation in the context of the simulation of molecular dynamics. The proposed method was linked to the nature of the model and the symplectic character of their algorithm is shown, which is an important property to verify for molecular dynamics.

In the context of our application to the thermalhydraulic code Cathare, we want to derive a multistep variant of parareal that will not be intrusive in the software. We seek to derive a consistent approximation of u(T n -δt). The only fine trajectory at our disposal is

F (T n-1 , ∆T, u n-2,N f -1 k , u n-1 k ).
Its final value at time T n is:

F (T n-1 , ∆T, u n-2,N f -1 k , u n-1 k )(T n
) from which we compute u n k+1 by the parareal correction. Hence, we translate the solution:

F (T n-1 , ∆T -δt, u n-2,N f -1 k , u n-1 k
)(T nδt) by the same correction:

u n k+1 -F (T n-1 , ∆T, u n-2,N f -1 k , u n-1 k
) and obtain the so called consistent approximation u n-1,N f -1 k+1 to initialise the fine propagation in [T n , T n+1 ]. We now detail our algorithm:

                       u n+1 0 = G(T n , ∆T, u n 0 ), 0 ≤ n ≤ N -1 u n-1,N f -1 0 = u n 0 u n+1 k+1 = G(T n , ∆T, u n k+1 ) + F (T n , ∆T, u n-1,N f -1 k , u n k ) -G(T n , ∆T, u n k ), 0 ≤ n ≤ N -1, k ≥ 0 u n,N f -1 k+1 = F (T n , ∆T -δt, u n-1,N f -1 k , u n k ) + u n+1 k+1 -F (T n , ∆T, u n-1,N f -1 k , u n k ), 0 ≤ n ≤ N -1, k ≥ 0 (2.3)
Another option to treat this issue is to use a one-step time scheme to initialise the fine computation. We will see from the numerical results that this choice generates an error greater than the target accuracy and prevents the parareal algorithm to converge towards the solution with the desired accuracy.

The algorithm (2.3) adds consistency with the fine scheme. Also, this strategy can be generalised to multi-step time schemes involving several fine time steps preceding the time T n by applying the same correction to terms taking the form:

u n,N f -i k+1 , i = 1, • • • , I. Iteration k T 0 T 1 T 2 T 3 G(T 1 , ∆T, U 1 k+1 ) -G(T 1 , ∆T, U 1 k ) F (T 1 , ∆T, U k (T 1 -δt), U 1 k ) U 1 k U 2 k+1 U k+1 (T 2 -δt) Figure 2.1: Correction of u 1,N f -1 k+1 at time T 2 -δt in [T 1 , T 2 ] for the initialisation of the fine propagation in [T 2 , T 3 ]

Case of two-step coarse and fine time schemes

In the sequel, we will consider the case where the fine solver and the coarse solver are both two-step time schemes. Hence we will use the following notation for the coarse G solver that takes two initial values:

G(t, s, x, y), for t ∈ [0, T ], s ∈ [0, T -t[ and x, y ∈ U
We propose to add a correction for the solution at time T n -δT , where δT is the coarse time step, to initialise the coarse propagation in each time window in a consistent way. The parareal solution at time T n -δT and k-th parareal iteration is:

u n-1,N f -R k u((n -1)∆T + (N f -R)δt) = u(T n -δT )
where R is the ratio between the coarse and the fine time steps: R = δT δt .

The full multi-step parareal algorithm (2.4) makes two additionnal corrections compared to the classical parareal algorithm when the coarse and fine propagators are based on one two-step time schemes: one at times T nδt (see figure 2.1) and the other at times T n -δT (see figure 2.2).

                       u n+1 k+1 = G(T n , ∆T, u n-1,N f -R k , u n k+1 ) + F (T n , ∆T, u n-1,N f -1 k , u n k ), -G(T n , ∆T, u n-1,N f -R k , u n k ), 0 ≤ n ≤ N -1, k ≥ 0 u n,N f -1 k+1 = F (T n , ∆T -δt, u n-1,N f -1 k , u n k ) + u n+1 k+1 -F (T n , ∆T, u n-1,N f -1 k , u n k ), 0 ≤ n ≤ N -1, k ≥ 0 u n,N f -R k+1 = G(T n , ∆T -δT, u n-1,N f -R k+1 , u n k+1 ) + F (T n , ∆T -δT, u n-1,N f -1 k , u n k ) -G(T n , ∆T -δT, u n-1,N f -R k , u n k ), 0 ≤ n ≤ N -1, k ≥ 0 (2.4) Iteration k T 0 T 1 T 2 -δT T 2 G(T 1 , ∆T, U k+1 (T 1 -δT ), U 1 k+1 ) F (T 1 , ∆T, U k (T 1 -δt), U 1 k ) G(T 1 , ∆T, U k (T 1 -δT ), U 1 k ) U 1 k U 1 k+1 U 2 k+1 U k+1 (T 2 -δT ) Figure 2.2: Correction of u 1,N f -R k+1 at time T 2 -δT in [T 1 , T 2 ] for the initialisation of the coarse propagation in [T 2 , T 3 ]
The full multi-step version (2.4) of the algorithm does not take into account the adaptive time stepping of the Cathare time scheme. A consequence of the adaptive time steppping can be seen on the update of the solution u n,N f -R k+1 . This quantity is corrected at each iteration as follows:

u k+1 (T n+1 -δT n k+1 ) = G k+1 (T n , ∆T -δT n k+1 , u n-1,N f -R k+1 , u n k+1 ) + F (T n , ∆T -δT n k , u n-1,N f -1 k , u n k ) -G k (T n , ∆T -δT n k , u n-1,N f -R k , u n k ) (2.5)
We distinguish the coarse solver G at iterations k and k + 1 by the subscript G k and G k+1 because the solver G does not have the same sequence of time steps at iterations k and k + 1. Hence: δT n k+1 = δT n k and the correction (2.5) combines a quantity at time T n+1 -δT n k+1 with other quantities at time T n+1 -δT n k . This correction becomes incoherent. To correct this inconsistency, we propose to store the times T n -δT n 0 for each time window, at the parareal initialisation. Then, we impose to the coarse solver to pass by the point T n -δT n 0 in [T n , T n+1 ], for every parareal iterations.

In the next section, we apply the multi-step parareal algorithm to two-phase test cases: the oscillating manometer and the simulation of a breach in the primary circuit of a nuclear reactor. The convergence analysis of this new algorithm is the subject of the next chapter.

Application to the Cathare code 2.3.1 The oscillating manometer

Here we apply the multi-step parareal algorithm to the resolution of an oscillating manometer ( [START_REF] Ait-Ameur | Multi-step variant of the parareal algorithm[END_REF][START_REF]Time-parallel algorithm for two phase flows simulation[END_REF]). This test case is proposed in [START_REF] Hewitt | [END_REF] for system codes to test the ability of each numerical scheme to preserve system mass and to retain the gas-liquid interface. In [START_REF] Hewitt | [END_REF], the oscillating manometer is proposed as a numerical benchmark test for system codes to test the ability of each numerical sheme to preserve system mass and to retain the gas-liquid interface. It consists in a U-shaped tube manometer which is connected at the top, so that a closed system is formed. The system contains initially gas and liquid with the liquid forming equal levels in each arm of the manometer. Further, all parts of the fluid system have initially a uniform non zero velocity, but zero acceleration. Under these initial conditions, a hydrostatic pressure hypothesis is made In this test case, the phases are separated and the interfacial friction term will be important in this configuration.

Model

The model used in Cathare is the 6 equation two-fluid model that considers a set of balance laws (mass, momentum and energy) for each phase, liquid and vapor. It assumes independent velocities and a pressure equilibrium. The unknowns are the volume fraction α k ∈ [0, 1], the pressure p ≥ 0, the velocity u k and the enthalpy H k of each phase. The subscript k stands for l if it is the liquid phase and g for the gas phase. For the sake of simplicity, we write the terms of the model involved in our test case, studied in section 2.3.1.

           ∂ t (α k ρ k ) + ∂ x (α k ρ k u k ) = 0 α k ρ k ∂ t u k + α k ρ k u k ∂ x u k + α k ∂ x p = α k ρ k g + F int k ∂ t α k ρ k H k + u 2 k 2 + ∂ x α k ρ k u k H k + u 2 k 2 = α k ∂ t p + α k ρ k u k g (2.6)
with α v + α l = 1 and the two equations of state :

ρ k = ρ k (p, H k ).
The interfacial forces F int k are of 2 types. The first ensures hyperbolicity of the system (see [START_REF] Ndjinga | Influence of interfacial pressure on the hyperbolicity of the two-fluid model[END_REF] for the well-posedness of the 6 equation model). The second is the interfacial friction term that has an important role for our test case. For the oscillating manometer, the phases are separated which means that one of the two phases vanishes in some parts of the domain. It is numerically challenging to compute the velocity of the vanishing phase (see [START_REF] Ndjinga | A 2 × 2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF]). For this reason, the Cathare scheme forces the two velocities to be equal with the use of a damping term: the interfacial friction term.

Numerical method

The Cathare scheme is based on a finite volume method on a staggered grid (MAC scheme) and on a two step time scheme. In a staggered scheme the i-th component of the velocity is located at the center of the edge orthogonal to the i-th unit vector. Pressures, void fractions and enthalpies are cell-centered. Given a time discretisation T 0 , T 1 , T 2 , • • • of the full time interval [0, T ), we use the following notations: (α k ρ k ) n is an approximation of (α k ρ k ) at time T n . Here, we write the time discretisation of the Cathare scheme:

                         (α k ρ k ) n+1 -(α k ρ k ) n ∆t + ∂ x (α k ρ k u k ) n+1 = 0 (α k ρ k ) n+1 u n+1 k -u n k ∆t + (α k ρ k u k ) n+1 ∂ x u n+1 k + α n+1 k ∂ x p n+1 = (α k ρ k ) n+1 g + F n,n+1 k 1 ∆t (α k ρ k ) n+1 H k + u 2 k 2 n,n+1 -(α k ρ k ) n H k + u 2 k 2 n-1,n +∂ x α k ρ k u k H k + u 2 k 2 n+1 = α n+1 k p n+1 -p n ∆t + (α k ρ k u k ) n+1 g (2.7)
Where the notation F n,n+1 k shows that the time discretisation of F k is a function of the numerical solution at times T n and T n+1 . After discretisation, the non linear system is solved by a Newton method. In this test case, a numerical difficulty of two-phase flows simulations arises, namely the vanishing phase. An important issue is to guarantee the positivity of the volume fraction. Many schemes were designed to ensure this property (like [START_REF] Ndjinga | A 2 × 2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF] for two incompressible phases). The Cathare code uses a high interfacial friction to deal numerically with this difficulty.

About the convergence

In this section, we apply the multi-step parareal algorithm (2.4) to the simulation of the oscillating manometer. We use the same physical model and the same mesh (110 cells) for both the coarse and the fine solvers: the only difference is the size of the time steps, δt for F and ∆T for G. All the calculations have been evaluated with a stopping criteria where the tolerance is fixed to the fine solver accuracy, = 5 • 10 -2 . With this threshold, parareal convergence is achieved after 2 or 3 iterations. In the sequel, after giving a numerical proof of the convergence of the parareal algorithm in our test case, some results about measured speed-up will be presented. Figure 2.4 illustrates that the multi-step parareal algorithm effectively converges when applied to the problem of the oscillating manometer. For a given time step T n and parareal iteration k, the relative error in L 2 norm between the parareal solution and the sequential fine solver decreases beyond our given convergence threshold . In the figure, the test case has been solved with the multi-step parareal algorithm when δt = 10 -5 and ∆T = 10δt. These results are obtained on 16 time windows.

Speed-up performances

In the following strong scaling tests, the same setting is used for the multi-step parareal algorithm. The test case has been solved on an increasing number N proc of processes N proc = 5, 10, 15, • • • , 70. In figure 2.5, with 25 processes, we obtain a speed up of 3.4 and of 3.7 with 50 processes. Here, we observe two global trends: For N proc = {5, 10, 15, 20, 25, 40, 50}, the speed up first monotonically increases until reaching 25 processes and then increase again with 40 and 50 processes. This is due to the number of parareal iterations that is equal to 2 in this case For N proc = {30, 35, 45, 55, 60, 65, 70}, the speed up is drastically reduced because the parareal algorithm converges in 3 iterations in this case In the sequel, we recall the well-known dependence of the computational cost of the parareal algorithm on the number of iterations. Let T f ine be the CPU time to run the fine solver in a sequential way on the whole time interval [0, T ). Since the coarse time step is ten times greater than the fine time step we suppose that the cpu time of the coarse solver T coarse = T f ine 10 . This ratio between coarse and fine solvers should be as high as possible to minimise the computational cost of the coarse solver which is launched in a sequential way. When the algorithm converges in N it iterations, the coarse solver is launched N it times and the fine solver N it -1 times in parallel over the number of processes N proc . Hence, we can write the cpu time in parallel T para in terms of T f ine :

T para = (N it -1) T f ine N proc + N it T coarse + τ = N it -1 N proc + N it 10 T f ine + τ
where τ contains the time of communication between processes and the cpu time for the computation of the parareal corrections and of the error. Now, we can deduce an upper bound of the speed up S when the parareal algorithm converges in 2 or 3 iterations by neglecting τ :

S = T f ine T para ≤ 1 N it -1 Nproc + N it 10
Example: On 25 processes, the algorithm converges in 2 iterations: S ≈ 4 when the measured speed up is 3.4. On 35 processes, the algorithm converges in 3 iterations: S ≈ 2.8 when the measured speed up is 2.3.

An industrial test case

Description of the test case

In this section, we report our efforts to apply the parareal algorithm to the time parallelisation of an industrial test case representative of the numerical difficulties met in nuclear safety studies. This test case simulates a breach in the primary circuit downstrean on the reactor core of a Pressurised Water Reactor. The size of this breach is 512mm 2 . This accidental scenario is studied to simulate numerically the behaviour of the nuclear reactor under a hypothetic break of a pipe welding and during the emergency procedure following the accident. In our test case, we only simulates the reactor core before the emergency procedure. Here, the system is composed of: the primary circuit the downcomer where the fluid transits from the primary circuit to the reactor core the reactor core, modeled by two sets of uranium rods

The boundary conditions replicate the effects of the breach on the system:

Inlet boundary conditions: at t = 0 there is a single-phase liquid flow then the liquid flowrate decreases from t = 2.8s, time of the breach appearance, until reaching 0 at time t = 7.9s. After the breach, the volume fraction of the vapour phase increases until reaching the maximal value 1 at time t = 19s.

Outlet boundary conditions: the pressure within the system decreases since there is a leak of liquid.

Concerning the power generated by the uranium rods in the reactor core, it is set to zero at time t = 7.9s to model the effect of the control rods on reducing the reactivity.

After the breach and the disappearance of liquid in the reactor core, the fuel rods cladding is not anymore in contact with the coolant. Hence, the heat exchanges are limited and this can drive to the meltdown of the rod cladding.

Application of the parareal algorithm

The Parareal library is an intermediate between the solutions computed by Cathare and the corrections made by the parareal algorithm. The library collects the fine and coarse approximations from Cathare, extracts from the Cathare arrays the principal unknows (liquid and vapour velocities, pressure, enthalpies, void fraction) and stores them in other arrays belonging to the Parareal library. From there, the local initial conditions are updated with the parareal correction by combining the solutions coming from the coarse and fine propagations. Then, the updated initial conditions are copied in Cathare arrays. Before making new coarse and fine propagations the Cathare code needs auxiliary values in addition to the principal unknowns. These auxiliary values are of different natures:

variables computed from the principal unknowns variables tracked during the simulation: for example, the water level in an element of the system boolean variables giving the flow regime in an area of the system: bubbly, annular, dispersed, seperated phases, etc.

For the first category of auxiliary variables, they can be recalculated using the updated initial conditions. However, the others are the last saved quantities before the copy of the corrected initial conditions. Hence, there is an inconsistency between the updated variables depending on the principal unknowns and the other variables that are unchanged. We show the consequences of this inconsistency by the following numerical experiment.

We denote F ref the target fine solution we seek to approximate with the parareal algorithm. This solution is computed in a sequential way with a finer time step than the one of the fine solver used within the parareal algorithm, denoted F para . The propagator F ref replaces the exact propagator since we do not have the expression of the exact solution for this test case and will be called the reference solution in this section. On the one hand, we initialise in the Parareal library every time window with the principal variables of the reference solution F ref at times T n , n = 0, • • • , N . Then we transfer these initial conditions to the Cathare code to make N parallel fine propagations

F ref |[T n ,T n+1 ] (T n , ∆T, F ref (T 0 , T n - T 0 , u 0 )(T n )), where F ref |[T n ,T n+1 ]
is the restriction of the reference solution F ref to the time interval [T n , T n+1 ]. From this information, the Cathare code builds the Cathare array composed of principal unknowns and auxiliary variables and propagates this initial state over the time window [T n , T n+1 ]. We plot in figure 2.6 the following quantity for every time windows [T n , T n+1 ]:

||F ref |[T n ,T n+1 ] (T n , ∆T, F ref (T 0 , T n -T 0 , u 0 )(T n )) -F ref (T 0 , T n -T 0 , u 0 )|| L 2 ||F ref (T 0 , T n -T 0 , u 0 )|| L 2 (T i ), i = 0, • • • , N f (2.8)
where N f = ∆T δt is the number of fine time steps in a time window of size ∆T . This numerical experiment will allow us to see the error made on the recalculation of the auxiliary variables. If there was no error at this stage, the quantity (2.8) would be equal to zero. In figure 2.6, we see that the inconsistency between principal variables and auxiliary variables in the Cathare array leads to an additionnal error. This error prevents the parareal algorithm to recover the fine sequential solution with the target accuracy. This target accuracy is computed in the following way:

max n=1,••• ,N ||F para (T 0 , T n -T 0 , u 0 ) -F ref (T 0 , T n -T 0 , u 0 )|| L 2 (Ω) max n=1,••• ,N ||F ref (T 0 , T n -T 0 , u 0 )|| L 2 (Ω)
where F para is the fine solver used within the parareal algorithm. This target accuracy is thus the accuracy F of the fine solver compared to the reference solution F ref . On the other hand, the Parareal library transfers to the Cathare code the whole initial state with principal and auxiliary variables computed by the reference solution F ref (T 0 , T n -T 0 , u 0 ). Then from these initial states, the Cathare code makes N parallel fine propagations

F ref |[T n ,T n+1 ] (T n , ∆T, F ref (T 0 , T n -T 0 , u 0 )(T n )) spread over the N time windows [T n , T n+1 ].
In figure 2.7, we plot the error (2.8) for this new configuration. In figure 2.7, we see that sending all the state to initialise the fine propagations in [T n , T n+1 ] improves the accuracy of the fine solutions computed by the Cathare code, compared to figure 2.6 where we only initialise the time propagations with the principal variables. However, we still observe a non negligible error of about 10 -4 , especially starting from time t = 5, after the appearance of the breach in the system that occurs at time t = 2.8. We propose a strategy to improve the transfer of the auxiliary variables between time windows. The method consists in correcting at each parareal iteration the vector of auxiliary values in a specific way. Firstly, we need to define V aux , the vector of auxiliary variables and V p the vector of principal variables. We distinguish:

V aux (V p (T n )) the auxiliary variables computed from the principal variables at time T n .

V aux (T n ) the auxiliary variables computed by the Cathare code at time T n .

Secondly, we denote T n

-the point belonging to the interval [T n-1 , T n ] and T n + the one belonging to [T n , T n+1 ]. We propose the following correction for the auxiliary variables at time T n and (k+1)-th parareal iteration:

V aux k+1 (T n + ) = V aux (V p k+1 (T n + )) -V aux (V p k (T n -)) + V aux k (T n -).
(2.9)

Since the Cathare time scheme is a two-step time scheme, we also apply a correction to the auxiliary variables at time T nδt.

V aux k+1 (T n -δt) = V aux (V p k+1 (T n -δt)) -V aux (V p k (T n -δt)) + V aux k (T n -δt). ( 2

.10)

This choice for the correction of the auxiliary values is motivated by the two numerical experiments 2.6 and 2.7. The figure 2.6 illustrates the error we make by computing the auxiliary variables with:

V aux (T n + ) = V aux (V p (T n + )).
In this configuration, the error is greater than the accuracy of the fine solver and is the main actual barrier to efficiently apply the parareal algorithm to the Cathare code in a non intrusive way. On the other hand, the figure 2.7 illustrates the error we make by computing the auxiliary variables with:

V aux (T n + ) = V aux (T n -).
In this case, the error is lower than the accuracy of the fine solver but is still non negligible. The correction we propose (2.9) combines the two approaches with a third term V aux (V p (T n )). It is important to mention that when the parareal algorithm converges, we obtain:

V aux (T n + ) = V aux (T n -),
since we have:

V p (T n + ) = V p (T n -)
, and then :

V aux (V p (T n + )) = V aux (V p (T n -))
We have not implemented this promising strategy yet to handle specifically the auxiliary variables of the Cathare code. The actual option to reconstruct the auxiliary variables is to compute them from the principal variables. In figures 2.6 and 2.7, we see that the reconstruction error is greater than the fine solver accuracy on 10 time windows: in figure 2.6, the reconstruction error is about 4 × 10 -2 while the fine solver accuracy is F = 5 × 10 -3 . Hence we can not reach the fine solver accuracy with the parareal algorithm since the reconstruction error of the auxiliary variables will dominate and pollute the simulation. However, the reconstruction error of the auxiliary variables is in the order of the fine solver accuracy F on 5 time windows. In figures 2.8 and 2.9, we make the same numerical experiment as in figures 2.6 and 2.7 on 5 time windows. We see in figure 2.8 that, in this particular case, the reconstruction of the auxiliary variables from the principal variables generates an error in the order of F at the initial times of the 5 time windows. In figure 2.9, we observe the same behaviour as on 10 processors: sending all the state to initialise the fine propagations in [T n , T n+1 ] improves the accuracy of the fine solutions computed by the Cathare code.

Hence, in the case of 5 processors, the reconstruction error of the auxiliary variables from the principal variables is in the order of the accuracy of the fine solver and may not pollute the parareal algorithm. In the sequel, we investigate the performances of the parareal algorithm when the reconstruction error is close to the fine solver accuracy to illustrate the behaviour of the parareal method when we will implement the specific treatment of the auxiliary values (2.9-2.10). Hence, we apply the parareal algorithm to the simulation of a breach in the primary circuit on 5 processors. The reference solution F ref is computed with a fine time step δt ref = 10 -4 . The fine solver F para within the parareal algorithm has a time step δt = 10 -4 and an accuracy F = 5 × 10 -3 , hence the target accuracy of the parareal algorithm is fixed to 5 × 10 -3 . The coarse solver has an accuracy G = 10 -1 with a coarse time step δT = 0.5 = 500δt. We apply the full multi-step parareal (2.4) that corrects both the coarse and the fine solvers to properly initialise the time propagations. The simulation lasts 25 seconds and we split the time interval over 5 processors. In figure 2.10, we plot the following error at each time T n = n∆T and each parareal iteration k, where n = 0, • • • , 5 and ∆T is the size of the time window :

||F para (T n , ∆T, u n k ) -F ref (T 0 , T n -T 0 , u 0 )|| L 2 (Ω) ||F ref (T 0 , T n -T 0 , u 0 )|| L 2 (Ω) (T n ), n = 0, • • • , 5
In this case, the parareal algorithm reaches the fine solver accuracy at the initial times of the 5 time windows after 3 iterations. Hence, when the reconstruction error is close to the fine solver accuracy the parareal algorithm can reach the target accuracy for the solution at times T n . In figures 2.11, 2.12 and 2.13, we see the volume fraction α v in one cell of the mesh computed by the fine solver in the multi-step parareal algorithm for the three first iterations, compared to the reference solution. After the first parareal iteration, we clearly distinguish the fine propagations made by the parareal algorithm and the reference solution, over the different time windows. Then for the second and the third iteration, the reference and the fine parareal solutions are very close. We capture the dynamic of the test case with α v 0 at time t = 0 since the system is initially filled with water. Then the vapour volume fraction increases from t = 2.8s, time of the breach appearance until reaching its maximal value 1 at time t = 19s.

Conclusion

In this chapter, we developed two strategies to apply the parareal algorithm to the Cathare code: a numerical clone of Cathare that is restricted to one test case and a Parareal library that uses the Cathare code in a non intrusive way. The main contribution of this work has been to adapt the parareal algorithm to the architecture of the software and to its time discretisation in a non intrusive way. The results obtained with the numerical clone on the oscillating manometer show that the parareal algorithm can effectively speed-up two-phase flows simulations. These preliminary results illustrate the behaviour of the multi-step parareal algorithm on a test case that is representative of the numerical challenges for two phase flows. However, the Parareal library that uses the Cathare The reference solution and the fine solution after 3 parareal iterations code as a black box suffers from a lack of accuracy. This is due to the data structure of the Cathare code that depends on principal and auxiliary variables. Hence, a reconstruction of the auxiliary values is necessary after the parareal update of the principal variables. The actual reconstruction error is greater than the target accuracy and prevents the parareal algorithm to converge towards the desired solution. We derived a strategy to accurately reconstruct the auxiliary variables by adding a parareal update specific to them and the efficiency of this method will be the subject of future works.

Chapter 3

Convergence analysis of the multi-step variant of the parareal algorithm In this paper, we consider the problem of accelerating the numerical simulation of time dependent problems involving a multi-step time scheme by the parareal algorithm. A multi-step time scheme can potentially bring higher approximation orders than plain one-step methods but the initialisation of each time window needs to be appropriately chosen. Our main contribution is the design and analysis of an algorithm adapted to this type of discretisation without being intrusive in the coarse or fine propagators. This property allows to apply this variant of the parareal algorithm on a software as a black box and ensures the portability of the method. The parareal method is based on combining predictions made by two propagators: an accurate and expensive one used in a parallel way over the time windows and a coarse and cheaper one used in a sequential way. At convergence, the parareal algorithm provides a solution that has the fine solver's accuracy. In the classical version of parareal, the local initial condition of each time window is corrected at every iteration. When the fine and/or coarse propagators is a multi-step time scheme, we need to choose a consistent approximation of the solutions involved in the initialisation of the fine solver at each time window. Otherwise, we could loose one of the well known property of the parareal method: to recover the fine solution at the machine precision after N iterations, where N is the number of time windows. In this paper, we develop a variant of the algorithm that overcomes this obstacle. Thanks to this, the parareal algorithm is more coherent with the underlying time scheme and we recover the properties of the original version. We show both theoretically and numerically that the accuracy and convergence of the multi-step variant of parareal algorithm are very competitive when we carefully choose the initialisation of each time window.

Introduction

Solving complex models with high accuracy and within a reasonable computing time has motivated the search for numerical schemes that exploit efficiently parallel computing architectures. In this paper, the model consists of a Partial Differential Equation (PDE) set on a domain D. In this context, one of the main ideas to parallelize a simulation is to break the problem into subproblems defined over subdomains of a partition of D. The domain can potentially have high dimensionality and be composed of different variables like space, time, velocity or even more specific variables for some problems. There exist algorithms with very good scalability properties for the decomposition of the spatial variable (see [START_REF] Quarteroni | Domain decomposition methods for partial differential equations[END_REF] or [START_REF] Toselli | Domain decomposition methods: algorithms and theory[END_REF] for an overview) and time domain decomposition is more and more considered to complement this strategy when the speed up performances stagnates despite remaining computing resources. Research on time parallel algorithms is currently very active and has by now a history of at least 50 years (back to at least [START_REF] Nievergelt | Parallel methods for integrating ordinary differential equations[END_REF]) during which several algorithms have been explored (see [START_REF] Gander | 50 years of time parallel time integration[END_REF] for an overview).

In this work, we report our recent effort to adapt one particular time-parallel algorithm: the parareal in time algorithm, to multi-step time schemes. The method was first introduced in [START_REF] Lions | Résolution par un schéma en temps "pararéel[END_REF] and has been well accepted by the community because it is easily applicable to a relatively large spectrum of problems (some specific difficulties are nevertheless encountered on certain types of PDEs as reported in [START_REF] Dai | Stable Parareal in time method for first-and second-order hyperbolic systems[END_REF][START_REF] Farhat | Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid-structure applications[END_REF] for hyperbolic systems or [START_REF] Dai | Symmetric parareal algorithms for hamiltonian systems[END_REF] for hamiltonian problems). Some limitations persist for the classical version of the parareal algorithm like the parallel efficiency that decreases with the final number of iterations K as 1/K. This limitation is addressed in [START_REF] Maday | An adaptive parareal algorithm[END_REF] that proposes an adaptive variant of the parareal method where the only remaining factor limiting high performance becomes the cost of the coarse solver. Without entering into very specific details of the algorithm at this stage, we can summarize the procedure by saying that we build iteratively a sequence to approximate the exact solution of the problem by a predictor-corrector algorithm. At every iteration, predictions are made by a solver which has to be as numerically inexpensive as possible since it is run on the full time interval. It usually involves coarse physics and/or coarse resolution. Corrections involve an expensive solver with high-fidelity physics and high resolution which is propagated in parallel over small time subdomains. In the classical version of parareal, the fine solver has a fixed high accuracy across all iterations. It is set to the one that we would use to solve the dynamics at the desired accuracy with a purely sequential solver. At each iteration, the local initial conditions are corrected for every time windows until convergence. Multi-step time schemes require several previous steps to compute the solution at a new point in time. When the fine and/or coarse propagators is a multi-step time scheme, we need to choose a consistent approximation of the solutions at previous steps involved in the initialisation of the fine and/or coarse solver at each time windows. Otherwise, the initialisation error will prevent the parareal algorithm to converge towards the solution with fine solver's accuracy. This point was addressed in the context of multigrid in time method in [START_REF] Falgout | A parallel-in-time algorithm for variable step multistep methods[END_REF][START_REF] Falgout | Multigrid methods with space-time concurrency[END_REF]. Here, the authors adapt the MGRIT algorithm framework to the use of multi-step time schemes, the BDF methods. In this paper, we propose a variant of the algorithm that overcome this obstacle. Thanks to this, the parareal algorithm is more coherent with the underlying time scheme and we recover the properties of the original version.

We present in section 2 the variant of the parareal algorithm adapted to multi-step time schemes. This method includes additional corrections at previous steps involved in the intialisation of the fine and/or coarse solver at each time window. This choice has the benefit to be non intrusive into the code we seek to parallelise by a time domain decomposition.

In the last section, we illustrate the performance of the algorithm on numerical examples: the damped oscillator and the Brusselator. We show that this variant allows the parareal algorithm to converge towards the solution with fine solver's accuracy.

A multi-step variant of the parareal algorithm

In this section, after introducing some preliminary notations in section 2.1, we formulate the new variant of the parareal algorithm adapted to multi-step time schemes (section 2.2). We then present the hypothesis we consider in this article and restrict ourselves to two-step time schemes for the convergence analysis (section 2.3). We prove that the multi-step variant converges with a convergence rate similar to that of the classical parareal algorithm. Finally, we discuss how the new paradigm can be generalised to multi-step time schemes, not only two-step times schemes, used in the fine and/or the coarse solver (section 2.4).

Setting and preliminary notations

Let U be a Banach space of functions defined over a domain

Ω ⊂ R d (d ≥ 1). Let S : [0, T ] × [0, T ] × U → U (3.1)
be a solver, that is, an operator such that, for any given time t ∈ [0, T ], s ∈ [0, Tt] and any function w ∈ U takes an initial value at time t and propagates it at time t + s. We further assume that S is defined through the discretisation of the time-dependent problem:

u (t) + A(t, u(t)) = 0, t ∈ [0, T ] u(0) ∈ U (3.2)
where A is a locally Lipschitz operator from [0, T ]×U. This ODE system can also be obtained from the discretisation of a PDE in space with A : R × R N → R N and N denotes the number of degrees of freedom. We seek to approximate the solution of problem (3.2) at a given target accuracy by a solver S. We denote ε(t, s, w) the propagator giving the exact solution of system (3.2), for any initial value w ∈ U, any t ∈ [0, T ] and any s ∈ [0, Tt]. Thus, S(t, s, w) approximates ε(t, s, w) with an accuracy η > 0 such that we have:

||δS(t, s, w)|| = ||ε(t, s, w) -S(t, s, w)|| ≤ ηs(1 + ||w||), ∀t ∈ [0, T ], s ∈ [0, T -t], w ∈ U (3.3)
where || • || denotes the norm in U.

The choice of the solver S determines the quality of the approximation and the computational cost of its implementation. One can potentially bring higher approximation orders than plain one-step methods by using a multi-step time discretisation method. Multi-step time schemes require several previous steps to compute the solution at a new point in time. Hence, the notation (3.1) does not hold when the solver is based on a multi-step time scheme. Here and in the following, we consider only two-step time schemes and we will use the following notation:

S : [0, T ] × [0, T ] × U × U → U (3.4)
such that S(t, s, w 1 , w 2 ) for any given time t ∈ [0, T ], s ∈ [0, Tt] and any functions w 1 , w 2 ∈ U takes two initial values at times t and tδt and propagates them at time t + s, where δt is the time step of the two-step time scheme underlying in S.

A multi-step variant of the parareal algorithm

We consider a given decomposition of the time interval

[0, T ] into N subintervals [T n , T n+1 ], n = 0, • • • , N -1.
Without loss of generality, we will take them of uniform size ∆T = T /N which means that T n = n∆T for n = 0, • • • , N . For a given target accuracy η > 0, the goal of the parareal algorithm is to accelerate the computation of an approximation ũ(T n ) of u(T n ) such that:

max 1≤n≤N ||u(T n ) -ũ(T n )|| ≤ η
The classical way to compute such an approximation is to set ũ(T n ) = S seq (0, T n , u(0)), 1 ≤ n ≤ N , where S seq is some sequential solver in [0, T ]. On the other hand, the strategy of the parareal algorithm follows the following steps, using two propagation operators:

G(T n , ∆T, u n ) computes a coarse approximation of u(T n+1 ) with initial condition u(T n ) u n .
The coarse propagation is sequential but have a low computational cost.

F (T n , ∆T, u n ). computes a more accurate approximation of u(T n+1 ) with initial condition u(T n ) u n . The action of F is distributed over N time windows and N processors solve over each interval [T n , T n+1 ] of size ∆T instead of solving over [0, T ].

In the sequel, we analyse the convergence rate of the multi-step variant of parareal algorithm when the coarse solver is a one-step time scheme and the fine one is a two-step time scheme.

Hypotheses (H):

There exists G , C d , C > 0 such that for any functions x, y ∈ U and for any t ∈ [0, T ] and s ∈ [0, Tt],

||ε(t, s, x) -G(t, s, x)|| ≤ s(1 + ||x||) G ⇔ ||δG(t, s, x)|| ≤ s G (1 + ||x||) (3.5) ||G(t, s, x) -G(t, s, y)|| ≤ (1 + Cs)||x -y|| (3.6) ||F (t, s, x 1 , y 1 ) -F (t, s, x 2 , y 2 )|| ≤ (1 + Cs) (||x 1 -x 2 || + ||y 1 -y 2 ||) (3.7) ||δG(t, s, x) -δG(t, s, y)|| ≤ C d s G ||x -y|| (3.8) || (F (t, s, ε(t, -δt, y 1 ), y 1 ) -ε(t, s, y 1 )) -(F (t, s, ε(t, -δt, y 2 ), y 2 ) -ε(t, s, y 2 )) || ≤ Csδt||y 1 -y 2 || (3.9) ||(F (t, s -δt, y 1 -δ 1 , y 1 ) -F (t, s -δt, y 2 -δ 2 , y 2 )) -(F (t, s, y 1 -δ 1 , y 1 ) -F (t, s, y 2 -δ 2 , y 2 ))|| ≤ (1 + Cs)δt||δ 1 -δ 2 || + (1 + Cs)δt||y 1 -y 2 || (3.10) ||F (t, s, ε(t, -δt, y), y) -ε(t, s, y)|| ≤ s F (1 + ||y||) (3.11)
Note that the hypothesis (3.5)-(3.8) are the classical properties of numerical schemes related to stability and accuracy. Hypotheses (3.6) and (3.7) are Lipschitz conditions and the quantity G is a small constant which, in the case of the explicit Euler scheme, would be proportionnal to the time step size. Hypotheses (3.9) and (3.10) are specific to two-step time schemes. There are two sources of error for two-step time schemes:

the error from the discretisation of the time derivative, common to one-step time schemes.

the error from the inconsistency between the two initial conditions x 1 and y 1 in F (t, s, x 1 , y 1 ).

In hypothesis (3.9), assuming the term x 1 = ε(t, -δt, y 1 ), y 1 ) is defined, there is no inconsistency between the two initial values since x 1 is computed with the exact propagator starting from y 1 . Hence, the only remaining errors are:

the difference between y 1 and y 2 .

the error from the time propagation over a time window of size s.

the error between the fine and the exact propagators that is proportionnal to the fine time step δt.

On the other hand, we assume hypothesis (3.10) holds for s ≥ ∆T , the time window size. This hypothesis includes the inconsistency between the two initial values and is denoted δ 1 and δ 2 . Hence, we describe here the errors coming from:

the inconsistency δ i between the two initial values of the fine solver.

the difference between the principal initial values y 1 and y 2 .

the time propagation over time windows of size s.

Example 3.1. Here, we illustrate the validity of hypothesis (3.10) on a simple linear ODE. The parameters involved in hypothesis (3.10) are: s, y 1 , y 2 , δ 1 and δ 2 . In the proof of convergence, we apply this hypothesis for y 1 = u n k-1 , the parareal solution at k -1 iteration and time T n , and y 2 = u(T n ), the exact solution at time T n , hence these two parameters are very close. On the other hand,

δ 1 = u n k-1 -u n-1,N f -1 k-1
, where

u n-1,N f -1 k-1
is the parareal solution at (k -1) iteration and time T nδt, and δ 2 = u(T n )u(T nδt), where u(T nδt) is the exact solution at time T nδt and s is equal to the time window size. y (t) = y(t), t ∈ [0, T ] y 0a , y 0b given where y 0a , y 0b are the two seed values to initialise the time propagation with the second-order BDF method:

y 0a = (y 1 -δ 1 ) -(y 2 -δ 2 ), y 0b = y 1 -y 2
We solve this system by a second-order BDF method:

3 2 y n+1 -2y n + 1 2 y n-1 = δty n+1 (3.12)
From (3.12), we have the expression of the numerical solution

y n for n = 0, • • • , N f with N f = T δt : y n = αr n 1 + βr n 2
such that:

r 1 = 2 + √ 1 + 2δt 3 -2δt = 1 + δt + O(δt 2 ) r 2 = 2 - √ 1 + 2δt 3 -2δt = 1 3 - δt 9 + O(δt 2 ).
In (3.12), the term r n 2 tends rapidly to zero when n goes to infinity. Thus we neglect its contribution.

α = r 2 (δ 1 -δ 2 ) + (1 -r 2 )(y 1 -y 2 ) r 1 -r 2 β = (r 1 -1)(y 1 -y 2 ) -r 1 (δ 1 -δ 2 ) r 1 -r 2
In the linear case, we can write hypothesis (3.10):

||y N +1 -y N || ≤ (1 + Cs)δt||y 0a -y 0b || + (1 + Cs)δt||y 0b ||
where: y N = F (0, sδt, y 0a , y 0b ), y N +1 = F (0, s, y 0a , y 0b ) and (N + 1) is the number of fine time steps in a time window of size s:

N + 1 = s δt = ∆T δt .
From the expression of y n in (3.12):

y N +1 -y N = αr N 1 (r 1 -1) + βr N 2 (r 2 -1)
Neglecting the term r N 2 , we obtain:

y N +1 -y N = r 2 r 1 -r 2 (1 + ∆T )δt(δ 1 -δ 2 ) + 1 -r 2 r 1 -r 2 (1 + ∆T )δt(y 1 -y 2 ) + O(δt 2 )
Hence, the second-order BDF method verify hypothesis (3.10).

In the following example, we explain the problem of initialising the fine solver in a time window when the fine propagator is a two-step time scheme:

Example 3.2. If one solves (3.
2) with a multi-step time scheme as fine propagator F like the second-order BDF method:

3 2 u j+1 -2u j + 1 2 u j-1 = -δtA(u j+1 , t j+1 ), j = 1, • • • , N f , t j+1 -t j = δt
Here the fine solver reads as: u j+1 = F (t j , δt, u j-1 , u j ). Now, we apply the parareal algorithm with a coarse grid: T 0 , • • • , T N where:

T n+1 -T n = ∆T = N f δt.
Then we can write:

u(T n + jδt) u n,j , j = 1, • • • , N f , n = 1, • • • , N c .
In order to perform the fine propagation, in a given time window [T n , T n+1 ], we only need the local initial condition u n k and a consistent approximation of u(T nδt).

In [START_REF] Audouze | Symplectic multi-time step parareal algorithms applied to molecular dynamics[END_REF], the authors propose a consistent approximation in the context of the simulation of molecular dynamics. The proposed method was linked to the nature of the model and the symplectic character of their algorithm is shown, which is an important property to verify for molecular dynamics. We now detail our algorithm:

                       u n+1 0 = G(T n , ∆T, u n 0 ), 0 ≤ n ≤ N -1 u n,N f -1 0 = u n+1 0 , 0 ≤ n ≤ N -1 u n+1 k+1 = G(T n , ∆T, u n k+1 ) + F (T n , ∆T, u n-1,N f -1 k , u n k ) -G(T n , ∆T, u n k ), 0 ≤ n ≤ N -1, k ≥ 0 u n,N f -1 k+1 = F (T n , ∆T -δt, u n-1,N f -1 k , u n k ) + u n+1 k+1 -F (T n , ∆T, u n-1,N f -1 k , u n k ), 0 ≤ n ≤ N -1, k ≥ 0 (3.13)
At this point, several comments are in order. To derive a consistent approximation of u(T nδt), we use the only fine trajectory at our disposal which is

F (T n-1 , ∆T, u n-2,N f -1 k , u n-1 k
). Its final value at T n is:

F (T n-1 , ∆T, u n-2,N f -1 k , u n-1 k )(T n
) from which we compute u n k+1 by the parareal correction. Hence, we translate the solution:

F (T n-1 , ∆T -δt, u n-2,N f -1 k , u n-1 k
)(T nδt) by the same correction:

u n k+1 -F (T n-1 , ∆T, u n-2,N f -1 k , u n-1 k
) and obtain the so called consistent approximation u n-1,N f -1 k+1 to initialize the fine propagation in [T n , T n+1 ].

Moreover, an important feature of this new algorithm is to preserve a well known property of the parareal algorithm:

u n k = F (T 0 , T n -T 0 , u 0 ), for: k ≥ n, n = 0, • • • , N (3.14)
This comes from the term:

G(T n , ∆T, u n k+1 ) -G(T n , ∆T, u n k ) that is equal to zero when k ≥ n, n = 0, • • • , N .
In our case, the multi-step variant of the parareal algorithm verifies (3.14) and the additionnal correction of the solution at time T nδt leads to:

u n,N f -1 k+1 = F (T n , ∆T -δt, u n-1,N f -1 k , u n k ) + u n+1 k+1 -F (T n , ∆T, u n-1,N f -1 k , u n k ), 0 ≤ n ≤ N -1, k ≥ 0 = F (T n , ∆T -δt, u n-1,N f -1 k , u n k ) + G(T n , ∆T, u n k+1 ) -G(T n , ∆T, u n k ), 0 ≤ n ≤ N -1, k ≥ 0
Hence, the multi-step parareal method satisfies the same property (3.14) at time T nδt:

u n,N f -1 k = F (T 0 , T n -δt -T 0 , u 0 ), for: k ≥ n, n = 0, • • • , N (3.15)
The convergence result of theorem 3.3 and its proof are helpful to understand the main mechanisms driving the convergence of the algorithm and explaining its behavior. To present it, we introduce the shorthand notation for the error norm:

E n k := u n k -ε(T 0 , T n -T 0 , u 0 ), k ≥ 0, 0 ≤ n ≤ N,
We introduce the following quantities:

             α := C d G ∆T µ := C∆T δt β := 1 + C c ∆T γ G := ∆T G max 0≤n≤N (1 + ||u(T n )||) γ F := ∆T F max 0≤n≤N (1 + ||u(T n )||) (3.16)
as shorthand notations for the proof of convergence. Theorem 3.3 (Convergence of the multi-step parareal algorithm). Let G, F and δG satisfy Hypotheses (3.5)- (3.11). Let k ≥ 0 be any given positive integer. If the time step δt of the fine solver verifies:

δt ≤ ∆T 2 2 G , (3.17)
then the (u n k ) n of the multi-step parareal scheme (3.13) satisfy:

       max 0≤n≤N ||u n 0 -u(T n )|| ≤ γG γ G e CT max n (1 + ||u(T n )||)T e -C∆T G , n ≥ 1 max 0≤n≤N ||u n k -u(T n )|| ≤ λ τ k+1 k + 1! f k 2 k+1 γG γ G α α k + k + 1 τ f k-1 2 k+1 γF γ G α α k-1 , n ≥ k + 1, k ≥ 1 (3.18) where: λ = e CT max 0≤n≤N (1 + ||u(T n )||) C d , τ = 2τ = 2C d T e -C∆T G , f k = (1 + √ 5) k+1 -(1 - √ 5) k+1 2 k+1 √ 5 .
and α, γG and γF are perturbations of the coefficients α, γ G and γ F respectively, such that:

γG γ G = 1 + O(∆T G ), α α = 1 + O(∆T G ), γF γ G = O(∆T G )
Let us make a couple of remarks before giving the proof of the theorem. First, the convergence rate of the multi-step parareal algorithm is similar to the one of the classical parareal algorithm with the factor τ k+1 k + 1! , since in the classical version the convergence rate is

τ k+1 k + 1! . The remaining factors γG γ G , α α k and γF γ G
are close to 1 and their contributions are negligible. The term f k is specific to the multi-step variant and have the following asymptotic behaviour when k tends to infinity:

f k 2 k+1 ∼ k→+∞ 1 + √ 5 4 k+1
Proof. The proof is in the spirit of existing results from the litterature [START_REF] Lions | Résolution par un schéma en temps "pararéel[END_REF][START_REF]Parallelization in time of (stochastic) ordinary differential equations[END_REF][START_REF] Gander | Nonlinear convergence analysis for the parareal algorithm[END_REF]. If k = 0, using definition (3.13) for u n 0 , we have for 0 ≤ n ≤ N -1,

E n+1 0 = u n+1 0 -ε(T 0 , T n+1 -T 0 , u 0 ) E n+1 0 = G(T n , ∆T, u n 0 ) -ε(T n , ∆T, u(T n )) ||E n+1 0 || ≤ ||G(T n , ∆T, u n 0 ) -G(T n , ∆T, u(T n ))|| + ||G(T n , ∆T, u(T n )) -ε(T n , ∆T, u(T n ))|| ≤ (1 + C∆T )||E n 0 || + ∆T G (1 + ||u(T n )||) ≤ β||E n 0 || + γ G
where we have used (3.5) and (3.6) to derive the second to last inequality.

For k ≥ 1, starting from (3.13), we have

E n+1 k = u n+1 k -ε(T 0 , T n+1 -T 0 , u 0 ) = G(T n , ∆T, u n k ) + F (T n , ∆T, u n-1,N f -1 k-1 , u n k-1 ) -G(T n , ∆T, u n k-1 ) -ε(T n , ∆T, u(T n ))
In the sequel, we add and substract the following quantites to E n+1 k :

G(T n , ∆T, u(T n )) and ε(T n , ∆T, u n k-1 ) ε(T n , ∆T, u(T n )) and F (T n , ∆T, u(T n ) -δx, u(T n )) F (T n , ∆T, u n k-1 -δx, u n k-1 ) and F (T n , ∆T, u n k-1 -δx, u n k-1 ) where: δx = u(T n ) -u(T n -δt), δx = u n k-1 -u n-1,N f -1 k-1 , δx = u n k-1 -ε(T n , -δt, u n k-1 ) E n+1 k = G(T n , ∆T, u n k ) -G(T n , ∆T, u(T n )) + δG(T n , ∆T, u(T n )) -δG(T n , ∆T, u n k-1 ) + F (T n , ∆T, ε(T n , -δt, u n k-1 ), u n k-1 ) -ε(T n , ∆T, u n k-1 ) -(F (T n , ∆T, ε(T n , -δt, u(T n )), u(T n )) -ε(T n , ∆T, u(T n ))) + F (T n , ∆T, u n k-1 -δx, u n k-1 ) -F (T n , ∆T, u n k-1 -δx, u n k-1 ) + F (T n , ∆T, u n k-1 -δx, u n k-1 ) -F (T n , ∆T, u n k-1 -δx, u n k-1 ) + F (T n , ∆T, u(T n -δt), u(T n )) -ε(T n , ∆T, u(T n ))
Taking norms and using (3.6), (3.7), (3.8), (3.9), (3.11), we derive:

||E n+1 k || ≤ (1 + C∆T )||E n k || + C∆T G ||E n k-1 || + C∆T δt||E n k-1 || + C|| δx -δx|| + C||δx -δx|| + ∆T F (1 + ||u(T n )||)
On the one hand, the term δx -δx becomes:

δx -δx = u n-1,N f -1 k-1 -u(T n -δt) -(u n k-1 -u(T n )) = E n,N f -1 k -E n+1 k = δE n+1 k
On the other hand, we derive a bound for the term: ||δx -δx||:

||δx -δx|| = ||u(T n -δt) -ε(T n , -δt, u n k-1 ) -(u(T n ) -u n k-1 )||
Writing the Taylor expansions of u(T nδt) and ε(T n , -δt, u n k-1 ) around T n and u n k-1 respectively, we obtain formally:

u(T n -δt) -u(T n ) = δtA(T n , u(T n )) + δt 2 2 ∂A ∂t + ∂A ∂u A (T n , u(T n )) + O(δt 3 ) ε(T n , -δt, u n k-1 ) -u n k-1 = δtA(T n , u n k-1 ) + δt 2 2 ∂A ∂t + ∂A ∂u A T n , u n k-1 + O(δt 3 )
Hence, assuming the operator A and its derivatives ∂A ∂t , ∂A ∂u are locally Lipschitz:

||δx -δx|| ≤ Cδt + Cδt 2 2 ||E n k-1 || + Cδt 3
We recall:

γ F = ∆T F max 0≤n≤N (1 + ||u(T n )||).
Since, the fine solver is based on a two-step time then F ≈ δt 2 . Hence, we neglect in the sequel the contribution Cδt 3 :

||E n+1 k || ≤ β||E n k || + (α + µ + Cδt + Cδt 2 2 )||E n k-1 || + C||δE n k-1 || + γ F
In the sequel, we derive an upper bound for the error terms

E n,N f -1 k and δE n+1 k = E n,N f -1 k -E n+1 k . δE n+1 k = E n,N f -1 k -E n+1 k δE n+1 k = u n,N f -1 k -ε(T 0 , T n+1 -δt -T 0 , u 0 ) -u n+1 k + ε(T 0 , T n+1 -T 0 , u 0 )
In the sequel, we add and substract the following quantites to δE n+1 k :

F (T n , ∆T -δt, u(T n ) -δx, u(T n )) F (T n , ∆T, u(T n ) -δx, u(T n )) δE n+1 k = F (T n , ∆T -δt, u(T n ) -δx, u(T n )) -ε(T n , ∆T -δt, u(T n )) -(F (T n , ∆T, u(T n ) -δx, u(T n )) -ε(T n , ∆T, u(T n ))) + F (T n , ∆T -δt, u n k-1 -δx, u n k-1 ) -F (T n , ∆T -δt, u(T n ) -δx, u(T n )) -F (T n , ∆T, u n k-1 -δx, u n k-1 ) -F (T n , ∆T, u(T n ) -δx, u(T n ))
Taking norms and using (3.10), (3.11), we derive:

||δE n+1 k || ≤ 2∆T F (1 + ||u(T n )||) + Cδt||u n-1,N f -1 k-1 -u(T n -δt) -(u n k-1 -u(T n ))|| + Cδt||u n k-1 -u(T n )|| ||δE n+1 k || ≤ Cδt||δE n k-1 || + Cδt||E n k-1 || + 2γ F E n,N f -1 k = u n,N f -1 k -ε(T 0 , T n+1 -δt -T 0 , u 0 ) = F (T n , ∆T -δt, u n-1,N f -1 k-1 , u n k-1 ) + u n+1 k -F (T n , ∆T, u n-1,N f -1 k-1 , u n k-1 ) -ε(T n , ∆T -δt, u(T n ))
In the sequel, we add and substract the same quantites to E n,N f -1 k as those for the term

δE n+1 k . E n,N f -1 k = u n+1 k -ε(T n , ∆T, u(T n )) + F (T n , ∆T -δt, u(T n ) -δx, u(T n )) -ε(T n , ∆T -δt, u(T n )) -(F (T n , ∆T, u(T n ) -δx, u(T n )) -ε(T n , ∆T, u(T n ))) + F (T n , ∆T -δt, u n k-1 -δx, u n k-1 ) -F (T n , ∆T -δt, u(T n ) -δx, u(T n )) -F (T n , ∆T, u n k-1 -δx, u n k-1 ) -F (T n , ∆T, u(T n ) -δx, u(T n ))
Taking norms and using (3.10), (3.11), we derive:

||E n,N f -1 k || ≤ β||E n k || + (C + Cδt)||δE n k-1 || + (α + µ + 2Cδt + Cδt 2 2 )||E n k-1 || + 3γ F
We summarise the obtained inequalities:

||E n+1 k || ≤ β||E n k || + (α + µ + Cδt + Cδt 2 2 )||E n k-1 || + C||δE n k-1 || + γ F (3.19) ||E n,N f -1 k || ≤ β||E n k || + (α + µ + 2Cδt + Cδt 2 2 )||E n k-1 || + C(1 + δt)||δE n k-1 || + 3γ F (3.20) ||δE n+1 k || ≤ Cδt||δE n k-1 || + Cδt||E n k-1 || + 2γ F (3.
||E n k || ≤ β||E n-1 k || + (α + µ + Cδt + Cδt 2 2 )||E n-1 k-1 || + C 2 δt k j=2 (Cδt) j-2 ||E n-j k-j || + C(Cδt) k-1 ||δE n-k 0 || + 1 + 2C 1 -(Cδt) k-1 1 -(Cδt) γ F (3.22)
The governing term in the sum

C 2 δt k j=2 (Cδt) j-2 ||E n-j k-j || is the term C 2 δt||E n-2 k-2 ||.
To ensure that it does not dominate the term (α + µ + Cδt + Cδt 2 2 )||E n-1 k-1 ||, we suppose that the fine time step verifies: δt ≤ ∆T 2 2 G (see hypothesis (3.17)). In the sequel, we show that the residual terms δE n-k

0 , 1 + 2C 1 -(Cδt) k-1 1 -(Cδt)
γ F and all the terms of the sum for j ≥ 3 can be distributed over the terms:

|| Ẽn k ||, || Ẽn-1 k ||, || Ẽn-1 k-1 || and || Ẽn-2 k-2 ||, where || Ẽn k || is a perturbation of ||E n k ||.
Setting the error perturbation to:

|| Ẽn k || = ||E n k || + ||δE n k || + Cδt ||E n-1 k-1 || + ||δE n-1 k-1 || + 2C + 3 + 2C 2 δt β -1 + α + 3µ + Cδt + Cδt 2 2 + C 2 δt γ F (3.23) The following inequality is satisfied by || Ẽn k ||: || Ẽn k || ≤ β|| Ẽn-1 k || + α|| Ẽn-1 k-1 || + C 2 δt|| Ẽn-2 k-2 || (3.24)
where the constant α is defined in (3.25). Hence, for k = 0, we have:

|| Ẽn 0 || ≤ β|| Ẽn-1 0 || + γG
Then, for k = 1, we have:

|| Ẽn 1 || ≤ β|| Ẽn-1 1 || + α|| Ẽn-1 0 || + γF
where: In the sequel, we use the following notations:

                 α = α + 3µ + 3Cδt + Cδt 2 2 γG = γ G + Cδt + 2C + 3 + 2C 2 δt β -1 + α + 3µ + Cδt + C δt 2 2 + C 2 δt γ F γF = 2C + 3 + 2C 2 δt β -1 + α + 3µ + Cδt + C δt 2 2 + C 2 δt + 3 γ F + Cδt + 2Cδt 2 (3.25)
Following [START_REF] Gander | Nonlinear convergence analysis for the parareal algorithm[END_REF], we consider the sequence (ẽ n k ) n≥0,k≥0 defined recursively as follows. For k = 0:

ẽn 0 = 0 , if n = 0 βẽ n-1 0 + γG , if n ≥ 1 (3.26)
For k = 1:

ẽn 1 = 0 , if n = 0, 1 βẽ n-1 1 + αẽ n-1 0 + γF , if n ≥ 2 (3.27)
For k ≥ 2:

ẽn k = 0 , if n = 0, 1, 2 βẽ n-1 k + αẽ n-1 k-1 + C 2 δtẽ n-2 k-2 , if n ≥ 3 (3.28) Since ||E n k || ≤ || Ẽn k || ≤ ẽn k , for k ≥ 0, n = 0, • • • , N
, we analyse the behavior of (ẽ n k ) to derive a bound for Ẽn k . For this, we consider the generating function:

ρk (ξ) = n≥0 ẽn k ξ n
From (3.26), (3.27) and (3.28) we get:

             ρ0 (ξ) = γG ξ (1 -βξ)(1 -ξ) ρ1 (ξ) = αξ 1 -βξ ρ0 (ξ) + γF ξ (1 -βξ)(1 -ξ) ρk (ξ) = αξ 1 -βξ ρk-1 (ξ) + C 2 δtξ 2 1 -βξ ρk-2 (ξ), k ≥ 2 (3.29)
From which we derive, for k ≥ 1:

ρk (ξ) = γG αk ξ k+1 (1 -ξ) [k/2] j=0 (C 2 δt) j α2j k -j j 1 (1 -βξ) k+1-j +γ F αk-1 ξ k (1 -ξ) [k-1/2] j=0 (C 2 δt) j α2j k -1 -j j 1 (1 -βξ) k-j (3.30) For k = 0, we have: ρ0 (ξ) = γG ξ   p≥0 ξ p     p≥0 β p ξ p   = γG p≥0 p l=0 β l ξ p+1
By a change of variable p = n -1, we obtain:

ẽn 0 = γG n-1 l=0 β l ≤ γG γ G e CT max n (1 + ||u(T n )||)T e -C∆T G , n ≥ 1
For k ≥ 1, using the binomial expansion in (3.30):

1 (1 -βξ) k+1-j = p≥0 k -j + p p β p ξ p
and by a change of variable, we obtain:

n≥0 ẽn k ξ n = γG αk n≥k+1 K n-k-1 ξ n + γF αk-1 n≥k K n-k ξ n
Identifying the term ξ k in the expansion yields to:

ẽk k = γF αk-1 K 0
This gives an upper bound for the error terms || Ẽk k ||, k ≥ 1. We do not use this estimate since the parareal algorithm ensures u n k = F (T 0 , T n -T 0 , u 0 ) for k ≥ n, which yields:

||E k k || = O( F ), k ≥ 1
In the sequel, we identify the terms ξ n for n ≥ k + 1 in the expansion:

ẽn k = γG αk K n-k-1 + γF αk-1 K n-k
We now compute the terms

K p = p l=0 [k/2] j=0 (C 2 δt) j (α + 3µ + Cδt) 2j k -j j k -j + l l β l and K p = p l=0 [k-1/2] j=0 (C 2 δt) j α2j k -1 -j j k -1 -j + l l β l Using: k -j + l l ≤ k + l l and: C 2 δt (α + 3µ + Cδt) 2 ≤ 1, from hypothesis (3.17): δt ≤ ∆T 2 2
G . We have:

K p ≤ p l=0 k + l l β l [k/2] j=0 k -j j ≤ f k k + 1 + p p β p K p ≤ p l=0 k -1 + l l β l [k-1/2] j=0 k -1 -j j ≤ f k-1 k + p p β p
where f k is the general term of the Fibonacci sequence defined by

f 0 = f 1 = 1 and f k+1 = f k + f k-1 , k ≥ 1: f k = [k/2] j=0 k -j j = (1 + √ 5) k+1 -(1 - √ 5) k+1 2 k+1 √ 5
Hence, we derive the bound:

||E n 0 || ≤ || Ẽn 0 || ≤ ẽn 0 ≤ γG γ G e CT max n (1 + ||u(T n )||)T e -C∆T G , n ≥ 1 ||E k k || = O( F ), k ≥ 1 ||E n k || ≤ || Ẽn k || ≤ ẽn k ≤ γG αk f k n k + 1 β n-k-1 + γF αk-1 f k-1 n k β n-k , n ≥ k + 1, k ≥ 1 (3.
31) which ends the proof of the theorem.

Advantages of the multi-step parareal algorithm

We proposed in the last section a new variant of the parareal algorithm with a consistent approximation of the solution at time T nδt in a non intrusive way. The initialisation of the fine propagation in each time window has te be appropriately chosen because an initialisation error would be propagated over the whole time interval and would prevent the parareal algorithm to converge towards the target solution. Another option to treat this issue is to use a one-step time scheme or a multi-stage Runge Kutta method to initialize the fine computation. This option is intrusive since we have to implement new time scheme for the initialisation. Morevover, we will see in section 3.4 that this strategy prevents the parareal to converge to the numerical solution with the target accuracy since the first-order scheme error will dominate. This method adds consistency with the fine scheme. Also, this strategy can be applied to multistep time schemes involving several fine time steps preceding the time T n by applying the same correction to terms taking the form:

u n,N f -i k+1 , i = 1, • • • , I.
We detail the algorithm for a multi-step time scheme involving more than one fine time step preceding the time T n .

                 u n+1 0 = G(T n , ∆T, u n 0 ), 0 ≤ n ≤ N -1 u n+1 k+1 = G(T n , ∆T, u n k+1 ) + F (T n , ∆T, u n-1,N f -1 k , u n k ) -G(T n , ∆T, u n k ), 0 ≤ n ≤ N -1, k ≥ 0 u n,N f -i k+1 = F (T n , ∆T -iδt, u n-1,N f -I k , u n-1,N f -I+1 k , • • • , u n k ) + u n+1 k+1 -F (T n , ∆T, u n-1,N f -I k , u n-1,N f -I+1 k , • • • , u n k ), i = 0, • • • , I, 0 ≤ n ≤ N -1, k ≥ 0 (3.32)
where we denote F (t, s, w 1 , w 2 , • • • , w I ) the multi-step propagator for any given time t ∈ [0, T ], s ∈ [0, Tt] and any function w 1 , • • • , w I ∈ U takes I initial values at times t and tiδt and propagates it at time t + s, where δt is the fine time step. We illustrate the good convegence properties in the next section by applying the parareal algorithm to an ODE system solved by a coarse solver based on a one-step time scheme and a fine solver based on a third-order BDF method.

When the coarse solver is a multi-step time scheme, there exists several options to initialise it on each time window:

If the coarse time step δT is equal to the size of the time window ∆T , there is no additional correction in the parareal algorithm since the solution at every coarse time step are updated If δT < ∆T , there are intermediate coarse time iterations in each time window. In [START_REF] Audouze | Symplectic multi-time step parareal algorithms applied to molecular dynamics[END_REF], the initialisation of the coarse solver is addressed and the authors propose a parareal-type correction at time T n -δT :

u n+1-N c int k+1 = G(T n , ∆T -δT, u n-N c int k+1 , u n k+1 ) + F (T n , ∆T -δT, u n-1,N f -1 k , u n k ) -G(T n , ∆T -δT, u n-N c int k , u n k ), N c int = δT δt 0 ≤ n ≤ N -1, k ≥ 0 (3.33)
We illustrate the behavior of the full multi-step parareal algorithm with specific initialisation of the fine and coarse solver in the next section where the two solvers involved in the parareal method are the second-order BDF method.

Numerical tests

We apply the multi-step parareal algorithm to a simple ODE firstly, the damped oscillator and then to a stiff problem, the Brusselator system. Our results illustrate that our approach improves the convergence properties with respect to the classical parareal algorithm. We also show that the generalisation of this approach to third-order time schemes holds and the convergence properties derived in Theorem (3.3) are preserved. Finally, we address the question of the parallel efficiency of the multi-step parareal. In the last section, we apply the adaptive parareal algorithm (see [START_REF] Maday | An adaptive parareal algorithm[END_REF]) where the accuracy of the fine solver is increased accross the iterations.

Numerical convergence results

The damped oscillator

We consider the damped oscillator system:

u (t) + 2λu (t) + ω 2 0 u(t) = 0, t ∈ (0, T ), with u(0) = u 0 , T = 10,
We rewrite it as a first order ODE system:

X (t) = AX(t), t ∈ (0, T ), with X(0) = u 0 u 0 , A = 0 1 -ω 2 0 -2λ
This system models the dynamic of a simple nonstiff harmonic oscillator under a frictional force.

For our tests, we set λ = 0.05, ω 0 = 1 and X 0 = 0.1 0.2 .

The coarse solver is a Backward Euler method with a coarse time step: ∆T = 0.1 which corresponds to 100 time windows and the fine solver is a second-order BDF method with a fine time step δt = 10 -4 (respecting hypothesis (3.17)). In figure 3.1, the fine solver is based on a two-step time scheme where the computation of the solution u n,j+1 at time T n + (j + 1)δt depends on the solutions u n,j and u n,j-1 at times T n + jδt and T n + (j -1)δt, respectively. We use the multi-step variant of parareal (3.13) to initialise the fine solver in each time window, starting from the parareal iteration k ≥ 2. At the parareal iteration k = 1, we use a Backward Euler method to initialise the fine solver since we did not use the fine propagator yet.

In this section, we analyse the evolution of two different errors accross the parareal iterations:

the error between the fine solution computed in a sequential way and the parareal solution in L ∞ (0, T ) norm,

max 1≤n≤N ||u n k -F (T 0 , T n -T 0 , u 0 )|| (3.34)
the error between the exact solution and the parareal solution in L ∞ (0, T ) norm

max 1≤n≤N ||u n k -u(T n )|| (3.35)
In all the figures of this section, we plot the evolution of errors (3.35-3.34) in the two following cases:

Without a multi-step adaptation (red curve): the error between the parareal solution where the Backward Euler method is used at each iteration for the initialisation of the fine solver and the fine solution computed in a sequential way for (3.34) (the exact solution for (3.35)), on one hand, With a multi-step adaptation (blue curve): the error between the solution given by the multistep parareal algorithm the fine solution computed in a sequential way for (3.34) (the exact solution for (3.35)), on the other hand.

In figure 3.1, we see that without the multi-step adaptation the error (3.34) stagnates around 10 -6 without recovering the fine solution at the machine precision, even after 100 iterations. On the other hand, using the multi-step parareal algorithm, the error continues to decrease until reaching the machine precision. Moreover, in the right figure, we see that the only way to recover the correct approximation of the exact solution is to use a multi-step adaptation, otherwise, without adaptation, the parareal algorithm will not reach the target accuracy. This result shows that making an initialisation error for a multi-step fine solver will prevent the parareal algorithm to obtain the approximation of the exact solution with the desired accuracy.

The convergence properties are illustrated in figure 3.1 on a fine solver based on the second-order BDF method with time step δt = 10 -4 .

In figure (3.2), we apply the extension of the multi-step parareal algorithm (3.32) to three-step time schemes by giving a consistent approximation of the solutions u(T nδt) and u(T n -2δt). We illustrate the convergence properties of this strategy by applying it on a fine solver based on the third-order BDF method with a time step δt = 10 -4 (see figure 3.2). We observe the same behaviour of the errors (3.35-3.34): without a multi-step adaptation, the fine propagation is initialised by two Backward Euler iterations and does not allow to recover the target approximation of the exact solution while the multi-step parareal converges to the exact solution with the desired accuracy.

The Brusselator system

We consider the Brusselator system:

x = A + x 2 y -(B + 1)x y = Bx -x 2 y
with initial condition x(0) = 0 and y(0) = 1. This is a stiff ODE that models a chain of chemical reactions. It was already studied in previous works on the parareal algorithm ( [START_REF] Gander | Nonlinear convergence analysis for the parareal algorithm[END_REF][START_REF] Maday | An adaptive parareal algorithm[END_REF]). The system has a fixed point at x = A and y = B A which becomes unstable when B > 1 + A 2 and leads to oscillations. We place ourselves in this oscillatory regime by setting A = 1 and B = 3. The dynamics present large velocity variations in some time sub-intervals, making the use of high order time schemes particularly desirable for an appropriate treatment of the transient. The coarse solver is a Backward Euler method with a coarse time step: ∆T = 0.1 which corresponds to 180 time windows since T = 18. The fine solver is a second-order BDF method with a fine time step δt = 10 -4 (respecting hypothesis (3.17)). In figure 3.3, the fine solver is based on a two-step time scheme. We use the multi-step parareal algorithm (3.13) to initialise the fine solver in each time window.

Likewise, we analyse the evolution of the errors (3.34) and (3.35) accross the parareal iterations. In Figure 3.3, we see that without the multi-step adaptation the error (3.34) stagnates around 10 -6 without recovering the fine solution at the machine precision, even after 180 iterations. On the other hand, using the multi-step parareal algorithm, the error continues to decrease until reaching the machine precision. Moreover, in the right figure, we see that the only way to recover the correct approximation of the exact solution is to use a multi-step adaptation, otherwise, without adaptation, the parareal algorithm will not reach the target accuracy. This result shows that making an initialisation error for a multi-step fine solver will prevent the parareal algorithm to obtain the approximation of the exact solution with the desired accuracy.

The convergence properties are illustrated in figure 3.3 on a fine solver based on the second-order BDF method with time step δt = 10 -4 . In the figure (3.4), we apply the extension of the multi-step parareal algorithm (3.32) to threestep time schemes by giving a consistent approximation of the solutions u(T n -δt) and u(T n -2δt). We illustrate the convergence properties of this strategy by applying it on a fine solver based on the third-order BDF method with time steps δt = 10 -4 (see figure 3.4). We observe the same behavior of the errors (3.35-3.34): without a multi-step adaptation, the fine propagation is initialised by two Backward Euler iterations and does not allow to recover the target approximation of the exact solution while the multi-step parareal converges to the exact solution with the desired accuracy.

Parallel efficiency

We address in this section the question of the speed up performances for the multi-step parareal algorithm. The only additional operations in the multi-step variant compared to the classical parareal are the corrections of solutions involved in the initialisation of the fine solver in each time window (update of u n,N f -1 k+1 in (3.13) for example). Hence, we consider that the computational cost of the multi-step variant is the same as the one of the classical parareal. In a recent work [START_REF] Maday | An adaptive parareal algorithm[END_REF], the authors propose a new method, the adaptive parareal algorithm, where the accuracy of the fine solver is increased across the iterations. This new point of view improves the speed up performances of the parareal method and the only remaining factor limiting high performance becomes the cost of the coarse solver. In this section, we seek to improve the parallel efficiency of the multi-step parareal method by increasing the accuracy of fine solver at each iteration. We first recall the parallel efficiency for the classical parareal (CP) and the adaptive parareal (AP) to obtain a solution with accuracy η and a propagation over [0, T ]:

ef f CP (η, [0, T ]) ∼ 1 K(η) ef f AP (η, [0, T ]) ∼ 1 1 + 1/α G
, under the hypothesis of Proposition 3.1 in [START_REF] Maday | An adaptive parareal algorithm[END_REF] where K(η) is the number of parareal iterations to obtain the approximation of the exact solution with the target accuracy η and α, the order of the fine time scheme. To apply this approach the method is then defined as the ratio of the above speed up with the number of processors which gives a target of 1 to any parallel solver:

ef f (η, [0, T ]) := S(η, [0, T ]) N .
To compare the speed-up of the multi-step and adaptive parareal algorithms, we use the number of fine and coarse propagations involved in the numerical solution and the computational cost of the coarse and fine propagations (communication delays have not been taken into account). For example, in table 3.1, the cost of the multi-step parareal algorithm is equal to the cost of 9 coarse propagations over [0, T ] plus 8 fine propagations over [T n , T n+1 ] with a fine time step δt = 10 -4 .

In [START_REF] Maday | An adaptive parareal algorithm[END_REF], the authors show that the main element affecting the performance of the adaptive parareal method is no longer the cost of the fine solver but the cost of the coarse solver. Hence, we compare the speed-up and efficiency when we count or do not count the cost of the coarse solver in table 3.2.

Obviously, when we do not count the cost of the coarse solver, the performance of both algorithms improves.

Conclusion

We have built a new variant of the parareal algorithm allowing to overcome the issue of initialising the fine and the coarse solvers when they are based on a multi-step time scheme ( [5]). The convergence properties of the multi-step parareal are very close to that of the classical parareal algorithm in the case of two-step time schemes. An extension of our approach to generic multi-step time schemes is proposed and validated numerically on a three-step time scheme. In addition, the accuracy of the multi-step parareal algorithm is illustrated on the numerical solution of a stiff ODE such as the Brusselator system. Finally, we address the question of the parallel efficiency of our strategy by coupling it with the adaptive parareal algorithm proposed in [START_REF] Maday | An adaptive parareal algorithm[END_REF]. The new adaptive formulation of the parareal algorithm opens the door to improve significantly the parallel efficiency of the method provided that the cost of the coarse solver is moderate.

Introduction

As an introduction to the issue, we consider a 1D conservative non linear hyperbolic system

∂ t U (x, t) + ∂ x F (U )(x, t) = 0, (4.1)
with unknown vector U ∈ R m and Lipschitz flux F : R m → R m with real-diagonalisable Jacobian matrix

A(U ) = ∇ U F (U ) ∈ R m×m .
In the rest of the introduction we review some notions about the numerical diffusion.

Consistency analysis

When approximating smooth solutions U of (4.1) by a consistent numerical method on a regular mesh with space step x, the semi-discrete equations approximate to the first order in x the following perturbed version of equation (4.1) :

∂ t U ∆x + ∂ x F (U ∆x ) = D(U ∆x , x) + o( x), (4.2)
where U ∆x is the numerical solution and D is a second order differential operator. When the flux function F is linear, linear numerical methods yield a linear diffusion operator D(U, x) = x∂ x (D∂ x U ). The matrix D comes from the upwind (off-centered) contributions of the discrete equations and gives a first insight into the scheme precision and stability. In the non linear case (F Lipschitz), the numerical diffusion operator can often be approximated to the first order by a non linear diffusion operator :

D(U, x) = x∂ x (D(U )∂ x U ) + o( x). (4.3)
This is the case for instance for colocated schemes based on characteristic upwinding such as Godunov [START_REF] Godunov | A difference scheme for numerical solution of discontinuous solution of hydrodynamic equations[END_REF], Roe [START_REF] Roe | Approximate riemann solvers, parameter vectors and difference schemes[END_REF], VFRoe [START_REF] Masella | On an approximate godunov scheme[END_REF] or VFFC [START_REF] Ghidaglia | Une méthode volumes finis à flux caractéristiques pour la résolution numérique des systèmes hyperboliques de lois de conservation[END_REF] schemes where the non linear numerical diffusion tensor is

D(U ) = |A(U )|.
We recall that in the case of symmetric hyperbolic systems ( t A(U ) = A(U )), any entropy solution to (4.1) preserves the L 2 norm (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] Example 3.2 in the Introduction chapter) and we would like the discrete L 2 norm of any scheme to be bounded as well. In the case of non-symmetric systems, one first symmetrises the system using entropic variables V (U ) = ∇s(U ) where s a strictly convex entropy of the system (4.1) is assumed to exist (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] Theorem 3.2 in the introduction chapter). The new symmetric system:

∂ t V + Ā(V )∂ x V = 0, with t Ā = Ā (4.4)
is linearly L 2 -stable. Any numerical scheme yields a numerical diffusion D(V ) in the symmetrised basis and we require that the diffusion operator D have positive symmetric part : t D + D ≥ 0.

The operator D gives a first insight into the scheme precision since the smaller the operator D, the closer the approximate solution U ∆x is to U . However the numerical diffusion operator gives also important informations about the scheme stability.

Stability analysis

In the case of symmetric hyperbolic systems ( t A(U ) = A(U )), the exact equation (4.1) yields the conservation of the L 2 norm (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] Example 3.2 in the Introduction chapter)

∀t ∈ R + , ∂ t R ||U || 2 2 (x, t)dx = 0, (4.5)
whilst the perturbed equation (4.2) yields most of the time the first order estimate

∀t ∈ R + , ∂ t R ||U ∆x || 2 2 (x, t)dx = x R t U ∆x ∂ x (D(U ∆x )∂ x U ∆x )dx + o( x) = -x R t (∂ x U ∆x )D(U ∆x )(∂ x U ∆x )dx + o( x), (4.6)
In order to obtain an L 2 stable scheme, it is therefore usual to require that the diffusion operator D have positive symmetric part : t D + D ≥ 0.

In the case of non symmetric systems, one first symmetrises the system using entropic variables ξ(U ) = ∇ U s(U ) where s is an entropy of the system. The new system is entropic and yields

∀t ∈ R + , ∂ t R ||ξ(U )|| 2 2 (x, t)dx = 0, (4.7)
We remark that a scheme that is entropic ∂ t R s(U ) ≤ 0 is not necessarily stable since s is not necessarily bounded below as is the case with the full Euler system. We recall the expression of the entropy for the full Euler system:

s = C v ln ρE - q 2 2ρ -γ ln ρ , (4.8)
where C v is the specific heat and γ the adiabatic constant such that: p = (γ -1)ρe, e is the internal energy and E, the total energy. If a numerical scheme applied to the full Euler system is entropic then the quantity (4.8) decays. This does not imply the boundedness of the unknowns ρ, u, E. For example, assuming ρ is constant, the variables E and q can grow infinitely while maintaining the difference ρE -q 2 2ρ constant.

We investigate in this chapter the L 2 -stability of staggered schemes. In this first account of our research, we investigate the isentropic Euler system which raises an issue that will remain with more complex fluid model : the numerical treatment of the mass balance equation and of the momentum equation yields a non classical diffusion operator. In order to obtain a straightforwardly stable scheme we propose a new discretisation with positive numerical diffusion. In section 4.2, we determine the numerical diffusion of the staggered schemes and show it does not straightforwardly yield a linear stability. We then present a new class of staggered schemes and prove their linear stability in section 4.3. Some numerical results are given in section 4.4.

The numerical diffusion of staggered schemes for the Euler system

We address here the following system, the isentropic Euler equations, written in the following conservative form :

∂ t ρ + ∂ x (ρu) = 0 ∂ t (ρu) + ∂ x (ρu 2 ) + ∂ x p = 0 (4.9)
This problem is posed over an open bounded connected subset Ω of R, with boundary ∂Ω, and a finite time interval (0, T ). The variable t stands for the time, ρ, u and p are the density, velocity and pressure in the flow.

The results can be extended to the multidimensional Euler system but the calculations are lengthier and do not help the intuition. The Euler system (4.9) can take the non conservative form

∂ t U + A(U )∂ x U = 0 (4.10) with c 2 = ∂p ∂ρ , u = q ρ and U = ρ q , A(U ) = 0 1 c 2 -u 2 2u . ( 4

.11)

Remark 4.1. In the sequel, we seek to show the upwind matrices for numerical schemes based on the principle of vector upwinding using the eigenbasis of the Jacobian matrix A such as: Godunov [START_REF] Godunov | A difference scheme for numerical solution of discontinuous solution of hydrodynamic equations[END_REF], Roe [START_REF] Roe | Approximate riemann solvers, parameter vectors and difference schemes[END_REF] and VFRoe [START_REF] Masella | On an approximate godunov scheme[END_REF]. The upwind matrices have two arguments: a left state U L and a right state U R . Here, we illustrate their behaviour with the Roe scheme. When the right and left states are equal, the upwind matrix has the following expression:

D upw (U, U ) = |A(U )| = 1 c c 2 -u 2 u u(c 2 -u 2 ) u 2 + c 2 .
(4.12)

If the flux F in 4.1 is linear, the upwind matrices of the Godunov and VFRoe schemes have the same form (4.12). We remark that by a change of basis, we obtain a system of decoupled transport equation where the transport speed is positive, hence ensuring the L 2 -stability for each equation and thus for the original problem.

Also for low Mach numbers ( |u| c 1), for a fixed velocity u and the sound speed c tending to infinity, the upwind matrix of the Roe scheme converges towards the identity matrix:

D upw = cId + O( 1 c ) (4.13)
The diffusion is evenly distributed on the mass and momentum equations. We see that there is problem in the order of magnitude:

Considering the mass equation of the system (4.9), the discretisation introduced a perturbation of the order ρc in the right hand side while the left hand side is of the order ρu:

∂ t ρ + ∂ x q = c x∂ xx ρ + o( x). (4.14)
Hence, the numerical scheme is too diffusive for the mass equation.

Considering the momentum equation of the system (4.9), the discretisation introduced a perturbation of the order ρuc in the right hand side while the left hand side is of the order ρc 2 :

∂ t q + ∂ x q 2 ρ + c 2 ∂ x ρ = c x∂ xx q + o( x). (4.15)
Hence, the numerical scheme is not diffusive enough for the mometum equation.

The upwind type schemes can be proven to be linearly L 2 stable ( [START_REF] Godunov | A difference scheme for numerical solution of discontinuous solution of hydrodynamic equations[END_REF], [START_REF] Roe | Approximate riemann solvers, parameter vectors and difference schemes[END_REF], [START_REF] Masella | On an approximate godunov scheme[END_REF]). However the amount of numerical diffusion is proportional to the sound speed c and for low Mach number flows, the schemes based on characteristics upwinding are not able to capture nearly incompressible solutions (see [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF][START_REF] Dellacherie | Analysis of godunov type schemes applied to the compressible euler system at low mach number[END_REF] for more details).

On the contrary, staggered schemes are known to be more precise for low Mach number flows in practice and are very popular in the thermal hydraulics community ( [START_REF] Prosperetti | Computational methods for multiphase flow[END_REF]). However their stability analysis is historically based on heuristics ( [START_REF] Hirt | Heuristic stability theory for finite difference equations[END_REF]). Yet the conservative staggered schemes presented in [START_REF]On some implicit and semi-implicit staggered schemes for the shallow water and euler equations[END_REF][START_REF] Herbin | Staggered schemes for all speed flows[END_REF] are proven to be entropic and to satisfy a kinetic energy preservation [START_REF] Herbin | A kinetic energy preserving convection operator for the mac discretization of compressible navier-stokes equations[END_REF]. Likewise in [START_REF] Berthelin | Kinetic schemes on staggered grids for barotropic euler models : entropy-stability analysis[END_REF], the authors present a kinetic scheme on staggered grids for the barotropic Euler equations, derive stability conditions which preserve both the positivity of the density and the decay of the discrete global entropy, and satisfy a kinetic energy preservation. Unfortunately the boundedness of the entropy does not necessarily imply the boundedness of the solution. Indeed a strictly convex function is not necessarily bounded below. This is in particular the case for the full Euler system since the entropy involves the functionln which is strictly convex but not bounded below (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] Example 3.3 in the Introduction chapter). In the next subsection we show that the first order perturbed equation (4.2) associated to staggered schemes yields not the classical diffusion operator (4.3) but instead a strongly nonlinear numerical diffusion operator.

The staggered scheme of Herbin et al.

Using staggered schemes, the density and pressure are located on cells and the velocity on faces (nodes in 1D) [START_REF]A numerical fluid dynamics calculation method for all flow speeds[END_REF][START_REF] Harlow | Numerical calculation of almost incompressible flow[END_REF]. The momentum variable is usually split as a product between the density and the velocity : q = ρ u. The main difference between the various staggered schemes is the treatment of the convection term ρ u ⊗ u in the momentum equation. We consider the staggered scheme of [START_REF] Herbin | A kinetic energy preserving convection operator for the mac discretization of compressible navier-stokes equations[END_REF] as a prototype of staggered schemes. Indeed, the main results of this section extends to other staggered schemes. This is because the mass discretisation is the same in all staggered schemes and therefore the mass diffusion operator will be non classical. In the past decade, Herbin, Latché and their coauthors have proposed a new approach with rigorous proofs of stability: discrete inequality for the kinetic energy and entropic character. They discretise the conservative form of the Euler equations (eq:euler syst cons 1D) with a conservative scheme. The different variants include one step ([65] section 2.1, [66] section 3.1) and prediction/correction steps ([65] section 2.2, [66] section 4.1) variants, fully implicit ([66] section 3, [START_REF] Herbin | Staggered schemes for all speed flows[END_REF] section 2.1), semi implicit and almost explicit [START_REF]On some implicit and semi-implicit staggered schemes for the shallow water and euler equations[END_REF] (all but the pressure gradient are explicit-in-time) variants.

For simplicity we present the discrete equation of the fully implicit variant ([66] section 3, [65] section 2.1) for the 1D isentropic Euler equations in conservative form.

ρ n+1 i -ρ n i t + ρ up,n+1 i+ 1 2 u n+1 i+ 1 2 -ρ up,n+1 i-1 2 u n+1 i-1 2 x = 0 (4.16) ρn+1 i+ 1 2 u n+1 i+ 1 2 -ρn i+ 1 2 u n i+ 1 2 t + ρu n+1 i+1 u up,n+1 i+1 -ρu n+1 i u up,n+1 i x + p n+1 i+1 -p n+1 i x = 0. ( 4.17)
The pressure p i and the density ρ i are located at the cell centers whereas the velocity u i+ 1 2 are located at the cell interfaces. The expression of the products ρu, ρ∂ t u and ρu 2 between the velocity located at cell interfaces and the density located at cell centers thus has to be defined through interpolation formula.

The mass flux ρu at the cell interfaces is defined using an upwind density ρ up i+ 1 2 defined as :

ρ up i+ 1 2 =    ρ i if u i+ 1 2 > 0 ρ i+1 if u i+ 1 2 ≤ 0 = ρ i + ρ i+1 2 + sign(u i+ 1 2 ) ρ i -ρ i+1 2 , (4.18)
which is the sum of a centered and an upwind terms. The expression of ρi+ 1 2 in the discrete momentum equation accounts for an average of the neighbouring densities

ρi+ 1 2 = 1 2 (ρ i + ρ i+1 ). (4.19)
The expression of ρu in the discrete momentum equation is

ρu i = 1 2 (ρ up i-1 2 u i-1 2 + ρ up i+ 1 2 u i+ 1 2
).

(4.20)

The upwind velocity u up i at cell centers is defined as :

u up i =    u i-1 2 if ρu i > 0 u i+ 1 2 if ρu i ≤ 0 = u i-1 2 + u i+ 1 2 2 + sign(ρu i ) u i-1 2 -u i+ 1 2 2 , (4.21)
which is the sum of a centered and an upwind terms. It is possible to use a centered velocity ū instead of the upwind velocity u up (see [START_REF]On some implicit and semi-implicit staggered schemes for the shallow water and euler equations[END_REF]).

The numerical diffusion of the scheme

The scheme presented in [START_REF]On some implicit and semi-implicit staggered schemes for the shallow water and euler equations[END_REF] section 3 and [65] section 2.1 is proven to be entropic. However the boundedness of the entropy does not necessarily imply the boundedness of the solution. Indeed a convex function is not necessarily bounded below. This is in particular the case for the full Euler system (not the isentropic one) because the entropy involves the functionln which is convex but not bounded below.

We would like to study the L 2 stability of the scheme by analysing its numerical diffusion operator. This would give a new insight into the scheme and prove at least a linear stability and a nonlinear stability in the case of almost constant initial data.

In this section, we assume that the exact solution is smooth and we determine the numerical diffusion of the scheme (4.16-4.17). In the context of staggered schemes, we have two meshes: one for the cell centered unknowns at points x i , i = 1, • • • , N and the other for the unknowns defined on edges at points x i+1/2 . Then we could write the consistency error for the momentum equation around x = x i or x i+1/2 . To stay in the spirit of the staggered schemes, we first develop the analysis around x = x i+1/2 then in a second time around x = x i to be closer to the classical analysis of the consistency error. To obtain the consistency error for the mass equation, we perform the following Taylor expansions around x = x i , assuming smooth solutions with u = 0 :

ρ i-1 = ρ(x i ) -∆x∂ x ρ(x i ) + 1 2 (∆x) 2 ∂ xx ρ(x i ) + O(∆x 3 ) (4.22) ρ i+1 = ρ(x i ) + ∆x∂ x ρ(x i ) + 1 2 (∆x) 2 ∂ xx ρ(x i ) + O(∆x 3 ) (4.23) u i-1 2 = u(x i ) - ∆x 2 ∂ x u(x i ) + 1 2 ∆x 2 2 ∂ xx u(x i ) + O(∆x 3 ) (4.24) u i+ 1 2 = u(x i ) + ∆x 2 ∂ x u(x i ) + 1 2 ∆x 2 2 ∂ xx u(x i ) + O(∆x 3 ) (4.25)
Mass numerical diffusion From (4.16), (A.1) and (A.2) the discrete mass flux is (we omit the time indices)

ρ up i+ 1 2 u i+ 1 2 -ρ up i-1 2 u i-1 2 = ∆x∂ x (ρu)(x i ) - 1 2 (∆x) 2 sign(u(x i ))∂ x (u(x i )∂ x ρ) + O(∆x 3 ).(4.26)
Hence the mass flux consistency is finally

ρ up i+ 1 2 u i+ 1 2 -ρ up i-1 2 u i-1 2 ∆x = ∂ x (ρu)(x i ) - ∆x 2 sign(u(x i ))∂ x (u∂ x ρ)(x i ) + O(∆x 2 ). (4.27)
and the numerical diffusion associated to th mass conservation law is the strongly non linear diffusion term sign(u)∂ x (u∂ x ρ). The linear stability analysis of such a strongly nonlinear diffusion is not classical and we are not aware of any reference.

If we assume that u does not change sign then the diffusion term simplifies to the weakly non linear diffusion term ∂ x (|u|∂ x ρ) which is positive diffusion coefficient |u|. The weakly non linear diffusion term ∂ x (|u|∂ x ρ) can be linearised around a constant state (ρ 0 , u 0 = 0 as:

∂ x (|u 0 |∂ x ρ) + ∂ x (|u|∂ x ρ 0 ) = |u 0 |∂ xx ρ).
Hence if u > 0 or u < 0 the mass equation has a positive contribution on the diagonal of the numerical diffusion tensor D and thus has a stabilising effect. If we allow u to change sign then the multiplication with sign(u) makes things more complicate and we can not rule out potential instabilities. The linearisation is not trivial since the smoothness of u does not imply even the continuity of sign(u). The consistency analysis is only a first step that requires smooth solutions but the final goal of capturing discontinuous weak solutions with velocity that change sign will raise even more questions.

Momentum numerical diffusion For the momentum numerical diffusion we compute the contribution from the pressure

p i+1 -p i = ∆x∂ x p(x i )+ 1 2 (∆x) 2 ∂ xx p(x i )+O(∆x 3 ).p i+1 -p i = ∆x∂ x p(x i+ 1
2

)+O(∆x 3 ). (4.28)

the contribution from the time evolution term (case of the conservative scheme (4.16-4.17))

∂ t (ρ i+ 1 2 u)(x i+ 1 2 , t) = ∂ t (ρu)(x i+ 1 2 , t) + ∂ t ρ(x i+1 ) -2ρ(x i+ 1 2 ) + ρ(x i ) 2 u (x i+ 1 2 , t) (4.29)
= ∂ t (ρu)(x i+ 1 2 , t) + 1 2 ∆x 2 2 ∂ t ((∂ xx ρ)u) (x i+ 1 2 , t) + O(∆x 3 )(4.30)
The time evolution term will bring a perturbation in (∆x) 2 that we can neglect since we are interested in first order error terms in (∆x).

To obtain the consistency error, we perform the following Taylor expansions around x = x i+ 1 2 , assuming smooth solutions with u = 0:

ρ i+1 = ρ(x i+ 1 2 ) + ∆x 2 ∂ x ρ(x i+ 1 2 ) + ∆x 2 8 ∂ xx ρ(x i+ 1 2 ) + O(∆x 3 ) ρ i+2 = ρ(x i+ 1 2 ) + 3∆x 2 ∂ x ρ(x i+ 1 2 ) + 9∆x 2 8 ∂ xx ρ(x i+ 1 2 ) + O(∆x 3 ) ρ i = ρ(x i+ 1 2 ) - ∆x 2 ∂ x ρ(x i+ 1 2 ) + ∆x 2 8 ∂ xx ρ(x i+ 1 2 ) + O(∆x 3 ) ρ i-1 = ρ(x i+ 1 2 ) - 3∆x 2 ∂ x ρ(x i+ 1 2 ) + 9∆x 2 8 ∂ xx ρ(x i+ 1 2 ) + O(∆x 3 ) u i+ 3 2 = u(x i+ 1 2 ) + ∆x∂ x u(x i+ 1 2 ) + 1 2 (∆x) 2 ∂ xx u(x i+ 1 2 ) + O(∆x 3 ) u i-1 2 = u(x i+ 1 2 ) -∆x∂ x u(x i+ 1 2 ) + 1 2 (∆x) 2 ∂ xx u(x i+ 1 2 ) + O(∆x 3 ) (4.31)
Momentum numerical diffusion From (4.17), (A.7) and (A.8), the discrete momentum flux yields (we omit the time indices)

ρu i+1 u up i+1 -ρu i u up i = ∆x∂ x (ρu 2 )(x i+ 1 2 ) - ∆x 2 2 sign(ρu(x i+ 1 2 ))∂ x (ρu∂ x u)(x i+ 1 2 ) - ∆x 2 2 sign(u(x i+ 1 2 ))∂ x (u 2 ∂ x ρ)(x i+ 1 2 ) + O(∆x 3 ) (4.32)
Hence, from (4.32) and (4.28), the momentum flux consistency is finally:

ρu i+1 u up i+1 -ρu i u up i ∆x + p i+1 -p i ∆x = ∂ x (ρu 2 )(x i+ 1 2 ) + ∂ x p(x i+ 1 2 ) - ∆x 2 sign(u(x i+ 1 2 ))∂ x (u 2 ∂ x ρ)(x i+ 1 2 ) + sign(ρu(x i+ 1 2 ))∂ x (ρu∂ x u)(x i+ 1 2 )
+ O(∆x 2 ) (4.33) For smooth solutions u and ρ:

∂ x (ρu∂ x u) = ∂ x (ρu∂ x q ρ ) = ∂ x (u∂ x q -u 2 ∂ x ρ)
We obtain:

ρu i+1 u up i+1 -ρu i u up i ∆x + p i+1 -p i ∆x = ∂ x (ρu 2 )(x i+ 1 2 ) + ∂ x p(x i+ 1 2 ) - ∆x 2 sign(u(x i+ 1 2 )) -sign(ρu(x i+ 1 2 )) ∂ x (u 2 ∂ x ρ)(x i+ 1 2 ) +sign(ρu(x i+ 1 2
))∂ x (u∂ x q)(x i+ 1 2

) + O(∆x 2 ) (4.34) Assuming that the scheme (4.16-4.17) preserves the positivity of ρ, we have:

sign(ρu(x i+ 1 2 )) = sign(u(x i+ 1 2 ))
Hence:

ρu i+1 u up i+1 -ρu i u up i ∆x + p i+1 -p i ∆x = ∂ x (ρu 2 )(x i+ 1 2 ) + ∂ x p(x i+ 1 2 ) - ∆x 2 sign(u(x i+ 1 2 ))∂ x (u∂ x q)(x i+ 1 2 ) + O(∆x 2 ) (4.35)
The numerical diffusion associated to the momentum conservation law is the term sign(u)∂ x (u∂ x q) and is decoupled from the mass diffusion. The linear stability analysis of such a strongly non linear diffusion is not classical and we are not aware of any reference in the litterature. If we assume that u does not change sign then the diffusion term simplifies to the weakly non linear diffusion term ∂ x (|u|∂ x q) which involves a positive diffusion coefficient |u|. The weakly non linear diffusion term ∂ x (|u|∂ x q) can be linearised around a constant state (q 0 , u 0 = 0 as:

∂ x (|u 0 |∂ x q) + ∂ x (|u|∂ x q 0 ) = |u 0 |∂ xx q).
Hence if u > 0 or u < 0 the momentum equation has a positive contribution on the diagonal of the numerical diffusion tensor D and thus has a stabilising effect. If we allow u to change sign then the multiplication with sign(u) makes things more complicated and we can not rule out potential instabilities. The linearisation is not trivial, even taking u smooth enough, since the function sign(u) is not continuous. The consistency analysis is only a first step that requires smooth solutions but the final goal of capturing discontinuous weak solutions with velocity that change sign will raise even more issues. Finally, the staggered scheme (4.16-4.17) have the following numerical diffusion operator:

∂ t U ∆x + ∂ x F (U ∆x ) = D(U ∆x , ∆x) + O(∆x 2 ) D(U ∆x , ∆x) = ∆x 2 sign(u(x i )) 0 0 sign(u(x i+ 1 2 )) ∂ x u(x i ) 0 0 u(x i+ 1 2 ) ∂ x ρ(x i ) q(x i+ 1 2 ) (4.36)
Now that we derived the consistency error around the point x = x i+1/2 for the momentum equation we turn to the classical way to analyse the consistency error around x = x i . the contribution from the pressure

p i+1 -p i = ∆x∂ x p(x i ) + 1 2 (∆x) 2 ∂ xx p(x i ) + O(∆x 3 ). (4.37)
the contribution from the time evolution term (case of the conservative scheme (4.16-4.17))

ρi+ 1 2 u i+ 1 2 = 1 2 (ρ i + ρ i+1 )u i+ 1 2 = (ρu)(x i ) + ∆x 2 ∂ x (ρu)(x i ) + ∆x 2 8 ρ∂ xx u(x i ) + ∆x 2 4 ∂ x (u∂ x ρ)(x i ) + O(∆x 3 ) (4.
38) Since we are interested in first order error terms in (∆x), we neglect the perturbation in (∆x) 2 and obtain:

∂ t (ρ i+ 1 2 u)(x i , t) = ∂ t (ρu)(x i , t) + ∆x 2 ∂ t (∂ x (ρu))(x i , t) + O(∆x 2 ) (4.39)
Since (ρu)(x i , t) is a solution of (4.9):

∂ t (∂ x (ρu))(x i , t) = ∂ x (∂ t (ρu))(x i , t) = -∂ xx (ρu 2 )(x i , t) -∂ xx p(x i , t) ∂ t (ρ i+ 1 2 u)(x i , t) = ∂ t (ρu)(x i , t) - ∆x 2 ∂ xx (ρu 2 )(x i , t) - ∆x 2 ∂ xx p(x i , t) + O(∆x 2 ) (4.40)
We seek to derive the numerical diffusion for the momentum equation around x = x i :

ρu i+1 u up i+1 -ρu i u up i = ∆x ∂ x (ρu 2 )(x i ) + ∆x 2 ∂ xx (ρu 2 )(x i ) - ∆x 2 2 sign(u(x i ))∂ x (u∂ x q)(x i ) + O(∆x 3 )
(4.41) Hence, from (4.41) and (4.37), the momentum flux consistency is finally:

∂ t (ρ i+ 1 2 u) + ρu i+1 u up i+1 -ρu i u up i ∆x + p i+1 -p i ∆x = ∂ t (ρu)(x i ) + ∂ x (ρu 2 )(x i ) + ∂ x p(x i ) - ∆x 2 {sign(u(x i ))∂ x (u∂ x q)(x i ) -∂ xx p(x i ) -∂ xx (ρu 2 )(x i ) -∂ t (∂ x (ρu))(x i ) + O(∆x 2 ) ( 4
.42) From (4.40), we obtain:

∂ t (ρ i+ 1 2 u) + ρu i+1 u up i+1 -ρu i u up i ∆x + p i+1 -p i ∆x = ∂ t (ρu)(x i ) + ∂ x (ρu 2 )(x i ) + ∂ x p(x i ) - ∆x 2 sign(u(x i ))∂ x (u∂ x q)(x i ) + O(∆x 2 )
(4.43) Finally, the staggered scheme (4.16-4.17) have the following numerical diffusion operator when the consistency error is derived around x = x i in the momentum equation:

∂ t U ∆x + ∂ x F (U ∆x ) = D(U ∆x , ∆x) + O(∆x 2 ) D(U ∆x , ∆x) = ∆x 2 sign(u(x i )) 0 0 sign(u(x i )) ∂ x u(x i ) 0 0 u(x i ) ∂ x ρ(x i ) q(x i ) (4.44)
Hence, we have the following result on the numerical diffusion operator of the staggered schemes.

Theorem 4.2 (Numerical diffusion of staggered schemes). The second order perturbation operator associated to the staggered scheme (4.16-4.17) on a 1D regular mesh with space step ∆x is the strongly nonlinear operator :

D(U, ∆x) = x sign(u) 0 0 sign(u) ∂ x u 0 c 2 u ∂ x ρ q + o( x). (4.45)
The numerical diffusion associated to the mass conservation law is the same for every staggered scheme since they all have the same discretisation of the mass equation. Hence all the staggered schemes have a non classical diffusion operator at least for the mass equation. Each staggered scheme differs from the others with the discretisation of the momentum equation, especially the convective term (non conservative or conservative scheme, implicit or semi implicit, ...). Hence, theorem 4.2 holds for both the mass and the momentum equations only for the Herbin et al staggered scheme (4.16-4.17).

In the context of staggered schemes, there are two meshes, one for the density ρ and one for the velocity u. Hence we have two options to write the consistency error for the momentum equation. This leads to two numerical diffusion operators (4.36) and (4.44). (4.36) writes the momentum consistency equation at point x i+1/2 whereas (4.44) writes it at point x i . Provided that the velocity u has a constant sign then these operators are equivalent and the diffusion operator (4.36) obtained around x i+1/2 = x i + ∆x 2 converges towards the diffusion operator (4.44) obtained around x i , when ∆x tends to zero. However, if the velocity changes sign, there is not enough regularity to have the equivalence of (4.36) and (4.44). In practice, the velocity u changes sign on a finite number of points in the domain, hence, the numerical diffusion operators (4.36) and (4.44) will differ only in some points.

In the next section, we propose a new class of staggered schemes whose numerical diffusion operator has a classical form (4.3) and ensures the L 2 -stability of the scheme.

A new class of staggered schemes for the Euler equations

The most advanced result regarding staggered schemes are the kinetic energy inequality and the entropic stability proved for the scheme (4.16-4.17) in [START_REF]On some implicit and semi-implicit staggered schemes for the shallow water and euler equations[END_REF]. These properties however do not guarantee the boundedness of the solution (see the discussion in section 4.1.2). In this section, we propose a class of staggered schemes for conservation laws (4.1) which are linearly L 2 -stable. Unlike classical staggered schemes which have a non classical diffusion operator (see the discussion at the end of section 4.2.1 and therorem 4.2). We impose a classical diffusion operator D (4.3), such that the diffusion tensor verifies: D + t D ≥ 0, where D is the matrix D in the basis that symmetrises the Euler system. We specify this new class in the particular case of the following 2D isentropic Euler equations in conservative form:

     ∂ t ρ + ∇ • q = 0 ∂ t q + ∇ • q ⊗ q ρ + ∇p = 0 . ( 4.46)
This can be written:

               ∂ t ρ + ∂q x + ∂q y = 0 ∂ t q x + ∂ x q 2 x ρ + ∂ y q x q y ρ + c 2 ∂ x ρ = 0 ∂ t q y + ∂ x q 2 y ρ + ∂ x q x q y ρ + c 2 ∂ y ρ = 0 . (4.47)
The isentropic Euler system takes the conservative form

∂ t U + ∇ • F (U ) = 0, (4.48)
where U = (ρ, q) and F is the flux matrix. The isentropic Euler system takes the non conservative form

∂ t U + A x (U )∂ x U + A y (U )∂ y U = 0. (4.49)
It is usual to define the Jacobian of F along vectors n = (n x , n y ) ∈ R 2 as

A(U, n) = n x A x (U ) + n y A y (U ) = 0 t n c 2 n -( u • n) u u ⊗ n + ( u • n)I 2 , (4.50)
with c 2 = ∂p ∂ρ assumed constant.

We consider the class Stag of discrete staggered conservative schemes of the form: and: (4.51)

U i,j (t) + F x + -F x - ∆x + F y + -F y - ∆y = 0, with: U i,j =    ρ i,j q x i+1/2,j q y i,j+1/2    ,
F x + = F x (U i,j ) + F x (U i+1,j ) 2 + D Stag (U i,j , U i+1,j , n x ) U i,j -U i+1,j 2 , (4.52)
F y + = F y (U i,j ) + F y (U i,j+1 ) 2 + D Stag (U i,j , U i,j+1 , n y ) U i,j -U i,j+1 2
D(U, ∆x, ∆ y ) = ∆x∂ x (D Stag (U, U, n x )∂ x U ) + ∆y∂ y (D Stag (U, U, n y )(∂ y U ) + o( x, y).
on a regular mesh with space steps ∆x and ∆y.

In the sequel, we analyse an element of the class of schemes Stag whose expression is:

                       ∂ t ρ i,j + qx i+1/2,j -qx i-1/2,j ∆x + qy i,j+1/2 -qy i,j-1/2 ∆y = 0 ∂ t q x i+1/2,j + ¯ q 2 x ρ i+1,j - ¯ q 2 x ρ i,j ∆x + ¯ qxqy ρ i+1/2,j+1/2 - ¯ qxqy ρ i+1/2,j-1/2 ∆y + c 2 ρ i+1,j -ρ i,j ∆x = 0 ∂ t q y i,j+1/2 + ¯ q 2 y ρ i,j+1
F x + =         qx i+1/2,j ¯ q 2 x ρ i+1,j + c 2 ρ i+1,j ˜ q x q y ρ i+1/2,j+1/2         , F y + =          qy i,j+1/2 ¯ q x q y ρ i+1/2,j+1/2 ¯ q 2 y ρ i,j+1 + c 2 ρ i,j+1          . ( 4.55)
Using the Roe average, the scheme takes a more compact form that follows:

qx i+1/2,j = q x i+1/2,j + (|u x | -u x ) ρ i,j -ρ i+1,j 2 ¯ q 2 x ρ i+1,j = (q x i+1/2,j ) 2 ρ i,j + (|u x | -u x ) q x i+1/2,j -q x i+3/2,j 2 ˜ q x q y ρ i+1/2,j+1/2 = q x i+1/2,j q y i,j+1/2 ρ i,j + (|u y | -u y ) q x i+1/2,j -q x i+1/2,j+1 2
Where u x is a Roe-type average:

u x = q x i+1/2,j √ ρ i,j + q x i+3/2,j √ ρ i+1,j √ ρ i,j + √ ρ i+1,j qy i,j+1/2 = (q y i,j+1/2 ) 2 ρ i,j + (|u y | -u y ) q y i,j+1/2 -q y i,j+3/2 2 ¯ q x q y ρ i+1/2,j+1/2 = q x i+1/2,j q y i,j+1/2 ρ i,j + (|u x | -u x ) q y i,j+1/2 -q y i+1,j+1/2 2
Where u y is a Roe-type average:

u y = q y i,j+1/2 √ ρ i,j + q y i,j+3/2 √ ρ i,j+1 √ ρ i,j + √ ρ i,j+1
Using the fact that the numerical scheme (4.54) has a numerical flux of the form (4.51-4.53), we can determine the scheme upwind operator: D Stag (U i,j , U i+1,j , n x ).

F x + - F x (U i,j ) + F x (U i+1,j ) 2 =         qx i+1/2,j ¯ q 2 x ρ i+1,j + c 2 ρ i+1,j ˜ q x q y ρ i+1/2,j+1/2         - 1 2         q x i+1/2,j + q x i+3/2,j (q x i+1/2,j ) 2 ρ i,j + (q x i+3/2,j ) 2 ρ i+1,j + c 2 (ρ i,j + ρ i+1,j ) q x i+1/2,j ρ i,j q y i,j+1/2 + q x i+3/2,j ρ i+1,j q y i+1,j+1/2         = - 1 2   |u x | -u x 1 0 -c 2 -u 2 x |u x | + u x 0 -u x u y u y |u x |      ρ i+1,j -ρ i,j q x i+3/2,j -q x i+1/2,j q y i+1,j+1/2 -q y i,j+1/2    (4.56)
We obtain D Stag (U i,j , U i+1,j , n x ) and D Stag (U i,j , U i,j+1 , n y ) (with the same approach on F + y ) the matrices coming from the upwind contributions of the discrete equations:

D Stag (U i,j , U i+1,j , n x ) =    |u x | -u x 1 0 -c 2 -u 2 x |u x | + u x 0 -u x u y u y |u x |    , D Stag (U i,j , U i,j+1 , n y ) =   |u y | -u y 0 1 -u x u y |u y | u x -c 2 -u 2 y 0 |u y | + u y   (4.57)

The linearised system

In order to simplify the stability analysis we consider the linearised Euler system around a state with density ρ 0 , momentum q 0 and velocity u 0 = 1 ρ 0 q 0 . From the identity:

∇ • q ⊗ q ρ = (∇ • q) q ρ + ( ∇ q) q ρ - q ⊗ q ρ 2 ∇ρ,
the linearised Euler system is obtained as the following constant coefficient PDE system

   ∂ t ρ + ∇ • q = 0 ∂ t q + (∇ • q)I d + ∇ q u 0 + (c 2 I d -u 0 ⊗ u 0 ) ∇ρ = 0 . (4.58)
This can be written:

       ∂ t ρ + ∂ x q x + ∂ y q y = 0 ∂ t q x + u x 0 (∇ • q + ∂ x q x ) + u y 0 ∂ y q x + (c 2 -(u x 0 ) 2 )∂ x ρ -u x 0 u y 0 ∂ y ρ = 0 ∂ t q y + u y 0 (∇ • q + ∂ y q y ) + u x 0 ∂ x q y + (c 2 -(u y 0 ) 2 )∂ y ρ -u x 0 u y 0 ∂ x ρ = 0 . ( 4

.59)

The linearised Euler system takes the form

∂ t U + ∇ • F (U ) = 0, (4.60)
where U = (ρ, q) and F is the linearisation of the matrix flux F , and satisfies F (U ) n = Ā( n)U where the jacobian matrix Ā( n) = n x Āx + n y Āy has expression

∀ n ∈ R 2 , Ā( u 0 , n) = 0 t n c 2 n -( u 0 • n) u 0 u 0 ⊗ n + ( u 0 • n)I d , (4.61)
The scheme (4.54) applied to the linearised Euler system writes:

                                 ∂ t ρ i,j + qx i+1/2,j -qx i-1/2,j ∆x + qy i,j+1/2 -qy i,j-1/2 ∆y = 0 ∂ t q x i+1/2,j + 2u x 0 qx i+1,j -qx i,j ∆x + u x 0 qy i,j+1 -qy i,j ∆y + u y 0 qx i+1/2,j+1/2 -qx i+1/2,j-1/2 ∆y +c 2 ρ i+1,j -ρ i,j ∆x -(u x 0 ) 2 ρ i,j -ρ i-1,j ∆x -u x 0 u y 0 ρ i,j -ρ i,j-1 ∆y = 0 ∂ t q y i,j+1/2 + 2u y 0 qy i,j+1 -qy i,j ∆y + u y 0 qx i+1,j -qx i,j ∆x + u x 0 qy i+1/2,j+1/2 -qy i-1/2,j+1/2 ∆x +c 2 ρ i,j+1 -ρ i,j ∆y -(u y 0 ) 2 ρ i,j -ρ i,j-1 ∆y -u x 0 u y 0 ρ i,j -ρ i-1,j ∆x = 0 . (4.62)
The numerical scheme (4.62) can be written in the form (4.51)-(4.53) as follows: and: (4.63)

U i,j (t) + F x + -F x - ∆x + F y + -F y - ∆y = 0, with: U i,j =    ρ i,j q x i+1/2,j q y i,j+1/2    ,
F x + = F x (U i,j ) + F x (U i+1,j ) 2 + DStag ( u 0 , u 0 , n x ) U i,j -U i+1,j 2 , (4.64)
F y + = F y (U i,j ) + F y (U i,j+1 ) 2 + DStag ( u 0 , u 0 , n y ) U i,j -U i,j+1 2 , (4.65)
based on the following numerical flux:

F x + =      qx i+1/2,j 2u x 0 qx i+1,j + c 2 ρ i+1,j -(u x 0 ) 2 ρ i,j u y 0 qi+1,j -u x 0 u y 0 ρ i,j + u x 0 qy i+1/2,j+1/2      , F y + =     qy i,j+1/2 u x 0 qy i,j+1 -u x 0 u y 0 ρ i,j + u y 0 qx i+1/2,j+1/2 2u y 0 qy i,j+1 + c 2 ρ i,j+1 -(u y 0 ) 2 ρ i,j     .
(4.66) Using the Roe average, the scheme takes a more compact form that follows:

qx i+1/2,j = q x i+1/2,j + (|u x 0 | -u x 0 ) ρ i,j -ρ i+1,j 2 qx i+1,j = q x i+1/2,j + q x i+3/2,j 2 + 1 + sign(u x 0 ) 2 q x i+1/2,j -q x i+3/2,j 2 qx i+1,j = q x i+1/2,j qy i+1/2,j+1/2 = q y i,j+1/2 + q y i+1,j+1/2 2 + sign(u x 0 ) q y i,j+1/2 -q y i+1,j+1/2 2 qy i,j+1/2 = q y i,j+1/2 + (|u y 0 | -u y 0 ) ρ i,j -ρ i,j+1 2
qy i,j+1 = q y i,j+1/2 + q y i,j+3/2 2 + 1 + sign(u y 0 ) 2 q y i,j+1/2 -q y i,j+3/2 2 qy i,j+1 = q y i,j+1/2 qx i+1/2,j+1/2 = q x i+1/2,j + q x i+1/2,j+1 2 + sign(u y 0 ) q x i+1/2,j -q x i+1/2,j+1 2
Using the fact that the numerical scheme (4.62) has a numerical flux of the form (4.63-4.65), we can determine the scheme upwind operator: DStag ( u 0 , u 0 , n x ).

F x + - F x (U i,j ) + F x (U i+1,j ) 2 =      qx i+1/2,j 2u x 0 qx i+1,j + c 2 ρ i+1,j -(u x 0 ) 2 ρ i,j u y 0 qi+1,j -u x 0 u y 0 ρ i,j + u x 0 qy i+1/2,j+1/2      - 1 2      q x i+1/2,j + q y i+3/2,j 2u x 0 (q x i+1/2,j + q x i+3/2,j ) + (c 2 -(u x 0 ) 2 )(ρ i,j + ρ i+1,j ) u y 0 (q x i+1/2,j + q x i+3/2,j ) -u x 0 u y 0 (ρ i,j + ρ i+1,j ) + u x 0 (q y i,j+1/2 + q y i+1,j+1/2 )      = - 1 2   |u x 0 | -u x 0 1 0 -c 2 -(u x 0 ) 2 |u x 0 | + u x 0 0 -u x 0 u y 0 u y 0 |u x 0 |      ρ i+1,j -ρ i,j q x i+3/2,j -q x i+1/2,j q y i+1,j+1/2 -q y i,j+1/2    (4.67)
We have:

(V i,j -V i+1,j ) • DStag (V 0 , V 0 , n x )(V i,j -V i+1,j ) = 1 2 (V i,j -V i+1,j )• DStag (V 0 , V 0 , n x ) + t DStag (V 0 , V 0 , n x ) (V i,j -V i+1,j ) 1 2 d||V || 2 2 dt = - ∆y 4 i j (V i,j -V i+1,j ) • DStag (V 0 , V 0 , n x ) + t DStag (V 0 , V 0 , n x ) (V i,j -V i+1,j ) - ∆x 4 i j (V i,j -V i,j+1 ) • DStag (V 0 , V 0 , n y ) + t DStag (V 0 , V 0 , n y ) (V i,j -V i,j+1 ) = - ∆y 4 DStag (V 0 , V 0 , n x ) + t DStag (V 0 , V 0 , n x ) i j ||V i,j -V i+1,j || 2 2 - ∆x 4 DStag (V 0 , V 0 , n y ) + t DStag (V 0 , V 0 , n y ) i j ||V i,j -V i,j+1 || 2 2 Since: DStag (V 0 , V 0 , n) + t DStag (V 0 , V 0 , n) ≥ 0, we obtain 1 2 d||V || 2 2 dt ≤ 0.
Corollary 4.5. The numerical scheme (4.54) is linearly L 2 stable.

Proof. We determine the upwinding matrices DStag (V 0 , V 0 , n x ) and DStag (V 0 , V 0 , n y ):

We have: V = P U with:

P = c 0 -u 0 I 2 , P -1 =    1 c 0 u 0 c I 2   .
Thus, we have:

DStag (V 0 , V 0 , n x ) = P DStag (V 0 , V 0 , n x )P -1 = | u 0 • n x | c t n x -c n x | u 0 • n x |I 2 DStag (V 0 , V 0 , n y ) = P DStag (V 0 , V 0 , n y )P -1 = | u 0 • n y | c t n y -c n y | u 0 • n y |I 2
We have:

1 2 d||V || 2 2 dt = - ∆y 4 i j (V i,j -V i+1,j ) • DStag (V 0 , V 0 , n x ) + t DStag (V 0 , V 0 , n x ) (V i,j -V i+1,j ) - ∆x 4 i j (V i,j -V i,j+1 ) • DStag (V 0 , V 0 , n y ) + t DStag (V 0 , V 0 , n y ) (V i,j -V i,j+1 ) Since: DStag (V 0 , V 0 , n x )+ t DStag (V 0 , V 0 , n x ) = | u 0 • n x |I 3 and: DStag (V 0 , V 0 , n y )+ t DStag (V 0 , V 0 , n y ) = | u 0 • n y |I 3 1 2 d||V || 2 2 dt = - ∆y 4 | u 0 • n x | i j ||V i,j -V i+1,j || 2 2 - ∆x 4 | u 0 • n y | i j ||V i,j -V i,j+1 || 2 2 1 2 d||V || 2 2 dt ≤ 0
Since DStag satisfies the hypotheses of Theorem 4.4, this numerical scheme is linearly L 2 stable.

D Stag = 0 1 -c 2 0 + O( 1 c ),
obtained for a fixed velocity u and the sound speed c tending to infinity in (4.57). We see that the diffusion is less important on the mass than on the momentum equation while it is evenly distributed for the Roe scheme (see 4.13). In the mass equation, the discretisation introduced a perturbation

x∂ xx q: ∂ t ρ + ∂ x q = x∂ xx q + o( x). ( 4

.69)

It can be expected that the scheme (4.54) yield a better accuracy with less numerical diffusion on the mass equation than the upwind scheme (4.14).

Now that we proved that the scheme (4.54) is L 2 -stable, we show in the next section some numerical results with an implicit version of the numerical scheme.

Numerical results

In this section, we assess the behaviour of our new staggered scheme ( [START_REF] Ndjinga | A new class of L 2 -stable schemes for the isentropic Euler equations on staggered grids[END_REF]) on a one dimensional Riemann problem. The robustness of the scheme is illustrated on a compressible fluid with isothermal equation of state p = ρc 2 where the sound speed is c = 300m/s. We choose initial conditions such that the structure of the solution consists in a rarefaction wave followed by a shock wave, with sufficiently strong shock to allow an easy discrimination of correct numerical solutions. These initial conditions are: left state: ρ lef t q lef t = 10 9 100ρ lef t , right state:

ρ right q right = ρ lef t 2 -100ρ lef t
We consider this Riemann problem for the isentropic Euler system (1.11). The problem is posed over Ω = (0, 1) and the discontinuity is initially located at x = 0.5. In figures 4.1 and 4.2, the solution displays a rarefaction (smooth) wave followed by a (discontinuous) shock wave. Our new method is able to capture both waves in a distinct and stable way.

Conclusion

In this chapter, we developed a rigorous framework for the L 2 -stability analysis of finite volume schemes on staggered grids. The derivation of the method has required the analysis of some theoretical aspects beforehand. We also presented some analytical numerical examples on the solution of the isentropic Euler equations with the purpose of illustrating the technique and its performances. The family of schemes presented here could be applied in the future to the Cathare code to handle the numerical difficulties specific to two-phase flows models. A major challenge for the simulation of two-phase flows is the configuration of the vanishing phase where one of the phases disappears in some parts of the domain. The prediction of this complex dynamic mainly relies on the capture of the void waves that appear in the two-fluid model. Since the void waves have a complex structure with a propagation speed that frequently change signs, it is important for the numerical scheme to be stable regardless of the velocity sign. Our new class of staggered schemes is a promising alternative to the actual numerical treatment of the vanishing phase implemented in the Cathare code. The actual strategy relies on an interfacial friction coefficient that becomes high when one of the phases vanishes and has reached its limits for some test cases like the one of the Water-packing, that is relevant for nuclear safety studies. A first step toward the implementation of a new staggered scheme in the Cathare code could be the application to the reduced system of [START_REF] Ndjinga | A 2 × 2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF] that focuses on the study of the void waves.

ρ up i-1 2 u i-1 2 = ρ i-1 + ρ i 2 + sign(u i-1 2 ) ρ i-1 -ρ i 2 u i-1 2 = ρ(x i ) + 1 2 -∆x∂ x ρ(x i ) + 1 2 (∆x) 2 ∂ xx ρ(x i ) + O(∆x 3 ) (1 + sign(u(x i ))) × u(x i ) - ∆x 2 ∂ x u(x i ) + 1 2 ∆x 2 2 ∂ xx u(x i ) + O(∆x 3 ) = ρ(x i )u(x i ) - ∆x 2 (ρ(x i )∂ x u(x i ) + (1 + sign(u(x i )))u(x i )∂ x ρ(x i )) + ∆x 2 4 ρ(x i ) 1 2 ∂ xx u(x i ) + (1 + sign(u(x i ))) (u(x i )∂ xx ρ(x i ) + (∂ x ρ)(x i )(∂ x u)(x i ))
+O(∆x 3 ) = ρ(x i )u(x i ) -∆x 2 (ρ(x i )∂ x u(x i ) + (1 + sign(u(x i )))u(x i )∂ x ρ(x i ))

Introduction

In the previous chapter we have derived a linearly L 2 stable class of staggered schemes for the isentropic Euler equations. This linear stability has been characterised by using the properties of the numerical diffusion matrix. A L 2 stable scheme can however capture non entropic weak solutions as is the case with the Roe scheme ( [START_REF] Desprès | Systèmes hyperboliques de lois de conservation[END_REF][START_REF] Roe | Approximate riemann solvers, parameter vectors and difference schemes[END_REF]).

In this section, we focus on some non linear property of the scheme, namely the entropy property. We will characterise the staggered schemes that are entropic by analysing the properties of the diffusion matrix. In section 5.1.1, we recall the derivation of the entropy-entropy flux pair for the isentropic Euler equations. The main result of this chapter is in section 5.2 with the definition of a new class of entropic staggered schemes. For the sake of simplicity, the analysis is performed in the one dimensional case with a linear state equation p = ρc 2 but the strategy can be extended to the multidimensional case with a general state law p(ρ). In the last section, we illustrate the performance of a prototype scheme belonging to the class of entropic staggered schemes on a one dimensional Riemann problem.

Entropy of the isentropic Euler system

The derivation of the entropy-entropy flux pair is recalled in this section for the particular case p = ρc 2 . The proof that the family of staggered schemes is entropic will follow the same lines albeit using discrete quantities.

We recall that the mathematical entropy function in this case is

s = ρ 1 2
|| u|| 2 + c 2 ln ρ .

(5.1)

s is strictly convex since its hessian matrix

H(s) =     || q|| 2 ρ 3 + c 2 ρ - t q ρ 2 - q ρ 2 1 ρ     (5.2)
is positive definite since c 2 > 0.

Let us assume that u and ρ are smooth functions solving

∂ t ρ + ∇ • (ρ u) = 0 ρ∂ t u + ρ u • ∇ u + ∇p = 0, (5.3)
we compute (5.10)

∂ t s = 1 2 || u|| 2 + c 2 ln ρ ∂ t ρ + ρ u • ∂ t u + c 2 ρ ∂ t ρ = 1 2 || u|| 2 + c 2 (1 + ln ρ) ∂ t ρ + ρ u • ∂ t u.

Entropy bound for a new class of staggered schemes

For the sake of pedagogy, we investigate the state law p = ρc 2 . For general 1D flows, we study the following scheme:

∂ t ρ i + qi+ 1 2
qi-1 2 x = 0, (5.11)

∂ t q i+ 1 2 + 1 ∆x q2 i+1 ρ i+1 - q2 i ρ i + p i+1 -p i x = 0, (5.12)
In the sequel, we analyse the entropic character of this class of staggered schemes (5.11)-(5.12) that can be written in a more compact form as follows:

U i (t) + F i+1 -F i ∆x = 0, with: U i = ρ i q i+1/2
, and:

(5.13)

F i+1 = F (U i ) + F (U i+1 ) 2 + D Stag (U i , U i+1 ) U i -U i+1 2 (5.14)
where D Stag is a 2 × 2 matrix valued function:

D Stag (U i , U i+1 ) = a i,i+1 b i,i+1 c i,i+1 d i,i+1
Hence, the interpolated quantities qi+1/2 , qi-1/2 , q2 i+1 ρ i+1 and q2 i ρ i have the following expressions: qi+1/2 = q i+1/2 + q i+3/2 2 + a i,i+1 ρ iρ i+1 2 + b i,i+1 q i+1/2q i+3/2 2 = q i+1/2 + a i,i+1 ρ iρ i+1 2 + (b i,i+1 -1) q i+1/2q i+3/2 2 q2 i+1

ρ i+1 = q 2 i+1/2 ρ i + q 2 i+3/2 ρ i+1 2 + c i,i+1 ρ i -ρ i+1 2 + d i,i+1 q i+1/2 -q i+3/2 2 = q 2 i+1/2 ρ i - 1 2
q 2 i+1/2 ρ i - q 2 i+3/2 ρ i+1 + c i,i+1 ρ i -ρ i+1 2 + d i,i+1 q i+1/2 -q i+3/2 2 = q 2 i+1/2 ρ i + (c i,i+1 + u 2 i,i+1 ) ρ i -ρ i+1 2 + (d i,i+1 -2u i,i+1 ) q i+1/2 -q i+3/2 2
where the Roe average is given by:

u i,i+1 = q i+1/2 √ ρ i + q i+3/2 √ ρ i+1 √ ρ i + √ ρ i+1 (5.15)
The coefficients a i,i+1 , b i,i+1 , c i,i+1 and d i,i+1 depend on the left state U i = ρ i q i+1/2 and the right state U i+1 = ρ i+1 q i+3/2 . In the sequel, we will derive constraints on these coefficients to ensure the entropic character of a new class of staggered schemes.

From now on, we set the coefficient b i,i+1 = 1. This choice is motivated by the remarks 4. In this section, we discuss the entropic character of the L 2 -stable staggered scheme introduced in Chapter 4, (equation 4.54) in the one dimensional case. This staggered scheme has the following upwinding matrix:

D Stag (U i , U i+1 ) = |u i,i+1 | -u i,i+1 1
-c 2 -u 2 i,i+1 |u i,i+1 | + u i,i+1
In practice, the lack of entropy decrease results in a shock wave being captured instead of a rarefaction wave. This is particularly the case for Riemann solvers ( [START_REF] Desprès | Systèmes hyperboliques de lois de conservation[END_REF][START_REF] Roe | Approximate riemann solvers, parameter vectors and difference schemes[END_REF]) where the violation of entropy decrease leads to the capture of a non entropic stationnary shock wave instead of a transonic rarefaction wave. Since the upwinding matrix of our class of staggered schemes is not based on the eigenvalues of the system, it is not straightforward to choose a numerical test that will show the entropy default. Hence, we need the entropy balance (5.33) to identify the configurations that could lead to the capture of non entropic solutions.

For the L 2 -stable staggered scheme (4.54), we have:

∂ t s i + 1 ∆x g q i+1/2 ρ i , ρ i -g q i-1/2 ρ i-1 , ρ i-1 = - 1 ∆x E(a i-1,i , c i-1,i , d i-1,i )
with:

E(a i-1,i , c i-1,i , d i-1,i ) = c 2 q i-1/2 ρ i ρ i-1 -1 -ln ρ i ρ i-1 + |u i-1,i | -u i-1,i 2 (ρ i-1 -ρ i ) ln ρ i-1 ρ i + q i+1/2 ρ i - q i-1/2 ρ i-1 2 1 2 q i-1/2 + |u i-1,i | -u i-1,i 4
(ρ i-1 + ρ i )

(5.33) From (5.33), we see that the scheme violates the entropy decrease property in the two following cases:

For constant densities, if q i-1/2 < 0 and u i-1,i = 0, then E(a i-1,i , c i-1,i , d i-1,i ) < 0. Hence, in this case, the L 2 -stable scheme locally generates a positive contribution of the entropy in the order of |q i-1/2 | 2 .

For variable densities and high negative velocities u i-1/2 , then E(a i-1,i , c i-1,i , d i-1,i ) < 0. Hence, in this case, the L 2 -stable scheme locally generates a positive contribution of the entropy in the order of c 2 .

The largest amount of entropy creation occurs in the second case, because the entropy creation is proportional to c 2 . In the sequel, we solve Riemann problems for the isentropic Euler system with initial conditions such that the scheme initially violates the entropy decrease.

Entropy default for negative high velocities

In the sequel, we assess the behaviour of the L 2 -stable staggered scheme (4.54) on a one dimensional Riemann problem. We choose initial conditions such that the structure of the solution consists in a shock followed by a rarefaction wave to allow an easy discrimination of correct numerical solutions. These initial conditions are: left state: ρ lef t q lef t = 1 ρ lef t u 0 , right state: ρ right q right = 2 ρ right u 0 (5.34) We consider this Riemann problem for the isentropic Euler system (1.11). The problem is posed over Ω = (0, 1) and the discontinuity is initially located at x = 0.5. This case illustrates the second configuration of entropy default for the L 2 -stable scheme with variable densities and an entropy creation proportional to c 2 . In figures 5.1 and 5.2, the L 2 -stable staggered scheme captures a shock wave followed by a rarefaction wave. Hence, the method captures the correct entropic solution despite the initial entropy creation proportional to c 2 . Hence for long time simulations the entropy default we observe in theory in (5.33) does not lead to the capture of non entropic solutions. In this section, the numerical tests we choose do not generate an entropy violation in the entropy balance (5.33) of the L 2 -stable staggered scheme. They illustrate the behaviour of the L 2 -stable scheme for high and low Mach number flows. Here, we solve the Riemann problem with the initial condition (5.34) where u 0 = 300 (see figures 5.3 and 5.4) and u 0 = -1 (see figures 5.5 and 5.6). In these two cases, our conservative staggered scheme is able to capture shock and rarefaction waves without any prior information on the characteristic fields. In this section, we consider the Riemann problem where the initial conditions are chosen such that the structure of the solution consists in two transonic rarefaction waves: left state: ρ lef t q lef t = 1 -300

, right state: ρ right q right = 1 300 (5.35) This case illustrates the first configuration of entropy default for the L 2 -stable scheme with constant densities and an entropy creation proportional to |q i-1/2 | 2 (see introduction of section 5.1.1). In figures 5.7 and 5.8, we see that despite the initial entropy creation the scheme captures the correct entropic solution of the Riemann problem with two transonic rarefaction waves. Hence for long time simulations the entropy default we observe in theory in (5.33) does not lead to the capture of non entropic solutions. In this section, we show some numerical results obtained with the entropic staggered scheme defined in Corollary 5.2. This staggered scheme has the following numerical diffusion matrix:

D Stag (U i , U i+1 ) = a i,i+1 1
-c 2u 2 i,i+1 + 2u i,i+1 ∆v i,i+1 a i,i+1 + 2u i,i+1

with ∆v i,i+1 = q i+1/2 ρ i -q i+3/2 ρ i+1 . The coefficient a i,i+1 of this diffusion matrix is chosen such that the entropic character of this numerical method is ensured. For this numerical scheme, we have:

∂ t s i + 1 ∆x g q i+1/2 ρ i , ρ i -g q i-1/2 ρ i-1 , ρ i-1 = - 1 ∆x E(a i-1,i , c i-1,i , d i-1,i )
with:

E(a i-1,i , c i-1,i , d i-1,i ) = c 2 q i-1/2 ρ i ρ i-1 -1 -ln ρ i ρ i-1 + a i-1,i 2 (ρ i-1 -ρ i ) ln ρ i-1 ρ i + q i+1/2 ρ i - q i-1/2 ρ i-1 2 1 2 q i-1/2 + a i-1,i 4 (ρ i-1 + ρ i ) + u i-1,i (ρ i-1 -ρ i )
(5.36) The coefficient a i-1,i is set to ensure E(a i-1,i , c i-1,i , d i-1,i ) ≥ 0 and hence the entropic character of this numerical method, according to the semi-discrete analysis of the entropy balance in section 5.2.

Case of negative high velocities

In the sequel, we assess the behaviour of the entropic staggered scheme (5.2) on a one dimensional Riemann problem. We choose initial conditions such that the structure of the solution consists in a shock followed by a rarefaction wave to allow an easy discrimination of correct numerical solutions. This initial condition is given by (5.34) with u 0 = -300. In figures 5.9 and 5.10, the numerical scheme (5.2) captures a shock wave followed by a rarefaction wave with spurious oscillations. Hence, the method captures the correct entropic solution of the Riemann problem but generates oscillations around the rarefaction wave. The entropic staggered scheme seems to have a lower numerical diffsion than the L 2 -stable staggered scheme. The semi-discrete analysis gives a theoretical lower bound for the coefficient a i-1,i to ensure the entropic character of the scheme but seems not sufficient to avoid spurious oscillations. In this section, we illustrate the behaviour of the entropic staggered (5.2) scheme for high and low Mach number flows. Here, we solve the Riemann problem with the initial condition (5.34) where u 0 = 300 (see figures 5.11 and 5.12) and u 0 = -1 (see figures 5.13 and 5.14). In these two cases, our conservative entropic staggered scheme (5.2) is able to capture shock and rarefaction waves without any prior information on the characteristic fields. In the case of high Mach number flow, the numerical method captures the correct entropic solution of the Riemann problem but generates oscillations around the stationnary shock wave. In the case of low Mach number flows, we do not observe spurious oscillations.

Conclusion

In this chapter, we developed a rigorous framework for the analysis of the entropic character of finite volume schemes on staggered grids. We also presented some analytical numerical examples on the solution of the isentropic Euler equations. These numerical results illustrate the behaviour of the L 2 -stable staggered scheme from Corollary 4.5 in configurations with a theoretical entropy violation, on the one hand. They also illustrate the behaviour of the entropic staggered scheme from Corollary (5.2) in the same configurations, on the other hand. The L 2 -stable staggered scheme captures the correct entropic solution despite the initial entropy default. In future works, we will analyse in detail the local entropy violation in the first time steps for the L 2 -stable scheme to have a better insight on the entropic character of this scheme. Depending on the conclusions from this first analysis, we will derive a fully discrete entropic analysis. On the other hand, the entropic staggered scheme from Corollary 5.2 captures the correct entropic solution but generates spurious oscillations. In future works, we will derive a strategy to limit the appearance of these spurious oscillations for the entropic staggered scheme. The multidimensionnal formulation of the new class of entropic schemes, introduced in this chapter, and its entropic stability analysis will be in a forthcoming paper. An extension of this strategy to the six-equation two-fluid model will be the subject of future works.

Conclusions and perspectives

In the present work, we have first of all developed acceleration techniques in the Cathare code in order to deal with the computational complexity with reasonable computing times. These methods take advantage of modern computer architectures by the use of a time domain decomposition method. The latter has been implemented with the parareal in time algorithm. A very special stress has been put on the adaptation of this algorithm to the Cathare code in a non intrusive way allowing to use the Cathare code as a black box. This development will be useful for several applications. First of all, the tool is important for safety calculations in the nuclear industry for the analysis of two-phase flows during accidental scenarii. A major challenge for the Cathare code is to produce real-time simulations when the software is governing a reactor simulator. A reactor simulator allows to reproduce the behaviour of a nuclear power plant under nominal or accidental conditions for the training of the operators and for the validation of the emergency procedures. Coupling the parareal algorithm with the actual acceleration techniques of the Cathare code represents a step toward a real-time response of the code.

Furthermore, since the time discretisation of the two-fluid model is done through a two-step time scheme, we designed a new variant of the parareal algorithm adapted to this family of methods. This work aims at solving the loss of accuracy in the parareal algorithm that can arise when an initialisation error is made at each time window. The parallel efficiency of this new variant is similar to the one achieved by the classical parareal algorithm and a way to improve these performances is to introduce adaptivity in the algorithm ( [START_REF] Maday | An adaptive parareal algorithm[END_REF]), by dynamically increasing the accuracy of the fine solver accross the parareal iterations. This adaptive formulation of the parareal algorithm offers new degrees of freedom to optimise the speed-up performances such as the choice of increasing target tolerances. A very interesting and challenging task would be the design of adaptive refinements based on a posteriori estimators. It would allow local time stepping adaptation in the parareal algorithm as well as spatial refinement if the problem involves also spatial variables. This strategy would be particularly interesting in the context of hyperbolic equations with whom the parareal algorithm may suffer from an instability. It has been shown ( [START_REF] Ruprecht | Convergence of Parareal with spatial coarsening[END_REF], [START_REF] Ruprecht | Wave propagation characteristics of Parareal[END_REF], [START_REF] Lunet | Time-parallel simulation of the decay of homogeneous turbulence using Parareal with spatial coarsening[END_REF]) that the numerical diffusion of the coarse and the fine solvers has an impact on the appearance of the parareal instability. Hence, introducing local refinements in the algorithm may reduce the instability we may observe when the parareal algorithm is applied to hyperbolic problems.

The second main contribution of this work has been devoted to the understanding of the theoretical properties of finite volume schemes on staggered grids such as the one used in the Cathare code. The idea consists in analysing the properties of the numerical diffusion operator. After showing that the staggered schemes do not straightforwardly yield a linear stability, we derive a linearly L 2 -stable class of staggered schemes for the isentropic Euler equations. We also implemented a scheme of the L 2 -stable class of schemes and perform a simulation of a Riemann problem with satisfactory results. Unlike classical staggered schemes that are L 2 -stable for constant sign velocities, the new class is L 2 -stable for variable sign velocities. This property is important in the context of the approximation of two-phase flows models since the phasic velocities frequently change signs. A major challenge for the simulation of two-phase flows is the configuration of the vanishing phase where one of the phases disappears in some parts of the domain. The prediction of this complex dynamic mainly relies on the capture of the void waves that appear in the two-fluid model. Since the void waves have a complex structure with a propagation speed that frequently change signs, it is important for the numerical scheme to be stable regardless of the velocity sign. Our new class of staggered schemes is a promising alternative to the actual numerical treatment of the vanishing phase implemented in the Cathare code. The actual strategy relies on an interfacial friction coefficient that becomes high when one of the phases vanishes and has reached its limits for some test cases like the one of the Water-packing, that is relevant for nuclear safety studies. A first step toward the implementation of a new staggered scheme in the Cathare code could be the application to the reduced system of [START_REF] Ndjinga | A 2 × 2 hyperbolic system modelling incompressible two phase flows: theory and numerics[END_REF] that focuses on the study of the void waves.

In a second time, we analyse a non linear property of the staggered schemes, namely the entropy property. We design a new class of entropic staggered schemes for the isentropic Euler equations by deriving conditions on the coefficients of the numerical diffusion operator. We implemented a scheme of this class of entropic schemes and perform with success a simulation of a Riemann problem that displays a transonic rarefaction wave with successful results. Hence, our conservative entropic staggered scheme is able to capture shock and rarefaction waves without any prior information on the characteristic fields. In the continuation of this analysis, there are several points that still need to be addressed for a better theoretical understanding of staggered schemes like the low Mach number accuracy. The multidimensionnal formulation of the new class of entropic schemes and its entropic stability analysis will be in a forthcoming work. An extension of this strategy to the six-equation two-fluid model will be the subject of future works.

1. 1

 1 Oscillating manometer test case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Convergence of the 1D Cathare scheme . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Setting of the Water-Packing benchmark in [68] . . . . . . . . . . . . . . . . . . . . 1.4 Behaviour of the liquid velocity during the Water-Packing phenomenon, from [86] . 1.5 Pressures spikes with or without Cathare anti-Water-Packing correction, from [86] 1.6 Example of a geometry in the Cathare code with elements and junctions . . . . . . 1.7 Performances with the actual parallelism in Cathare for two industrial test cases on 12 processors, [101] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 Load balancing between two threads in a Cathare simulation . . . . . . . . . . . .2.1 Correction of u1,N f -1 k+1 at time T 2δt in [T 1 , T 2 ] for the initialisation of the fine propagation in [T 2 , T 3 ] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Correction of u 1,N f -R k+1 at time T 2 -δT in [T 1 , T 2 ]for the initialisation of the coarse propagation in [T 2 , T 3 ] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Oscillating manometer test case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Convergence of the multi-step parareal algorithm when δt = 10 -5 and ∆T = 10δt . 2.5 Strong scaling results with the multi-step variant of the parareal algorithm . . . . 2.6 Sending only the principal variables to the Cathare code . . . . . . . . . . . . . . . 2.7 Sending all the state to the Cathare code on 10 time windows . . . . . . . . . . . . 2.8 Sending only the principal variables to the Cathare code on 5 time windows . . . . 2.9 Sending all the state to the Cathare code on 5 time windows . . . . . . . . . . . . 2.10 Convergence of the multi-step parareal algorithm when δt = 10 -4 and ∆T = 0.5 for an industrial test case on 5 processors . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 The reference solution and the fine solution after 1 parareal iteration . . . . . . . . 2.12 The reference solution and the fine solution after 2 parareal iterations . . . . . . . 2.13 The reference solution and the fine solution after 3 parareal iterations . . . . . . . 3.1 Convergence of the multi-step parareal for the second-order BDF method, δt = 10 -4 (left: error (3.34), right: error (3.35)) . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Convergence of the multi-step parareal for the third-order BDF method, δt = 10 -4 (left: error (3.34), right: error (3.35)) . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Convergence of the multi-step parareal for the second-order BDF method, δt = 10 -4 (left: error (3.34), right: error (3.35)) . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Convergence of the multi-step parareal for the third-order BDF method, δt = 10 -4 (left: error (3.34), right: error (3.35)) . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Density at time t = 0.001 with ∆x = 0.002 and CFL= 0.99 . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Oscillating manometer test case

5 Figure 1 . 2 :

 512 Figure 1.2: Convergence of the 1D Cathare scheme

Figure 1 . 4 :

 14 Figure 1.4: Behaviour of the liquid velocity during the Water-Packing phenomenon, from[START_REF] Mohaupt | Étude du problème numérique de Water-Packing[END_REF]

Figure 1 . 6 :

 16 Figure 1.6: Example of a geometry in the Cathare code with elements and junctions

Figure 1 . 7 :

 17 Figure 1.7: Performances with the actual parallelism in Cathare for two industrial test cases on 12 processors, [101]

Figure 1 . 8 :

 18 Figure 1.8: Load balancing between two threads in a Cathare simulation

  Parareal library for the Cathare code . . . . . . . . . . . . . . . . . . . 48 2.1.1 Obstructions linked to the data structure of the Cathare code . . . . . . . . . . 48 2.1.2 Obstructions linked to the time discretisation . . . . . . . . . . . . . . . . . . . 50 2.1.3 A numerical clone of the Cathare code . . . . . . . . . . . . . . . . . . . . . . . 50 2.2 Multi-step variant of the parareal algorithm . . . . . . . . . . . . . . . . . . 51 2.2.1 Original parareal algorithm and notations . . . . . . . . . . . . . . . . . . . . . 52 2.2.2 Adaptation to multi-step time schemes . . . . . . . . . . . . . . . . . . . . . . . 52 2.3 Application to the Cathare code . . . . . . . . . . . . . . . . . . . . . . . . . 55 2.3.1 The oscillating manometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 2.3.2 An industrial test case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

2. 1 . 1

 11 Obstructions linked to the data structure of the Cathare code 2.1.1.1 Data structure

Figure 2 . 3 :

 23 Figure 2.3: Oscillating manometer test case

Figure 2 . 4 :

 24 Figure 2.4: Convergence of the multi-step parareal algorithm when δt = 10 -5 and ∆T = 10δt

Figure 2 . 5 :

 25 Figure 2.5: Strong scaling results with the multi-step variant of the parareal algorithm

Figure 2 . 6 :

 26 Figure 2.6: Sending only the principal variables to the Cathare code

Figure 2 . 7 :

 27 Figure 2.7: Sending all the state to the Cathare code on 10 time windows

Figure 2 . 8 :Figure 2 . 9 :

 2829 Figure 2.8: Sending only the principal variables to the Cathare code on 5 time windows

Figure 2 . 10 :

 210 Figure 2.10: Convergence of the multi-step parareal algorithm when δt = 10 -4 and ∆T = 0.5 for an industrial test case on 5 processors

Figure 2 . 11 :

 211 Figure 2.11: The reference solution and the fine solution after 1 parareal iteration

Figure 2 . 12 :

 212 Figure 2.12: The reference solution and the fine solution after 2 parareal iterations

Figure 2 .

 2 Figure 2.13: The reference solution and the fine solution after 3 parareal iterations

Contents 3 . 1 70 3. 3 79 3. 4 80 3. 4 . 2

 317037948042 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.2 A multi-step variant of the parareal algorithm . . . . . . . . . . . . . . . . . 69 3.2.1 Setting and preliminary notations . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.2.2 A multi-step variant of the parareal algorithm . . . . . . . . . . . . . . . . . . . Advantages of the multi-step parareal algorithm . . . . . . . . . . . . . . . . Numerical tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 3.4.1 Numerical convergence results . . . . . . . . . . . . . . . . . . . . . . . . . . . . Parallel efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

21 ) 1 k

 211 Since the upper bound of the error term ||E n,N f -|| depends on the error terms ||E n+1 k || and ||δE n+1 k || we focus on the inequalities (3.19)-(3.21). Hence, we can write by induction:

Figure 3 . 1 :

 31 Figure 3.1: Convergence of the multi-step parareal for the second-order BDF method, δt = 10 -4 (left: error (3.34), right: error (3.35))

Figure 3 . 2 :

 32 Figure 3.2: Convergence of the multi-step parareal for the third-order BDF method, δt = 10 -4 (left: error (3.34), right: error (3.35))

Figure 3 . 3 :

 33 Figure 3.3: Convergence of the multi-step parareal for the second-order BDF method, δt = 10 -4 (left: error (3.34), right: error (3.35))

Figure 3 . 4 :

 34 Figure 3.4: Convergence of the multi-step parareal for the third-order BDF method, δt = 10 -4 (left: error (3.34), right: error (3.35))

Remark 4 . 6 .

 46 For one dimensionnal flows and low Mach numbers ( |u| c 1), the numerical scheme (4.54) of the class Stag has the following numerical diffusion matrix:

Figure 4 . 1 :

 41 Figure 4.1: Density at time t = 0.001 with ∆x = 0.002 and CFL= 0.99

Figure 4 . 2 :

 42 Figure 4.2: Momentum at time t = 0.001 with ∆x = 0.002 and CFL= 0.99

2 ∂ 2 ∂ 2 ∂ 2 ∂

 2222 xx ρ(x i ) + O(∆x 3 ) (1sign(u(x i ))) xx u(x i ) + O(∆x 3 ) = ρ(x i )u(x i ) + ∆x 2 (ρ(x i )∂ x u(x i ) + (1sign(u(x i )))u(x i )∂ x ρ(x i )) xx u(x i ) + (1sign(u(x i ))) (u(x i )∂ xx ρ(x i ) + (∂ x ρ)(x i )(∂ x u)(x i )) +O(∆x 3 ) = ρ(x i )u(x i ) + ∆x 2 (ρ(x i )∂ x u(x i ) + (1sign(u(x i )))u(x i )∂ x ρ(x i )) xx u(x i ) + (1sign(u(x i )))∂ x (u∂ x ρ)(x i ) +O(∆x 3 ) (A.1)

(5. 4 ) 2 || u|| 2 + c 2 ( 1 + 1 2( 5 . 5 )

 42221155 Using the fact that u and ρ solve (5.3), we obtain-∂ t s = 1 ln ρ) ∇ • (ρ u) + ρ u • ( u • ∇u) + u • ∇p = || u|| 2 ∇ • (ρ u) + c 2 ln ρ∇ • (ρ u) + c 2 ρ∇ • u + c 2 u • ∇ρ + ρ u • ( u • ∇u) + u • ∇p. Since p = ρc 2 and c 2 ρ∇ • u + u • ∇p = ∇ • (p u) we obtain -∂ t s = 1 2 || u|| 2 ∇ • (ρ u) + c 2 ln ρ∇ • (ρ u) + c 2 u • ∇ρ + ρ u • ( u • ∇u) + ∇ • (p u). (5.6) Since u • ∇u = 1 2 ∇|| u|| 2 + ( ∇ × u) × u we obtain -∂ t s = 1 2 || u|| 2 ∇ • (ρ u) + c 2 ln ρ∇ • (ρ u) + c 2 u • ∇ρ + 1 2 ρ u • ∇|| u|| 2 + ∇ • (p u) || u|| 2 ρ u) + c 2 ln ρ∇ • (ρ u) + c 2 u • ∇ρ + ∇ • (p u).

(5. 7 )

 7 Since ln ρ∇ • (ρ u) + u • ∇ρ = ∇ • (ρ ln ρ u) we finally obtain-∂ t s = ∇ • ( 1 2 || u|| 2 ρ u + c 2 ρ ln ρ u + p u).(5.8)and finally ∂ t s + ∇ • ((s + p) u) = 0, (5.9)and s is an entropy with the associated entropy flux g( u, ρ) = ρ 1 2 || u|| 2 + c 2 ln ρ + p u.

  1 and 4.6 made in the previous chapter about the low Mach number accuracy of the upwind-type schemes and 5.3 Numerical results 5.3.1 Entropy default of the L 2 -stable staggered scheme from Corollary 4.5

Figure 5 . 1 :Figure 5 . 2 : 5 . 3 . 1 . 2

 51525312 Figure 5.1: Density at time t = 5 × 10 -4 with ∆x = 1 × 10 -3 and CFL= 0.99 for u 0 = -300

Figure 5 . 3 :Figure 5 . 4 :Figure 5 . 5 : 1 Figure 5 . 6 : 1 5. 3 . 1 . 3

 5354551561313 Figure 5.3: Density at time t = 5 × 10 -4 with ∆x = 1 × 10 -3 and CFL= 0.99 for u 0 = 300

Figure 5 . 7 : 99 Figure 5 . 8 :

 579958 Figure 5.7: Density at time t = 5 × 10 -4 with ∆x = 1 × 10 -3 and CFL= 0.99

Figure 5 . 9 :Figure 5 . 10 :

 59510 Figure 5.9: Density at time t = 5 × 10 -4 with ∆x = 1 × 10 -3 and CFL= 0.99 for u 0 = -300

Figure 5 . 11 :Figure 5 . 13 : 1 Figure 5 . 14 :

 5115131514 Figure 5.11: Density at time t = 5 × 10 -4 with ∆x = 1 × 10 -3 and CFL= 0.99 for u 0 = 300

  3.3.1), simulations made with the Cathare scheme on test cases showing vanishing phase (see section 1.3.3.2,

strategies of discretisation for discontinuous source terms occuring in two-fluid models (see section 1.3.3.3).

  Stag is a 3 × 3 matrix valued function. An example of scheme in the class Stag is the staggered centered scheme, which correspond to the case D Stag = 0. Schemes of the class Stag admit a classical diffusion operator (4.3). (Classical diffusion of Stag schemes). Let D Stag : R 3 → R 3×3 be a matrix valued Lipshitz function. A staggered conservative scheme (4.51) with a numerical flux F x + and F x -of the form (4.53) admits the following classical diffusion operator

	,	(4.53)
	where D Theorem 4.3

1.2. Hyperbolicity of the model

1.3. Discretisation of the two-fluid model

Acknowledgements

Acknowledgements v Introduction (Version française) Introduction (English version)

Part II

Analysis of finite volume schemes on staggered grids

on the multi-step variant, we need to carefully initialise each time window. If the fine scheme is the second-order BDF method, the computation of u n+1 depends on u n and u n-1 and with the adaptive paradigm we have:

t nt n-1 = t n+1t n Hence, we initialise the fine solver with one variable step-size BDF method.

We apply this strategy to the damped oscillator system with the Backward Euler method as a coarse solver (∆T = 0.1) and the second-order BDF method as a fine solver with the sequence of time steps indicated in table 3.1. Table 3.1: Convergence of the adaptive parareal and the multi-step parareal with a target accuracy η = 3 × 10 -8

Multi-step parareal

The multi-step parareal algorithm with adaptivity converges to the exact solution with an accuracy obtained by a sequential fine solution with time step δt = 10 -4 after 8 iterations such as the multi-step method without adaptivity (see table 3.1). With the sequence of fine time steps used in the adaptive parareal method, convergence is reached with the same number of iterations as the multi-step variant. The adaptive algorithm allows to obtain better speed-up performances compared to the nonadaptive version since the fine solver (δt = 10 -4 ) is used only one time instead of 8 times in the multi-step variant. In table 3 adaptive and multi-step parareal algorithms applied to the damped oscillator. The speed-up is defined as the ratio:

between the cost to run a sequential fine solver achieving a target accuracy η with the cost to run a parareal algorithm providing at the end the same target accuracy η. The parallel efficiency of We obtain DStag ( u 0 , u 0 , n x ) and DStag ( u 0 , u 0 , n y ) (with the same approach on F + y ) the matrices coming from the upwind contributions of the discrete equations:

Linear stability of the class Stag

In order to study the stability of the Stag class of staggered schemes, the first step is to use the variables that symmetrise the continuous system and the second to prove that the energy (L 2 norm) of the new variables decreases with time. The reason why symmetrising is important is that the contribution of the first order derivatives to the energy balance:

which is the divergence of the vector field 1 2

The linearised Euler system can be symmetrised using the variable V = cρ qρ u 0 . We obtain the following system that is equivalent to the linearised Euler system:

The symmetrisation of the linearised Euler system therefore takes the form:

where F is the symmetrisation of the linearised matrix flux F . The Jacobian matrix associated to F is the symmetric operator Ā with expression:

From the symmetrised Euler system, we have the following property of L 2 stability for the schemes of the class Stag. Proof. After the linearisation around the state V 0 R 3 and symmetrisation, the Euler system takes the form:

The values V i,j of V in the cell C i,j are solutions of:

with fluxes:

We compute the evolution in time of

and the positiveness of DStag :

Then:

The numerical diffusion of the Herbin et al staggered scheme

Preliminary calculations for the mass equation From (4.18) and (4.21) we have

Preliminary calculations for the momentum equation

3) and (A.5), we have:

From (A.4) and (A.6), we have:

Chapter 5

A new class of entropic staggered schemes for the Euler equations Contents the staggered schemes. Setting b i,i+1 = 1, we obtain the following expressions for the interpolated quantities:

Hence, the class of staggered schemes we study in this section takes the form (5.13)- (5.14) where the matrix diffusion has the following generic expresssion:

where the term -c 2 comes from the discretisation of the pressure gradient.

Theorem 5.1 (A class of entropic staggered schemes). A staggered conservative scheme (5.13) with a numerical flux F i,i+1 of the form (5.14) such that, the coefficients a i-1,i , c i-1,i and d i-1,i of the diffusion operator D Stag satisfy:

(5.17) is entropic and the following discrete entropy dissipation estimate holds:

where:

(5.19) is the numerical entropy flux that is consistent with (5.10).

Proof. Firstly, we seek to derive the discrete analog of the velocity evolution equation:

From (5.12), we derive:

Assuming ρi+1/2 = ρ i , and from (5.11), we obtain:

From the expression of the interpolated quantities qi+1/2 and q2 i+1 ρ i+1 , this yields to:

||u|| 2 + c 2 ln ρ , we have:

From (5.21) and (5.12), we obtain:

), and:

We obtain:

In order to deal with the term (q i+1/2q i-1/2 ) ln ρ i , we add and substract the quantity

In order to absorb the terms containing the variables ρ i+1 and q i+3/2 into flux differences, we add and substract the quantities:

(5.23)

We obtain:

which becomes:

Corollary 5.2. The numerical scheme (5.13) with a numerical flux F i,i+1 of the form (5.14) such that the coefficients a i-1,i , c i-1,i and d i-1,i of the diffusion operator D Stag satisfy:

where u i-1,i is the Roe average (5.15) and:

is entropic and the discrete entropy dissipation estimate (5.18) holds. The numerical entropy flux is given by (5.19) and is consistent with (5.10).

Proof. Let us recall the entropy balance of our class of entropic staggered schemes:

(5.27) We define α i-1,i as:

Hence the coefficients a i-1,i and d i-1,i are defined as:

The entropy balance (5.27) becomes:

(5.28) From hypothesis (5.25), the coefficient c i-1,i is defined as:

and we finally obtain:

This term is positive provided the upwinding coefficient a i-1,i is large enough. The threshold values are

(ρ i-1ρ i ) ln ρ i-1 ρ i

(5.31) and

q i-1/2 + u i-1,i (ρ i-1ρ i ) (5.32)

Remark 5.3. Particular case: ρ iρ i-1 → 0:

Here, we give the asymptotic behaviour, when (ρ iρ i-1 ) tends to zero, of the coefficient a i-1,i that satisfies the condition (5.26). When ρ i ρ i-1 -1 → 0, we have:

, where C is a constant.

Hence, the first lower bound (5.31) of the coefficient a i-1,i in (5.26) is asymptotically bounded when (ρ iρ i-1 ) tends to zero.

Case of constant densities

In this section, we consider the Riemann problem where the initial conditions are chosen such that the structure of the solution consists in two transonic rarefaction waves: left state: ρ lef t q lef t = 1 -300

, right state: ρ right q right = 1 300

(5.37)

In figures 5.15 and 5.16, we see that the entropic staggered scheme (5.2) captures the correct entropic solution of the Riemann problem with two transonic rarefaction waves. In this case, we do not observe spurious oscillations.