domain: phys spi

text:
Keywords: NUMERICAL METHODS, ENERGETICS, MULTIPHYSICS, STEADY

Même si trois années passent vite, une thèse est une véritable période de vie, durant laquelle beaucoup de personnes se retrouvent investies. Je tiens à remercier ici ces personnes qui m'auront accompagné, aidé et soutenu durant cette période.

Je commence par remercier les membres du jury qui ont accepté leur rôle. Je remercie en particulier les rapporteurs. Comme j'ai eu l'occasion de leur dire, leurs rapports sont d'un grand intérêt car ils permettent de prendre du recul sur son propre travail, en donnant un point de vue critique mais surtout extérieur sur la thèse. Enfin, je remercie Xavier Vasseur qui n'était pas à priori lié à cette thèse mais qui y a contribué, en répondant à mes questions, me fournissant des références et enfin en acceptant le rôle de juré.

Ma thèse s'étant déroulée initialement à Toulouse puis à Châtillon, je n'ai jamais été physiquement avec mes deux encadrants en même temps, à part de rares occasions. Cependant, je ne me suis jamais senti seul dans mon travail. Que ce soit au téléphone, par chat, par mail ou lors du point hebdomadaire, vous étiez toujours disponible. Pour ceci, Guillaume et Lionel, je vous suis extrêmement reconnaissant et je ne peux que recommander un tel encadrement de thèse. Je pense aussi que vos qualités d'encadrants sont liées à vos personnalités : en dehors du contexte professionnel, nous nous sommes toujours bien entendus, ce qui facilite d'autant plus nos échanges. Merci à tous les deux, à Guillaume pour ton point de vue optimiste lorsque je n'en avais pas, à Lionel de m'avoir laissé te harceler par chat et au téléphone. Je suis convaincu que vous avez contribué au bon déroulement de cette aventure, et pour cela merci encore.

J'ai donc commencé ma thèse dans l'unité HEAT avant de déménager dans l'unité PLM. Ce fut pour moi l'occasion de rencontrer davantage de personnes qui ont ajouté leur pierre à l'édifice. Certes, cela veut en fait dire que j'avais d'avantage de personnes à aller embêter avec mes questions, mais j'avais toujours de l'aide en retour. Je remercie donc les membres de ces deux unités, et même plus largement les membres du DMPE, qui m'ont aidé durant ces trois ans.

Parmi eux, je remercie en particulier Jean-Michel Lamet. C'est parce qu'il s'est intéressé à ma thèse que nous avons pu trouver un avantage immédiat à l'utilisation de mes travaux. Merci pour le temps que tu m'as consacré, à re-travailler ton modèle, échanger avec moi, et même relire une partie de ce mémoire ! Plus largement, je tiens à remercier l'ONERA et son personnel. Je m'adresse ici aux personnes qui ont n'ont pas contribué à ma thèse d'un point de vue scientifique, mais qui ont tout de même apporté leur aide.

ii iii En-dehors du personnel ONERA, je remercie les différentes personnes qui ont fait partie du support CEDRE. Je pense notamment à Julien R. qui m'impressionna par ses compétences et son investissement. Je pense également à Alexandre et Christ, avec qui j'ai pu râler sur les problèmes et les bugs que je rencontrais. À ces personnes s'ajoutent Julien V., Yacine, ou encore mes collègues doctorants. Merci à vous tous pour ces moments d'échanges parfois scientifique, parfois moins.

Il ne faut pas oublier le plus grand contributeur de cette thèse, sans qui ce mémoire n'existerait pas aujourd'hui. Merci Zhao pour ton ordi ! Je remercie enfin ma famille pour son soutien, et en particulier mes parents. Vous dites être fiers de moi, mais je suis comme je suis de par votre éducation et votre amour, merci. Merci Manon, de partager ma vie, et d'en faire un bonheur.

Pour mamie Raymonde, qui disait que plus tard, je serai ingénieur.

Introduction

General context

This work is rooted in the field of numerical simulation of fluid dynamics, applied to aeronautics, defense, and space industries. This field gathers many industrial actors (the French Direction Générale de l'Armement, ArianeGroup, Safran, Airbus, etc) as well as academics (ON-ERA, CERFACS, DLR, VKI, universities, etc). It is interested in how to simulate a fluid flow with a computer, trying to represent faithfully the physical reality. The various actors in this field need to be able to access certain physical quantities associated with specific phenomena and operating regimes. These regimes are often not feasible on our scale, due to material or financial limitations. Examples include the study of icing on the wing of an aircraft, which is experimentally feasible but represents an imposing budget for the aircraft manufacturer, or the study of heat transfer in an atmospheric reentry capsule, which is much more difficult to achieve experimentally. To overcome these limitations, numerical simulation is the best option, as it allows such a case study to be modelled by the execution of a computer program, and to obtain a large set of data that will be analysed afterwards to answer the desired questions.

The analysis of physics usually produces a set of equations, often partial differential equations, representing the real system that one wishes to study. Algorithms are then required to determine the fluid flow from these equations, in the working domain, as a function of time. Thus, to obtain the desired quantities, the physical system is integrated using mathematical algorithms to obtain its evolution in time. Often, the expected result is not the complete temporal evolution of the system, but only its equilibrium state. This is called a steady numerical simulation or steady computation. Steady computations are opposed to unsteady computations that aim to accurately describe the system's temporal evolution.

For all computational fluid dynamics users, the most crucial issue is to achieve a satisfying compromise between computational cost and the accuracy of the results. Indeed, a quick and inexpensive computation tends to be not very faithful to physics, while an accurate computation tends to use more computational resources and time. For steady computations, speed corresponds to getting the final steady state at a low time cost, for both computational time and the time it actually took to a user. It results in a compromise between methods that are expensive in terms of computational resources and take a long time but give accurate results that are close to the physical reality, and faster methods that save resources but give lesser quality results. A software developer working on a numerical simulation tool needs to choose methods and algorithms to obtain a compromise that is considered satisfactory. The final interest for a player in computational fluid dynamics is therefore to have a result that is sufficiently precise and inexpensive enough to obtain. An accurate result is needed to answer the questions that required 1 INTRODUCTION the simulation. The search for a result that is inexpensive to obtain is motivated by questions of savings in computational cost. This cost is applied to the user, while he waits for the simulation to end, and to the company through the cost of computer resources, electricity, investments in more efficient machines, etc.

Depending on the problem to be solved, there are more or less suitable algorithms and methods. Performant methods were originally developed to answer the needs of the aerodynamics community. In the case of energetics and multiphysics problems, the methods issued from the more traditional aerodynamics are limited by the coupling between the different physics which have their distinct characteristic times. Thus, the algorithms used for the numerical simulation of classical fluid dynamics are focused on the resolution of the Navier-Stokes equations, and are not necessarily the most adapted to a simulation in the field of energetics. The involvement of several distinct physical phenomena imposes constraints on the choice and use of algorithms. Consequently, it would be advisable to adapt or replace the algorithms involved in the time integration for multiphysics problems.

Context of this thesis

This thesis was conducted at ONERA, the French Aerospace Lab, in the multiphysics for energetics department. Although numerical simulation has been greatly developed in the aeronautical field, it has not been so well adapted to the multiphysics field, and many industrial codes are content to reuse the same algorithms. This is the example of the software system CEDRE, developed at ONERA by the multiphysics for energetics department [START_REF] Refloch | CEDRE software[END_REF]. This software system constitutes a platform grouping several solvers to integrate several physics: each solver is dedicated to its physical model. There is a solver for the resolution of compressible, multifluid, reactive, and turbulent flows, two solvers for the calculation of dispersed phase (drops, crystals, particles) in the Eulerian and Lagrangian approach respectively, a solver dedicated to the calculation of liquid films, a solver dedicated to radiation, etc. Thus, CEDRE is in fact a global platform adapted to multiphysics problems. Older work has been done to set up a time integration adapted to the problems solved by CEDRE [START_REF] Selva | Méthodes itératives pour l'intégration implicite des équations de l'aérothermochimie sur des maillages non-structurés[END_REF]. Work on time integration has led to the development of implicit integration methods for the integration of steady problems, and to the development of a GMRES method for the solution of linear systems. Users can choose from a panel of integration methods to obtain methods adapted to their problems. However, the weak coupling between solvers hampers the convergence to the steady state and the choice in integration methods is limited, compared to what can be found in the literature. This is at least the opinion of the actors of the CEDRE code, i.e. its developers and users, who would like more robust methods to use CEDRE on steeper problems and converge faster to save on computational costs.

On the research side, however, many efforts have been made in the direction of multiphysics, but have not yet left the academic framework. This is for example the case of [START_REF] Wong | Sequential-implicit Newton method for multiphysics simulation[END_REF], who were interested in the temporal integration of coupled equations. They designed an integration method adapted to coupled equations, which is an evolution of a fixed point method with the addition of a step of a Newton method. They then compared this method to the standard fixed point method, which is more commonly used for coupled computations. By implementing their method on two simple coupling problems, they finally showed the value of their method compared to the more standard method. If this method is well fitted for their calculations, it is however only implemented for simple problems, less complex than the multiphysics problems that CEDRE wishes to solve. Moreover, the tests carried out are on problems on an academic scale, and not on an industrial scale. The same issue concerns [START_REF] Coulette | Implicit time schemes for compressible fluid models based on relaxation methods[END_REF], in which authors used relaxation techniques to transform the original equation into a simpler one. Furthermore, this seems to be specific to a fluid model, and would not scale well in CEDRE where many different models coexist. A similar idea uses splitting of the multiphysics parts to use IMplicit-Explicit methods [START_REF] Huang | High-order partitioned spectral deferred correction solvers for multiphysics problems[END_REF]. It seems well fitted for out multiphysics applications, however we assumed that having a better implicit solver was a priority.

At the same time, tools used in aerodynamic simulation could prove interesting for multiphysics problems. The Jacobian-Free Newton-Krylov method, or JFNK, is already well used in the numerical simulation of the Navier-Stokes equations. In [START_REF] Park | On physics-based preconditioning of the Navier-Stokes equations[END_REF], a time integration mechanism is implemented around a JFNK formulation. A physics-based preconditioning is developed to solve the linear problem more precisely. This time integration is then tested on a thermally driven cavity problem. However, this work is only concerned with the Navier-Stokes equations and is not extended to more general energetics problems. Finally, it is only tested on a two-dimensional academic problem. In [START_REF] Content | Coupled/uncoupled solutions of RANS equations using a Jacobian-free Newton-Krylov method[END_REF], it was also used on academic problems for the turbulent Navier-Stokes equations. The Jacobian-Free Newton-Krylov method seems to be well fitted to problems that couple multiple physics. It was used in [START_REF] Turpault | Modélisation, approximation numérique et applications du transfert radiatif en déséquilibre spectral couplé avec l'hydrodynamique[END_REF] on applications coupling radiative transfer with hydrodynamics. It could prove interesting as a time integration method for CEDRE multiphysics problems.

This thesis was co-funded by ONERA and AID, the French agency for innovation and defence.

General overview and motivations

It can therefore be seen that numerical methods are already available to solve computational fluid dynamics problems. In particular, there are already solvers capable of solving steady multiphysics problems of industrial scales. However, such solver use methods inherited from the aerodynamic context, that are not always adapted to handle the difficulties arising from the multiphysics properties of our problems. It is at least the case for our solver CEDRE. On the other hand, other methods developed in an academic context have shown their interest in simple multiphysics problems. It would now seem interesting to adapt these new methods for the resolution of multiphysics problems in an industrial code. This is the reason for this study. It consists in improving the convergence, speed and robustness of the time integration of the CEDRE platform on steady multiphysics problems by adding numerical methods not used in industrial numerical simulation. The aim of this work is to improve the current performances of our solver by using such more recent methods.

INTRODUCTION

Outline

This document is divided into two parts. The first part, consisting of Chapters 1, 2, and 3, describes the work we did with the solver CEDRE, around the Jacobian-Free Newton-Krylov method. Some of it was presented at the 2022 Aviation Forum [START_REF] Seize | Towards multi-physics simulations using Jacobianfree Newton-Krylov method[END_REF]. The second part, consisting of Chapters 4 and 5, focuses on exponential time integration methods, initially investigated in CEDRE and then in another solver: JAGUAR. The exponential time integrator shares many ingredients with the JFNK method but is dedicated to unsteady flows. The novelty of the present part is the association of exponential integration and an advanced discretization technique called the Spectral Difference method. We plan to submit an article just after the defense on this research line.

The objective of Chapter 1 is to identify, from the literature, methods that we think can improve CEDRE performances, meaning speed, accuracy and robustness. We will introduce time integration methods and the mathematical tools we need to analyse them. We will go into more detail about implicit time integration methods by describing each step involved in its implementation. Those details will lead us from the differential equation to the algebraic solver used to solve linear systems. This chapter also describes what already existed in CEDRE before this thesis. It shows the strengths and the weaknesses of the traditional time integration methodology. The idea is then to work on those weaknesses by using previously identified methods. Indeed, we will look at the literature and choose new methods to implement that go well with what we want to keep, in place of what we want to discard. At this stage, we will have identified methods that we hope will improve our solver performances.

Chapter 2 aims to present the implementation of the new methods in CEDRE. To do so, we will briefly describe the structure of the solver and its different parts. We will then go into more detail and explain some implementation choices. Finally, we will justify those choices and show they are indeed valid to accomplish what set out to do. After this, we have a method that we need to test to see if it improves our solver performances.

The goal of Chapter 3 is to evaluate the robustness, accuracy and speed of the previously implemented method, and compare those performances to those of previously existing methods. The idea is to select representative test cases, that correspond to typical CEDRE applications. We will present at first simple aerodynamic test cases, and we will increase the complexity of the physical model as we go. For each of those test cases, we will compare the newly implemented method to the older one. This will justify the choices made in Chapter 1 and their implementation in Chapter 2. Finally, we will use the newly implemented method on an also newly implemented fluid model and show that this lead to some benefits for the solver users.

In Chapter 4, we will step out of the context of steady problems and extend our work to unsteady computations that use larger time steps. The objective of this chapter is to introduce a new category of time integration methods and give some of their characteristics. Then, we will show on a simple analytic test case why those new methods may be of interest to our computational fluid dynamics applications. We will also see that we need to use another solver to study properly those methods.

Chapter 5 aims to show the potential of this new class of time integration methods in our field. To do so, we use another solver, as was recommended in the previous chapter. After introducing the novelties of this other solver, we will define a set of methods we want to compare to one another, made of previously existing and newly added methods. Then, we describe several test cases to characterise and compare those methods.

Finally, we will conclude this thesis and give perspectives that can continue this work further.

Part I

Solving efficiently multiphysics steady problems

Chapter 1

Analysis of existing methods

Résumé du chapitre : Analyse des méthodes existantes Le but de ce chapitre est d'identifier les méthodes qui serviront pour la suite. Le contexte mathématique et ses notations sont introduits en premier. L'objectif est de donner les étapes permettant de passer d'une équation aux dérivées partielles, issue du modèle physique, à une équation différentielle ordinaire, purement numérique. Pour fermer le système d'équations, les méthodes de discrétisation spatiale sont brièvement introduites et en particulier les deux principales du solveur CEDRE. Grâce à ces méthodes de discrétisation appliquées à un maillage, l'équation différentielle ordinaire finale est donc obtenue et il faut alors l'intégrer avec des méthodes d'intégration temporelle.

La stabilité et la notion d'ordre sont les deux principales caractéristiques nécessaires à de telles méthodes. Puis, les méthodes d'intégration que l'on trouve classiquement dans la littérature sont rappelées, en soulignant la dichotomie explicite/implicite.

Après avoir détaillé les raisons qui nous poussent à analyser les méthodes d'intégration implicite, plusieurs éléments interviennent dans la mise en place : le solveur non-linéaire, le solveur linéaire, et en particulier le passage de l'un à l'autre. Ce passage étant jugé insatisfaisant pour notre solveur, nous présentons la méthode JFNK, ou plus généralement une formulation sans matrice, qui s'appuie sur des parties déjà présentes de CEDRE et corrige, a priori, ce défaut.

As we aim to enhance the performances of the solver CEDRE on steady problems, we first have to define some notions. In particular, we need to explain what performance means. In this chapter, we will introduce the general problem, and the tools used to solve it. We will confirm some choices already made in the solver, and identify the features we would like to modify or replace.

CHAPTER 1. ANALYSIS OF EXISTING METHODS

Problem setup

In order to set the mathematical framework for this study, let us start with a partial differential equation arising from the physical model, in the form of

∂ξ ∂t + F (ξ) = 0 (1.1)
where the function F uses some space derivatives of the state variable ξ. This equation then describes the temporal evolution of the state variables ξ.

A particular class of such partial differential equations are conservative equations. They correspond to the case where the function F can be written as a divergence term. Finally, with a source term S, those equations look like:

∂ξ ∂t + ∇ • f (ξ) = S . (1.2)
One might notice that equation (1.2) is indeed a particularisation of equation (1.1), with F (ξ) = ∇ • f (ξ) -S. Those conservative equations are the ones we will focus on in this study, as they describe the physical systems we are interested in.

In this work, attention is paid to computational fluid dynamics and in particular to solving the Navier-Stokes equation and its variants: the reactive Navier-Stokes equation, the Reynoldaveraged Navier-Stokes equation, etc. A simple form of this equation can be:

     ∂ t (ρ) + ∇ • (ρu) = 0 ∂ t (ρu) + ∇ • ρu ⊗ u + pId = ∇ • τ ∂ t (ρE) + ∇ • ((ρE + p) u) = ∇ • τ • u -q (1.3)
with the closing relation ρE

= p γ-1 + ρ u•u
2 . This relation is quite simple, as it is the one that corresponds to ideal gas with constant heat capacity. In the solver CEDRE, this relation is usually more complex as we are working with multifluid flows and with variable heat capacities. The deviatoric stress tensor τ accounts for the fluid's viscosity, and its computation depends on the model used. The heat flux term q enables heat exchange in the fluid, through the heat diffusivity, and it is generally related to the temperature gradient. Without those last two terms, one recovers the Euler equations. To this simple form can be added source terms from the reactive model, source terms from the turbulence model, divergence terms from a molecular diffusion model, etc. Yet it is clear that with a bit of rewriting, it is possible to get back to the starting form (1.1) and even the conservative form with source terms (1.2). The quantity ξ is no longer a scalar but a vector with the density ρ, each component of the momentum ρu and the energy ρE as its components. Apart from this small change, the idea is the same.

When solving equations like (1.1) numerically, the domain of interest is first bounded since the infinite domain cannot be represented by a discrete set of variables. Let D be this computational domain. For numerical computation, the different quantities, such as the state variable ξ, are represented discretely over the domain D and are stored in the memory of a computer. To do so, the domain D is split into a finite number of cells, or elements, associated with the degrees of freedom of the solution: this is the mesh. Those cells are small disjoints volumes in 3D, faces in the two-dimensional space or segments in 1D, such as their union recovers the original domain.

BRIEF INTRODUCTION TO THE SPATIAL INTEGRATION SCHEMES

Interest quantities, such as the fluid velocity, density, . . . , are then stored at each node, averaged at the center of each cell or sometimes in a more complex fashion depending on the method. They are no longer mathematically represented by a function of the continuous physical domain ξ : D → R but by a finite-sized vector Ξ gathering all the information across the discretised domain. For some simple discretisation methods, this vector consists of the quantity evaluated at the mesh nodes or averaged at the center of the cells. For more complex methods, this vector consists of information used to construct the solution over the domain: polynomial coefficients, spectral decomposition coefficients, etc. Anyway, a discretised domain is used rather than D, the continuous one.

The partial differential equation (1.1) transforms then into an ordinary differential equation:

dΞ dt -G (Ξ) = 0. (1.4)
The difference here is that the function G is a function of a discrete vector whereas F was a function of continuous state variables that are in turn functions of the space variable, and therefore G does not use any spatial derivatives. The negative sign is here so that this new equation is similar to the original partial differential equation (1.1) and the new equation can be written in the form that corresponds to most of the literature:

dΞ dt = G (Ξ) . ( 1.5)
Thanks to the spatial discretisation method, the only derivative remaining is with regard to time.

The rest is then up to the temporal integration method, which is the main topic of this thesis. This integration method will works on equation (1.5) no matter where the function G comes from, but sometimes understanding the origin of this function can help so we will now introduce the spatial discretisation method used in our solver.

Brief introduction to the spatial integration schemes

A spatial discretisation method is what tells how to represent a quantity over a discretised domain, and how to compute the spatial derivative of this quantity from this representation. Indeed, before solving equation (1.1), one must decide how to transform the continuous model into a discretised one. We also have to look at how the spatial derivatives arising from equation (1.1) translate in the discretised model.

The Finite Volume method

The spatial discretisation method used in the fluid solver CHARME is called the Finite Volume method [START_REF] Eymard | Finite volume methods[END_REF][START_REF] Leterrier | Discrétisation spatiale en maillage non-structuré de type général[END_REF]. This method is particularly well fitted for conservative equations such as equation (1.2). A such equation has the property that the quantity ξ is conserved: without source terms, the variation of the total quantity ξ over the domain D is equal to the flux f (ξ) coming through the boundary ∂D. In the case of the Navier-Stokes equations (1.3), the density, the momentum and the energy are conserved throughout time, apart from fluxes going through domain boundaries. In a close domain where nothing comes in or out, they are indeed conserved. The main interest of the Finite Volume method is that this property stays true through the spatial discretisation step.

CHAPTER 1. ANALYSIS OF EXISTING METHODS

The Finite Volume method consists in integrating the partial differential equation over each cell of the mesh. Writing V i the volume of the ith cell:

Vi ∂ξ ∂t dv + Vi ∇ • f (ξ) dv = Vi S dv. (1.6)
Then the Green-Ostrogradski theorem transforms the flux divergence into a surface integral:

d dt Vi ξdv + ∂Vi f (ξ) • ds = Vi S dv. (1.7)
By writing i = 1

Vi

Vi dv the average in the ith cell, it comes:

dξ i dt + 1 V i ∂Vi f (ξ) • ds = S i . (1.8)
As stated before, the spatial discretisation method does transform the partial differential equation into an ordinary differential equation. It tells to store quantities as their averaged values represented at the center of gravity of each cell as the vector Ξ. An important aspect concerns the transformation of the divergence from equation (1.2) into a flux surface average, and the flux computation is still here an open question. The mesh cells are supposed to be polyhedrons and they have a finite number of planar faces. The integral over the boundary of the cell can be decomposed by the faces, to get the approximation:

∂Vi f (ξ) • ds ≈ j neighbor of i f ij • s ij (1.9)
where s ij is the product of the face area and the unit normal vector directed outwards volume i and f ij • s ij is an approximation of the flux going through the face between cells i and j. This approximation is a key element of the Finite Volume method and will be discussed later. The function G from equation (1.5) can now be computed: for each face of the mesh, one computes f ij • s ij , adds this value to the ith component and removes it from the jth component of the new vector. Then, after adding the source terms, this yields the vector containing the result of G (Ξ). As can be seen, every contribution of the flux added in a cell is removed from another, and therefore this spatial discretisation method preserves the conservation of the underlying equation.

The Riemann problem

The last remaining problem with this presentation of the Finite Volume method is how to compute the flux going through cell interfaces. On the interfaces between two cells, the left and right quantities ξ L and ξ R are known and are used to compute the corresponding flux. It is possible here to use a reconstruction method to get a better approximation of the quantities at the left and right sides of the interface, and therefore the left and right quantities at the face ξ * L and ξ * R are used, rather than the quantities at the cell center. The idea is now to compute the flux going through the face as a function of ξ * L , ξ * R and the surface vector s. From the interface point of view, there are two possible different states, one from each side: this is what is called a Riemann problem. A Riemann problem is an initial value problem applied to a conservation equation, where the initial solution is piecewise constant with a single possible discontinuity. By working with the equation and deriving the jump condition, it is possible to

BRIEF INTRODUCTION TO THE SPATIAL INTEGRATION SCHEMES

compute the quantity from a possibly discontinuous state at the interface. Then it is possible to evaluate the flux associated with this state going through the surface. This approach can be called the exact Riemann solver as it uses the exact solution of the Riemann problem. But the drawback of this approach usually is the computational cost required to find this exact solution and people generally revert to approximate Riemann solvers, compromising between speed and accuracy. Several approximate Riemann solvers are available to the user in our solver, such as the well-known Roe, HLLC or AUSM+ schemes [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF].

Gradient reconstruction methods

The standard Finite Volume method represents quantities with the averaged value in each cell. This corresponds to a first-order discretisation method. Simply put, it means that it can represent quantities exactly as 0-order polynomials locally to each cell. There are ways to achieve higher-order representation such as with the MUSCL approach [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF]. It consists in handling the surface flux evaluation as explained in the Riemann problem part on one hand and deciding what left and right quantities to feed to this flux computation on the other hand. In our solver, there are two ways to construct high-order states to give to the flux computation method. They are described in the following parts. For both of them, the idea is to use neighbouring data to enhance the order of the local representation.

The k-exact method

The first method used to reconstruct high-order quantities is called the k-exact method with successive corrections. The idea is to construct iteratively an order k representation of the quantity using the neighbouring order k-1 representation [START_REF] Haider | Stability analysis of the cell centered finite-volume MUSCL method on unstructured grids[END_REF]. This method is often employed to achieve a second-order reconstruction, as choose most users, but it can also achieve higher-order reconstructions [START_REF] Haider | Efficient implementation of high order reconstruction in finite volume methods[END_REF][START_REF] Haider | A high-order interpolation for the finite volume method: the coupled least squares reconstruction[END_REF][START_REF] Pont | Multiple-correction hybrid k-exact schemes for high-order compressible RANS-LES simulations on general unstructured grids[END_REF]. At each step, while increasing the order of the representation, it is important not to create a local maximum or minimum. This might happen close to discontinuities in the solution, or near rapidly varying spots. It is indeed common, when interpolating, to create local overshoot or undershoot. One might think here about Gibbs or Runge's phenomena, and despite the problem here being a different one, the idea is the same. Creating local extrema in the solution can be troublesome, and so the k-exact method includes limiters to limit reconstructed data.

The Multislope method

The second method used to reconstruct high-order quantities is called the Multislope method. This method is a direct adaptation of the structured directional method to unstructured grids. It consists in defining local directions for the interpolation of the unknown on the face. For a given cell and one of its faces, one can draw the line from their centers. Information is interpolated on this line using neighbouring cells to get the local slope. The mechanism to define the slopes and the geometric treatment is quite complex and it won't be discussed here. The details are provided in [START_REF] Le Touze | Multislope MUSCL method for general unstructured meshes[END_REF]. Finally, the Multislope method gives a second-order reconstruction. Once again, this reconstruction might create local maxima or minima, and therefore it uses slope limiters [START_REF] Venkatakrishnan | On the accuracy of limiters and convergence to steady state solutions[END_REF][START_REF] Berger | Analysis of slope limiters on irregular grids[END_REF] to prevent it.

We briefly explained the spatial discretisation method used in our solver, as it might help the analysis of the time integration part. The Finite Volume method averages the partial differential equation over each cell of the mesh, which transforms the flux divergence into a flux surface average. A reconstruction method, the k-exact method or the Multislope method, is then used CHAPTER 1. ANALYSIS OF EXISTING METHODS to get a higher-order representation of the solution so that the surface flux can be computed at each face. Slope limiters are used to prevent the formation of local extrema, which can be harmful to the computation. The choices made for the spatial discretisation methods are motivated by the fact that CEDRE and in particular the fluid solver CHARME aim to work with general unstructured meshes. It restricts the choices of algorithms, as handling general unstructured meshes can prove difficult for simple spatial discretisation methods. It means that the methods must be more robust, which often means less precise.

Introduction to time integration methods

With the help of a spatial discretisation method, the equation to solve is now an ordinary differential equation. The main objective of this thesis is focused on the resolution of steady problems. The steady solution of equation (1.5) is given by G (Ξ) = 0. To get the solution, one might then try to find a root of the function G. Unfortunately, with our typical applications, this function G has got bad mathematical properties, such as its stiffness, arising from the nonlinearities of the underlying equations. Therefore, algorithms that try to find a root of G struggle and usually fail. Another approach is to take an initial value Ξ 0 , and to solve the equation (1.5) for this initial value. After a long enough time, one hopes that Ξ will reach the desired steady solution.

The idea is now to solve the temporal evolution of Ξ to get the solution after a long time when it approaches the steady solution. The equation is solved numerically, which means the next solution is computed iteratively after a given time step, knowing the current one. It is also possible to modify the equation, as the interest is in the final state, not in the transient one. It is possible, for example, to use local time-stepping, which consists in having each cell of the mesh move forward in time with its own time step. The resulting transient states do not make sense from a physical point of view, as the equation solved is not the initial one, but it converges to the same steady solution. Therefore, it is alright to change the equation as long as it gives the same steady solution. Finally, this way of finding a converged steady solution is what is called a Pseudo-Transient Continuation method [START_REF] Kelley | Convergence analysis of Pseudo-Transient Continuation[END_REF].

After deciding on an initial value, the equation to solve is:

   dΞ dt = G (Ξ) Ξ (t 0 ) = Ξ 0 .
(1.10)

A time integration method is going to produce a succession of solutions: starting from Ξ 0 , it produces Ξ 1 at t 1 , then Ξ 2 at t 2 , and Ξ n at t n , etc. Let ∆t n be the time step t n+1 -t n . The following is related to a single step of the time integration method, so the subscript on the time step can be dropped as it is not meaningful. For a steady problem, the evolution of the solution has no interest and it seems reasonable to want to "go fast" to the steady state, meaning to use as large a time step as possible. Unfortunately, not every time integration method allows large time steps due to the well-known Courant-Friedrichs-Lewy condition [START_REF] Courant | On the partial difference equations of mathematical physics[END_REF]. Some tools must be defined in order to help decide on the method to use.

INTRODUCTION TO TIME INTEGRATION METHODS

15

Analysis of time integration methods

Consistency and order

A time integration method must respect some properties to be "well-behaved". For instance, it has to be consistent. To define the consistency, let us look at equation (1.10). After one step, a numerical method gives a value Ξ 1 , believed to be near the exact value Ξ (t 0 + ∆t). A numerical time integration method is said to be consistent if:

lim ∆t→0 Ξ 1 -Ξ (t 0 + ∆t) ∆t = 0. (1.11)
Also, the method is of order p if the local error is in ∆t p+1 [START_REF] Iserles | A first course in the numerical analysis of differential equations[END_REF]:

Ξ 1 -Ξ (t 0 + ∆t) = O ∆t p+1 . (1.12)
This means a p-order method can recover exactly a solution that is a polynomial function of time with an order less or equal to p.

Note:

The order of a time integration method reflects its "local" behaviour, meaning on a single given time step, provided it is small enough. In the field of spatial discretisation of partial differential equations, the order p of a method is such as:

Ξ -Ξ exact = O (h p ) (1.13)
where h is the spatial discretisation parameter. There is a difference between the two definitions: the error order of magnitude is p + 1 for the temporal method and p for the spatial one. This is due to the fact that the error in the spatial case is global: it sums the error over the whole domain. It would correspond to counting the error on each step for the temporal integration. To convince oneself, one could say that when solving the differential equation on an interval [0, T ] with a fixed T , the global error of a p-order method would behave as O (∆t p ) as it amounts to summing T /∆t local errors of O ∆t p+1 . The coherency with the definition of the order for spatial discretisation methods is now clear. If this trick can help understand the difference between the two definitions, this is indeed just a mental trick and not a rigorous mathematical proof. To get this proof, more hypotheses on the method are required [START_REF] Iserles | A first course in the numerical analysis of differential equations[END_REF].

Stability

A meaningful criterion in the choice of a time integration method is stability. Depending on the application, different levels of stability are expected in order to avoid a numerically induced divergence of the computation.

The stability of a time integration method is usually analysed on the ordinary differential equation with a linear right-hand side [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF]. The reason is that if Ξ is solution of equation ( 1 This new equation used to analyse the stability of time integration methods is the one called the Dahlquist test equation. This equation is studied in C, so that eigenvalues and eigenvectors of the matrix J exist.

CHAPTER 1. ANALYSIS OF EXISTING METHODS

Note:

When looking at a method applied to the Dahlquist test equation (1.14), it is assumed that the real parts of the eigenvalues of J are all negative. This choice may seem arbitrary but can be understood with the following example. Let us work in C 2 , with:

J = -1 0 0 10 3 , y 0 = 1 0 . (1.15)
The solution of the equation is then:

y (t) = e -t 0 . (1.16)
As the equation is solved numerically, the floating point representation introduces some roundoff error. The initial condition may then be

y 0 = 1 (1.17)
instead of the exact one y 0 , with a typical = 10 -15 for double precision. Let us suppose that we have an exact time integration method that gives the exact solution at each time step t n = n∆t.

The solution computed by this method will be:

y n = e -n∆t e 10 3 n∆t (1.18)
that gives an error of e 10 3 n∆t . For the numerical values suggested here, this amount to an error as large as 10 6 for n = 5 and 10 28 for n = 10. The explosion of the error comes from the fact that the positive eigenvalue of J amplifies the roundoff error. This phenomenon has nothing to do with the time integration scheme, but with the equation. Finally, this is why equation (1.14) is studied assuming the eigenvalues of J are negative.

Single-step methods

To compute the solution at the next time step, some methods need only to know the current solution. Such methods are called single-step methods. A single-step method applied to the Dahlquist test equation (1.14) leads to the following relationship between the current solution and the next one:

y n+1 = g (∆tJ) y n . (1.19)
For most time integration methods, g is an analytic function. The initial value y 0 can be decomposed on a basis of eigenvectors of J, v 1 , . . . , v N , associated with the eigenvalues α 1 , . . . , α N :

y 0 = N i=1 λ i v i .
Because g is an analytic function, it comes immediately:

y n = N i=1 λ i g (∆tα i ) n v i . (1.20)
It is now straightforward to deduce a stability condition for the single-step method: if for any i there is |g (∆tα i )| < 1, then y n converges to 0. This is how the stability region of a time integration method is defined:

{ z ∈ C | |g (z)| < 1 } . (1.21)
When each eigenvalue of ∆tJ falls in the stability region, then the method is stable.

INTRODUCTION TO TIME INTEGRATION METHODS

17

If an eigenvalue is not in the stability region, the associated eigendirection will be amplified and a numerical instability will lead to the divergence of the computation. We can see that the argument of the function g is not J but ∆tJ. This means that stability can be achieved by choosing wisely the time step: with a small enough ∆t one can ensure that each eigenvalue falls into the stability region. Unfortunately, this often forces the user to set a relatively small time step, which is an issue for steady computations.

Multistep methods Some methods do not fall into the previous framework. The multistep methods, in particular, cannot be written under the form of equation (1.19). These methods use not only y n to find y n+1 , but the k previous steps. Applied to equation (1.14), they can be written in the form:

y n+1 = k i=1 g k-i (∆tJ) y n+1-i . (1.22)
When looking for y i under the form y i ∝ µ i , one has:

µ k = k i=1 g k-i (∆tJ) µ k-i (1.23)
which leads to identifying the polynomial

g ∆tJ (µ) = µ k - k-1 i=0 g i (∆tJ) µ i .
If each root of this polynomial is of modulus less than 1, the solution converges to 0. This is how the stability region for multistep methods is defined [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF]:

z ∈ C | all roots of X k - k-1 i=0 g i (z) X i
are of modulus less or equal to 1, strictly less to 1 for roots with multiplicity .

(1.24)

Once again, if all eigenvalue of ∆tJ falls in the stability region, then the method is stable.

The key property resulting in the stability analysis of time integration methods is the Astability [START_REF] Dahlquist | A special stability problem for linear multistep methods[END_REF]. A time integration method is A-stable if its stability region contains the left half complex plane. Simply put, a method is A-stable if it converges to 0 when it should, and does not diverge due to numerical errors. The A-stability is interesting, as an A-stable method is also said to be unconditionally stable, whereas a method that is not is conditionally stable. This other characterisation comes from the fact that a non-A-stable method needs to respect some additional criteria to be stable, on the time step it uses for example, whereas an A-stable method is stable no matter the time step. As we said before, we would like to use large time steps to quickly find the steady state of our applications, and that is why we look for A-stability in our methods.

Having defined the concepts of stability and order of time integration methods, we are now equipped to analyse them. We see from the literature that time integration methods are usually split into two groups: explicit and implicit methods. We can now look at methods from this literature, starting with the simpler explicit methods.

CHAPTER 1. ANALYSIS OF EXISTING METHODS

Explicit methods

Explicit methods are called this way because, at each step, the computation of the next solution is straightforward: they give it explicitly as a function of currently available data. They are largely used in unsteady computational fluid dynamics simulations. Their strength comes from the fact that they are usually simple and therefore easy to implement in a solver, and computationally inexpensive compared to non-explicit methods.

Explicit Euler method

The explicit Euler method is the most simple time integration method. It consists in integrating equation (1.10) between t n and t n+1 assuming the function G stays constant, equal to G (Ξ n ). Equivalently, it consists in replacing the time derivative dΞ dt by a finite difference Ξn+1-Ξn ∆t and evaluating G in Ξ n . Then, the method gives:

Ξ n+1 = Ξ n + ∆t G (Ξ n ) . (1.25)
After verifying that this is a first-order method, a quick stability analysis gives a stability region equal to the open unity disk centred in -1. Practically, this stability region is often deemed unsatisfactory as it forces the use of small time steps. Yet this method is a classic that must be introduced before talking about more complex methods.

Runge-Kutta methods

Instead of making one step forward in time as the explicit Euler method, a Runge-Kutta method will make a set of intermediate steps, and find the final solution as a combination of those intermediate steps. It starts from the value of Ξ n known in t n and given intermediates steps t n,i = t n + c i ∆t for 1 ≤ i ≤ k, with a fixed k. An exact integration between t n and t n,i of equation (1.10) leads to:

Ξ (t n,i ) = Ξ n + ∆t tn,i tn G (Ξ (t)) dt.
(1.26)

The integral on the right-hand side is then approximated with a quadrature using the previously computed intermediate steps:

tn,i tn G (Ξ (t)) dt ≈ i-1 j=1 a ij G (Ξ (t n,j )) . (1.27)
Once each intermediate step is known, equation (1.10) is integrated between t n and t n+1 and the resulting integral is approximated using the intermediate steps.

To sum up, a Runge-Kutta method iterates in the following way:

             Ξ n+1 = Ξ n + ∆t k i=1 b i G (Ξ n,i ) with Ξ n,i = Ξ n + ∆t i-1 j=1 a ij G (Ξ n,j ) .
(1.28) A Runge-Kutta method is characterised by its size k and by the quadrature coefficients: a ij,1≤j<i≤k , b i,1≤i≤k and c i,1≤i≤k . There are as many Runge-Kutta methods as there are choices in the quadrature coefficients, but not all choices give good methods. There are criteria that the coefficients must follow to ensure the consistency of the method, and then criteria with more complexity as the order increases. The quadrature coefficients are often arranged in the Butcher tableau:

c A b T = 0 c 2 a 21 c 3 a 31 a 32 . . . . . . . . . c k a k1 a k2 . . . a k,k-1 b 1 b 2 . . . b k-1 b k . (1.29)
The Butcher tableau of some well-known Runge-Kutta methods are shown in table 1.1. The RK1 method is in fact equivalent to the explicit Euler method. The RK2 method is one of the 2 steps second-order Runge-Kutta method. This one is also called the Midpoint method. The RK4 method is a fourth-order method with 4 steps. This is the most famous Runge-Kutta method, vastly used for explicit time integration of ordinary differential equations.

We note that the coefficients c i do not appear in the definition from equation (1.28). It is because the differential equation (1.5) is autonomous: the function G does not depend explicitly on the time. If it did, the equation would then be non-autonomous, and G would be a function of two variables: the state vector Ξ and the time t. Equation (1.28) becomes then:

             Ξ n+1 = Ξ n + ∆t k i=1 b i G (Ξ n,i , t n + c i ∆t) with Ξ n,i = Ξ n + ∆t i-1 j=1 a ij G (Ξ n,j , t n + c j ∆t) . (1.30)
It can be shown that the order p of the method is less or equal to the number of steps k. Up to 4 steps, it is possible to choose the quadrature coefficients to have p = k. Above that, getting a lower bound on the order depending on the number of steps is still an open problem as of today. For a Runge-Kutta method of order p, the corresponding function used for the stability analysis is [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF]:

g (z) = 1 + z + z 2 2 + • • • + z p p! + O z p+1 .
(1.31) When p = k, the last term O z p+1 is in fact null. Figure 1.1 shows the stability region of the Runge-Kutta methods of orders up to 4.

We can see that Runge-Kutta methods are not A-stable, and they do not have a large stability region. Increasing the order of the method does increase the stability region, but not enough for our applications. Practically, this imposes the use of small time steps, which is agreeable for unsteady computations but not for our steady ones. If one iteration of the method is inexpensive, the total number of iterations needed to reach the steady solution will make the overall computation too costly.

Adams-Bashforth methods

We could try to use other explicit methods, such as multistep Adams-Bashforth methods. The kth order Adams-Bashforth method uses the last k computed steps to find the next one. One can indeed check that the index k designating the method also corresponds to its order [START_REF] Hairer | Solving ordinary differential equations I. Nonstiff problems[END_REF].

The idea is to apply a Lagrange interpolation of the function G from equation (1.10) in the k last computed points, and then replace G with the interpolation polynomial when integrating from t n to t n+1 . Contrary to the Runge-Kutta methods, as this method reuses previous information, a single G evaluation is required at each step. Because of that, the cost of one iteration of the Adams-Bashforth method is quite inexpensive. The stability analysis for such methods is a bit more complex [START_REF] Hairer | Solving ordinary differential equations I. Nonstiff problems[END_REF][START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF], and so the result obtained numerically is shown in figure 1.2 without giving the details. The conclusion is even worse than with Runge-Kutta methods, as the stability region decreases as the order increases. The Adams-Bashforth methods can reach a high order of accuracy while staying computationally inexpensive, but they drastically lack stability, and that is why they are often not used in computational fluid dynamics computations. More generally, an explicit multistep method cannot be A-stable [START_REF] Dahlquist | A special stability problem for linear multistep methods[END_REF].

Implicit methods

We explained in the previous section why explicit time integration methods are not suited for our applications. It is then natural to look at implicit methods. Contrary to the explicit methods, implicit methods do not give the looked-for solution right away, but as the solution of a specific equation. The name is appropriate: the next state is not given explicitly but implicitly.

Implicit Euler method

The implicit Euler method is the equivalent of the explicit Euler method but on the implicit side. It is quite similar, as it consists in integrating equation (1.10) assuming the function G is constant, but this time equal to G (Ξ n+1 ):

Ξ n+1 = Ξ n + ∆t G (Ξ n+1 ) ⇔ Ξ n+1 -Ξ n -∆t G (Ξ n+1 ) = 0 (1.32)
and the next state Ξ n+1 is given as the solution of a nonlinear problem, the root of a nonlinear function.

After checking that this is a first-order method, the stability analysis gives the corresponding function g (z) = (1 -z)

-1 , which gives a stability region equal to the whole complex plane minus the closed unity disk centered in 1. Therefore this method is A-stable.

CHAPTER 1. ANALYSIS OF EXISTING METHODS

Implicit Runge-Kutta methods

Runge-Kutta methods can also be implicit methods. This happens when the quadratures use points that have not already been computed. In other words, this corresponds to a full A matrix in the Butcher tableau, where it is strictly lower triangular for explicit Runge-Kutta methods. This also means that any step of the Runge-Kutta method may be used in any other step, and therefore one may have to simultaneously solve an implicit system of equations. This may lead to awfully expensive methods, and therefore users tend to restrict themselves to some particular methods.

Despite their cost, implicit Runge-Kutta methods can be appealing. The main reason is that they can easily achieve a high order of accuracy. For example, the methods based on Gauss-Legendre quadratures achieve an order 2k with k steps, and they are all A-stable [START_REF] Iserles | A first course in the numerical analysis of differential equations[END_REF]. Theoretically, this means that an arbitrarily high order can be achieved while keeping the stability quality with these methods. However, users usually stop at the 3-stage 6th order method, as the computational cost tends to be too much for higher order methods.

When applying the stability analysis to a Runge-Kutta method, one can get the function g using the matrix A and the array b:

g (z) = 1 + zb T (Id -zA) -1 (1, . . . , 1) T .
(1.33)

We will not try to show the corresponding stability regions as there are too many methods possible. Instead, we will review some of the most frequent from the literature.

Diagonally Implicit Runge-Kutta methods

The Diagonally Implicit Runge-Kutta methods [START_REF] Alexander | Diagonally implicit Runge-Kutta methods for stiff O.D.E.'s[END_REF], or DIRK methods, are Runge-Kutta methods with a lower triangular matrix A. Then, each step is given as an implicit problem using the already known steps and the next one. The difference is that instead of solving a full implicit system, one just needs to solve each implicit step successively. This help drastically reduces the cost of the method. Furthermore, if all quadrature coefficient a ii are all equals, as each step requires the inversion of the matrix .34) this can help solve the linear problem. This variant is called Singly Diagonally Implicit Runge-Kutta methods (SDIRK) [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF].

Id -a ii ∆t d G dΞ (Ξ n ) , ( 1
Rosenbrock methods Rosenbrock methods are also called linearly implicit Runge-Kutta methods [START_REF] Rosenbrock | Some general implicit processes for the numerical solution of differential equations[END_REF]. They start from the recurrence relation of a DIRK method:

Ξ n,i = Ξ n + ∆t i-1 j=1 a ij G (Ξ n,j ) + ∆ta ii G (Ξ n,i ) (1.35)
and then take its image by G. Instead of solving this nonlinear equation with a nonlinear solver as the DIRK method would do, they linearize the result:

G n,i = G   Ξ n + ∆t i-1 j=1 a ij G n,j   + ∆ta ii d G dΞ   Ξ n + ∆t i-1 j=1 a ij G n,j   G n,i (1.36)
by noting G n,j = G (Ξ n,i ). The recurrence equation for the integration scheme is no longer on the intermediate states but on the intermediate slopes. Usually, the Jacobian matrix is kept constant throughout the stages, equal to the Jacobian matrix evaluated in Ξ n . Additional parameters γ ij are often to add more flexibility to the method [START_REF] Wanner | On the integration of stiff differential equations[END_REF]. It gives the Rosenbrock-Wanner method:

               Ξ n+1 = Ξ n + ∆t k i=1 b i G n,i with G n,i = G   Ξ n + ∆t i-1 j=1 a ij G n,j   + ∆t d G dΞ (Ξ n ) i j=1 γ ij G n,j .
(1.37)

This corresponds to classical Rosenbrock methods when all γ ij are null, except for γ ii that correspond to the a ii . However, the distinction is not made very often, and Rosenbrock-Wanner methods are simply called Rosenbrock methods. The method is defined by the coefficients a ij,1≤j<i≤k , b i,1≤i≤k and γ ij,1≤j≤i≤k . If the problem is non-autonomous, the method is then:

                                       Ξ n+1 = Ξ n + ∆t k i=1 b i G n,i with G n,i = G   Ξ n + ∆t i-1 j=1 a ij G n,j , t n + c i ∆t   + ∆t ∂ G ∂Ξ (Ξ n , t n ) i j=1 γ ij G n,j + ∆t ∂ G ∂t (Ξ n , t n ) i j=1 γ ij (1.38) with c i = i-1
j=1 a ij as for other Runge-Kutta methods. If one takes all the γ ii equals to one another, then the matrix to invert at each step of the Rosenbrock method is the same. This can help reduce the cost of the method by doing a single matrix inversion, or by recycling information between Krylov subspace methods [START_REF] Gaul | Recycling Krylov subspace methods for sequences of linear systems -analysis and applications[END_REF]. Overall, even if they can achieve high-order and A-stability, Rosenbrock methods are still computationally expensive.

Backward differentiation formula

As explicit Runge-Kutta methods are extended to implicit methods, Adams-Bashforth methods also have their implicit counterparts. The idea is to still compute the Lagrange interpolation of the function G but with one additional point: the point at the next time step. This new method is called the Adams-Moulton method. However, the stability region of Adams-Moulton is quite narrow, as they were originally not made for stiff equations [START_REF] Iserles | A first course in the numerical analysis of differential equations[END_REF]. This is why the Backward Differentiation Formula methods, or BDF methods, were introduced. Contrary to Adams-Bashforth and Adams-Moulton methods, they use the Lagrange interpolation of the solution Ξ instead of the interpolation of the function G. Then, the time derivative of the solution can be replaced by the time derivative of the interpolating polynomial, and the resulting equality is evaluated at the next time step t n+1 . This gives an implicit equation that is solved to get Ξ n+1 . The name of those methods comes from the fact that if it uses a constant time step between iterations, this equation can be written using the differentiating operator defined by ∇ 0 i = i and Figure 1.3: Stability regions (in colour) of the BDF methods. Zoom is not constant between plots.

∇ j+1 i = ∇ j i -∇ j i-1 : k i=1 1 i ∇ i Ξ n+1 = ∆t G (Ξ n+1 ) . (1.39)
The order of the method is equal to the index of the method, corresponding to the number of previous states needed to compute the next one. These methods allow for an arbitrarily high order without increasing the cost, as the implicit equation is not harder to solve as the order increases. However, the stability analysis limits the higher order achievable. The stability analysis of the BDF methods can be done numerically. At this stage, one can note that the first-order BDF method is in fact the implicit Euler method. Also, methods of order 7 or higher are unstable, so we can limit our analysis to methods with orders from 1 to 6. The corresponding stability regions are satisfying, as can be seen in figure 1.3. Particularly, the first and second methods are A-stable. We see when looking at the scales in figure 1.3 that the complement of the stability region grows with the index of the method. More generally, there are no A-stable multistep methods with orders higher than 2 [START_REF] Dahlquist | A special stability problem for linear multistep methods[END_REF][START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF]. This introduction to the classic time integration methods helps us to decide what to do for our solver. As we said earlier, we want to solve an ordinary differential equation in order to recover the steady solution, reached after a long time. Then, explicit methods that constrain the time step are not well fitted. Even if they cost more, algorithmically speaking, implicit methods are indeed the best choice for stiff equations. It is often better to do a single expensive iteration over a large time step with an implicit method than make a lot of inexpensive iterations over 1.4. IMPLICIT METHODS FRAMEWORK 25 small time steps with an explicit method. This is finally why we will continue working with the A-stable implicit Euler method. This method is already used in our solver as the base implicit method, for the same reasons we want to use it. Furthermore, it is at the base of all other implicit methods, as they can be seen as small variations of the implicit Euler methods. Choosing it is also smart as updating it into another implicit time integration method will not require too much work.

Implicit methods framework 1.4.1 Methodology of the implicit time integration

As we explained in the previous section, implicit time integration methods give the next state as the solution of a nonlinear equation or a system of nonlinear equations. If we want to use implicit methods, we then need to be able to solve a nonlinear problem. We will continue our discussion using the implicit Euler method, but everything can easily be adapted for any other implicit method. For Diagonally Implicit Runge-Kutta methods, we apply the nonlinear solver successively for each step. For BDF methods, we apply the nonlinear solver to a slightly different equation.

From a nonlinear problem . . . Newton's method can be used to solve a nonlinear problem of the form:

f (Ξ n+1 ) = 0.
(

Equivalently, as it works for a single n at a time, this equation can be expressed in terms of the increment x = Ξ n+1 -Ξ n without writing the subscript:

f (x) = 0. (1.41)
Newton's method starts from an initial guess x 0 . This initial guess is often equal to zero, as it is equivalent to taking Ξ n as an initial guess for Ξ n+1 . The method will then iterate to approximate the solution of equation (1.41).

At each step of Newton's method, the problem is linearised in the current estimation x i and the new equation is evaluated at the next iteration x i+1 :

f (x i+1 ) ≈ f (x i ) + f (x i ) (x i+1 -x i ) = 0 (1.42)
This gives a linear problem in which x i+1 is the solution.

. . . to a linear one

Equation (1.42) can be rewritten into the classic linear problem:

Ax = b (1.43)
where

A = f (x i ), b = -f (x i ) and x = (x i+1 -x i ).
The x notation is reused as it is often used in the literature to name the unknown in linear problems. Such linear problems are found during the iterations of Newton's method, but each is solved for a given iteration number so there is no ambiguity in this notation. A lot of methods were conceived to solve such linear problems, and we will explain later how to handle them. To sum up, we started from the partial differential equation (1.1) arising from the physical model. With a spatial discretisation method, and after choosing an initial value, this transforms into an ordinary differential equation (1.10). For stability reasons, we decided to use implicit time integration methods. Iteratively, such methods are going to produce one or several nonlinear problems in the form of equation (1.41). Newton's method used to solve such problems is going to produce a succession of linear problems in the form of equation (1.43). This complicated sequence of operations describes the time integration procedure used in our solver to find the solution to steady problems. It is schematised in figure 1.4, where the blue circles correspond to the problems being solved, while green circles correspond to the methods used to solve them. In this thesis, we are interested in the non-hatched parts.

Nonlinear solver

As we said, we solve the nonlinear problem with Newton's method. What was done in our solver was in fact a single step of Newton's method, which means a single linearization of the nonlinear problem. Since we could benefit from an actual Newton's method, it was implemented. Contrary to the single linearisation that was done before, several linear problems need to be solved during one single time step, so several Jacobian matrices are needed.

The issue when using Newton's method is that it should be accompanied by a Line Search algorithm to determine the step size α so that the next iterate of the method is:

x i+1 = x i -αf (x i ) -1 f (x i ) . (1.44)
During this thesis, we did some work towards using a complete Newton's method, but we did not have time to develop a Line Search algorithm, so we ended up using the standard method from our solver: a single linearization of the nonlinear problem.

Linear solver

We want to be able to solve efficiently linear systems like (1.43). We additionally assume that A is an invertible matrix. Let us note the size of the linear system N . This size is quite large in our typical applications, but the linear system is sparse. This means that the coefficients of the matrix A are mostly zeros. This is because a coefficient in the matrix A corresponds to a link between two degrees of freedom. For the considered spatial discretisation methods, a cell only depends on a small number of neighbouring cells, and therefore a degree of freedom is not linked with most of the others. In practice, the stencil is generally reduced to the current cell and those that share a face with the current cell. This transcribes into lots of zeros in the matrix, and therefore its sparsity. Such sparse matrices can be seen in figure 1.5. Using the sparsity of the matrix is essential, as it would not be possible to store it as a dense matrix in the memory of a computer. Instead, some clever formats are used, such as the Compressed Sparse Row format [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. Some operations such as matrix-vector products are more efficiently done using such formats.

Many methods exist to solve linear problems. Some are even taught in school, such as the Gaussian elimination. Such methods are called direct methods, as they first do some work and then find directly the exact solution to the problem. But for problems with huge sizes such as the ones encountered in computational fluid dynamics, the amount of work is too much to be computed with today's means. It would take too much time as well as too much memory. Instead, one can use iterative methods [START_REF] Olshanskii | Iterative methods for linear systems[END_REF]. An iterative method starts from an initial guess x 0 and produces as it iterates a supposedly better estimation of the solution x n . The subscript n has nothing to do with the subscript used to identify the iterates of the time integration method: this happens at a given fixed time step. We can then decide when to stop the method, whether the solution estimate is good enough or the resolution is taking too much time. Iterative methods are the most well-fitted to solve large sparse linear problems, and they are the preferred solution in computational fluid dynamics.

CHAPTER 1. ANALYSIS OF EXISTING METHODS

Among the iterative methods are what could be called the "classic" methods, or relaxation. They are the Jacobi, Gauss-Seidel or Successive Over Relaxation methods. They decompose the matrix into A = M -N with M a matrix easily invertible. The choice of M depends on the method. Then, starting from a given x 0 , each step computes x n+1 = M -1 (N x n + b). Those methods are often deemed not efficient enough for computational fluid dynamics applications. They are often used, however, as preconditioning methods, as we will see later.

Another class of iterative methods is becoming the standard for computational fluid dynamics: the Krylov subspace methods. A Krylov subspace method projects the linear problem (1.43) on a smaller linear subspace called Krylov subspace. The obtained smaller system is then solved, much more easily. The cleverness resides in the fact that those Krylov subspaces are nested, and each iteration reuses the information obtained in the previous ones.

In the following, we note at iteration n the residual r n = b -Ax n . Practically, we keep n N so that the cost of the method stays reasonable, but this does not change the following. The corresponding Krylov subspaces are defined as:

K n (A, r 0 ) = Vect r 0 , Ar 0 , . . . , A n-1 r 0 .
(1.45)

Then, the next iterate in K n must satisfy a Petrov-Galerkin condition [START_REF] Simoncini | Recent computational developments in Krylov subspace methods for linear systems[END_REF]:

x n ∈ x 0 + K n (A, r 0 ) such as r n ⊥ L n (1.46)
where L n is a linear subspace with dimension n. For example, L n = K n corresponds to a Galerkin condition, and L n = AK n is a minimum residual condition.

To construct the growing Krylov subspace, one can use the Arnoldi iteration [START_REF] Trefethen | Numerical linear algebra[END_REF]. As a result, the matrix A is only used through matrix-vector products. Indeed, Krylov subspace methods do not require the matrix A to solve the linear system (1.43), just to know how to compute matrix-vector products. We will use this property later on.

There are many Krylov subspace methods that can solve a linear problem, characterised by the choice of L n . The minimal residual condition gives the Generalized Minimal Residual method or GMRES [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF], vastly used in computational fluid dynamics [START_REF] Franco | High-order matrix-free incompressible flow solvers with GPU acceleration and low-order refined preconditioners[END_REF] and other fields of numerical simulations [START_REF] Ernst | Why it is difficult to solve Helmholtz problems with classical iterative methods[END_REF][START_REF] Mercier | Fast nonlinear solvers in solid mechanics[END_REF]. Its main advantage is that even if its iterations may be slightly more expensive than other methods such as the Bi-CGSTAB method [START_REF] Van Der | Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems[END_REF][START_REF] Trefethen | Numerical linear algebra[END_REF], the residual norm is minimised and the error is then decreasing as the method iterates. This allows for some control of the residual norm. With the Bi-CGSTAB method for example, also available in our solver, the residual norms do not have to be decreasing, which can lead to chaotic convergence. The GMRES method already exists in our solver and is often used. Because it is discussed a lot in the literature, many variants and enhancements were developed [START_REF] Coulaud | Deflation and augmentation techniques in Krylov subspace methods for the solution of linear systems[END_REF][START_REF] Vasseur | Contribution to the study of efficient iterative methods for the numerical solution of partial differential equations[END_REF][START_REF] Jolivet | Block iterative methods and recycling for improved scalability of linear solvers[END_REF]. For those reasons, we decided to keep the GMRES method as the base of our linear solver.

The convergence of Krylov subspace methods, and more generally of iterative methods, depends on the linear system matrix. It is common knowledge that the convergence is linked to the condition number of the matrix [START_REF] Nevanlinna | How fast can iterative methods be[END_REF]. The condition number of an invertible matrix A is defined as:

κ (A) = A A -1 (1.47)
and so it depends on the chosen norm. A problem is said to be well-conditioned when κ = O (1), and ill-conditioned when κ 1. For the Euclidean norm • 2 ,

κ (A) = σ max σ min (1.48)
where σ min and σ max are the smallest and largest singular values of the matrix A. As a reminder, the singular values of A are the eigenvalues of A * A where A * notes the conjugate transpose of A. The relation (1.48) is often simplified as

κ (A) = |λ max | |λ min | (1.49)
with λ min and λ max the eigenvalue of A with the smallest and largest modulus. If this helps picture what the condition number stands for, this is only true for normal matrices: matrices that commute with their conjugate transpose. However, the matrices we face in our field have no reasons to be normal matrices. An intriguing result is found in [START_REF] Greenbaum | Any nonincreasing convergence curve is possible for GMRES[END_REF]. For a given decreasing convergence curve and a given spectrum, it is possible to construct a linear problem such as the convergence of GMRES follows the given curve and the matrix has the given spectrum. In other words, there are ill-conditioned matrices for which GMRES converges quickly and wellconditioned matrices for which GMRES will be slow. In particular, one can find a matrix with the best condition number possible, meaning 1, on which GMRES will make no progress until the last iteration. This is because of the nonnormality of the matrix [START_REF] Greenbaum | Matrices that generate the same Krylov residual spaces[END_REF][START_REF] Greenbaum | Any nonincreasing convergence curve is possible for GMRES[END_REF]. To handle it in the analysis of the convergence, one must not look at the spectrum but the pseudospectrum [START_REF] Trefethen | Pseudospectra of matrices[END_REF][START_REF] Trefethen | The Graduate Student's Guide to Numerical Analysis[END_REF]. This analysis is quite complex and is mostly done in an analytical context, not an industrial one.

Using Krylov subspace methods with nonnormal matrices and the study of their convergence is no simple task [START_REF] Liesen | Convergence analysis of Krylov subspace methods[END_REF][START_REF] Huhtanen | Aspects of nonnormality for iterative methods[END_REF]. As the matrices we encounter in our computational fluid dynamics problems are arbitrary, meaning mostly nonnormal, we decided not to study finely the spectrum of the operator.

Even if we will not look deeply into the spectrum of the matrix, we still help the linear solver with preconditioning. Preconditioning consists in multiplying the equation (1.43) with a preconditioning operator P to the left for left preconditioning:

P Ax = P b (1.50)
and to the right for right preconditioning:

AP x = b P x = x. (1.51)
The idea is to transform A into a matrix that is more easily invertible. The matrix to invert is now P A for left preconditioning and AP for right preconditioning. If the preconditioning operator is close to A -1 and is cheap to compute, this new matrix is close to Id and so is easily invertible. If P is too close to A -1 , computing the preconditioning will be as expensive as solving the original linear problem. If P is too easy to compute, meaning close to Id, the preconditioning will be useless. Choosing a preconditioner means making a compromise between those two extrema.

Just as Krylov subspace methods only need to compute matrix-vector products, and not to know the matrix, we do not need to know the preconditioning matrix, but only to know how to apply it on any given vector. To better understand how preconditioning works, let us take an example, inspired by example 35.2 of [START_REF] Trefethen | Numerical linear algebra[END_REF]. We take a square matrix of size 200 by 200, which is the sum of a diagonal part and a random perturbations part:

CHAPTER 1. ANALYSIS OF EXISTING METHODS

A = A diag + 1 2 √ N A rand where A diag,k = 2 sin 2kπ N -1 -1 + i cos 2kπ N -1 and A rand → N (0, 1) . (1.52)
The vector b is a vector of ones. The spectrum of A is displayed in figure 1.6 in blue. It is spread out around the origin. When we apply the Jacobi preconditioning, meaning the preconditioning matrix is the invert of the diagonal part, the spectrum is gathered around 1. This helps a lot GMRES, as can be seen in figure 1.6: it struggles to reduce the residual norm on the standard problem but does it easily on the preconditioned problem. This is only a dummy example, but it helps visualise the importance of preconditioning.

We need to use preconditioners when solving the linear system (1.43) with GMRES. One of the less sophisticated preconditioners is the Jacobi one. It is the one used in the example above: the preconditioning matrix is the invert of the diagonal part of A. Its advantage resides in its simplicity. It is inexpensive and easily scalable. On the other hand, it lacks efficiency, and that is why many preconditioners were developed in the literature. A slight variation consists in taking the diagonal blocs instead of the diagonal elements. This makes sense for our matrices, as they can be divided into blocs, each block corresponding to the degrees of freedom in a cell. The resulting preconditioner, called Block Jacobi, already exists in our solver and is often used as the standard preconditioner.

Others preconditioners exist in our solver. One consists in using the invert of the cell volumes as a diagonal matrix: the preconditioning matrix is a diagonal matrix, where each diagonal element is the inverse of the corresponding mesh cell volume. This was originally conceived to preserve the conservative properties through the GMRES solver. Indeed, as GMRES produces an approximation of the linear problem solution, the overall scheme may not be conservative. This preconditioner makes sure it is. More explanations on this statement can be found in appendix A. As it is even simpler than the Jacobi preconditioner, it shows poor performance. Some work was done towards Polynomial preconditioning: the preconditioning matrix P is equal to the partial Neumann series of Id -A, as it approximates A -1 [START_REF] Dubois | Approximating the inverse of a matrix for use in iterative algorithms on vector processors[END_REF]. Its major drawback is its computational cost, as it requires multiple matrix-vector products. For this reason, we decided not to look into it.

We were at first interested in some other preconditioners found in the literature, such as the Incomplete LU preconditioner. When taking the LU factorisation of the sparse matrix A, the triangular matrices L and U lose the sparsity pattern of A. This is troublesome as it is not possible to store dense matrices of the size of A. The Incomplete LU preconditioner, characterised by its index k, computes only the coefficients of the triangular factors that belong to the sparsity pattern of A k+1 . This way, the factorisation is incomplete in the sense that it does not recreate A, but it is still sparse and can be used as a preconditioner. The Incomplete LU factorisation is often used in computational fluid dynamics [START_REF] Liu | CFD high-order accurate scheme Jacobian-Free Newton Krylov method[END_REF][START_REF] Ahrabi | An implicit block ILU smoother for preconditioning of Newton-Krylov solvers with application in high-order stabilized finite-element methods[END_REF], but it was not a good fit for our solver, due to the difficulty in developing the method in an industrial-size solver, and because it relies directly on the matrix. Indeed, the matrix available to us is not precise, as we will explain later, so using it as a preconditioner may not work as well as expected.

As the preconditioning matrix P is not explicitly needed, and its effect should be near A -1 , one could take as a preconditioning procedure another Krylov subspace method: to apply the preconditioning to any vector v would mean solving Ax = v with a Krylov subspace method. However, applying a Krylov subspace method is not a linear operator. To handle this case, the GMRES algorithm must be modified into the Flexible Generalized minimal residual method, or FGMRES [START_REF] Saad | A flexible inner-outer preconditioned GMRES algorithm[END_REF][START_REF] Simoncini | Flexible inner-outer Krylov subspace methods[END_REF]. GMRES can be preconditioned with an inner GMRES, and as GMRES was already written in our solver this idea was simple to implement. This preconditioning does not need more information on the matrix A than the outer GMRES does: in particular it does not require knowing the matrix coefficients. This is interesting for some reasons that are discussed later. Finally, this method interests many scientists [START_REF] Coulaud | Deflation and augmentation techniques in Krylov subspace methods for the solution of linear systems[END_REF][START_REF] Vasseur | Contribution to the study of efficient iterative methods for the numerical solution of partial differential equations[END_REF] and has shown promising results in numerical simulation [START_REF] Pinel | A perturbed two-level preconditioner for the solution of three-dimensional heterogeneous Helmholtz problems with applications to geophysics[END_REF]. For all those reasons, we decided to add this method to our solver.

The GMRES method is often used with restarting: instead of keeping on iterating the method, one could stop it and start again. This is the Restarted GMRES method. The cost of one iteration increases as the method iterates. Restarting allows for some control of the Krylov subspace dimension, and therefore on computational cost and memory usage. The issue with restarting is that it can be harmful to the convergence. When the Krylov subspace generated on the restart is too close to the previously generated Krylov subspace, the residual norm may stagnate [START_REF] Simoncini | A new variant of restarted GMRES[END_REF]. A solution to this problem would be to reuse information between restarts. This is done with the Loose GMRES method [START_REF] Baker | A technique for accelerating the convergence of restarted GMRES[END_REF] or with Augmentation or Deflation techniques [START_REF] Chapman | Deflated and augmented Krylov subspace techniques[END_REF][START_REF] Morgan | GMRES with deflation restarting[END_REF][START_REF] Ramos | Projections, deflation and multigrid for nonsymmetric matrices[END_REF]. Augmentation and Deflation give in fact equivalent results [START_REF] Coulaud | Deflation and augmentation techniques in Krylov subspace methods for the solution of linear systems[END_REF]. Their idea is to recycle the spectral information acquired at the end of the cycle onto the next cycle. Recycling spectral information could even be done throughout multiple linear solve, during a nonlinear solve for example [START_REF] Gaul | Recycling Krylov subspace methods for sequences of linear systems -analysis and applications[END_REF]. Those techniques look promising but we did not have time to try them in our solver, as we focused on other points.

Evaluating the Jacobian matrices

If we can solve precisely linear problems, but the problem is not the one the nonlinear solver wants, it may hurt the nonlinear solver convergence. We then need to be able to get the right linear problem from the nonlinear one. At this point, we know how to evaluate the function f from equation (1.41), which is often expressed as a linear combination of previous states and G evaluations, where G is the function introduced in equation (1.10), the function from the starting ordinary differential equation. On the other hand, computing its derivative with regard to x is a different story. As this derivative is the matrix used when solving the linear problem, it is crucial to have a good representation of it. This derivative f (x) uses the Jacobian matrix of the function G. As the function G comes from the spatial discretisation method applied to the original partial differential equations (1.1), it is quite complex to evaluate this Jacobian. The underlying algorithm is hard to fully understand and write, and the numerical evaluation is usually expensive. Furthermore, our software keeps evolving, and models are constantly added and modified. It would require constant work to maintain the Jacobian matrix computation. Computing by hand the exact Jacobian matrix would amount to too much work in our industrial software. We must use other alternatives to get the Jacobian matrix we need for Newton's method.

An idea would be to use Automatic Differentiation [START_REF] Griewank | Evaluating derivatives[END_REF]. This means to give the source code of the G function to some software [START_REF] Hascoët | The Tapenade automatic differentiation tool: principles, model, and specification[END_REF], which gives in return a way to compute its Jacobian matrix. One advantage of this method is that the cost of derivating the function is done only once, at software compilation. After that, computing the Jacobian matrix amount to calling a function. Another advantage is that the given Jacobian is supposedly exact, contrary to some alternatives we will discuss later. For those reasons, Automatic Differentiation is today being used in actual computational fluid dynamics software [START_REF] Bilanceri | Comparison of explicit and implicit time advancing in the simulation of a 2D sediment transport problem[END_REF][START_REF] Kenway | Effective adjoint approaches for Computational Fluid Dynamics[END_REF]. In our software, using Automatic Differentiation did seem too hard and therefore we looked at other methods instead.

A possible idea is to use an approximation of the Jacobian matrix that is inexpensive to compute. We said that the function G comes from the second-order or higher-order Finite Volume method used as the spatial discretisation method. One can then take the function G 1 given by the first-order corresponding method, and use its Jacobian matrix instead. Indeed, computing G 1 is inexpensive in comparison to G, and the same goes for the corresponding Jacobian matrices. Going further, one could also approximate the Jacobian matrix of G 1 : when G uses some complex models such as turbulence models, it is often decided to not include those models' contributions to the Jacobian matrix, for complexity and stability reasons [START_REF] Content | Coupled/uncoupled solutions of RANS equations using a Jacobian-free Newton-Krylov method[END_REF]. This is what is done originally in our in-house solver: using a cheap low-order approximation of the Jacobian matrix.

When we decided to use Krylov subspace methods as our linear solver, we highlighted the fact that they only need to know how to compute matrix-vector products. Using this, and the fact that the matrix is a Jacobian matrix, one could approximate the matrix-vector product of f (x) with a vector v by:

f (x) v ≈ f (x + εv) -f (x) ε (1.53)
with a scalar parameter ε discussed below. This is the finite difference approximation of the Jacobian matrix-vector product. A quick analysis shows that the error on this approximation is o (ε). Usually, the value of f (x) is previously computed, so one evaluation of f is required for one matrix-vector product. One could also use the centered finite difference approximation:

f (x) v ≈ f (x + εv) -f (x -εv) 2ε (1.54)
that gives a smaller error of o ε 2 , but it requires two f evaluations, so it is twice as expensive as the first-order approximation. Therefore, we will use the first-order approximation.

Using this approximation, it is possible to recover the full Jacobian matrix. If the approximation (1.53) is applied to each vector of the canonical basis as the vector v, one can gather each column of the Jacobian matrix. This gives what is often called the finite difference Jacobian matrix. Instead of computing a matrix-vector product for each direction, a technique consists of computing independent directions beforehand to get a colouring of the matrix: two directions are of the same colour if they are independent through the matrix. Then, a smaller number of matrix-vector product approximations are required: as many as there are colours. This is often called the finite difference Jacobian matrix with colouring [START_REF] Gebremedhin | What color is your Jacobian? Graph coloring for computing derivatives[END_REF]. Those techniques are still not well fitted for our solver. As we work with large dimensions, computing the Jacobian matrix takes time. Furthermore, we explained that we do not need to know explicitly the matrix. That is why we will use the approximation (1.53) each time we need a matrix-vector product instead of computing the Jacobian matrix first.

We still have not talked about the parameter ε introduced in equation (1.53). As it represents the size of the step made to approximate a derivative, it should be small. But taking it too small leads to roundoff errors so it must be chosen carefully [START_REF] Knoll | Jacobian-free Newton-Krylov methods: a survey of approaches and applications[END_REF]. We decided to look at some strategies in the choice of ε, and this work will be presented in a later part. Now the choices we made for our nonlinear and linear solve strategies are starting to make sense. Working with different Jacobian matrices is not an issue, as they are not computed. The Krylov subspace method GMRES does not need the Jacobian matrix, only to compute matrixvector products, contrary to other iterative methods such as the Gauss-Seidel method. The flexible preconditioning also does not need the Jacobian matrix for the same reason, contrary to other types of preconditioners such as ILU preconditioners.

Overall, using Newton's method to solve the nonlinear problem and a Krylov subspace method to solve the linear problem while using a matrix-free method is known as the Jacobian-Free Newton-Krylov method, or JFNK. This methods and its variants are as of today still discussed [START_REF] An | On finite difference approximation of a matrix-vector product in the Jacobian-free Newton-Krylov method[END_REF][START_REF] Turpault | Modélisation, approximation numérique et applications du transfert radiatif en déséquilibre spectral couplé avec l'hydrodynamique[END_REF] and used on actual computational fluid dynamics solvers [START_REF] Liu | CFD high-order accurate scheme Jacobian-Free Newton Krylov method[END_REF][START_REF] Franco | High-order matrix-free incompressible flow solvers with GPU acceleration and low-order refined preconditioners[END_REF]. We are interested in the JFNK method as we hope it will give a better representation of the Jacobian matrix than the low-order one already in use in our solver. We hope that a better Jacobian matrix will lead to a better convergence, in particular when some models are ignored in the classic Jacobian matrix, such as turbulence models.

Furthermore, CEDRE works with multiple distinct solvers. When a computation uses multiple solvers, two, for example, the function f from equation (1.41) is in fact (f 1 , f 2 ) and x is (x 1 , x 2 ). The Jacobian matrix is then:

f (x) =   ∂f1 ∂x1 ∂f1 ∂x2 ∂f2 ∂x1 ∂f2 ∂x2   .
(1.55)

CHAPTER 1. ANALYSIS OF EXISTING METHODS

The cross-term Jacobian matrices, ∂fi ∂xj with i = j, are usually hard to compute, both analytically and computationally, and are the main obstacle to a fully implicit solver. Instead, in CEDRE, each solver iterates independently to the others, with some exchanges between iterations. Using the matrix-free approximation and the JFNK method, those cross-solver terms are taken into account. This choice is then perfectly aligned with the direction our solver is aiming: towards a better implicit cross-solver mechanism. Some specific situations are known to be troublesome for the JFNK method, such as shocks, reaction fronts or discontinuities arising from high-order advection schemes [START_REF] Knoll | Jacobian-free Newton-Krylov methods: a survey of approaches and applications[END_REF]. Unfortunately, those features are precisely at the center of our solver. We are hoping that the difficulties of the method will be overbalanced by its advantages: getting a better convergence, allowing a fully coupled implicit solve and allowing an implicit integration method for new models. It means that we do not expect the method to enhance the robustness of the solver, but we hope that because it uses a better Jacobian matrix than the one already in use it will help the convergence. And once again, the second advantage is that the matrix-free method would allow us to do coupled implicit time integration. We will come back to the third advantage later in the application's part.

In this section, we described how an implicit method operates: the method itself produces one or several nonlinear problems at each time step. Those problems are solved with a nonlinear solver that uses linearisation to transform them into linear systems, which are then solved with a linear solver like a preconditioned GMRES method in the case of our solver. We explained that this linearisation in CEDRE is done by using Jacobian matrices that are poorly approximated. It means that we do not solve the correct linear systems, which hinders the nonlinear resolution and the time integration method ends up finding an incorrect increment. The Jacobian-Free Newton-Krylov method, which is compatible with other algorithms already available in CEDRE such as the GMRES method, should rectify this issue. For this reason, we will continue this work with this method for our time integration method.

Chapter 2

Description of the development workflow in CEDRE

Résumé du chapitre : Description des développements dans CEDRE Le but de ce chapitre est de décrire les développements réalisés dans CEDRE afin d'implémenter le mécanisme "sans matrice" de la méthode JFNK. Il s'agit de développer astucieusement ce mécanisme afin d'obtenir la méthode souhaitée tout en garantissant l'efficacité et la portabilité du code produit.

Dans ce chapitre, nous décrivons un peu plus en détail les différentes briques qui composent CEDRE. Nous donnons ensuite les idées utilisées afin de réaliser les développements de la manière la plus générique possible. La nécessité de donner ces détails vient du fait que CEDRE est un code particulièrement massif, qui se prête mal aux développements. Afin d'obtenir le résultat souhaité, nous avons donc utilisé des notions modernes du langage Fortran que nous détaillons dans ce chapitre. Nous précisons ensuite un point précis de l'implémentation de la méthode JFNK, qui est le choix du paramètre de différenciation ε, en apportant une correction par rapport aux choix principalement utilisés par le reste de la communauté. Nous montrons enfin comment ces nouveaux développements s'articulent autour des mécanismes déjà existants, pour vérifier la compatibilité de l'ensemble des méthodes.

In this brief chapter, we will discuss the details of how we implemented selected methods in CEDRE. This does not constitute research work, but it ended up being a large part of the work done during this thesis. It is also not without interest, as we sometimes used advanced features in order to implement what we set out to do.

CHAPTER 2. DESCRIPTION OF THE DEVELOPMENT WORKFLOW IN CEDRE

General description of CEDRE

The software system CEDRE gather several solvers to solve problems in the field of multiphysics [START_REF] Refloch | CEDRE software[END_REF]. Each solver is dedicated to a given model. As of today, there are seven solvers embedded in CEDRE:

• CHARME, the fluid solver, for compressible multifluid and reactive flow, with RANS or LES turbulence models

• SPIREE, the dispersed phase solver using an Eulerian framework

• SPARTE, the dispersed phase solver using a Lagrangian framework

• ASTRE, the radiation solver using a Monte Carlo method

• REA, the radiation solver using a discrete ordinates method

• FILM, for shallow water equations, used to model ice accretion

• ACACIA, the conduction solver, for heat transfer in solids.

Combining different solvers, CEDRE is able to simulate multiphysics phenomena numerically. Using those solvers, CEDRE applications go from aerodynamics to aeroacoustics, aerothermodynamics, combustion, icing, etc. The solver is coupled either through boundary conditions for example in a thermal interaction at a fluid-structure contact or inside the computational domain for example in the case of mass and energy transfer between dispersed phases and the main flow. The coupling can either be one-way or two-way, depending on the user's choice. Each solver is integrated in time separately, and the coupling consists of some data exchange between iterations: it is an explicit coupling.

Some functionalities common to multiple solvers exist outside the solver in helper libraries. For instance:

• ASSEMBLAGE acts as the conductor by handling the overall simulation, telling the solvers what to do and when to do it, when to exchange data and with which other coupled solver

• BIBCEDRE contains tools for geometrical operations, linear algebra methods, mesh handling, parallel communications and other general functionalities

• THERMOLIB is used to compute the different thermophysical properties such as heat capacities or to compute the closure relations that can be complex in our applications.

Those libraries are the main engines of the solver. CEDRE does a computation with CHARME, it is in fact ASSEMBLAGE that tells CHARME when to initialise itself, to do a step of the time integration method, etc. When CHARME needs to compute the ordinary differential equation function, it is BIBCEDRE that does the higher-order polynomial reconstruction of the solution so that the Riemann solver can compute the interface flux. When CHARME needs to compute the primitive variables such as the pressure and temperature from the conservative variables, it is THERMOLIB that solves a nonlinear problem with Newton's method to compute the result. Now that we defined the existing framework for the software system CEDRE, we can develop new functionalities. Despite allowing some flexibility in the programming language, most of CEDRE is written in Fortran. We decided to keep working with Fortran to help with the integration of our work. As CEDRE is used by industrial clients, and as they rely on their own supercomputers, we need to limit ourselves to Fortran 2003 standards, so as to ensure compatibility.

Implementation details

In this thesis, we focused on the most used solver: CHARME. Indeed, not only is it the most used, but other solvers use it as a base for many applications. For example, when simulating ice accretion around a wing profile, a standard methodology with CEDRE is to first get the base aerodynamic flow with CHARME and then compute the ice particles with SPIREE or SPARTE. Working with the solver CHARME was the way to benefit the most from our work. Even if during this thesis we only worked on CHARME, we always kept in mind that the finality was multiphysics simulations using multiple solvers. That is why we tried to develop generic functionalities so that they could be easily imported to other solvers, provided the developers of said solvers wanted to use them. The same reason was also used as a criterium in our choices, as was explained previously. Choosing the Jacobian-Free Newton-Krylov method goes towards fully implicit coupling between solvers, instead of the explicit coupling existing today.

FGMRES

In order for our work to be usable in every other solver, we had to work on the common library BIBCEDRE. When the implicit Euler method of CHARME needs to solve a linear problem, it uses BIBCEDRE. It contains everything needed to solve linear problems, such as GMRES and preconditioners. A linear problem is stored in BIBCEDRE as the Fortran derived type type_sys. In order to add Flexible preconditioning to the existing GMRES, we added a pointer to an inner instance of type_sys inside of type_sys, so that the linear system and its corresponding solver may use an inner solver for an inner problem: type type_sys ! Inner linear system and solver type(type_sys), pointer :: sys_int => null()

... ! Additional data end type type_sys

This way, when we need to apply the preconditioner during a GMRES iteration, we can use the inner type_sys instance to call the inner GMRES. Furthermore, having a pointer to an inner instance allows for more freedom for the inner solver. One could for instance use multiple depths of preconditioning and have the inner GMRES also be an FGMRES method, preconditioned by another GMRES, etc.

Matrix-free

Sparse matrices are stored in an in-house format, using an array for the diagonal blocks, another one for the extra-diagonal blocks and a third one to index the extra-diagonal blocks. Matrix-vector products are made inside BIBCEDRE to handle this matrix format, with the routine: This way, when the client solver creates an instance sys of type_sys, it can choose how to evaluate matrix-vector products by setting the procedure pointer sys%gmvec_ctx%gmvec. It can for example point to the already existing routine to use the classical matrix-vector product, but it also can use a custom routine that implements the approximation (1.53). Furthermore, the client solver can use polymorphism and create an extended type of type_gmvec_ctx in order to store additional data in the context. This is what is done by the solver CHARME, as it does need additional data to approximate the Jacobian matrix-vector product. Finally, with this implementation, any solver that wants to use the approximation (1.53) just needs to write the corresponding routine and set the context accordingly. Then, a user can choose at execution time whether to use the standard Jacobian matrix or the matrix-free method.

subroutine gmvec(sys,
The computational cost of the matrix-vector approximation is approximately the cost of one function evaluation. It requires some additional operations, such as vector subtractions, and vector scaling, but those operations are way less expensive. Using a matrix-free algorithm is clearly less expensive in terms of memory, as the matrix does not need to be stored at all. Some even claim that it can lead to a less expensive method in terms of time [START_REF] An | On finite difference approximation of a matrix-vector product in the Jacobian-free Newton-Krylov method[END_REF]. This happens when the alternative, which is computing the Jacobian matrix, is an expensive task. In CHARME, the Jacobian matrix that is traditionally used is inexpensive to compute, with a cost similar to a function evaluation. The issue is that the function evaluation itself is expensive in our solver. This is troublesome, and the issue is currently addressed by the team of developers in charge of the solver. In the meantime, it means that the cost of the matrix-free approximation will suffer even more. Practically, using it increases the computational time of one iteration. We now hope that as it will not improve the iteration speed, it will improve the convergence of our solver. Also, as we said the team responsible for CEDRE is working to improve the speed of the function evaluation. If they succeed, our method will then benefit greatly from their work.

Strategy for the choice of ε

When we introduced the approximation (1.53) it introduced a new parameter ε. It is easy to check that the truncation error on this approximation decreases linearly with regard to ε. It is natural to take a small value for ε. But unfortunately, when working with floating-point arithmetic, dividing by a small ε introduces roundoff errors [START_REF] Martel | Semantics of roundoff error propagation in finite precision calculations[END_REF]. This parameter needs to balance truncation and roundoff errors. We need to decide on a strategy for the choice of epsilon. One could take for example ε = √ ε mach where ε mach is the machine epsilon: around 10 -6 for single precision and 10 -15 for double precision. This choice is often discarded as it is deemed too simplistic. Instead, works from the literature tend to use the same few options [START_REF] Park | On physics-based preconditioning of the Navier-Stokes equations[END_REF][START_REF] Liu | CFD high-order accurate scheme Jacobian-Free Newton Krylov method[END_REF][START_REF] Abhyankar | PETSc/TS: a modern scalable ODE/DAE solver library[END_REF]] that come from [START_REF] Pernice | NITSOL: a Newton iterative solver for nonlinear systems[END_REF] and [START_REF] Dennis | Numerical methods for unconstrained optimization and nonlinear equations[END_REF]. Those options are well described in [START_REF] Knoll | Jacobian-free Newton-Krylov methods: a survey of approaches and applications[END_REF]. In particular, the one that we encounter the most is the one from [START_REF] Pernice | NITSOL: a Newton iterative solver for nonlinear systems[END_REF]:

ε wp = ε 0 (1 + x 2 ) v 2 (2.1)
using the same x and v as in equation (1.53). Here ε 0 is the estimated relative error in function evaluation. A reason this choice is so popular is that apart from this ε 0 value, it does not require any user input. Furthermore, ε 0 is often simply set to machine epsilon ε mach .

To analyse this strategy in the choice of ε, we do the following numerical experiment. We consider the one-dimensional Burgers' equation over a regular periodic mesh made of 10 cells (or segments in the one-dimensional space). The function f is taken as the right-hand side of equation (1.5) when solving the Burgers' equation, using a first-order Finite Volume method as the spatial discretisation method, itself using an exact Riemann solver (Godunov's scheme). This function was chosen because it is a nonlinear conservative equation, often viewed in computational fluid dynamics as a simplified version of Euler equations and it contains most aspects of nonlinear hyperbolic equations. The vectors x and v are 10 + r 1 and 0.1 (2r 2 -1) where r 1 and r 2 are random vectors following a uniform distribution on [0, 1[, obtained with the Python package NumPy: numpy.random.random.

The experiment result can be seen in figure 2.1. The relative error in the approximation is shown as a function of ε with small grey crosses. On the right part of the figure, the error decreases linearly with regard to ε. This corresponds to a truncation error-dominated region. On the left part, the error increases as ε decreases. This corresponds to a roundoff error-dominated region. Also, this increase is not smooth as the linear decrease, because roundoff errors tend to produce more chaotic results. The choice of ε from [START_REF] Pernice | NITSOL: a Newton iterative solver for nonlinear systems[END_REF] is also shown in the figure, as the large blue dot. As can be seen, it falls into a low error region, not too much on the left, not too much on the right.

The issue we noticed is the following. Some references in the literature do not specify the norm used in equation (2.1). As it is most of the time the Euclidian norm • 2 , or 2-norm, we can assume that it is also the case when it is not specified. However, this means that in our Figure 2.1: Error in the Jacobian matrix-vector product approximation as a function of ε, and in particular for the most popular choice ε wp . The function is the right-hand side of equation (1.5) from a Finite Volume method using an exact Riemann solver for Burgers' equation, over a 10-cell one-dimensional regular mesh. example the value of ε wp will depend on the vector size. With our example, if one increases the number of cells in the mesh the typical size of the vector components stays roughly the same, so the shape of the error as a function of ε should not change much from the one in figure 2.1. But if the vector components size stays the same, the 2-norm does increases with the vector dimension. One can even see that when the dimension is N 1, ε wp ∼ N -1/4 . Having ε to depend on N did not make sense to us, as we are working with possibly large vectors in our industrial applications. That is why we decided to use a variation from what is currently found in the literature: we will use the norm • 2 / √ N that can be seen as a scaled 2-norm. We then get a new strategy for the choice of ε that we note ε wp,N :

ε wp,N = ε 0 1 + x 2 / √ N v 2 / √ N . (2.2)
We compared the two strategies on the same test case while increasing the number of cells from 10 to 100, 1000, ..., up to 10 7 . Figure 2.2 shows the result of this experiment. On the left, we see that the shape of the relative error as a function of ε when N = 10 7 is similar to when N = 10, albeit the roundoff error-dominated region is more regular. In particular, the ideal trade-off between the two error types did not move a lot. The value of ε wp decreases as expected: this translates as the fact that the blue dot moved to the left. When looking at the error as a function of the dimension N , one can even recover the -1/4 expected slope. In the right figure, there is in fact a 1/4 slope but the error is inversely proportional to the value of ε so if the error behaves as N 1/4 , then ε wp behaves as N -1/4 . Because of the shape of the error as a function of ε, if ε wp changes with N it will eventually end up increasing the error. This is an undesired feature. Our new choice ε wp,N on the other hand does not change a lot when N increases. Therefore, the error level stays the same no matter the dimension.

The same numerical experiment was made on a more complex case: the one-dimensional Euler equations. Here, the function f corresponds to the right-hand side of equation (1.5) given by a centred Finite Volume method used as the spatial discretisation method, in which the interface flux is defined as an average of left and right fluxes. This means the Riemann solver just averages the left and right fluxes. This spatial method is known for its instability when used in an actual solver, but it is useful in this experiment as the analytic Jacobian matrix is easy to derive, contrary to methods using more complex Riemann solvers. This physical model assigns 3 degrees of freedom in each cell. We will use primitive variables. The vector x corresponds to a uniform density of 1kg m -3 , a velocity equal to a sine making one period over the mesh and of amplitude 10m s -1 , and uniform pressure of 10 5 Pa. The vector v is a random vector in [0, 1] as before, where the first component is scaled by 10 -3 , the second by 10 -2 , and the third by 10 2 , in order to impose 10 -3 relative perturbations. Figure 2.3 shows the corresponding results. As for the previous experiment, ε wp depends on the dimension N , and then the associated error ends up increasing as N grows. Our correction ε wp,N does not exhibit the same drawback. Even if it does not fall at the bottom of the error curve, it at least does not lead to a larger error.

In this last case, we see that at the beginning the two errors were close with ε wp,N being better than ε wp . In the first one, their starting positions were a bit different and this time ε wp was the better one. This is largely due to the randomness of this analysis: the one used to construct the x and v vectors. In fact, it would be unwise to draw a conclusion from the position of the points in figures 2.2 and 2.3. Changing the choice of vectors, of the function f , or even the randomness in the vectors is enough to modify their position. For example, we can not conclude anything from the fact that ε wp,N is almost at the bottom of the error curve in figure 2.2. It may as well be a bit more on the side like in figure 2.3. What we can use from this analysis, however, is the tendencies that choices of ε show. The main conclusion is then that with our strategy we removed a dependency of the relative error on the dimension, which makes sense as this dependency is not expected a priori.

The issue with this analysis is that it relies on extremely simple examples. The function f in our actual applications is in fact much more complicated than the two we used here. But in order to do this analysis we need to be able to compute the relative error, and this means we need to be able to compute a Jacobian matrix-vector product analytically. Unfortunately, this is not possible for our applications. If it was, we would not even need to introduce the approximation (1.53) and therefore do this analysis. Using one-dimensional Burgers' equation and one-dimensional Euler equations with a simple Riemann solver was then a good compromise. Those equations share similarities with the ones we work with in our solver, as they are nonlinear hyperbolic equations, but they are still simple enough so that they can be used here.

Linear system normalisation

An available option in BIBCEDRE for the linear solver is normalisation. It consists in taking an invertible diagonal matrix N and rewriting the linear system (1.43) as:

Ãx = b with à = N -1 AN , x = N -1 x and b = N -1 b (2.3)
One can check that it is equivalent to the original linear system. The diagonal coefficients of the matrix N are estimates of the order of magnitude of the corresponding degrees of freedom. This estimation can be computed at each time step depending on the current state, or can simply be constant, equal to some reference values, depending on the choice of the user. When working with Euler equations, for example, the numerical values corresponding to the energy degrees of freedom are usually many orders of magnitudes over other degrees of freedom. The goal of this normalisation is to bring back each vector component at the same level in the linear system. It can also be seen as a combination of a left preconditioning by N -1 and a right preconditioning by N . This normalisation is often used by our solver's users, so we need to ensure compatibility with our matrix-free method. We can adapt the approximation (1.53) to account for the normalisation as:

N -1 f (x) N v ≈ N -1 f (x + εN v) -f (x) ε . (2.4)
We also need to adapt our strategy for the choice of ε. With the normalisation, we see that εN sort of plays the role of ε. It then makes sense to take εN equal to the original choice of ε, and therefore we take:

ε = ε 0 (1 + x ) N v . ( 2.5)
We take as a matrix norm the one induced by the vector norm. This means that:

N = sup v =0 N v v . ( 2.6)
When using the 2-norm, the corresponding norm is N = max N ii as N is diagonal with strictly positive coefficients. When using the scaled 2-norm discussed in the previous part, the corresponding norm is also the same. This is nice because whatever strategy we use in the choice of ε, it does not interfere with the handling of normalisation.

Preconditioning

Our matrix-free implementation also needs to be compatible with the already existing preconditioning from BIBCEDRE. There are two ways to precondition the linear system in BIBCEDRE.

• The outer left preconditioning is applied once before the linear solve. It changes the numerical values of the matrix coefficients and the right-hand side. Once it is applied, the resolution can continue as if nothing had happened.

• The inner right preconditioning is applied before each matrix-vector product, and once more at the end of the solve.

Handling the inner preconditioning with our matrix-free implementation is easy as those two operations are independent of each other. The outer preconditioning relies on the fact that the preconditioners available are simple, and therefore it is inexpensive to apply it to the whole matrix, and less expensive than applying after each matrix-vector product. This advantage is lost with the matrix-free method. The left preconditioning is still applied once to the right-hand side before the solve, but it is also applied after each matrix-vector approximation.

Local time-stepping

We said earlier that local time-stepping is used in CEDRE to accelerate the convergence on steady problems. What it does is that each cell of the mesh moves forward in time with its own time step, independently of other cells. This local time step is expressed as a fraction of the global time step. There are several options a user can choose from to set the time-stepping strategy. It can use a CFL number based option and compute the corresponding local time step so that the CFL number local to the cell is not too large. The preferred and default option is based on the explicit increment and some reference values. It first estimates what the solution increment will be by using the explicit increment and possibly the reference values, depending on the choice of method, and compute a local time step so that the increment local to the cell is not too large. This local time-stepping strategy adds various parameters that need to be tuned, in the form of scalar coefficients, whether or not to apply the local time-stepping to the boundaries, whether to use conservative or primitive variables, etc. In the end, what local time-stepping methods do is artificially modify the volume of the cells: increasing a cell volume corresponds to increasing its inertia, and therefore scaling the time step down. The volumes of the cells appear in the linear system matrix, in the form of a diagonal "mass" matrix. By using the scaled mass matrix, the matrix-free approximation does account for local time-stepping. With this last feature, we ensured that the matrix-free approximation is compatible with the already existing methods. This chapter explained in a simplified manner some implementation detail of our solver. It is important to remind oneself that CEDRE is an industrial-sized software system, compared to smaller solvers developed for academic purposes. The scale of the solver was a limiting factor in this thesis. We often wanted to try some new scientific idea but were faced with many troubles that are inherent to the solver size. At some point at the early stage of this thesis, for example, we wanted to improve the quality of the linear solve by using advanced solvers and preconditioners, with the help of the PETSc library [START_REF] Balay | PETSc Web page[END_REF][START_REF] Balay | PETSc/TAO users manual[END_REF][START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF]. Using the library in CEDRE ended up being impractical, or at least would have taken too much time that we could not afford, and so we had to change direction. In the end, a solver of this size is not an ideal sandbox to try methods that are too innovative. A good amount of experience is required to master the various models and methods, an amount I did not have at the beginning of this thesis. For example, the spatial discretisation methods come with a lot of parameters: the kind of method, the slope limiter, the Riemann solver, and many more that are specific to CEDRE. The same goes for the time integration with for instance the local time-stepping or the linear system normalisation. A single parameter may be deterministic in the success of a computation. This forces us to allocate a quite large amount of time to every computation. Moreover, because of the large number of developers that are working on the solver, and the fact that many developments are small tweaks from occasional developers such as interns and PhD students, fixing bugs also took a decent amount of time. Overall, this explains why we decided to implement the solutions previously discussed. They are already well known from the literature and feel like a natural step up from the already existing methods of CEDRE. It is good to mention that many difficulties arise from the solver, and are inherent to it, mostly due to its overall size. We would probably have taken a different route if we had worked with a different solver.

Without getting into too much detail, this chapter explained how the selected methods and modifications were implemented into the solver CEDRE. Besides some minor developments, the main feature is the matrix-free implementation. As it was explained, this implementation is generic so that it can be used with any solver that wants to solve a linear problem with the library BIBCEDRE. The idea was also that if multiple solvers use this implementation, it becomes easy to do a coupled implicit solve. Without it, because of the structure of CEDRE with the library ASSEMBLAGE conducting the solvers, a coupled implicit solve is not possible without reworking the basic structure of CEDRE. During this thesis, however, because we were limited in the time allowed to us, the implementation was only used for the reactive Navier-Stokes solver CHARME. This is why we had to compromise and develop only for the CHARME solver but in a generic way.

Chapter 3

Analysis of the JFNK method in CEDRE

Résumé du chapitre : Analyse de la méthode JFNK dans CEDRE Le but de ce chapitre est de vérifier les gains apportés par l'ajout de la méthode JFNK dans le solveur CEDRE. Il s'agit de comparer les performances de cette nouvelle méthode avec les méthodes préexistantes en matière de stabilité, robustesse et rapidité. Nous adaptons pour cela le point de vue d'un utilisateur, en considérant ce qui l'intéresse lorsqu'il réalise une simulation stationnaire.

Pour cela, nous comparons tout d'abord la méthode JFNK avec la méthode utilisant explicitement la Jacobienne. Cette comparaison se traduit principalement par la comparaison de la convergence des deux méthodes à travers la norme des résidus. Ainsi, nous réalisons cette comparaison sur une succession de cas tests, de complexité croissante, sélectionnés car ils représentent différentes facettes du solveur.

Nous commençons par un calcul purement aérodynamique d'un profil d'aile en deux dimensions dans un écoulement turbulent, qui utilise une modélisation des effets turbulents de la couche limite. Puis, nous utilisons le même profil mais cette fois en incidence, avec un maillage bien plus fin qui capture les effets de couche limite turbulente, plutôt que de la modéliser.

Nous regardons ensuite un cas de rentrée atmosphérique avec une sphère placée dans un écoulement à haute énergie. Cela se traduit par un très fort choc, ainsi qu'une zone de déséquilibre thermodynamique où ont lieu d'intenses réactions chimiques. Dans un premier temps, nous regardons comme précédemment les deux méthodes pour comparer leur convergence. Pour finir, nous utilisons un modèle de fluide récemment ajouté au solveur afin de représenter plus finement le déséquilibre thermodynamique en caractérisant l'écoulement par plusieurs températures distinctes. Avec les développements réalisés actuellement sur ce nouveau modèle, il n'est pas possible d'utiliser la méthode implicite préexistante, et les utilisateurs sont forcés d'employer des méthodes explicites. Nous montrons donc que la méthode JFNK peut elle être utilisée en l'état, sans nécessiter de développements supplémentaires, et permet d'accélérer les calculs ainsi que d'obtenir une meilleure convergence.

CHAPTER 3. ANALYSIS OF THE JFNK METHOD IN CEDRE

In a previous chapter, we identified some methods from the literature that we wanted to use in our solver CEDRE, and some others from CEDRE that we wanted to improve. In another chapter, we discussed the practical implementation of said methods in the solver. The goal of this chapter is now to test those methods on several applications to comment on the choices we made. We need to define test cases that represent well enough target applications so we can comment on the performances of our choices.

Comparison between matrix-free formulation and Jacobian matrix approximation

From the previous analysis and implementation, the main addition to the solver is the Jacobian-Free Newton-Krylov method, and in particular the matrix-free approach. In this part, the new method that uses the matrix-free approach will be compared with the implicit Euler method as the interest is focused on implicit time integration schemes. The traditional method linearises the equation that the implicit Euler method produces, approximates the Jacobian matrix using the Jacobian matrix of the corresponding first-order scheme, and then solves the linear system with the Krylov subspace method GMRES. In order to understand the impact of a better Jacobian matrix, the only difference is how the Jacobian matrix is handled. The new method will then work in a similar way, except the matrix used in the linear solver is not actually computed, but the matrix-vector products are approximated using equation (1.53).

Turbulent transonic airfoil

Definition of the test case

The first application is a typical aerodynamics test case. It is a two-dimensional simulation of the flow around an RAE 2822 wing profile. The fluid is standard air assumed to be a perfect gas. The Mach number is taken equal to 0.75, the chord is equal to 1m, the angle of attack is 0°and we use the atmospheric conditions at 10km. This gives a laminar Reynolds number of 6.5 × 10 6 . This first test case is chosen for multiple reasons. Firstly, it is a simple case in the field of computational fluid dynamics. It is a standard aerodynamics case, with a small mesh in comparison to many other three-dimensional cases. This allows testing our methods inexpensively. Secondly, this case belongs to the tutorial suite of our solver. It means that it is already well mastered by the team. Thirdly, even if it is only a standard aerodynamics case it still has some stiff features, such as turbulence modelling and a shock. Finally, it is a standard test case for turbulence modelling validation. Therefore there are many references in the literature using this case.

The mesh used for this simulation is an unstructured hybrid mesh made of triangles and quadrangles. Parts of it can be seen in figure 3.1. Cell sizes range from 2.5m far from the airfoil and 100µm at the wall. At the wall, there is a C-shaped layer of regular cells. This helps better capture boundary layer effects near the profile and the wake. Also, under those conditions, a shock is expected to develop on the upper part. Special treatment such as refinement was applied to the mesh at the expected shock location. The model used is the Reynolds Averaged Navier-Stokes equations, or RANS. Simply put, every scale of the turbulence is modelled. The Spalart-Allmaras turbulence model is chosen to close the RANS equation. It is well known for being one of the simplest turbulence models, which is fine to set up a simple first test case. Figure 3.1 shows that the mesh near the wall boundary is not particularly fine, and therefore it is not fine enough to compute the boundary layer. It is indeed because this computation uses a turbulence wall model that replaces the standard no-slip boundary condition with a more complex relationship between the variables and their derivatives. The physics of the boundary layer is introduced in the model, which limits the stiffness of the system. Such models are known to be troublesome with some more complex turbulence models such as the k -model but are used with others such as the one we use here. Turbulence wall modelling may be nonconventional in aerodynamic simulations but it is a commonly used feature of our solver, so it is interesting to see how our new method behaves on such test cases. The spatial discretisation method is a second-order Finite Volume method, using the HLLC Riemann solver and Multislope method [START_REF] Le Touze | Multislope MUSCL method for general unstructured meshes[END_REF] with a Van-Leer slope limiter. Local time-stepping is used to speed up the convergence.

COMPARISON BETWEEN MATRIX-FREE AND MATRIX APPROXIMATION

Analysis of the results

We are now going to compare the results from the two different simulations. It is first worth noting that before trying the JFNK method, the computation needs to start using the traditional method. As explained before, the matrix-free approximation is not ideal to handle discontinuities such as shocks, and a shock is indeed present in this computation. Only some iterations are needed, just to start the development of the boundary layer and approximately place the shock. After that, the computation can be continued, with the traditional method on one side and with the new method on the other.

The global flow is shown in figure 3.2 through the pressure field near the airfoil. The shock is clearly seen on the upper part. The two computations give similar results: they are indistinguishable by just looking at the pressure field. This is to be expected as the traditional method gives already satisfying results on such applications. The pressure coefficient around the airfoil is given in figure 3.3. Once again, the two curves are indistinguishable. Finally, the aerodynamic coefficients, the drag and lift coefficients, are given throughout the simulation in figure 3.4. This time, the JFNK method gives an improvement: the coefficients are much more stable, which means the convergence is better. The blue curve corresponding to the traditional method struggles to converge and seems to oscillate periodically. In comparison, the oscillation of the orange curve is dampened: this corresponds to convergence. The scale of the oscillation seen in figure 3.4 is not significant physically speaking, as it is negligible. It shows an issue with the traditional implicit method of the solver CEDRE: it fails to reach proper convergence but it oscillates around the solution. The new JFNK method however achieves proper convergence. Even if the figure shows an advantage of the new method, the result is almost meaningless as the oscillation of the blue curve is numerically insignificant. Looking at other data is necessary to properly conclude on the convergence of the two methods.

The residual is the best way to measure if a method can help the convergence of the solver. To define the residuals, let us use the partial differential equation (1.1). The residual is in fact the value of the function F from this equation. The name is quite adequate: for steady problems, the residual is what is left and still needs to be removed to find a steady solution. To decide on the convergence of a steady simulation, this residual norm is measured throughout the computational domain D. More precisely, this residual is analysed component by component. For example, the residual 1-norm associated with the ith component is:

F i 1 = D |F i | dv (3.1)
and the ∞-norm is:

F i ∞ = max D |F i | . (3.2)
In this fluid dynamics application, the residual norm is associated with the conservative variables: the density, the momentum components, the energy and the turbulent eddy viscosity ν t . With the Spalart-Allmaras compressible model, the conservative variable is in fact the product of the density and the Spalart-Allmaras variable ρν.

Figure 3.5 shows the residual norms from both computations. The gain from using the JFNK method appears clearly. The final residual norm is much smaller for each and every conservative component. This means that the Jacobian-Free Newton-Krylov method is able to reach a better convergence level than the traditional method. In other words, the better Jacobian matrix approximation leads to a better convergence than the older poor approximation.

In this first test case, we looked at the residual norms as the fields computed by both methods were almost identical. Moreover, our method was even able to reach better convergence levels. This validates our method on a first simple case, albeit with turbulence modelling and a shock.

COMPARISON BETWEEN MATRIX-FREE AND MATRIX APPROXIMATION

Turbulent transonic airfoil with boundary layer resolution

Definition of the test case

The previous case uses a mesh that is too coarse to resolve the boundary layer but instead uses a turbulent wall model. It is representative of CEDRE applications, but not of actual applications from the aerodynamic community. The RAE 2822 airfoil is often used for turbulence model validations, but computations in the field prefer to use a finer mesh to resolve the boundary layer. This is what we present next. The new mesh, shown in figure 3.6 is clearly different. It is a three-dimensional mesh, which is in fact a two-dimensional mesh that was extruded on 6 cells in the y direction. Periodicity is set between the two parallel sides of the domain. The mesh is made of 198 144 cells whereas the previous RAE 2822 mesh was only made of 36 541 cells. With these additional cells, this mesh is able to be fine on the wall boundary condition as seen in figure 3.6. This means that the use of the turbulent wall model is no longer necessary, unlike in the previous test case. Turbulence closure is still provided by the Spalart-Allmaras one equation turbulence model. The scale of the mesh is the same: the chord is still equal to 1m.

The flight conditions are also modified to match some existing experimental conditions: the AGARD-AR-138 experimental database. The Mach number is then 0.734 and the angle of attack is 2.54°, with a Reynolds number of 6.5 × 10 6 . This corresponds to Case 9 from the experimental conditions. The parameters are not exactly the same as they are corrected to account for the effects of the wind tunnel [START_REF] Hellstrom | Reynolds stress transport modelling of transonic flow around the RAE2822 airfoil[END_REF].

Analysis of the results

As before, the computation starts from an initial state constant in the domain. Then the computation is continued for a few more iterations with the traditional implicit Euler method. This is to reach a state close enough to the steady solution more quickly before analysing the convergence of our methods. After that, the computations are restarted, with the traditional method on one side and the matrix-free method on the other. The residuals are given in figures 3.7 and 3.8. Overall, it takes much more iterations to decrease residual norms for this test case compared to the previous one. It is due to the larger number of cells, and the presence of small flattened cells in the boundary layer. Figure 3.7 shows that the reduction in residual 1-norms with the standard method is not satisfying. Users usually expect the solver to reduce the order of magnitude of the residuals by 4 or 5, which is more than what our solver does. In contrast, the matrix-free method reduces residual norms much more efficiently and converges to a more precise solution. The quality of the solution computed by the matrix-free method is then better.

The ∞-norm of the residual is given in figure 3.8 for the same simulations. This time, the traditional method is not able at all to reduce residual norms. It means there are some cells in the computational domain for which the method is not able to reduce the local residual. The JFNK method in orange does not have this issue and converges to a better solution. The width of the curve is only due to data noise, as all lines use the same width.

The value of the turbulent variable ν in the first cell above the wing profile is provided for both computations in figure 3.9. The top views show the upper part of the airfoil, and the discontinuity due to the shock is visible, as expected. The bottom views show the lower part of the airfoil. The views on the right correspond to the JFNK computation, whereas the ones on the left to the computation with the traditional method. This traditional method using an approximated Jacobian matrix produces a solution of lower quality, as the pattern is not physical. In particular, there are cells for which ν is null in the computational domain, which should not happen with the standard Spalart-Allmaras model. Figure 3.8 shows that the ∞-norm of the residual on the conservative variable ρν normalised by the initial residual is near 10 0 , so figure 3.9 highlights cells for which the residual is higher than 10 -0.1 . Those cells correspond to the ones for which ν = 0. It shows a flaw in the standard time integration method, as it has non-physical features in the solution that prevent residual convergence. In contrast, it shows the quality of the JFNK method on this application.

Hypersonic reactive sphere

The second application we selected to compare the Jacobian-Free Newton-Krylov with the traditional one is the computation of the flow around a hypersonic solid sphere. Because of the high energy of the surrounding flow, the air molecules can separate and even form a plasma. The hypersonic reactive sphere is a well-known case, for both experimental [START_REF] Lobb | Experimental measurement of shock detachment distance on spheres fired in air at hypervelocities[END_REF] and numerical studies [START_REF] Dobrov | Numerical simulation of hypersonic flow with non-equilibrium chemical reactions around sphere[END_REF]. This is a simple yet representative test case of CEDRE applications. It then makes a lot of sense to analyse our new method on it.

Definition of the test case

The two-dimensional mesh is shown in figure 3.10. It is a regular mesh made of quadrangles, with refinement in the radial direction at the wall boundary and the expected shock location. The refinement at the wall boundary helps to compute more precisely the physical phenomena that happen at said boundary. The refinement at the shock is a requirement of the spatial discretisation methods. In order to get a clean and slim shock, the cells near the shock must have a high aspect ratio. The mesh takes this into account and is made according to CEDRE best practices. This expected shock location is obtained from a previous computation, made with a coarser mesh. This mesh is used in a two-dimensional axisymmetric computation to get the flow around a three-dimensional solid sphere. Figure 3.9: Spalart-Allmaras variable ν on the RAE 2822 wing profile for the traditional method (left) and the JFNK method (right). The top views show the profile from above, and the bottom views from below. The highlighted region in green corresponds to cells where the absolute value of the residual on ρν is higher than 10 -0.1 .

The solid sphere is modelled by an isothermal wall boundary condition. This is a representative choice as well, as the heat flux going through the wall is one of the main interests of such computation. Moreover, as it depends on the derivatives of the flow variables, it is usually harder to get a correct value for the heat transfer. A better convergence will lead to better derivatives, which will lead to better physical results for such case users. No turbulence model is used, as the flow is mostly laminar.

The main feature of this test case is that it simulates a hypersonic flow. At the left, the input boundary condition feeds air at Mach 15. This will induce a strong shock, meaning a strong discontinuity in the flow. It is indeed seen in figure 3.10 on the pressure and temperature. As this is a typical application of our solver, it is of importance to us that the new method behaves well with such flow features. When going through a strong shock, the temperature of the flow will increase a lot. The flow is made of air, or a mixture of 77% N 2 and 23% O 2 . At the high temperature they reach after the shock, the molecules can decompose into N, O, and NO. They can even get ionised. The choosen model uses 11 possible species: N 2 , O 2 , N, O and NO, the corresponding cations N 2 + , O 2 + , N + , O + and NO + and the electrons e -.

We argued previously that our solver is quite complex and that it makes it hard to do correct computations for non-experimented users. This hypersonic reactive sphere is a perfect example of this statement. Even if this case looks simplistic, it should not be taken lightly. When using unfit spatial discretisation methods, we ended up with a well-known problem of a such simulation called carbuncle [START_REF] Maccormack | Carbuncle computational fluid dynamics problem for blunt-body flows[END_REF]. This phenomenon does not come from physics but is solely due to numerical methods. It is shown in figure 3.11 from some earlier computations that were using a different mesh. Simply put, it creates a recirculation bubble downstream of the shock, near the stagnation point. This bubble tends to push the shock farther from the sphere, and modify greatly the heat flux at this location. The choice of the Riemann solver ended up being the key element to getting rid of this undesirable effect: there is no recirculation with the AUSM+ scheme. There are other pitfalls regarding this test case in CEDRE. For example, the user can choose whether to interpolate the mass fraction y i for the species at the cell interfaces for the MUSCL scheme or the mass concentration ρy i . Choosing the default option leads to nonconverging residuals no matter the time integration method. Indeed, this default option is recommended when running simulations with multiple phases, with significant variations in the density. Because CEDRE is made to solve a large variety of problems, it has a lot of methods that come with a lot of parameters, so even a simple simulation such as this one requires a lot of knowledge if a good convergence is required. This test case is a typical application of CEDRE, as opposed to the previous one which focuses more on the aerodynamic properties. Indeed CEDRE does not try to be a pure aerodynamic solver but a multiphysics one, equipped to solve high-energy problems. Reentry phenomena are therefore in the scope of our solver. Any improvements coming from our method on such applications would benefit a lot of CEDRE users on their applications.

COMPARISON BETWEEN MATRIX-FREE AND MATRIX APPROXIMATION

Analysis of the results

Once again, a first computation is done using the more robust traditional method. Indeed, the flow in the first cell against the sphere and the symmetry axis is initialised with a Mach number of 15 going straight into the wall. The matrix-free approximation struggles with such stiff initialisations, which is why the computation must start with another method, just until the shock has started to detach from the wall. Practically this amount to just a few iterations, compared to the number required to achieve convergence.

We first look at the different fields from the two computations. Once again they are similar and the residual norms are used to compare them, and this shows the simulation accounts for the features of interest. In figure 3.10, we see that the shock is present and that it falls as expected in the refined region. We see in the same figure the mass fraction of the various species downstream of the shock, which means the simulation is actively computing the chemical features of the flow, as expected.

To analyse the result from this case, one might once again look at the residuals. The first result corresponds to a computation without a reactive model. This simplification makes sense as the case to better understand the behaviour of the methods on an easier problem in terms of numerical complexity. It also means a different mesh is used, as not accounting for ionisation means the shock location is slightly different, but the two meshes are built the same way. The residual 1-norms are shown in figure 3.12. It shows that the new method gives a better convergence. The difference between the traditional method and the new one is that the traditional method uses the Jacobian matrix of the first-order spatial discretisation method, whereas the new one approximates the true Jacobian matrix using the approximation (1.53). This can indeed be verified. Because this simulation does not use complex models, the poor approximation used in the traditional model should use the Jacobian matrix of the first-order scheme. If the approximation (1.53) uses the first-order evaluation of f , then it should approximate this same Jacobian matrix. The evolutions of the residual norms from this method match the ones from the traditional method in figure 3.12. It first validates the development of the matrix-free approximation, and that the unbridle method should then use the actual second-order Jacobian matrix. It also confirms that the traditional method uses the Jacobian matrix of the first-order spatial discretisation method. Finally, it validates the choice of ε and the matrix-free approximation, as it is able to give the same result as when we use the Jacobian matrix.

With the new method, using the true function f with a second-order lead to much smaller residual norms. The only difference between this method and the previously existing one is that the latter uses the Jacobian matrix of the first-order spatial discretisation method, but the JFNK method takes into account the MUSCL reconstruction when using the Jacobian matrix. Because it uses a better Jacobian matrix in the linear problem, it gives a better solution to the nonlinear solver, which gives in turn a better time integration. It confirms that using a better Jacobian matrix helps the overall convergence when solving steady problems and validates the choice of the Jacobian-Free Newton-Krylov method.

When adding the reactive model, the conclusion is not in favour of the JFNK method anymore. The evolution of the residual norms is shown in figure 3.13. This time, the JFNK method does not longer find smaller residuals than the traditional method. In fact, when comparing figures 3.12 and 3.13, it appears that the traditional method converges better, enough so that it converges better than the JFNK method. This appears to happen because the mesh quality is better in this computation, compared to the mesh used without a reactive model. However, understanding the mesh quality is hard even on this simple test case for the reconstruction methods used in CEDRE. This is under current investigation by the responsible team. Still, the difference in the convergence of both methods is relatively small. A possible interpretation of this result is that when the mesh is not ideal for the spatial discretisation methods, using the JFNK method leads to a better convergence as it uses more precise Jacobian matrices for the linear system. When the mesh quality is better, the JFNK method is slightly less converged. Overall, it is still an improvement in the convergence on such test cases.

In a previous part, we suggested that the poor quality of the Jacobian matrix used in our solver was an obstacle to achieving good convergence. We then proposed an improvement: the Jacobian-Free Newton-Krylov method using the matrix-free approximation (1.53). We compared this new method with the traditional one on simple yet representative applications of our solver. This comparison showed that indeed, using a better Jacobian matrix leads to a lower residual norm. The new method can be useful when looking at quantities that require precise convergence, for instance, the derivatives of the flow variables. However, the main drawback of this new method

USING THE MATRIX-FREE FORMULATION WITH A NEW FLUID MODEL

63

is its computational cost. This cost is not inherent to the method but comes from our solver properties as was discussed previously. The recommended usage is then to start a computation using the inexpensive traditional method, and then achieve better convergence using the more precise yet more expensive new method.

Using the matrix-free formulation with a new fluid model

In the previous section, the new Jacobian-Free Newton-Krylov method was compared to the traditional method that uses the first-order Jacobian matrix. This first-order matrix is easily computed for the usual Navier-Stokes equations. But as CEDRE is under constant development, new models are added in order to handle more finely various multiphysics phenomena. For example, a new model was recently added to better handle multiphasic flows [START_REF] Cordesse | Contribution to the study of combustion instabilities in cryotechnic rocket engines : coupling diffuse interface models with kinetic-based moment methods for primary atomization simulations[END_REF]. Another model that we will investigate in the following was added to better account for thermodynamic disequilibriums. Adding a new model amount to writing the function F from equation (1.1), or equivalently G from equation (1.5). Usually, the development of this function is straightforward as it is given explicitly from the equations. However, the development of the Jacobian matrix is much more difficult. This is why the new models can not yet compute Jacobian matrices, and therefore can not use implicit methods as of today. The Jacobian-Free Newton-Krylov method is then a nice way to use implicit methods as it does not require the Jacobian matrix. Because the matrix-free approximation is written in a generic fashion, it can be used on any fluid model. In the following, we will use the JFNK method on a new fluid model that does not have an available Jacobian matrix and therefore can not use traditional implicit methods.

Multi-TEmperature model

The new model we will use is called the Multi-TEmperature model, or MTE. Its difference from the traditional Navier-Stokes model is that it allows for a thermodynamic non-equilibrium of some flow components. A particle has multiple degrees of freedom: translation, rotation and vibration. In more traditional models, it is assumed that all modes are at equilibrium, and they are grouped in what we call energy. One can define a time constant for each energy mode that corresponds to the number of collisions that are required to get to the equilibrium. A few collisions are required for translation modes, and about ten for rotation modes. This gives a small time constant of order 10 -9 s. For vibration modes, however, up to 20 000 collisions are required. The corresponding time constant is significantly larger, and there may be regions where there can be a disequilibrium between vibrational energy on one side and translation and rotation energy on the other. The energy of such components can no longer be described with a single temperature. As electrons are much lighter, they move more than the other heavier flow components. Their transitional energy and the transitional energy of other components can not be described with the same temperature. With the new Multi-TEmperature model, the flow components may be divided into three classes:

• the ones that always are at the equilibrium • the ones that may be at vibrational disequilibrium • the electrons that are handled separately from the other heavy component.

To account for the disequilibrium, the Navier-Stokes equations (1.3) are modified into the following:

                                                                         ∂ t (ρ s ) + ∇ • (ρ s u) = ∇ • ρD s ∇y s + ωs ∂ t (ρu) + ∇ • ρu ⊗ u + (p + p e ) Id = ∇ • µ ∇u + ∇u T - 2 3 µ (∇ • u) Id ∂ t (ρ m e v,m ) + ∇ • (ρ m e v,m u) = ∇ • λ v,m ∇T v,m + ρe v,m D m ∇y m + S v-t m + S v-v m + S v-e m + ωm
) u) = ∇ • λ eq ∇T + m λ m ∇T m + λ e ∇T e + µ ∇u + ∇u T u - 2 3 µ (∇ • u) u + ρ s h s D s ∇y s (3.3)
The first noticeable feature is that the equations of this physical model are more complex compared to the standard Navier-Stokes equations. Furthermore, there are multiple energies corresponding to the multiple temperatures: e e the energy of the electronic gas, e v,m the vibrational energy for each flow component that may be at disequilibrium, and E the total energy. It increases the number of conservative variables. The variables y correspond to the mass fraction of the different flow components. The source terms ω correspond to chemical production or decay.

The source terms S correspond to the energy transfers between the different energy modes. The details of those terms can be complex and are not discussed here as it is not the subject of this thesis, but can be found in [START_REF] Soubrié | Prise en compte de l'ionisation et du rayonnement dans la modélisation des écoulements de rentrée terrestre et martienne[END_REF]. A noticeable feature of this model is the presence of a nonconservative term in the conservation equation of the electrons' energy: u • ∇p e . This term is due to the effect of the electric field. This is a known issue of this model, as the Finite Volume method handles better conservative terms. There are two ways of handling this issue. The first one uses electronic entropy instead of energy as a conservative variable [START_REF] Coquel | A Roe-type linearization for the Euler equations for weakly ionized multi-component and multi-temperature gas[END_REF]. To do so, it makes some simplifications concerning electronic dissipation terms. However, this simplification can lead to inaccuracy in the physical result, and it appears that using the nonconservative formulation is necessary to get physically accurate results [START_REF] Soubrié | Prise en compte de l'ionisation et du rayonnement dans la modélisation des écoulements de rentrée terrestre et martienne[END_REF][START_REF] Kim | Modeling of electron energy phenomena in hypersonic flows[END_REF]. The other approach is the one used in our model, which handles the nonconservative term from equation (3.3). Not to go into too much detail, such models that account for thermodynamic disequilibrium are still under discussion as of today [START_REF] Blanco | Numerical modeling of hypersonic weakly ionized external flowfields with Poisson's equation[END_REF].

The Multi-TEmperature model allows for non-equilibrium between the vibrational mode of the flow components, the electronic energy and the total energy. Taking this non-equilibrium into account helps to compute more precisely physical quantities such as the electronic density, which can be useful for magnetohydrodynamics applications for example. It gives a better prediction of shock layer radiation with a better representation of the population of energy states which helps the simulation of radiation cooling or the computation of wall heat fluxes. This model is most interesting on reentry problems, as there is indeed some non-equilibrium downstream of the strong shocks that appear on such problems. Taking the non-equilibrium into account can change the flow downstream of the shock and on the solid wall, and give more precise temperatures and chemical compositions. Generally speaking, the Multi-TEmperature model gives more precise results than the standard Navier-Stokes model on fast flow with low relative density where there may be some thermodynamic non-equilibrium, such as reentry problems or problems with the strong expansion of a plume.

Hypersonic reactive sphere

The hypersonic sphere test case is a perfect fit for the Multi-TEmperature model. There is some thermodynamic non-equilibrium in the region downstream of the shock, because of its intensity. Unfortunately, the Jacobian matrix for the new MTE model is not currently available. This limits the MTE model users to explicit time integration methods, and therefore small time steps for stability reasons. As they deemed explicit time integration too slow, they wanted to use the newly implemented Jacobian-Free Newton-Krylov method for their simulations. We will show in the following another interest of the matrix-free approximation: it allows for implicit time integration methods despite having no Jacobian matrix. The goal here is to be able to get a steady solution to the problem as fast as possible for the user. Initially, the MTE model users were using a second-order Runge-Kutta method: the Midpoint method. We will compare the JFNK method to this reference method, and we will look in particular at how much time is required from the user to get to a physically satisfying result.

Full reactive model

For the first step of our analysis, we will use the most complete physical model. We use the same 11 flow components as before: N 2 , O 2 , N, O and NO, the corresponding cations N 2 + , O 2 + , N + , O + and NO + and the electrons e -. We of course use the Multi-TEmperature model, and we consider that N 2 and O 2 may be at vibrational disequilibrium. Other components may in fact also be at disequilibrium and we could take them into account but the computational cost increase is not worth it. As they are less present, and they quickly get to the equilibrium relatively to N 2 and O 2 , accounting for their disequilibrium would not be significant in the results [START_REF] Park | Thermochemical relaxation in shock tunnels[END_REF]. This amount to having four temperatures to describe the flow: T e for the electrons, T v,N2 and T v,O2 for the vibration modes of N 2 and O 2 , and T for the total energy.

As before, the Jacobian-Free Newton-Krylov method needs some help at the beginning. This is why the computation starts with some iterations of the Midpoint method and uses the firstorder Finite Volume method. Then, it continues with the JFNK method, still using the first-order Finite Volume method, to quickly get a good enough approximation of the expected solution. Finally, to get a finer result, the simulation switches the second-order spatial discretisation method. Local time-stepping based on the CFL number is used with the implicit simulation using the JFNK method. The GMRES method can not use preconditioners that require the matrix, so the only one available is the diagonal preconditioner based on the cell volumes. It is not ideal as it is extremely simple, and some physics-based preconditioner would be preferable [START_REF] Knoll | Jacobian-free Newton-Krylov methods: a survey of approaches and applications[END_REF] but none were available to us at the time. The spatial discretisation method used is the Multislope method, as it appears more robust than the k-exact method with the Multi-TEmperature model. The goal of this thesis is not the physical analysis of this test case so it will be left out, but according to the MTE model users, the results are in agreement with the literature. Figure 3.14 the various temperatures along the symmetry axis. This justifies the use of the Multi-Temperature model: the region downstream of the shock is indeed in a thermodynamic nonequilibrium state.

We will also look at the residuals for this analysis. But this time, using the number of iterations as the x-axis would be unfair to the explicit method. Indeed the cost of one iteration is a lot smaller for the explicit method than for the implicit JFNK one. What matters to us and a typical user is not the number of iterations but the time spent waiting for the result. It is called the elapsed real time or wall-clock time. As the computations run each time in the same parallel environment, with a fixed number of CPU cores, the elapsed time is proportional to the CPU time. The x-axis will then be the elapsed real time: the actual time it took to get to the current residual.

Figure 3.15 shows the residual norms for the two computations. The leftmost black curve corresponds to the initialisation: a few iterations of the explicit time integration method. From then the computation can continue with the Midpoint method on one side and the Jacobian-Free Newton-Krylov method on the other. We see that the residual norm obtained with the JFNK method is lower than its explicit counterpart. Indeed, the explicit method is limited to small time steps because of the shock and the stiff reactive model. The implicit method however is not. The fact that the implicit method can use larger time steps is not new, but with this specific model, it was impossible to use implicit methods until the matrix-free approximation was added. The JFNK method even ends up being cheaper in terms of CPU time. This means that for a user, it is faster to use the newly implemented Jacobian-Free Newton-Krylov method than to use explicit methods. Our method gives a better alternative to users working on problems with the Multi-TEmperature model.

Simplified reactive model

The equations concerning electrons in the Multi-TEmperature model bring a lot of issues. It is because the quantity of free electrons is usually relatively small when compared to other flow components. It can lead to numerical instabilities or cause other problems. As the equations hold everywhere in the numerical domain, they hold in particular in areas without free electrons, upstream of the shock in our hypersonic sphere for example. The question one must ask is then how to compute quantities such as the electronic temperature T e where there are no free electrons. We will not go into such implementation details as it goes beyond the scope of this work, but the point is that having electrons is quite troublesome for numerical methods. According to users that work on such reentry applications, it is not always useful to compute fluid ionisation. We explained why it might be crucial for some simulations, but the impact of ionisation on other quantities is almost negligible when the interest is not on the electrons. Sometimes, disregarding degrees of freedom corresponding to the electrons leads to simpler yet meaningful computations. The corresponding computation is drastically simplified as the contribution of the electrons is difficult to handle. The result however still gives interesting information, when one is not interested in the effects of ionisation. It is done by considering only the chemical components N 2 , O 2 , N, O and NO.

We decided to do the same comparison as before but this time with the simplified reactive model. The computation starts once again with a few iterations of the Midpoint method. Then, we compare the explicit method with the Jacobian-Free Newton-Krylov method. We show the corresponding residual norms in figure 3.16, while still using the elapsed real time as the xaxis. The result is similar to the previous one, except it is even more in favour of the JFNK method. The implicit method reaches much lower residual norms than the explicit one. With the full reactive model, the limiting factor for the implicit method came from the stiffness due to the electrons. Here, the term stiffness corresponds to both stiffness from the equation and difficulties arising from the implementation, that uses tools such as min/max functions that add non-differentiability. Without them, the convergence is significantly faster. However, the time step for the explicit method is under the same limitations with both reactive models. The gain in CPU time is then huge, in favour of the JFNK method.

Conclusion

This chapter compared the performances of the Jacobian-Free Newton-Krylov method with already existing methods of CEDRE on typical applications. When compared to the traditional implicit method on turbulent RANS computations, it showed that the JFNK achieves better convergence in the residual norms, when the convergence of the traditional method is unsatisfactory. The conclusion is similar when working with typical reentry applications. When the traditional method can converge well, however, the result is no longer in favour of the JFNK method, although the difference between the two methods is not significant. Finally, the matrix-free method was used with a newly implemented fluid model that does not give access to a Jacobian matrix yet. For this reason, using the JFNK method is the only way of using an implicit method. It was compared to the explicit Midpoint method that is currently used for computations with this fluid model. The results are once again in favour of the JFNK method, both in terms of convergence and speed. The superiority of the matrix-free method depends on the model's complexity level: the required time to reach the same convergence in residual norms is divided by 2 with the full model, and by 4 without ionisation.

Chapter 4

Introduction to exponential integration methods

Résumé du chapitre : Introduction aux méthodes d'intégration exponentielles

Le but de ce chapitre est de d'introduire un nouveau type de méthodes d'intégration temporelle.

Les méthodes classiquement utilisées dans la résolution des problèmes stationnaires sont souvent réutilisées après être adaptées dans la résolution de problèmes instationnaires intégrés avec de grands pas de temps pour s'affranchir des critères de stabilité de l'intégration explicite. Nous présentons dans ce chapitre le principe des méthodes d'intégration exponentielles, qui sortent de la dichotomie explicite/implicite. Nous décrivons ensuite plus en détail la famille des méthodes Rosenbrock exponentielles, qui seront utilisées par la suite. Nous réalisons également une brève analyse théorique de ces méthodes.

Comme leur nom l'indique, ces méthodes reposent sur l'exponentielle. Leur difficulté vient du fait qu'elles nécessitent de calculer des exponentielles de matrices, de très grande taille dans nos applications. Nous détaillons donc comment ce calcul est réalisé en utilisant des sous-espaces de Krylov.

Finalement, nous essayons une méthode exponentielle développée dans CEDRE sur un cas analytique simple. Nous montrons alors la pertinence de cette nouvelle méthode au vu de ses performances sur ce cas.

In this thesis, we were interested in finding solutions to steady problems. As we explained in the previous part, we use implicit time integration methods to efficiently get solutions to such problems. This is a pretty standard choice: most computational fluid dynamics solvers use implicit time integration methods to solve steady problems. The reason is that they can use larger time steps than their explicit counterparts. Because of this advantage, they are also often used when solving unsteady problems but with large time steps. Indeed, solving an unsteady problem with large time steps is similar to solving a steady problem. Most of what is done towards solving the steady problem can therefore be reused for unsteady problems with large time steps.

73

Even if implicit time integration methods are quite standard when solving problems with large time steps, there are other less conventional methods we could choose from. We could step out from the explicit implicit dichotomy, and decide to use IMplicit-EXplicit, or IMEX, methods. They split the function from the ordinary differential equation (1.5) into two parts: the stiff part that is integrated by an implicit method, because of its stiffness, and the other part that can just be integrated with an explicit method. The Additive Semi-Implicit Runge-Kutta methods, or ASIRK methods, are such methods [START_REF] Zhong | Additive semi-implicit Runge-Kutta methods for computing high-speed nonequilibrium reactive flows[END_REF]. They are already in use in computational fluid dynamics problems, such as fluid-structure interaction problems [START_REF] Huang | High-order, linearly stable, partitioned solvers for general multiphysics problems based on implicit-explicit Runge-Kutta schemes[END_REF]. Previous work already implemented ASIRK methods in our solver CEDRE for specific multiphysics applications. Some methods are even less common and correspond to a total paradigm shift: the parallel time integration methods [START_REF] Nievergelt | Parallel methods for integrating ordinary differential equations[END_REF][START_REF] Lions | Résolution d'EDP par un schéma en temps «pararéel[END_REF]. Just as we classically split the computational domain over processes and compute the spatial discretisation method in parallel, parallel time integration methods decompose the time integration interval into subintervals. They then solve the ordinary differential equation on each interval concurrently then ensure the continuity between the subintervals. They then iterate with Newton's method to find the solution over the whole time integration interval. As a result, they can approximate accurately the solution at a later time without knowing accurately the solution at a previous time. Despite being nontraditional, they were successfully used to solve fluidstructure interaction problems or Navier-Stokes problems [START_REF] Gander | Analysis of the Parareal time-parallel time-integration method[END_REF]. They were even more recently used to solve simple turbulent flow problems [START_REF] Lunet | Stratégies de parallélisation espace-temps pour la simulation numérique des écoulements turbulents[END_REF]. They can even be used in a more convoluted way, with for instance exponential methods [START_REF] Gander | PARAEXP: a parallel integrator for linear initial-value problems[END_REF]. Just as parallel time integration is inspired by parallel spatial discretisation, some other time integration methods are inspired by spectral discretisation methods. Time spectral methods, that were originally used for fluid dynamic time-periodic problems [START_REF] Gopinath | Time spectral method for periodic unsteady computations over two and three dimensional bodies[END_REF][START_REF] Gopinath | Application of the time spectral method to periodic unsteady vortex shedding[END_REF], are now used on non-periodic problems [START_REF] Ekici | Modeling periodic and non-periodic response of dynamical systems using an efficient Chebyshev-based time-spectral approach[END_REF]. But as both time parallel integration methods and time spectral methods are considered highly unconventional they require a major refactoring of the solver, which should be avoided for large-scale industrial solvers. In this part, however, we prefer to focus on another class of time integration methods.

Exponential integration methods

Let us focus on methods known as exponential methods. Despite being known for a long time [START_REF] Pope | An exponential method of numerical integration of ordinary differential equations[END_REF], they were not widely used in computational fluid dynamics, because of some difficulties that are explained later, but still interested some scientists [START_REF] Edwards | Krylov methods for the incompressible Navier-Stokes equations[END_REF]. They started to come back in the literature [START_REF] Hochbruck | Explicit exponential Runge-Kutta methods for semilinear parabolic problems[END_REF][START_REF] Hochbruck | Exponential integrators[END_REF] and are now being used on applications similar to ours [START_REF] Nie | Efficient semi-implicit schemes for stiff systems[END_REF][START_REF] Bhatt | Efficient Krylov-based exponential time differencing method in application to 3D advection-diffusion-reaction systems[END_REF]. They are still actively studied [START_REF] Narayanamurthi | Efficient implementation of partitioned stiff exponential Runge-Kutta methods[END_REF][START_REF] Caliari | BAMPHI: matrix-free and transpose-free action of linear combinations of ϕ-functions from exponential integrators[END_REF].

The ordinary differential equation

dy dt = f (y) (4.1)
is the same as the ordinary differential equation (1.5) from the previous part but with different notations. The main idea of exponential integration methods is to start from the ordinary differential equation (4.1) and split the function f into a linear part and a nonlinear part:

dy dt = Ly + N (y) . (4.2)
This decomposition is sometimes natural for some particular equations, but in the more general case it is always possible: it consists in choosing a linear part L and then setting the nonlinear part N (y) = f (y) -Ly. There are then an infinite number of decompositions, but we will see later that some are more interesting.

EXPONENTIAL INTEGRATION METHODS

75

To define the time integration step, we start from the current estimate of the solution y n that we assume to be exact: y n = y (t n ). We note ∆t = t n+1 -t n as we work here with a fixed n. We can integrate equation (4.1) using the variation of constants formula to get:

y (t n+1 ) = e ∆tL y n + tn+1 tn e (tn+1-t)L N (y (t)) dt. (4.
3)

The exponential integration methods then approximate the integral to compute the next value y n+1 . What defines the method is how it approximates this integral. For instances, methods that use the Taylor expansion of the nonlinear term are called exponential Taylor methods [START_REF] Koskela | Exponential Taylor methods: analysis and implementation[END_REF].

Here we focus on methods that use Runge-Kutta like quadratures. As we can see in equation ( 4.3), the linear part is treated exactly and the nonlinear one is approximated. If there is no nonlinear part, with N = 0, then the solution y n+1 is exact. If there is no linear part, with L = 0, then equation ( 4.3) transforms into

y (t n+1 ) = y n + tn+1 tn N (y (t)) dt (4.4)
and the exponential integration method behaves like a standard time integration method. This is the advantage of exponential integration methods: at best they are exact methods, and at worst, they are equivalent to traditional methods.

Similarly to classic methods, there are explicit [START_REF] Bhatt | Efficient Krylov-based exponential time differencing method in application to 3D advection-diffusion-reaction systems[END_REF] and implicit [START_REF] Nie | Efficient semi-implicit schemes for stiff systems[END_REF] exponential integration methods. It depends on whether it uses y n+1 or not to compute the integral from equation (4.3).

Before continuing, let us define some functions that are going to be convenient later on. Let the functions ϕ k be defined by:

     ϕ 0 : z → e z ∀k ∈ N * , ϕ k : z → 1 0 e (1-θ)z θ k-1 (k -1)! dθ. (4.5)
They could also be defined using the recurrence relation:

   ϕ 0 : z → e z ∀k ∈ N, ϕ k+1 : z → ϕ k (z) -ϕ k (0) z , ( 4.6)
or directly using their analytic formula:

∀k ∈ N, ϕ k : z → +∞ i=0 z i (i + k)! . (4.7)
This analytic formula ensures it is well-defined for squared matrices.

Exponential Rosenbrock-Euler method 4.2.1 Definition of the exponential Rosenbrock-Euler method

The exponential Euler method is the most basic one, and it can be used to understand exponential integration methods. Just as the Euler method assumes that N (y) is constant and equal to N (y n ) in equation (4.4), its exponential counterpart makes the same assumption but in equation (4.3). It then gives:

y n+1 = e ∆tL y n + ∆tϕ 1 (∆tL) N (y n ) (4.8)
using the ϕ 1 function defined above, and finally:

y n+1 = y n + ∆tϕ 1 (∆tL) f (y n ) . ( 4.9)
Note the difference with the corresponding standard Euler method, for which the ϕ 1 function is replaced with the identity function. If there is no linear part in the decomposition, with L = 0, then nothing is treated differently by the exponential method, and the standard explicit Euler method is recovered.

In this method, the linearisation L = f (y n ) is assumed. Because of this choice, the method is called the exponential Rosenbrock-Euler method. We will discuss the nomenclature below. This choice feels natural and minimises the error of the method as is shown in the following. This choice is frequent among exponential integration methods, despite some methods using a fixed linearisation L = f (y 0 ). The difference is that a new Jacobian matrix is required at each iteration of the time integration method, but it is something we are used to with our traditional implicit methods. It is worth noting that we might want to try the implicit equivalent of this method by taking N (y) = N (y n+1 ) in equation (4.3). However, because we took L = f (y n ), N (y) = 0 and after a linearisation the two variants are equivalent.

Analysis of the exponential Rosenbrock-Euler method

From the definition of the method and equation (4.3), we have that the error made after one step is:

y (t n+1 ) -y n+1 = tn+1 tn e (tn+1-t)L N (y (t)) dt -∆tϕ 1 (∆tL) N (y n ) = tn+1 tn e (tn+1-t)L (N (y (t)) -N (y n )) (4.10)
We take the Taylor series of the nonlinear part:

N (y) = +∞ i=0 N (i) (y n ) i! (y -y n ) i (4.11)
and in particular, using the fact that

y (t) = y n + y (t n ) (t -t n ) + O (t -t n ) 2 = y n + f (y n ) (t -t n ) + O (t -t n ) 2 (4.12)

EXPONENTIAL RUNGE-KUTTA AND ROSENBROCK METHODS

77

the partial series of order 2 is:

N (y (t)) = N (y n ) + N (y n ) f (y n ) (t -t n ) + O (t -t n ) 2 . (4.13)
The error of the method is then:

y (t n+1 ) -y n+1 = tn+1 tn e (tn+1-t)L N (y n ) f (y n ) (t -t n ) + O (t -t n ) 2 dt = ∆t 2 ϕ 2 (∆tL) N (y n ) f (y n ) + O ∆t 3 . (4.14)
The decomposition chosen earlier is now justified: if L = f (y n ) and therefore N (y n ) = 0 then y (t n+1 ) -y n+1 = O ∆t 3 and the method is of order 2. This is a noticeable property of this method: despite being a single-step method it has a second-order of accuracy. A more rigorous analysis can be found in [START_REF] Hochbruck | Exponential Rosenbrock-type methods[END_REF].

It is much harder to analyse the stability of this method, and more generally of exponential methods. Indeed, the stability analysis is done on the Dahlquist test equation (1.14). But as this equation is linear, exponential methods can solve them exactly no matter the time step size. We can still write the single step equation (1.19) y n+1 = g (∆tJ) y n with g (z) = e z , and get the corresponding stability region which is the left half complex plane. We could then conclude that the method is A-stable. The issue with exponential methods is that the standard stability analysis is no longer pertinent. We tried to do a better stability analysis for simple exponential methods such as the exponential Rosenbrock-Euler method, but we did not succeed. Some work in the literature does define some stability notions for exponential methods [START_REF] Du | Stability analysis and application of the exponential time differencing schemes[END_REF], but usually takes a linear N . It goes against the idea that all the linear part of f is in L and N is purely nonlinear. We did not find a stability analysis that was satisfying to us.

Exponential Runge-Kutta and Rosenbrock methods

When we introduced explicit methods, we used more convoluted approximations than the explicit Euler methods for the integral in equation (4.4). This gave the Runge-Kutta methods. The same can be done for the integral in equation (4.3) to get exponential Runge-Kutta methods. The difference with a standard Runge-Kutta method is that the quadrature coefficients are now functions of the matrix L:

             y n+1 = e ∆tL y n + ∆t k i=1 b i (∆tL) N (y n,i ) with y n,i = e ci∆tL y n + ∆t i-1 j=1 a ij (∆tL) N (y n,j ) . (4.15)
The Butcher tableau that we used to represent the Runge-Kutta methods is also used to represent the exponential Runge-Kutta methods. For instance, the Butcher tableau of the exponential Rosenbrock-Euler method, which is a special simple case of exponential Runge-Kutta methods, is shown in table 4.1.

Until now, we have not been exact when naming the methods we constructed here. From the literature, exponential Runge-Kutta methods are defined just as we did above, except they use They are used when the remainder r is small or at least bounded in terms of A. This is not the case in our applications, and in many others so that is why another type of exponential method was introduced: exponential Rosenbrock methods [START_REF] Hochbruck | Explicit integrators of Rosenbrocktype[END_REF]. They can be seen as a variation of the exponential Runge-Kutta methods where the linearisation is done at each time step. The issue is that a new Jacobian matrix is needed at each iteration, but once again we are used to this with implicit methods.

The nomenclature is a bit blurry to us, as the distinction between exponential Runge-Kutta and exponential Rosenbrock methods is not the same as the one between Runge-Kutta methods and Rosenbrock methods. First of all, the exponential Rosenbrock method is not implicit contrary to its standard counterpart. Furthermore, it would make sense to us to name an exponential integration method base on the underlying method used to approximate the integral from equation (4.3), and exponential Rosenbrock methods do not use Rosenbrock methods to do it. The only parallel we found was that the Rosenbrock methods do a linearisation at each inner stage, and the exponential Rosenbrock methods do a linearisation at each step of the time integration. It just seems that in the literature, the term Rosenbrock indicates that the linearisation is computed at each step instead of once at the beginning, with no apparent link with Rosenbrock methods.

Just as the Runge-Kutta method has order conditions for its quadrature coefficients, the exponential Rosenbrock method quadrature functions must follow some rules to ensure the expected order. To be consistent, it must have:

k i=1 b i = ϕ 1 .
(4.17)

It is reasonable to want multiple-stage methods to have a higher order than the single-stage exponential Euler method, and therefore it must also have

i-1 j=1 a ij = z → c i ϕ 1 (c i z) , 1 ≤ i ≤ k (4.18)
to achieve a second order. Using these two conditions, exponential Rosenbrock methods can be written in the form:

             y n+1 = y n + ∆tϕ 1 (∆tL) f (y n ) + ∆t k i=2 b i (∆tL) D n,i with y n,i = y n + c i ∆tϕ 1 (c i ∆tL) f (y n ) + ∆t i-1 j=2 a ij (∆tL) D n,j (4.19)
using D n,i = N (y n,i ) -N (y n ). This form shows that exponential Rosenbrock methods are variations of the exponential Euler method. The defects D n,i are small in size, so less computational effort can be spent when computing their contribution.

0

ϕ 1 0 1 ϕ 1 ϕ 1 -2ϕ 3 2ϕ 3 0 3 4 3 4 ϕ 1 3 4 • ϕ 1 -32 9 ϕ 3 32 9 ϕ 3 Exponential Rosenbrock-Euler ExpRB32 ExpRB42
Table 4.1: Butcher tableau for some exponential Rosenbrock integration methods Among exponential Rosenbrock methods, we looked at the ExpRB32 method from [START_REF] Hochbruck | Exponential Rosenbrock-type methods[END_REF] and ExpRB42 methods from [START_REF] Luan | Fourth-order two-stage explicit exponential integrators for time-dependent PDEs[END_REF]. They were developed as adaptive methods, but we will use them without the error estimate as standard methods. They are both 2-stage exponential Rosenbrock methods and achieve respectively a third and fourth order of accuracy. This shows the quality of exponential Rosenbrock methods. Their Butcher tableaux are in table 4.1. For both methods, the coefficient that is the hardest to get is b 1 , however equation (4. [START_REF] Le Touze | Multislope MUSCL method for general unstructured meshes[END_REF]) tells us that it will not be computed. We decided to work with those two methods, along with the Rosenbrock-Euler method, as they are simple among exponential methods and are supposedly well suited for our problems. They give us methods with orders of accuracy from 2 to 4. But what is interesting to us is that on linear problems those methods are exact, or with an infinite order of accuracy. As our problems are not purely linear, it means that the linear part will be solved exactly while the nonlinear part will be solved with a still good order of accuracy.

Evaluating matrix functions

At first glance, exponential methods may seem too good to be true. Indeed, they give explicit methods with few stages but high order of accuracy, and they have the formidable quality that they solve exactly the linear part of the ordinary differential equation. In other words, the exponential Euler method should be able to get the solution of Poisson's equation after any time with a single step. Indeed, Poisson's equation gives often a linear ordinary differential equation after the spatial discretisation, with the Finite Volume method on a regular mesh for example. The exponential Rosenbrock-Euler method can solve this linear differential equation exactly, for any time step, as large as it may be.

The reason why exponential methods are not that good is the reason why they were not used for a long time after being introduced. It is because evaluating matrix exponential, or more generally evaluating the ϕ-functions on matrices, is not trivial. Indeed, the power series definition of the matrix exponential is not well suited for numerical evaluations. In particular, in our context, computing the powers of the matrix is completely out of the question. The dimension of our matrices can get large, and the powers of the matrices would lose the sparsity. In fact, any approximation working directly with the matrix would not be satisfying. The question of how to evaluate ϕ-functions in the context of exponential methods is a subject of active research, as in this extremely recent reference [START_REF] Caliari | BAMPHI: matrix-free and transpose-free action of linear combinations of ϕ-functions from exponential integrators[END_REF]. We will however use more standard well-established algorithms.

The other issue with exponential methods is the question of the Jacobian matrix. Exponential Runge-Kutta reuse the same matrix throughout the computation but it is not a good idea in our applications since the nonlinear part can change. Exponential Rosenbrock methods were

CHAPTER 4. INTRODUCTION TO EXPONENTIAL INTEGRATION METHODS

introduced for this reason but they require the computation of a new Jacobian matrix at each iteration.

Evaluating matrix functions with Krylov subspace methods

The solution for both issues is the same: Krylov subspace methods. The Cayley-Hamilton theorem ensure that the exponential of a matrix is a polynomial of degree less or equal to N where N is the dimension of the matrix. Instead of looking for the matrix exponential as an infinite series, we can look for an N -order polynomial. Moreover, exponential methods do not need to evaluate the exponential or the ϕ-functions of the matrix, but their effect when applied to a vector. For instance, the exponential Rosenbrock-Euler method from equation (4.9) does not need to compute ϕ 1 (∆tL) but ϕ 1 (∆tL) f (y n ). It means we do not need to evaluate a given function h on the matrix A, but to compute the effect of h (A) on the vector b. For exponential methods, h is a linear combination of ϕ-functions, and is therefore analytic. The idea, originally developed to approximate eigenvalues of a large sparse matrix, is to use a Krylov subspace method to approximate h (A) b [START_REF] Saad | Analysis of some Krylov subspace approximations to the matrix exponential operator[END_REF][START_REF] Sidje | Expokit: a software package for computing matrix exponentials[END_REF]. In fact, it is similar to using a Krylov subspace method to solve the linear system Ax = b: it means taking h = z → z -1 . With the Arnoldi iteration, the Krylov subspace method produces at the m-th iteration the relation:

AV m = V m+1 Hm = V m H m + h m+1,m v m+1 e T m . (4.20)
The columns of V m , v 1 , . . . , v m , form an orthonormal basis of the Krylov subspace:

K m (A, b) = Vect b, Ab, . . . , A m-1 b (4.21) with v 1 = b/ b 2 .
The matrix Hm is a Hessenberg matrix and H m is the square matrix equal to Hm without its last line:

Hm =        h 1,1 h 1,2 • • • h 1,m h 2,1 h 2,2 • • • h 2,m . . . . . . . . . . . . 0 • • • h m,m-1 h m,m 0 • • • 0 h m+1,m        , H m =      h 1,1 h 1,2 • • • h 1,m h 2,1 h 2,2 • • • h 2,m . . . . . . . . . . . . 0 • • • h m,m-1 h m,m      (4.22)
The vector v m+1 is the last column vector of V m+1 , added to the column vectors of V m to get an orthonormal basis of K m+1 (A, b), and e m is the m-th canonical vector (0, . . . , 0, 1) T . Since H m represent the compression of A on K m (A, b) with respect to basis V m and b = b 2 V m e 1 with e 1 = (1, 0, . . . , 0), we make the approximation [START_REF] Eiermann | A restarted Krylov subspace method for the evaluation of matrix functions[END_REF]:

h (A) b ≈ b 2 V m h (H m ) e 1 . ( 4.23)
With this approximation, we can evaluate the function h by applying it to a smaller matrix. We reduced the original function evaluation to the much smaller Krylov subspace. Since the method does not need to know the matrix but only to be able to apply it as a linear operator, we can use the matrix-free approximation (1.53) to approximate the Jacobian matrix L = f (y n ). This way, having to compute a new Jacobian matrix at each iteration of the exponential Rosenbrock method is not an issue. Just as Krylov subspace methods for solving linear problems, this method can be restarted to limit the size of the Krylov subspace and therefore control the memory consumption of the algorithm. It is even possible to define an error estimate in the approximation to dynamically control the number of iterations needed to reach a given tolerance, albeit it is more convoluted than the residual used for solving linear systems with Krylov subspace methods [START_REF] Eiermann | A restarted Krylov subspace method for the evaluation of matrix functions[END_REF].

EVALUATING MATRIX FUNCTIONS 81

Note: Since the matrix used to evaluate the quadrature functions of the exponential Rosenbrock methods is always the same, L, only one Arnoldi iteration is required for one step of the method. This could lead to a great cost reduction, however, we did not try it in our developments. Indeed, our work on exponential integration methods was more to demonstrate their feasibility in our applications, but it would be interesting to use this idea for future performance enhancements.

Evaluating matrix exponentials

With the Krylov subspace method just described, the next step is the computation of the exponential of the smaller Hessenberg matrix H m . Attention is now paid to how to compute the exponential of relatively small dense matrices. Many methods use special properties of the matrices but our matrices do not share the same properties. Instead of using a truncated Taylor series, using the Padé approximant seems to be the best option [START_REF] Higham | Computing matrix functions[END_REF]. The exponential function is approximated by the [k/m] Padé approximant:

e z ≈ r k,m (z) = p k,m (z) q k,m (z) with            p k,m = k i=0 (k + m -i)!k! (k + m)! (k -i)! z i i! q k,m = m i=0 (k + m -i)!m! (k + m)! (m -i)! (-z) i i! . ( 4.24)
It is best to use diagonal approximants, with k = m, as they are more accurate [START_REF] Higham | Computing matrix functions[END_REF].

The final method used to compute the exponential of a dense matrix is the Scaling and Squaring method [START_REF] Higham | The Scaling and Squaring method for the matrix exponential revisited[END_REF]. It computes an optimal parameter s and takes:

e z ≈ r m,m (z/2 s ) 2 s . ( 4.25)
The goal is to choose s such as the error on the computation of the exponential is bounded, with minimal cost. It then decides on an optimal order m for the Padé approximant. The matrix is scaled by 2 s , the Padé approximant is used to compute the exponential, and the result is then squared back to the original scale. It is what gives its name to the method.

Evaluating the ϕ-functions

The exponential time integration methods do not only use matrix exponential, the ϕ 0 function, but may use any ϕ k function. There is a nice relation that helps compute those functions. We first define the augmented squared Hessenberg matrix:

Hm+p =          H m c 0 • • • 0 0 1 . . . . . . 0 . . . 0 . . . 1 0 0          (4.26)
with a vector c of dimension m. Then, the following relation holds [START_REF] Sidje | Expokit: a software package for computing matrix exponentials[END_REF]:

exp τ Hm+p =             exp (τ H m ) τ ϕ 1 (τ H m ) c τ 2 ϕ 2 (τ H m ) c • • • τ p ϕ p (τ H m ) c 1 τ 1! . . . τ p-1 (p-1)! 1 . . . . . . . . . τ 1! 0 1             . ( 4.27)
The additional cost of using this augmented matrix is negligible in comparison with the cost of the rest of the methods. Furthermore, the methods we introduced before only need to augment with p ≤ 3. By adjusting the scalar coefficient τ and reading the result in the appropriate column, we can finally compute the quadrature coefficients of our exponential Rosenbrock methods.

We explained how to compute one step of an exponential integration method. From what we have seen, the computational cost of one step of an exponential method is roughly the same as the cost of one step of an equivalent implicit method. For example, the computational cost of the exponential Rosenbrock-Euler method is the cost of one Jacobian matrix computation, one Arnoldi decomposition and then one function evaluation on the Hessenberg matrix. The implicit Euler method, which is also a one-stage method, also computes a Jacobian matrix and does an Arnoldi decomposition. It then inverts the Hessenberg matrix with a QR decomposition which is equivalent to applying a function to the Hessenberg matrix. This is why the computational costs of the implicit Euler method and the exponential Rosenbrock-Euler method are similar.

When we wrote the general formula for exponential Rosenbrock methods in equation (4.19), we used the defects D n,i . The point of this formulation is to write the exponential Rosenbrock methods as corrections of the Rosenbrock-Euler method. Because the defects are usually small in size, we can spend less computational power to compute their contribution to the final step. Practically, it means we can use smaller Krylov subspaces to account for the defects. It means shorter Arnoldi iterations, which means smaller Hessenberg matrices, and overall a much smaller computational cost. We explained why the one-stage exponential method cost is equivalent to the one-stage implicit method cost. Because we use the formulation (4.19) and the fact that the defects are small, adding more stages to the exponential method will be less expensive than adding stages to the implicit method.

First developments in CEDRE

We first tried exponential methods in the solver CEDRE. The development was fast as it reuses a lot of existing parts from the implicit time integration: principally the Krylov subspace methods tools. The matrix used in the decomposition (4.2) can either be the traditional loworder Jacobian matrix or can use the matrix-free approximation (1.53). It shows that this work about exponential integration methods is compatible with our work about the Jacobian matrix, presented before. The Arnoldi decomposition (4.20) is computed by the same routines as the ones used for the GMRES algorithm. We only had to develop the Scaling and Squaring algorithm, with the computation of the Padé approximant to be able to use exponential Rosenbrock methods. Because it only required this development we decided to deviate from the original framework of this thesis and try this new type of time integration method. Now that everything was accessible, we started with the exponential Rosenbrock-Euler method, as it is the basic one.

To test the method, we decided to use a simple academic test case, adapted from the one proposed at the International Workshop on High Order CFD Methods (HiOCFD). It is a twodimensional inviscid isentropic vortex that is convected in a periodic box. This vortex is a known exact solution of the unsteady Euler equations and is generally considered a good test case to analyse the accuracy of both spatial discretisation and time integration methods. It is a non-dimensionalised computation, and the mesh is regular over a square of length L = 20, with 100 2 cells. The vortex is convected in a flow at a Mach number of Ma = 0.5, with the dimensioning parameters P ∞ = 1 and T ∞ = 1. The undisturbed velocity is then U ∞ e x with U ∞ = Ma γr gas T ∞ and the vortex is defined using the polar coordinates system (e r , e θ ):

                       u = U ∞ e x + βU ∞ r R 0 e -r 2 /2R 2 0 e θ T = T ∞ -β 2 U 2 ∞ γ -1 2γr gas e -r 2 /R 2 0 P = P ∞ 1 - β 2 U 2 ∞ T ∞ γ -1 2γr gas e -r 2 /R 2 0 γ/(γ-1)
.

(4.28)

The heat capacity ratio γ is equal to 1.4, and the specific gas constant is r gas = 1. The vortex is defined by its characteristic radius R 0 = 1 and its intensity β = 0.2. This initial condition is shown in figure 4.1. As we solve inviscid Navier-Stokes equations, or Euler equations, the vortex moves in the e x direction undisturbed, and as the domain is periodic the solution after one period ∆T = L/U ∞ = 33.8 is the same as the initial solution.

We compare the exponential Rosenbrock-Euler method to the standard implicit Euler method and the classic fourth-order explicit Runge-Kutta method. Both the implicit and the exponential methods need to compute Jacobian matrices. To reduce the number of variables in this comparison, both will use the traditional Jacobian matrices which is the low-order approximation. We then end up with two equivalent methods in terms of computational cost. We simulate the vortex advection throughout one period ∆T . The solution for each method is shown in figure 4.1. We first look at the bottom center and right simulations. We see that the explicit Runge-Kutta method preserves well enough the vortex compared to the implicit method. It is because the RK4 method is a fourth-order method whereas the implicit method is a first-order method. However, the simulations are made at a CFL number of 0.9. With a higher CFL number of 1.5, the explicit Figure 4.1: Initial solution (top left) and solution after one period for the exponential Rosenbrock-Euler method (top right), the implicit Euler method (bottom left and center) and the explicit RK4 method (bottom right). The first two simulations correspond to a CFL number of 1.5 whereas the last two to a CFL number of 0.9.

Chapter 5

Analysis of exponential integration methods in JAGUAR

Résumé du chapitre : Analyse des méthodes d'intégration exponentielles dans JAGUAR Le but de ce chapitre est de comparer les méthodes exponentielles aux méthodes d'intégration classiques de JAGUAR pour valider leur intérêt.

Nous avons montré dans le chapitre précédent que CEDRE n'était pas le contexte idéal afin d'analyser des méthodes d'intégration temporelles très précises comme celles qui nous intéressent. Nous utilisons donc le solveur JAGUAR, basé sur la méthode des Differences Spectrales, que nous présentons au début de ce chapitre.

Nous analysons ensuite les méthodes exponentielles ajoutées à JAGUAR grâce à l'utilisation de la librairie SLEPc sur un cas simple, afin de constater l'ordre de ces méthodes à travers l'erreur globale en fin de simulation. Sur un autre cas académique, nous comparons ensuite notre ensemble de méthodes en regardant le temps d'exécution qu'il leur est nécessaire afin de réaliser une même simulation. Pour finir, nous utilisons les méthodes exponentielles sur un cas plus significatif car proche d'une application réelle du solveur: une aube de turbine.

We want to analyse time exponential integration methods in a framework that gives us access to high-order spatial discretisation methods. It is possible to use high-order Finite Volume methods, but this means using larger stencils which hurts parallelism. In CEDRE, users often stop at second-order methods, so we need to use another solver. As the spatial discretisation method does not play a direct role in the time integration and the work done in this thesis, we can step out from the Finite Volume framework. Furthermore, using a less complex solver will also help us try and develop new methods more easily than we already did with CEDRE. This is why we decided to accomplish our analysis of exponential time integration methods with the solver JAGUAR.

CHAPTER 5. ANALYSIS OF EXPONENTIAL METHODS IN JAGUAR

JAGUAR: a Spectral Difference solver

JAGUAR means proJect of an Aerodynamic solver using General Unstructured grids And high ordeR schemes. It is a reactive Navier-Stokes solver originally developed at CERFACS and has been jointly owned by ONERA and CERFACS for four years. It is made for unstructured grids and uses a Large Eddy Simulation model to solve unsteady turbulence effects. It means that the large turbulence scales are computed explicitly while the smaller ones are modelled. Its particularity is that it uses a spectral method as a spatial discretisation scheme: the Spectral Difference method. The Spectral Difference method follows some principles of the Discontinuous Galerkin formulation, in the sense that both are families of schemes of arbitrary order using polynomial approximations, but their mathematical foundations differ.

The Spectral Difference method

With the Spectral Difference method, the solution is represented by a polynomial of degree p inside each cell. It means that in the partial differential equation

∂u ∂t + ∇ • F (u) = 0, (5.1)
u is a p-degree polynomial of the coordinate variables, where the polynomial coefficients are functions of time. Then F (u) has to be a p+1-order polynomial of the coordinate variables in order to have a p-order polynomial of the flux divergence. The key to the Spectral Difference method is how to compute a p + 1-order F (u) from a p-order u. This method uses key elements that were first mentioned by [START_REF] Kopriva | A conservative staggered-grid Chebyshev multidomain method for compressible flows. II. a semi-structured method[END_REF], and were later developed by [START_REF] Liu | Spectral difference method for unstructured grids I: basic formulation[END_REF].

Because the solution is represented in each cell by a polynomial, there is no reason for the solution to be continuous throughout cell interfaces. Indeed, the Spectral Difference method uses a piecewise representation of the solution. For example, let us consider the two-dimensional regular mesh, made of 2 × 2 cells over [0, 2]

2 . The exact analytical function u over the mesh in the left part of figure 5.1:

u : [0, 2] 2 → R (x, y) → cos (5x) tanh (5 (1 -y)) (5.2)
cannot be represented by polynomial exactly and the right picture of figure 5.1 represents the solution associated with the Spectral Difference method. The colour mapping corresponds to the output of the function but is not given as this function is of no particular interest: it is a rough example only. This function is interpolated with a second-order method using Lagrange polynomials over each cell. The result is shown in the right part of figure 5.1, where each colour corresponds to one cell. There are discontinuities at the cell interfaces, which shows a specificity of the Spectral Difference methods. If the function on the left part of figure 5.1 was used as an initial condition for a simulation, the Spectral Difference method would in fact use the polynomial interpolation on the right part. However, to respect the underlying conservation property of equation (5.1), it is necessary for F to be continuous over the whole computational domain. The Spectral Difference method will then make sure that the polynomial representations of the flux density are continuous throughout cell interfaces. To better understand how this method works, let us take a one-dimensional cell: the segment [0, 1]. Because the following is done at a fixed time, the dependency on the time is dropped, but the solution and the coefficients are in reality functions of the time, not scalars. The solution u inside this segment is then u (x) = p i=0 a i x i . Using Lagrange interpolation polynomials, it is equivalent to use the set of the p + 1 coefficients a i or the set of the p + 1 values u (x i ) computed at the distinct points x i ∈ [0, 1] called solution points. The solution as a p-order polynomial is represented by either one of those two sets. As the flux density F must be a p + 1-order polynomial, it can also be represented by its values in the p + 2 distinct flux points. To ensure that F is continuous at the segment end points, we take 0 and 1 as flux points. The choice of the rest of the flux points will be discussed later, but let us say for now that they are staggered with the solution points: each flux point, apart from the segment end points, is between two solution points and vice versa. The solution points are marked by red dots. Parts of the solution from the left and right neighbouring cells can be seen. As the figure shows, the piecewise local solutions inside each cell may be discontinuous at the mesh interface. This is the starting point of the Spectral Difference method and will be called the 0-th step.

1. In the first step, the method computes the value of the solution in the p + 2 flux points, marked by blue dots. As the polynomial representing the solution is known, this step just consists in evaluating it at flux points. 3. Because there are different values of the flux density at the segment end points, the method would not preserve the conservative property of the partial differential equation. This is why in the third step, the method computes a unique interface flux density for both cells at each segment end point. The problem is to find the interface flux at the discontinuity of a piecewise solution. In other words, this is a Riemann problem. Once again, this is solved with a Riemann solver, exact or approximate, to get in the end a single value for the left and right parts of the discontinuity.

4. An interpolation from the value of the flux density at the p + 2 flux points using Lagrange polynomials enables to get a p+1-order F as expected in the fourth step, represented with a blue line in the figure. Therefore, one can end up with a continuous representation of F , differentiable everywhere except at cell interfaces.

5. In the last step, the representation of F is differentiated to get the flux divergence, represented by a green line in the figure. Finally, a time integration method can use a p-order representation of ∇ • F to compute the solution at the next time step.

To work with any segment, not only [0, 1], an isoparametric transformation is introduced to get back to this unity segment. It is the same when working with multiple dimensions, where the isoparametric transformation gets the cell back to the tensor powers of this segment. The placement of the solution points does not seem to matter much, but the placement of the flux points does [START_REF] Van Den Abeele | On the stability and accuracy of the spectral difference method[END_REF]. The p + 1 solution points we will use are defined in the [0, 1] segment as the Chebyshev roots: As a side note, the Spectral Difference method with p = 0 corresponds roughly to the firstorder Finite Volume method. Indeed, the solution is assumed constant in the cell, represented by the value at its barycenter. The flux balance is made at the cell interfaces with a Riemann solver. This corresponds to the placement of the solution and flux points in figure 5.3 when p = 0. For any polynomial degree p of the solution, the proposed Spectral Difference method is naturally of order p + 1. Indeed, the method can represent exactly any polynomial of degree p and, as a consequence, the error term is of order p + 1. In JAGUAR, the polynomial order p ranges from 2 to 10, but 2 ≤ p ≤ 6 for practical applications.

x i = 1 2 1 -cos 2i + 1 2p + 2 π , 0 ≤ i ≤ p. ( 5
The demonstration of accuracy of the Spectral Difference method is proposed in a reference paper published recently [START_REF] Vanharen | Revisiting the spectral analysis for high-order spectral discontinuous methods[END_REF]. JAGUAR was optimised to be efficient for parallel computations [START_REF] Cassagne | Implémentation multi GPU de la méthode spectral differences pour un code de CFD[END_REF][START_REF] Cassagne | High-order method for a new generation of large eddy simulation[END_REF][START_REF] Marait | Performances du code de mécanique des fluides JAGUAR sur CPU et accélérateur[END_REF]. To perform Large Eddy Simulation, it is of paramount importance to implement an unsteady characteristic boundary condition: Navier-Stokes Characteristic Boundary Condition, or NSCBC. The initial formulation [START_REF] Fiévet | Strong compact formalism for characteristic boundary conditions with discontinuous spectral methods[END_REF] was extended to deal with acoustic conditions and liner optimisation [START_REF] Cardesa | Optimizing an acoustic liner by automatic differentiation of a compressible flow solver[END_REF]. Recently, the solver was extended to handle h -p adaptation [START_REF] Hartmann | Space adaptive methods/meshing[END_REF] and the schemes were extended to triangles and tetrahedrons [START_REF] Veilleux | Stable spectral difference approach using Raviart-Thomas elements for 3D computations on tetrahedral grids[END_REF][START_REF] Veilleux | A stable spectral difference approach for computations with triangular and hybrid grids up to the 6th order of accuracy[END_REF]. Another improvement concerns the treatment of combustion, with a PhD thesis funded by Safran [START_REF] Marchal | Extension of the spectral difference method to combustion[END_REF]. Finally, the last extension concerns the treatment of shocks [START_REF] Deniau | Entropy stable spectral difference scheme[END_REF]. This review shows the strong involvement of researchers in the solver.

Exponential integration methods in JAGUAR

Using a high-order spatial discretisation method and a refined mesh is the best-suited situation to test high-order time integration methods: the resulting error will come mostly from the time integration and not the spatial discretisation. The development of an exponential time integration method was inexpensive within CEDRE as the method reuses lots of already existing parts. For JAGUAR which only has explicit methods, the Arnoldi iteration and the exponential computation functions must be implemented to use exponential methods. They could be developed as internal stand-alone routines to produce a finely tuned method for JAGUAR. However, since the goal is to analyse exponential methods scientifically, it was decided to rely on the SLEPc library [START_REF] Hernandez | SLEPc: a scalable and flexible toolkit for the solution of eigenvalue problems[END_REF]. The Scalable Library for Eigenvalue Problem Computation, or SLEPc, is an extension of the software library PETSc [START_REF] Balay | PETSc Web page[END_REF][START_REF] Balay | PETSc/TAO users manual[END_REF][START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF]. Instead of rewriting the needed algorithms, it was chosen to use their SLEPc implementation. In particular, SLEPc handles what it calls Matrix Function objects or MFN: "Given a matrix A and a vector b, the call MFNSolve(mfn, b, x) computes x = f (A) b, where f is a function such as the exponential."

This is precisely what is needed to develop exponential methods: as it handles the previously introduced ϕ-functions with its MFN objects, SLEPc is then an obvious choice of a library for implementing exponential integration methods. Furthermore, because it relies on PETSc, it is efficiently scalable to fit our multiprocessing needs.

Exponential methods are based on a decomposition of the ordinary differential equation such as equation (4.2). However, JAGUAR uses only explicit time integration methods, so no Jacobian matrix is available. Computing analytically the Jacobian matrix would prove challenging as any ingredient should be differentiated, such as the Spectral Difference scheme, the Riemann solvers, the diffusion scheme, etc. It is not insurmountable but would amount to more work than what could be afforded during this thesis. It was decided instead to reuse the work from the previous part: the Jacobian matrix will not be formed, but its effect will be computed by a finite difference approximation. Another reason to work with the SLEPc library is that using a finite difference approximation is easy with it. More precisely, it is PETSc that handles this approximation. To sum up, to compute ϕ k (L) b that is required for the exponential methods with L = f (y), L is created with the appropriate PETSc data type. As it is a Jacobian matrix of a function, PETSc only needs to know the function to compute matrix-vector products. Then, after setting the MFN object of SLEPc, the library can compute the desired result. The Arnoldi iteration, the Scaling and Squaring algorithm and the exponential evaluation are all handled internally by the library. Overall, this procedure is extremely simple from our perspective and requires only a few lines of code to implement. This highlights the relevance of using the SLEPc library to easily test procedures and acquire the associated knowledge.

Other people that also work with JAGUAR did develop implicit time integration methods using PETSc. To continue in this direction, it was decided to add a time integration method based on the generic PETSc Time Stepper [START_REF] Abhyankar | PETSc/TS: a modern scalable ODE/DAE solver library[END_REF]. This way, a user may use most methods from PETSc in JAGUAR with no additional developments. It includes any explicit and implicit time integration methods, but also nonlinear solvers such as Newton's methods with various line search algorithms, linear solvers with various Krylov subspace methods and finally various preconditioning. However, this generic PETSc Time Stepper method was added on the side of this thesis, and therefore will not be discussed here, along with the previous work on implicit methods.

In the end, the high-order spatial discretisation method called the Spectral Difference method allows us to analyse and compare exponential methods that are available through the SLEPc library. In the following, attention is mainly focused on the three methods that were presented earlier in table 4.1: the Rosenbrock-Euler, ExpRB32 and ExpRB42 methods.

Analysis of exponential time integration methods

There are in JAGUAR all the tools necessary to analyse exponential time integration methods and several test cases must be chosen. During the analysis, the exponential methods and explicit Runge-Kutta methods will be compared, as the latest ones are the only available methods in JAGUAR and represent the state-of-the-art. Those methods are:

• RK2, the Midpoint method • RK4, the classical four-stage fourth-order method • RKo6s, a low-storage, low-dissipation and low-dispersion six-stage second-order method from [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for flow and noise computations[END_REF] dedicated to unsteady computations

• TVDRK(3, 3), a three-stage third-order TVD method from [START_REF] Shu | Efficient implementation of essentially non-oscillatory shockcapturing schemes[END_REF][START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF] • SSPRK [START_REF] Huang | High-order partitioned spectral deferred correction solvers for multiphysics problems[END_REF][START_REF] Coulette | Implicit time schemes for compressible fluid models based on relaxation methods[END_REF], an optimal SSP five-stage fourth-order method from [START_REF] Spiteri | A new class of optimal high-order strong-stability-preserving time discretization methods[END_REF].

Low-storage methods are Runge-Kutta methods for which the matrix A in their Butcher tableau is lower diagonal, or in other words all a ij are null except when j = i -1. Also, all b i are null except the last one that is equal to 1. It means that the final stage and all intermediate ones use only the previous stage. Therefore, there is no need to keep track of intermediate stages in memory, only to update the current value, which gives the name low-storage method. The last two methods are Runge-Kutta methods with an additional property originally called TVD [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF] and more recently SSP [START_REF] Gottlieb | Strong stability-preserving high-order time discretization methods[END_REF]. Without going into too much detail, it means that they are methods that are convex combinations of explicit Euler steps so that their stability is guaranteed for sufficiently small time steps. To us, it means they are defined differently than other methods, using the coefficients α 0<i≤k,0≤j<i and β 0<i≤k,0≤j<i as:

             y n+1 = y (k)
with

y (i) = i-1 j=0 α ij y (j) + ∆tβ ij f y (j) , 1 ≤ i ≤ k
and y (0) = y n .

(5.4) However, they are equivalent to traditional Runge-Kutta methods in the sense that they can be represented by the standard Butcher tableau, with:

             a ij = i-1 l=j+1 α i-1,l-1 a l,j + β i-1,j-1 b i = k l=i+1 α k,l-1 a l,i + β k,i-1 .
(5.5)

We can finally say that we are indeed working with explicit Runge-Kutta methods, as they were introduced initially. The point of this analysis is to see how the global truncation error depends on the time step ∆t that is kept constant throughout the computations. Because some methods are more accurate than others, several computations with different mesh resolutions and Spectral Difference method orders were made. Let us start with the Midpoint method, associated with the order of the spatial discretisation method p = 4. The computation is stopped after 20 periods and the L 2 error for the pressure is analysed for different time steps. Results are shown using the CFL number as it is proportional to the time step: This analysis is repeated for other methods, and we look now at the RK4 method. However, figure 5.5 shows that the error is constant no matter the value of N . It is because the error due to the RK4 method is so small that it is dominated by the error due to the spatial discretisation scheme. It is useless to do the same for lower time steps, as the error would not change, therefore the RK4 error lines are continued by black dashed lines. As the time step decreases, the curves from the Midpoint scheme tend to the same dashed lines: this is proof that this lower bound for the error corresponds to the error of the spatial discretisation method. It shows once again that a spatial discretisation method with small errors is necessary in order to analyse time integration methods. Increasing N up to 128 and then 256 was unsuccessful: the error was still constant. Instead of refining the mesh, an analysis playing with the Spectral Difference method can be performed by increasing the spatial discretisation order. As the theory behind this analysis does not directly depend on N and p, we can adjust them to better fit the method we want to analyse. This is done in figure 5.6 where the global truncation error is shown as a function of the time step, or CFL number, for the RK4, RKo6s and TVDRK(3, 3) methods. The expected slopes for the second and third-order methods are obtained but the result is more subtle for the fourth-order method. The error is once again dominated by the error of the spatial discretisation scheme. On the right part of the curve, one can guess the correct slope but it is necessary to plot the error for larger time step values to see it more clearly. However, those larger time step values are outside the stability domain of the method. In fact, all curves end on the right at the largest time step for which the computation did not fail. It means from figure 5.6 that the RKo6s method is more stable than both the RK4 and TVDRK(3, 3) methods. This is the issue with this analysis: for too small time steps, the error from the spatial discretisation method is dominant and the expected slope is not recovered. On the other hand, with too large time steps, computations fall outside the stability region of the method. The difficulty is to find the correct window that allows us to observe the desired slope by choosing suitable values for p and N . Because the error of the SSPRK(5, 4) is particularly small, it was not possible to find reasonable values below p = 8 and N = 64. Above those values, the computations get rather long so we skipped this method.

The analysis of traditional explicit Runge-Kutta methods can now be applied to the newly added exponential Rosenbrock methods. Their corresponding error curves are shown in figure 5.7. As expected, the respective slope values are 2, 3 and 4. Furthermore, by looking at the abscissa ranges, they work correctly with higher CFL numbers than all the explicit methods that were tested. The goal of this first analysis of the exponential integration methods is not to analyse their robustness, but it already seems better than with explicit methods. What is interesting from this analysis is that one can use higher-order methods with fewer Runge-Kutta stages. Furthermore, the additional computational cost of the additional stages is low. As discussed earlier, exponential Rosenbrock methods can be reformulated as in equation (4.19) using the defects D n,i . The first stage of an exponential Rosenbrock method is then an exponential Rosenbrock-Euler method, for which the contribution of the "full" right-hand side f (y n ) is used. Here, this is done with the algorithms previously described, using a Krylov subspace of dimension 20. Then, because the defects are relatively smaller, the other stages compute their contribution with only 5 Krylov basis vectors. The two exponential Rosenbrock methods with two stages are not a lot more expensive than the Rosenbrock-Euler method. In comparison, a 2-stage explicit Runge-Kutta method will be twice as computationally expensive as the Euler method. This analysis was also performed with larger Krylov subspace methods. In this other configuration, the first stage uses 40 Krylov basis vectors for all three exponential methods, and the second stage uses 10 Krylov basis vectors for the ExpRB32 and ExpRB42 methods. However, using a larger subspace did not have any impact on the result of this analysis. For all (N, p) configurations, the error curve as a function of the time step was almost the same as with smaller Krylov subspace, such as they were almost indistinguishable when plotted in figures like figure 5.7.

ANALYSIS OF EXPONENTIAL TIME INTEGRATION METHODS

We choose this first test case as it is widespread in the computational fluid dynamics community, as a tool to analyse the order of integration schemes. It often concerns spatial discretisation methods, but it can be considered for the time integration method as it gives access to the global truncation error, which is the error that interests the solver users. First, this analysis recovers the correct slope on the error curve as a function of the time step for the already existing explicit Runge-Kutta integration methods. Indeed, the error converges, as the time step decreases, to a minimum error value that comes from the spatial discretisation scheme. By refining the mesh and using higher-order spatial discretisation methods, we are able to reduce this minimal error and observe the expected curve for almost all methods. The same analysis but with the newly implemented exponential Rosenbrock methods showed the expected order. It also showed that the additional cost of using more stages with exponential methods is relatively small compared to the same additional cost for explicit Runge-Kutta methods. Finally, for this convected inviscid isentropic vortex case, the exponential Rosenbrock methods were much more robust, as they work with higher CFL numbers. The analysis of their stability is the topic of the next section.

Robustness analysis: Taylor-Green vortex

The Taylor-Green vortex consists in following the evolution of a set of vortices in a periodic domain. It is a three-dimensional case, solution of the Navier-Stokes equations in a periodic box of length 2πL centered around the origin. The initial flow is given by:

                               u = Ma γr gas T ∞      sin x L cos y L cos z L cos x L sin y L cos z L 0      T = T ∞ P = P ∞ 1 + γ Ma 2 16 cos 2x L + cos 2y L cos 2z L + 2
(5.10) over the domain (x, y, z) ∈ [-πL, πL] 3 , as seen in figure 5.8. Even if u z = 0 at the initialisation, it becomes non-null afterwards and the problem is fully three-dimensional. The initial flow transitions to turbulence: the initial large scales decay into smaller ones, that end up dissipated. Figure 5.8 shows contours for which the second invariant of the velocity gradient, also called Qcriterion, is equal to 0.1. This was chosen as it roughly indicates the sizes of turbulent structures.

The numerical values are non dimensionalised: L = 1, Ma = 0.1, r gas = 71.4284, T ∞ = 1 and P ∞ = 71.4288. The goal is to recover the Reynolds number Re = ρ ∞ U ∞L/µ = 1600 with the associated constant fluid viscosity is µ = 6.25 × 10 -4 . The Prandtl number, the ratio of momentum diffusivity to thermal diffusivity, is equal to 0.71. Using the convective time t c = L/U ∞ , it is known that the maximum dissipation happens around 8t c , and after 15t c the flow is fully turbulent without any trace of the initial structures. The definition of the test case is part of the International Workshop on High-Order CFD Methods. A reference solution is also provided by the workshop, obtained using a spectral solver on a very refined mesh. The Spectral difference method uses the approximated Riemann solver of Roe.

The values of interest for this test case are the temporal evolution of the mean turbulent kinetic energy over the computational domain Ω:

E k = 1 ρ ∞ |Ω| Ω ρ u 2 2 dΩ, (5.11)
the mean enstrophy, defined by The mesh used here is the same for all computations: a regular Cartesian mesh made of 80 3 cells. This choice is not natural from the workshop but previous experiments have led to the conclusion that such a mesh gives the best solution at the lowest CPU cost, for the Spectral Difference method with the order p = 4. Each computation is run on 308 CPU cores. The goal is to see if the result from the simulation matches reference data up to 20t c , and how much CPU time the computation took. The RKo6s, TVDRK(3, 3) and SSPRK [START_REF] Huang | High-order partitioned spectral deferred correction solvers for multiphysics problems[END_REF][START_REF] Coulette | Implicit time schemes for compressible fluid models based on relaxation methods[END_REF] methods are the already existing reference methods, and the exponential Rosenbrock-Euler, ExpRB32 and ExpRB42 methods are the newly implemented methods. For the exponential Rosenbrock methods, Krylov subspaces of dimension 20 are used for the first stage, and dimension 5 for the additional stages for the ExpRB32 and ExpRB42 methods. The time step is constant in a simulation, and some effort was done to use the highest time step compatible with each method. Above this upper bound on the time step, either the computation fails or the quality of the solution is not satisfying anymore. However, this higher time step can be increased for exponential methods by increasing the dimension of the Krylov subspaces used to compute ϕ-functions. By doing so, the computational cost of a single iteration is higher, but it reduces the error in the ϕ-functions evaluations. Symmetrically, it was also tried to reduce the dimension of the Krylov subspaces. It means that the method can not work with the same time step as before as it is now too high, but this reduces the cost of the iteration. This is to see if making more iterations that are cheaper, or in the opposite less that are more expensive can be a good idea. From now on, ExpEuler(d) method refers to the exponential Rosenbrock-Euler method that uses a Krylov subspace of dimension d.

E = 1 ρ ∞ |Ω| Ω ρ ∇ × u
The result of all simulations is shown in figure 5.9. As a simulation corresponds to the computation with the largest time step for which the results are satisfying, this figure does not give much information. At first sight, all methods produce results that are close to reference data. When looking more closely, there are some differences with the reference, but all JAGUAR computations are almost indistinguishable from one another. As a consequence, one can argue that the error from the spatial discretisation method is the origin of the difference between the reference and JAGUAR computations. It is then reasonable to assume that all curves correspond to valid computations. They all give the same result, which is taken as the reference solution for the considered mesh and the set of parameters of the study. Finally, as those methods give similar physical results, the correct way to compare them is to analse their statistics. Statistics for each time integration method are provided in table 5.1. As the cost of the algorithms used in the time integration methods does not change between iterations, it makes sense to compute the average elapsed real time per iteration. Among the explicit methods, we see that the RKo6s and SSPRK(5, 4) methods run at a higher CFL number than the TVDRK(3, 3) method. However, they have more stages per iteration so each iteration takes longer, and the total elapsed real time is rather similar between all three methods. We note that the cost of one step of an explicit Runge-Kutta method is not proportional to the number of stages. Indeed, the implementation of the RKo6s method uses the fact that it is a low storage method to reduce the computational cost in terms of both time and memory. Nevertheless, the cost of an iteration of those explicit methods is much less than exponential methods. As expected, the cost of the exponential Rosenbrock-Euler method is similar to d 2 where d is the dimension of the Krylov subspace used for the Arnoldi iteration. As explained, the ExpRB32 method does not cost a lot more than the exponential Rosenbrock-Euler method that uses Krylov subspace with the same dimension. This is because the second stage uses only a small Krylov subspace, and this is what is seen from the results in table 5.1. We can draw two important conclusions from this analysis. First, we see that the exponential Rosenbrock methods are faster than the explicit Runge-Kutta methods available in JAGUAR. For some of them, they are almost twice as fast. Secondly, we see that the exponential Rosenbrock-Euler method can accurately simulate the Taylor-Green vortex at high CFL numbers. Accurate explicit methods would not work with such high CFL numbers, and even if implicit methods could they would not be accurate enough to capture and preserve the small scales of the flow. To reach high CFL numbers, we have to use larger Krylov subspaces, which makes the method a lot slower. The method used here should not be used for actual applications: the goal was just to show feasibility. It seems that with exponential methods, the robustness can be increased by paying the price in terms of the Krylov subspace dimension. We can think of ways to compute more accurately the ϕ-functions less expensively: by using restarted methods, for instance, we can control the maximal dimension of the Krylov subspace. However, this was not tried during this thesis. This unsteady test case shows that exponential methods are interesting for our applications. We found a method that halves the elapsed real time, which is a great improvement for an unsteady simulation. We also showed that we could accurately simulate this flow while using high CFL numbers. It means that we have added more robust methods than the already existing ones. They could prove useful in applications where we are limited by the lack of stability of the explicit Runge-Kutta methods.

Industrial application: LS89

We decided to try exponential methods on a final test case. Contrary to the previous two, we stepped out of the academic context with an industrial application. It is a simulation of the flow inside the LS89 turbine blade cascade that was studied experimentally by the von Karman Institute [START_REF] Arts | Aero-thermal performance of a two-dimensional highly loaded transonic turbine nozzle guide vane: a test case for inviscid and viscous flow computations[END_REF]. We choose this test case as it has already been done with several solvers, including JAGUAR recently [START_REF] Brunet | Comparison of various CFD vodes for LES simulations of turbomachinery: from inviscid vortex convection to multi-stage compressor[END_REF]. The computation uses a two-dimensional unstructured mesh, made of 2374 cells, that was extruded over 20 cells for a length of 0.15c with c the chord of the turbine blade. The final three-dimensional mesh made of 47 480 cells is shown in figure 5.10. The mesh is periodic in the extruded direction and the vertical direction. The boundary condition at the wall is isothermal, with T wall = 297.75K. At inlet, the stagnation conditions are P 0 = 184 900Pa and T 0 = 409.2K. The outlet pressure is set to 116 487Pa. The upper and lower boundaries, as well as the front and rear planes, correspond to periodic boundary conditions. The Spectral difference method uses the HLLC approximate Riemann solver. This is an unsteady LES simulation: the From the previous results, we decide to use the exponential Rosenbrock-Euler method for the time integration. Indeed, it seems to be the quickest method, and we wish here to reduce the computational time. The other exponential methods would improve the accuracy of the physical results, so we keep only the Rosenbrock-Euler method to simplify the analysis.

The reason why we decided to use exponential methods for this numerical application is that in previous works, the explicit time step was limited by the viscous effect that prevails in the small cells of the boundary layer. Indeed, this test case was used to compare several solvers, but as JAGUAR only had explicit methods it took longer to find a solution [START_REF] Brunet | Comparison of various CFD vodes for LES simulations of turbomachinery: from inviscid vortex convection to multi-stage compressor[END_REF]. Because exponential methods handle well the linear part of equations, and viscous effects are mostly linear, we thought that exponential methods would allow JAGUAR to break free from the constraint on the time step, and improve the overall speed of the solver. In figure 5.11, the Mach number is shown at the end of the computation. It is in line with what is expected, and this proves the capacity of exponential methods to work as time integration schemes for such applications.

The physical value investigated in [START_REF] Brunet | Comparison of various CFD vodes for LES simulations of turbomachinery: from inviscid vortex convection to multi-stage compressor[END_REF] is the heat transfer coefficient. The reader is invited to look at this reference for the interpretation of the physical results. Here, we assume that JAGUAR produces a satisfying solution with its classical explicit TVD(3, 3) method, and we compare it with the exponential Rosenbrock-Euler method. The heat transfer coefficient is shown in figure 5.12, and we see that the results are indistinguishable between the two methods. In this figure, the origin of the curvilinear abscissa corresponds to the leftmost point on the blade, and positive (respectively negative) abscissa corresponds to the upper (respectively lower) part of the blade. This result validates the capability of exponential methods to be used as time integration methods in computational fluid dynamics.

ANALYSIS OF EXPONENTIAL METHODS IN JAGUAR

Our logic was that if the time step for this application is limited by viscous effects, and viscous effects are mainly linear, exponential methods would accept higher time steps, and in the end, the computation would be faster. However, our initial hypothesis was wrong: the time step is limited by advection and not diffusion. It was originally true at the time of some previous investigations [START_REF] Brunet | Comparison of various CFD vodes for LES simulations of turbomachinery: from inviscid vortex convection to multi-stage compressor[END_REF], but as the diffusive scheme was recently greatly improved [START_REF] Marchal | Extension of the spectral difference method to premixed laminar and turbulent combustion[END_REF], it is not true anymore. This is why in the end, there are no significant improvements in terms of computational speed when using exponential integration methods. Even if we did not find the improvements we were originally looking for, this numerical experiment did show that exponential time integration methods can be used in actual computational fluid applications.

Conclusion

The work made towards exponential methods during this thesis shows their interest in the time integration of computational fluid dynamics problems. It is now clear that they are precise methods, and using them makes a lot of sense when the spatial discretisation method is accurate as well, with high-order methods for instance. Then, some numerical experiments with the JAGUAR solver showed their quality. For one application, they were particularly quicker and more stable than classical explicit integration methods. For the last application, they did not improve the solver performances as expected, but it could be argued that we expected too much because of some wrong hypothesis on our part. Overall, this work showed that exponential methods are worth investigating.

Conclusion and perspectives

Jacobian-Free Newton-Krylov for implicit time integration

This thesis was interested in solving efficiently steady problems for multiphysics applications. It means improving the convergence, stability and speed of already existing methods, particularly for the software system CEDRE. A preliminary analysis of existing methods within CEDRE and the literature focussed our effort on improving the nonlinear and linear resolutions required by implicit time integration methods. In particular, the Jacobian-Free Newton-Krylov method caught our interest. This is first because it reuses many algorithms that were already available in CEDRE. The other reason is that it is a better alternative to what was currently done, which is using a poor approximation of the Jacobian matrix. Indeed, using the JFNK method amounts to using a full Jacobian matrix that accounts for all numerical models. We hoped that using a more accurate matrix for the linear problem would improve the quality of the nonlinear solution, which would in turn improve the overall quality of the implicit time integration method. With a careful implementation, it was added to the software system to ensure compatibility with other solvers and existing algorithms.

To check the quality of the Jacobian-Free Newton-Krylov method, it was tested on several applications. Those applications were chosen to represent typical CEDRE computations, to evaluate if the method is interesting for the solver. The first application was the simulation of the turbulent flow around a wing profile in two dimensions. When compared with the traditional method, the new one that is matrix-free gives similar results when looking at flow data. The improvement is solely on the convergence, as the JFNK method converges. The second application is the same wing profile, but this time using a much finer mesh to finely represent the boundary layer. This time, the traditional method lacks the ability to converge satisfyingly when the JFNK method, however, does. It shows that the JFNK method improves the convergence of the solver on such applications.

As CEDRE aims to solve problems in the field of multiphysics, it is not enough to test the new method on aerodynamics test cases. This is why the next chosen application was the simulation of the hypersonic reactive flow around a solid sphere. It is a simplified typical reentry application. Similar results are found when looking at the non-reactive equivalent of this test case: the JFNK method convergence is better. However, when using a fine mesh with the reactive model, the conclusion is no longer in favour of the matrix-free method as the older method converges better. Still, the difference in the convergence between the two methods is small compared to the same difference in the non-reactive case. A possible conclusion is that in the best-case scenario, when the traditional method converges well, the JFNK does not improve the quality of the solution. When the traditional method struggles to find a steady solution, however, the JFNK method becomes interesting.

Finally, the Jacobian-Free Newton-Krylov method was tested on a newly implemented fluid model. Indeed, as CEDRE is under current development, new fluid models are sometimes added to represent flow features with higher fidelity. For instance, the new model used here accounts more precisely for thermodynamic disequilibrium. A newly implemented model often does not give access to their Jacobian matrices, as it would require much more work from developers and is not a priority. It means that users must restrict themselves to using explicit time integration methods. As the JFNK method does not require Jacobian matrices, it is a good candidate to be used as a time integration method with such models. It was then compared with the explicit Midpoint method. In the end, the JFNK method converges better than the explicit method and is also quicker in terms of CPU time. The result is even more in favour of the implicit method when using a simplified version of the reactive model that disregards ionisation.

One key argument of the Jacobian-Free Newton-Krylov method is the possibility to handle the Jacobian matrix approximation in the same order of accuracy as the spatial discretization. We showed that accuracy can lead to improved convergence. Rough approximations of the Jacobian matrices should be avoided in the future, or considered only for tests or as preconditioners. If the choice of CEDRE developers is finally to revert to a Jacobian matrix-based formulation, we can assume that the computation of the Jacobian matrix will be a key ingredient. Doing this by hand is time-consuming and can lead to implementation issues. An alternative could be considering code differentiation. But in the context of CEDRE, it means revisiting the full structure of the whole set of solvers.

The downfall of the Jacobian-Free Newton-Krylov method is its slowness. Indeed, it requires multiple right-hand side evaluations, and such evaluations are quite expensive in our solver. It means that for a user, there is no interest to use the JFNK methods on problems for which the traditional method gives satisfying results and convergence. It is still useful to achieve precise convergence when the traditional method fails to do so. However, improving the efficiency of right-hand side evaluations is an active research topic within the team in charge of the solver. Such improvement will benefit the JFNK method a lot more than it will benefit the traditional method, so the interest in the JFNK method should increase.

In light of this work, it seems natural to want to try the Jacobian-Free Newton-Krylov method with other fluid models that also do not have yet an available Jacobian matrix. Such models already exist in CEDRE, and they are restricted to the explicit time integration method, as was the model detailed in this thesis before it used the matrix-free method. It is reasonable to assume that the JFNK method will yield similar results, and will unlock implicit methods for those models.

Going further, one might want to try the same method but with another solver from the platform that is CEDRE. Indeed, this would amount to very few developments which is why this idea is interesting. Finally, if multiple solvers use the matrix-free method, it means that CEDRE would be able to integrate implicitly those solvers in a coupled fashion. Because of the structural choices made in CEDRE, this is currently impossible unless refactoring a significant amount of code.

As the Jacobian-Free Newton-Krylov method does not use an actual matrix, it is limited in terms of available preconditioners. More generally, preconditioners available in CEDRE are quite simple and the linear solver would benefit from better preconditioners from the literature that already works well in other computational fluid dynamics solvers. Physics-based preconditioner could improve the linear solver performances which would improve robustness, convergence and speed of the overall integration method [START_REF] Park | On physics-based preconditioning of the Navier-Stokes equations[END_REF][START_REF] Liu | CFD high-order accurate scheme Jacobian-Free Newton Krylov method[END_REF]. In particular, for the Multi-TEmperature model one that was used in this thesis, a physics-based preconditioner would help handle the stiff part corresponding to electrons and ionisation.

Exponential integration methods for unsteady time integration with large time steps

Another topic of interest in this thesis was exponential time integration methods. Indeed, such methods reuse many parts and algorithms that are classicly used for implicit time integration. Furthermore, methods used for implicit time integration are often adapted to solve unsteady problems with large time steps. A preliminary analysis showed that the accuracy of exponential methods was similar to the one of explicit methods while being able to use relatively large time steps as implicit methods. As exponential methods are quite precise, they require spatial discretisation methods that are at least as accurate. This is why exponential methods were then studied with the JAGUAR solver that uses a Spectral Difference method.

A first numerical experiment showed that the newly added exponential methods had the expected order of accuracy in JAGUAR, and prove more stable than explicit methods already available in the solver. Then, the same methods were used in another test case. It showed that exponential methods can reduce significantly the CPU time required to fulfil the computation. Furthermore, some methods were able to use time steps 40 times higher than explicit methods and were still quicker to fulfil the computation. This highlight the quality of exponential methods for unsteady time integration when using large time steps. In a final test case, it was shown that exponential methods still work fine on an actual LES computation, with fine cells in a wall boundary condition. Our first hope was that exponential methods would prove faster, but in the end, their performances were similar to the ones of explicit methods. This preliminary work did demonstrate the feasibility of exponential methods in JAGUAR and their suitability for unsteady time integration. Rewriting the exponential computation routines within JAGUAR instead of using the external SLEPc library would improve their quickness even more. Indeed, using a well-written library is often a good idea to benefit from its quality as its developers are experts in their fields, but values are constantly copied back and forth between JAGUAR and SLEPc data types with the current implementation. The solution is either to modify the solver so that it uses SLEPc data types globally, or to rewrite SLEPc methods suited for JAGUAR data types. The latter would also allow for specific optimisations. For instance, we mentioned earlier that exponential methods often compute ϕ-functions of the same matrix. One could reuse spectral information between each of those computations to speed them up and get more accurate results.

As JAGUAR had only explicit integration methods before exponential methods were added, no attention was paid to preconditioning. Now that there are algorithms that use Krylov subspace methods, preconditioners could prove beneficial for the time integration [START_REF] Castillo | Preconditioning the matrix exponential operator with applications[END_REF].

Abstracts Résumé de la thèse en français

Cette thèse s'intéresse aux performances de l'intégration temporelle du code CEDRE sur des problèmes stationnaires. CEDRE est une plateforme logicielle visant la résolution des problèmes multi-fluides pour des applications en énergétique à échelle industrielle. Elle est composée de plusieurs solveurs, chacun dédié à un ensemble de phénomènes physiques. Plus précisément, nous regardons comment améliorer la rapidité, robustesse et convergence de l'intégration temporelle. Pour des raisons de stabilité nous nous intéressons à des méthodes implicites, en particulier à la méthode d'Euler implicite. Ces méthodes nécessitent la résolution de problèmes non-linéaires, qui nécessitent à leur tour la résolution de problèmes linéaires. Le passage de l'un à l'autre se fait par la présence de la Jacobienne des fonctions du problème non-linéaire. Une méthode de Krylov est déjà existante dans CEDRE pour l'inversion de systèmes linéaires: la méthode GMRES. Nous utilisons le fait qu'elle n'a pas explicitement besoin de la matrice pour l'inverser et mettons en place une méthode JFNK. Le but est d'améliorer la précision de la matrice Jacobienne utilisée, en espérant que cela améliorera la précision globale de l'intégration temporelle. Ceci est justifié par le fait qu'avant cette thèse la Jacobienne utilisée est très approximée, notamment en ce qui concerne les modélisations fines des solveurs, comme les termes sources turbulents, et les méthodes de reconstruction, comme les méthodes MUSCL. Une implémentation d'une méthode sans-matrice est mise en place de manière générique de sorte que tout solveur de CEDRE puisse utiliser cette formulation. Cela ouvre de plus la porte à une résolution implicite couplée des solveurs, chose non permise avec la structure actuelle de CEDRE. La méthode JFNK est comparée aux méthodes préexistantes de CEDRE sur des applications typiques de complexité croissantes choisies afin de représenter les fonctionnalités du solveur.

Dans un second temps, nous élargissons le contexte d'étude en nous intéressant aux méthodes d'intégration exponentielles, cette fois avec le solveur JAGUAR. Ce changement de solveur est justifié par la plus grande précision apportée par la méthode des Différences Spectrales utilisée comme schéma de discrétisation spatiale, précision nécessaire à l'analyse de ces nouveaux schémas temporels très précis. Le choix de s'intéresser aux méthodes exponentielle est justifié par le fait que ces méthodes réutilisent beaucoup des ingrédients de l'approche JFNK précédente. Nous choisissons, implémentons et analysons un ensemble de méthodes exponentielles, en comparaison à des méthodes déjà présentes, sur plusieurs cas pour montrer leur intérêt.

128

Thesis abstract in English

This thesis focuses on the performance of the time integration of the CEDRE code on steady problems. CEDRE is a software platform aimed at solving multi-fluid problems for industrialscale energetics applications. It is composed of several solvers, each dedicated to a set of physical phenomena. More precisely, we are looking at improving speed, robustness and convergence of the time integration. For stability reasons, we are interested in implicit methods, in particular the implicit Euler method. These methods require solving nonlinear problems, which in turn requires solving linear problems. The transition from one to the other corresponds to the presence of the Jacobian matrix of functions from the non-linear problem. A Krylov method already exists in CEDRE for the inversion of linear systems: the GMRES method. We use the fact that it does not explicitly need the matrix to invert it and implement a JFNK method. The aim is to improve the accuracy of the Jacobian matrix used, in hope that this will improve the overall accuracy of the time integration. This is justified by the fact that prior to this thesis the Jacobian matrix used is very approximate, especially with respect to many of the fine modelling features of the solver, such as turbulent source terms, as well as reconstruction methods, such as the MUSCL one. A formulation of the matrix-free method is generically implemented so that any CEDRE solver can use this formulation. This also opens the door to coupled implicit solving, something not allowed with the current CEDRE structure. The JFNK method is compared to pre-existing CEDRE methods on typical applications of increasing complexity chosen to represent the functionality of the solver.

In a second step, we broaden the context by looking at exponential integration methods, this time with the JAGUAR solver. This change of solver is justified by the higher accuracy brought by the Spectral Difference method that JAGUAR uses as a spatial discretization scheme, necessary to analyze these new very accurate temporal methods. The choice to focus on exponential methods is justified by the fact that these methods reuse many of the ingredients of the previous JFNK approach. We select, implement and analyse a set of exponential methods, in comparison to existing methods, on several test cases to show their interest.

Méthodologies permettant l'obtention efficace de solutions multi-physiques stationnaires pour des applications en énergétique

Cette thèse s'intéresse aux performances de l'intégration temporelle du code CEDRE sur des problèmes stationnaires. CEDRE est une plateforme logicielle visant la résolution des problèmes multi-fluides pour des applications en énergétique à échelle industrielle. Elle est composée de plusieurs solveurs, chacun dédié à un ensemble de phénomènes physiques. Plus précisément, nous regardons comment améliorer la rapidité, robustesse et convergence de l'intégration temporelle. Pour des raisons de stabilité nous nous intéressons à des méthodes implicites, en particulier à la méthode d'Euler implicite. Ces méthodes nécessitent la résolution de problèmes non-linéaires, qui nécessitent à leur tour la résolution de problèmes linéaires.

Le passage de l'un à l'autre se fait par la présence de la Jacobienne des fonctions du problème non-linéaire. Une méthode de Krylov est déjà existante dans CEDRE pour l'inversion de systèmes linéaires: la méthode GMRES. Nous utilisons le fait qu'elle n'a pas explicitement besoin de la matrice pour l'inverser et mettons en place une méthode JFNK. Le but est d'améliorer la précision de la matrice Jacobienne utilisée, en espérant que cela améliorera la précision globale de l'intégration temporelle. Ceci est justifié par le fait qu'avant cette thèse la Jacobienne utilisée est très approximée, notamment en ce qui concerne les modélisations fines des solveurs, comme les termes sources turbulents, et les méthodes de reconstruction, comme les méthodes MUSCL. Une implémentation d'une méthode sans-matrice est mise en place de manière générique de sorte que tout solveur de CEDRE puisse utiliser cette formulation. Cela ouvre de plus la porte à une résolution implicite couplée des solveurs, chose non permise avec la structure actuelle de CEDRE. La méthode JFNK est comparée aux méthodes préexistantes de CEDRE sur des applications typiques de complexité croissantes choisies afin de représenter les fonctionnalités du solveur.

Dans un second temps, nous élargissons le contexte en nous intéressant aux méthodes d'intégration exponentielles, cette fois avec le solveur JAGUAR. Ce changement de solveur est justifié par la plus grande précision apportée par la méthode des Différences Spectrales utilisée comme schéma de discrétisation spatiale, précision nécessaire à l'analyse de ces nouveaux schémas temporels très précis. Le choix de s'intéresser aux méthodes exponentielle est justifié par le fait que ces méthodes réutilisent beaucoup des ingrédients de l'approche JFNK précédente. Nous choisissons, implémentons et analysons un ensemble de méthodes exponentielles, en comparaison à des méthodes déjà présentes, sur plusieurs cas pour montrer leur intérêt.

Mots-clés : METHODE NUMERIQUE ; ENERGETIQUE ; MULTIPHYSIQUE ; STATIONNAIRE

Efficient methods for steady multi-physics simulations

This thesis focuses on the performance of the time integration of the CEDRE code on steady problems. CEDRE is a software platform aimed at solving multi-fluid problems for industrial-scale energetics applications. It is composed of several solvers, each dedicated to a set of physical phenomena. More precisely, we are looking at improving speed, robustness and convergence of the time integration. For stability reasons, we are interested in implicit methods, in particular the implicit Euler method. These methods require solving nonlinear problems, which in turn requires solving linear problems. The transition from one to the other corresponds to the presence of the Jacobian matrix of functions from the non-linear problem. A Krylov method already exists in CEDRE for the inversion of linear systems: the GMRES method. We use the fact that it does not explicitly need the matrix to invert it and implement a JFNK method. The aim is to improve the accuracy of the Jacobian matrix used, in hope that this will improve the overall accuracy of the time integration. This is justified by the fact that prior to this thesis the Jacobian matrix used is very approximate, especially with respect to many of the fine modelling features of the solver, such as turbulent source terms, as well as reconstruction methods, such as the MUSCL one. A formulation of the matrix-free method is generically implemented so that any CEDRE solver can use this formulation. This also opens the door to coupled implicit solving, something not allowed with the current CEDRE structure. The JFNK method is compared to pre-existing CEDRE methods on typical applications of increasing complexity chosen to represent the functionality of the solver.

  .10), G can be linearised in Ξ. With y = Ξ -Ξ and J = ∂ G ∂Ξ Ξ , assumed constant, it comes:

Figure 1 . 1 :

 11 Figure 1.1: Stability regions (in colour) of the first fourth-orders Runge-Kutta methods.

1. 3 . 21 Figure 1 . 2 :

 32112 Figure 1.2: Stability regions (in colour) of the first four Adams-Bashforth methods.

Figure 1 . 4 :

 14 Figure 1.4: Procedure used to find the solution of steady problems.

  (a) GT01R: 2D inviscid flow in the inter-blade channel of a linear cascade turbine. (b) HV15R: 3D RANS simulation of an engine fan. (c) RM07R: 3D viscous flow in a jet engine compressor.

Figure 1 . 5 :

 15 Figure 1.5: Matrices from various CFD simulations [33] using a Finite Volume method. Coloured points correspond to nonzero values.

Figure 1 . 6 :

 16 Figure 1.6: Spectrum of A in blue and A pre = AD -1 in orange (left). GMRES convergence for those matrices (right).

Figure 2 . 2 :

 22 Figure 2.2: Error in the Jacobian matrix-vector product approximation as a function of ε (left) and of the dimension N (right). On the left figure, coloured cross markers correspond to the values computed on a mesh of 10 1 , 10 2 , ... cells, and circle markers to the last value with 10 7 cells. Grey cross markers on the left figure correspond to a mesh of 10 7 cells.

Figure 2 . 3 :

 23 Figure 2.3: Error in the Jacobian matrix-vector product approximation as a function on ε (left) and of the dimension N (right). The function is the right-hand side of equation (1.5) from a Finite Volume method using a centred scheme for the Euler equations over a regular one-dimensional mesh. On the left figure, coloured cross markers correspond to the values computed on a mesh of 10 1 , 10 2 , ... cells, and circle markers to the last value with 10 7 cells. Grey cross markers on the left figure correspond to a mesh of 10 7 cells.

Figure 3 . 1 :

 31 Figure 3.1: Mesh for the RAE 2822 test case. Close up on the leading edge and on the expected shock location.

Figure 3 . 2 :

 32 Figure 3.2: Pressure around the RAE 2822 airfoil and sonic Ma = 1 contour.

Figure 3 . 3 :

 33 Figure 3.3: Pressure coefficient around the airfoil for the traditional method (blue) and the JFNK method (orange).

Figure 3 . 4 :

 34 Figure 3.4: Aerodynamic coefficients (lift and drag) for the profile throughout the computation.

Figure 3 . 5 :

 35 Figure 3.5: Residual 1-norms for the two momentum components, turbulent viscosity and energy throughout the computation for the turbulent transonic airfoil case. Values are normalised by the initial residual.

Figure 3 . 6 :

 36 Figure 3.6: Finer mesh for the RAE 2822 test case. Close up on the leading and trailing edges.

Figure 3 . 7 :

 37 Figure 3.7: Residual 1-norms for the horizontal and vertical momentum components, turbulent viscosity and energy throughout the computation for the turbulent transonic airfoil with boundary layer resolution case. Values are normalised by the initial residual.

3. 1 . 57 Figure 3 . 8 :

 15738 Figure 3.8: Residual ∞-norms for the horizontal and vertical momentum components, turbulent viscosity and energy throughout the computation for the turbulent transonic airfoil with boundary layer resolution case. Values are normalised by the initial residual.

Figure 3 . 10 :

 310 Figure 3.10: Mesh along with final pressure, temperature, and NO and e mass fractions for the hypersonic reactive sphere test case.

Figure 3 . 11 :

 311 Figure 3.11: Streamlines near the stagnation point for two Riemann solvers.

Figure 3 . 12 :

 312 Figure 3.12: Residual 1-norms for the horizontal and vertical momentum components, N 2 volumic mass and energy throughout the computation for the hypersonic non-reactive sphere case. Values are normalised by the initial residual.

Figure 3 . 13 :

 313 Figure 3.13: Residual 1-norms for the horizontal and vertical momentum components, N 2 volumic mass and energy throughout the computation for the hypersonic reactive sphere case. Values are normalised by the initial residual.

Figure 3 . 14 :

 314 Figure 3.14: Temperature on the symmetry axis for the Multi-TEmperature model.

Figure 3 . 15 :

 315 Figure 3.15: Residual 1-norms for the electron mass fraction, vertical momentum, N 2 vibrational energy and total energy throughout the computation for the hypersonic reactive sphere case using the Multi-TEmperature model. Values are normalised by the initial residual.

Figure 3 . 16 :

 316 Figure 3.16: Residual 1-norms for the electron mass fraction, vertical momentum, N 2 vibrational energy and total energy throughout the computation for the hypersonic reactive sphere case using the Multi-TEmperature model without ionisation. Values are normalised by the initial residual.

CHAPTER 4 .

 4 INTRODUCTION TO EXPONENTIAL INTEGRATION METHODSa fixed linearisation: L = f (y 0 )[START_REF] Hochbruck | Explicit exponential Runge-Kutta methods for semilinear parabolic problems[END_REF]. It is because they were developed for semilinear parabolic problems: dy dt + Ay = r (y) .(4.16)

Figure 5 . 1 :

 51 Figure 5.1: Continuous function (left) and discontinuous representation made by the Spectral Difference method (right).

Figure 5 .

 5 Figure 5.2 shows the steps the Spectral Difference make to compute the flux divergence ∇•F (u) from the solution u, with the specific choice p = 2. 0. At first, the solution is represented with a red line by its value in the p + 1 solution points.The solution points are marked by red dots. Parts of the solution from the left and right neighbouring cells can be seen. As the figure shows, the piecewise local solutions inside each cell may be discontinuous at the mesh interface. This is the starting point of the Spectral Difference method and will be called the 0-th step.

Figure 5 . 2 :

 52 Figure 5.2: Steps of the Spectral Difference methods with p = 2.

Figure 5 . 3 :

 53 Figure 5.3: Solution points and flux points in the [0, 1] segment used in the Spectral Difference method for 0 ≤ p ≤ 11.

. 3 )

 3 They are traditionally defined inside the [-1, 1] segment but are scaled into [0, 1]. They are the roots of the Chebyshev polynomials of the first kind and are often used in polynomial interpolation as they tend to minimise Runge's phenomenon. The p + 2 flux points are the p roots of the p-th Legendre polynomials to which are added the two segment end points. This choice has the property that there is a flux point between each contiguous solution point. There is no explicit formula for the Legendre polynomial roots as there is one for the Chebyshev polynomial roots. They are also defined in the [-1, 1] segment and are scaled into [0, 1]. Finally, the solution and flux points can be seen in figure5.3.

Figure 5 . 4 :

 54 Figure 5.4: Global truncation error for the RK2 Midpoint method.

9 )

 9 Using the CFL number instead of the time step allows comparing results with different N or even p. It can be seen here as a nondimensional time. Computations are made for multiple values of N : 16, 32 and 64, and the corresponding results are shown in figure5.4. The 16 × 16 mesh does not allow for a good analysis: the error due to the spatial discretisation method is too important and prevents seeing the expected dependency between the global truncation error and the time step. The two other meshes however show the expected relation: err = O ∆t 2 .

Figure 5 . 5 :

 55 Figure 5.5: Global truncation error for the RK2 and RK4 method.

Figure 5 . 6 :

 56 Figure 5.6: Global truncation error for the RK4, RKo6s and TVDRK(3, 3) methods.

Figure 5 . 7 :

 57 Figure 5.7: Global truncation error for the exponential Rosenbrock-Euler, ExpRB32 and Ex-pRB42 methods.

  by the present simulation are compared with the reference values.

Figure 5 . 8 :

 58 Figure 5.8: Isocontour (0.1) of the Q-criterion coloured by velocity norm of the Taylor-Green vortex at various times.

Figure 5 . 9 :

 59 Figure 5.9: Kinetic energy and enstrophy for various time integration methods compared to reference data.

Figure 5 . 10 :

 510 Figure 5.10: Mesh for the LS89 turbine blade. Close up on the leading and trailing edges.

Figure 5 . 11 :

 511 Figure 5.11: Mach number in the LS89 turbine blade cascade.

Figure 5 . 12 :

 512 Figure 5.12: Heat transfer coefficient from the TVDRK(3, 3) and exponetial Rosenbrock-Euler methods. Negative abscissa corresponds to the lower part of the blade and positive abscissa to the upper part.

  arguments: the type_sys instance, an index identifying the vector to multiply and an index identifying the vector where to put the result. As we explained when we introducedCHAPTER 2. DESCRIPTION OF THE DEVELOPMENT WORKFLOW IN CEDREKrylov subspace methods, GMRES uses the linear system matrix through this matrix-vector product routine. Classically, a client solver such as CHARME fills the matrix coefficients and then lets BIBCEDRE solve the linear system. In order to use the matrix-free approximation from equation(1.53), one only needs to replace this routine with a new one that computes the approximation. Unfortunately, the approximation uses a function that belongs to the client solver. The library BIBCEDRE does not know this function and how to compute it, as it is part of CHARME or any other client solver. As we said, we want to write generic solutions, so merging the library BIBCEDRE with the solver CHARME is not a good way to proceed. What we need here is to allow BIBCEDRE to use a callback from the client solver. We did that using the Fortran 2003 feature: polymorphism. The Fortran 2003 standard is almost twenty years old today, and many may consider that the functionalities added to the previous Fortran 95 standard are therefore not that innovative. But in the context of a project as large as CEDRE, developed by researchers and used by industrial partners, even the Fortran 2003 standard and polymorphism in particular are innovative. Without going into too much detail, we added a member to the type type_sys that contains the context to evaluate a matrix-vector product:

		i_p, i_ap)
	type(type_sys), intent(inout) :: sys
	integer,	intent(in)	:: i_p
	integer,	intent(in)	:: i_ap
	that takes three

type type_sys ! Matrix-vector product context class(type_gmvec_ctx), pointer :: gmvec_ctx => null() ... ! Additional data end type type_sys with: type type_gmvec_ctx procedure(interface_gmvec), pointer, nopass :: gmvec ... ! Additional data end type type_gmvec_ctx

  e v,m ∂ t (ρ e e e ) + ∇ • ((ρ e e e + p e ) u) = ∇ • λ e ∇T e + ρh e D e ∇y e + u • ∇p e

+ S e-t + S e-r + S e-v + ωe e e ∂ t (ρE) + ∇ • ((ρE + p + p e

Table 5 .

 5 1: Statistics for various time integration methods. The time corresponds to elapsed real time, or wall-clock time. A computation corresponds to the simulation of the Taylor-Green vortex from equation (5.10) on time interval [0, 20t c ].

Remerciements

Part II

Unsteady time integration using large time steps

method is unstable and the computation fails. The implicit method is able to handle higher CFL numbers, of course as it is the reason implicit methods are used. But unsurprisingly we see in the bottom left figure that the solution is even worse. Finally, the exponential method is able to handle the higher CFL number and preserve the vortex well enough. It has the robustness of the implicit method and the precision of the higher-order explicit method.

The analysis of this simple test case is quite rudimentary. Indeed, we did not look finely at the error introduced by the time integration methods, but we visually compared the results. This work showed that at the cost of very few developments we have a method that could prove interesting for time integration with large time steps. Its purpose was not to precisely analyse exponential methods but to see their feasibility in our solver. Now that this preliminary test showed the quality of exponential time integration methods, we need to perform an actual analysis.

We see in figure 4.1 that even the methods that preserve the vortex, the exponential and the RK4 methods, do introduce some error. Indeed, the vortex after one period is not exactly the same as at the beginning, as it should be. We see that it loses in intensity, and the contour lines are not exactly concentric. Those symptoms are seen in both computations with the exponential and the higher-order explicit methods. Furthermore, the error in the result of those simulations looks similar. It is because this error does not come from the time integration method but from the spatial discretisation method. Indeed, the spatial discretisation method induces some error when evaluating the right-hand side of equation (4.1). This error is accumulated throughout the simulation, and because of it, we could not get the exact result no matter the time integration method. What we can do, however, is take a spatial discretisation method such as the error it adds is small compared to the error from the time integration scheme. This way, the error we get at the end of the computation is mostly due to the time integration method, and it allows us to compare methods. Here the spatial discretisation method used in CEDRE is a second-order method. We need to use a higher-order spatial discretisation method to perform the analysis of exponential time integration methods.

CHAPTER 5. ANALYSIS OF EXPONENTIAL METHODS IN JAGUAR

Order analysis: convected inviscid isentropic vortex

This first case is designed to analyse the order of accuracy of the proposed methods and to validate their implementation. It is the simple case of the two-dimensional inviscid isentropic vortex, the same as the one considered with CEDRE. The case is still described by equation 4.28, but the numerical values are different. The convective flow Mach number is still 0.5, but this time P ∞ = 1 × 10 5 Pa and T ∞ = 300K. The heat capacity ratio is γ = 1.4 and the specific gas constant is r gas = 287.058J kg -1 K -1 . Those values correspond to standard values for dry air. The mesh represent the two-dimensional box [0, L]

2 with L = 0.1m, and is made of

The vortex is defined by its characteristic radius R 0 = 0.005m and its intensity β = 0.2. It corresponds in fact to the same vortex as the one used with CEDRE, in size and in intensity, but it has been scaled down into a smaller box and dimensioned to some given pressure and temperature. The present choice is in agreement with the prescription of the International Workshop on High Order CFD Methods. For this application, the Spectral difference method uses the approximated Riemann solver of Roe.

The period for those numerical values is T = L/U ∞ = 5.760 × 10 -4 s. After 20 periods, the vortex should recover the initial position and the 2-norm of the error between initial and computed solutions can be defined for any scalar variable u (x, t):

Once the vortex moves in the domain, the numerical solution is subject to two schemes and errors are a consequence of both. Once the spatial scheme is defined, the mesh is sufficiently refined in order to transport the vortex accurately and the goal is to have a much lower influence of the spatial scheme than of the time integration scheme on the total error. The order of our time integration methods can be determined by looking at this error as a function of the time step.

The error at the end of the simulation is the global truncation error: the error the method makes while getting to a fixed time, no matter how many iterations it took. However, the order of a time integration method was previously defined using the local truncation error: the error it makes after a single small step. There is a link between global and local truncation errors for our single-step methods. Let us call τ n the local truncation error and e n the global truncation error at step n. Because single-step methods can be written as:

with an increment function g that is K-Lipschitz continuous in the y variable, the global truncation error is bounded [START_REF] Süli | An introduction to numerical analysis[END_REF]:

A method is a p-order method if τ n = O ∆t p+1 . Therefore, a p-order method verify e n = O (∆t p ). However, the reciprocal is not true: observing a global truncation error e n = O (∆t p ) does not mean the order of the method is p. As a consequence, attention is focused on the expected order for the global truncation error at the end of the computation, as it is the error users look at. The proof of the order for a method is analytical, by writing the partial Taylor series of the exact solution, as we did earlier for the exponential Rosenbrock-Euler method.

CONCLUSION AND PERSPECTIVES

The last JAGUAR application that was presented here was originally chosen as it was limited in the time step size by the diffusive part of the equations, which is mostly linear. Exponential methods would have then improved the time integration as they deal efficiently with the linear parts. However, this initial statement was not true, and therefore the gain from using exponential methods did not appear. It would be interesting to see what would happen on a test case for which this statement holds.

Going outside the scope of JAGUAR, exponential methods could prove to be great candidates as time integration methods on other solvers. For instance, solvers that simulate plasmas and arc lightning are significantly limited in the time step size by viscous effects. As such effects are mostly linear, exponential integration methods would allow breaking free from this limitation. The same idea applies to the field of thermic for which large physical times need to be simulated. Exponential methods should be compatible with larger time steps, which would reduce the overall computational time.

Exponential methods are becoming more and more active in computational fluid dynamics. For the same reason we decided to use them along a Spectral Difference method, they are often used with high-order spatial discretisation methods. For instance, the Discontinuous-Galerkin method from [START_REF] Li | A parallel discontinuous Galerkin method with physical orthogonal basis on curved elements[END_REF] is used with an exponential Rosenbrock method described in [START_REF] Li | An exponential time-integrator scheme for steady and unsteady inviscid flows[END_REF][START_REF] Li | Time advancement of the Navier-Stokes equations: p-adaptive exponential methods[END_REF]. In [START_REF] Narayanamurthi | Partitioned exponential methods for coupled multiphysics systems[END_REF], exponential integrators are used to solve multiphysics problems. Many other recent references from the literature use exponential methods as time integrators for similar computational fluid dynamics problems. This shows that exponential methods are worth investigating.

Appendix A

Conservative preconditioning

Let us consider the ordinary differential equation:

This equation comes from the Finite Volume discretisation, where Q is the vector of the conservative variables and F a conservative right-hand side. To keep this simple, we consider that there is a single conservative variable, so that the dimension of Q is equal to the number of cells. This is only to simplify the notations, but everything can be adapted to use more conservative variables. The mass matrix M is a diagonal matrix where the ith diagonal coefficient corresponds to the volume of the ith cell. We note by S the vector of the same dimension as Q made of ones and transposed. This way, multiplying on the left a vector by S amount to sum a vector components. The conservation property of the equation means that S F (Q) = 0 or equivalently that the sum of the conservative variable over the domain SM Q is constant.

When solving this equation with the explicit Euler method, we have at the nth step that:

and so the explicit Euler method preserves the conservation property.

When using the implicit Euler method with a single linearisation, we need the Jacobian matrix J n of the function F evaluated in Q n . Since J n = F (Q n ) and S F (Q n ) = 0, we have SJ n = 0 also. The method gives the increment δQ n as the solution of:

Multiplying this relation on the left by S and using that SJ n = S F (Q n ) = 0, we have SM δQ n = 0, which means this method also preserves the conservation property.

However, we do not usually use the exact solution of the linear problem but the solution of a subspace Krylov method. Let us consider that we use k steps of a Krylov subspace method with a zero initial guess. Then, the increment belongs to the corresponding Krylov subspace:

Now, there are no reasons for SM δQ n to be equal to zero. For example, the solution computed with the smallest Krylov subspace dimension (k = 1) is parallel to F (Q n ), and then SM δQ n ∝ SM F (Q n ) and therefore is not null.

With preconditioning however, we can recover the conservation property. If we use the invert of the mass matrix as a preconditioner, either a left or a right one, we have that:

We can verify by recurrence that SM v i = 0:

Then, we have that SM δQ n = 0 which means the preconditionned method preserves the conservation property.