Tout d'abord je tiens à remercier mes encadrants Edouard Audit, Pierre Kestener, Samuel Kokh et Pascal Tremblin. En seulement trois années vous m'avez permis d'apprendre énormément de choses sur des domaines très variés. Je remercie tout particulièrement Pascal Tremblin pour sa très grande disponibilité et son état d'esprit positif à toute épreuve. J'ai beaucoup appris grâce à toi et à travers les nombreuses discussions sur les schémas numériques, le calcul intensif et la physique. Je veillerai à appliquer du mieux possible le théorème d'Astérix. Je remercie Samuel Kokh pour son encadrement, notamment sur les modèles diphasiques et pour sa rigueur de relecture. Je remercie également Pierre Kestener pour avoir partagé ses connaissances sur les modèles de programmation parallèle et sur les bonnes pratiques du C++et d'avoir mis à disposition ses nombreuses maquettes "Euler2d". J'aimerais ensuite remercier mon jury de thèse. Merci à Ann Almgren et Philippe Helluy d'avoir accepté d'être rapporteurs de mon manuscrit, en espérant qu'ils aient apprécié la lecture. Je remercie également Isabelle Baraffe, Caroline Nore et Carlos Parès d'avoir participé au jury en tant qu'examinateurs. Je les remercie pour l'intérêt qu'ils ont porté à mon travail ainsi qu'aux échanges lors de la soutenance.

Je remercie chaleureusement les permanents et non-permanents de mon laboratoire d'accueil, la Maison de la Simulation. Les nombreuses discussions autour des pauses café ou déjeuner ont notablement participé au bon déroulement de la thèse. Je remercie mes anciens collègues thésards Arnaud, Hélène, Ksander, Maxime, Vitaly et tout particulièrement mon collègue de bureau Ziqiang pour toutes les discussions, presque philosophiques, autour des modèles diphasiques. Enfin je souhaite une bonne continuation à Rémi pour cette aventure qu'est la thèse.

Je remercie également mon comité de suivi de thèse Julien Bigot et Patrick Hennebelle pour leurs conseils et pour s'être assurés du bon déroulement de ma thèse.

Je remercie mes amis Nantais de longue date Anthony, Fadel et Tito pour tous ces bons moments et échanges qui ont

confinement.

Enfin je tiens à remercier toute ma famille, mes parents, mes frères ainsi que mes belles-soeurs pour m'avoir permis et encouragé à poursuivre mes études lorsque la motivation se faisait plus rare. v List of Figures Graphical representation of different heat transfer mecanisms from the core of a star to its surface. At first, the heat generated by nuclear reactions from the core is diffused by radiation. Then the heat transfer is achieved by the advection of the gas up to the surface of the star. Source https://solarscience.msfc.nasa. 2 Comparison of performance on different architectures: Intel KNL, Intel Skylake (one socket), NVIDIA K80, NVIDIA P100 and NVIDIA V100. Measures on Intel KNL and Intel Skylake were performed on Joliot-Curie's supercomputer at TGCC using the same code. In our case we obtain better results with the Intel Skylake than the Intel KNL due to a lack of vectorization. Going to a GP-GPU we have a speed-up around five with a NVIDIA K80 compared to multi-core architecture and seven between NVIDIA K80 and V100. . . . . . . . . . . . . . . . . . . . . . 1.3 Weak scaling results obtained on Joliot Curie's Intel Skylake partition at TGCC.

We use a hybrid MPI-OpenMP configuration in which one MPI task is bound to a socket. Simulations run for 1000 time steps and each MPI process treats 128 3 cells. We see that the efficiency reaches a plateau of 85%. . . . . . . . . . . . . . 1.4 Sod's test case simulations. Figure shows a snapshot of the density profile ρ for the All-Regime scheme, with and without the low-Mach correction, a first order Godunov-type scheme (HLLC) and the exact solution. Spatial resolution is nx = 100. We see that the All-Regime scheme gives results close to the Godunov-type scheme around discontinuities but is more diffusive in the rarefaction wave. . . . 1.5 Gresho vortex simulations. Snapshots of the magnitude of the velocity field at time t f = 10 -3 , for a resolution of 512 2 and for different Mach numbers. First line shows results where the low-Mach correction is disabled and second line where it is enabled. We see that without the low-Mach correction the scheme fails at simulating low-Mach flows. . two in yellow at time t = 12.4 and for a resolution of 200 × 600. First line show results with the the all-regime scheme, where on the left the low-Mach correction is disabled and is enabled on the right. Second line shows results with a Godunovtype scheme, on the left it is first order, on the right it is second order using a Muscl-Hancock scheme. We see that with the low-Mach correction we recover features only present at second order for a standard Godunov-type scheme. . . . 30 1.9 Rayleigh-Bénard instability simulations in 2D. Figure shows the time evolution of the mean absolute velocity for different ratios of Rayleigh number over critical Rayleigh number (see legend). Blue points show the case where the low-Mach correction is enabled and orange ones where it is disabled. We observe that when the low-Mach correction is enabled the onset of convection is closer to the expected critical Rayleigh number. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 1.10 Rayleigh-Bénard instability simulations. Snapshot of the local Mach number field and the velocity field. We see that in the strong stratification case, there is a large range of Mach, near zero at the center of rolls up to half at the upper boundary. 32 1.11 Rayleigh-Bénard instability simulations in 3D. Figure shows the kinetic energy spectrum of the horizontal middle plane. The blue line corresponds to the scheme with low-Mach correction and the orange one without the low-Mach correction. We see more kinetic energy at all scales in the case of the low-Mach correction. . 32 1.12 Rayleigh-Bénard instability simulations in 3D. The figure represents the time evolution of the moving average maximum pressure at the wall for simulations at high resolution 1024 3 and 2048 3 . We add the ideal wall as a reference. The black vertical line indicates the impact of the droplet whereas the horizontal line indicates the water hammer pressure. . . . . . . . . . 2.8 Droplet impact at resolution 1024 3 . The liquid density is represented in red.

The white zone at the back of the droplet represents the vapor volume fraction. Velocity field lines are displayed near the droplet. . . . . . . . . . . . . . . . . . .

3.1

Example of two 2D-grids used to represent the numerical solution. The mesh on the left is a standard Finite Volume mesh, each cell approximating the mean value of the exact solution. The mesh on the right is the corresponding spectral mesh with spectral cells highlighted by a blue border. Inside each Spectral Cell are represented Control Volumes, also called degrees of freedom, colored alternatively in grey and white. Both Spectral Cells and Control Volumes exactly match a set of Finite Volume Cells for a given pattern; [START_REF] Uvsq | Maison de la Simulation[END_REF]3,[START_REF] Uvsq | Maison de la Simulation[END_REF] in each direction in this example. If Finite Volume Cells are interpreted as degrees of freedom, we see that the spectral scheme can be interpreted as a compression of information when a spectral representation is adapted to the solution. . . . . . . . . . . . . . . . . . . 3.2 Example of a 2D Spectral Cell ω S is,js for a third order Spectral volume method, i.e. it has 3 Control Volumes in each direction. On the right figure the vertical green lines (resp. red lines) represent the external faces (resp. internal faces) in the X-direction where the polynomial flux F x is reconstructed. The Lagrange interpolation is realized at the center of Control Volumes represented by the black points. . We can see that due to interfaces instabilities, the number of valid Spectral Cells keeps decreasing. . . . 3.9 Snapshot of a jet simulation at Mach 80 with SVM (4) at a Finite Volume resolution 500 × 250. The top panel shows the density in a log scale. The bottom panel shows valid cells when the value is 1, i.e. when the update is done with the spectral scheme. We see that the fallback scheme is essentially enabled near discontinuities. 3.10 Snapshot of a jet simulation at Mach 80 with SVM (4) at a Finite Volume resolution 3200 × 1600. The top panel shows the density in a log scale. The bottom panel shows valid cells when the value is 1, i.e. when the update is done with the spectral scheme. We see that the fallback scheme is essentially enabled near discontinuities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Two snapshots of a Rayleigh-Taylor simulation at resolution 500 × 1500 and times t = 8.1 and t = 12.4. We observe dense material going down due to instability. . List of Tables 1.1 Isothermal atmospheres at rest. Table shows for different spatial resolutions the maximum velocity in the domain. We see that the velocity is maintained around zero up to the machine precision, thus illustrating the well-balanced property. .

This table gives the distribution of Control

Volumes in terms of Finite Volume Cells. 3.2 Summary of L 1 errors using SVM [START_REF] Uvsq | Maison de la Simulation[END_REF] on the density variable and corresponding effective orders obtained on the advection of sinusoide. . . . . . . . . . . . . . . . 3.3 Summary of L 1 errors using SVM (2) on the density variable and corresponding effective orders obtained on the advection of sinusoide. . . . . . . . . . . . . . . . 3.4 Summary of L 1 errors using SVM (3) on the density variable and corresponding effective orders obtained on the advection of sinusoide. . . . . . . . . . . . . . . . 3.5 Summary of L 1 errors using SVM (4) on the density variable and corresponding effective orders obtained on the advection of sinusoide. . . . . . . . . . . . . . . . 3.6 Summary of L 1 errors using SVM [START_REF] Uvsq | Maison de la Simulation[END_REF] on the density variable and corresponding effective orders obtained for the advection of the isentropic vortex. . . . . . . . . 3.7 Summary of L 1 errors using SVM (2) on the density variable and corresponding effective orders obtained for the advection of the isentropic vortex. . . . . . . . . 3.8 Summary of L 1 errors using SVM (3) on the density variable and corresponding effective orders obtained for the advection of the isentropic vortex. . . . . . . . . 3.9 Summary of L 1 errors using SVM (4) on the density variable and corresponding effective orders obtained for the advection of the isentropic vortex. . . . . . . . .

Introduction générale version française

Contexte général

La modélisation consiste à donner un cadre et un ensemble d'outils permettant de décrire une partie de la réalité. Ce cadre est formé à partir d'un ensemble d'hypothèses réfutables, au sens de Karl Popper, et se traduisent souvent sous forme d'équations dont la solution fournit une description approchée de phénomènes observés. L'expérimentation permet ensuite de tester ces hypothèses en comparant les résultats de manipulations ou d'observations aux solutions des équations du modèle. Avec un point de vue tourné vers les applications industrielles, la connaissance des solutions d'un modèle, alors supposé pertinent, permet de faciliter le dimensionnement de nombreux systèmes.

Parmi les différents modèles utilisés en astrophysique et dans le domaine des écoulements diphasiques se trouvent des modèles décrivant la dynamique d'un ou plusieurs fluides qui appartiennent à la catégorie des modèles de milieux continus. Ces modèles prennent la forme d'équations aux dérivées partielles et s'appuient souvent sur des principes communs sous-jacent comme l'expression de lois de conservation, par exemple la conservation de la masse, de la quantité de mouvement ou encore de l'énergie auxquelles peuvent s'ajouter d'autres équations. Hormis quelques cas simples (exploitant de nombreuses symétries du problème) les solutions de ces équations ne peuvent pas être obtenues analytiquement.

On a alors recours aux simulations numériques qui consistent à utiliser les capacités de calcul d'un ordinateur pour calculer une approximation de cette solution. Afin de pouvoir réaliser ces simulations numériques, une étape supplémentaire est nécessaire: la discrétisation. La discrétisation consiste à traduire le modèle mathématique en un algorithme composé d'un nombre fini d'opérations élémentaires et réalisables par un ordinateur par exemple: les additions, et les multiplications. L'algorithme obtenu est un nouveau modèle, discret également appelé schéma numérique, qui en un certain sens est une approximation du modèle d'origine qualifié, par opposition, de modèle continu. L'étude du lien entre le modèle continu et un schéma numérique est le domaine de l'analyse numérique au sein des mathématiques appliquées. Plusieurs schémas numériques pouvant approcher un même modèle continu, des critères d'évaluation des schémas numériques sont donc nécessaires. Ces critères peuvent prendre différentes formes par exemple: le temps de restitution de la solution numérique ou encore la précision de la solution.

Dans le cas de l'astrophysique, parmi les différents mécanismes de transfert de chaleur, l'étude de la convection naturelle est un phénomène important, au coeur de la dynamique interne des étoiles [START_REF] Spruit | Solar Convection[END_REF][START_REF] Pinsonneault | MIXING IN STARS[END_REF]) et des exoplanètes gazeuses [START_REF] Tremblin | CLOUDLESS ATMOSPHERES FOR L/T DWARFS AND EXTRASOLAR GIANT PLANETS[END_REF]), comme illustré Figure 1. Pour des raisons de temps de restitution raisonnable, l'étude par simulation numérique de la convection est souvent réalisée à partir d'approximations supplémentaires connues sous le nom d'approximation de Boussinesq ou anélastique. Ces différentes approximations permettent de filtrer les ondes sonores en supposant un écoulement convectif à faible vitesse devant la vitesse du son (régime bas Mach) et sur de petites échelles de hauteur. Cependant Figure 1: Représentation graphique des différents régimes de transfert de chaleur depuis le coeur d'une étoile vers sa surface. La chaleur émise depuis le coeur de l'étoile par les réactions nucléaires est dans un premier temps diffusé par radiation, sans mouvement du gaz. Le transfert de chaleur est ensuite effectué à partir de mouvements macroscopiques du gaz vers la surface de l'étoile. Source https://solarscience.msfc.nasa.gov/interior.shtml l'étude de régimes intermédiaires où la compressibilité du fluide joue un rôle sort du cadre de ces modèles et ne peut donc pas être capturé numériquement. L'utilisation du modèle hydrodynamique compressible est alors nécessaire mais pose également des difficultés numériques lorsque l'écoulement est proche de l'équilibre hydrostatique.

En ce qui concerne le domaine des écoulements diphasiques, dans de nombreuses situations (souvent d'origines industrielles) le liquide se trouve être en régime bas Mach. Une approximation proche du modèle incompressible est alors utilisée, présentant des similarités avec l'approximation de Boussinesq. Cependant dans des régimes intermédiaires, notamment dans les circuits hydrauliques de réacteurs nucléaires, les variations de pression au sein de gouttes d'eau liquide peuvent être importantes et sortent du cadre de ce modèle. On est alors amené à utiliser des modèles tenant compte des phénomènes compressibles. A terme, ces pics de pression peuvent endommager la paroi du circuit, comme illustré sur la Figure 2. L'étude par simulation numérique des pics de pression au niveau de la paroi fournit une approche complémentaire aux expériences.

La réalisation de ces simulations peut demander de nombreuses heures de calcul (jusqu'à plusieurs dizaines de millions d'heures de calcul). Afin de pouvoir réaliser ces simulations en des temps raisonnables l'utilisation de supercalculateurs, des ordinateurs dédiés à ce type de calcul, est nécessaire. A partir de 2021, les nouveaux supercalculateurs atteindront la puissance de calcul permettant de réaliser jusqu'à 10 18 opérations à la seconde, ce que l'on appelle couramment l'ère du calcul exaflopique ("exascale" en anglais). Afin d'obtenir cette puissance de calcul, les constructeurs s'appuient sur des architectures matérielles variées et ce qui par conséquent cela requiert l'utilisation de modèles de programmation adaptés à chacune d'entre elles pour un même schéma numérique. Cette duplication de code peut, à terme, poser des problèmes de maintenance. C'est pourquoi différentes initiatives ont donné lieu à des bibliothèques telles que Kokkos, à travers une abstraction des architectures matérielles, permettant de retrouver une

Objectifs de la thèse

Cette thèse s'inscrit à l'interface entre les aspects de modélisation astrophysique et diphasique, de développements de schémas numériques ainsi que de leur implémentation parallèle. Elle s'est déroulée au sein du laboratoire Maison de la Simulation, laboratoire multidisciplinaire dont les activités concernent des aspects de modélisation et de calcul haute performance.

Dans cette thèse nous nous intéressons au développement de schémas numériques dits "tout régime", à partir de modèles hydrodynamiques, ainsi qu'à leur implémentation dans un objectif de calcul haute performance. Ces schémas s'appuient sur un découplage des phénomènes acoustiques, a priori rapides, et des phénomènes de transport, a priori lents. Ce découplage permet entre autre de développer des schémas numériques avec diverses propriétés en facilitant leur analyse.

En ce qui concerne l'application astrophysique, nous nous intéressons à des instabilités de convection. Ces instabilités ont lieu autour de l'équilibre hydrostatique et généralement à bas nombre de Mach. Nous cherchons alors à obtenir un schéma numérique capable de simuler des écoulements convectifs à tout nombre de Mach. Nous cherchons donc à maintenir exactement, ou à précision machine près, l'équivalent discret des équilibres hydrostatiques. Les caractères bas nombre de Mach et proche équilibre hydrostatique des écoulements convectifs nécessitent un traitement particulier de la discrétisation des gradients de pression liés aux phénomènes acoustiques. Le caractère haut Mach est garanti par l'utilisation d'un modèle compressible et d'un schéma conservatif.

L'application diphasique de la thèse concerne l'impact d'une goutte d'eau sur une paroi au sein du circuit hydraulique d'un réacteur nucléaire. Les liquides et les gaz ont des propriétés thermodynamiques très différentes. La localisation précise de ces deux fluides au sein du domaine de simulation est donc importante afin de pouvoir calculer correctement les quantités thermodynamiques. De part la raideur de l'équation d'état du liquide, de fortes tensions peuvent conduire le liquide à caviter. Nous nous intéressons alors au Modèle Homogène à l'Equilibre (HEM pour Homogeneous Equilibrium Model) qui suppose, à chaque instant, l'équilibre thermodynamique entre les phases.

En lien avec ces applications nous nous intéressons également à la problématique de performance des schémas numériques et de leur implémentation. L'accélération du calcul peut être envisagée à la fois de manière algorithmique et de manière informatique. Du point de vue algo-rithmique nous nous intéressons à un schéma d'ordre élevé afin de tirer parti d'une représentation plus compacte et donc a priori moins coûteuse de la solution numérique. Du point de vue informatique, nous abordons l'aspect portabilité de performance à travers l'utilisation des bibliothèques Kokkos et Trilinos. Ces deux bibliothèques fournissent des algorithmes élémentaires, parallèles et portables, qu'il s'agit d'exploiter dans l'implémentation de nos schémas.

Description des travaux

Chapitre 1

Dans ce chapitre nous nous inscrivons dans la continuité des travaux de Chalons et al. 2016et de Chalons et al. 2017 pour le développement d'un schéma "tout-régime" et équilibre des équations d'Euler avec terme source de gravité statique. Nous cherchons à capturer les états d'équilibre hydrostatique de ces équations. La séparation des phénomènes acoustiques et de transport en deux sous-systèmes d'équations nous a conduit à placer le terme source de gravité dans le système acoustique. Le système acoustique est discrétisé à partir d'un schéma de type Godunov [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF]. Le solveur de Riemann approché associé à ce schéma est obtenu à partir d'une méthode de relaxation [START_REF] Suliciu | On the thermodynamics of rate-type fluids and phase transitions. I. Rate-type fluids[END_REF][START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF][START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF] des équations permettant une linéarisation des termes associés à l'équation d'état. Le système de transport est quant à lui discrétisé à partir d'un schéma décentré amont tout en garantissant une mise à jour globale conservative. L'observation qu'en gravité statique, par opposition à l'auto-gravité, l'énergie tenant compte du potentiel gravitationnel est également une quantité conservée, nous proposons également un schéma numérique permettant la conservation de cette quantité. De manière similaire au travail de Chalons et al. 2016 et faisant suite des analyses bas Mach des schémas colocalisés [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF][START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF], une correction de flux est apportée au gradient de pression afin d'améliorer la précision du schéma dans un régime bas Mach.

Ces schémas numériques ont ensuite été utilisés pour des applications d'instabilité de convection dans différents régimes d'écoulement. Il a été montré l'importance de la correction de flux dans des régimes bas Mach. Ils ont notamment été utilisés afin de proposer une théorie visant à unifier différents types de phénomènes convectifs.

Chapitre 2

Ce chapitre porte sur l'application diphasique de la thèse. Nous nous intéressons à la simulation d'impacts d'une goutte d'eau contre une paroi. Lorsque la goutte d'eau impacte la paroi, un important saut de pression se développe connu sous le nom de coup de bélier [START_REF] Ghidaoui | A Review of Water Hammer Theory and Practice[END_REF]) égal au produit de l'impédance acoustique de l'eau liquide par le saut de vitesse. Afin de capturer ces sauts de pression nous utilisons un modèle compressible: un Modèle Homogène à l'Equilibre [START_REF] Helluy | Relaxation models of phase transition flows[END_REF] qui suppose qu'à chaque instant l'équilibre thermodynamique est atteint. Ce modèle est approché par une relaxation instantanée et à partir d'un modèle dit "modèle à 5 équations" [START_REF] Allaire | A Five-Equation Model for the Simulation of Interfaces between Compressible Fluids[END_REF]. Ce dernier est alors discrétisé en séparant également les phénomènes acoustiques et les phénomènes de transport. Le système acoustique est discrétisé à partir d'une méthode de relaxation [START_REF] Suliciu | On the thermodynamics of rate-type fluids and phase transitions. I. Rate-type fluids[END_REF][START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF][START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF] adaptée au vide [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws[END_REF]). Le système de transport utilise un schéma anti-diffusif [START_REF] Kokh | An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model[END_REF][START_REF] Lagoutière | Modelisation mathematique et resolution numerique de problemes de fluides compressibles a plusieurs constituants[END_REF]) afin de préserver du mieux possible les interfaces présentes dans la situation initiale. Des simulations d'impact de goutte ont été menées sur différents types de paroi avec une étude de convergence. Une simulation Grand Challenge a été menée sur la moitié du supercalculateur Joliot-Curie au Très Grand Centre de Calcul (TGCC), Saclay.

Chapitre 3

Dans ce chapitre nous explorons la possibilité de diminuer le temps de restitution des simulations numériques en adaptant la représentation numérique et le schéma numérique aux propriétés locales de la solution [START_REF] Schaal | Astrophysical hydrodynamics with a high-order discontinuous Galerkin scheme and adaptive mesh refinement[END_REF][START_REF] Sonntag | Efficient Parallelization of a Shock Capturing for Discontinuous Galerkin Methods using Finite Volume Sub-cells[END_REF]. Dans une région où la solution est discontinue nous souhaitons utiliser les méthodes traditionnelles Volumes Finis par exemple à partir de schémas de type Godunov [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF]. A l'inverse, dans une région où la solution est régulière nous souhaitons compresser sa représentation numérique à l'aide de polynômes et de schémas spectraux. Nous avons donc développer une méthode hybride à partir d'un schéma MUSCL-Hancock [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]) et d'un schéma Volume Spectral [START_REF] Wang | Spectral (Finite) Volume Method for Conservation Laws on Unstructured Grids. Basic Formulation[END_REF][START_REF] Wang | Spectral (Finite) Volume Method for Conservation Laws on Unstructured Grids: II. Extension to Two-Dimensional Scalar Equation[END_REF]Wang and Liu 2004;Wang et al. 2004;[START_REF] Liu | Spectral (finite) volume method for conservation laws on unstructured grids V: Extension to three-dimensional systems[END_REF]. Ce dernier utilise une formulation intégrale des équations afin de fournir une évolution temporelle des coefficients du polynôme. Cette mise à jour est effectuée par l'intermédiaire des valeurs moyennes dans des volumes de contrôle. Cela permet une écriture très proche des schémas standards, non spectraux, et par conséquent facilite le passage d'une représentation à l'autre. Enfin des tests d'ordre de convergence et de performance ont été menés dans le cas des équations d'Euler et une extension au terme source de gravité est proposée.

Chapitre 4

Dans ce chapitre nous abordons différents outils numériques utilisés au cours de la thèse permettant de réaliser des simulations haute performance. Dans un premier temps nous présentons la problématique de portabilité de performance qui se pose avec l'arrivée des calculateurs exaflopiques. Nous présentons ensuite brièvement la bibliothèque Kokkos [START_REF] Edwards | Kokkos: Enabling manycore performance portability through polymorphic memory access patterns[END_REF] ainsi que les stratégies de parallélisation des codes s'appuyant sur ARK implémentant les schémas numériques des Chapitres 1 et 2. Nous proposons une adaptation du schéma convectif "tout-régime" aux maillages adaptatifs dans un code 2D. Enfin nous présentons les efforts de portage GPU d'un schéma implicite résolvant une équation de diffusion modélisant le transfert radiatif [START_REF] Commerçon | A fast, robust, and simple implicit method for adaptive time-stepping on adaptive mesh-refinement grids[END_REF]) dans le code RAMSES [START_REF] Teyssier | Cosmological hydrodynamics with adaptive mesh refinement: A new high resolution code called RAMSES[END_REF], code utilisé parmi la communauté astrophysique.

Publications et communications

Les travaux effectués durant cette thèse ont fait l'objet de publications et ont été présentés à plusieurs conférences internationales, à savoir

General introduction english version

Modelling consists in giving a framework and a set of tools allowing to describe part of the reality. This framework is made from a set of falsifiable assumptions, in the sense of Karl Popper, and are often represented as equations providing an approximation of observed phenomena. Then experimentation allows to test these assumptions by comparing the results of manipulations or observations to the solutions of the equations of the model. From the point of view of industrial applications, the knowledge of the solutions of a model, then assumed to be relevant, makes it possible to dimension many systems.

Among the different models used in astrophysics and in the field of two-phase flows, there are models describing the dynamics of one or more fluids that belong to the category of continuous media models. These models take the form of partial differential equations and are often based on common underlying principles such as the expression of conservation laws, for instance conservation of mass, conservation of momentum or energy to which can be added other equations. Apart from a few simple cases (that use multiple symmetries of the considered problem) the solutions of these equations cannot be obtained analytically.

Therefore we resort to use numerical simulations that consist in using the computation capabilities of a computer to calculate an approximation of this solution. In order to run these numerical simulations, a additional step is required: the discretization. This step consists in translating the equations of the model into an algorithm made of a finite number of elementary operations that a computer can realize such as additions and multiplications. The obtained algorithm is a new, discrete model also called numerical scheme, that in a sense to be defined is an approximation the first model. The study of the link between the continuous model and a numerical scheme is the field of numerical analysis within applied mathematics. Multiple numerical schemes can approximate the same continuous model, therefore criteria of evaluation are necessary. These criteria can take different forms for example: the time to numerical solution or the precision of the solution.

In the case of astrophysics, among the different heat transfer mecanisms, natural convection is an important phenomenon, at the basis of the dynamics of stars [START_REF] Spruit | Solar Convection[END_REF][START_REF] Pinsonneault | MIXING IN STARS[END_REF] and gaseous giant exoplanets [START_REF] Tremblin | CLOUDLESS ATMOSPHERES FOR L/T DWARFS AND EXTRASOLAR GIANT PLANETS[END_REF], see Figure 3. In order to have a reasonable time to solution, the numerical study of convection is often carried out from additional approximations known as the Boussinesq approximation or the anelastic approximation. These different approximations filter the sound waves assuming that the convective flow is slow compared to the speed of sound (low Mach regime) and on small scale heights. However the study of intermediate regimes where the compressibility of the fluid plays a role cannot be performed using this type of model and thus cannot be captured numerically. The use of the compressible hydrodynamics model is then required although numerical difficulties arise near the hydrostatic balance.

Regarding the field of two-phase flows, in many situations (often of industrial origin) the liquid is found to be in a low Mach regime. An approximation close to the incompressible model Figure 3: Graphical representation of different heat transfer mecanisms from the core of a star to its surface. At first, the heat generated by nuclear reactions from the core is diffused by radiation. Then the heat transfer is achieved by the advection of the gas up to the surface of the star. Source https://solarscience.msfc.nasa.gov/interior.shtml Figure 4: Formation of microstructures after 5 (left) and 10 (right) millions of liquid droplets impacts at high velocity [START_REF] Luiset | Mechanisms of stainless steels erosion by water droplets[END_REF].

is then used, showing similarities with the Boussinesq approximation. However in intermediate regimes, particularly in the hydraulic circuits of nuclear reactors, pressure variations within water liquid droplets can be significant and fall out of the scope of this model. We are then brought to use models taking into account compressible phenomena. Ultimately, these pressure peaks can damage the circuit wall, see Figure 4. The study by numerical simulations of pressure peaks at the wall provides a complementary approach to real experiments.

Performing these simulations can take many hours of calculation (up to tens of millions of hours of calculation). In order to be able to carry out these simulations in a reasonable amount of time the use of supercomputers, computers dedicated to this type of calculation, is necessary. Starting from 2021, the new supercomputers will reach a computing power allowing to perform up to 10 18 operations per second, what is commonly called the exascale computing. In order to obtain this computing power, manufacturers rely on various hardware architectures and thus it requires the use of programming models adapted to each of them for the same numerical scheme. This duplication of code may ultimately become problematic in terms of maintenance. This is why different initiatives have led to the development of libraries such as Kokkos, through an abstraction of material architectures, allowing to achieve the performance portability while writing a single version of the algorithm.

Aim of the thesis

This thesis is at the interface between the aspects of modeling astrophysics and two-phase flow, development of numerical schemes as well as their parallel implementation. It took place at the Maison de la Simulation laboratory, a multidisciplinary laboratory whose activities concern aspects of modeling and high performance computing.

In this thesis we are interested in the development of numerical schemes so-called "all regime", from hydrodynamic models, as well as their implementation for the purpose of high computing performance. These numerical schemes are based on a decoupling of the acoustic phenomena, a priori fast, and transport phenomena, a priori slow. This decoupling allows, among other things, to develop numerical schemes with various properties facilitating their analysis.

Regarding the astrophysics application, we are interested in convection instabilities. These instabilities take place near the hydrostatic equilibrium and generally at low Mach number. We then seek to obtain a numerical scheme able of simulating convective flows at any Mach number. We therefore seek to maintain exactly, or up to the machine precision, the discrete equivalent of hydrostatic equilibrium. Both the low Mach and near hydrostatic equilibrium regimes require a particular treatment of the discretization of pressure gradients associated with acoustic phenomena. The high Mach regime is guaranteed by the use of a compressible model and of a conservative scheme.

The two-phase flow application of the thesis concerns the impact of a water droplet on a wall within the hydraulic circuit of a nuclear reactor. Liquids and gases have very different thermodynamic properties. The precise location of these two fluids within the simulation domain is therefore important in order to be able to calculate correctly thermodynamic quantities. Due to the stiffness of the liquid equation of state, strong tensions can lead to cavitation within the liquid phase. We are then interested in the Homogeneous Equilibrium Model (HEM) which assumes, at each instant, the thermodynamic equilibrium between the phases.

Along with these applications we are also interested in the performance issue of numerical schemes and their implementation. The acceleration of computation can be considered in two manners. From the algorithmic point of view we are interested in a high-order scheme in order to take advantage of a more compact representation and therefore a priori less expensive numerical solution. From a computer science point of view, we address the aspect of performance portability by using the Kokkos and Trilinos libraries. These two libraries provide elementary algorithms, parallel and portable, that shall be used in the implementation of our numerical schemes. 0.6 Description of the work 0.6.1 Chapter 1

In this chapter we are in the continuity of the works initiated by [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF][START_REF] Chalons | A large timestep and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF] for the development of an "all-regime" and well-balanced scheme of Euler equations with the static gravity source term. We seek to capture the hydrostatic equilibrium states of these equations. The splitting between acoustic and transport phenomena in two subsystems of equations led us to put the gravity source term in the acoustic system. The acoustic system is discretized from a Godunov-type scheme [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF]. The approximate Riemann solver associated with this numerical scheme is obtained from a relaxation method [START_REF] Suliciu | On the thermodynamics of rate-type fluids and phase transitions. I. Rate-type fluids[END_REF][START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF][START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF] of the equations allowing the linearization of terms associated with the equation of state. The transport system is discretized using an upwind numerical scheme while ensuring an overall conservative update. Noticing that in the static gravity case, by opposition to the self-gravity, an energy taking into account the gravitational potential is satisfying a conservation law, we also propose a numerical scheme allowing its conservation. Similar to the work of [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF] and following low Mach analyses of co-located schemes [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF][START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF], a flux correction is made to the pressure gradient in order to improve the precision of the numerical scheme in a low Mach regime. These numerical schemes were then used for convection instability applications in different regimes of the flow. It has been shown the importance of the flux correction in low Mach regimes. They were also used to propose a theory aiming at unifying different types of convective phenomena.

Chapter 2

This chapter deals with the two-phase application of the thesis. We are interested in the simulation of a liquid water drop on a wall. When the water droplet hits the wall, a large jump pressure develops known as the water hammer [START_REF] Ghidaoui | A Review of Water Hammer Theory and Practice[END_REF]) that is equal to the product of the acoustic impedance of the liquid water by the velocity jump. In order to capture these pressure jumps we use a compressible model: a Homogeneous Equilibrium Model [START_REF] Helluy | Relaxation models of phase transition flows[END_REF]) that assumes at each instant an thermodynamic equilibrium. This model is approached by an instantaneous relaxation from a five-equation model [START_REF] Allaire | A Five-Equation Model for the Simulation of Interfaces between Compressible Fluids[END_REF]. The latter is then also discretized by separating the acoustic phenomena and transport phenomena. The acoustic system is discretized from a relaxation method [START_REF] Suliciu | On the thermodynamics of rate-type fluids and phase transitions. I. Rate-type fluids[END_REF][START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF][START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF] adapted to the vacuum [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws[END_REF]). The transport system uses an anti-diffusive scheme [START_REF] Kokh | An anti-diffusive numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model[END_REF][START_REF] Lagoutière | Modelisation mathematique et resolution numerique de problemes de fluides compressibles a plusieurs constituants[END_REF] in order to preserve as much as possible the interfaces present in the initial situation. Simulations of droplet impact were carried out on different types of wall along with a convergence study. A Grand Challenge simulation was carried out on half of the Joliot-Curie supercomputer at the Très Grand Centre de Calcul (TGCC), Saclay.

Chapter 3

In this chapter we explore the possibility of reducing the time to solution of numerical simulations by adapting the numerical representation and the numerical scheme to local properties of the solution [START_REF] Schaal | Astrophysical hydrodynamics with a high-order discontinuous Galerkin scheme and adaptive mesh refinement[END_REF][START_REF] Sonntag | Efficient Parallelization of a Shock Capturing for Discontinuous Galerkin Methods using Finite Volume Sub-cells[END_REF]. In a region where the solution is discontinuous we want to use standard Finite Volume methods such as Godunov-type methods [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF]. Conversely, in a region where the solution is regular we want to compress its numerical representation using polynomials and spectral schemes. We have therefore developed a hybrid method using a MUSCL-Hancock scheme [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]) and a Spectral Volume Method [START_REF] Wang | Spectral (Finite) Volume Method for Conservation Laws on Unstructured Grids. Basic Formulation[END_REF][START_REF] Wang | Spectral (Finite) Volume Method for Conservation Laws on Unstructured Grids: II. Extension to Two-Dimensional Scalar Equation[END_REF]Wang and Liu 2004;Wang et al. 2004;[START_REF] Liu | Spectral (finite) volume method for conservation laws on unstructured grids V: Extension to three-dimensional systems[END_REF]. The latter uses an integral formulation of the equations in order to provide a time evolution of the polynomial coefficients. This update is carried out through the mean value associated with Control Volumes. This allows a formulation close to standard, non-spectral diagrams, and therefore facilitates the change of one representation to the other. Finally tests of order of convergence and of performance were carried out in the case of the Euler equations and an extension to the source term of gravity is proposed.

Chapter 4

In this chapter we discuss different numerical tools used during the thesis allowing to carry out high performance simulations. First, we present the problem of performance portability that arises with the arrival of exaflop supercomputers. We then briefly present the Kokkos library (Carter [START_REF] Edwards | Kokkos: Enabling manycore performance portability through polymorphic memory access patterns[END_REF] as well as the parallelization strategies of codes based on ARK that implement the numerical schemes presented in chapters 1 and 2. We propose an extension of the convective "all-regime" scheme to the Adaptive Mesh Refinement (AMR) technique in a 2D code. Finally we present the GPU porting efforts of a implicit numerical scheme solving a diffusion equation that models the radiative transfer [START_REF] Commerçon | A fast, robust, and simple implicit method for adaptive time-stepping on adaptive mesh-refinement grids[END_REF] in the code RAMSES [START_REF] Teyssier | Cosmological hydrodynamics with adaptive mesh refinement: A new high resolution code called RAMSES[END_REF], code used among the astrophysical community.

Publications et communications

The work carried out during the thesis made the object of publications and have been presented in international conferences, namely • invited webinar, EUROfusion, July 2020, Introduction to Kokkos.

In addition, the work of Chapter 2 and of Chapter 3 are the subject of articles in preparation.

Introduction

The study of convection is an active topic of research in the astrophysics community because of its major role in different mechanisms such as heat transport in solar and stellar interiors [START_REF] Spruit | Solar Convection[END_REF]), mixing of elements [START_REF] Pinsonneault | MIXING IN STARS[END_REF]) and dynamo [START_REF] Charbonneau | Solar Dynamo Theory[END_REF]. As these mechanisms play a role in the estimation of the lifetime of these objects it is of great importance for stellar evolution theory. Different approximations have been developed to ease the study of convection. The Boussinesq and the anelastic approximations simplify the Navier-Stokes system by getting rid of acoustic waves and keeping buoyancy effects. In practice these approximations are derived by looking at the equations satisfied by small perturbations near a reference state [START_REF] Spiegel | On the Boussinesq Approximation for a Compressible Fluid[END_REF]. The Boussinesq approximation is quite restrictive as it is valid for a small layer of the reference state, such that the flow can be considered incompressible. On the other hand the anelastic approach allows to have a larger scale height by keeping the density stratification of the reference state [START_REF] Gilman | Compressible convection in a rotating spherical shell. I -Anelastic equations. II -A linear anelastic model. III -Analytic model for compressible vorticity waves[END_REF]. Another way to understand these approximations is to consider the flow regime in terms of the Mach number Ma. As it is shown in Mentrelli 2018, these approximations can be recovered by considering low-Mach asymptotic limits of the Navier-Stokes system. The Froude number, defined as the non-dimensional ratio of kinetic energy to gravitational energy, characterizes the influence of gravity in the flow. By taking into account different Froude regimes, they recover the incompressible, the Boussinesq and the anelastic models. From a numerical point of view the removal of the acoustics waves in these models is quite attractive because it allows to have larger time steps. The anelastic model has been successfully implemented in different codes like Rayleigh [START_REF] Featherstone | THE SPECTRAL AMPLITUDE OF STELLAR CONVECTION AND ITS SCALING IN THE HIGH-RAYLEIGH-NUMBER REGIME[END_REF] or Magic [START_REF] Gastine | Effects of compressibility on driving zonal flow in gas giants[END_REF] and it is widely used in the community (see [START_REF] Glatzmaier | Introduction to Modeling Convection in Planets and Stars[END_REF]. We can also mention the MAESTRO code (see [START_REF] Nonaka | MAESTRO: AN ADAPTIVE LOW MACH NUMBER HYDRODYNAMICS ALGORITHM FOR STELLAR FLOWS[END_REF] which uses an extended version of the anelastic model. The velocity constraint takes into account the time variation of pressure. However these approaches present some drawbacks. The addition of new physics and source terms to the model is difficult, one has to derive another asymptotic model to take the new physics into account in the anelastic regime (see [START_REF] Mentrelli | Modelling of the convective plasma dynamics in the Sun: anelastic and Boussinesq MHD systems[END_REF]. Furthermore one has to be careful that the simulation stays in the regime of validity of the model (especially in the Boussinesq regime). Finally a numerical difficulty is the parallelization of those codes. They usually use pseudo-spectral methods for which it is more difficult to achieve a good scalability (e.g. need to use pencil-type domain decomposition [START_REF] Featherstone | THE SPECTRAL AMPLITUDE OF STELLAR CONVECTION AND ITS SCALING IN THE HIGH-RAYLEIGH-NUMBER REGIME[END_REF].

We chose to take a more flexible approach by solving the full compressible Navier-Stokes system, as in the MUSIC code [START_REF] Viallet | Towards a new generation of multi-dimensional stellar evolution models: development of an implicit hydrodynamic code[END_REF][START_REF] Goffrey | Benchmarking the Multidimensional Stellar Implicit Code MUSIC[END_REF]) but with a collocated finite volume solver instead of using a staggered grid. Different discretization techniques of the Euler system are used in the astrophysics community. We can classify them in various ways. One way is to separate SPH techniques from grid-based techniques. Furthermore grid-based approaches can be divided in different families, finite difference, finite element and finite volume. The finite volume method is of particular interest because of its natural property of being conservative and to capture shocks and discontinuities. Designing a finite volume scheme essentially resides in the definition of a numerical flux, numerical counterpart of the physical flux. A widely used family of fluxes is the Godunov (see [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF] flux which is the flux of the -usually approximate -Riemann problem between two neighbour cells.

However we have to face multiple numerical difficulties with this approach. Compressible solvers and mainly Godunov-type solvers are known to have an excessive amount of numerical diffusion in the low-Mach regime which make them unusable in this regime (see [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF][START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF][START_REF] Miczek | New numerical solver for flows at various Mach numbers[END_REF][START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF][START_REF] Barsukow | A Numerical Scheme for the Compressible Low-Mach Number Regime of Ideal Fluid Dynamics[END_REF]. In this regime, in which flows are smoother, considering Riemann problems at interfaces is not adapted. Indeed in the work of [START_REF] Miczek | New numerical solver for flows at various Mach numbers[END_REF] they show that part of the kinetic energy is dissipated into internal energy whereas it should be conserved. To tackle this issue they propose a preconditionned Roe scheme to remove the numerical diffusion. Secondly, hydrodynamics and gravity are usually discretized independently from each other. In the case of highly stratified medium, the numerical scheme does not maintain the hydrostatic equilibrium and produces spurious flows that pollutes the simulation. Different approaches have been investigated to solve this issue both for the Euler and the shallow water equations. In [START_REF] Leroux | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF], they rewrite the Euler system as a fully conservative system by defining an hydrostatic pressure satisfying a conservation law. In Chandrashekar and Klingenberg 2015 they use a variable reconstruction by taking advantage of the equilibrium profile. In [START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter dependent source. The case of Euler system with friction[END_REF][START_REF] Vides | A Godunov-Type Solver for the Numerical Approximation of Gravitational Flows[END_REF][START_REF] Chalons | A large timestep and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF] they incorporate the source term in the Riemann problem itself allowing to compensate pressure gradients at the interface. As in [START_REF] Leroux | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF]Cargo 1994, authors from Chertock et al. 2018 also propose to discretize the Euler system with gravity as a fully conservative system but using global fluxes and a reconstruction on equilibrium variables. Finally the last numerical difficulty is the time step in the low-Mach regime. Because of the stability condition involving the fast acoustic waves, the time step becomes very small compared to the material transport timescale. It can either be resolved using a full implicit approach as in the MUSIC code [START_REF] Viallet | Towards a new generation of multi-dimensional stellar evolution models: development of an implicit hydrodynamic code[END_REF][START_REF] Goffrey | Benchmarking the Multidimensional Stellar Implicit Code MUSIC[END_REF], or by using an implicit-explicit (IM-EX) approach in which only the system with fast acoustic waves is solved implicitly [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF][START_REF] Chalons | A large timestep and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF].

Following the original work of [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF][START_REF] Chalons | A large timestep and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF] we use an acoustictransport splitting. In [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF] they derive a finite volume scheme of the Euler system on unstructured mesh. This scheme uses an acoustic splitting to separate acoustic waves from material ones. In the low Mach regime, this translates to a splitting between fast waves and slow waves. In the low Mach regime, the fast waves can be treated with an implicit solver to get rid of the restrictive stability condition. Then in the work of [START_REF] Chalons | A large timestep and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF], the scheme has been adapted to shallow water equations with a source term which is the topography. This source term is added in the equivalent acoustic subsystem to obtain a well-balanced scheme. In this paper we adapt their approach for the Euler system by taking care of the discretization of the energy equation.

The paper is organized as follows. In Section 1.1 we briefly recall the compressible model we use to study convection, i.e. the Navier-Stokes equations with gravity. In Section 1.2 we present the derivation of the well-balanced and all-regime numerical scheme using a splitting approach between an acoustic step and a transport step both solved explicitly in this work. In Section 1.3 we present some implementation features about the "ARK"1 code in particular the Kokkos library used for the shared memory parallelization. We also give some performance results. Finally in Section 1.4 we present different numerical test cases illustrating the importance of the low-Mach correction and the well-balanced discretization of gravity.

Navier-Stokes equations

We want to solve Navier-Stokes equations expressing conservation of mass, balance of momentum and balance of energy, respectively written as follows

∂ t ρ + div (ρu) = 0, ∂ t (ρu) + div (ρu ⊗ u + pI -τ visc ) = ρg, ∂ t (ρE) + div ((ρE + p) u -τ visc u + q heat ) = ρg • u, (1.1)
where ρ is the density, u the material velocity, p the pressure, g the external gravitational field, ρE = ρe + 1 2 ρu 2 the density of total energy with e the specific internal energy, q heat the heat flux and τ visc the viscous tensor satisfying

τ visc = µ ∇u + ∇u T + η (divu) I, (1.2)
where µ is the dynamic viscosity and η the bulk viscosity. We use • as a scalar product and thus div represents the divergence operator. In order to close Navier-Stokes system (1.1) we add constitutive equations namely a pressure law p EOS (1.3a), the Fourier's law (1.3b) and the Stokes hypothesis (1.3c)

p = p EOS (ρ, e) , (1.3a
)

q heat = -κ∇T, (1.3b) η = - 2 3 µ (1.3c)
We recall that the gravitational field is derived from a gravitational potential Φ for which g = -∇Φ. Dealing with a constant in time external gravity field, ∂ t Φ = 0 and using the conservation of mass we get (1.4)

∂ t (ρΦ) + div (ρΦu) = ρu • ∇Φ. (1.4)
Let us emphasize that in this equation, the gravitational energy ρ(x, t)Φ(x) is time dependent only through the density ρ(x, t). Hence the energy equation (1.4) can be rewritten in the following conservative form

∂ t (ρE) + div (ρEu -σ stress u + q heat ) = 0 (1.5)
where we define ρE = ρe + 1 2 ρu 2 + ρΦ. Equation (1.5) expresses the local conversion between three different energy reservoirs, as depicted in figure 1.1: internal, kinetic and gravitational. There can be a direct transfer between gravitational energy and kinetic energy through the

Gravitational energy

Internal energy Kinetic energy Work of gravitational forces Dissipation work of gravitational forces, from kinetic energy to internal energy because of the second law of thermodynamics but no direct transfer between gravitational energy and internal energy, see also Section 5 of Springel 2010 and Section 2.2 of Marcello and Tohline 2012 for a discussion on energy conservation for both external and self-gravity cases.

Because of this conservation of energy including gravitational energy we will use the formulation (1.5) of the energy equation and we will use the gravitational potential instead of the usual gravitational field g. To our knowledge this approach is quite rare, see [START_REF] Graham | Numerical simulation of two-dimensional compressible convection[END_REF][START_REF] Chertock | Well-balanced schemes for the Euler equations with gravitation: Conservative formulation using global fluxes[END_REF] where they use global fluxes to have a well-balanced and conservative scheme.

An important steady state solution of this system for stratified objects is the hydrostatic balance. The flow is static and the gravitational force is balanced by the pressure forces, i.e. following equation (1.6) ∇p = -ρ∇Φ, u = 0, (1.6)

As we mentioned in the introduction, convective flows can be considered as a perturbation flow of the hydrostatic equilibrium. Thus this steady state is particularly important in order to study convection problems in stratified flows.

1.2 Numerical scheme

Euler system -Hyperbolic system

Before going into the derivation of the scheme we introduce the notations. We define by ∆x (resp. ∆y and ∆z) the step along the x-direction (resp. the y and z-direction). We note by ∆t the time interval between current time t n and t n+1 . We use the notation q n i (resp. q n i,j,k ) to represent the averaged quantity associated to the field q at time t n and in the cell i (resp. i, j, k) in the onedimensional case (resp. the three-dimensional case). We use the notation q n i+1/2 (resp. q n i+1/2,j,k ) to represent the quantity associated to the field q at time t n and at the interface between cells i and i + 1 (resp. i, j, k and i + 1, j, k) in the one-dimensional case (resp. the three-dimensional case). Finally we define the notation [q] i = q i+1/2 -q i-1/2 in the one-dimensional case.

Acoustic-Transport splitting approach

Following [START_REF] Chalons | A large timestep and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF] we use a splitting strategy that separates acoustic terms and transport terms and we choose to add the gravitational source terms to the acoustic part. This way, pressure gradient can be balanced by the gravity source term.

However we have another equation compared to the shallow water system that is the energy equation. As in [START_REF] Chalons | A large timestep and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF], we want an isentropic acoustic step for smooth solutions. Thereby we choose to solve the equation on the gravitational energy,

∂ t ρ + div (ρu) = 0, ∂ t (ρu) + div (ρu ⊗ u + pI) = -ρ∇Φ, ∂ t (ρE) + div ((ρE + p) u) = 0, ∂ t (ρΦ) + div (ρΦu) = ρu • ∇Φ, (1.7) ρE = ρe + 1 2 ρu 2 + ρΦ.
However, this leads to a non constant gravitational potential in the acoustic step whose time variations are compensated in the transport step in order to have a constant potential in the full step. The potential is constant in the full step at the continuous level, but discretization errors with the splitting can lead to a non-constant discretized potential. Thus we choose to introduce an approximation of the gravitational called Ψ ≈ Φ and a relaxation parameter λ

∂ t ρ + div (ρu) = 0, ∂ t (ρu) + div (ρu ⊗ u + pI) = -ρ∇Φ, ∂ t (ρE) + div ((ρE + p) u) = 0, ∂ t (ρΨ) + div (ρΨu) = ρu • ∇Φ + ρ λ (Φ -Ψ) , (1.8
)

ρE = ρe + 1 2 ρu 2 + ρΨ.
We consider the relaxation system (1.8) to be an approximation of the original system (1.7) that we formally recover in the limit λ → 0. System (1.8) is solved by first solving the system in the limit λ → ∞ and then in the limit λ → 0 in which Ψ is projected onto Φ, the initial condition. This way, the evolution of the gravitational potential Ψ, consistent with zero, is forced to be constant. The relaxation technic used here for the gravitational potential is similar to what is done for pressure relaxation in many approximate Riemann solvers and we emphasize that Ψ is just an intermediate used to design the scheme and can be removed when writing the final scheme (see 1.2.1).

We now turn to the discretization of the system (1.8) in the limit λ → ∞. Transport phenomena of the form u • ∇ are separated from the other terms to give two subsystems, first the acoustic subsystem (1.9)

∂ t ρ + ρdivu = 0, ∂ t (ρu)+ρudivu + ∇p = -ρ∇Φ, ∂ t (ρE) + ρEdivu + div (pu) = 0, ∂ t (ρΨ)+ρΨdivu = ρu • ∇Φ, (1.9)
then the transport subsystem (1.10)

∂ t ρ + u • ∇ρ = 0, ∂ t (ρu) + u • ∇ (ρu) = 0, ∂ t (ρE) + u • ∇ (ρE) = 0, ∂ t (ρΨ) + u • ∇ (ρΨ) = 0.
(1.10)

We now briefly study the eigenstructure of systems (1.9)-(1.10). Let n be any unit normal vector, the acoustic system (1.9) involves seven eigenvalues: -c, 0, c. The fields associated with 0 (resp. ±c) are linearly degenerate (resp. genuinely nonlinear), see Appendix 1.A for more details. The eigenvalues for transport system (1.10) are given by u • n. Both systems (1.9)-(1.10) are hyperbolic. We emphasize here that the choice of using a relaxation procedure for the gravitational potential by introducing the equation on the gravitational potential energy ρΨ has been made to obtain this simple wave pattern for the splitted Euler system with gravity. (i.e. the same pattern as without gravity). Other choices for the relaxation procedure (e.g. ∂ t Ψ = 0 in both steps) would either lead to the introduction of u • n in the eigenvalues of the acoustic subsystem or would significantly complexify the relaxation procedure for the pressure.

To summarize our numerical procedure, we propose to define a flux interface by approximating system (1.7) with a three-step procedure that involves solving the acoustic system (1.9) (acoustic step), the transport system (1.10) (acoustic step) and finally project Ψ onto Φ (relaxation step). We detail each step in the next sections using the one-dimensional equations.

Acoustic step

We follow the idea of [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF] to discretize the acoustic subsystem. They introduce a pressure relaxation Π ≈ p, an acoustic impedance a ≈ ρc and a relaxation parameter ν to get a fully linearly degenerated system. It is then written using Lagrangian variables (τ, u, v, E, Ψ) where u represents the normal velocity component at an interface and v a transverse component. We also use a mass variable dm = ρ(t n , x)dx where time is frozen at instant t n

∂ t τ -∂ m u = 0, ∂ t u + ∂ m Π = - 1 τ ∂ m Φ, ∂ t v = 0, ∂ t E + ∂ m (Πu) = 0, ∂ t Π + a 2 ∂ m u = 1 ν (Π -p) , ∂ t Ψ = u τ ∂ m Φ,
where

E = e + 1 2 (u 2 + v 2 ) + Ψ.
The discretization of this system is realized with an approximate Riemann solver that accounts for the source term by means of integral consistency and composed by three waves -a, 0, a, see [START_REF] Gallice | Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source[END_REF][START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF][START_REF] Chalons | A large timestep and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF]. After the relaxation, in which ν → 0, it gives

τ i = τ n i + ∆t ∆m i [u * ] i , u i = u n i - ∆t ∆m i [Π * ] i + ∆t ∆m i S n i , v i = v n i , E i = E n i - ∆t ∆m i [Π * u * ] i , Π i = p EOS 1 τ i , e i , Ψ i = Ψ n i - ∆t ∆m i (uS) n i ,
where

u * i+1/2 = 1 2 (u n i+1 + u n i ) - 1 2a Π n i+1 -Π n i -S n i+1/2 , Π * i+1/2 = 1 2 Π n i+1 + Π n i - a n i+1/2 2 u n i+1 -u n i , a n i+1/2 ≥ max ρ n i c n i , ρ n i+1 c n i+1 , S n i = 1 2 S n i+1/2 + S n i-1/2 , (uS)
n i = 1 2 (u * i+1/2 S n i+1/2 + u * i-1/2 S n i-1/2 ), S n i+1/2 = - 1 2 1 τ n i + 1 τ n i+1 Φ n i+1 -Φ n i .
and a n i+1/2 ≥ max ρ n i c n i , ρ n i+1 c n i+1 which is the so-called sub-characteristic condition (see [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF].

The update of the conservative variables is then

L i ρ i = ρ n i , L i (ρu) i = (ρu) n i - ∆t ∆x [Π * ] i + ∆t ∆x S n i , L i (ρv) i = (ρv) n i , L i (ρE) i = (ρE) n i - ∆t ∆x [Π * u * ] i , L i (ρΨ) i = (ρΨ) n i - ∆t ∆x (uS) n i
where

L i = 1 + ∆t ∆x [u * ] i .

Transport step

The transport subsystem can be written in the following form, for b ∈ {ρ, ρu, ρv, ρE, ρΨ}

∂ t b + ∂ x (bu) -b∂ x u = 0, that is discretized as follows b n+1 i = b i - ∆t ∆x bu * i + b i ∆t ∆x [u * ] i .
The interface term b i+1/2 is defined by the upwind choice with respect to the velocity

u * i+1/2 b i+1/2 = b i if u * i+1/2 ≥ 0 b i+1 if u * i+1/2 ≤ 0

Relaxation step

At this stage, the relaxed gravitational potential Ψ still evolves in time. So we perform the relaxation λ → 0 that boils down to set

Ψ n+1 i = Φ i .

Overall algorithm

Gathering the previous steps and intermediate variables, the overall scheme reads

ρ n+1 i = ρ n i - ∆t ∆x [ ρu * ] i , (ρu)
n+1 i = (ρu) n i - ∆t ∆x (ρu)u * + Π * i + ∆t ∆x S n i , (ρv)
n+1 i = (ρv) n i - ∆t ∆x (ρv)u * i , (ρE)
n+1 i = (ρE) n i - ∆t ∆x (ρE) + Π * u * i (1.11)
It may also be expressed as a first-order classic finite-volume scheme involving flux terms for the conservative part for energy ρE = ρe + 1 2 ρu 2 and source terms for gravity

ρ n+1 i = ρ n i - ∆t ∆x [ ρu * ] i , (ρu)
n+1 i = (ρu) n i - ∆t ∆x (ρu)u * + Π * i -∆t{ρ∂ x Φ} i , (ρv)
n+1 i = (ρv) n i - ∆t ∆x (ρv)u * i , (ρE)
n+1 i = (ρE) n i - ∆t ∆x (ρE) N G + Π * u * i -∆t{ρu∂ x Φ} i , (1.12)
where

∆x{ρu∂ x Φ} i = [ ρu * Φ] i -[ ρu * ] i Φ i , ∆x{ρ∂ x Φ} i = -S n i , (ρE)
N G i = (ρE) n i - ∆t ∆x [Π * u * ] i .
We emphasize that both formulations are equivalent and conservative for the energy ρE. A non-conservative energy approach is also detailed in Appendix 1.B. We can notice that in the case of a constant gravitational potential, we recover the original scheme derived in [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF].

On the low-Mach correction

As for the scheme of [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF] and as explained in Dellacherie 2010, the numerical scheme defined by (1.11) poorly performs in the low Mach regime due to truncature error of magnitude ∆x Ma that comes from the term Π * i+1/2 . To tackle this issue, following Chalons et al. 2016 we modify the upwinding part of Π * i+1/2 thanks to an extra parameter θ i+1/2 by setting

Π * i+1/2 = 1 2 Π n i+1 + Π n i - a n i+1/2 θ i+1/2 2 u n i+1 -u n i , θ i+1/2 = min Ma i+1/2 , 1 , Ma i+1/2 = |u * i+1/2 | max c n i , c n i+1 . (1.13)
Using a truncation analysis in dimensionless form it can be shown that this correction acts like a rescaling of the numerical diffusion induced by the pressure discretization (see [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF]).

As we can see, the low-Mach correction does not directly come from the derivation of the numerical scheme 1.11. Some ongoing works are trying to derive directly all-Mach schemes using more sophisticated relaxation schemes (see [START_REF] Bouchut | An entropy satisfying twospeed relaxation system for the barotropic Euler equations[END_REF]).

On the well-balanced property

A numerical scheme is said to be well-balanced for equilibrium states satisfying equation (1.6), if it exists a discrete counterpart of equation (1.6) in which solutions are preserved by the numerical scheme.

The discrete counterpart of equation (1.6) for scheme (1.11) is given by

u n i = 0, v n i = 0, Π n i+1 -Π n i = - 1 2 ρ n i + ρ n i+1 (Φ i+1 -Φ i ) , (1.14)
Let us now verify that we have obtained a well-balanced scheme. If at time t n , for some density profile the initial state reads as in (1.14) then fluxes from the acoustic step reduce to

u * i-1/2 = u * i+1/2 = 0 [Π * ] i = 1 2 Π n i+1 -Π n i + 1 2 Π n i -Π n j-1 + S n i .
Then we have for the acoustic step

u i = u n i , v i = v n i , ρ i = ρ n i , E i = E n i .
Finally, because u * i+1/2 vanishes, transport step is trivial and the initial state remains unchanged. Once we have made the appropriate choice for the discretization of the gravitational source term in the acoustic step, the well-balanced property is automatically verified without the need to introduce an other algorithmic correction.

Dissipative fluxes -Parabolic system

We now turn to the discretization of dissipative fluxes (1.2)- (1.3b). They are discretized using second order discrete fluxes

divf dissipative i,j,k = (f x,i+1/2,j,k -f x,i-1/2,j,k ) ∆x + (f y,i,j+1/2,k -f y,i,j-1/2,k ) ∆y + (f z,i,j,k+1/2 -f z,i,j,k-1/2 ) ∆z KNL Skylake K80 P100 V100 10 1 10 2 10 3 Performance [Mcell-updates/s] Figure 1
.2: Comparison of performance on different architectures: Intel KNL, Intel Skylake (one socket), NVIDIA K80, NVIDIA P100 and NVIDIA V100. Measures on Intel KNL and Intel Skylake were performed on Joliot-Curie's supercomputer at TGCC using the same code.

In our case we obtain better results with the Intel Skylake than the Intel KNL due to a lack of vectorization. Going to a GP-GPU we have a speed-up around five with a NVIDIA K80 compared to multi-core architecture and seven between NVIDIA K80 and V100.

where f dissipative is either the heat flux q heat or the viscous flux τ visc . In the case of the heat flux we have

q x,i+1/2,j,k = -κ (T i+1,j,k -T i,j,k ) ∆x q y,i,j+1/2,k = -κ (T i,j+1,k -T i,j,k ) ∆y q z,i,j,k+1/2 = -κ (T i,j,k+1 -T i,j,k ) ∆z
With the addition of the viscous terms and the heat flux, this all-regime well-balanced scheme is now well-suited for the study of convection problems in highly stratified flows in both low Mach and high Mach regimes. Before showing validating numerical tests, we present some specificities about the numerical implementation and parallelization used in this work.

Implementation and parallelization

In this section we describe the implementation of the scheme using Kokkos library. We begin by giving a brief overview of the Kokkos library.

Exascale computing

To reach the exascale, the distributed memory model is not sufficient to take advantage of all the computing power of new architectures. There are mainly two reasons for this. First, nodes of supercomputers tend to grow more and more and hence are more suited to a shared memory model [START_REF] Sunderland | An Overview of Performance Portability in the Uintah Runtime System through the Use of Kokkos[END_REF]. Secondly, nodes tend be more and more heterogeneous by using multi-core, many-core and/or accelerators like GP-GPUs. So it means that even if shared memory is exposed, it needs to be handled differently from one architecture to another. For example we can think of OpenMP or C++11 threads for multi-core and many-core processors, and CUDA or OpenACC for GP-GPUs. Moreover this architecture heterogeneity raises a performance portability issue. Currently, many HPC codes are optimized for some specific architectures to get the maximum computing power. However this optimization process couples the numerical scheme to its implementation details like the memory management, the loop ordering, cache blocking and so on. Hence running a code on a different architecture results in bad performance.

We propose to use the recent C++ library Kokkos (see Carter [START_REF] Edwards | Kokkos: Enabling manycore performance portability through polymorphic memory access patterns[END_REF]) that implements a new shared memory model. Using abstract concepts such as execution spaces (where a function is executed), data spaces (where data resides) and execution policies (how the function is executed) the library is able to efficiently take advantage of multi-core many-core processors and GP-GPUs. This way the portability relies on the library and no more on the numerical code.

Implementation

Following the work of Kestener 2017, the code is then organized with computation kernels:

• Acoustic and transport kernels,

• Viscous and heat diffusion operator kernels,

• Conservative variables to primitive variables kernel,

• Time step kernel.

Each kernel is a C++ functor. They are given to Kokkos through the function Kokkos::parallel_for. Internally, depending on the device chosen at compile-time, it hides a parallelized one-dimensional loop where the current index is given as an argument to the functor. This index is then interpreted as a cell index in the domain.

Kokkos only deals with shared memory systems. We use the Message Passing Interface (MPI) programming model with a regular domain decomposition to take advantage of distributed memory machines across multiples computing nodes. Kokkos is then used as a shared memory programming model inside each node. These domains are endowed with ghost zones which are used to both implement physical boundary conditions and to contain values from neighbour domains. Communications are handled through the ghost cell pattern (see [START_REF] Kjolstad | Ghost Cell Pattern[END_REF]. Thus for a given direction X, Y (or Z) and a given side, left or right, one MPI process sends data from its domain to its neighbour's ghost zone and receives data into its own ghost zone.

Performance results

Thanks to Kokkos, we were able to use the same code on different architectures like Intel Skylake, Intel Knights Landing (KNL) and NVIDIA GP-GPUs (K80, P100, V100). We measured performance on the Intel Skylake and the Intel KNL partition of the Joliot Curie machine at TGCC. Figure 1.2 shows the results. We see that the Kokkos library is able to provide good performance on the different tested architectures. Nevertheless, even if the peak performance of the Intel KNL architecture is higher than the Intel Skylake one we have better performance on We use a hybrid MPI-OpenMP configuration in which one MPI task is bound to a socket. Simulations run for 1000 time steps and each MPI process treats 128 3 cells. We see that the efficiency reaches a plateau of 85%.

the Intel Skylake architecture. We also notice the important speed-up (around five) between the Intel Skylake architecture and the NVIDIA V100 GP-GPU.

Figure 1.3 shows a weak scaling test performed with a hybrid configuration OpenMP/MPI. We went up to 512 MPI processes, one MPI process per Intel Skylake socket to avoid NUMA effects. It results in a total of 12288 cores at 512 MPI processes. Each MPI process is getting a piece of the whole domain of 128 3 , so a domain of 44 3 per core. We can see that we obtain a plateau of 85% of maximum performance from 128 MPI processes.

The performances obtained with the use of the Kokkos library are encouraging for the study of convection problems with the ARK code on massively parallel present and future architectures. In the next section, we use several numerical tests to show that the numerical scheme used in the ARK code is indeed very well suited for the study of convection.

Numerical results

In this section we specialize the equation of state 1.3a. We will use an ideal gas satisfying p EOS (ρ, e) = (γ -1) ρe where γ is the adiabatic index of the gas. The speed of sound satifies the following simple relation

c 2 = γ p ρ 0.0 0.2 0.4 0.6 0.8 1.0 x 0.2 0.4 0.6 0.8 1.0 ρ Exact solution All Regime [θ = Ma] All Regime [θ = 1] Godunov-type [order 1]
Figure 1.4: Sod's test case simulations. Figure shows a snapshot of the density profile ρ for the All-Regime scheme, with and without the low-Mach correction, a first order Godunov-type scheme (HLLC) and the exact solution. Spatial resolution is nx = 100. We see that the All-Regime scheme gives results close to the Godunov-type scheme around discontinuities but is more diffusive in the rarefaction wave.

We emphasize that it is possible to use a different equation of state with the all-regime wellbalanced numerical scheme. Moreover we consider two versions of the all-regime scheme depending on the low-Mach correction. We will refer to the disabled low-Mach correction scheme when θ = 1 and to the enabled one when θ follows equation 1.13.

We will test different properties of the scheme with different test cases: wave speeds with the Sod test (no gravity), low-Mach accuracy with the Gresho vortex test (no gravity), hydrostatic balance with the test of an atmosphere at rest and out of equilibrium behavior with the Rayleigh-Taylor test. We then use the ARK code for the study of Rayleigh-Bénard convection.

Shock tube test

The Sod shock tube [START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF]) is a classical test for compressible solvers. It tests the ability of the solver to have correct wave speeds and its numerical diffusion near discontinuities.

The computational domain is the interval [0,1], the initial condition is defined by

(ρ, p, u) = (1, 1, 0) if x < 0.5, (0.125, 0.1, 0) if x ≥ 0.5.
Results are shown in figure 1.4 for simulations with nx = 100. First we can observe that the solver is as good as a first order Godunov-type scheme with a HLLC approximate Riemann solver around the contact discontinuity and the shock. However the rarefaction wave is a bit more diffused. We also notice that the low-Mach correction does not influence the behavior of the scheme for this test case. However we want to stress out some instability near discontinuities, as shown in [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF]. This can also be seen in a double shock waves test case.

Gresho vortex test case

The Gresho vortex [START_REF] Gresho | On the theory of semi-implicit projection methods for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix. Part 2: Implementation[END_REF][START_REF] Miczek | New numerical solver for flows at various Mach numbers[END_REF]) is a test case that has already been used to test numerical schemes in the low Mach regime. It is a two dimensional stationary test case that can be parameterized by the maximum value of the Mach number. It is thus well-suited to study the behavior of the scheme in the low Mach regime. We recall that the test case is defined using polar coordinates (r, θ) defined with respect to the center of the vortex as follows

ρ = ρ 0 , (u r , u θ ) =      (0, 5r) 0 ≤ r < 0.2, (0, 2 -5r) , 0.2 ≤ r < 0.4, (0, 0) , 0.4 ≤ r p =      p 0 + 12.5r 2 , 0 ≤ r < 0.2, p 0 + 12.5r 2 + 4 -20r + 4 ln(5r), 0.2 ≤ r < 0.4, p 0 -2 + 4 ln 2, 0.4 ≤ r.
where p 0 satisfies p 0 = 1 γMa 2 . In this case Ma is a parameter and γ is the adiabatic index of the ideal gas. The velocity is normalized so a particle placed at the peak of velocity (u = 1.0 at location r = 0.2) make a full rotation in ∆t = 2 5 π ≈ 1.26. We ran a serie of simulations with different solvers where we explored parameter space nx and Ma from 32 to 2048 and from 1.0 to 1.0 × 10 -5 respectively. Final time is set to t f = 1.0 × 10 -3 , which has been chosen sufficiently small such that the error doesn't saturate.

Figure 1.5 shows snapshots of the the magnitude of the velocity field at the final time and at resolution 512 2 . We see that when the Mach number decreases the velocity field becomes more and more degraded when the low-Mach correction is disabled. At Ma = 10 -5 , the vortex has completely disappeared. Figures 1.6 and 1.7 show more quantitative results where we show absolute L 1 error on velocity in function of the Mach number Ma and the spatial resolution dx respectively. Figure 1.6 shows that L 1 error on velocity depends on the Mach number. More precisely we measure a slope of -1 on schemes or order 1 and a slope of -0.5 on scheme of order 2. On the other hand the low Mach correction of the all-regime scheme gives a uniform error with respect to the Mach number. see that both Godunov-type and all-regime without the low Mach correction converge at order 1 as expected. Nevertheless Godunov-type with Muscl-Hancock reconstruction converges only at order 1.5. It may be due to the lack of regularity of the velocity field as it can be observed in the case of a contact discontinuity (see [START_REF] Springel | E pur si muove: Galilean-invariant cosmological hydrodynamical simulations on a moving mesh[END_REF]. All-Regime scheme shows two different behaviors, at first it converges at order 1.5 then around nx = 1024 the slope changes and it converges at order 1.2. We assume that at higher resolution we would recover order 1. We see that at low Mach number the precision, independently of the order, is better than the one of a Godunov-type scheme.

Well-balanced test case

The well-balanced test case is a simple isothermal column of atmosphere at equilibrium. This column of atmosphere is in a stable equilibrium state. The test allows us to measure the ability of the scheme to preserve this equilibrium. After normalization, it is given by

p(z) = ρ(z) = e -z T = 1
which is the solution of the following system

dp dz = -ρ dΦ dz T = 1 p = ρT
We take advantage of the formula (1.14) and we initialize the test case with the following formula

p i+1 -p i ∆z = - ρ i + ρ i+1 2 Φ i+1 -Φ i ∆z T i = 1 p i = ρ i T i
The computational domain used is the interval [0, 3]. Results are displayed in table 1.1 at time t = 10, more than three times the sound crossing time in the box. We see that we stay near machine precision at the end of the simulation. We see a shift of two orders of magnitude in the error when using the low-Mach correction. The reason of this shift is not entirely clear and is difficult to interpret as it involves truncature errors. Looking at the spatial pattern of the error in the simulation, it does seem to come from the boundary conditions (extrapolation of the hydrostatic balance for pressure and density and reflexive conditions for the velocity) with the use of the low-Mach correction. A more appropriate boundary condition might remove this shift in the error (which is in any case sufficiently small and stable to allow the use of controlled seeded perturbations).

ρ (x, y) = 1 for y < 0 2 for y >= 0 p (x, y) = ρgy u (x, y) = 0 v (x, y) = C 4 1 + cos 2πx L x 1 + cos 2πy L y
Where C = 0.01 is the magnitude of the velocity perturbation, L x = 0.5 and L y = 1.5 are the size of the domain in each direction. We do not need to use the well-balanced formula (1.14), the equilibrium is preserved in the case A = 0. Figure 1.8 shows two simulations of the Rayleigh-Taylor test case, one with the low Mach correction and the other without it (θ n i+1/2 = 1). Both simulations are at the same time t = 12.4 and the same resolution 200 × 600. The yellow part is at density 2 and the purple is at density 1. We see that we recover the classical linear growing mode. Moreover the simulation with the low Mach correction is able to capture secondary instabilities in the non linear regime. They are closer to the second order Godunov-type simulation than the order one. However the low Mach correction does not help on the interface diffusion between the two mediums. It also shows a peak that is not present without the correction at the same resolution. This spurious behavior is therefore caused by the low Mach correction that removes some numerical diffusion in the scheme. By looking at higher resolutions, we identify that this peak is a grid-seeded secondary RT unstable mode that appears at the top of the large scale seeded mode. This type of secondary modes are not unexpected and can be seen for example in Fig. 9 of [START_REF] Almgren | CASTRO: A NEW COMPRESSIBLE ASTROPHYSICAL SOLVER. I. HYDRODYNAMICS AND SELF-GRAVITY[END_REF]. This peak disappears with the addition of some physical viscosity in the simulation.

Rayleigh-Bénard instability test case

This last test case is about compressible convection simulations both in 2D and 3D. In this test case there are different important parameters. First, from stability analysis we know that the Rayleigh number Ra is an important non-dimensional number. Beyond a threshold, called First line show results with the the all-regime scheme, where on the left the low-Mach correction is disabled and is enabled on the right. Second line shows results with a Godunov-type scheme, on the left it is first order, on the right it is second order using a Muscl-Hancock scheme. We see that with the low-Mach correction we recover features only present at second order for a standard Godunov-type scheme.

the critical Rayleigh Ra c , the convection process starts and efficiently transports the heat (see Figures 1 and 3 in Hurlburt et al. 1984). Below this threshold, diffusion processes are sufficient to transport heat and no material displacement is necessary. Then another important parameter is the density stratification χ which the ratio between the density at the bottom of the domain and the density at the top. In the highly stratified case, study of convection becomes more difficult as there is not a unique Rayleigh number but more a whole range of values extending on the scale height. Notice that when χ → 1 we recover the Boussinesq-like situation. Finally the last parameter is the polytropic index m which is a measure of how close is the initial temperature gradient from the adiabatic gradient. One can show that the Schwarzschild criterion writes

m + 1 < γ γ -1 = 2.5, γ = 5 3
The initial setup is inspired from [START_REF] Hurlburt | Two-dimensional compressible convection extending over multiple scale heights[END_REF][START_REF] Toomre | Three-dimensional compressible convection at low prandtl numbers[END_REF]. Following their notation, the initial state is given by a polytropic profile of polytropic index m

T = z, ρ = z m , p = z m+1
where z is the vertical variable. It is initialized using to the recursive formula (1.14). So we begin with a hydrostatic equilibrium that we destabilize whether with a velocity mode perturbation or with a temperature random perturbation.

2D case

We begin with 2D simulations in a weak stratification setup where χ = 1.1 and m = 1.3 in order to be close to the adiabatic gradient. The initial perturbation is close to the fundamental velocity mode. The spatial resolution is set to 128 2 , and we impose the temperature flux on the bottom boundary. We then obtain stationary symmetric convective rolls. We study the effect of the low-Mach correction on the onset of the Rayleigh-Bénard instability by varying the initial Rayleigh number. Figure 1.9 shows the evolution of the mean absolute velocity. The linear phase, in logarithmic scale, corresponds to the exponential growth of modes. We can see that If we now turn to a stronger stratification, the convective rolls pattern change. We increase the density ratio to χ = 21. Figure 1.10 shows a snapshot of the local Mach number field with the velocity field, low-Mach correction enabled. As observed in [START_REF] Hurlburt | Two-dimensional compressible convection extending over multiple scale heights[END_REF] we see a downward shift of the center of mass of convective rolls compared to the weak stratification case. By conservation of mass, the upper part of the convective roll has to be larger. The strong stratification case also exhibits higher Mach flows, around Ma ≈ 0.5 at the top of the box due to the low density. The all-regime well balanced scheme is indeed able to capture properly convection in highly stratified and high Mach flows.

3D case

We now turn to 3D simulations in a weak stratification situation. In this setup we want to look at the effect of the low-Mach correction on the kinetic energy spectrum in a more turbulent situation. So we change the polytropic index to m = 0.1 and increase the initial Rayleigh number to Ra ≈ 650000. We also change the boundary condition to a fixed temperature for both boundaries in order to continuously force a large Rayleigh number in the simulation. By using different upscaling, from 128 3 to 512 3 we reach a stationary state2 . Figure 1.12 shows a snapshot of the velocity in the box. We see large and structured vertical flows whereas in horizontal plans the flow is more turbulent. In order to study the different scales and the energy in this turbulent state we compute power spectrum of the kinetic energy of the horizontal middle plane. correction and the blue one without it. We see a net difference in the overall kinetic energy due to a lower dissipation into the internal energy. We notice that we recover higher kinetic energies at all scales showing that the low Mach correction is important to properly capture the power spectrum of turbulent convection.

Conclusion

We have presented a new numerical code that is able to perform simulations of convection without any approximation of Boussinesq nor anelastic type. To do so we have adapted an all-Mach number scheme into a well-balanced scheme for gravity. We have been able to show that it preserves arbitrary discrete equilibrium states up to the machine precision. Moreover the low-Mach correction in the numerical flux allows to be more precised in the low-Mach regime. This new scheme is well suited to properly study highly stratified and high Mach convective flows.

The low Mach correction is important to properly capture convection modes in the laminar low Mach regime and the kinetic energy power spectrum in the turbulent regime. This code has been parallelized using a hybrid approach MPI+Kokkos in order to be well prepared for running on forthcoming exascale machines. Further work will consist in using the implicit-explicit approach to reach very low Mach number simulations, see [START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF], and still keeping the well-balanced property for the gravity source term. Indeed by solving the acoustic part implicitly we avoid the restrictive CFL condition due to the fast acoustic waves. With both the explicit-explicit and implicit-explicit approach, this numerical scheme will be able to efficiently study convection problems in all regimes, low Mach and high Mach on the largest next generation massively parallel architectures.

Appendix

1.A Eigenstructure of the acoustic system

For the sake of simplicity, the eigenstructure analysis of the acoustic system (1.9) is made in the one-dimensional case. We use the following change of variables, valid for smooth flows

(ρ, ρu, ρE, ρΨ, Φ) → (ρ, u, s, Ψ, Φ),
where s is the specific entropy. By using equation of mass, one obtains

∂ t ρ + ρ∂ x u = 0, ∂ t u + 1 ρ ∂ x p EOS + ∂ x Φ = 0, ∂ t e - p ρ ∂ x u = 0, ∂ t Ψ -u∂ x Φ = 0, ∂ t Φ = 0.
By using the second law of Thermodynamics and the equation on the specific internal energy, one can show that ∂ t s = 0 (see [START_REF] Godlewski | Hyperbolic systems of conservation laws[END_REF]. Thus the acoustic system (1.9) writes equivalently

∂ t ρ + ρ∂ x u = 0, ∂ t u + 1 ρ ∂ x p EOS + ∂ x Φ = 0, ∂ t s = 0, ∂ t Ψ -u∂ x Φ = 0, ∂ t Φ = 0.
(1.15)

The Jacobian matrix associated to the quasi-linear system 1.15 involves five eigenvalues: -c < 0 < c where 0 is degenerated three times and c satifies c 2 = ∂ ρ p EOS (ρ, s). It is then hyperbolic.

The four eigenvectors are given by

r 0 0 = ∂ s p, 0, -c 2 , 0, 0 T , r 1 0 = (0, 0, 0, 1, 0) T , r ±c = (ρ, ±c, 0, 0, 0) T .
Clearly the field associated to the stationary wave is linearly degenerated. The fields associated to ±c are genuinely non-linear under the condition that the following quantity does not vanish

±∇c (ρ, s) • r ±c = ±ρ∂ ρ c = ± ρ 2c ∂ 2 ρρ p EOS .

1.B Non-conservative energy scheme

To obtain the non-conservative scheme, we do not need anymore the relaxation on the gravitational potential. This scheme is then obtained through the following splitting, for the acoustic subsystem

∂ t ρ + ρdivu = 0, ∂ t (ρu) + ρudivu + ∇p = -ρ∇Φ,
∂ t (ρE)+ρEdivu + div (pu) = -ρu • ∇Φ, followed by the transport subsystem

∂ t ρ + u • ∇ρ = 0, ∂ t (ρu) + u • ∇ (ρu) = 0, ∂ t (ρE) + u • ∇ (ρE)
= 0, then we use the same techniques for the acoustic system as in 1.2.1, in other words the use of the mass variable and the Lagrangian variables. The acoustic system in these variables writes

∂ t τ -∂ m u = 0, ∂ t u + ∂ m p = - 1 τ ∂ m Φ, ∂ t v = 0, ∂ t E + ∂ m (pu) = - u τ ∂ m Φ, E = e + 1 2 (u 2 + v 2 ).
Using a pressure relaxation, an approximate Riemann solver with source term, see [START_REF] Gallice | Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source[END_REF], and the same upwind scheme for the transport system as in 1.2.1 we obtain the following nonconservative counterpart scheme

ρ n+1 i = ρ n i - ∆t ∆x [ ρu * ] i , (ρu)
n+1 i = (ρu) n i - ∆t ∆x (ρu)u * + Π * i + ∆t ∆x S n i , (ρv)
n+1 i = (ρv) n i - ∆t ∆x (ρv)u * i , (ρE)
n+1 i = (ρE) n i - ∆t ∆x (ρE) + Π * u * i + ∆t ∆x (uS) n i , where u * i+1/2 = 1 2 (u n i+1 + u n i ) - 1 2a Π n i+1 -Π n i -S n i+1/2 , Π * i+1/2 = 1 2 Π n i+1 + Π n i - a 2 u n i+1 -u n i , a n i+1/2 ≥ max ρ n i c n i , ρ n i+1 c n i+1 , S n i = 1 2 S n i+1/2 + S n i-1/2 , (uS)
n i = 1 2 (u * i+1/2 S n i+1/2 + u * i-1/2 S n i-1/2 ), S n i+1/2 = - 1 2 1 τ n i + 1 τ i+1 (Φ i+1 -Φ i ) .
This scheme is not conservative for the whole energy but is closer to the scheme proposed for the shallow water equations in [START_REF] Chalons | A large timestep and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF], for which the authors have obtained a discrete entropy inequality. It seems therefore possible to obtain a similar inequality for this non-conservative scheme, but this demonstration is beyond the scope of this paper.

1.C Generalized convection, diabatic convection

Convective instabilities such as standard adiabatic convection, moist convection or even doublediffusive convection are usually studied independently. In this section we present a unified derivation of the instability criteria by introducing a general source term on the energy equation of Navier-Stokes equations and a fraction that parameterizes the equation of state. As it has been presented in [START_REF] Tremblin | Thermo-compositional Diabatic Convection in the Atmospheres of Brown Dwarfs and in Earth's Atmosphere and Oceans[END_REF], this allows to consider other convective instabilities such as the radiative convective instability in the CO/CH 4 transition in the atmosphere of brown dwarfs and exoplanets. Even though all these convective instabilities involve different mechanisms they have in common to either influence local density or temperature which are the key ingredients to convection movements. By opposition to adiabatic convection, we qualify the convection involving such mechanisms as diabatic convection. Omitting the momentum viscosity, the Navier-Stokes system writes

∂ t ρ + div (ρu) = 0, ∂ t (ρu) + div (ρu ⊗ u + pI) = ρg, ∂ t (ρE) + div ((ρE + p) u) = ρc p H (C, p, T ) , ∂ t (ρC) + div (ρCu) = ρR (C, p, T ) , (1.16)
p = p eos (ρ, e, C) .

(1.17)

in which H is a heat source and R is a composition source. Both source terms and the equation of state depend on a fraction C. In Sections 1.C.1 and 1.C.2, we derive the Boussinesq system from a rescaled Navier-Stokes system, omitting the source terms H and R. This system consists in studying Navier-Stokes equations on a small length scale in which pressure variations are considered to be small thus removing pressure dependence on density variations. Then in Section 1.C.3 we derive the instability criteria which are obtained by studying stability of eigen modes after linearization of the Boussinesq model.

1.C.1 Non-dimensional form

Before deriving the Boussinesq system, we write the Euler equations in a non-dimensional form. We choose l 0 , u 0 , ρ 0 , p 0 , g 0 to be some reference quantities. For the reference speed of sound we choose c 2 0 = p0 ρ0 . We also choose the internal energy as total energy reference ρ 0 E 0 = ρ 0 e 0 = p 0 . By rescaling variables as f = f 0 f , the Euler system writes

ρ 0 t 0 ∂t ρ + ρ 0 u 0 l 0 div (ρ ū) = 0 ρ 0 u 0 t 0 ∂t (ρ ū) + 1 l 0 div ρ 0 u 2 0 ρ ū ⊗ ū + p 0 pI = ρ 0 g 0 ρḡ ρ 0 E 0 t 0 ∂t ρ Ē + 1 l 0 div ρ 0 E 0 ρ Ē + p 0 p u 0 ū = ρ 0 g 0 u 0 ρḡ • ū
This system can be rewritten as follows:

Sr∂t ρ + div (ρ ū) = 0 Sr∂t (ρ ū) + div ρ ū ⊗ ū + 1 Ma 2 pI = 1 Fr 2 ρḡ Sr∂t ρ Ē + div ρ Ē + p ū = Ma 2 Fr 2 ρ ū • ḡ ρ Ē = ρē + Ma 2 ρ 2 ū 2
where we introduced the Strouhal number, the Mach number and the Froude number:

Sr = l 0 u 0 t 0 Fr 2 = u 2 0 g 0 l 0 = ρ 0 u 2 0 ρ 0 g 0 l 0 = 2 E kinetic E gravitational Ma 2 = u 2 0 c 2 0 = u 2 0 p 0 /ρ 0 = ρ 0 u 2 0 ρ 0 e 0 = 2
E kinetic E internal that are respectively the ratio of kinetic energy over gravitational energy and the ratio of kinetic energy over internal energy. In the following we consider that Sr = 1 and we omit the bar to simplify notations.

1.C.2 Boussinesq regime

In this section we derive the Boussinesq system from Euler system with gravity. We follow the approach in Mentrelli 2018; Feireisl and Novotný 2017 using an Hilbert expansion in terms of a non-dimensional small parameter . However we choose a different set of variables for the Hilbert expansion. Moreover we set H = 0 and R = 0.

We start the analysis from the rescaled Euler system

∂ t ρ + div (ρu) = 0 ∂ t (ρu) + div ρu ⊗ u + 1 Ma 2 pI = 1 Fr 2 ρg ∂ t (ρE) + div (u (ρE + p)) = Ma 2 Fr 2 ρg • u p = p eos (ρ, e) ρE = ρe + Ma 2 ρ 2 u 2 , ρE = ρE + Ma 2 Fr 2 ρφ (1.18)
We are interested in the Boussinesq regime, which we define to be the limit of the latter system when Ma = , Fr = √ and Bo = Ma Fr = √ for → 0. We refer the reader to Feireisl and Novotný 2017 for the scaling leading to the anelastic system. In order to derive the limit we use an Hilbert expansion, variables are then written as f = n n f n . Because we are interested in smooth solutions, System 1.18 can be written in non-conservative variables u, p and T

ρ (∂ t u + u • ∇u) + 1 2 ∇p = 1 ρg ∂ t p + u • ∇p + ρc 2 divu = 0 ρc p (∂ t T + u • ∇T ) = αT (∂ t p + u • ∇p) ρ = ρ eos (p, T ) (1.19) in which α (p, T ) = -1 ρ ∂ T ρ eos (p, T )
is the isobaric thermal expansion coefficient. We also introduce the isothermal compressibility coefficient χ (p, T ) = 1 ρ ∂ p ρ eos (p, T ). Temperature equation is obtained from entropy equation and using Maxwell relations. Hilbert expansion for orders O -2 , O -1 and O 0 reads

∇p 0 = 0 ∇p 1 = ρ 0 g ∂ t p 0 + u 0 • ∇p 0 + ρ 0 c 2 0 divu 0 = 0 ρ 0 (∂ t u 0 + u 0 • ∇u 0 ) + ∇p 2 = ρ 1 g ρ 0 c p 0 (∂ t T 0 + u 0 • ∇T 0 ) = α 0 T 0 (∂ t p 0 + u 0 • ∇p 0 ) ρ 0 = ρ eos (p 0 , T 0 ) (1.20)
First equation imposes uniform pressure p 0 (t, x) = p 0 (t), using compatible boundary conditions we also have a constant pressure p 0 (t) = p 0 (0). From pressure equation, this imposes a divergence-free velocity field divu 0 = 0. Hence, temperature equation reduces to a transport equation. Assuming uniform initial temperature profile and compatible boundary conditions we deduce that temperature is uniform in space and constant in time, T 0 (t, x) = T 0 (t) = T 0 (0). Using equation of state, ρ eos we also deduce that density is uniform and constant. Density ρ 1 is unknown, we need to go one order further, it writes

ρ 0 c p 0 (∂ t T 1 + u 0 • ∇T 1 ) = α 0 T 0 (∂ t p 1 + u 0 • ∇p 1 ) ρ 1 = ∂ p ρ (p 0 , T 0 ) p 1 + ∂ T ρ (p 0 , T 0 ) T 1 (1.21)
Akin to p 0 , assuming that p 1 is constant in time we obtain

ρ 0 c p 0 (∂ t T 1 + u 0 • ∇T 1 ) = α 0 T 0 u 0 • ∇p 1 ρ 1 = -ρ 0 α 0 T 1 + ρ 0 χ 0 p 1 (1.22)
We introduce the so-called adiabatic temperature gradient, in hydrostatic equilibrium

∇T ad = α 0 T 0 ρ 0 c p 0 ∇p 1 = α 0 T 0 c p 0 g (1.23)
Boussinesq system then writes

divu 0 = 0 ∂ t u 0 + u 0 • ∇u 0 + 1 ρ 0 ∇p 2 = ρ 1 ρ 0 g ∂ t T 1 + (∇T 1 -∇T ad ) • u 0 = 0 ρ 1 = -ρ 0 α 0 T 1 + ρ 0 χ 0 p 1 (1.24)
It is standard to write the system in variables relative to the hydrostatic equilibrium, we emphasize that this is no more an approximation. Hydrostatic solution follows

u h 0 = 0 ∇p h 2 = ρ h 1 g ρ h 1 = -ρ 0 α 0 T h 1 + ρ 0 χ 0 p 1 (1.25)
40 CHAPTER 1. AN ALL-REGIME HYDRODYNAMIC SOLVER FOR STRATIFIED FLOWS Noting f = f -f h variables relative to the hydrostatic equilibrium, Boussinesq system 1.24 also write divu 0 = 0

∂ t u 0 + u 0 • ∇u 0 + 1 ρ 0 ∇p 2 = ρ1 g ∂ t T1 + u 0 • ∇ T1 + ∇T h 1 -∇T ad • u 0 = 0 ρ1 = -ρ 0 α 0 T1 (1.26)
To summarize the expansion of variables we have

T = T 0 + T h 1 + T1 + 2 T 2 + . . . , T 1 = T h 1 + T1 ρ = ρ 0 + ρ h 1 + ρ1 + 2 ρ 2 + . . . , ρ 1 = ρ h 1 + ρ1 p = p 0 + p h 1 + 0 + 2 p h 2 + p2 + . . . , p h 1 = p 1 , p 2 = p h 2 + p2
(1.27)

1.C.3 Diabatic convection instability analysis

To fix ideas the equation of state is chosen to be an ideal gas for which mean molecular weight m depends on a fraction / concentration

C p = ρ k B m (C) T (1.28)
As an example of source term in this section we consider thermal conduction given by Fourier's law H (C, p, T ) = H (T ) = σ T ∆T and molecular diffusion given by Fick's law R (C, p, T ) = R (C) = σ C ∆C. Both source terms have the same structure with σ T being the thermal diffusion coefficient and σ C the molecular diffusion coefficient. These source terms correspond to the double diffusive convection. Other source terms are developed in [START_REF] Tremblin | Thermo-compositional Diabatic Convection in the Atmospheres of Brown Dwarfs and in Earth's Atmosphere and Oceans[END_REF]. Following notations introduced in Section 1.C.2, the Boussinesq system writes divu 0 = 0

∂ t u 0 + u 0 • ∇u 0 + 1 ρ 0 ∇p 2 = ρ1 g ∂ t T1 + u 0 • ∇ T1 + ∇T h 1 -∇T ad • u 0 = H C h 1 + C1 , p 0 , T h 1 + T1 -H C h 1 , p 0 , T h 1 ∂ t C1 + u 0 • ∇ C1 + ∇C h 1 • u 0 = R C h 1 + C1 , p 0 , T h 1 + T1 -R C h 1 , p 0 , T h 1 ρ1 = -ρ 0 α 0 T1 + ρ 0 β 0 C1 (1.29) with β 0 = 1 ρ0 ∂ C ρ eos (C 0 , p 0 , T 0 )
and in which we neglect pressure variations of source terms H and R. By linearizing the system near equilibrium we have

divu 0 = 0, ∂ t u 0 + 1 ρ 0 ∇p 2 = ρ1 g, ∂ t T1 -w 0 ∇T h 1 -∇T ad = H C C1 + H T T1 , ∂ t C1 -w 0 ∇C h 1 = R C C1 + R T T1 , ρ1 = -ρ 0 α 0 T1 + ρ 0 β 0 C1 . (1.30) in which ∇T h 1 = -∇T h 1 • e z , ∇T h ad = -∇T h ad • e z and ∇C h 1 = -∇C h 1 • e z .
In order to study stability of perturbations, a standard method is to first derive dispersion relation by seeking for non-vanishing modes of the form f (t, x) = f e ωt+ik•x in a periodic box. Solutions of this form are said to be unstable if the real part of the dispersion relation is positive, i.e. (ω) = (ω (k)) > 0, the case (ω) = 0 defines the marginal stability (which defines the critical Rayleigh number in the adiabatic situation). It can be noticed that other boundary conditions in the z-direction generally lead to a differential system in variable z resulting in a shift of the marginal stability. Even though this analysis needs to be adapted to each situation, it gives qualitative results regarding physical mechanisms. It results in a linear system given by

k • û0 = 0, ρ 0 ω û0 + ik p2 = ρ1 g, ω T1 -ŵ0 ∇T h 1 -∇T h ad -ĤC Ĉ1 -ĤT T1 = 0, ω Ĉ1 -ŵ0 ∇C h 1 -RC Ĉ1 -RT T1 = 0, ρ1 = -ρ 0 α 0 T1 + ρ 0 β 0 C1 . (1.31)
In the particular case of thermal conduction we have ĤT T1 = -σ T k 2 T1 and for molecular diffusion RC Ĉ1 = -σ C k 2 Ĉ1 . For a given wave-number k, we look for non-zero solutions of the linear system 1.31 which are then characterized by a non-zero determinant of the associated matrix. The problem now reduces to find values for ω which both cancel out the determinant and have negative real part. It can be shown that the problem reduces to solving a third degree polynomial. Finally unstable solutions follow one of the two inequalities

∇T h 1 -∇T h ad α 0 -∇C h 1 β 0 > 0, (1.32) ∇T h 1 -∇T h ad α 0 RC + β 0 RT -∇C h 1 β 0 ĤT + α 0 ĤC < 0. (1.33)
The first inequality 1.32 corresponds to the known Ledoux criterion for the convective region of stars. The second criterion involves the coupling of source terms and is referred as diabatic criterion. Again in the particular case of thermal conduction and molecular diffusion, inequality 1.33 reduces to the standard double-diffusive criterion 1.34

∇T h 1 -∇T h ad α 0 σ T -∇C h 1 β 0 σ C > 0.
(1.34)

In addition to double diffusion convection Equation 1.33 can also describe the moist convective instability in earth atmospheres when source terms are replaced by condensation/evaporation of liquid water. Note that this general theory has been derived thanks to the simulations that were performed with code ARK. Other applications such as the CO/CO 2 radiative convection in the atmosphere of hot rocky exoplanets are under development (Daley-yates et al in prep., see Figure 1.13). We also point out that the diabatic theory with condensation/evaporation of steam/liquid water could also be of interest for two-phase flow applications in nuclear engineering, the subject of the next chapter.

Chapter 2

Droplet impact

Introduction

In the context of a nuclear reactor, the increase of the power leads to the decrease of the pressure in the secondary circuit and to a decrease of the vapor mass fraction down to 97 %. Liquid droplets can form from the steam matrix. The velocity of the droplets is estimated between 30 m s -1 and 100 m s -1 . When impacting the pipes, for example during a change of direction, liquid droplets can damage them. An estimation of the pressure developed at the moment of the impact can help in the design of the circuitry to better cope with potential high local pressure values.

A few liquid droplet impact experiments have been realized, as for example in the work of [START_REF] Field | The effects of target compliance on liquid drop impact[END_REF]. They could record snapshots of the density gradient during the impact of a high velocity two-dimensional droplet, see snapshots below 2.2. From the numerical point of view, some simulations investigated the liquid impacts. In the work of [START_REF] Wu | On high-speed impingement of cylindrical droplets upon solid wall considering cavitation effects[END_REF], authors explore cylindrical droplets impacts at very high velocities, up to 200 m s -1 . Three phases are involved, liquid water, steam water and air. They use a pressure threshold to trigger cavitation along with an instantaneous relaxation. They find good agreement in comparing Schlieren images. In [START_REF] Kyriazis | Modelling cavitation during drop impact on solid surfaces[END_REF], authors model the phases as a mixture at mechanical and thermal equilibrium. They have found a large cavitation area at the back of the droplet.

In this work we investigate droplet impacts by means of a Homogeneous Equilibrium Model. We begin by introducing this model and how we relate it to the five-equation model using a relaxation. The second section presents numerical discretization of these models. The third section presents the numerical validation of the model and numerical simulations of droplet impacts at high resolution.

Model

Before discussing models and discretization we first introduce necessary notations used in the following. Vectors are noted v and matrices T, in particular identity matrix is noted I.

We study a particular case of two-phase flows that involves water in two phases: liquid and gas. In the following we use as synonyms the terms gas, vapor and steam. We use the index α to generically refer to the phases: α = 0 (resp. α = 1) will refer to the gas (resp. liquid). Phasic quantities are indexed by α, for example ρ α . We call pure state a state which is only composed of either gas or liquid and a mixture state otherwise.

Following notations from Chapter 1, each phase is characterized by a density ρ α = 1 τα , a specific internal energy e α and a velocity u α . Moreover we note s α and T α , respectively, the specific phasic entropy and the phasic temperature. We assume that both phases, gas and liquid, are described with given equations of state (ρ α , s α ) → p eos α (ρ α , s α ). These equations of state are such that we can define a real valued speed of sound

c α 2 = ∂p eos α ∂ρ α (ρ α , s α ) > 0. (2.1)
We also recall the definition of the phasic specific chemical potential

g α = e α + p α ρ α -T α s α .
(2.2)

Description of the mixture

As in the mono-fluid case, the mixture is characterized by a density ρ, a specific internal energy e, a velocity u and fractions that link the phasic quantities to the mixture quantities.

The volume fraction z α is the volume occupied by one phase in a given volume. It is used to link the mixture densities and the specific energies to their phasic equivalents

ρ = α ρ α z α , ρe = α ρ α z α e α .
(2.3)

These two equations come from the extensivity properties of mass and energy. Furthermore we consider that the two phases are immiscible. This directly translates to

α z α = 1. (2.4)
This allows us to note z = z 1 = 1 -z 0 .

The mass fraction, y α is the fraction of mass of one phase in a given volume, i.e.

y α = ρ α ρ z α .
(2.5)

By definition it satisfies α y α = 1, thus we also note y = y 1 = 1 -y 0 . We choose a mixture equation of state of the form p = p eos (ρ, e, z, y) which depends on the two fractions z and y. In the following we consider the thermodynamics equilibrium equation of state and the isobaric equation of state.

Different definitions can be chosen for p that will characterize the mixture. There exists other closure models related to mixtures such as isothermal and isochoric, see [START_REF] Lagoutière | Modelisation mathematique et resolution numerique de problemes de fluides compressibles a plusieurs constituants[END_REF][START_REF] Allaire | A Five-Equation Model for the Simulation of Interfaces between Compressible Fluids[END_REF][START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF]; Helluy and Seguin 2006.

Isobaric mixture, p eos isob

We consider an isobaric closure that simply writes p eos isob (ρ, e, z, y) = p eos 1 (ρ 1 , e 1 ) = p eos 0 (ρ 0 , e 0 ) .

(2.6)

For a given tuple (ρ, e, z, y) and together with the specific internal energy relation e = α e α , the pressure equality allows to give a value to phasic specific internal energies e α , i.e. to the phasic temperatures T α .

Remark 1 (Domain of validity) It follows that the liquid pressure must remain positive because of the gas phase. Then this doesn't allow to take into account meta-stable states of liquid like stretching that can lead to negative pressures. However, depending on phasic equations of state, this formula may allow to derive an effective equation of state with a larger domain of validity as it is the case with stiffened gases. • the equilibrium that corresponds to a pure state,

• the equilibrium that corresponds to a mixture state.

Pure states

In the first case, the equilibrium equation of state reduces to the phasic equation of state, that is to say

p eos isob (ρ, e, z, y) = p eos 1 (ρ, e) , if z 1 = y 1 = 1, p eos 0 (ρ, e) , if z 0 = y 0 = 1.
(2.7)

Mixture states

In the second case, we consider that the two phases have the same pressure, temperature and chemical potential that is to say

τ = (1 -y) τ 0 + yτ 1 , e = (1 -y) e 0 + ye 1 ,
p eos 0 (ρ 0 , e 0 ) = p eos 1 (ρ 1 , e 1 ), T eos 0 (ρ 0 , e 0 ) = T eos 1 (ρ 1 , e 1 ), g eos 0 (p 0 , T 0 ) = g eos 1 (p 1 , T 1 ).

(2.8)

The equality of chemical potentials allows to define a saturation curve T sat (p). By inverting the system of Equations 2.8 both in τ and e, it writes (2.10)

y = τ -τ 0 τ 1 -τ 0 = e -e 0 e 1 -e 0 . ( 2
Knowing the equilibrium pressure p we can then recover all necessary thermodynamics quantities at equilibrium T eq = T sat (p) , (τ α ) eq = τ eos α (p, T ) , (e α ) eq = τ eos α (p, T ) ,

y eq = τ -τ 0 τ 1 -τ 0 , z eq = τ 1 τ y.
(2.11)

Finally, we notice that the equation of state at thermodynamics equilibrium is only determined by the mixture density ρ and the specific mixture energy e. Fractions y = y eq (ρ, e) and z = z eq (ρ, e) can be deduced from the equilibrium. This is similar to having only one equation of state used for both phases, as in IAPWS [START_REF] Wagner | The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use[END_REF] Remark 2 IAPWS-95, which stands for International Association for the Properties of Water and Steam, is a standard that formulates properties of water using a serie expansion of the Helmholtz free energy in terms of variables (ρ, T ). This serie involves up to 64 terms which are obtained from experimental measures. Other thermodynamics quantities can be recovered from the Helmholtz free energy.

The Homogeneous Equations Model

The velocity of the water liquid droplet in the pipe is estimated between 30 and 100 m/s. Because of this high velocity, we have to take into account the compressibility of liquid water. We choose an Homogeneous Equations Model, also shortly referred as HEM, that assumes that both phases are at mechanical and thermodynamics equilibrium at each time t. This translates to a velocity equilibrium and to the usage of the thermodynamics equilibrium equation of state u = u 1 = u 0 , p = p eos eq (ρ, e).

(2.12)

This model simply takes the form of a single fluid Euler equations

∂ t ρ + div (ρu) = 0,
∂ t (ρu) + div ρuu + p eos eq I = 0, ∂ t (ρE) + div ρE + p eos eq u = 0, ρE = ρe + 1 2 ρu 2 , y = y eq (ρ, e) , z = z eq (ρ, e) ,

(2.13) that expresses the conservation of mixture quantities: the mixture mass ρ, the mixture momentum ρu and the mixture total energy ρE. The pressure p eos eq = p eos eq (ρ, e) is the thermodynamics equilibrium pressure between liquid and gas. We can distinguish two cases in this model, pure phase regions (for which mass fraction is exclusively 0 or 1) and mixture regions.

First in the pure phase case, then Equations 2.13 locally degenerate to the Euler system with p = p eos α (ρ, e), in which mixture variables degenerate to phasic variables. The speed of sound also reduces to the phasic speed of sound c α . On the other hand, in the mixture case, mass transfer between the two phases occurs as mass fraction instantaneously adapts to local variations in density and internal energy. As already mentioned time evolution of the mass fraction solely depends on the time evolution of mass and internal energy equations, formally ∂ t y eq (ρ, e) = ∂ ρ y eq (ρ, e) ∂ t ρ + ∂ e y eq (ρ, e) ∂ t e.

(2.14)

This allows to potentially create new phase inside region of pure phase, in case of strong pressure drop for instance.

Remark 3 The kinematic properties of a mixture lead to different physics. For example the Baer and Nunziato model (for example see [START_REF] Baer | A two-phase mixture theory for the deflagration-todetonation transition (ddt) in reactive granular materials[END_REF] allows to model a sub-scale velocity drift between the two phases.

Remark 4 In terms of mass transfer we can distinguish models in which it occurs at the interface between two fluids or using a local relaxation, possibly instantaneous, to thermodynamics equilibrium.

In the first case, mass transfer takes the form of a source term at the interface. Velocity and pressure at interface are not anymore continuous, see [START_REF] Fechter | A sharp interface method for compressible liquid-vapor flow with phase transition and surface tension[END_REF]. The interface is moving at its own velocity that is determined by jump relations. The difficulty related to such models is the generation of new interfaces.

In the second case, mass transfer takes the form of a local source term. One assumes that out of equilibrium, the evolution is relaxing towards thermodynamics equilibrium at some rate λ, which may be instantaneous. The difficulty with such model is the choice of the relaxation parameter. The particular choice of λ = ∞ falls back to using models known as Homogeneous Equations Models, HEM.

The relaxed five-equation model

In order to solve the system of Equations 2.13 we notice that we can formally rewrite the system of Equations 2.13, in the limit of λ → ∞, with the following system

∂ t z + u • ∇z = λ (z eq -z) , ∂ t (ρy) + div (ρyu) = λρ (y eq -y) , ∂ t ρ + div (ρu) = 0, ∂ t (ρu) + div (ρuu + p eos isob I) = 0, ∂ t (ρE) + div (ρE + p eos isob ) u = 0, ρE = ρe + 1 2 ρu 2 ,
(2.15) in which we introduce two new equations on the volume fraction z and the mass fraction y and the isobaric mixture equation of state p eos isob . On the one hand, when λ → ∞, we formally recover the former HEM system with p eos isob → p eos eq along with the following compatibility relation p eos isob (ρ, e, z eq (ρ, e), y eq (ρ, e)) = p eos eq (ρ, e) .

(2.16)

On the other hand, when λ = 0, the model reduces to the following five-equation model

∂ t z + u • ∇z = 0, ∂ t (ρy) + div (ρyu) = 0, ∂ t ρ + div (ρu) = 0, ∂ t (ρu) + div (ρuu + p eos isob I) = 0, ∂ t (ρE) + div (ρE + p eos isob ) u = 0, ρE = ρe + 1 2 ρu 2 .
(2.17)

It is an hyperbolic system with five wavespeeds given by, in the direction n,

{u • n -c eos isob , c eos isob , c eos isob , c eos isob , u • n + c eos isob } , (2.18)
with c eos isob satisfying

ξ α = ρ α (∂ pα e eos α ) ρα , ξ = α z α ξ α , ρξ(c eos isob ) 2 = α ρ α z α ξ α (c eos α ) 2 .
(2.19)

The relaxation 2.15 is motivated by numerical reasons. Indeed as it is presented in Kokh and Lagoutière 2010, along with an anti-diffusive numerical scheme and an isobaric equation of state, the five-equation model is able to transport interfaces with low numerical diffusion without generating spurious pressure oscillations. This property is important to avoid smearing of the initial interface between the liquid droplet and the steam matrix, see also Section 2.A for a connection between interface models and mixture models.

Other relaxations can be considered as in the work of [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF]; Helluy and Seguin 2006 in which they derive relaxation models based on entropy considerations. These models allow to consider non-instantaneous relaxation.

Discretization

We start this section by introducing some notations required for the discretization. Let Ω ⊂ R 3 be a Cartesian domain

Ω = [x min , x max ] × [y min , y max ] × [z min , z max ].
It is discretized using a Cartesian regular mesh. We define by ∆x (resp. ∆y and ∆z) the step along the X-direction (resp. the Y and Z-direction). In the X-direction (resp. Y, Z-direction) interfaces are located at x i+1/2 = (i + 1)∆x (resp. y j+1/2 = (j + 1)∆y, z k+1/2 = (k + 1)∆z). Time is discretized into ordered instants t n , separated by intervals ∆t and we identify t n and n.

Let φ : (x, t) → φ(x, t) be some scalar field, we note φ n i,j,k the Finite Volume approximation of φ in a given cell (i, j, k) at instant t n . Finally we note φ n i+1/2,j,k (resp. φ n i,j+1/2,k , φ n i,j,k+1/2 ) the quantity associated to the field φ at the time t n and at the interface between cells (i, j, k) and (i + 1, j, k) (resp. (i, j, k) and (i, j + 1, k), (i, j, k) and (i, j, k + 1)).

As presented above we use a two-step procedure to discretize 2.15

• Step 1 discretizes hydrodynamic evolution. The fluid parameter φ will update from φ n i,j,k to φi,j,k . Anticipating further notations we introduce sub-steps.

-Sub-step a, an acoustic sub-step that will update quantities from φ n i,j,k to φi,j,k .

-Sub-step b, a transport sub-step that will update quantities from φi,j,k to φi,j,k .

•

Step 2 deals with phase change using instantaneous relaxation. The fluid parameter φ will evolve from φi,j,k to φ n+1 i,j,k .

Step 1, hydrodynamic evolution

We begin by recalling the system that is discretized

∂ t z + u • ∇z = 0, ∂ t (ρy) + div (ρyu) = 0, ∂ t ρ + div (ρu) = 0, ∂ t (ρu) + div (ρuu + p eos isob I) = 0, ∂ t (ρE) + div ((ρE + p eos isob ) u) = 0.
(2.20)

The discretization of this system is also handled using a two-step procedure by splitting acoustic and transport phenomena leading to two systems.

Remark 5 This operator splitting procedure is similar to the discretization used in Chapter 1.

It involves separation of acoustic and transport waves. In Chapter 1 special treatment in the acoustic step was needed to have a well-balanced numerical for hydrostatic equilibrium. Here we need to pay attention to the transport equation which can lead to numerical, undesired, mixture of the two phases.

We first consider an acoustic system that reads

∂ t z = 0, ∂ t y = 0, ∂ t ρ + ρdivu = 0, ∂ t (ρu) + ρudivu + ∇p eos isob = 0, ∂ t (ρE) + ρEdivu + div (p eos isob u) = 0,
(2.21) and the transport system writes

∂ t z + u • ∇z = 0, ∂ t (ρy) + u • ∇ (ρy) = 0, ∂ t ρ + u • ∇ρ = 0, ∂ t (ρu) + u • ∇ (ρu) = 0, ∂ t (ρE) + u • ∇ (ρE) = 0.
(2.22)

In the transport system we can see that each conservative variable is advected at material velocity u. This operator splitting procedure allows to treat independently each system.

Sub-step a, acoustic

In order to discretize the system of Equations 2.21 we propose to use a relaxation for the acoustic system that is reminiscent of the Lagrangian Gas dynamics equations (see Després 2010). To do so, we first rewrite the system using the variables τ = 1 ρ , u and E. Then we proceed to a linearization of the equations and then discretize them.

It then writes

∂ t V + τ div (F (V )) = 0, V = (z, y, τ, u, E) T , F (V ) =       0 T 0 T -u T p eos isob I p eos isob u T       . (2.23)
We use a relaxation based on ideas introduced by Bouchut 2004 for the Euler system that extends the Suliciu-based relaxation (see [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF][START_REF] Chalons | Modified Suliciu relaxation system and exact resolution of isolated shock waves[END_REF][START_REF] Suliciu | On the thermodynamics of rate-type fluids and phase transitions. I. Rate-type fluids[END_REF]. This relaxation allows to remove some non-linearities due to the equation of state. The relaxation introduces two variables Π and a which are approximations of the pressure p eos isob and the Lagrangian speed of sound ρc eos isob , respectively. The augmented acoustic system writes

∂ t W + τ div (F (W )) = 0, W = z, y, τ, u, E, a, Π a 2 T , F (W ) =           0 T 0 T -u T ΠI Πu T 0 T u T          
.

(2.24)

The connection with respect to the equation of state is restored at each time step by initializing Π = p eos isob (ρ, e, z, y) and a = ρc eos isob (ρ, e, z, y). We recover the initial variable V by a projection P from the relaxation variables

P : W → V = (z, y, τ, u, E)
T .

(2.25)

This system is close to have a conservative form, up to the variable τ in front of the divergence. Then we freeze this variable at discrete time t n , which can also be interpreted as a relaxation procedure, it then becomes

∂ t W + τ n (x) div (F (W )) = 0. (2.26)
This system is then discretized using a Finite Volume method

Wi,j,k = W n i,j,k -τ n i,j,k ∆t ∆x F i+1/2,j,k -F i-1/2,j,k -τ n i,j,k ∆t ∆y F i,j+1/2,k -F i,j-1/2,k -τ n i,j,k ∆t ∆z F i,j,k+1/2 -F i,j,k-1/2 ,
(2.27) in which we have F i+1/2,j,k = F ∆ W n i,j,k , W n i+1,j,k ; e x and for e ∈ {e x , e y , e z },

F ∆ (W L , W R ; e) = (0, 0, -u * , Π * e, Π * u * , 0, u * ) T , u * = a L u L • e + a R u R • e a L + a R - Π R -Π L a L + a R , Π * = a R Π L + a L Π R a L + a R - a L a R (u R • e -u L • e) a L + a R ,
(2.28)

where the velocity waves a L , a R are given by

a L = ρ L c L + 1 2 (1 + max (γ 0 , γ 1 )) Π R -Π L ρ R c R + u L • e -u R • e + a R = ρ R c R + 1 2 (1 + max (γ 0 , γ 1 )) Π L -Π R a L + u L • e -u R • e +        if Π R > Π L , (2.29) and a R = ρ R c R + 1 2 (1 + max (γ 0 , γ 1 )) Π L -Π R ρ L c L + u L • e -u R • e + a L = ρ L c L + 1 2 (1 + max (γ 0 , γ 1 )) Π L -Π R a R + u L • e -u R • e +        if Π R ≤ Π L . (2.30)
Remark 6 This scheme falls back to a regular Suliciu-type relaxation when a L = a R = a, indeed we get

u * = 1 2 (u L • e + u R • e) - 1 2a (Π R -Π L ) , Π * = 1 2 (Π L + Π R ) - a 2 (u R • e -u L • e) .
(2.31)

Thus a low-Mach correction can also be used to reduce numerical diffusion in the pressure Π * , see Chapter 1. We have noticed that using both a low-Mach correction and a directional splitting leads to an unstable numerical scheme.

After relaxation towards equilibrium, it simply results in Ṽi,j,k = P Wi,j,k .

(2.32)

It can be noticed that the numerical scheme also writes

L i,j,k ρi,j,k Ṽi,j,k = ρ n i,j,k V n i,j,k - ∆t ∆x P F i+1/2,j,k -P F i-1/2,j,k - ∆t ∆y P F i,j+1/2,k -P F i,j-1/2,k - ∆t ∆z P F i,j,k+1/2 -P F i,j,k-1/2 , (2.33)
with

L i,j,k = 1 + ∆t ∆x u * i+1/2,j,k -u * i-1/2,j,k + ∆t ∆y u * i,j+1/2,k -u * i,j-1/2,k + ∆t ∆z u * i,j,k+1/2 -u * i,j,k-1/2 .
(2.34)

Up to the compression factor L i,j,k , the variable ρi,j,k Ṽi,j,k admits a quasi-conservative update that will be used in the next section.

Sub-step b, transport

We first recall the equation we want to discretize, for φ ∈ {z, ρy, ρ, ρu, ρE}

∂ t φ + u • ∇φ = 0. (2.35)
To take advantage of the first acoustic step, this equation is simply rewritten as

∂ t φ + div (φu) -φdivu = 0, (2.36)
and is partially discretized as follows φi,j,k = φi,j,k -

∆t ∆x φ * i+1/2,j,k u * i+1/2,j,k -φ * i-1/2,j,k u * i-1/2,j,k - ∆t ∆y φ * i,j+1/2,k u * i,j+1/2,k -φ * i,j-1/2,k u * i,j-1/2,k - ∆t ∆z φ * i,j,k+1/2 u * i,j,k+1/2 -φ * i,j,k-1/2 u * i,j,k-1/2 + ∆t ∆x φi,j,k u * i+1/2,j,k -u * i-1/2,j,k + ∆t ∆y φi,j,k u * i,j+1/2,k -u * i,j-1/2,k + ∆t ∆z φi,j,k u * i,j,k+1/2 -u * i,j,k-1/2 ,
(2.37)

for which u * has been defined in Equations 2.24. Using the factor L i,j,k , it equivalently writes

φi,j,k = L i,j,k φi,j,k - ∆t ∆x φ * i+1/2,j,k u * i+1/2,j,k -φ * i-1/2,j,k u * i-1/2,j,k - ∆t ∆y φ * i,j+1/2,k u * i,j+1/2,k -φ * i,j-1/2,k u * i,j-1/2,k - ∆t ∆z φ * i,j,k+1/2 u * i,j,k+1/2 -φ * i,j,k-1/2 u * i,j,k-1/2
(2.38)

In this way, the discretization is conservative with respect to variables ρy, ρ, ρu, ρE. It remains to define the flux φ * to have a well-defined numerical scheme.

We want to minimize numerical diffusion in the case of an interface. Different numerical techniques can be used to define the flux φ * . A first approach would be to use a standard upwind scheme. It is able to transport accurately interfaces when time step is adapted to the transport velocity and with a CFL condition exactly to 1. This constraint is too strong as we wish to treat shocks, it is thus not possible to use such a scheme. Hence we choose to use an anti-diffusive approach, originally developed in Després and Lagoutière 2001 and adapted to the five-equation model in Kokh and Lagoutière 2010. The main idea of the latter scheme is to be the most downwind possible but while ensuring stability. We start by defining fluxes relative to the volume fraction z * using standard upwind scheme relative to u * from the acoustic step

u * i+1/2,j,k = ũi,j,k if u * i+1/2,j,k > 0, ũi+1,j,k otherwise. (ρ α ) * i+1/2,j,k = (ρ α ) i,j,k if u * i+1/2,j,k > 0, (ρ α ) i+1,j,k otherwise. (e α ) * i+1/2,j,k = (ẽ α ) i,j,k if u * i+1/2,j,k > 0, (ẽ α ) i+1,j,k otherwise. (ρy) * i+1/2,j,k = (ρ 1 ) * i+1/2,j,k z * i+1/2,j,k , ρ * i+1/2,j,k = z * i+1/2,j,k (ρ 1 ) * i+1/2,j,k + (1 -z * i+1/2,j,k )(ρ 0 ) * i+1/2,j,k .
(2.39) Thus, these fluxes are consistent regarded that the flux z * i+1/2,j,k is consistent. Then it suffices to give z * i+1/2,j,k a definition. As it is detailed in Appendix 2.B and in Kokh and Lagoutière 2010, we express both consistency and stability in one direction to define this latter flux. To treat the multi-dimensional case we can either apply a directional splitting or use a mean decomposition as in Peluchon 2017. Multi-dimensional update can be rewritten as zi,j,k = 1 3 zx i,j,k + zy i,j,k + zz i,j,k , (2.40) in which zx i,j,k , zy i,j,k , zz i,j,k are one dimensional updates along the X, Y and Z directions, respectively. As an example, developing zx i,j,k gives

zx i,j,k = zx i,j,k -3 ∆t ∆x z x i+1/2,j,k u * i+1/2,j,k -z x i-1/2,j,k u * i-1/2,j,k + 3 ∆t ∆x zx i,j,k u * i+1/2,j,k -u * i-1/2,j,k .
(2.41)

Hence a sufficient condition for stability and consistency for update 2.40 is that each individual update is stable and consistent. It results that, using similar stability and consistency developments as in Appendix 2.B, we adapt the trust intervals such that, along the X-direction it gives

z * i+1/2,j,k = P Ĩi,j,k zi ⇒ ,j,k Ĩi,j,k = mi+1/2,j,k , Mi+1/2,j,k ∩ ãi+1/2,j,k , Ãi+1/2,j,k , ãi+1/2,j,k = zi ,j,k + M(i+1/2) ,j,k -zi ,j,k u * (i+1/2) ,j,k u * i+1/2,j,k - ∆x 3∆t 1 u * i+1/2,j,k , Ãi+1/2,j,k = zi ,j,k + m(i+1/2) ,j,k -zi ,j,k u * (i+1/2) ,j,k u * i+1/2,j,k - ∆x 3∆t 1 u * i+1/2,j,k . (2.42)
This final volume fraction flux is used in fluxes for other variables.

Step 2, thermodynamics relaxation

Finally the last step concerns relaxation to thermodynamics equilibrium. We recall that the equations formally write, in the limit λ → +∞,

∂ t z = λ (z eq -z) , ∂ t y = λ (y eq -y) , ∂ t ρ = 0, ∂ t (ρu) = 0, ∂ t (ρE) = 0.
(2.43) Numerical discretization is straightforward and only consists in setting

z n+1 i,j,k = (z eq ) i,j,k , y n+1 i,j,k = (y eq ) i,j,k , ρ n+1 i,j,k = ρ i,j,k , (ρu) n+1 i,j,k = ( ρu) i,j,k , (ρE) n+1 i,j,k = ρE i,j,k . (2.44)
We need to find (z eq ) i,j,k and (y eq ) i,j,k knowing ρi,j,k and ēi,j,k . As the process is local to each cell (i, j, k), in the following we omit indices to simplify notations. We solve the non-linear Equations 2.10 in pressure. Knowing p eq we can then recover all necessary thermodynamics quantities at equilibrium T eq = T sat (p eq ) , (τ α ) eq = τ α (p eq , T eq ) , (e α ) eq = τ α (p eq , T eq ) , y eq = τ -(τ 0 ) eq (τ 1 ) eq -(τ 0 ) eq , z eq = (τ 1 ) eq y eq τ .

(2.45)

For a given couple (τ , ē), if Equations 2.10 have a solution then it is a mixture state, otherwise it is a pure state. As there is no analytical expression for the function T sat , it needs to be discretized.

Saturation curve T sat (p)

Both phases have their own equation of state. At thermodynamics equilibrium they obey to g 0 (p eq , T eq ) = g 1 (p eq , T eq ) ,

By enforcing equality between the two fluids, one is able to invert the relation

∀p ∈ [p min , p max ] , find T sat such that g 0 (p, T sat ) = g 1 (p, T sat ) , (2.47)
We use this non-linear equation to build a polynomial interpolation curve

T sat (p) = n i a i (log (p)) i .
(2.48)

Coefficients are computed using the numpy.polyfit function.

Overall scheme

We summarize the numerical scheme in this section

z n+1 i,j,k = (z eq ) i,j,k , y n+1 i,j,k = (y eq ) i,j,k , ρ n+1 i,j,k = ρ n i,j,k - ∆t ∆x ρ * i+1/2,j,k u * i+1/2,j,k -ρ * i-1/2,j,k u * i-1/2,j,k - ∆t ∆y ρ * i,j+1/2,k u * i,j+1/2,k -ρ * i,j-1/2,k u * i,j-1/2,k - ∆t ∆z ρ * i,j,k+1/2 u * i,j,k+1/2 -ρ * i,j,k-1/2 u * i,j,k-1/2 , (ρu) n+1 i,j,k = (ρu) n i,j,k - ∆t ∆x ρ * i+1/2,j,k u * i+1/2,j,k u * i+1/2,j,k -ρ * i-1/2,j,k u * i-1/2,j,k u * i-1/2,j,k - ∆t ∆y ρ * i,j+1/2,k u * i,j+1/2,k u * i,j+1/2,k -ρ * i,j-1/2,k u * i,j-1/2,k u * i,j-1/2,k - ∆t ∆z ρ * i,j,k+1/2 u * i,j,k+1/2 u * i,j,k+1/2 -ρ * i,j,k-1/2 u * i,j,k-1/2 u * i,j,k-1/2 - ∆t ∆x Π * i+1/2,j,k -Π * i-1/2,j,k - ∆t ∆y Π * i,j+1/2,k -Π * i,j-1/2,k - ∆t ∆z Π * i,j,k+1/2 -Π * i,j,k-1/2 , (ρE) n+1 i,j,k = (ρE) n i,j,k - ∆t ∆x ρ * i+1/2,j,k E * i+1/2,j,k u * i+1/2,j,k -ρ * i-1/2,j,k E * i-1/2,j,k u * i-1/2,j,k - ∆t ∆y ρ * i,j+1/2,k E * i,j+1/2,k u * i,j+1/2,k -ρ * i,j-1/2,k E * i,j-1/2,k u * i,j-1/2,k - ∆t ∆z ρ * i,j,k+1/2 E * i,j,k+1/2 u * i,j,k+1/2 -ρ * i,j,k-1/2 E * i,j,k-1/2 u * i,j,k-1/2 - ∆t ∆x Π * i+1/2,j,k u * i+1/2,j,k -Π * i-1/2,j,k u * i-1/2,j,k - ∆t ∆y Π * i,j+1/2,k u * i,j+1/2,k -Π * i,j-1/2,k u * i,j-1/2,k - ∆t ∆z Π * i,j,k+1/2 u * i,j,k+1/2 -Π * i,j,k-1/2 u * i,j,k-1/2
, in which star quantities are defined in Sections 2.2.1 and 2.2.1. We notice that, as expected, the final form of the scheme is conservative for mixture density, momentum and total energy. Volume and mass fractions are at equilibrium even though they temporary evolved out of equilibrium during acoustic and transport steps.

In the case where the thermodynamics relaxation is omitted, the numerical scheme simply falls back to the discretization of the five-equation model using anti-diffusive scheme.

Numerical experiments

Global setup

We use stiffened gases as equations of state for both phases. These equations of state write

p eos α (ρ α , e α ) = (γ α -1) ρ α (e α -e 0α ) -γ α Π α , T eos α (ρ α , p α ) = p α + Π α ρ α c v α (γ α -1) , s eos α (ρ α , e α ) = c v α log e α -e 0α - Π α ρ α -c v α (γ α -1) log (ρ α ) + s 0α , c 2 α (ρ α , p α ) = γ α p α + Π α ρ α .
(2.49)

The pressure shift Π α of the stiffened gas can be interpreted as the macroscopic manifestation of attractive interactions between water molecules (see [START_REF] Métayer | The Noble-Abel Stiffened-Gas equation of state[END_REF]. This term is responsible for negative pressures, even though the speed of sound remains positive. This setup uses the values in Table 2.2 for the thermodynamics, which are taken from Hurisse 2017 (Table B.3), along with values from Table 2.1 for the interpolation of the saturation curve.

In order to avoid spurious phase transition in the bulk of the liquid droplet after the first time step we need to pay attention to initialize it under the saturation curve. The liquid droplet size is 1 mm, the density is initialized at 828.34 kg/m 3 and the pressure at 30 bar. Steam is initialized at 15.76 kg/m 3 . Both phases are initialized at 50 m s -1 .

Droplet test case, 2D

We start with 2D simulations to provide qualitative comparison with similar experiments from [START_REF] Field | The effects of target compliance on liquid drop impact[END_REF]. In this experiment, the droplet size is 10 mm, the initial velocity is at 110 m s -1 and at ambient pressure, see Figure 2.2 for the snapshots from this experiment. We realized a twodimensional simulation, see 2.1.

Even though the experiment is led in different conditions from the nuclear reactor conditions, we recover the main stages. A shock wave propagates from the point of impact towards the back of the droplet. It reflects on the interface, inside the droplet, and gets focused, we refer the reader to [START_REF] Wu | On high-speed impingement of cylindrical droplets upon solid wall considering cavitation effects[END_REF] for a detailed analysis on the wave propagation inside the droplet.

In our simulation the liquid jets at the wall are less developed because the droplet velocity is lower. Starting from time t = 3 × 10 -6 s we notice the development of cavitation regions (z = 0, 1), which do not seem to be present on the experiment 2.2, at least on such long time scale. The cavitation region is in agreement with results from [START_REF] Kyriazis | Modelling cavitation during drop impact on solid surfaces[END_REF].

Droplet test case, 3D

Ideal wall

We begin with a slippery plane wall, we consider it ideal. Figure 2.3 represents the time evolution of the maximum pressure at the position of the wall. We observe a peak pressure that develops when the droplet impacts the wall. The maximum pressure is p max ≈ 586 bar. This pressure can be compared to the analytical water hammer pressure (see [START_REF] Ghidaoui | A Review of Water Hammer Theory and Practice[END_REF], for a review) ∆p = ρ 1 c 1 ∆u.

(2.50)

In our configuration we have ρ 1 = 828.34 kg/m 3 , c 1 ≈ 1161 m s -1 and ∆u = u 1 = 50 m s -1 . It results in a water hammer pressure p hammer ≈ 30 + 481 = 511 bar. Thus the numerical measured value presents a pressure excess of 15% compared to the water hammer pressure. After the impact, the pressure globally decreases with time. However we notice oscillations right after the initial peak, they seem to be related to cells undergoing phase transition at the interface between liquid and vapor and in contact with the wall. Thus, these local peaks may not be physically relevant as such, and a mean profile should be considered.

Perturbation of the wall

In this section we explore the impact of a droplet on a perturbed wall. This non-ideal wall represents a damaged wall because of, for instance, previous droplet impacts that have formed cavities. We choose to perturbe the wall using 2D Gaussian obstacles, the wall position then writes

x wall (y, z) = x 0 + i=N i=1 a i e -1
2σ 2 i ((y-yi) 2 +(z-zi) 2 ) .

(2.51) Coefficients a i , σ i , y i and z i are randomly chosen and represent respectively the amplitude, the width and the position each Gaussian perturbation. Figure 2.5 represents the wall used for the following simulations.

Numerical results

As part of a Grand Challenge on the AMD partition of the Joliot-Curie supercomputer at TGCC, Saclay, we have realized two simulations at resolutions 1024 3 and 2048 3 , for a total of 13 million CPU hours. Figure 2.4 shows isopressure surfaces after the droplet has impacted the wall. We can see pressure waves in the liquid droplet propagating from the wall towards the droplet's end. Moreover we see a peak pressure which is developed inside a local cavity. Interested in the pressure developed near the wall, we have extracted the pressure field on the first cell layer after the wall. Thus the maximum pressure evolution at the wall is shown in Figure 2.6 for the two simulations. Both time series exhibit high frequencies in pressure. They can be related to a threshold effect due to the phase transition and the first order discretization of the wall. Indeed, there is no hysteresis in the phase transition, thus if a fluid particle is in a state close to saturation it can oscillate between the two phases and produce pressure oscillations. As we can see on Figure 2.6, the simulation at resolution 2048 3 does not have higher pressures than simulation at lower resolution. In term of peak pressure we can assume that convergence is achieved.

On both figures, we distinguish two stages. In a first stage we see an increase in pressure up until the droplet impacts the wall, represented by the black vertical line. After the impact we see a global decrease of the maximum pressure to which is superposed rapid peaks. These peaks are due to the inhomogeneities in the wall. Early peaks even surpass the initial impact pressure. This may be due to small but rapid jets that impact the wall on the perturbations.

Cavitation

Figure 2.8 shows the liquid droplet after the impact onto the wall. We can see in the white area inside the droplet, vapor formation, void fraction near 0.01, due to cavitation, that is to say a pressure drop which leads to formation of vapor. This cavitation zone appears after the reflection of the shock on the back of the droplet.

Cavitation refers to the creation of vapor bubbles inside liquids. This phase transition occurs due to a pressure drop in the liquid, a pressure lower than the saturation pressure p sat (T ). The standard picture follows with a vapor bubble which collapses usually in an asymmetric way resulting in high velocity and localized liquid jets. As a result strong pressure can develop.

Following the work of Hantke and Thein 2019, the previous description of cavitation can be referred as strong cavitation because it produces pure vapor phase. It is opposed to weak cavitation which produces wet steam. As it is explained in Hantke and Thein 2019, an adiabatic expansion in equilibrium models can only produce wet steam. For that reason equilibrium models cannot capture correctly cavitation effects from nucleation to the bubble collapse. Only the final phase of collapsing can be captured.

Discussion and interpretation of models

Even though the driving physical test formulation is simple, a liquid droplet impacting a wall, we have experienced multiple modeling and interpretation difficulties. The two-phase flow literature is rich of many models. Many of them consider different types of fractions, either volume, mass, energy. However their interpretation is rather difficult as they seem to introduce cut-off length scale or time scale. This avoids to capture phenomena such as nucleation.

We also chose to represent vapor-liquid as two phases with two different equations of state. In the initial situation this makes sense as both phases are under the critical point. However when the liquid droplet impacts the wall, pressure rises sufficiently high to reach the critical region, i.e. above the critical pressure. In this region the fluid is neither steam neither liquid and is usually referred as critical fluid. This leads to an interpretation difficulty using a color function as it either represents a liquid or vapor.

As the droplet impact leads to high pressure, mainly following the water hammer pressure, the liquid phase cannot be considered anymore as incompressible. In the compressible model, pressure is then evaluated using the couple density and internal energy. This results in a stiff formulation of thermodynamics compared to the formulation couple pressure-temperature.

The modeling of phase change, in particular driven by pressure changes, is still a challenging question. Contrary to temperature driven phase transition which can be treated in a sharp interface method as a source term at the interface, cavitation may appear in the bulk of the liquid leading to out of equilibrium liquid which may translate to temporary negative pressures.

Finally negative pressures, also called tension, is the ability of liquids to support stretching. Contrary to gases, in the liquid phase inter-molecular interactions are important. From statistical mechanics it can be shown that liquid pressure is composed of two terms: a thermal pressure due to shock interactions, as in gases, and a term representing short range attractive interactions which is negative. Experiments, see Caupin and Herbert 2006 for a review, and molecular dynamics simulations, for example see [START_REF] Tan | Molecular dynamics study of a simple liquid at negative pressures[END_REF] show that in wellprepared conditions a liquid may resist to stretching instead of undergoing phase transition. In a sufficiently strong stretching, water liquid "breaks" and steam is formed: this is called homogeneous cavitation. In the case of stiffened gases they are known to lead to negative pressures while keeping a positive speed of sound. The interpretation of negative pressures appearing in such a dynamical context is not clear as it is rather different from experiments. Further modeling along with numerical simulations is thus needed to properly establish if negative pressures arise for physical or numerical reasons.

Conclusion

As a conclusion we have used an Homogeneous Equilibrium Model to explore liquid droplet impacts onto a wall. We have used successive operator splittings to discretize the equations. A first splitting decouples the flow dynamic from phase transition. As a result, the flow dynamic can be temporarily out of thermodynamics equilibrium. Then a second splitting allows to decouple acoustic waves from transport waves. Each splitting stage allows to use specific numerical algorithms. In particular, the acoustic step is solved using a solver allowing large density ratios to avoid strong stability constraints on the time step. The transport system is solved using an anti-diffusive solver in order to keep a sharp interface.

As a result we have performed high resolution three-dimensional simulations as part of a Grand Challenge. We have been able to reproduce in a full 3D simulation the peak pressure at the impact onto the wall, compatible with the water hammer pressure. Adding roughness to the wall exhibits secondary peak pressures even higher than the initial peak.

Considering this model, a parametric study, for different initial speed, droplet radius would be easier using AMR. Indeed AMR would be of great value to help in reducing computation in the gas phase and in sharpening pressure gradients in the liquid droplet. However as it has been emphasized, AMR would not help in capturing cavitation zones in the strong sense as nucleation does not seem to be part of the model, even though they may also lead to strong pressure peaks. It can be argued that even though the cavitation mechanism is not well captured, the region where it appears is far from the wall. However we have assumed homogeneous liquid droplet at initial time, if the droplet presents initial steam cavities, this could also lead to high pressures when they collapse.

Because the interface is only advected with the flow velocity u, the interface model 2.52 equivalently reads in Ω ∂ t z + u • ∇z = 0,

∂ t ρ + div (ρu) = 0, ∂ t (ρu) + div (ρuu + p eos I) = 0, ∂ t (ρE) + div ((ρE + p eos ) u) = 0, ρE = ρe + 1 2 ρu 2 ,
(2.56) in which the pressure takes the form p eos (ρ, e, z) = (1 -z)p eos 0 (ρ, e) + zp eos 1 (ρ, e) .

(2.57)

Let us emphasize that there is no new assumption in Equation 2.57, only one equation of state is used at a time, either p eos 0 or p eos 1 . Indeed if the initialization of the color function z lies in the ensemble {0, 1}, then we have z (t, x) ∈ {0, 1} due to the properties of the advection equation. From a numerical point of view, it is standard and straightforward to use a grid-based method to discretize the interface model 2.56. As a drawback, such a discretization does not respect the condition z ∈ {0, 1} near the interface due to the numerical diffusion. Even though this numerical diffusion can be made small using adapted techniques such as an anti-diffusive scheme it cannot be entirely removed. As a consequence a region in which z / ∈ {0, 1} needs a proper definition for the quantities required by the numerical scheme such as the pressure and the speed of sound. We then choose to model these regions using a mixture.

Let us now outline some reasonable properties that a candidate mixture model should fulfill for approximating a sharp interface flow. Let ω be a small volume which contains a portion of the interface, then both phases are present in this volume. For a sufficiently small volume ω and based on the advection properties of the interface, the mixture should satisfy the following properties

• the components of the mixture are immiscible in the sense of the Section 2.1.1,

• the components of the mixture have the same pressure.

These properties motivate the usage of the five-equation model along with the isobaric equation of state in the modelling of such two-phase flows, see Section 2.1.5. Indeed for a sufficiently thin isobaric mixture zone it can be interpreted as an interface, with the volume fraction tending to the color function. The isobaric closure allows any temperature jump as in the interface model. However the five-equation model imposes that the mixture components are in mechanical equilibrium. The thin mixture zone then locally enforces adherence of the two phases.

2.B Anti-diffusive flux

In this appendix we recall derivation of the anti-diffusive flux in the 1D context, see [START_REF] Després | Contact Discontinuity Capturing Schemes for Linear Advection and Compressible Gas Dynamics[END_REF]; Kokh and Lagoutière 2010 for more details. We follow notations previously introduced. The equation to be solved writes (2.58) in which u = u (t, x) is a given function. We assume to be known a set of discrete values u * i+1/2 consistent with u associated with the position x i+1/2 as defined in the Section 2.2. We further introduce min and max values of z at interfaces

∂ t z + u∂ x z = 0,
m n i+1/2 = min z n i , z n i+1 , M n i+1/2 = max z n i , z n i+1 .
(2.59)

Discretization is carried out in the equivalent semi-conservative following form

∂ t z + ∂ x (zu) -z∂ x u = 0. (2.60)
We choose a first order discretization of the form

z n+1 i = z n i - ∆t ∆x z * i+1/2 u * i+1/2 -z * i-1/2 u * i-1/2 + z n i ∆t ∆x u * i+1/2 -u * i-1/2 . (2.61)
This discretization is consistent as long as z * i+1/2 is a consistent flux. The main idea of the derivation is to have the most downwind flux with a stability constraint. This is achieved in two steps using stability and consistency criteria.

2.B.1 Flux consistency

A sufficient condition for the flux consistency is given by

z * i+1/2 ∈ m n i+1/2 , M n i+1/2 . (2.62) Indeed if (z n i , u n i ) = z n i+1 , u n i+1 = (z, u) then z * i+1/2 = z and the flux satisfies z * i+1/2 u * i+1/2 =
zu which is the definition of the flux consistency.

2.B.2 Stability

Stability is analyzed depending on local velocities. We distinguish four cases to analyze the stability depending on the sign of u * k+1/2 for k

= i -1, i, i + 1 1. u * i+1/2 > 0 and u * i-1/2 > 0 2. u * i+1/2 > 0 and u * i-1/2 < 0 3. u * i+1/2 < 0 and u * i+3/2 > 0 4. u *
i+1/2 < 0 and u * i+3/2 < 0 In the following we develop only the two first cases.

Case u *

i+1/2 > 0, u * i-1/2 > 0 In the case of u * i+1/2 > 0 and u * i-1/2 > 0, a sufficient condition for obtaining local stability within cell i is given by, see [START_REF] Després | Contact Discontinuity Capturing Schemes for Linear Advection and Compressible Gas Dynamics[END_REF],

z n+1 i ∈ m i-1/2 , M i-1/2 .
(2.63) Thus developing update z n+1 i from Equation 2.61, Equation 2.63 writes

m i-1/2 ≤ z n i - ∆t ∆x z * i+1/2 u * i+1/2 -z * i-1/2 u * i-1/2 + z n i ∆t ∆x u * i+1/2 -u * i-1/2 ≤ M i-1/2 . (2.64)
It equivalently writes

-z n i u * i+1/2 + z n i u * i-1/2 - ∆x ∆t + ∆x ∆t m i-1/2 -z * i-1/2 u * i-1/2 ≤ -z * i+1/2 u * i+1/2 -z * i+1/2 u * i+1/2 ≤ -z n i u * i+1/2 + z n i u * i-1/2 - ∆x ∆t + ∆x ∆t M i-1/2 -z * i-1/2 u * i-1/2
(2.65)

By consistency of fluxes, see 2.62 and hypothesis u * i-1/2 > 0, we have a sufficient condition for stability

-z n i + z n i -m i-1/2 u * i-1/2 u * i+1/2 - ∆x ∆t 1 u * i+1/2 ≤ -z * i+1/2 -z * i+1/2 ≤ -z n i + z n i -M i-1/2 u * i-1/2 u * i+1/2 - ∆x ∆t 1 u * i+1/2
(2.66) Using hypothesis u * i+1/2 > 0, we obtain the following result

                   z * i+1/2 ∈ a n i+1/2 , A n i+1/2 a n i+1/2 = z n i + M i-1/2 -z n i u * i-1/2 u * i+1/2 - ∆x ∆t 1 u * i+1/2 A n i+1/2 = z n i + m i-1/2 -z n i u * i-1/2 u * i+1/2 - ∆x ∆t 1 u * i+1/2 (2.67) Case u * i+1/2 > 0, u * i-1/2 < 0 In the case of u *
i+1/2 > 0 and u * i-1/2 < 0 we use a standard upwind flux given by

z * i+1/2 = z n i (2.68)
Case u * i+1/2 < 0, u * i+3/2 < 0 In the case of u * i+1/2 < 0 and u * i-1/2 < 0 we have a similar result as for the case u * i+1/2 > 0,

u * i-1/2 > 0                    z * i+1/2 ∈ a n i+1/2 , A n i+1/2 a n i+1/2 = z n i+1 + M i+3/2 -z n i+1 u * i+3/2 u * i+1/2 - ∆x ∆t 1 u * i+1/2 A n i+1/2 = z n i+1 + m i+3/2 -z n i+1 u * i+3/2 u * i+1/2 - ∆x ∆t 1 u * i+1/2 (2.69) Case u * i+1/2 < 0, u * i+3/2 > 0 In the case of u *
i+1/2 > 0 and u * i+3/2 < 0 we use a standard upwind flux given by

z * i+1/2 = z n i+1
(2.70)

2.B.3 Overall algorithm

Thus in order to have both consistency and stability, we can choose the flux in the interval I n i+1/2 defined by

I n i+1/2 = ω n i+1/2 , Ω n i+1/2 = a n i+1/2 , A n i+1/2 ∩ m n i+1/2 , M n i+1/2 .
(2.71)

One can show that I n i+1/2 = ∅ as z n i ∈ I n i+1/2 (resp. z n i+1 ∈ I n i+1/2 ) when u i+1/2 > 0 (resp. u i-1/2 ). The flux value is then chosen to be the most downwind possible value that lies within I n i+1/2 . This can be formulated as the projection of the downwind value onto the trust interval I n i+1/2 , noted P I n i+1/2 , defined by

P I n i+1/2 (x) =      ω n i+1/2 if x < ω n i+1/2 , x if ω n i+1/2 ≤ x ≤ Ω n i+1/2 , Ω n i+1/2 if Ω n i+1/2 < x.
(2.72)

We define the downwind and upwind indices relative to the velocity on the right face u * i+1/2 . We define the downwind indices as

i ⇒ = i + 1 if u * i+1/2 > 0 i else (i + 1/2) ⇒ = i + 3/2 if u * i+1/2 > 0 i -1/2 else (2.73)
and the upwind index as

i = i if u * i+1/2 > 0 i + 1 else (i + 1/2) = i -1/2 if u * i+1/2 > 0 i + 3/2 else (2.74) Thus, when u * i+1/2 u * (i+1/2)
> 0 the anti-diffusive flux writes

                               z * i+1/2 = P I n i+1/2 z n i ⇒ , m n i+1/2 = min z n i , z n i+1 , M n i+1/2 = max z n i , z n i+1 , a n i+1/2 = z n i + M n (i+1/2) -z n i u * (i+1/2) u * i+1/2 - ∆x ∆t 1 u * i+1/2 , A n i+1/2 = z n i + m n (i+1/2) -z n i u * (i+1/2) u * i+1/2 - ∆x ∆t 1 u * i+1/2 , (2.75) otherwise we set z * i+1/2 = z n i .
(2.76)

Chapter 3

Spectral Volume method

Introduction

In the field of hyperbolic conservation laws, a lot of work has been dedicated to improve precision of simulations since the developments of the first order Godunov method [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF]) and subsequent Godunov-like methods. This method has provided a general framework in which to derive upwind numerical schemes i.e. taking into account the direction of the flow. However this stability can be interpreted into larger numerical diffusion not only near discontinuities but also in smooth regions. This can lead to precision difficulties when dealing with smooth solutions. When the solution is smooth enough, the truncation error of a numerical scheme can formally be estimated at Ch p These three letters offer three options to improve precision

• decreasing grid spacing with parameter h also referred as h-convergence

• increasing order of convergence of the numerical method with parameter p also referred as p-convergence

• lowering the value of the coefficient C: this can happen by adapting the numerical scheme to better fit specific simulation cases.

When a numerical scheme has a convergence order p ≥ 2, it is qualified as high-order. A first category of robust and high-order methods is based on the reconstruction of the solution using neighbours of Finite Volume Cells. The general idea is to follow Godunov's method but using reconstructed states at faces of cells instead of cell centered values. Such a well-known method is called MUSCL-Hancock scheme which stands for "Monotonic Upstream-centered Scheme for Conservation Laws" along with a predictor-corrector for time integration, see Toro 2009.

A second category is called Spectral Methods, inspired by Finite Elements, in which one adds degrees of freedom inside cells in order to reconstruct a polynomial, such cells are then called Spectral Cells. Because of the polynomial reconstruction, the solution is smooth inside a Spectral Cell and this can be exploited to derive a numerical scheme. Again at interfaces between Spectral Cells, the solution is discontinuous and a Riemann Solver is used as in the original Godunov's method. There are essentially three families of Spectral Methods: Discontinuous Galerkin method, Spectral Difference method and Spectral Volume method. We can distinguish them based on the equation formulation used to derive the numerical scheme; strong, integral and weak formulation. In the case of the Discontinuous Galerkin method is based on a weak formulation of the conservation laws and leads to different numerical schemes depending on the definition of the test function space, for example see [START_REF] Cockburn | Discontinuous Galerkin Methods: Theory, Computation and Applications[END_REF]. It has recently gained popularity because of its ease to adapt to different types of equations and meshes. Spectral Difference method relies on the reconstruction of both smooth solutions and flux functions. Evolution of degrees of freedom is then based on a strong formulation of the equations and finite differences [START_REF] Liang | Spectral difference method for compressible flow on unstructured grids with mixed elements[END_REF]. Its strength is an easy formulation. Finally Spectral Volume method has been introduced in a serie of papers [START_REF] Wang | Spectral (Finite) Volume Method for Conservation Laws on Unstructured Grids. Basic Formulation[END_REF][START_REF] Wang | Spectral (Finite) Volume Method for Conservation Laws on Unstructured Grids: II. Extension to Two-Dimensional Scalar Equation[END_REF]Wang and Liu 2004;Wang et al. 2004;[START_REF] Liu | Spectral (finite) volume method for conservation laws on unstructured grids V: Extension to three-dimensional systems[END_REF] and is close to the original Finite Volume method. Indeed it is based on the reconstruction of smooth solutions but degrees of freedom are evolved using the integral formulation of the equations. In the latter, the conservation property is not only expressed at the Spectral Cell level but also for degrees of freedom. All of them are here presented as spatial reconstructions and can be used along with a Runge-Kutta time integrator.

Because of the presence of discontinuities in conservation laws, a limitation process has to be used on the reconstructed solution. They play an important role in high-order schemes in order to avoid spurious oscillations near discontinuities by removing artificial extrema generated by the reconstruction step. Limiters force the numerical scheme to locally fall back to first order in accordance with the local regularity of the solution. As a consequence a too sensitive limiter forces the numerical scheme to behave like a first order scheme. Examples of well-known limiters are minmod, gminmod, superbee, TVD, . . . . This limiting process can be referred as a priori in comparison to an a posteriori method like MOOD which stands for "Multi-dimensional Optimal Order Detection", see [START_REF] Clain | A high-order finite volume method for systems of conservation laws-Multi-dimensional Optimal Order Detection (MOOD)[END_REF]. In such approach, different numerical schemes, which may be of different order of convergence are tried for a given cell update and the best solution is selected based on predefined criteria.

Another aspect is the computational cost of a numerical scheme. High-order methods usually become more expensive for a given mesh because they involve more computation per cell. However an other point of view consists in evaluating efficiency by comparing computational costs at a given error. In that sense it has been successfully shown that high-order methods are more efficient than a second order MUSCL-Hancock method on smooth examples (see [START_REF] Schaal | Astrophysical hydrodynamics with a high-order discontinuous Galerkin scheme and adaptive mesh refinement[END_REF]. Different arguments have been invoked

• more efficient p-convergence over h-convergence on smooth test cases leading to coarser grid hence less time iterations,

• less expensive non-linear computation such as Riemann problems and limiters,

• easier data locality coming from Spectral Cells making them good candidates to target new supercomputers in which cache plays an important role.

In the case of discontinuous solutions, p-convergence is however not suitable. Indeed discontinuous solutions benefit from h-convergence rather than p-convergence. Moreover as it has been said, computation cost per cell is more expensive when using a high-order scheme. This usually results in using an expensive high-order numerical scheme in discontinuous regions. Hence we propose to solve this problem by using cheaper standard Finite Volume scheme near discontinuities and accelerate computation in smooth regions by taking advantage of p-convergence. Acceleration can be achieved by building Spectral Cells from Finite Volume Cells resulting in less cells to update in smooth regions. The Spectral Volume method is used as it provides an easy transition between Finite Volume Cells and its degrees of freedom. The chapter is organized as follows, we start by defining some notation in Section 3.1. Section 3.2 recalls the main idea of the Spectral Volume Method using a polynomial flux reconstruction for Cartesian meshes to improve efficiency. Then in Section 3.3 we present our numerical method. In Section 3.4 we give a matrix formulation of the Spectral Volume Method. In Section 3.5 we present numerical experiments for common test cases found in the literature.

Notations

We begin by defining some notations. Let

Ω = [x min , x max ] × [y min , y max ] × [z min , z max ]
be a Cartesian subdomain of R 3 . We want to solve an Initial-Boundary Value Problem (IBVP) in the case of a hyperbolic scalar conservation law [START_REF] Godlewski | Hyperbolic systems of conservation laws[END_REF][START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws[END_REF][START_REF] Dafermos | Hyperbolic Conservation Laws in Continuum Physics[END_REF], that takes the form

u(0, x) = u 0 (x), ∀t > 0, ∀x ∈ Ω ∂ ∂t u(t, x) + div (f (u(t, x))) = 0, B.C. on ∂Ω, (3.1)
in which div (.) is the 3D spatial divergence operator and f = (f x , f y , f z ) is called the flux function. Equation 3.1 also admits an integral formulation. Let u be a solution of Equation 3.1. By integrating the Equation 3.1 over a given volume ω ⊂ Ω and by applying the flux-divergence theorem the solution u satisfies

ū(t) = 1 |ω| ω u(t, x)dx, d dt ū(t) + 1 |ω| ∂ω f (u (t, x)) • n (x) dσ(x) = 0. (3.2)
This integral formulation will be used to derive a conservative numerical scheme. Conservative numerical schemes are known to be able to capture weak solutions, see [START_REF] Lax | Systems of conservation laws[END_REF]. As it will be developed in following sections, we will define two meshes M and M S over Ω in order to use two different representations of the numerical solution, see Figure 3.1.

First we define a Cartesian mesh, referred to as Finite Volume mesh noted M: it is composed of N x , N y and N z cells in each dimension and ∆x, ∆y, ∆z are the corresponding mesh space steps. In the X-direction (resp. Y and Z-direction) interfaces are located at x i+1/2 = (i + 1)∆x, (resp. y j+1/2 = (j + 1)∆y, z k+1/2 = (k + 1)∆z) with (x min , x max ) = (x -1/2 , x Nx+1/2 ) (resp. (y min , y max ) = (y -1/2 , y Ny+1/2 ) and (z min , z max ) = (z -1/2 , z Nz+1/2 )). Hence the domain Ω is partitioned into non-overlapping regular cells noted ω i,j,k with

ω i,j,k = x i-1/2 , x i+1/2 × y j-1/2 , y j+1/2 × z k-1/2 , z k+1/2 , Ω = i,j,k ω i,j,k , (3.3)
and |ω i,j,k | = ∆x∆y∆z. We define x i , y j and z k by setting

x i = 1 2 x i-1/2 + x i+1/2 , y j = 1 2 y j-1/2 + y j+1/2 , z k = 1 2 z k-1/2 + z k+1/2 . (3.4)
A second Cartesian mesh, referred to as the Spectral Mesh noted M S , is built by aggregating Finite Volume Cells within a set of Spectral Cells ω S is,js,ks that also forms a new Cartesian grid over Ω, see Figure 3.1. More precisely: a Spectral Cell ω S is,js,ks is an aggregate of Figure 3.1: Example of two 2D-grids used to represent the numerical solution. The mesh on the left is a standard Finite Volume mesh, each cell approximating the mean value of the exact solution. The mesh on the right is the corresponding spectral mesh with spectral cells highlighted by a blue border. Inside each Spectral Cell are represented Control Volumes, also called degrees of freedom, colored alternatively in grey and white. Both Spectral Cells and Control Volumes exactly match a set of Finite Volume Cells for a given pattern; [START_REF] Uvsq | Maison de la Simulation[END_REF]3,[START_REF] Uvsq | Maison de la Simulation[END_REF] in each direction in this example. If Finite Volume Cells are interpreted as degrees of freedom, we see that the spectral scheme can be interpreted as a compression of information when a spectral representation is adapted to the solution.

= x S is-1/2 , x S is+1/2 × y S js-1/2 , y S js+1/2 × z S ks-1/2 , z S ks+1/2 .
(3.5) Thus, for a Spectral Cell ω S is,js,ks ∈ M S we introduce the following set of indices

F V is,js,ks = {(i, j, k) such that ω i,j,k ⊂ ω S is,js,ks }.
This way a Spectral Cell ω S is,js,ks ∈ M S can be written as

ω S is,js,ks = (i,j,k)∈F V is,js,ks ω i,j,k .
The number of Finite Volume Cells inside each Spectral Cell is a free parameter and is discussed in Section 3.3.

Finally each Spectral Cell is composed of degrees of freedom, also called Control Volumes in the case of the Spectral Volume Method, noted ω CV l,m,n . Let us consider three sets of p ≥ 1 integers n x r , n y r and n z r following

p-1 r=0 n r = N S,loc x , p-1 r=0 n y r = N S,loc y , p-1 r=0 n z r = N S,loc z . (3.6)
We consider the sets of positions defined by

x S is-1/2 + n x r ∆x, 0 ≤ r ≤ p -1 , y S js-1/2 + n y r ∆x, 0 ≤ r ≤ p -1 , z S ks-1/2 + n z r ∆x, 0 ≤ r ≤ p -1 , (3.7)
that are respectively spanned by the strictly increasing sequences x CV l+1/2 , y CV m+1/2 and z CV n+1/2 . The Control Volumes are then defined by

ω CV l,m,n = x CV l-1/2 , x CV l+1/2 × y CV m-1/2 , y CV m+1/2 × z CV n-1/2 , z CV n+1/2 . (3.8)
We also define x CV l , y CV m and z CV n by setting

x CV l = 1 2 x CV l-1/2 + x CV l+1/2 , y CV m = 1 2 y CV m-1/2 + y CV m+1/2 , z CV n = 1 2 z CV n-1/2 + z CV n+1/2 . (3.9)
We also introduce two new sets of indices

CV is,js,ks = {(l, m, n) such that ω CV l,m,n ⊂ ω S is,js,ks }, F V l,m,n = {(i, j, k) such that ω i,j,k ⊂ ω CV l,m,n }.
In the following we identify the triplet of indices (i, j, k) (resp. (i s , j s , k s ), (l, m, n)) to the corresponding cell ω i,j,k (resp. ω S is,js,ks , ω CV l,m,n ). We define u i,j,k (resp. u is,js,ks , u l,m,n ) to be a value associated with a Finite Volume Cell (i, j, k) (resp. Spectral Cell (i s , j s , k s ), Control Volume (l, m, n)).

Let us note the difference of construction of the two meshes compared to classical high-order numerical methods. In a classical high-order method one defines a unique spectral grid first and then defines degrees of freedom inside each Spectral Cell. This difference is related to the point of view of accelerating computation using a spectral method rather than increasing precision.

For the sake of clarity we will first present our discretization strategy in the case of a scalar conservation law 3.1. Then we will present the extension to the full Euler system in Section 3.5.1 consisting in applying the same procedure to each equation.

Spectral Volume method

Semi-discrete scheme

In this section we only deal with the spectral mesh M S and forget about the initial Finite Volume mesh M. The main idea of the Spectral Volume Method is to enforce Spectral Cell's polynomial to follow a conservation principle for each Control Volume. As it has already been mentioned, in the Spectral Volume Method, degrees of freedom can be interpreted as a Finite Volume sub-grid within the Spectral Cell. This method only deals with spatial discretization and does not provide any insight into the time integration, we are free to choose one's preferred method. Therefore we shall present the method as a semi-discrete numerical scheme.

Let u S is,js,ks be a polynomial associated with the Spectral Cell ω S is,js,ks of the form u S is,js,ks

u S is,js,ks (t) q,r,s x q y r z s .

(3.10)

The number of coefficients of this polynomial matches the number of Control Volumes per Spectral Cell and is of degree p -1 per direction. At each time t, u S is,js,ks represents a polynomial approximation of the exact solution u within the Spectral Cell ω S is,js,ks . In order to derive a numerical scheme we take advantage of the following property: the polynomial u S is,js,ks is characterized by its mean values in each Control Volume ω CV l,m,n of the Spectral Cell ω S is,js,ks . This allows us to seek for a numerical scheme evolving the polynomial u S is,js,ks based on the discretization of the Equation 3.2

u l,m,n (t) = 1 |ω CV l,m,n | ω CV l,m,n u S is,js,ks (t, x)dx d dt u l,m,n (t) + 1 |ω CV l,m,n | ∂ω CV l,m,n F (t, x) • n (x) dσ(x) = 0 (3.11)
where F = (F x , F y , F z ) is a polynomial numerical flux that will be defined in the next section (for the X-direction).

Polynomial flux reconstruction F x

We need a high-order estimation of the numerical flux F x in the Equation 3.11. To this point there are two types of faces for the Control Volumes: internal faces (black lines in Figure 3.2) and external faces (blue lines in Figure 3.2) relative to the Spectral Cell ω S is,js,ks . Even though for internal faces, the numerical flux

f x u S is,js,ks (t, x) (3.12)
can be analytically computed for some simple flux functions, in the general case the flux function f x is highly non-linear so that the computation of the flux integral can be difficult. The standard choice for Spectral Methods is to use quadrature formulas for both internal and external faces. These quadrature rules are chosen with a sufficient order of accuracy preserve the order of accuracy of the Spectral Method. So for both external and internal faces, the evaluation of the flux integral falls back to provide an approximation of the flux function f x at quadrature points.

In the original work of Wang 2002, they used a quadrature formula per Control Volume. This choice had the advantage to work for any unstructured meshes in which distribution of Control Volumes is more complex. The inconvenient is that it is more expensive than other spectral methods such as Discontinuous Galerkin or Spectral Difference methods.

In our case we aim at taking advantage of the Cartesian structure of the mesh M S to reduce the number of quadrature points per Spectral Cell. To do so we consider a layered polynomial flux reconstruction using a Lagrange interpolation inside a Spectral Cell at centers of the Control Volumes, the points of coordinates y CV j , z CV k 1≤j≤p,1≤k≤p , it writes

F x (t, x, y, z) = m,n F x m,n (t, x)L m,n (y, z), (3.13) in which L m,n is a normalized Lagrange interpolation polynomial L m,n (y, z) = (j,k) =(m,n) y -y CV j y CV m -y CV j z -z CV k z CV n -z CV k , (3.14)
and interpolating the flux values F x m,n (t, x) 1≤j≤p,1≤k≤p defined as in Equation 3.15 and Equation 3.16. Figure 3.2 shows an example of such a distribution of point interpolations in the case of a third order method.

For internal faces, we take advantage of the polynomial representation u S is,js,ks to directly evaluate the flux function, we set

F x m,n t, x CV l+1/2 = f x u S is,js,ks (t, x CV l+1/2 , y CV m , z CV n ) . (3.15)
In the case of external faces, such that x CV l+1/2 = x S is+1/2 , the polynomials of the two adjacent Spectral Cells are a priori discontinuous. We then use a Riemann Solver F RS to have a unique flux at the face

F x m,n (t, x CV l+1/2 ) = F RS (u L , u R ), u L = u S is,js,ks (t, x CV l+1/2 , y CV m , z CV n ), u R = u S is+1,js,ks (t, x CV l+1/2 , y CV m , z CV n ).
(3.16)

The numerical flux F RS can be an Exact Riemann Solver [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF] or any Approximate Riemann Solver, we refer the reader to the large literature for this topic, local Lax-Friedrichs Lax 1954, HLL Harten et al. 1983, HLLC Toro and Chakraborty 1994, relaxation, Roe Roe 1981 solvers. . . From the semi-discrete numerical scheme 3.11, the interpolated flux functions (F x , F y , F z ) are then integrated over the faces of Control Volumes. These integrals can be exactly computed and reformulated using a matrix-vector product, see Section 3.4.

If this method is used as such, one needs to add a limitation procedure. In our case, no limiter is used at this stage as it will be coupled to a Finite Volume scheme using the MOOD procedure.

Figure 3.2: Example of a 2D Spectral Cell ω S is,js for a third order Spectral volume method, i.e. it has 3 Control Volumes in each direction. On the right figure the vertical green lines (resp. red lines) represent the external faces (resp. internal faces) in the X-direction where the polynomial flux F x is reconstructed. The Lagrange interpolation is realized at the center of Control Volumes represented by the black points.

Presentation of an hybrid FV-SV method

In this section we present an hybrid FV-SV scheme. We begin by defining two kinds of representation of the solution according to the local regularity of the solution. Then we define numerical schemes adapted to each representation. Finally we define how to select the representation and how to switch from one to the other.

Discussion on the solution representation

In a classical Finite Volume method, being first order or second order like in a MUSCL-Hancock reconstruction, cell values represent an integral approximation of the exact solution u

u i,j,k (t) ≈ 1 |ω i,j,k | ω i,j,k u(t, x)dx.
(3.17)

Increasing the order of accuracy of the scheme is synonym to increasing the precision for this unique value.

On the other hand, in the case of spectral schemes, cells hold high-order representation of the solution. A cell does not hold anymore only a single value but degrees of freedom akin to the Finite Element method. These degrees of freedom are used to represent solution which becomes then a point-wise approximation of the solution.

As far as we are concern, we keep the grid M defined earlier as in the Finite Volume method. Each value in the cell is an approximation of the mean value of the exact solution. We then use the second mesh M S when solution is considered smooth enough.Thus with these two grids we are able to have two types of representation for the solution. First, we have a piecewise constant representation as in a standard Finite Volume scheme which formally writes

u is,js,ks (t, x) = i,j,k∈F V is,js,ks φ i,j,k (t)1 ω i,j,k (x), (3.18)
1 ω i,j,k being the characteristic function of cell ω i,j,k and φ i,j,k (t) being coefficients at time t. This representation is composed of mean values in FVC. It is a priori discontinuous and is well suited when the solution is presenting discontinuities such as contacts or shocks. When the solution is smooth enough a spectral representation is used with a basis of functions which again formally writes on a smooth function basis (f α ) α to be specified

u is,js,ks (t, x) = α φ α is,js,ks (t)f α (x) (3.19)
φ α is,js,ks are again coefficients at time t. To be as close as possible to the Finite Volume representation we choose to use the Spectral Volume Method as a high-order method in a Spectral Cell.

As it has been recalled in Section 3.2, degrees of freedom are averaged values in sub-volumes of the Spectral Cell called Control Volumes, and these Control Volumes are cells of different sizes which also have a time evolution made with fluxes. In the case of the Spectral Volume Method reconstruction, the definition of (f α ) α falls back to the definition of a partition of Control Volumes inside a Spectral Cell. Once again, Control Volumes are chosen to follow the original Cartesian structure of the Finite Volume Grid of Ω. Such an example is given in Figure 3.1 in a 2D mesh. Comparing both Equations 3.18 and 3.19, we see that if there are less Control Volumes than Finite Volume Cells we have a compression of information in regions where the solution is smooth enough. This leads to less computational needs in smooth regions.

Choice of Control Volumes

It remains to define distribution of Control Volumes. The number of Control Volumes is determined by the order of accuracy we wish to achieve. On the one hand, sizes of Control Volumes determines the stability of the method. As it has already been explored in Wang and Liu 2004, some distributions are invalid for a stable numerical scheme. This has to be related to the socalled Runge phenomenon, see [START_REF] Boyd | Exponentially-Convergent Strategies for Defeating the Runge Phenomenon for the Approximation of Non-Periodic Functions, Part I: Single-Interval Schemes[END_REF]. It refers to the oscillations of polynomials near the end of the interpolation intervals in the case of an equidistant point distribution. A simple way to avoid these spurious oscillations is to use closer points near end of interpolation interval. On the other hand, it also determines the precision of the numerical scheme. Larger Control Volumes will lead to a coarser spectral mesh thus less computations. We see that a trade off has to be found between the precision and the cost of the computation. In practice we use a distribution close to Gauss-Lobatto points. Table 3.1 gives 1D distribution of Control Volumes in terms of Finite Volume Cells up to fourth order.

Time integration

In the case of Finite Volume Cells we wish to use the cheapest robust scheme able to deal correctly with discontinuities either shocks or contacts. A first order scheme is cheap and able to deal with shocks. However contact discontinuities are rapidly regularized. We have chosen a second order MUSCL-Hancock which is still able to deal with shocks using a slope limiter and contact discontinuities. Again a minmod slope limiter tend to be too much diffusive and we have chosen the generalized minmod limiter, see Equation 3.20 6,1) 8 4 (1,4,4,1) 10 Table 3.1: This table gives the distribution of Control Volumes in terms of Finite Volume Cells.

minmod θ (a, b) = minmod 1 2 (a + b), θa, θb , θ ∈ [1, 2] (3.20) order distribution N S,loc x = N S,loc y = N S,loc z 1 (1) 1 2 (1, 1) 2 3 (1,
which reduces to minmod when θ = 1. It has been used with θ = 2 to keep contact discontinuities as sharp as possible.

In the case of Spectral Volume Method, we have a semi-discrete scheme which writes

d dt u S = L SV ∆ u S
operator L SV ∆ being our Spectral Volume semi-discrete scheme. We then use a classical Runge-Kutta time integration of the same order of accuracy. A Strong Stability Preserving Runge-Kutta see [START_REF] Gottlieb | Strong Stability-Preserving High-Order Time Discretization Methods[END_REF], referred as SSP-RK, is usually preferred in the literature for conservation laws to avoid non-admissible states. It has been tested but did not show any major difference for our test cases.

Limitation procedure

As presented above we have two types of solution representation but only one of them can be used at a time. We need to define when to use them. Moreover we need a mapping to go from one representation to the other. We choose to use the Multi-dimensional Optimal Order Detection method, referred as MOOD see [START_REF] Clain | A high-order finite volume method for systems of conservation laws-Multi-dimensional Optimal Order Detection (MOOD)[END_REF]. The general philosophy of this method could be summarized as a "try-catch" algorithm. It means that a high-order update candidate is realized, if it doesn't satisfy some predefined criteria it is discarded and a fallback scheme is used. Criteria used are observations of NaN, negative pressure or density and a relaxed discrete maximum principle. The result may fail to satisfy these criteria because of trying to construct a polynomial when the solution is locally discontinuous. The polynomial then presents overshoots that may go into negative densities or pressures, or violates the stability condition. In our case, the fallback scheme is already determined by the numerical scheme on the Finite Volume mesh M. When the Spectral Volume scheme fails to update a Spectral Cell, the solution on the Finite Volume mesh is updated using MUSCL-Hancock scheme. The idea of using sub-grids inside a Spectral Cell is not new and has already been explored in different ways [START_REF] Dumbser | A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolic conservation laws[END_REF][START_REF] Sonntag | Efficient Parallelization of a Shock Capturing for Discontinuous Galerkin Methods using Finite Volume Sub-cells[END_REF] To our understanding they both have drawbacks regarding discontinuities. In the work of Dumbser et al. 2014, using a DG scheme at order n, a uniform Finite Volume sub-grid of 2n+1 cells is used along with operators of projection P and reconstruction R. This choice of sub-grid is made to match the time step in the DG scheme. The major inconvenient is the need for a reconstruction operator to compute the solution on a Spectral Cell. This means that if a Spectral Cell is detected two times in a row, because P • R = I even though they follow R • P = I, there is a loss of information at the sub-grid level. In the paper of Sonntag and Munz 2017 they use a nodal DG scheme at order n. This scheme is closer to what we propose. A non-uniform sub-grid with n + 1 Finite Volume Cells, centered on Gauss points, is used to apply a standard Finite Volume scheme. There is no more loss of information between two time steps. However non-uniform sub-grid does not make it adapted to shocks. Indeed, at a fixed number of degrees of freedom, when the order increases the size ratio between the smallest and largest cell also increases. In the scheme we propose, we try to solve both problems by using uniform Finite Volume grid as a basis. If the fallback scheme is used on Finite Volume Cells, it is applied to solution at time n from the Finite Volume mesh without modification. This means that up to the quality of the detection criteria, shocks should be as precise as a fallback scheme simulation on this regular grid.

Time step computation

We suppose that both Finite Volume and Spectral Volume methods are constrained by a CFL time relation of the form

∆t = C CFL h. (3.21)
where h is the space step involved with the grid of each discretization. Because both schemes can be used to advance cells between time n and n + 1, both schemes have to use the same time step.

In our case to have compatibility between spectral and Finite Volume representations, we choose to have at least one Control Volume to match exactly a Finite Volume Cell, see Section 3.3.2. This enforces that both schemes have a time step computed with the same, minimal space step.

Then the CFL factor C CFL has to be chosen so that both schemes are stable.

Algorithm summary

We finally give a summary of the different identified steps in our algorithm:

1. Initialization on the Finite Volume mesh M.

2. Time loop, while t < t max , n < n max do:

(a) Construction of the Spectral representation on the Spectral mesh M S .

(b) Apply the MOOD procedure:

i. Try a high-order time evolution of the Spectral representation. ii. Projection of the Spectral representation onto the Finite Volume mesh M. iii. Catch bad evolution of Spectral Cells. iv. Apply the MUSCL-Hancock scheme on Finite Volume Cells of detected Spectral Cells and change fluxes between detected and non-detected Spectral Cells.

Implementation details and matrix formulation

Implementation is realized using two global grids. A grid is dedicated to hold Finite Volume representation and the other grid to hold Control Volume values. We emphasize that latter grid is small compared to the Finite Volume grid. When high-order time evolution is realized, these fluxes are also stored in a global array. This way, fluxes can be further modified if a Spectral Cell is detected. If a Spectral Cell is not detected but one neighbour is, it also needs to be projected on the Finite Volume grid again because of corrected fluxes at edges. In the following we omit the time dependency in the Spectral representation and for the corresponding discrete values. By a translation, the following calculations can always be done from the Spectral Cell (i s , j s , k s ) = (0, 0, 0). Thus we write CV = CV is,js,ks and F V = F V is,js,ks . This allows to interpret triplets of indices (l, m, n) and (i, j, k) as local indices relative to this Spectral Cell.

Let us also note U S = u S q,r,s 0≤q,r,s≤p-1 (resp. U CV = (u l,m,n ) (l,m,n)∈CV and U F V = (u i,j,k ) (i,j,k)∈F V ) the vector of coefficients of this polynomial (resp. the vector of values associated to Control Volume and the vector of values associated to the Finite Volume Cells).

Regarding the flux, we only consider a face in the X-direction, located at x CV l+1/2 , as an example. We note F x l+1/2 the vector of flux integrals and

F x,pt l+1/2 = F x m,n x CV l+1/2 0≤m,n≤p-1
the vector of flux points associated with the Lagrange interpolation 3.13.

Finite Volume Cells to Control Volumes

This step is a simple projection onto a coarser grid. Finite Volume Cells that belong to the same Control Volume simply add their content

∀(l, m, n) ∈ CV, u l,m,n def = 1 |ω CV l,m,n | (i,j,k)∈F V l,m,n |ω i,j,k |u i,j,k . (3.22)

Control Volumes to Finite Volume Cells

This step is the counterpart of Finite Volume Cells to Control Volumes. For a valid Spectral Cell it consists in projecting the polynomial u S representation onto the Finite Volume Cells.

∀ (i, j, k) ∈ F V, u i,j,k def = 1 |ω i,j,k | ω i,j,k u S (x, y, z)dxdydz. (3.23)
This integral can be explicitly computed and only consists in a matrix-vector product, more details in the following.

From polynomial coefficients

We want to give a matrix formulation of Equation 3.23, i.e. to compute from the values associated with the Control Volumes to the values associated with the Finite Volume Cells,

U F V = (u i,j,k ) (i,j,k)∈F V . We have ω i,j,k u S (x, y, z)dxdydz = q,r,s u S q,r,s ω i,j,k
x q y r z s dxdydz (3.24)

= |ω i,j,k | q,r,s
M ((i,j,k),(q,r,s)) U S (q,r,s)

(3.25)
in which we define the matrix M by

M ((i,j,k),(q,r,s)) = 1 |ω i,j,k | ω i,j,k x q y r z s dxdydz (3.26)
The integral in the Equation 3.26 is only geometrical and thus only depends on the mesh. It is then straightforward that we have a matrix formulation

U F V = M U S (3.27)

From Control Volumes

We need to define how to obtain the polynomial coefficients from the Control Volumes data. We start by computing the vector of coefficients U S from the vector of Control Volumes U CV

1 |ω CV l,m,n | ω CV l,m,n u S (x, y, z)dxdydz = u l,m,n (3.28)
Developing the polynomial u S we obtain

ω CV l,m,n u S (x, y, z)dxdydz = q,r,s u S q,r,s ω CV l,m,n
x q y r z s dxdydz (3.29)

= |ω CV l,m,n | q,r,s N ((l,m,n),(q,r,s)) U S (q,r,s) , (3.30)
with

N ((l,m,n),(q,r,s)) = 1 |ω CV l,m,n | ω CV l,m,n
x q y r z s dxdydz.

(3.31)

Because there are as much coefficients as Control Volumes we can invert this system and obtain the polynomial coefficients with

U S = N -1 U CV (3.32)
If we combine both results we are able to compute the values of the Finite Volume Cells from the values of Control Volumes without having to compute the polynomial coefficients

U F V = M N -1 U CV . (3.33)
This formulation has the main advantage to be independent of the location of Spectral Cell (i s , j s , k s ) in the mesh M S . This implies that matrix M N -1 can be computed before the time loop on a Spectral Cell of reference.

Flux integral

In this section we aim at formulating flux integrals as a matrix-vector product using the flux points F x,pt l+1/2 computed at the center of Control Volumes. As it has been presented in Section 3.2.2 the integral of the polynomial flux reconstruction 3.13 writes

F x l+1/2,m,n = 1 |Γ l+1/2,m,n | Γ l+1/2,m,n F x (x CV l+1/2 , y, z)dydz, Γ l+1/2,m,n = ω CV l,m,n ∩ ω CV l+1,m,n , F x (x, y, z) = m,n F x m,n (x)L m,n (y, z), L m,n (y, z) = (j,k) =(m,n) y -y CV j y CV m -y CV j z -z CV k z CV n -z CV k .
(3.34)

To have a matrix-vector product we first find the polynomial coefficients of F x x CV l+1/2 , y, z = 0≤r,s≤p-1 c r,s y r z s by solving the following linear system

N f C = F x,pt l+1/2 , (3.35) in which N f ((m,n),(r,s)) = y CV m r z CV n
s and C r,s = c r,s . As in the previous section it results that

F x l+1/2 = M f N f -1 F x,pt l+1/2 , M f ((m,n)(r,s)) = 1 |Γ l+1/2,m,n | Γ l+1/2,m,n y r z s dydz.
(3.36)

Numerical experiments

In this section we present numerical experiments both in term of precision and efficiency. The tests are done with order ranging from 1 to 4. For a given variable φ, the L 1 error on φ is measured as

φ num -φ sol L 1 = 1 N x N y i,j |φ num i,j -φsol i,j | in which φsol i,j = 1 |C i,j | Ci,j φ(x, y)dxdy
In the following we refer to "CV perf" as the number of Control Volumes updated per second whereas "FV perf" refers to the number of Finite Volume Cells updated per second which are standard metrics in grid based algorithms. Of course one has to be careful when comparing such numbers as they depend on the algorithm, their implementation and the hardware architecture.

In our case all performance numbers are presented on a single Intel Sandy-Bridge E5-2670 processor. The first order simulation is implemented with the SVM algorithm which may lower performances due to unnecessary computations. We chose to test our numerical method on the well-known adiabatic Euler system and selected test cases showing different aspects of numerical schemes such as transport, acoustic, smooth, discontinuous or high Mach flows.

Adiabatic Euler system

The adiabatic Euler system is used to model dynamics of inviscid fluids, see [START_REF] Godlewski | Hyperbolic systems of conservation laws[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] for basic properties on this system. It is a system of conservation laws for the mass density ρ, the momentum density m = ρu and the energy density ρE

∂ ∂t ρ + div (ρu) = 0 ∂ ∂t (ρu) + div (ρu ⊗ u + pI) = 0 ∂ ∂t (ρE) + div ((ρE + p) u) = 0 ρE = ρe + 1 2 ρu 2 , p = p EOS (ρ, e) (3.37)
where u is the material velocity, p EOS is a given pressure law and e is the specific internal energy of the fluid. Bold expressions represent vectors in 2 or 3D, product u ⊗ u denotes tensor product, i.e. (u ⊗ u) i,j = u i u j and I is the identity matrix in 3D. The pressure law in the following is chosen to be a perfect gas where γ > 1 is the adiabatic index of the gas. Even though the numerical scheme in Section 3.3 is presented for a scalar conservation law it can be directly applied to a system of conservation by using the same procedure for each equation in the system.

p EOS (ρ, e) = (

Smooth transport test case

This test consists in the advection of a sinusoide density profile at constant velocity. Initial conditions are given at t = 0 by, ∀x, y ∈ [0, 1] ρ(x, y) = 2 + sin(2πx) sin( 2πy)

u(x, y) = 10, v(x, y) = 10 p(x, y) = 1
This test is used for measuring the effective order of the spectral volume method of order p when there is no limitation procedure. Indeed this test case is smooth which makes it well-suited for a convergence study. We report L 1 errors on density variable in tables 3.3, 3.4 and 3.5 for orders between 2 and 4 and for the same number of Finite Volume Cells. Effective order tends to match formal accuracy by lower values when the number of Spectral Cells increases. Figures 3.3a and 3.3b summarize the results from the tables. We notice that the plots of orders 2 and 3 cross at 160 fine cells. This can be understood by the way the numerical scheme is constructed. Assuming that the error decreases as C∆x p for a smooth solution, the number of Control Volumes specifies the order p of convergence but the number of fine cells inside a Control Volumes specifies the constant C. The more fine cells are aggregated to form a Control Volume the bigger C is. We also show performance results in the tables 3.2, 3.3, 3.4 and 3.5. In a smooth test case, with MOOD disabled, we only measure performance relative to the Spectral Volume method stage. We first notice a decrease in performance when going to higher orders. This is expected as more computation is done per Control Volume. However for a given error we can see that a high-order scheme helps in allowing to use coarser grids hence less time iterations. We also give FV performances as an indicator to the maximal performance that can be achieved when when the fallback scheme is enabled. Of course for first and second order, there is no benefit as Control Volumes exactly map Finite Volume Cells. However we see a boost of performances in the third and fourth orders.

Discontinuous transport test case

This test case is a transport test with a contact discontinuity. This test is meant to measure the smearing of the discontinuity in time and to detect spurious oscillations. Even in the case of the Godunov scheme that involves the exact solution of the Riemann problem, the contact discontinuity will be smeared by the averaging steps. We consider the initial conditions

ρ(x, y) = 1 if 0.1 < x < 0.3, 0.1 < y < 0.3, 0.1 otherwise, p(x, y) = 1, u(x, y) = 1, v(x, y) = 1.
The initialization is performed on the fine grid in order to avoid early smearing of the interface. We use periodic boundary conditions. Results are shown in Figure 3.4. All four simulations involve the same number of fine cells 320, which implies that there are 320, 160, 80, 40 Spectral Cells for orders 1, 2, 3 and 4 respectively.

Results are plotted after one advection tour. Top left quadrant of Figure 3.4a shows the result of the standard finite volume scheme at order 1. The contact discontinuity is highly smeared up to the point of loosing track of the initial shape. Other quadrants of Figures 3.4a show results at orders 2, 3 and 4 respectively. We notice that at order 3, there is a larger smearing of the interface compared to the second order case. This can be related to the convergence study in the sinusoide test case 3.5.2 and the observation of a larger error at order 3 compared to order 2 when the mesh is coarse. The more the discontinuity is smeared in time, the less cells are detected. Figure 3.4b shows the profile of density along the first diagonal at the final time for all orders.

The exact solution is also represented in dash grey line. First we notice no spurious oscillations due to the use of the fallback scheme that also involves a slope limiter. Even though SVM (2) to SVM (4) use the same fallback scheme in the discontinuous region we notice a slightly less smearing of the discontinuity.

Sod shock tube test case

The Sod test case is a Riemann problem where initial conditions are given by following two states

ρ L = 1, p L = 1, u L = 0 ρ R = 0.125, p R = 0.1, u R = 0 (3.39)
The exact solution is composed of three elementary waves, an expansion wave, a contact wave and a shock wave. As the Spectral Volume scheme is only used in smooth regions, we cannot expect much improvement issued from high order method in this test that involves two discontinuities. Figure 3.5 shows result at time t = 0.2 for a fixed Finite Volume mesh of N x = 320 and different orders. On the left panel we display density profiles and on the right panel detected bad cells. We see that detection mainly happens near discontinuities, or at a loss of regularity near the ending of the expansion wave. At order 3 and 4 we notice that a lot of Spectral Cells are detected. It is expected because we are at low resolution.

We do not observe a major improvement with methods of order higher than 2. However for the shock wave, we notice that it can be localized up to the size of a Finite Volume Cell and not to the size of a Spectral Cell. Small oscillations near the contact discontinuity are related to the slope limiter 3.20. Table 3.6: Summary of L 1 errors using SVM (1) on the density variable and corresponding effective orders obtained for the advection of the isentropic vortex.

Isentropic vortex test case

The isentropic vortex test involves the transport of a vortex by the velocity field (u 0 , v 0 ) = [START_REF] Uvsq | Maison de la Simulation[END_REF][START_REF] Uvsq | Maison de la Simulation[END_REF]. The vortex center is located at (x c , y c ) = (5, 5) in the Cartesian domain [0, 10] 2 . The solution is given by

ρ(x, y) = 1 - (γ -1)β 2 8γπ 2 exp 1 2 (1 -r 2 ) 1 γ-1 u(x, y) = u 0 - β 2π exp 1 2 1 -r 2 (y -y c ) v(x, y) = v 0 + β 2π exp 1 2 1 -r 2 (x -x c ) p(x, y) = ρ(x, y) γ (3.40)
We use periodic boundary conditions. The initialization of the Finite Volume Cells is achieved through the initialization of the spectral grid using conservative variables. It has been noticed that the initialization of the Finite Volume grid with high order quadrature can produce errors of the order of 10 -8 . This has resulted in a failure of the measurement of the order of convergence at time t = 10. However the order of convergence can be shown to be recovered on a longer time scale. In consequence, it has been decided to initialize conservative variables on the spectral grid in order to avoid this spurious error. Tables 3.6, 3.7, 3.8 and 3.9 give results for different meshes. For orders up to 3, we obtain a convergence rate close to the formal order. However in the case of SVM (4) we see a higher convergence rate when increasing number of cells that has also been observed in the case of discontinuous Galerkin method see [START_REF] Schaal | Astrophysical hydrodynamics with a high-order discontinuous Galerkin scheme and adaptive mesh refinement[END_REF].

Kelvin-Helmholtz instability test case

The Kelving-Helmholtz instability test simulates the development of an instability at the interface between two layers of different densities that undergoes a velocity shear. This density contrast can have different origins such as a gradient of temperature, a gradient of mean molecular weight, a phase difference between gas and liquid. Because this flow is unstable, a small perturbation can break the stationary fluid layering state and develop vortices. In the context of the Euler equations, the density ratio is due to a temperature gradient. is no temperature diffusion in the Euler equations, a temperature jump should be indefinitely maintained. The interface between the two layers keeps growing in time within the computational domain.

The initial situation is given by a small regularization version defined as follows, for x, y ∈ [0, 1]

ρ(x, y) = ρ 1 + (ρ 2 -ρ 1 )Θ(y; y 1 , y 2 ) u x (x, y) = u 1 x + (u 2 x -u 1 x )Θ(y; y 1 , y 2 ) u y (x, y) = u 0 y sin(2πx) Θ(y; y 1 , y 2 ) = 1 + exp 2 y -y 1 δ -1 1 + exp 2 y 2 -y δ -1 (3.41)
with parameter values displayed in Table 3.10 and periodic boundary conditions. A simulation at resolution 1280 2 was performed and the result is shown in Figure 3.7 for different times. We can see a well-defined interface in the first two snapshots that develops two spirals. In the last two snapshots, other Kelvin-Helmholtz instability modes have developed at the interface leading to the mixing of the two layers. Figure 3.8 shows the time evolution of valid Spectral Cells. As we can see, from time t = 1.2, the number of valid Spectral Cells keeps decreasing. This time corresponds to the development of new Kelvin-Helmholtz instability modes as shown in figure 3.8. It should be emphasized that no diffusion process is at stake, hence there is no limitation in the development of new Kelvin-Helmholtz modes. Such a process would define a length scale which, if resolved, could limit the use of the fallback scheme.

Mach 80 jet test case

We now turn to a more demanding test coming from the astrophysics field [START_REF] Ha | Numerical Simulation of High Mach Number Astrophysical Jets with Radiative Cooling[END_REF][START_REF] Gardner | Numerical simulation of the XZ Tauri supersonic astrophysical jet[END_REF] It consists in a high mach jet flow entering into a lower density material at rest. High mach flows are demanding for numerical schemes when solving the Euler system because of the total energy conservation. In these flows, the total energy is dominated by kinetic energy which make it difficult to be accurate for internal energy and thus may lead to negative energies. The jet enters with the following characteristics has a width of 0.1. In order to avoid to rely on a boundary related problem, we set only the first layer of Finite Volume Cells to jet values and we set Neumann boundary conditions for all faces. We observed that using only negative energy and density as MOOD detection criteria is not enough. Indeed it may lead to cells in a admissible but still non-physical state and thus to a near zero time step. We show results on Figures 3.9 and 3.10 at time t = 0.7. Figure 3.9 is a low-resolution simulation on a Finite Volume mesh 500 × 250 similar to results given in [START_REF] Ha | Numerical Simulation of High Mach Number Astrophysical Jets with Radiative Cooling[END_REF]. We notice a numerical shock instability at the head of the jet. It recalls a carbuncle instability, see Quirk 1994 which is usually observed for stationary shocks aligned to the grid. Moreover we can see that because of the low resolution and large Spectral Cells (fourth order), there is a large ratio of cells using the fallback scheme. We increased the resolution up to 3200 × 1600, see Figure 3.10 where we can see a Kelvin-Helmholtz instability developing near the head of the jet. This simulation with a grid spacing of ∆x = 6.2510 -4 was performed at a mean speed of 3.7M Finite Volume Cells update per second.The performance is no longer constant but depends on the number of detected cellsreaching 12% of total Spectral Cells at the end of the simulation, less than in the low resolution simulation. This number is similar to a sequential second order MUSCL-Hancock scheme on a regular grid.

Conclusions

In this chapter we have presented a new hybrid Spectral Volume and Finite Volume method. It aims at accelerating computation using high-order schemes in smooth regions to benefit from efficient p-convergence and standard Finite Volume schemes in discontinuous regions to benefit from standard h-convergence. We have shown results validating the order of convergence and a good behavior in presence of discontinuities. We have obtained good performance results near the performances of the second order MUSCL-Hancock scheme for third and fourth orders numerical schemes. Future works on the optimization of the scheme might lead to better performances. We also point out that a factor 2 can be easily obtained by simply switching from RK4 to RK2 without much impact on the simulation results for the KH and jet tests.

To fully take advantage of new supercomputers, numerical schemes have to present high degree of parallelism. Spectral schemes are known to ease data locality thanks to Spectral Cells. In the case of regular meshes numerical schemes usually have data parallelism. Thus domain decomposition is simple and there is no load-balancing issue. In our case the use of two schemes may break this load-balancing. Indeed if only one domain presents discontinuities, only this domain needs to use a fallback scheme, thus extra computation is needed only for some domains. when the value is 1, i.e. when the update is done with the spectral scheme. We see that the fallback scheme is essentially enabled near discontinuities.

Appendix

3.A Extension to source term

One can be interested to study hyperbolic conservation laws involving source terms, also referred to as balance laws, Equation 3.1 now takes the following form

∂ ∂t u + div (f (u)) = s(u) (3.44)
In integral form over a Spectral Cell, keeping the notations of Section 3.2, we have

u l,m,n (t) = 1 |ω CV l,m,n | ω CV l,m,n u S is,js,ks (t, x)dx, s l,m,n (t) = 1 |ω CV l,m,n | ω CV l,m,n s(u S is,js,ks (t, x))dx, d dt u l,m,n (t) + 1 |ω CV l,m,n | ∂ω CV l,m,n F (t, x) • n (x) dσ(x) = s l,m,n .
(3.45)

The source term s l,m,n is simply approximated by

s l,m,n (t) ≈ s (u l,m,n (t)) (3.46)
We use a Runge-Kutta time integrator as in the conservative case. In the MOOD strategy, if the Spectral Cell is detected, we use the MUSCL-Hancock numerical scheme on the Finite Volume Cells with a source term that is evaluated by the unique value in the cell.

3.A.1 Rayleigh-Taylor instability

We now consider the Euler system, see Section 3.5.1, with gravity

∂ ∂t ρ + div (ρu) = 0 ∂ ∂t (ρu) + div (ρu ⊗ u + pI) = ρg ∂ ∂t (ρE) + div ((ρE + p) u) = ρu • g ρE = ρe + 1 2 ρu 2 , p = p EOS (ρ, e) (3.47)
The gravity source term has the particularity of being linear in conservative variables so that the approximation in the previous section is exact. Chapter 4

High Performance Computing tools

In the previous chapter, we have explored acceleration by using high order methods in order to adapt the representation of the solution depending on the local smoothness of the solution. In this chapter we continue by exploring other ways of accelerating the time to solution. We mainly explore two approaches: software and hardware acceleration by exploiting parallelism, implicit algorithms, Adaptive Mesh Refinements and a mix of them. We begin by presenting the performance portability problematic, the Kokkos and Trilinos libraries. Then we present the different parallelism strategies explored in the codes ARK and ARK 2 used in previous chapters using Kokkos and MPI. Then we explore acceleration by using the acoustic implicit solver in low-Mach regimes and by using the Trilinos library. We explore AMR which allows to have a sparse representation of the numerical solution on the splitted numerical scheme. Finally we explore the porting to GPU of an implicit solver for radiative transfer from RAMSES using Trilinos.

The performance portability problematic

When going to large simulations one needs to use dedicated computers called supercomputers. By definition a supercomputer is a computer with high-end performance capabilities at a given time. To reach best performances, these computers use the state of the art technologies in terms of processors, networks, storage system, cooling system and software stack. Supercomputers are also called clusters because they are built upon aggregating nodes, connected together by a fast network. An example of such a supercomputer is shown on Figure 4.1 which is the Joliot-Curie machine hold at TGCC, Saclay, France.

In order to run on such clusters, a popular and standardized way is to use the Message Passing Interface, MPI. This programming model, qualified as distributed memory model, allows a process to access data from an other process' memory by explicitly engaging a communication. Hence the programmer has to explicitly split the work and distribute it between processes. MPI provides essential and basic communication patterns such as one-to-one, gather, scatter, all-toall. . . Then the main difficulty is to have correct load balancing between processes. Nowadays this model is mainly used to communicate between nodes or inside a node between Non-Uniform Memory Access regions, also called NUMA regions. At the NUMA level, the shared memory model is better suited even though MPI can still be used.

Technologies evolve rapidly and more specifically processors. Moore observed this evolution for the density of transistors which has become a law named after him. To be more precise Moore's law is an economic law predicting the doubling of the density of transistors on a CPU 4.2 shows the Moore's law in orange dots. It also shows other characteristics of processors such as frequency and the number of logical cores. We can see that up until the beginning of the years 2000, computation power was increased by the increase of frequency. The same software could then run faster on the next generation of processors, this is usually referred as "free lunch", see Sutter 2013. After reaching the "heat wall", CPU manufacturers changed strategy and started to develop multi-core CPUs. Such processors are able to run at the same time two or more threads.

With the advent of multi-core processors, the Open Multi-Processing programming model has been used to exploit shared memory parallelism. In this model, work is again splitted among threads which are mapped to physical cores. However all threads share the same memory address space and thus can access any allocated data in the memory at any time without a request to an other thread. It is of the responsibility of the programmer to avoid conflicting accesses between threads. This type of conflict is called a race condition. In a multi-core processor, threads are independent and can execute different instructions on different data. The computation power of cores comes from both the usage of hierarchical memory, called caches, and complex control units able to do pipelining of instructions, executing instructions out-of-order. . . In addition to the core complexity, vector instructions allow to reduce the number of cycles necessary to realize a computation. For example some architectures are able to realize vector additions on 8 double precision floating point in the same number of cycles as the scalar equivalent operation.

Originally used for graphics computation GPUs, Graphics Processing Units, are architectures dedicated to massive parallelism. Thus, near 2010, with the release of the programming model CUDA from NVIDIA, which stands for Compute Unified Devices Architectures, GP-GPU started to be used as accelerators for scientific applications. Contrary to a processor, a GP-GPU is composed of thousands of simple cores, called CUDA cores, grouped into Streaming Multiprocessors. From the programming point of view threads are logically grouped into CUDA blocks. Each block is assigned to a Streaming Multiprocessor. The Streaming Multiprocessor is then responsible to schedule groups of 32 threads, called warps, onto the CUDA cores. Warps have the particularity to execute the same instruction on different data. If a branch divergence occurs, the warp executes both branches resulting in a loss of performance. As we can see

The Kokkos library

The Kokkos abstract machine model

The Kokkos programming model, as of version 2.9, is based on an abstract machine model representing an exascale node, see Figure 4.3. Resources from this machine can be classified into two categories: execution spaces and memory spaces. First an execution space represents resources able to do computation. From Figure 4.3 this corresponds to CPU cores (multi-core and many-core) to which are added accelerators (GPUs). Then a memory space represents resources able to store data. From Figure 4.3 different memories are available such as on-chip memory, central memory and memory on accelerators.

The Kokkos programming model allows the user to manage which execution space to use, when computation takes place and to control data transfer between memory spaces which may be costly. Kokkos is built upon existing programming models such as OpenMP or CUDA. We refer to these internal programming models as backends.

The model is based on different key concepts

• memory space: specifies where data are stored.

• memory layout: specifies the multidimensional data layout.

• execution pattern: basic parallel algorithm to execute.

• execution space: specifies where computation takes place.

• execution policy: specifies how computation is executed.

Similar to the CUDA programming model, it defines host to be the CPU which dispatches work on devices. In the CUDA programming model, the device corresponds to the GPU, in the Kokkos context it corresponds to any backend, which can be the CPU itself.

The Kokkos execution patterns

Execution patterns are the building blocks for more complex parallel algorithms. Kokkos provides three of them • Kokkos::parallel_for,

• Kokkos::parallel_reduce,

• Kokkos::parallel_scan, First the Kokkos::parallel_for is a basic loop pattern. Each loop iteration calls the functor passed as an argument. The Kokkos::parallel_reduce pattern realizes a reduction such as a summation. It allows to do some computation before the reduction. The Kokkos::parallel_scan pattern is also called a cumulative reduction such as a cumulative summation.

Similarly to CUDA, the code executed in an execution space is called a Kokkos kernel. The Kokkos kernel is a C++ functor. A functor is an object that can be called like a function thanks to the function call class operator, operator (). This operator takes usually one argument that is the index, see example 4. An execution pattern is called from the host and executed on the device. As in the CUDA programming model, when such a pattern is called we say that a kernel is launched.

Kokkos does not guarantee any order of execution but only that the amount of work, specified by the execution policy, will be executed.

The Kokkos execution space

The Kokkos execution space directly maps to the machine execution space defined earlier. It specifies where the code will be executed. There are four possibilities:

• Kokkos::Serial represents the host and code is executed sequentially,

• Kokkos::Threads represents the host and uses POSIX threads,

• Kokkos::OpenMP represents the host and uses OpenMP,

• Kokkos::Cuda represents the device GPU and uses CUDA library.

At compile time, one matches the execution spaces to the architectures available on the target node. If multiple devices are enabled at compile time, one can choose at run time the execution space to use when launching a kernel. If no execution space is provided to the execution pattern, Kokkos relies on a default ranking of available execution spaces.

For a given execution space, kernels are executed sequentially. Between execution spaces Kokkos does not guarantee any order. To enforce an explicit barrier, Kokkos::fence waits until all kernels are finished to continue.

The Kokkos memory space and memory layout

The Kokkos memory space is also close to the memory space defined in the abstract machine model. It defines where the data will be stored. There are two main possibilities:

• Kokkos::HostSpace represents the memory on the host

• Kokkos::CudaSpace, Kokkos::CudaHostPinnedSpace,Kokkos::CudaUVMSpace represent different memory spaces on the GPU.

To remedy the lack of native multi-dimensional data container in C++, Kokkos provides its own called Kokkos::View. It is a multi-dimensional container for which the memory layout can be specified at compile-time. This feature is fundamental for both vectorization on CPU and coalescent memory access on GPUs.

The Kokkos execution policies

An execution policy specify how to execute an algorithm. Kokkos provides four of them The Range policy is the simplest policy that Kokkos provides. It gives to the user an integer work index. This work index spans a range given by the user.

The MDRange policy is the multi-dimensional equivalent of the Range policy. Kokkos provides a tuple of integers which spans a Cartesian range.

In the Team policy, one exposes different levels of parallelism in an algorithm. In this policy, one does not deal with an integer as work index but with a team of Kokkos threads, which are mapped to backend OpenMP and CUDA threads. Further details are given in the next section.

Finally the Task policy provides a way to express dependencies between unit of works.

Parallelism in ARK and ARK 2

ARK, which stands for All-Regime Kokkos, is a fork from an existing code originally developed by P. Kestener at Maison de la Simulation. It is a C++ code using both MPI and Kokkos for respectively distributed and shared memory parallelism.

Shared memory parallelism using Kokkos

We begin with the shared memory parallelism and the usage of the Kokkos library. In the first version of the code ARK, we use the flat parallelism from the Range policy. The index range used corresponds to the local mesh size, including boundary layers. The parallel index is then unlinearized to test if it corresponds to an internal cell. Because it is close to a raw

CUDA implementation this makes it efficient on GPUs. On CPUs it works well for distribution of work among threads. However the unlinearization along with the test usually prevents desired auto-vectorization from the compiler. As mentioned in the introduction CPU vectorization is the ability of a CPU to use particular instructions which are able to realize operations on multiple floating points in the same number of cycles that the scalar version. When relying on autovectorization, the usage of vector instructions is a shared responsibility between the compiler and the programmer. Depending on annotations from the programmer and loop complexity, the compiler determines if vector instructions can be used and if it provides a performance gain. As Kokkos 2.9 relies on auto-vectorization in all execution policies, we follow general guidelines from Intel to achieve vectorization, see Programming Guidelines for Vectorization 2020, namely

• avoid branches,

• functions called should be inline,

• favor unit stride memory accesses.

To avoid branches, we remove the unlinearization at the inner loop level by using a two-level parallelism strategy from the three-level hierarchical parallelism from the Team policy.

Our parallelism strategy, implemented in ARK 2 is then as follows

• work assignment, see also Figure 4.4:

a Kokkos team maps to a line of cells along the X direction, -a Kokkos thread maps to a chunk of cells, -a Kokkos lane maps to a cell update,

• functions called inside the inner loops are either marked KOKKOS_INLINE_FUNCTION or KOKKOS_FORCE_INLINE_FUNCTION,

• replace small loops by multiple function calls,

• the layout is fixed to the left layout.

An example of such kernel is given by listing 4.3.1.

Hybrid parallelism using MPI and Kokkos

In the MPI standard, one has to pass raw pointers to communicate data with the corresponding data type description. We use the fact that Kokkos exposes internal raw pointers to the user. However depending on the Kokkos backend used, the raw pointer may not be handled by MPI and extra copy to the host is needed. In order to avoid unnecessary copies between the host and the device we use the decision diagram represented in Figure 4.8. This diagram helps to determine if one can pass the raw pointer from the Kokkos::View or if copy to the host is needed. We use Kokkos utilities to determine if the device space is accessible from the host. If it is not we need to copy data to the host. One exception corresponds to the case in which the device is CUDA and the MPI implementation, for example Open MPI, has the non-standard ability to handle CUDA pointers, usually referred as the CUDA-aware feature.

Parallelism using Trilinos

As it has been presented in Chapter 1, we can use implicit solvers in a low-Mach regime. As a portable implementation of linear system solvers is a work beyond this thesis we use the Trilinos library to solve the acoustic system. The advantage of using Trilinos, over other widely used libraries such as PETSc, is that it uses Kokkos at its core to achieve performance portability. Hence this minimizes changes in terms of data structures in ARK. However a major difference with Kokkos is that Trilinos also handles the distributed memory parallelism using MPI.

As it has been presented in Section 4.3.2 we use a ghost cell pattern to deal with the MPI domain decomposition. Hence one does not need an explicit global numbering of cells because stencils from explicit numerical schemes connect only a few cells between two time steps. In the case of an implicit numerical scheme, all cells are in interaction in one time step. Moreover

Trilinos vector and matrix data structures are distributed and thus require a global numbering. We take advantage of the Cartesian structure and use the canonical global coordinates of cells within the mesh. These coordinates are then linearized to have a global index.

Even though matrix values evolve at each time step, the matrix keeps a static symmetric profile due to the uniform mesh. The matrix has a maximum of 1 + 5d non-zero entries per line, d corresponding to the spatial dimension because of the pressure equation. The assembling of the matrix is distributed however there is no shared memory parallelism: the assembling of the matrix is realized sequentially from the host.

Numerical experiments

We use a GMRES iterative solver along with a Schwarz preconditioner to take into account MPI domain decomposition. Each MPI process also uses an inner preconditioner that we choose to be an incomplete factorization, ILUT. We use the Gresho vortex as a test case to control the flow regime.

Table 4.2 compares the time to solution between the explicit and implicit solvers for different Mach regimes. We see that from Mach 10 -3 and lower Mach regimes there is a great benefit in using an implicit solver which is near 10 -3 Ma -1 . We emphasize that, even though the solution is captured faster, the solver presents the same drawbacks as in the explicit case. The solution diffuses and the vortex is lost.

As we can see from table 4.3 the solver shows a different behavior in time to solution when enabling the low-Mach correction. Indeed we notice difficulties in the convergence of the linear system. The solution exhibits checkerboard modes in density which are not damped due to the removal of numerical viscosity. Thus the linear solver error reaches a plateau higher than the threshold and stops on the maximal number of iterations. However, as in the explicit case, the vortex is preserved.

In table 4.4 we vary the size of the problem to be solved and compare it to the explicit case. We see a scaling as in the explicit numerical scheme. Indeed the numerical scheme is IM-EX, hence in a problem in dimension d, if the size is multiplied by 2 the cost is multiplied by 2 d+1 , which includes the CFL condition.

Towards an Adaptive Mesh Refinement All-Regime

The motivation of using AMR algorithms is to have a dynamic sparse representation of the numerical solution. This sparse representation can be used to focus computational efforts only in some parts of the domain like discontinuities and to capture strong gradients that would be too costly on a uniform mesh. AMR algorithms usually fall into two opposite categories: cell based and patch based, see Figure 4.9. In the former approach, cells can be refined independently from each other. This is a flexible approach which allows to have very localized and deep refinement. The mesh is represented using a tree, an octree in 3D, in which AMR leaves represent one cell of the mesh. As a drawback of this data structure the looping pattern is complex. Moreover this approach tends to couple the numerical scheme to the AMR algorithm when looping over neighbours.

In the latter approach refinement occurs by patch of cells, potentially including cells that do not need refinement. This allows to have regular patches thus simplifying the looping pattern inside a patch. Patches are then nested to allow deeper refinement.

In our case we are interested in an intermediate approach, called the block based AMR, which allows to keep the flexibility of the cell based AMR and the regularity from the patch based AMR. In this method an octree is again used to represent the mesh. However, the major change is that leaves are not mesh cells anymore but regular block of cells. This work has been done in collaboration with O. Abramkina 1 .

PABLO and Kokkos

As it has been presented in the sections above, the regularity of the computation is necessary to efficiently exploit new architectures through data parallelism. The introduction of block of cells inside leaves in the block based AMR reduces the problem of parallelization to the uniform mesh case.

So far three aspects have been mentioned

• the AMR algorithm,

• the numerical scheme,

• the performance portability.

In order to decouple them, the strategy adopted is to use a dedicated library PABLO to handle AMR and the Kokkos library for performance portability. In order to avoid to manage the combination of neighbours when implementing a numerical solver, a ghost layer is added to each AMR leaf. The costly combination of possibilities is then factorized and executed once per time step. Numerical solvers can then be implemented as in the uniform mesh case. However one needs to be careful about fluxes at jumps of levels. Because of the non-linearity of fluxes, two neighbours do not compute the same flux. This problem is addressed in the next section.

Because ghost layers are added to each AMR leaf, this increases memory footprint. Thus to reduce it, it has been chosen to have two types of data structures. A first data structure contains all leaves without ghost layers and a second data structure (called group) contains only a part of leaves along with ghost layer.

All-Regime numerical scheme adaptation

As presented in Chapter 1, the numerical scheme uses an operator splitting. This extends the stencil compared to a regular first order numerical scheme. There are two alternatives of the algorithm adaptation. Either we apply the numerical scheme for one group of leaves at once, thus requiring two ghost layers, either we apply the acoustic step for all groups and then apply the transport step, thus requiring only one ghost layer. Constraints with the PABLO library lead us to choose the second option.

In order to deal with the conservation issue, we use averaging. If l and l + 1 denote the levels between the two leaves, we average cells from leaves from fine level l + 1, see Figure 4.10. This enforces ghost cells from level l to match values of cells from level l + 1. The averaging is done before and after the acoustic step. The drawback is to increase numerical dissipation but as a consequence it provides smooth transition between levels. Figure 4.11 shows the result of a Blast wave simulation generated by a high energy deposit at the center of the box, see [START_REF] Sedov | Similarity and Dimensional Methods in Mechanics[END_REF] We compare an AMR simulation using the previous strategy (left hand panel) with a MOOD simulation as in Chapter 3 (right panel). On the bottom panel we show AMR levels and valid cells for the AMR and the SVM simulations respectively. We can see that the regions where cells are detected as invalid in the MOOD procedure corresponds to the regions where the mesh is refined. • radiative transfer equations,

Adaptive Mesh Refinement and implicit

• Poisson equation for self-gravity,

• particles of matter, sinks.

RAMSES uses Adaptive Mesh Refinement. AMR implementation is based on a tree also called octree to emphasize the number of potential sub-nodes within a node. As memory is organized as a one-dimensional array, a mapping is necessary between the three-dimensional mesh and the one-dimensional array. The Peano-Hilbert space filling curve is then used to realize this mapping.

In terms of parallelization, RAMSES uses MPI. The parallelization consists in chunking the space filling curve and assigning each chunk to a MPI process. Each MPI process is then responsible to apply the different numerical schemes. As the mesh refinement occurs at runtime, the initial work distribution may become unbalanced between MPI processes. Hence a load balancing procedure is applied consisting in a new chunking of the space filling curve.

Radiative transfer module

The radiative transfer module is used to model the propagation of radiative energy and its interaction with matter. It exists different moment models for the radiative transfer depending on different approximations. In this work we are interested in the so-called grey Flux Limited Diffusion, also known as FLD model. As the name suggests, the equation on the radiative energy E r is similar to a diffusion equation with source terms, it writes

∂ t E r -div cλ κ R ρ ∇E r = κ P ρc aT 4 -E r (4.1)
in which a is the Stefan-Boltzmann constant, c is the light speed, κ R and κ P are given opacity functions that solely depend on matter temperature T , ρ the density of matter and finally λ a flux limiter. The model is qualified as grey because the energy E r takes into account all frequencies at once. An implicit time discretization is used in order to avoid severe stability condition restriction when coupled to hydrodynamics due to the parabolic nature of Equation 4.1. The numerical scheme has the particularity to solve the radiative transfer equations per level of the AMR tree: one level uses other levels as boundary conditions, see [START_REF] Commerçon | A fast, robust, and simple implicit method for adaptive time-stepping on adaptive mesh-refinement grids[END_REF] along with boundary conditions, either physical or virtual between two AMR levels. Dirichlet virtual boundary conditions are used at interface between two levels. Moreover the radiative energy gradient is modified to account for the level difference, see [START_REF] Commerçon | A fast, robust, and simple implicit method for adaptive time-stepping on adaptive mesh-refinement grids[END_REF]. We note that this numerical scheme can have a different interpretation, it can be seen as the first iteration in a Schwarz domain decomposition. We recall that the non-overlapping Schwarz method consists in solving a global boundary value problem by solving smaller boundary value problems iteratively using appropriate virtual boundary conditions. If we consider AMR levels as the subdomains, we recover the above numerical scheme at first iteration.

RAMSES uses its own implementation of a BiCGSTAB solver along with a diagonal Jacobi preconditioner. A standard approach in such a code is to use a matrix-free implementation of the BiCGSTAB. However in the case of the radiative transfer module, the matrix is stored to avoid to compute too many times the costly coefficients.

Porting to GPU

Given the current state of the RAMSES code, it is not possible to run on GPUs. There exists different approaches to port existing Fortran codes to GPUs such OpenMP or OpenACC. However these approaches would require modifications of the code and possible loop reordering.

As it has been presented in Section 4.4, Trilinos implements linear solvers and preconditioners which can run on different architectures. Hence the objective of this work is to port the radiative transfer solver module to GPUs using Trilinos capabilities. This work has been done in collaboration with H. Bloch2 and S. Donfack3 .

We try to minimize modifications on the Fortran code by taking advantage of the assembling of the matrix. Matrix elements are computed per matrix line on the Fortran side and sent to the Trilinos data structure on the C++ side. To ease the assembling we use our own global numbering of AMR leaves.

Once the matrix is assembled for one AMR level, the Trilinos linear problem is prepared and solved. The solution is then copied back to RAMSES data structures.

In order to manage the interoperability between the Trilinos solver, which is C++, and the RAMSES code, which is Fortran code, we call C functions from the Fortran code, mainly from the diffusion_cg routine.

Results and discussion

Validation is realized using the "Dirac test case" from the RAMSES test suite. This test case consists in having an energy impulse at the center of the box and let it diffuse. Figure 4.12 shows strong scaling results in a three-dimensional Dirac simulation with AMR fixed at level 8. Simulations have run on the Jean-Zay cluster at IDRIS, France. Nodes are made of Intel Skylake with 40 cores accelerated with 4 NVIDIA V100 GPUs. We compare the GPU version, using Trilinos, and the original version, using only MPI. We see that 8 MPI processes are equivalent to 1 GPU. This result is due to the time spent in assembling the matrix which encompasses memory transfer between host and device. For instance in the case of using 4 GPUs, the total time spent to solve the linear systems is only of 4s which is only 12% of the total time.

Conclusion

Even though vectorization on Intel CPUs has been achieved using Kokkos hierarchical parallelism, auto-vectorization is not satisfactory as it relies both on the compiler and the respect of some vectorization rules. This approach makes it difficult to anticipate successful vectorization as it can only be checked after compilation. For instance in our case the GNU 8 compiler is not able to vectorize ARK 2. This is a known problem and the alternative to auto-vectorization is explicit vectorization by using architecture intrinsics. By definition intrinsics are language-bindings to architecture specific assembler instructions. By using them we enforce vectorization but we loose portability. As a remedy to this problem C++ wrappers have been developed relying on C++ features. The wrapper consists in defining a class which stores an intrinsic object. Operations on these objects are then vectorized thanks to internal calls to intrinsics operations. Portability between CPUs is recovered by choosing at compile time the correct architecture to use. This solution has been adopted in different libraries such as Boost.SIMD, xsimd and will be part of C++23 with std::simd. A future version of Kokkos (>3.1) should include such a solution. As of today, the Kokkos wrapper of intrinsics is also able to deal with GPUs, contrary to other known wrappers. Indeed, based on the hierarchical parallelism model the wrapper maps work to the first dimension of a CUDA block. Hence, by choosing at compile-time a vector length to be the size of warp, that is 32, this enforces a warp to execute the same instruction and to have coalescent memory access. As a drawback this implies some refactoring of the code. Scalar functions that are called in vectorized loops need to be redefined using the wrapper type.

Algorithmic acceleration using an implicit acoustic solver is promising. It revealed convergence difficulties when modifying the fluxes with the low Mach correction. The convergence is slow and shows that further work on the preconditioner is needed to be used in production.

AMR is a powerful tool to avoid a waste of computation resources in non-essential regions. We explored a simple extension to the uniform mesh formulation of the all-regime numerical scheme. There is still work to provide a well-balanced formulation on an AMR mesh. The difficulty may reside in the dynamic refinement. The initial one-dimensional profile can still be used to initialize the simulation. However the behavior of refinement on the profile may not be able to keep both equilibrium and conservation property.

Finally we explored the usage of Trilinos, an external library, to provide linear solvers in the AMR code RAMSES. Even though it successfully allowed the usage of GPU using the portability of Trilinos, at this stage the benefit is not clear. A lot of time is spent in the assembling phase mixing both computation of matrix coefficients and data transfer to the GPU.

Conclusion and perspectives

The aim of the thesis is to provide "all-regime" numerical schemes for convection simulations and liquid droplet impacts. We used schemes based on a splitting of operator between acoustic and transport phenomena.

In the Chapter 1, we have been able to adapt this splitting of operator with a gravity source term. Accouting for the source term in the acoustic system this allowed to derive two wellbalanced numerical schemes. A first version of the scheme is a sequel to the work of [START_REF] Chalons | A large timestep and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF]. A second version (more complex) allows to comply with the conservation of a new energy that includes the gravitational energy. Moreover it has been shown that the flux correction based on the Mach number, referred as low-Mach correction, allows to remove unnecessary diffusion in the low Mach regime. In order to test a generalized theory of convection a Grand Challenge simulation on the Jean-Zay supercomputer has been performed using the code ARK that implements this numerical scheme.

In the Chapter 2 we have applied this splitting of operator to a Homogeneous Equilibrium Model for simulating liquid droplet impacts. In order to keep a sharp interface between the liquid droplet and the gas, we have used an instantaneous relaxation from a five-equation model. The splitting of operator was then applied to the five-equation model allowing to use an anti-diffusive scheme for the transport step. We have performed different simulations of liquid droplet impacts including a Grand Challenge simulation on the Joliot-Curie supercomputer using a fork of ARK.

In the Chapter 3 we have explored the possibility of reducing the time to simulation using high-order numerical schemes on regular grids. By adapting the solution representation to the regularity of the solution we sought to reduce the number of expensive, non-linear computations. It has shown promising results using a sequential implementation.

Finally in the Chapter 4 we explored different strategies of accelerating the time to solution using different parallel programming models. In particular we have been able to obtain a single source code able to run efficiently on both GPUs and CPUs (with vectorization) by using the Kokkos library.

Different perspectives arise from the present work. Even though the numerical scheme for convective flows enables the use of time implicit numerical methods, performance and scalability also require an appropriate preconditioner for solving the linear systems involved in the simulation in order to obtain an efficient discretization that can be used in production, see Chapter 4.

We have worked on the building blocks of an AMR implementation of this numerical scheme without the gravity source term. It is not clear that the refinement (and derefinement) process is compatible with the well-balanced property. Future works could consist in studying if the well-balanced feature of the numerical scheme is impacted by the AMR context.

The hybrid high-order numerical scheme has shown promising results and suggests different further developments. A first strategy could be to use this numerical scheme for each subsystem of the splitting of operator. A second strategy could be to replace the MUSCL-Hancock fallback scheme by a numerical scheme focusing on contact discontinuities such as an anti-diffusive solver [START_REF] Després | Contact Discontinuity Capturing Schemes for Linear Advection and Compressible Gas Dynamics[END_REF] or the THINC method (a similar strategy is developed by M. Kucharik, R. Loubère). A last strategy would be to couple the hybrid method with AMR. The result would benefit from both the p-convergence in smooth regions and h-convergence from AMR in discontinuous regions.

Finally, a Kokkos SIMD wrapper is in development as an alternative to the auto-vectorization strategy used inside Kokkos. The usage of such a wrapper would enforce to think in terms of a chunk of cells. The vectorization would no longer rely on the ability of the compiler to find vectorizable loops, without loosing performance with GPUs. Titre: Développement de méthodes de simulation AMR "tout régime" pour la dynamique des fluides, applications en astrophysique et aux écoulements diphasiques Mots clés: hydrodynamique compressible, diphasique compressible, calcul haute performance, nombre de Mach, schéma ordre élevé, convection Résumé: Bien que performantes pour la capture des chocs, la plupart des méthodes de simulation standards ne sont pas adaptées à des régimes de Mach variés. Des méthodes numériques innovantes, utilisant des schémas de type Volumes Finis, robustes et précises uniformément selon le nombre de Mach (dites "tout régime") ont été récemment élaborées au CEA. Ces méthodes permettent de résoudre les équations de la mécanique des fluides compressibles pour capturer des chocs, mais aussi pour simuler des écoulements à très faible vitesse. Fort de ces résultats prometteurs, nous proposons dans cette thèse de mettre à l'épreuve ces nouvelles méthodes dans deux domaines applicatifs différents: les écoulements diphasiques à petit échelle et les écoulements compressibles en astrophysique. Pour ces deux domaines la simulation multi-régime est un point difficile. En effet, ces deux contextes d'applications ont pour coeur une modélisation d'écoulement compressible mais mettent en jeu des phénomènes de convection et de compressibilité à des régimes de Mach très variés. L'approche "tout régime" permettra de capturer des phénomènes très compressibles tout en gardant la précision sur les écoulements basse vitesse.

Title: Development of "all-regime" AMR simulation methods for fluid dynamics, application in astrophysics and two-phase flows Keywords: compressible hydrodynamics, compressible two-phase flow, high performance computing, Mach number, high order scheme, convection Abstract: Although classic simulation methods for compressible flow are efficient for shock capturing, they are not adapted to variable Mach regimes. Innovative methods using Finite Volume numerical schemes, robust and uniformly accurate with respect to the Mach number (so-called "all-regime"), were recently developed at CEA. These methods allow to solve the equations of compressible flows for both shocks capturing and flows involving very low material speed. Using the ground of these promising results, we propose within this thesis to challenge these new methods in two different application areas: small scale two-phase flows and compressible flows in astrophysics. For both contexts the multiregime simulation is a key issue: they both rely on a compressible flow modeling but involve convection and compressibility in highly-variable Mach regimes. The "all-regime" approach is a good candidate for capturing highly compressible phenomena while preserving the accuracy in the low speed flows.

Université Paris-Saclay Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

  gov/interior.shtml . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Formation of microstructures after 5 (left) and 10 (right) millions of liquid droplets impacts at high velocity (Luiset et al. 2013). . . . . . . . . . . . . . . . . . . . . . 1.1 Diagram representing energy transfers between energy reservoirs. . . . . . . . . . 1.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Gresho vortex simulations. L 1 error on the velocity in function of the Mach number at a fixed number of points of nx = 2048 . . . . . . . . . . . . . . . . . . 27 1.7 Gresho vortex simulations. L 1 error on the velocity in function of the spatial resolution, at a fixed Mach number of Ma = 10 -5 . . . . . . . . . . . . . . . . . . 28 1.8 Rayleigh-Taylor simulations. Figure shows snapshots of density, one in purple and

  Figure shows the velocity field in the box. The length of an arrow is scaled using the magnitude of the local velocity. The colorbar represents the vertical component of the velocity showing the direction of the flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 1.13 Simulation box of an atmosphere, composed of CO/CO 2 molecules, which is unstable to convection. The heating source term depends both on radiative energy and molecular composition; driving the convective instability. The colorbar represents the density. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.1 Time sequence of numerical Schlieren images from a 256 2 simulation droplet impact onto a plane wall. Time is in s. . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.2 Time sequence of density contrast from Field et al. 1989. The liquid droplet size is 10 mm at ambient pressure. At point F, cavitation occurs. . . . . . . . . . . . 58 2.3 The figure represents the time evolution of the maximum pressure at the wall for a simulation at coarse resolution 64 3 . The black vertical line indicates the impact of the droplet whereas the horizontal line indicates the water hammer pressure. . 59 2.4 Droplet impact at resolution 2048 3 . It is represented contour surfaces of the pressure field inside the liquid droplet. We see a peak pressure localized on a Gaussian perturbation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 2.5 This represents the perturbed wall. Color shows the height of Gaussian perturbations. Perturbations are removed near borders of the box. . . . . . . . . . . . . . 61 2.6 The figure represents the time evolution of the maximum pressure at the wall for simulations at high resolution 1024 3 and 2048 3 . The black vertical line indicates the impact of the droplet whereas the horizontal line indicates the water hammer pressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 The convergence rates of the Spectral Volume Method for the advection of a sinusoide test case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Advection of a square simulations, at constant velocity. . . . . . . . . . . . . . . . 3.5 Sod shock tube simulations at time t = 0.2 and for a fixed Finite Volume grid of N x = 320 cells. We show density profiles for different orders. The detected cells are colored in red whereas the valid, spectral, cells are colored in blue. . . . . . . 3.6 Convergence rates of the Spectral Volume Method for the isentropic vortex test case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Kelvin-Helmholtz instability test case. We show four snapshots at different times t = 1.2, t = 1.8, t = 2.4 and t f = 3. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Time evolution of the fraction of valid Spectral Cells.

4. 1

 1 Example of a supercomputer, here Joliot-Curie at TGCC, Bruyères-le-Châtel, France108 4.2 The figure represents the evolution over time of different characteristics for chips. 109 4.3 Figure represents an asbtract machine model of an exascale node. . . . . . . . . . 110 4.4 Graphical representation of the work assignment using hierarchical parallelism on a 2D mesh. A line of cells is assigned to a team of threads. Within the team, threads 0 to 3 split total work into contiguous chunks. Each lane in a thread is responsible of a cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 4.5 This figure shows difference of performance between initial range policy in ARK and team policy in ARK2 for different architectures. . . . . . . . . . . . . . . . . 116 4.6 Strong scaling test on the Intel Skylake partition of the Joliot-Curie supercomputer, TGCC. The time measured is the global execution time without IO operations.117 4.7 Weak scaling test on the NVIDIA V100 GPU partition of the Jean-Zay supercomputer, IDRIS. The time measured is the global execution time without IO operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 4.8 Decision diagram to determine whether MPI is device-aware. This decision can be made at compile time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 4.9 Different AMR types. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 4.10 Figure shows the averaging procedure to recover conservation between two levels l and l + 1. Both panels represent the same leaves, on the right panel ghost layer is added. In grey backgrounds, cells are averaged. . . . . . . . . . . . . . . . . . . 123 4.11 Blast wave simulation results. The left hand side shows an AMR simulation, with AMR levels between 5 and 8 along with 4 2 patches inside an AMR leaf. The right hand side shows a Spectral Volume simulation at order 4 with 64 2 spectral cells. The limitation is realized with MOOD using a second order MUSCL-Hancock scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.12 Strong scaling test. The time is measured inside the Fortran diffusion_cg routine.126

2. 1

 1 Parameters used for droplet simulations. These are interpolation coefficients for the saturation curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Parameters used for droplet simulations, taken from Hurisse 2017. . . . . . . . .

  3.10 Parameter values for the Kelvin-Helmholtz instability. . . . . . . . . . . . . . . . 4.1 TOP500 of June 2020 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 MPI processes = 1, Threads=8, N x × N y = 128 2 . . . . . . . . . . . . . . . . . . 4.3 MPI processes = 1, Threads=8, N x × N y = 128 2 . . . . . . . . . . . . . . . . . . 4.4 MPI processes = 1, Threads=8, Mach = 10 -4 . . . . . . . . . . . . . . . . . . .

Figure 2 :

 2 Figure 2: Formation de microstructures après 5 (gauche) et 10 (droite) millions d'impacts de gouttes d'eau à haute vitesse (Luiset et al. 2013).

Figure 1 . 1 :

 11 Figure 1.1: Diagram representing energy transfers between energy reservoirs.

Figure 1 . 3 :

 13 Figure 1.3: Weak scaling results obtained on Joliot Curie's Intel Skylake partition at TGCC. We use a hybrid MPI-OpenMP configuration in which one MPI task is bound to a socket. Simulations run for 1000 time steps and each MPI process treats 128 3 cells. We see that the efficiency reaches a plateau of 85%.

Figure 1 . 5 :

 15 Figure 1.5: Gresho vortex simulations. Snapshots of the magnitude of the velocity field at time t f = 10 -3 , for a resolution of 512 2 and for different Mach numbers. First line shows results where the low-Mach correction is disabled and second line where it is enabled. We see that without the low-Mach correction the scheme fails at simulating low-Mach flows.

Figure 1 .Figure 1 . 6 :

 116 Figure 1.6: Gresho vortex simulations. L 1 error on the velocity in function of the Mach number at a fixed number of points of nx = 2048

Figure 1

 1 Figure 1.7: Gresho vortex simulations. L 1 error on the velocity in function of the spatial resolution, at a fixed Mach number of Ma = 10 -5 .

Figure 1

 1 Figure 1.8: Rayleigh-Taylor simulations. Figure shows snapshots of density, one in purple and two in yellow at time t = 12.4 and for a resolution of 200 × 600.First line show results with the the all-regime scheme, where on the left the low-Mach correction is disabled and is enabled on the right. Second line shows results with a Godunov-type scheme, on the left it is first order, on the right it is second order using a Muscl-Hancock scheme. We see that with the low-Mach correction we recover features only present at second order for a standard Godunov-type scheme.

Figure 1

 1 Figure1.9: Rayleigh-Bénard instability simulations in 2D. Figure shows the time evolution of the mean absolute velocity for different ratios of Rayleigh number over critical Rayleigh number (see legend). Blue points show the case where the low-Mach correction is enabled and orange ones where it is disabled. We observe that when the low-Mach correction is enabled the onset of convection is closer to the expected critical Rayleigh number.

Figure 1

 1 Figure 1.10: Rayleigh-Bénard instability simulations. Snapshot of the local Mach number field and the velocity field. We see that in the strong stratification case, there is a large range of Mach, near zero at the center of rolls up to half at the upper boundary.

Figure 1 .

 1 Figure 1.11: Rayleigh-Bénard instability simulations in 3D. Figure shows the kinetic energy spectrum of the horizontal middle plane. The blue line corresponds to the scheme with low-Mach correction and the orange one without the low-Mach correction. We see more kinetic energy at all scales in the case of the low-Mach correction.

Figure 1 .

 1 Figure 1.12: Rayleigh-Bénard instability simulations in 3D. Figure shows the velocity field in the box. The length of an arrow is scaled using the magnitude of the local velocity. The colorbar represents the vertical component of the velocity showing the direction of the flow.

2. 1 . 3

 13 Thermodynamics equilibrium mixture, p eos eq Following the lines from Helluy and Seguin 2006; Barberon and Helluy 2005; Faccanoni et al. 2008, we distinguish two cases

Figure 2 . 1 :

 21 Figure 2.1: Time sequence of numerical Schlieren images from a 256 2 simulation droplet impact onto a plane wall. Time is in s.

Figure 2 Figure 2

 22 Figure 2.2: Time sequence of density contrast from Field et al. 1989. The liquid droplet size is 10 mm at ambient pressure. At point F, cavitation occurs.

Figure 2

 2 Figure 2.4: Droplet impact at resolution 2048 3 . It is represented contour surfaces of the pressure field inside the liquid droplet. We see a peak pressure localized on a Gaussian perturbation.

Figure 2

 2 Figure 2.5: This represents the perturbed wall. Color shows the height of Gaussian perturbations. Perturbations are removed near borders of the box.

Figure 2 . 6 :Figure 2 . 7 :

 2627 Figure 2.6: The figure represents the time evolution of the maximum pressure at the wall for simulations at high resolution 1024 3 and 2048 3 . The black vertical line indicates the impact of the droplet whereas the horizontal line indicates the water hammer pressure.

Figure 2

 2 Figure 2.8: Droplet impact at resolution 1024 3 . The liquid density is represented in red. The white zone at the back of the droplet represents the vapor volume fraction. Velocity field lines are displayed near the droplet.

  L 1 errors in function of the number of Spectral Cells.

  L 1 errors in function of the number of Spectral Cells.

Figure 3

 3 Figure 3.3: The convergence rates of the Spectral Volume Method for the advection of a sinusoide test case.

  Density color maps for different orders: 1, 2, 3 and 4 correspond respectively to the top left, top right, bottom left and bottom right quadrants. We can see a major improvement between the orders 1 and 2, however from the order 2 the smearing is not reduced. This is due to the common fallback MUSCL-Hancock scheme, that controls the smearing. Density profiles along the diagonal for different orders.

Figure 3 . 4 :

 34 Figure 3.4: Advection of a square simulations, at constant velocity.

Figure 3 . 5 :

 35 Figure3.5: Sod shock tube simulations at time t = 0.2 and for a fixed Finite Volume grid of N x = 320 cells. We show density profiles for different orders. The detected cells are colored in red whereas the valid, spectral, cells are colored in blue.

  L 1 errors in function of the number of Spectral Cells.

  L 1 errors in function of the number of Finite Volume Cells.

Figure 3

 3 Figure 3.6: Convergence rates of the Spectral Volume Method for the isentropic vortex test case.

Figure 3 Figure 3

 33 Figure 3.7: Kelvin-Helmholtz instability test case. We show four snapshots at different times t = 1.2, t = 1.8, t = 2.4 and t f = 3.

Figure 3 Figure 3

 33 Figure3.9: Snapshot of a jet simulation at Mach 80 with SVM (4) at a Finite Volume resolution 500 × 250. The top panel shows the density in a log scale. The bottom panel shows valid cells when the value is 1, i.e. when the update is done with the spectral scheme. We see that the fallback scheme is essentially enabled near discontinuities.

Figure 3

 3 Figure 3.11: Two snapshots of a Rayleigh-Taylor simulation at resolution 500 × 1500 and times t = 8.1 and t = 12.4. We observe dense material going down due to instability.

Figure 4 . 1 :

 41 Figure 4.1: Example of a supercomputer, here Joliot-Curie at TGCC, Bruyères-le-Châtel, France

  Figure

Figure 4 . 3 :

 43 Figure 4.3: Figure represents an asbtract machine model of an exascale node.

  template < int dim > K O K K O S _ I N L IN E_F UN CTI ON void Kernel < dim >:: operator ()( const Team & team ) const {// From given team compute coordinates of cell (0 , j , k

Figure 4 . 4 :

 44 Figure 4.4: Graphical representation of the work assignment using hierarchical parallelism on a 2D mesh. A line of cells is assigned to a team of threads. Within the team, threads 0 to 3 split total work into contiguous chunks. Each lane in a thread is responsible of a cell.

Figure 4

 4 Figure 4.5: This figure shows difference of performance between initial range policy in ARK and team policy in ARK2 for different architectures.

Figures 4 .

 4 Figures 4.6 and 4.7 show strong and weak scaling tests respectively. The weak scaling test was realized on a GPU partition of the Jean-Zay supercomputer, without the MPI CUDA-aware feature.

Figure 4

 4 Figure 4.6: Strong scaling test on the Intel Skylake partition of the Joliot-Curie supercomputer, TGCC. The time measured is the global execution time without IO operations.

Figure 4

 4 Figure 4.7: Weak scaling test on the NVIDIA V100 GPU partition of the Jean-Zay supercomputer, IDRIS. The time measured is the global execution time without IO operations

4. 6 . 1

 61 Figure 4.10: Figure shows the averaging procedure to recover conservation between two levels l and l + 1. Both panels represent the same leaves, on the right panel ghost layer is added. In grey backgrounds, cells are averaged.

Figure 4 .

 4 Figure 4.11: Blast wave simulation results. The left hand side shows an AMR simulation, with AMR levels between 5 and 8 along with 4 2 patches inside an AMR leaf. The right hand side shows a Spectral Volume simulation at order 4 with 64 2 spectral cells. The limitation is realized with MOOD using a second order MUSCL-Hancock scheme.

Figure 4 .

 4 Figure 4.12: Strong scaling test. The time is measured inside the Fortran diffusion_cg routine.

1

  Représentation graphique des différents régimes de transfert de chaleur depuis le coeur d'une étoile vers sa surface. La chaleur émise depuis le coeur de l'étoile par les réactions nucléaires est dans un premier temps diffusé par radiation, sans mouvement du gaz. Le transfert de chaleur est ensuite effectué à partir de mouvements macroscopiques du gaz vers la surface de l'étoile. Source https:

//solarscience.msfc.nasa.gov/interior.shtml . . . . . . . . . . . . . . . . . 2 Formation de microstructures après 5 (gauche) et 10 (droite) millions d'impacts de gouttes d'eau à haute vitesse (Luiset et al. 2013). . . . . . . . . . . . . . . . . 3

  • mini-symposium, CANUM, Cap d'Agde, France, Juin 2018, A well-balanced scheme for compressible flows with gravity at all Mach number, Introduction to Kokkos. De plus, les travaux du Chapitre 2 et du Chapitre 3 sont l'objet d'articles en cours de préparation.

	• mini-symposium, ICIAM, Valencia, Espagne, Juillet 2019, Compressible two-phase flow
	simulations of liquid droplet impacts,
	• webinar invité, EUROfusion, Juillet 2020,

• Thomas Padioleau et al. 2019. "A High-performance and Portable All-Mach Regime Flow Solver Code with Well-balanced Gravity. Application to Compressible Convection." The Astrophysical Journal 875, no. 2 (April): 128

Table 1 .

 1 1: Isothermal atmospheres at rest. Tableshowsfor different spatial resolutions the maximum velocity in the domain. We see that the velocity is maintained around zero up to the machine precision, thus illustrating the well-balanced property.

	Number of cells velocity (θ = 1) velocity (θ = M a) 128 2.910 -15 1.410 -13 256 8.110 -15 5.710 -13 512 1.510 -14 1.110 -12 1024 2.210 -14 2.210 -12 2048 4.710 -14 1.610 -12 4096 1.110 -13 4.010 -12
	1.4.4 Rayleigh-Taylor instability test case

The Rayleigh-Taylor test case is a two dimensional test case where two fluids of different densities are superposed and are at equilibrium. The denser one is on top. A small perturbation is introduced to break equilibrium.

The full setup is as follow, for a domain [-0.25, 0.25] × [-0.75, 0.75]:

  equation of state, see De Lorenzo et al. 2017 for this equation of state formulated in the (ρ, e) plane.

Table 2

 2 Table 2.2: Parameters used for droplet simulations, taken from Hurisse 2017.

	.1: Parameters used for droplet simulations. These are interpolation coefficients for the saturation curve.
	Phase	Gas (α = 0)	Liquid (α = 1)
	γ α c v α J K -1 kg -1 p 0α (Pa) e 0α J kg -1 s 0α J K -1 kg -1	1.135020250284312 2.803 × 10 3 0 9.991 579 270 858 164 × 10 5 -2.104 807 279 208 952 × 10 3	1.549523809523810 3.15 × 10 3 7.18 × 10 8 -1.4665 × 10 6 0

  Table3.2: Summary of L 1 errors using SVM (1) on the density variable and corresponding effective orders obtained on the advection of sinusoide.

	spectral	N x	L 1 error on ρ effective order CV perf FV perf
	40	40	2.2213 × 10 -1	-	3.7	3.7
	80	80	1.3576 × 10 -1	0.7104	4.2	4.2
	160	160 7.5613 × 10 -2	0.8444	4.2	4.2
	320	320 3.9979 × 10 -2	0.9194	4.1	4.1
	640	640 2.0566 × 10 -2	0.9590	3.9	3.9
	1280	1280 1.0432 × 10 -2	0.9793	3.9	3.9
					γ -1) ρe	(3.38)

  1. Contrary to a function, the advantage of a functor is that it can store data.

	Listing 4.1: Example of a functor
	struct Functor
	{
	K O K K O S _ INL IN E_F UN CTI ON
	void operator ()( int index ) const
	{
	// unit of work
	}
	};

Table 4 .

 4 3: MPI processes = 1, Threads=8, N x × N y = 128 2

	Mach number implicit time (s) iterations explicit time (s) acceleration
	10 -3 10 -4 10 -5	19 17 17		14 197 1970 (estimated)	0.74 11 115
	Table 4.2: MPI processes = 1, Threads=8, N x × N y = 128 2
	Mach number implicit no correction time (s) iterations	with correction time (s) iterations
	10 -3 10 -4 10 -5	19 17 17	39 74 74	200 200 200

  . Hence, at a given AMR level, the numerical discretization of Equation4.1 writes as in the uniform mesh case

	E r	n+1 i,j,k = E r	n i,j,k +	∆t ∆x		σ n i+1/2,j,k (∇E r ) n+1 i+1/2,j,k -σ n i-1/2,j,k (∇E r ) n+1 i-1/2,j,k
			+	∆t ∆y	σ n i,j+1/2,k (∇E r ) n+1 i,j+1/2,k -σ n i,j-1/2,k (∇E r )	n+1 i,j-1/2,k
			+	∆t ∆z	σ n i,j,k+1/2 (∇E r ) n+1 i,j,k+1/2 -σ n i,j,k-1/2 (∇E r ) n+1 i,j,k-1/2
			+κ P	n i,j,k ρ n i,j,k c a T n i,j,k	4 -E r	n+1 i,j,k
						σ n i+1/2,j,k =	cλ n i+1/2,j,k i+1/2,j, ρ n κ R n i+1/2,j,k
			(∇E r )	n+1 i+1/2,j,k =	1 ∆x	E r	n+1 i+1,j,k -E r	n+1 i,j,k

https://gitlab.erc-atmo.eu/erc-atmo/ark, version v1.0.0

The simulation outputs are available at http://opendata.erc-atmo.eu

Université Paris-Saclay, UVSQ, CNRS, CEA, Maison de la Simulation, 91191, Gif-sur-Yvette, France.

Université Paris-Saclay, UVSQ, CNRS, CEA, Maison de la Simulation, 91191, Gif-sur-Yvette, France. Université Paris-Saclay, CNRS, Institut du développement et des ressources en informatique scientifique, 91403, Orsay, France.

Remerciements

km

(g/cm 3 ) Figure 1.13: Simulation box of an atmosphere, composed of CO/CO 2 molecules, which is unstable to convection. The heating source term depends both on radiative energy and molecular composition; driving the convective instability. The colorbar represents the density.

Appendix

2.A Interface models and mixture models

In this section we make a connection between interface models and mixture models.

We define Ω ⊂ R 3 a domain which contains two fluids that are separated by an interface. We note Ω 0 (t) ⊂ R 3 the domain occupied by the gas and Ω 1 (t) ⊂ R 3 the domain occupied by the liquid. As the two phases are separated by an interface we have Ω 0 (t) ∩ Ω 1 (t) = ∅. We note Γ(t) = Ω 0 (t) ∩ Ω 1 (t) the location of the interface at a time t and for x ∈ Γ(t) we note n(t, x) the normal to the interface, directed towards the liquid phase. We also assume that they fill the whole domain i.e. Ω 0 (t) ∪ Ω 1 (t) = Ω. In the bulk of the phases, i.e. in Ω 0 (t) ∪ Ω 1 (t), the interface model writes ∂ t ρ + div (ρu) = 0, ∂ t (ρu) + div (ρuu + p eos I) = 0, ∂ t (ρE) + div ((ρE + p eos ) u) = 0,

(2.52) in which

(2.53)

Neglecting the surface tension and any other source term at the location of the interface, the two phases are connected with jump relations of Rankine-Hugoniot type at the interface velocity σn

(2.54)

The pressure jump involves the two equations of state. Neglecting the mass transfer, noted

between the two phases we have j = 0. Thus both the normal component of the velocity and the pressure fields are continuous across the interface. The interface velocity σ reduces to the normal component of the material velocity u • n. The interface behaves as a contact wave in a single phase Euler system, there is no constraint on the density jump and the tangential component of the velocity field.

Let us introduce a color function z defined by with values p 0 = 2.5, g y = -0.1 and A = 0.01. A velocity perturbation is added so that initial equilibrium cannot be sustained.

Results of simulation are given in Figure 3.11 at resolution 500 × 1500. We can see dense layer going down and light layer going up. As in the test case 3.5.6, a Kelvin-Helmholtz instability is developing at the interface between both layers leading to a numerical mixture at time t = 12.4. The behavior of the numerical scheme is really close to the Kelvin-Helmholtz simulation. The dynamics is focused near the interface between two large and smooth layers. The effect of the Spectral scheme is not really important in terms of accuracy gain, however it is significant in terms of computational cost as it allows to decrease the number of Riemann problems to be computed. Table 4.1: TOP500 of June 2020 maximum performance is reached when threads within a warp agree on the flow of instructions. As a result, the diversity in the main architectures of supercomputers increased as we can see from Table 4.1 which shows the first positions from the TOP500 ranking of the most powerful supercomputers in the world in June 2020. This diversity raises two related difficulties, portability and performance portability. Portability relates to the ability to run on different computers. The difficulty is due to differences in build systems, software dependencies, operating systems. . . Whereas performance portability is the ability to use efficiently different hardware architectures. The latter implies portability. As presented above, architectures are different and thus need different programming models to achieve good performances. This may lead to implement the same algorithm multiple times. In addition to the diversity of architectures, node size tends to increase. Multiple architectures can be on the same node thus requiring memory management and thus a flexible programming model to use all computation resources.

To address the aforementioned problems, a typical solution is to use OpenACC or OpenMP which avoids to refactor an existing code. The advantage is to add instructions to the compiler to generate parallel loops and kernels for GPUs. However memory management is difficult and can be a bottleneck on some applications as we will see further. We choose an other path which is to use a dedicated C++ library called Kokkos. We emphasize that it exists alternatives to This strategy mainly follows the aforementioned guidelines. Although explicit loop unrolling is not advised, we noticed better results in term of desired vectorization. The small loops are mainly coming from the way of writing kernels for multiple dimensions at a time through template parameter. The chosen alternative has been to develop a small template library for basic linear algebra operations and to manually unroll these loops. In order to check status of auto-vectorization we had to relied on vectorization reports from Intel compiler, see 4.A Results of this strategy is shown on Figure 4.5. We compare the initial strategy using range policy in ARK and the second strategy in ARK 2 for two architectures. The first architecture Intel KNL heavily relies on vectorization. We see that in the range policy case, whether vectorization is enabled or not does not impact performance. Whereas with the team policy we see a gain between 3 and 4. The second architecture is a NVIDIA V100 GPU. We don't notice loss in performance. Indeed a Kokkos thread corresponds to the first dimension of a CUDA block, this also ensures coalescent memory access on GPUs.

Distributed memory parallelism using MPI

In the distributed memory model, we use an explicit domain decomposition also referred as MPI domain decomposition. We recall that we work on a regular domain hence, we use a Cartesian decomposition. We use a one-to-one mapping between MPI processes and domains.

We then use a ghost cell pattern. Each MPI process has extra layer of cells, called ghost cells, which are used to implement physical boundary conditions or exchange data with MPI neighbours. All MPI processes advance in a synchronized way due to the computation of the time step through a reduction and data exchange with neighbours. Then processes are independent and apply the different stages of the numerical scheme.

Appendix

4.A Compilation time considerations

This section deals with different strategies of compilation when using Kokkos, however it may also apply to other header oriented libraries.

The C++ language is a strongly typed language meaning that the type of a variable is fixed at compile time. The language also offers generic programming, also referred as template metaprogramming, meaning that one can write at once multiple functions, for different types. The snippet 4.2 gives an example of a generic function returning the addition of objects of the same type class T. Of course this function can only be called if the body makes sense for the given type.

Listing 4.2: Template example template < class T > T plus ( T && t1 , T && t2 ) { return t1 + t2 ; }

The template parameter types of such generic function are then fixed when the function is called, either explicitly or implicitly by following deduction rules. Template code is automatically considered inline to avoid a violation of the one definition rule, keyword inline can be omitted.

From the compilation point of view this implies that binary code is generated only when these generic functions are called in a translation unit. This also means that if a generic function is called in multiple translation units, the compiler has to generate multiple times binary code related to this function. Because a generic function is inline, the compiler may choose to replace the function call by its body where it is called. Listing 4.3: Kokkos example // main . cpp # include " KernelA . hpp " # include " KernelB . hpp " int main ( int argc , char ** argv ) { // ... Kokkos :: parallel_for ( n , KernelA ( view1 , view2 ) ); Kokkos :: parallel_for ( n , KernelB ( view1 , view2 ) ); Kokkos :: deep_copy ( view1 , view2 ); // ...

}

We consider a first example with the snippet 4.3. In this example two kernels are executed in a parallel_for evolving view2 from view1 and then copy view2 back to view1. Three computation kernels are called, two from the user's code and one from Kokkos. To avoid a loss of performance kernel's function call operator void operator () (/*...*/) const is inline and thus defined in a header file. From this example, we notice some drawbacks of generic code about compilation

• any modification to KernelA body has the side effect to trigger recompilation of unchanged KernelB body,

• if the KernelA is called in an other source file, with the same call to parallel_for its body is recompiled,

• increases size of a translation thus exposing less parallelism,

• chaining kernels may reach rapidly inline limits of compilers, for example Intel compiler,

• analysis of inlining and vectorization of functors KernelA, KernelB and deep_copy are all placed in the same report. Hence it results in long and obscure optimization reports when they are generated on a per translation unit basis.

Hence we have tried to reduce these side effects in ARK. We consider then two types of inline code

• critical inline functions for optimization from the compiler such as auto-vectorization,

• kernel inline functions that should not be optimized in the context where they are called such as Kokkos kernels.

Thus we consider that optimization from the compiler should be dedicated to a given kernel. This prevents for example the compiler to do optimizations such as loop fusion between kernels that may benefit from cache reuse on CPU.

4.A.1 Solution 1

Listing 4.4: Solution 1 In snippet 4.4 we introduce a static method apply which is responsible to create the functor and launch the parallel execution. This answers all aforementioned issues. In the case of a template functor, one also uses explicit template instantiation for a set of template parameters commonly used.

Listing 4.5: Solution 1 However this introduces a minor issue. As we can see in snippet 4.5, it is a valid code but compiler is not able to inline function call as the definition is not visible in this translation unit. Hence inter-procedural optimization (IPO) can be used to trigger inlining. A related issue happens if the declaration of the function call operator is marked to be forced to be inlined. Hence snippet 4.5 will not even compile.

4.A.2 Solution 2

A more flexible approach would involve a third file.