# Contributions to the numerical simulation of some fluid dynamics models

Stéphane del Pino

► To cite this version:

Stéphane del Pino. Contributions to the numerical simulation of some fluid dynamics models. Mathematics [math]. Sorbone Université, 2023. tel-04264442v2

HAL Id: tel-04264442

https://cea.hal.science/tel-04264442v2

Submitted on 12 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Habilitation à diriger des recherches en mathématiques

# Contributions to the numerical simulation of some fluid dynamics models

Stéphane del Pino

Soutenue le 14 Novembre 2023 après avis des rapporteurs

**Jean-Luc Guermond**

Professeur, Texas A&M University, Department of Mathematics

**Claus-Dieter Munz**

Professeur, Universität Stuttgart, Institut für Aerodynamik und Gasdynamik

**Nicolas Seguin**

Directeur de Recherche, INRIA, Centre de l'Université Côte d'Azur, antenne de Montpellier

devant le jury

**Daniel Bouche**

Directeur de Recherche, CEA, DAM, DIF

**Bruno Després**

Professeur, Sorbonne Université, Laboratoire Jacques-Louis Lions

**Jean-Luc Guermond**

Professeur, Texas A&M University, Department of Mathematics

**Claus-Dieter Munz**

Professeur, Universität Stuttgart, Institut für Aerodynamik und Gasdynamik

**Frédéric Nataf**

Directeur de Recherche, Sorbonne Université, Laboratoire Jacques-Louis Lions

**Nicolas Seguin**

Directeur de Recherche, INRIA, Centre de l'Université Côte d'Azur, antenne de Montpellier
<figure><img src="image_1.png" /></figure>
Acknowledgements

First of all, I would like to express my warmest thanks to J.-L. Guermond, C.-D. Munz and N. Seguin who have accepted to review this manuscript and for the glowing reports they wrote. I would like to thank them again, and to thank D. Bouche, B. Després and F. Nataf for agreeing to be part of my jury.

Lots of thanks to Bruno for all the help you gave me in preparing and organizing this defense. Thank you!

I also warmly thank my hierarchy at CEA who encouraged me to prepare this HDR for so long. Chronologically, Bruno S., Hervé J., Pierre-Franck P., Stéphane J. and Cédric E., thank you. I want to thank Laurence B. and Daniel B. for the organization of the Séminaire HDR which triggered the writing of this manuscript.

Obviously, this kind of work is not a lonely adventure, so I want to thank all my coauthors and especially Silvia B., Gilles C., Tomás C., Bruno D., Pascal H., Hervé J., Emmanuel L., Frédéric L., Bertrand M., Isabelle M., Olivier P., Alexiane P., Ulrich R. and Driss Y.

My colleagues play a special role in these achievements. They prove me every day that one can work seriously in a friendly (and funny) atmosphere. Thank you Christophe B., Patricia C., Gilles C., Teddy C., Jean-François C., Frédéric D., David D., Nicolas G., Mathieu G., Agnès G., François H., Philippe H., Stéphane J., Pierre K., Mohammed K., Emmanuel L., Franck L., Benoît L., Isabelle M., Jérôme M., Marie-Pierre O., Tristan P., Alain P., Céline P. and Xavier V.

I am also grateful to my more distant friends and colleagues at CEA. Let me thank Céline B. (remember "*le petit pont de bois*"), Sonia B., Jean-Philippe B., Jérôme B., Alexandra C., Samuel K., Pierre-Henri M. and Alexiane P. It is always a pleasure to have scientific (or not) discussions with you.

I want to thank all my friends and colleagues from the outside world and especially Rémi A., Xavier B., Laurent B., Frédéric C., Bruno D., Laurent D., Pascal H., Frédéric H., Mourad I., François J., Frédéric L., Raphaël L., Hélène M., Bertrand M., Olivier P., Ulrich R., Nicolas S., Stéphanie S., Nicolas V., Driss Y.

I have a special thought for Cécile D. who left us far too soon and who inspired some of the works presented in this memoir.

Thanks to the PhD students at Ter@tec. It is always nice to have people that kindly smile at my worst jokes. Thank you Axelle D., Victor F. and Julie P. (ok you are not really there anymore). I want to especially thank the PhD students that I co-supervise(d) Axelle D., Alexiane P. and Driss Y. It has been (and it is) a pleasure to work with you.

Many thanks to the CEA "football" team. It is always refreshing and nice to play with you.

Benoît and Christophe, it is not easy to manage careers as both an engineer and a break manager. Thank you for your help and support! And more seriously, just thank you!

I shall thank "mes amis du Barrault", always Philippe, often Jérôme and Stéphane, and sometimes François, Pascal and Xavier. We spend some incredible moments "à refaire les matches" without any judgment, and discussing of all the rest (even if mathematics are never far away).

I would like to thank you Patricia, for all these scientific discussions full of humility and knowledge. Thank you also for your kindness and your laughter. We have to work together. I'm a little sad to no longer be your git reference and I'm still waiting to taste your homemade *cannelés*.

Thank you Gilles for your friendship, your energy, your kindness and your determination. I hope we'll be working together again soon. If I put aside the hairstyle, you are an example for young people.

A little word to my "conscrit" Emmanuel. As time goes, I have the strange feeling that we have more and more projects to go. Thank you for everything and understanding my wildness.

Philippe, my lifelong friend, thank you for everything. Even if it is not likely to happen, do not change. We have so much to do and so much to laugh about.

I cannot end these words of thanks without adding a final one. Thank you, Isabelle. We have been sharing the office for a few years now. Thank you for all these moments of hard work and laughter. We still have a long way to go in these two areas. I cannot wait!
# Contents

Contents ii

Notations iii

Introduction 1

1 Elliptic problems and incompressible flows 3

1.1 Fictitious Domain Methods 3

1.2 Incompressible flows 6

2 Hyperbolic problems and compressible flows 9

2.1 High-order schemes for the acoustic wave equation and applications to aeroacoustics 9

2.2 Multi-dimensional finite-volume methods for compressible flows in semi-Lagrangian coordinates 13

2.3 Mesh adaptation for semi-Lagrangian compressible flows 34

Conclusion and perspectives 41

Bibliography 43

Publications 43

References 44
# Notations

For the sake of uniformity the following notations are used in this manuscript.

## Bibliographic citations

All along the manuscript, I use the following convention to distinguish between my publications and other references.

* Publications I co-authored are cited using numeric reference and regrouped in the bibliographic section starting from page 43.

* Other references use an alphabetic style and are regrouped from page 44.

## Vectors and matrices

To ease the reading, we use the following classical typographic conventions.

* Vector quantities are denoted using bold faces.

* Matrix quantities are written using capital letters.

Let $V$ be a vector space, for readability and depending on the context, the dot product of two vectors $\mathbf{u}, \mathbf{v} \in V$ is either denoted by $\mathbf{u} \cdot \mathbf{v}$, $(\mathbf{u}, \mathbf{v})$ or $\mathbf{u}^T \mathbf{v}$.

## Integrals and volumes

In order to simplify the notations, when there is no ambiguity, the measure is omitted in integrals. For instance

$$ \int_{\Omega} f = \int_{\Omega} f d\mu, $$

where $(\Omega, \mu)$ is a measurable space.

Also, in order to simplify the discourse, we will refer as

$$ |\Omega| = \int_{\Omega} d\mu, $$

the volume of $\Omega \in \mathbb{R}^d$, whatever the value of the dimension $d$ is, and $|\partial\Omega|$ is called a surface.
## Mesh notations

With regard to meshes and their elements, the following notations are used.

• $\mathcal{J}$: the set of all the cells of a mesh $\mathcal{M}$,

• $\mathcal{L}$: the set of all the edges of $\mathcal{M}$,

• $\mathcal{R}$: the set of all the nodes of $\mathcal{M}$.

Also, to ease the reading, to denote generic elements of these sets, we write

• $j \in \mathcal{J}$: a cell of $\mathcal{M}$ (k might also be used when needed),

• $l \in \mathcal{L}$: an edge,

• $r \in \mathcal{R}$: a node or a vertex of $\mathcal{M}$ (s is sometimes used).

Mixing these notations allows to define some useful subsets of $\mathcal{J}$, $\mathcal{L}$ or $\mathcal{R}$. Let us give some examples

• $\mathcal{J}_r$: the set of cells of $\mathcal{M}$ that are connected to the vertex $r$,

• $\mathcal{L}_j \cap \mathcal{L}_r$: the set of edges that are connected to the cell $j$ and to the node $r$.

• $\mathcal{R}_j \cap \mathcal{R}_k$: the set of vertices that are connected to cell $j$ and to cell $k$.

Finally, in this document, we use the following convenient abuse of notation. When referring to the cell $j$ we do not distinguish between

• the domain of $\mathbb{R}^d$ that is defined by the cell,

• the mesh connectivity element,

• or the number (index) associated to the cell.

Thus, we write the volume of a cell $j$ as $V_j := \int_j 1$, for instance.

The same kind of abuse of notations is used for edges and vertices.
Introduction

This manuscript is a summary of the research works that I have done since my PhD Thesis. It is composed of two main chapters of non-equal sizes.

The first and shorter one summarizes my work in the context of elliptic problems and incompressible flows. The second and larger chapter presents more precisely my work dealing with hyperbolic problems and compressible flows. Actually, there is no direct link between these two topics, but the influence of the techniques that I used for elliptic problems is not negligible with respect to some of the works presented in the second chapter. For instance, the semi-Lagrangian finite-volume solver designed to treat curvilinear cells [9] is based on a finite-element approximation of the fluxes.

Before describing into more details the content of this document, let me summarize a few works that are not discussed in the following. The first one, see [12], was published during my PhD Thesis. In this work, with E. Heikkola, O. Pironneau and J. Toivanen, we propose a methodology to solve Helmholtz equation in 3D in geometries provided by Constructive Solid Geometry. The second work that is not described in this manuscript has been published with J.-L. Lions and O. Pironneau in [15] during my PhD Thesis. Here, we study a domain decomposition technique, a modified Schwarz algorithm, that is used to analyze the Chimera method [Ste91]. An application to geological flows is presented. The last work that is not described in this manuscript was published in [8] as a result of my first postdoctoral position at CEA. It provides an efficient method to visualize Cartesian tree-based AMR data in 2D or 3D.

The first chapter deals with elliptic problems and incompressible flows. It is composed of two main sections. The first section is dedicated to some contributions to Fictitious Domain Methods (FDM). After a brief introduction that recalls the interest of FDM in the elliptic context, I remind some results that I provided during my PhD Thesis. I also describe rapidly some related works performed with O. Pironneau. Then, I recall a simple application of the penalty method that we proposed with B. Maury to deal with fluid-structure interaction. My PhD work was actually the starting point of an extension to spectral methods that we performed with D. Yakoubi during his PhD Thesis. The second section is dedicated to incompressible flows. I summarize a theoretical joint work with D. Yakoubi and U. Razafison, where we provide a lower bound to the inf – sup constant for the divergence operator. Finally, I present the numerical analysis of a non-linear iterative method designed to solve an ocean-atmosphere turbulent coupling model. This is the result of a collaboration with T. Chacón Rebollo and D. Yakoubi.

The second chapter recalls my contributions to the numerical analysis and to the design of schemes to approximate hyperbolic problems and compressible flows. These works are more related to my activities at CEA. Firstly, I present the research I made during my second postdoctoral position under the supervision of H. Jourdren. Actually, we proposed an arbitrary high-order finite-difference scheme to approximate the
advection equation at constant velocity. We adapted it to define an Eulerian scheme (Lagrange+remap) that behaves well in transport or acoustic regimes. Then I recall the application of this arbitrary high-order scheme, that I performed with B. Després, P. Havé, H. Jourdren and P.-F. Piserchia, to aeroacoustics. Then, I switch to one of my main topics: the numerical resolution of Euler equations in semi-Lagrangian coordinates. I recall the 3D and second-order accurate extensions of the finite-volume scheme proposed by B. Després and C. Mazeran. This work was published with G. Carré, B. Després and E. Labourasse. Then I present some works dedicated to the extension of these semi-Lagrangian schemes to higher-order than two. A first contribution is the result of researches conducted with G. Carré and E. Labourasse. Then I recall a numerical method that I proposed to deal with curvilinear cells in the context of finite-volume semi-Lagrangian coordinates. I show how we applied it, with P. Hoch and E. Labourasse, to obtain promising results with regard to defining a very high-order semi-Lagrangian scheme. Then using the same ideas, I recall the numerical method that we proposed to deal with perfect sliding, with S. Bertoluzza and E. Labourasse. Then the work performed with E. Labourasse and G. Morel is summarized. It consists in the definition of an asymptotic preserving scheme designed to approximate the Scannapieco-Cheng model using a multidimensional indirect ALE framework. In the next section, I depict the PhD Thesis work of A. Plessier, that I co-supervise with B. Després. In this work, we propose and analyze an implicit semi-Lagrangian scheme to approximate gas dynamics in 1D and 2D. The scheme is unconditionally stable. This chapter is concluded by presenting a triangular metric-based mesh adaptation method that we proposed with I. Marmajou. The method is quite efficient, since it yields a quasi-Lagrangian method<fnref n="1" /> that preserves an adapted mesh all along the calculation.

<footnote n="1">¹At each time step, most of the cells are preserved by the mesh adaptation procedure, and are Lagrangian.</footnote>
# Elliptic problems and incompressible flows

In this chapter, I present some works dedicated to the approximation of elliptic problems and incompressible flows. These have been done mostly in parallel with my researches at CEA.

## 1.1 Fictitious Domain Methods

Fictitious Domain Methods were first introduced by C. Peskin [Pes72] (Immersed boundary method) in the seventies to compute viscous incompressible flows around heart valves using finite-differences. This approach consists in solving the flow on a rectangular grid (which does not fit the geometry of the computational domain) and to take into account the boundary conditions using a well chosen set of Dirac masses inside the domain. This leads to an efficient treatment of the evolution of the domain's geometry that does not require the use of any kind of remeshing technique.

Since then, many Fictitious Domain Methods have been proposed (the most popular one being probably the one by V. Girault and R. Glowinski [GG95] which imposes Dirichlet boundary conditions by means of Lagrange multipliers).

As one could expect, Fictitious Domain Methods have generally a low order of approximation (namely $O(h^{1/2})$ in the case of [GG95]), however this is balanced by the fact that the rectangular background meshes allows the use of fast solvers in $D$ (Fast Fourier Transform [CT65], cyclic-reduction [Swa77], PSCR [RT99],...) which makes the approach competitive especially when dealing with moving boundaries.

Let us describe the model problem that will be used in the following. Let $\Omega$ be a connected bounded domain of $\mathbb{R}^d$ (with $d \in \{1, 2, 3\}$) such that its boundary $\partial\Omega$ is regular enough. Let $f \in L^2(\Omega)$ and let $u$ be the unique solution of

$$ \left| \begin{array}{l} \text{Find } u \in H_0^1(\Omega), \\ -\Delta u = f. \end{array} \right. \qquad (1.1) $$

<figure><img src="image_2.png" /><figcaption>Figure 1.1: Example of a fictitious domain geometry. The domain $D$ embeds the computational domain $\Omega$.</figcaption></figure>
Let now $D \supset \Omega$ (see Figure 1.1) and let $f^D \in L^2(D)$ be an $L^2$ extension of $f$ ($f^D|_{\Omega} = f$ a.e. in $\Omega$). Then the “fictitious domain problem”, find $u^D \in H^1(D)$ such that

$$ \left| \begin{array}{l} -\Delta u^D = f^D, \\ u^D|_{\partial\Omega} = 0, \end{array} \right. \qquad (1.2) $$

admits a unique solution which satisfies $u^D|_{\Omega} = u$ in $H^1(\Omega)$. Fictitious Domain Methods mostly reduce to formulating precisely (1.2) (defining $f^D$ outside of $\Omega$ and how to impose $u^D|_{\partial\Omega} = 0$) and to provide an efficient discretization of it.

Many fictitious domain approaches have been developed, interested reader can check [PO01; Mau01; BIM11; GK98; Hei+98; GPP94] and the references therein.

## Finite-Element Methods

When O. Pironneau proposed to supervise my PhD Thesis [7], he had two main ideas in mind: study a simple fictitious domain-like method and use this method to build a general 3D PDE solver (for elliptic and parabolic problems) à la freefem<fnref n="1" />. Actually, freefem was already a successful tool used not only for teaching purpose but also for research. The choice of a fictitious domain based method was actually quite important in that context, since it allowed to avoid the complex task of 3D mesh generation<fnref n="2" />.

As a Fictitious Domain Method we decided to use a penalty approximation of Dirichlet boundary conditions, so that in the simple case of the Poisson problem, it reads find $u_h \in V_h$ such that

$$ \forall v \in V_h, \quad \int_D \mathbf{1}_\Omega \nabla u_h \cdot \nabla v_h + \frac{1}{\epsilon} \int_{\partial\Omega} u_h v_h = \int_D \mathbf{1}_\Omega f v_h, \qquad (1.3) $$

where $V_h$ is a finite-element subspace of $H^1(D)$. The main difficulty when computing (1.3) is the calculation of $\int_{\partial\Omega} u_h v_h$ which requires a surface mesh to define the numerical integration. This kind of mesh is actually much simpler to obtain than a conformal 3D finite-element one, since one can use a marching cube-like algorithm [LC87] to define it<fnref n="3" />. One also needs to take care of the integration of the volume terms in the cells that contain a portion of $\partial\Omega$'s mesh.

The main theoretical result I provided is an error estimate which states that if the solution $u$ of (1.1) belongs to $H^{k+1}(\Omega)$, one has

$$ \|u - u_h\|_{1,\Omega} \le C_1 (1 + C_2 h^k) \sqrt{\epsilon} \left\| \frac{\partial u}{\partial n} \right\|_{0,\Gamma} + C_3 (h^k \|u\|_{k+1,\Omega}), \qquad (1.4) $$

where $C_1, C_2$ and $C_3$ are positive constants that do not depend on $h$ and $\epsilon$. One should observe that (1.4) shows that the method has the same order of accuracy as classical finite-element approximation in the case of an exact quadrature. A sketch of the proof is given in [18] (see [7] for details), it strongly relies on [Bab73], where I. Babuška studies the approximation of Dirichlet boundary conditions by penalty. In [Mau08], B. Maury establishes the convergence of the solution of a penalized abstract problem to the solution of the Dirichlet problem. In this case of a distributed penalty<fnref n="4" />, the convergence rate is proved to be $O(h^{1/2}) + O(\epsilon^{1/2})$.

<footnote n="1">¹freefem (https://freefem.org) is a successful PDE solver whose C++ kernel is driven by a flexible DSL (domain specific language) which allows to implement with ease complex algorithms using a strong numerical toolbox.</footnote>

<footnote n="2">²In 1999, NETGEN[Sch97] and TetGen[Si00] were very new tools and not yet easy to use. The now popular Gmsh[GR09] was not released yet.</footnote>

<footnote n="3">³Since the mesh is only used to compute quadrature, its generally bad quality does not deteriorate the conditioning of the obtained linear system</footnote>

<footnote n="4">⁴Distributed penalty Poisson problem: find $u^\epsilon \in H^1(D)$ s.t. $\forall v \in H^1(D), \int_\Omega \nabla u^\epsilon \cdot \nabla v + \frac{1}{\epsilon} \int_{D \setminus \Omega} \nabla u^\epsilon \cdot \nabla v = \int_\Omega f v$.</footnote>
In a first proceeding [6], I described the method and some choices I made in order to implement FreeFEM3D. In a second proceeding [18], written with O. Pironneau, we published the analysis of the method (discussing also the easier case of Neumann or Fourier boundary conditions).

We later published an application of the method [19]. It was used to perform the *couplex exercice*: a simulation challenge designed by ANDRA<fnref n="5" /> to check the ability of the simulation of nuclear waste storage. The difficulty of this test resides in the treatment of the different scales (for instance, the computational domain is 25km × 695m × 300m, and porosity coefficients may vary by a factor 10⁷ alongside of the materials,...). We perform an asymptotic analysis of the problem which leads to a domain decomposition technique that splits the multi-scale problem into a collection of simpler ones.

I will close this section by describing a work that was performed with B. Maury and published on the occasion of O. Pironneau's 60th birthday [17]. The aim of this work is to propose a simple method to compute 2D/3D turbines. Let $\mathcal{O} \subset \Omega$ (the turbine domain). Let us denote by $\mathbf{U}$, a rigid body velocity field. Actually, the solution in the sense of distributions, $\mathbf{u}^\epsilon$ of the penalty Stokes problem

$$ \left| -\mu \Delta \mathbf{u}^{\epsilon} + \nabla p^{\epsilon} = \mathbf{f} + \xi^{\epsilon}, \quad \nabla \cdot \mathbf{u}^{\epsilon} = 0, \right. \quad \text{where} \quad \xi^{\epsilon} = \frac{1}{\epsilon} \mathbf{1}_{\mathcal{O}}(\mathbf{u}^{\epsilon} - \mathbf{U}) $$

converges, when $\epsilon \to 0$, to the solution of the Stokes problem

$$ \left| -\mu \Delta \mathbf{u} + \nabla p = \mathbf{f} + \xi, \quad \nabla \cdot \mathbf{u} = 0, \right. \quad \text{where } \xi \text{ is the force required to impose } \mathbf{u} = \mathbf{U} \text{ on } \partial \mathcal{O}. $$

Starting from this remark, we define an algorithm such that, given a rotation axis, one computes the angular velocity of an obstacle $\mathcal{O}$. We then solve the penalty Navier-Stokes problem discretizing the convective terms by means of the method of characteristics (as introduced by O. Pironneau [PLT92]). Simulations were performed with both freefem++ and FreeFEM3D.

## Spectral Methods

In 2007, I started to co-supervise, with C. Bernardi, the PhD Thesis of D. Yakoubi [Yak07]. Due to the optimal quality of approximation of the Fictitious Domain Method (1.3), I suggested to D. Yakoubi that similar results could be obtained in the case of spectral methods [BM97]. The objective was to benefit from spectral approximation in non-tensorial domains<fnref n="6" />.

In the case of the Poisson problem, the method reads exactly as (1.3), but the space of approximation $V_h$ is spanned by a tensorial Legendre basis. We showed similarly to the finite-element case, optimal approximation order (with a more elegant formulation than the one I produced during my PhD Thesis). Let $u \in H^m(\Omega)$ solution of the Poisson problem, if $\partial\Omega$ is regular enough,

$$ \begin{align} \|u - u_{\epsilon}^{\delta}\|_{H^1(\Omega)} &\le c \left( N^{1-m} \|f\|_{H^{m-2}(\Omega)} + \sqrt{\epsilon} \|\frac{\partial u}{\partial n}\|_{L^2(\partial\Omega)} \right), \tag{1.5} \\ \|u - u_{\epsilon}^{\delta}\|_{L^2(\Omega)} &\le c \left( N^{-m} \|f\|_{H^{m-2}(\Omega)} + \epsilon \|\frac{\partial u}{\partial n}\|_{L^2(\partial\Omega)} \right), \notag \end{align} $$

<footnote n="5">⁵Agence nationale pour la gestion des déchets radioactifs (French national radioactive waste management agency).</footnote>

<footnote n="6">⁶We wanted to avoid the use of conformal transformations, which limits accessible geometries, or the complexity of spectral element-like methods.</footnote>
where N is the degree of the polynomial approximation in each direction.

Unfortunately, we never published this work. The reason for that is that the numerical experiments were almost impossible to produce: actually the estimates holds for $\partial\Omega \in C^{m-1,1}$ and in the case of exact integration. This conducted us to define an octree to compute the integrals on $\Omega$. Obviously the size of this 3D octree grows exponentially which makes the method far too expensive$^7$. In my experience, this was an educative example of a numerical method that has very good properties on the paper but that is not practicable.

## 1.2 Incompressible flows

During D. Yakoubi's PhD Thesis [Yak07], we produced two additional works. The first one is the result of a collaboration with U. Razafison [20] where we provide a new estimate for the inf – sup condition's constant. The second one, in collaboration with T. Chacón Rebollo [5], studies a resolution method for a model that describes the coupling of ocean-atmosphere flows.

### A lower bound for the inf – sup condition's constant for the divergence operator

In the short note [20], published with U. Razafison and D. Yakoubi, we provide a lower bound to the inf – sup condition's constant for the divergence operator. The evaluation of such a constant is of practical interest since it may be useful for a posteriori error evaluations (see [HSV12] for instance). Let $\omega \in \mathbb{R}^d$, $d = 2,3$ and let us denote by $b_\omega(\cdot, \cdot)$ the bilinear form

$$ \forall (\mathbf{u}, p) \in H^1(\omega)^d \times L^2(\omega), \quad b_\omega(\mathbf{u}, p) = - \int_\omega p \nabla \cdot \mathbf{u}. $$

Denoting by $L_0^2(\omega) := \{q \in L^2(\omega) \text{ s.t. } \int_\omega q = 0\}$, we recall the definition of the inf – sup condition constant $\beta(\omega) > 0$:

$$ \beta(\omega) := \inf_{q \in L_0^2(\omega)} \sup_{\mathbf{v} \in H_0^1(\omega)^d} \frac{b_\omega(\mathbf{v}, q)}{\|\mathbf{v}\|_{H^1(\omega)^d} \|q\|_{L^2(\omega)}}. $$

We showed that, $\forall \Omega \supset \bar{\omega}$ connected open set of $\mathbb{R}^d$ with a Lipschitz-continuous boundary, the following inequality holds:

$$ \beta(\omega) \geq \frac{\beta(\Omega)}{(1 + P_{\Omega})(1 + \|R_{\omega}\|)} $$

where $P_{\Omega} > 1$ is Poincaré's constant on $\Omega$ and $R_{\omega}$ is the harmonic trace lifting operator on $\omega$. One should observe that Poincaré's constant can be difficult to estimate, and thus should be of concern in the choice of a suitable $\Omega$.

### An iterative procedure to solve a coupled two-fluids turbulence model

In this paper [5], published with T. Chacón and D. Yakoubi, we propose and analyze an iterative method to solve a system of two stationary turbulent flows (namely the ocean and the atmosphere) coupled at the interface by means of turbulent kinetic energy. This kind of model is a simplified version$^8$ of models often used in geophysics [Lew97; LTW93].

<footnote>$^7$Actually it was so expensive that we were not able to conduct a reasonable convergence study.</footnote>

<footnote>$^8$Convection is neglected.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1.2: Computational domain for the coupled two-fluids turbulence model.</figcaption></figure>

The computational domain $\Omega \in \mathbb{R}^d$, $d \in \{2,3\}$ is partitioned into two subdomains $\Omega_i$, $i \in \{1,2\}$ that share the interface $\Gamma = \partial\Omega_1 \cap \partial\Omega_2$. One denotes $\Gamma_i = \delta\Omega \setminus \delta\Omega_i$, see Figure 1.2.

The model reads

$$ \forall i \in \{1,2\}, \quad \left| \begin{array}{l} -\nabla \cdot (\alpha_i(k_i)\nabla \mathbf{u}_i) + \nabla p_i = \mathbf{f}_i, \\ \nabla \cdot \mathbf{u}_i = 0, \quad \text{in } \Omega_i, \\ -\nabla \cdot (\gamma_i(k_i)\nabla k_i) = \alpha_i(k_i)|\nabla \mathbf{u}_i|^2, \end{array} \right. \quad (1.6a) $$

with the following boundary conditions

$$ \forall i \in \{1,2\}, \quad \left| \begin{array}{ll} \mathbf{u}_i = \mathbf{0}, & \text{on } \Gamma_i, \\ k_i = 0, & \text{on } \Gamma_i, \end{array} \right. \quad (1.6b) $$

and on the fluids interface $\Gamma$, one imposes

$$ \forall i,j \in \{1,2\} \text{ s.t. } i \neq j, \quad \left| \begin{array}{l} \alpha_i(k_i) \partial_{\mathbf{n}_i} \mathbf{u}_i - p_i \mathbf{n}_i + \kappa_i (\mathbf{u}_i - \mathbf{u}_j) |\mathbf{u}_i - \mathbf{u}_j| = \mathbf{0}, \\ k_i = \lambda |\mathbf{u}_1 - \mathbf{u}_2|^2. \end{array} \right. \quad (1.6c) $$

In this system, $\mathbf{u}_i, p_i$ and $k_i$ denote respectively, the velocity, the pressure and the turbulent kinetic energy of fluid $i$. The quantity $\alpha_i(k_i)$ is its turbulent viscosity and $\gamma_i(k_i)$ its turbulent diffusion. Finally $\mathbf{f}_i$ denotes the external force applied to the fluid $i$. The parameters $\kappa_i > 0$ and $\lambda > 0$ are assumed to be constant.

This stationary problem is highly non-linear. On the one hand, this is due to the functions $\alpha_i, \gamma_i$ and to the kinetic energy production term $\alpha_i(k_i)|\nabla \mathbf{u}_i|^2$. On the other hand, the interface boundary conditions are themselves non-linear.

To approximate solutions of (1.6) we use the following iterative procedure that consists in separating the velocities and pressures calculations from the turbulent kinetic energies ones. Thus one solves at each step $n$

$$ \left| \begin{array}{ll} -\nabla \cdot (\alpha_i(k_i^n) \nabla \mathbf{u}_i^{n+1}) + \nabla p_i^{n+1} = \mathbf{f}_i, & \text{in } \Omega_i, \\ \nabla \cdot \mathbf{u}_i^{n+1} = 0, & \text{in } \Omega_i, \\ \mathbf{u}_i^{n+1} = 0, & \text{on } \Gamma_i, \\ \alpha_i(k_i^n) \partial_{\mathbf{n}_i} \mathbf{u}_i^{n+1} - p_i^{n+1} \mathbf{n}_i + \kappa_i (\mathbf{u}_i^{n+1} - \mathbf{u}_j^{n+1}) |\mathbf{u}_i^{n+1} - \mathbf{u}_j^{n+1}| = \mathbf{0}, & \text{on } \Gamma, \end{array} \right. \quad (1.7a) $$

followed by

$$ \left| \begin{array}{ll} -\nabla \cdot (\gamma_i(k_i^n) \nabla k_i^{n+1}) = \alpha_i(k_i^n) |\nabla \mathbf{u}_i^{n+1}|^2, & \text{in } \Omega_i, \\ k_i^{n+1} = 0, & \text{on } \Gamma_i, \\ k_i^{n+1} = \lambda |\mathbf{u}_1^{n+1} - \mathbf{u}_2^{n+1}|^2, & \text{on } \Gamma. \end{array} \right. \quad (1.7b) $$
We show that under classical assumptions of regularity on the domains $\Omega_i$ ($\Omega_i$ is convex or $\partial\Omega_i \in C^{1,1}$) and for small enough $f_i$ that the problem (1.6) admits a unique smooth solution. In that case, the iterative procedure described by (1.7) converges and $\lim_{n \to +\infty} (\mathbf{u}_i^n, p_i^n, k_i^n)$ is the unique solution of (1.6).

In practice, the non-linear boundary condition in (1.7a) is replaced by

$$ \alpha_i(k_i^n) \partial_{n_i} \mathbf{u}_i^n - p_i^n \mathbf{n}_i + \kappa_i (\mathbf{u}_i^{n+1} - \mathbf{u}_j^n) |\mathbf{u}_i^n - \mathbf{u}_j^n| = \mathbf{0}, \quad \text{on } \Gamma. $$

The proof of convergence of this algorithm is much more difficult<fnref n="9" /> to establish, but one can show that if it converges, it converges to the solution of (1.6).

I will end this paragraph by noticing that T. Chacón Rebollo and D. Yakoubi continued recently this work in [CY18] where they regularized the interface conditions to overcome some lack of regularity.

<footnote n="9">⁹The proof of convergence of the iterative procedure (1.7) is already quite technical.</footnote>
# Hyperbolic problems and compressible flows

In this chapter, I summarize my researches with regard to the approximation of solutions of hyperbolic systems of conservation laws and in particular to compressible gas dynamics.

## 2.1 High-order schemes for the acoustic wave equation and applications to aeroacoustics

In 2004, I started my second postdoctoral position at CEA under the supervision of H. Jourdren. We published [13] an arbitrary high-order scheme for the approximation of solutions of the linear advection equation in 1D. These fluxes were then used to build a scheme for Euler's equations in 1D, that is very efficient in acoustics and transport regimes.

### Arbitrary high-order scheme for strictly hyperbolic linear systems in dimension 1

The initial objective of this study was to derive an accurate scheme to deal with acoustic waves propagating on very long distances. H. Jourdren proposed that a good starting point was the paper of V. Daru and C. Tenaud [DT04]. Indeed, in this paper, the authors provide a collection of one step schemes up to the seventh-order of accuracy (both in space and time) to approximate linear advection at constant velocity,

$$ \left| \begin{array}{l} \partial_t u + a \partial_x u = 0, \quad \text{with } a \in \mathbb{R}, \\ \text{with } u(\cdot, 0) = u_0(\cdot), \quad \text{the initial condition.} \end{array} \right. $$

The scheme construction is quite simple, it relies on the calculation of the equivalent equation of a $n^{th}$-order scheme and on the discretization of the error term to reach order $n + 1$. The procedure is initiated by the upwind scheme<fnref n="1" />. The high-order in time discretization is obtained through the Cauchy-Kovalevskaïa procedure which consists, in that case, in using

$$ \partial_t^k u = (-a)^k \partial_x^k u, \quad \forall k \in \{1, \dots, p\} \quad \text{if } u_0 \in C^p(\mathbb{R}), $$

to substitute time derivatives into space derivatives in the equivalent equation.

<footnote n="1">¹The second-order scheme is nothing else but the Lax-Wendroff [LW60] scheme in the case of linear advection</footnote>
Using the procedure (by upwinding corrections for odd orders and centering corrections for even orders), I noticed a pattern and showed that the following one step scheme is Nth-order accurate in both space and time in the case a > 0 on uniform grids:

$$u_j^{n+1} = u_j^n - \nu (F_{j+1/2}^N - F_{j-1/2}^N), \qquad (2.1a)$$

where $\nu = a \frac{\Delta t}{\Delta x}$ is the Courant number, and where the fluxes $F_{j+1/2}^N$ are defined recursively by

$$\left| \begin{aligned} F_{j+1/2}^1 &= u_j^n, \\ F_{j+1/2}^N &= F_{j+1/2}^{N-1} - \frac{1}{N!} \left( \prod_{i=-m, i \neq 0}^{M} \nu + i \right) \left( \sum_{k=0}^{N-1} (-1)^{k+N} \binom{N-1}{k} u_{j+m-k}^n \right), \end{aligned} \right. \qquad (2.1b)$$

with $m = [\frac{N}{2}]$, $M = [\frac{N-1}{2}]$ and $\binom{n}{p} = \frac{n!}{p!(n-p)!}$. One obtains a similar formula in the case $a < 0$.

Actually this scheme was already defined by B. P. Leonard in [Leo91], but if he provides the procedure to explain its derivation (and the schemes up to order 7²), he does not write the recursive formula (2.1b), which is a novelty of our work. We also show that the scheme is stable in the sense of von Neumann for standard CFL condition $\nu < 1$, but the proof was so technical that we did not published it. However in [Des08], B. Després shows the same result in an elegant way<fnref n="3" />.

Obviously, the schemes defined by (2.1) provide straightforwardly Nth-order schemes for strictly hyperbolic linear systems in dimension 1. Indeed, this kind of systems reads

$$\partial_t \mathbf{u} + A \partial_x \mathbf{u} = \mathbf{0}, \qquad (2.2)$$

where $A \in \mathbb{R}^{q \times q}$ is diagonalizable with distinct real eigenvalues $\lambda_1 < \dots < \lambda_q$. Denoting $\mathbf{r}_k$, a kth right-eigenvector associated to $\lambda_k$ (i.e. $\forall k \in \{1, \dots, q\}$, $A\mathbf{r}_k = \lambda_k \mathbf{r}_k$), one classically rewrites (2.2) as a system of decoupled linear advection equations

$$\forall k \in \{1, \dots, q\}, \quad \partial_t w_k + \lambda_k \partial_x w_k = 0, \qquad (2.3)$$

where

$$\mathbf{u} = \sum_{k=1}^{q} w_k \mathbf{r}_k. \qquad (2.4)$$

A simple linearity argument allows to approximate the solution of (2.2) to the Nth-order: one uses (2.1) to approximate the q independent advection problems (2.3) at order N and simply reconstructs u by (2.4).

Actually, the acoustic wave equation is just a special case of (2.2)

$$\left| \begin{aligned} \partial_t u + \frac{1}{\rho_0} \partial_x p &= 0, \\ \partial_t p + \rho_0 c_0^2 \partial_x u &= 0, \end{aligned} \right.$$

where $\rho_0, c_0 \in \mathbb{R}_+^*$ are the density and the sound velocity of the propagating medium. In that case, the two Riemann invariants associated to $(u, p)$ are $w_\pm = u \pm \frac{1}{\rho_0 c_0} p$. They satisfy the two decoupled advection equations $\partial_t w_\pm \pm c_0 \partial_x w_\pm = 0$.

<footnote n="2">²Actually, due to a different decentering choice, the scheme given by V. Daru and C. Tenaud [DT04] matches (2.1) up to the fourth-order but not beyond that.</footnote>

<footnote n="3">³In this paper, he also showed asymptotic $L^1$ and $L^\infty$ stability, proving then the convergence of these schemes for BV initial data.</footnote>
The scheme (2.1) is very efficient. It is arbitrary high-order accurate (dealing with smooth solutions) and very cheap. However, it is limited to strictly hyperbolic problems with constant coefficients.

Thus, in [13], we also propose a scheme to approximate the compressible gas dynamics system in dimension 1. It is based on the following discussion. In both the acoustic and transport regimes, it can be relevant to produce a scheme that behaves like (2.1). To achieve this, we solve the Euler system of equations using a Lagrange plus remap strategy. Aside from the splitting of the Euler system, the isentropic Riemann invariants $J_{\pm} = u \pm \int \frac{dp}{\rho c}$ are used. In the case of polytropic gases they simplify to $J_{\pm} = u \pm \frac{2}{\gamma-1}c$.

I do not detail the scheme in this document and I invite the interested reader to consult [13].

It is worth noting that, even if the obtained scheme is not arbitrary high-order, it behaves quite well in the transport and acoustic regimes. This scheme is used in the HERA [Jou05] code with the acoustic invariants ($u \pm \frac{1}{\rho c}p$ which are a "linearized" versions of $J_{\pm}$). It has also been used successfully in [HJJ09] to study the convergence of high-order schemes when dealing with non-convex equations of state. Finally, this work has been extended by F. Duboc et al. [Dub+10] to the 6th-order of approximation of the non-linear Euler equations using a finite-volume formulation.

### Application to aeroacoustics

As stated in the beginning of this paragraph, the aim of the study was to propose a numerical method to compute accurately acoustic waves through long distances. In a joint work with B. Després, P. Havé, H. Jourdren and P.-F. Piserchia, we published in [11] an illustration of the efficiency of the scheme we developed.

**The Attenborough test.** The first test of interest has been defined by K. Attenborough et al [Att+95], where the Authors provide an analytical solution to a wave propagation problem in a domain with varying sound speed. We adapted the problem to the 2D case. It is defined as follows. Let $\Omega = [0, 5000[ \times [0, 4000[$ be the computational domain. The initial pressure and velocity is set to 0. In $\Omega$, the density is constant $\rho = 1.205$ and the sound speed is a function of $y$, $c(y) = 343.23 + 0.1y$. Symmetry conditions are imposed to the boundaries $x = 0$ and $y = 0$. A pressure source $p_s(t) = \sin(20\pi t)$ is set at position $(0, 5)$, see Figure 2.1.

The aim of this test is to compute the pressure absorption along the axis defined by $y = 1$.

The numerical setting consists in using directional splitting and to use in each direction the scheme (2.1). Since the sound velocity is constant in each layer of the x direction, the scheme is high-order in x. However, it is only second-order accurate in the y direction since the sound velocity growth linearly with y. Moreover, the use of a Strang directional splitting [Str68] limits the global order of accuracy to 2.

The numerical results are shown in Figure 2.2. On the left part (Figure 2.2a), we compare the results obtained for different schemes. The schemes denoted by "linearized Lax-Wendroff" correspond to (2.1). The "GAD" and the "Hybrid-Godunov" (GoHy) schemes are presented in [HJJ09]. The scheme denoted by "SDP-HP" consists of a third-order

<figure><img src="image_4.png" /><figcaption>Figure 2.1: Attenborough test case configuration.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>(a) Comparison of various schemes on a fine grid (15 cells per wavelength)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Results for various orders (2.1) on a coarse grid (5 cells per wavelength)</figcaption></figure>

Figure 2.2: Comparison of the obtained pressure absorption (dB) for various schemes along the y = 1 axis.

scheme defined for the non-constant in space advection equation

$$ \partial_t u + c(x) \partial_x u = 0. $$

This scheme was not published for two reasons: the fluxes were too complex and if the scheme improved slightly the results obtained by the third-order linearized one, it was very expensive.

Calculations are run using 15 cells per wavelength. The best result is given by the linearized 17th-order⁴ scheme (qualitatively and quantitatively quite close to the analytical solution). One notices that the second-order GAD scheme is too dissipative. However, if the level provided by the SDP-HP scheme is bad, the frequency of the signal is quite accurate. Finally, if the signal level obtained by the GoHy scheme is quite good<fnref n="5" />, a phase problem is observed.

On Figure 2.2b, the effect of the linearized scheme order is depicted using a coarser grid (5 cells per wavelength). One can observe that the 17th-order scheme gives pretty good results even on such a coarse grid.

Aside from the good results provided by this family of schemes it is important to notice that these high-order schemes are very cheap: the CPU cost of the 17th-order scheme is only five times the one of the first-order scheme. This is due to cache effects: many calculations are done with a few data.

**The Misty Picture experiment.** The second valuable test is the simulation of the Misty Picture event [Leh87] that held in the desert of New Mexico during the 80's. A chemical explosion of approximately 4kT of TNT was performed. It generated a long range acoustic wave around 0.1 Hz. For the simulation, the sound speed and the winds are modeled as functions of the altitude. The effects of the winds are taken into account by adding a convective part to the acoustic wave equation

$$ \left| \partial_t p + \mathbf{a} \cdot \nabla p + \rho c \nabla \cdot \mathbf{u} = 0, \\ \rho (\partial_t \mathbf{u} + \nabla \cdot (\mathbf{u} \otimes \mathbf{a})) + \rho c \nabla \cdot \mathbf{u} = 0, $$

<footnote n="411">⁴¹¹ was chosen arbitrarily by Pascal Havé: “One has to stop at some point!”.</footnote>

<footnote n="5">⁵My personal analysis is that the GoHy scheme is not stable (limitation was necessary to run the simulations). In my view, it can be an explanation to the good absorption level.</footnote>
where

$$ \mathbf{a}(x,y,z) = \begin{pmatrix} a_x(z) \\ a_y(z) \\ 0 \end{pmatrix}, \quad \rho(x,y,z) = \rho(z), \quad \text{and} \quad c(x,y,z) = c(z), $$

are constant in time functions. The numerical treatment follows exactly what was done for the previous test. In the splitting, the convective terms are discretized using the scheme (2.1). The computational domain is very large $] -1000, 1000]^2 \times ]0, 200[$ (in kilometers). Except at the top of the domain (which uses absorbing boundary conditions to simulate the end of the atmosphere), symmetry boundary conditions are imposed. On a 2-dimensional experiment, we showed that 15 cells per wavelength were enough to reach a satisfactory convergence using the 17th-order scheme. At this time, the simulation

<figure><img src="image_7.png" /><figcaption>(a) Pressure coloring at time 2520s. The zoom illustrates the complexity of the pressure field on the ground.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Comparison of experimental measures (black) and our results (in blue). The red curves have been computed with a ray tracing technique.</figcaption></figure>

<caption>Figure 2.3: Results obtained for the Misty-Picture experimentation.</caption>

was very large: 4092 processors of the TERA-10 supercomputer [Lah06] were used, the calculation itself involved several billions of cells.

## 2.2 Multi-dimensional finite-volume methods for compressible flows in semi-Lagrangian coordinates

In 2005, I got a permanent position at CEA. I integrated the development team of the massively parallel multi-physics 3D code Troll [Lef+18]. As many radiation-hydrodynamics codes dedicated to ICF<fnref n="6" />, it is based on a semi-Lagrangian<fnref n="7" /> discretization of the equations. From the numerical point of view, it means that the mesh is displaced at each time step to follow the fluid flow.

One year later, B. Scheurer and B. Després wanted to implement a 3D version of the Glace [DM05; Maz07] scheme<fnref n="8" /> into Troll.

<footnote n="6">⁶Inertial Confinement Fusion</footnote>

<footnote n="7">⁷Also referred as updated-Lagrangian in the literature. Actually, the semi-Lagrangian formulation is nothing else than the Eulerian formulation written on arbitrary domains moving at the flow velocity. The Lagrangian formulation on the other hand formulates the Euler equations on the initial domain.</footnote>

<footnote n="8">⁸The Glace scheme is a semi-Lagrangian finite-volume scheme that approximates Euler equations in arbitrary dimension. It will be presented in the following.</footnote>
Their idea was to evaluate the behavior of this new scheme for realistic applications. Thus, with G. Carré, B. Després and E. Labourasse, we implemented the scheme, extending the PhD Thesis work of C. Mazeran [DM05; Maz07] to the third dimension and to the second-order of accuracy.

Since then, Lagrangian hydrodynamics has been one of my main topics of research. This is the subject of this section.

In order to fix notations, let us first recall Euler equations in dimension $d$, which have to be taken in the sense of distributions on an open set $\Omega \subset \mathbb{R}^d$

$$ \left| \begin{array}{l} \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, \\ \partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla p = \mathbf{0}, \\ \partial_t (\rho E) + \nabla \cdot (\rho E \mathbf{u}) + \nabla \cdot (p \mathbf{u}) = 0, \end{array} \right. \qquad (2.5) $$

where $\rho$ is the mass density, $\mathbf{u}$ is the fluid velocity and $E$ is the specific total energy. The system is closed thanks to the equation of state that defines the pressure $p = p(\rho, \epsilon)$, where $\epsilon$ denotes the specific internal energy ($\epsilon = E - \frac{1}{2}\|\mathbf{u}\|^2$). Moreover in order to select physical weak solutions of (2.5), one must ensure physical entropy production (in the sense of distributions)

$$ \partial_t (\rho \eta) + \nabla \cdot (\rho \eta \mathbf{u}) \geq 0, \qquad (2.6) $$

where $\eta$ is the physical entropy defined by the Gibbs relation $Td\eta = d\epsilon + pd\tau$. Here we introduced $T := \frac{\partial \epsilon}{\partial \eta}|_{\tau}$, the fluid temperature and the specific volume $\tau = \rho^{-1}$. For the sake of simplicity, we deliberately omitted boundary conditions.

The system (2.5) is **strictly hyperbolic** as soon as the equation of state satisfies

$$ c^2 := \left. \frac{\partial p}{\partial \rho} \right|_{\eta} > 0, $$

which defines the sound speed $c$.

## Euler equations in semi-Lagrangian coordinates

The use of Lagrangian methods to simulate hydrodynamics is not new. Actually, the first method defined to compute hydrodynamics shocks has been published by J. von Neumann and R. D. Richtmyer in 1950 [NR50]. This scheme, generally referred as the vNR scheme, is a Lagrangian method. It has been successfully used to compute 1D flows since then, and remains a reference method that is still a source of inspiration. Nowadays, it enters in the family of staggered in space and time schemes. It is morally second-order accurate in both space and time since all derivatives are almost centered. Actually, the scheme is constructed on the discretization of isentropic Euler equations and the Authors introduced the notion of artificial viscosity to impose entropy production in shocks<fnref n="9" />. For these reasons, the scheme is very efficient, especially on coarse grids.

A few years later, S. K. Godunov developed another method [God59], based on conservation principles, to simulate gas dynamics in 1D. This method is by construction first-order accurate and the unknowns are the mean values of conservative variables (even the total energy) in each cell. Historically, this is the first finite-volumes scheme. The numerical fluxes are given by the exact resolution of the Riemann problems that are defined at the interfaces between connected cells. Finally, the Godunov scheme ensures naturally the growth of physical entropy. Observe that, if in its first version the scheme was Lagrangian, it became very popular as an Eulerian method.

<footnote n="9">⁹The scheme is not conservative in total energy and unfortunately the conservation error does not vanish at convergence for classical CFL. However, the continuous in time scheme is conservative.</footnote>
Lagrangian numerical methods present a lot of advantages (especially in 1D). Since the convective terms are not discretized, the numerical dissipation is generally smaller compared to their Eulerian counterparts. Also, the absence of mass fluxes allows the natural treatment of multi-material flows: there is no need to consider mixed cell closure since each cell can be attached to a single material. Boundary conditions are generally easy to implement and Lagrangian methods are naturally adapted to free-surface flows. Finally, another important feature is that these numerical methods are by construction Galilean invariant.

For all these reasons, a lot of research has been devoted to the improvement of La-grangian methods for fluids. Starting from vNR and Godunov schemes, the most chal-lenging extension was probably their generalizations to higher dimensions of space.

Actually, in the case of *approximate* Godunov solvers, the first valid semi-Lagrangian multidimensional extension was defined by C. Mazeran and B. Després in [DM05; Maz07], thus almost 50 years after the original PhD work of Godunov. Before, in the early 90s, a promising attempt was implemented in the CAVEAT code [Add+90]. Similarly to the multidimensional Eulerian extensions of the Godunov scheme, it treats Riemann prob-lems at the faces of the mesh. However, the mesh displacement (or node velocity) is not defined naturally. It leads to compromises between accuracy and stability that are not driven by Numerical Analysis<fnref n="10" /> and then it is very complex to cure its flaws.

On the other hand multidimensional extensions of the vNR scheme have been devel-oped quite early. M. L. Wilkins published in 1964 a foundation paper [Wil64], where he extends the vNR scheme to elastic-plastic flows and proposes a treatment of 2D cylindri-cal geometries, improving W. B. Goad's approach [Goa60] with regard to the symmetry preservation of radial flows. In 1D, the conservation defects of the vNR scheme have been cured in the early 60s [TT61], and much later (in the early 90s) in dimension 2 [Bur90]. In the end of the 90s, a conservative staggered scheme for 2D-cylindrical geometries was proposed [Car+98]. It is not staggered in time, but uses a predictor-corrector approach to achieve second-order in time. Research around staggered Lagrangian schemes is still an active topic. For instance, a very high-order staggered Lagrangian scheme was published in [DKR12], or a conservative fully staggered scheme was proposed in [LCF16].

Despite the fact that one must tune artificial viscosity to get entropic numerical so-lutions<fnref n="11" />, staggered methods are very efficient. However, they are not intrinsically com-patible with ALE<fnref n="12" /> treatment. To be more precise, it is quite difficult to build an *indirect* ALE<fnref n="13" /> method that is conservative in mass, momentum and total energy [LS05].

On the other hand, since all conservative variables are defined at the same location, indirect ALE is naturally achieved when using finite-volume methods. Conservation is not an issue and one can show (see [BHS20]) in dimension *d* that first-order remapping ensures Maximum Principles on $\rho$, $(\rho u_i)_{1\le i\le d}$, $\rho E$, $(u_i)_{1\le i\le d}$, *E* and even for the specific internal energy $\epsilon$. This also makes it natural to combine AMR-like techniques with semi-Lagrangian methods, see Section 2.3.

Let us finally write Euler equations (2.5) and the physical entropy production (2.6) using semi-Lagrangian coordinates. Here, for two reasons, we use the integral form. Firstly it is a convenient formulation to derive finite-volume methods and secondly it is a conservative form. Let us finally remark that this form is only defined and equivalent to (2.5) and (2.6) in the case of smooth solutions. It reads $\forall t > 0$, and for any Lagrangian

<footnote n="10">¹⁰Due to its structure, numerical analysis of the CAVEAT scheme seems out of reach.</footnote>

<footnote n="11">¹¹In dimension greater than 1, defining properly artificial viscosity sensors is not that easy.</footnote>

<footnote n="12">¹²Arbitrary Lagrangian-Eulerian.</footnote>

<footnote n="13">¹³Indirect ALE consists in a splitting of the ALE formulation of Euler equations into three phases: a (semi-)Lagrangian phase, a grid rezoning phase and finally a remapping phase on the new grid.</footnote>
subdomain $\omega(t) \subset \Omega(t)^{14}$ (i.e. moving at the fluid velocity),

$$ \left\{ \begin{aligned} & \frac{d}{dt} \int_{\omega(t)} 1 = \int_{\partial\omega(t)} \mathbf{u} \cdot \mathbf{n}, && \left( = \int_{\omega(t)} \nabla \cdot \mathbf{u} \right) \\ & \frac{d}{dt} \int_{\omega(t)} \rho = 0, \\ & \frac{d}{dt} \int_{\omega(t)} \rho \mathbf{u} = - \int_{\partial\omega(t)} p \mathbf{n}, && \left( = - \int_{\omega(t)} \nabla p \right) \\ & \frac{d}{dt} \int_{\omega(t)} \rho E = - \int_{\partial\omega(t)} p \mathbf{u} \cdot \mathbf{n}, && \left( = - \int_{\omega(t)} \nabla \cdot p \mathbf{u} \right) \end{aligned} \right. \qquad (2.7) $$

and the physical entropy satisfies

$$ \frac{d}{dt} \int_{\omega(t)} \rho \eta \geq 0. \qquad (2.8) $$

The derivation of (2.7)–(2.8) is a straightforward application of the Reynolds Transport Theorem<fnref n="15" /> for smooth flows. For a rigorous definition of the (total-)Lagrangian formulation of Euler equations, one can refer to [Wag87; Wag96; Maz07]. The first equation in (2.7), the volume conservation, is a key ingredient with regard to the discretization on moving grids, see for instance [FGG01] in the case of ALE methods, it is generally designated as the Geometric Conservation Law (GCL).

## Extension of the acoustic solver to arbitrary dimension

As mentioned previously, this work was done with G. Carré, B. Després and E. Labou-rasse. In the very beginning, B. Després proposed to give a name to the scheme he devel-oped with C. Mazeran. He chose Glace for Godunov LAgrangian Conservative for the total Energy variable.

Actually, P.-H. Maire and B. Nkonga published a similar paper [MN08] that was sub-mitted simultaneously with our work [3]. The main differences between these two con-tributions are that, on the one hand they address the tri-dimensional extension of differ-ent schemes, respectively the Eucclhyd scheme [Mai+07] and the Glace scheme [DM05; Maz07]; and on the other hand the construction of the scheme in [MN08] relies on a more geometric point of view. This leads to an arbitrary decomposition of the cells into tetra-hedrons which can break the flow symmetry (this was addressed later in [GBM16]). In our derivation, the construction of the scheme relies on the so called “corner vectors” $C_{jr}$, which are defined here after. These vectors actually contain the geometrical informa-tion<fnref n="16" /> and play a crucial role in the semi-Lagrangian schemes construction. Also, in our paper, we briefly enlighten the natural ability of the cell-centered schemes to deal with ALE and AMR.

## The corner vectors $C_{jr}$

In order to fix ideas, let us consider a polyhedral cell $j$ and let us assume that the vol-ume $V_j$ of the cell $j$ is defined by the positions $x_s$ of its vertices $s \in \mathcal{R}_j$, that is $V_j =$

<footnote n="14">¹⁴Here we added a time dependency to the computational domain $\Omega$ since it is very natural using a semi-Lagrangian formulation, but it is not necessary.</footnote>

<footnote n="15">¹⁵It is better known as the Leibniz integral rule in Mathematics.</footnote>

<footnote n="16">¹⁶Changing the geometry, actually the definition of the sub-domains that define the cells, just consists in computing the associated set of $C_{jr}$ vectors using (2.9). For admissible choices of geometries, these sets of $C_{jr}$ vectors will share the same abstract properties (2.10)–(2.13).</footnote>
$V_j((\mathbf{x}_s)_{s \in \mathcal{R}_j})$. Then, the corner vectors $\mathbf{C}_{jr}$ are defined in dimension $d$ by

$$ \forall j \in \mathcal{J}, \forall r \in \mathcal{R}_j, \quad \mathbf{C}_{jr} := \nabla_{\mathbf{x}_r} V_j. \qquad (2.9) $$

Since these vectors of $\mathbb{R}^d$ actually measure the rate of change in volume of a cell $j$ according to the change in position of the vertex $r$, it has been used intensively to derive multi-dimensional semi-Lagrangian schemes (see for instance [Wha96; Car+98] and references therein).

Let us just recall the fundamental properties of $\mathbf{C}_{jr}$ vectors:

$$ \forall j \in \mathcal{J}, \quad V'_j = \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \cdot \mathbf{u}_r, \qquad (2.10) $$

$$ \forall j \in \mathcal{J}, \quad V_j = \frac{1}{d} \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \cdot \mathbf{x}_r, \qquad (2.11) $$

$$ \forall j \in \mathcal{J}, \quad \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} = \mathbf{0}, \qquad (2.12) $$

$$ \forall r \in \mathcal{R}, \quad \sum_{j \in \mathcal{J}_r} \mathbf{C}_{jr} = \mathbf{0}, \qquad (2.13) $$

where in (2.10), $\mathbf{u}_r = \frac{d}{dt}\mathbf{x}_r$ is the velocity of the vertex $r$. The relation (2.12) only occurs in planar geometries (i.e. not in cylindrical or spherical geometries). Finally (2.13) accounts for the local volume conservation.

### Scheme structure

The Glace scheme has the following structure

$$ \forall j \in \mathcal{J}, \quad \left| \begin{array}{l} m_j \tau'_j = \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \cdot \mathbf{u}_r, \\ m'_j = 0, \\ m_j \mathbf{u}'_j = - \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} p_{jr}, \\ m_j E'_j = - \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \cdot (p_{jr} \mathbf{u}_r), \end{array} \right. \qquad (2.14) $$

where $m_j = \rho_j(t)V_j(t)$, is the lagrangian mass of cell $j$, and $V_j(t)$ is the cell volume at time $t$. The mean cell velocity is $\mathbf{u}_j(t)$ and the mean specific total energy is denoted by $E_j(t)$. Finally the specific volume is $\tau_j(t) = \frac{V_j(t)}{m_j} = \frac{1}{\rho_j(t)}$. Setting $p_j = p(\rho_j, E_j - \frac{1}{2}\|\mathbf{u}_j\|^2)$ and $c_j$ the mean sound speed, the scheme is completely defined by using the acoustic Riemann invariants to link the mean cell values to the nodal values

$$ \forall j \in \mathcal{J}, \forall r \in \mathcal{R}_j, \quad p_{jr} - p_j + (\rho c)_j (\mathbf{u}_r - \mathbf{u}_j) \cdot \frac{\mathbf{C}_{jr}}{\|\mathbf{C}_{jr}\|} = 0, \qquad (2.15) $$

and setting the conservation constraint

$$ \forall r \in \mathcal{R}, \quad \sum_{j \in \mathcal{J}_r} \mathbf{C}_{jr} p_{jr} = \mathbf{0}. \qquad (2.16) $$

It is easy to show (see [3]) that the scheme ensures local conservation in volume, mass, momentum and total energy. Also one can check that the semi-discrete scheme (2.14)–(2.16) is entropy stable.
As it was observed in [Klu08; KD10; Mai11] the scheme can be written in a more general form

$$ \forall j \in \mathcal{J}, \quad \left| \begin{array}{l} m_j \tau'_j = \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \cdot \mathbf{u}_r, \\ m'_j = 0, \\ m_j \mathbf{u}'_j = - \sum_{r \in \mathcal{R}_j} \mathbf{F}_{jr}, \\ m_j E'_j = - \sum_{r \in \mathcal{R}_j} \mathbf{F}_{jr} \cdot \mathbf{u}_r, \end{array} \right. \qquad (2.17) $$

with

$$ \forall j \in \mathcal{J}, \forall r \in \mathcal{R}_j, \quad \mathbf{F}_{jr} = \mathbf{C}_{jr} p_j + A_{jr} (\mathbf{u}_j - \mathbf{u}_r), \qquad (2.18) $$

$$ \text{and} \quad \forall r \in \mathcal{R}, \quad \sum_{j \in \mathcal{J}_r} \mathbf{F}_{jr} = \mathbf{0}. \qquad (2.19) $$

Actually this form defines a family of schemes that are well defined, conservative and entropy stable, as soon as the matrices $A_{jr}$ are non-negative and such that $\forall r$, $\sum_{j \in \mathcal{J}_r} A_{jr}$ are invertible. In order to fix ideas we recall that choosing

$$ \forall j \in \mathcal{J}, \forall r \in \mathcal{R}_j, \quad A_{jr} = (\rho c)_j \frac{\mathbf{C}_{jr} \otimes \mathbf{C}_{jr}}{\|\mathbf{C}_{jr}\|}, \qquad (2.20) $$

defines the Glace scheme [DM05],[3], and that the Eucclhyd scheme is obtained by setting

$$ \forall j \in \mathcal{J}, \forall r \in \mathcal{R}_j, \quad A_{jr} = (\rho c)_j \sum_{l \in \mathcal{L}_j \cap \mathcal{L}_r} \frac{\mathbf{N}_{jlr} \otimes \mathbf{N}_{jlr}}{\|\mathbf{N}_{jlr}\|}, \qquad (2.21) $$

where $\mathbf{N}_{jlr}$ is the “facet normal”, see [Mai+07; MN08; GBM16] for details<fnref n="17" />.

In [3], we also cover the boundary conditions treatment which is not recalled here.

## Second-order extension

Let us finally discuss the very natural second-order extension of the scheme (2.17)–(2.19). Actually, since it is a finite-volume scheme, one can use all the classical recipes that have been investigated in the Eulerian framework. In other words, to achieve second-order in space it is enough to consider linear reconstructions of $\mathbf{u}_j$ and $p_j$. Thus assuming that in each cell $j$, the following linear reconstructions are defined

$$ \forall j \in \mathcal{J}, \quad \left| \begin{array}{l} \mathbf{x} \in j \mapsto \bar{\mathbf{u}}_j(\mathbf{x}), \\ \bar{\mathbf{x}} \in j \mapsto \bar{p}_j(\mathbf{x}), \end{array} \right. $$

then replacing the fluxes $\mathbf{F}_{jr}$ in (2.18) by

$$ \forall j \in \mathcal{J}, \forall r \in \mathcal{R}_j, \quad \mathbf{F}_{jr} = \mathbf{C}_{jr} \bar{p}_j(\mathbf{x}_r) + A_{jr}(\bar{\mathbf{u}}_j(\mathbf{x}_r) - \mathbf{u}_r), \qquad (2.22) $$

defines a second-order scheme. This is quite easy to figure out. Indeed, assuming that the reconstruction is exact, that is

$$ \forall j \in \mathcal{J}, \quad \left| \begin{array}{l} \bar{\mathbf{u}}_j = \mathbf{u}, \quad \text{and,} \\ \bar{p}_j = p, \end{array} \right. \quad \text{with } \forall \mathbf{x}, \quad \left| \begin{array}{l} \mathbf{u}(\mathbf{x}) = \mathbf{u}_0 + G_{\mathbf{x}} \mathbf{x}, \quad \text{and,} \\ \mathbf{p}(\mathbf{x}) = p_0 + \mathbf{g}_p \cdot \mathbf{x}, \end{array} \right. $$

<footnote n="17">¹⁷The expression is given here in dimension 2, since we use the edges sets $\mathcal{L}_j$ and $\mathcal{L}_r$. The 3D version which is algebraically similar, substituting edges sets by faces sets, is not given here for the sake of simplicity.</footnote>
where $p_0 \in \mathbb{R}$, $\mathbf{u}_0 \in \mathbb{R}^d$, $\mathbf{g}_p \in \mathbb{R}^d$ and $G_u \in \mathbb{R}^{d \times d}$, then injecting (2.22) into (2.19), one gets

$$ \forall r \in \mathcal{R}, \quad \underbrace{\sum_{j \in \mathcal{J}_r} \mathbf{C}_{jr} p(\mathbf{x}_r) + \sum_{j \in \mathcal{J}_r} A_{jr} (\mathbf{u}(\mathbf{x}_r) - \mathbf{u}_r)}_{=0 \text{ by (2.13)}} = \mathbf{0}, $$

which implies that $\mathbf{u}_r = \mathbf{u}(\mathbf{x}_r)$ and finally $\mathbf{F}_{jr} = \mathbf{C}_{jr} p(\mathbf{x}_r)$. Actually, since the implicit quadrature formulas defined by $\mathbf{C}_{jr}$ vectors are exact for affine functions (it is a trapezium formula), the scheme (2.17),(2.22),(2.19) provides an exact approximation for the three first equations of (2.5) (namely volume, mass and momentum conservation laws). It is not the case for the total energy conservation equation, but it is second-order accurate by means of trapezium quadrature formula.

In order to deal with shocked solutions, these reconstructions are limited. Let us just remark that when we published [3], the limitation of the velocity remained an issue and we needed to use specific limitation procedures. Since then, we benefited from the excellent work of G. Luttwak and J. Falcovitz, who propose in [LF11] a multi-dimensional limiter for vectors (VIP<fnref n="18" />). I did not specifically published on this treatment but contributed in its implementation. The important improvement it produces is that the scheme remains Galilean invariant when using VIP for the velocity.

I do not comment the second-order extension in time since it is straightforward, let me just remark that in [3], we use a one-step scheme that is not second-order in time but mimics the Lax-Wendroff procedure. In practice it is very accurate and cheap. The interested Reader is invited to check the details in [3].

### Very high-order extensions of semi-Lagrangian cell-centered schemes

When dealing with new schemes, a natural question is: how to define a very high-order version of the scheme? Quite early following the second-order implementation of the Glace scheme in the Troll code [Lef+18], we began to think about it with G. Carré and E. Labourasse.

Actually, one identifies easily the main difficulties. First of all, the volume conservation equation indicates clearly that reaching higher-order than 2 requires the use of bendable edges: cells cannot remain polygonal. The second difficulty is more classical: how to extend the flux integration to higher-order? Finally, how to reconstruct properly velocities and pressures, and how to limit them to treat discontinuous cases?

Meanwhile, J. Cheng and C.-W. Shu published in [CS07] a third-order semi-Lagrangian scheme. It is based on an ENO conservative reconstruction. However the grid velocity calculation is not really satisfactory in dimension 2 and requires the use of logical grids.

This is why, with G. Carré and E. Labourasse, we proposed at CEMRACS'08 a project to study two difficulties: numerical fluxes integration, and high-order reconstruction and limitation of cell pressures and velocities. In this project, we worked with K. P. Gostaf and A. V. Shapeev [4].

### High-order reconstruction and limitation

Since the known values are the conservative quantities, in order to achieve higher-order than 2, one must use reconstructions of $\rho$, $\rho\mathbf{u}$ and $\rho E$ in each cell $j$, namely $\bar{\rho}_j$, $(\rho\mathbf{u})_j$ and $(\rho E)_j$. Then one deduces the specific quantities as

$$ \forall j \in \mathcal{J}, \quad \bar{\mathbf{u}}_j = \frac{\overline{(\rho \mathbf{u})}_j}{\bar{\rho}_j}, \quad \bar{\epsilon}_j = \frac{\overline{(\rho E)}_j}{\bar{\rho}_j} - \frac{1}{2} \bar{\mathbf{u}}_j \cdot \bar{\mathbf{u}}_j \quad \text{and} \quad \bar{\rho}_j = p(\bar{\rho}_j, \bar{\epsilon}_j). $$

<footnote n="18">¹⁸Vector Image Polygon.</footnote>
The polynomial reconstruction itself is performed using the least-squares ENO method that has been proposed by C. F. Ollivier-Gooch in [Oll96a; Oll96b] and that we adapt in order to improve the condition number of the least-square matrices, following [Abg94]. In order to limit (or reduce the reconstruction degree), as described in [Oll96a; Oll96b], one can use the residual of the least-square approximation as an oscillation indicator. However, this strategy does not apply straightforwardly to the reconstruction of the conservative variables since $\bar{u}_j$ and $\bar{\epsilon}_j$ are rational fractions and $\bar{p}_j$ are nonlinear functions of them. In the paper we define some kind of limiters that aim at preserving the Galilean invariant property, but the limitation problem is not solved, especially in the case of the velocity. With that regard, a promising direction of improvement is to adapt VIP [LF11] to our method.

## High-order fluxes integration

The second aspect of the high-order extension that we study in [4] is the effects of the numerical integration of the fluxes. Actually, without any change in the structure of the scheme (2.17), (2.22), (2.19), one cannot expect higher-order than 2: the $C_{jr}$ vectors implicitly define trapezium formula. Our aim in [4] is simply to illustrate these effects on relevant test cases. Actually the Kidder test [Kid74] is well suited for this.

* Solving the 1D Kidder problem on an initially 2D aligned grid, the quadrature formulas associated to the $C_{jr}$ are exact, thus one must observe high-order when using reconstructions of conservative variables. See Figures 2.4a and 2.4b.

* Solving the 2D Kidder problem on an initially 2D Cartesian grid, since the velocity is a linear function, the grid remains rectilinear all along the calculation. In other words, their is no geometrical error in considering straight edges. The order loss, that should be observed, is only related to the numerical integration of the fluxes. See Figure 2.4c.

The expected results are displayed on Figure 2.4.

Actually, the treatment of high-order numerical integration of the fluxes is somehow related to the extension to high-order geometry of the cells. In other words, defining a proper way to bend edges can provide the tools to compute high-order fluxes. This is discussed below.

## A curvilinear extension of cell-centered schemes for semi-Lagrangian flows

Starting from dimension 2, the definition of very high-order semi-Lagrangian schemes requires the treatment of curvilinear cells (see Figure 2.5a for an illustration). The edges of the cells must bend, in a compatible way with the GCL, during the calculation. This problem has been addressed using isoparametric-$\Pi^k$ finite-elements in the case of staggered schemes [DKR12]. In the finite-volume context, F. Vilar proposed a total-Lagrangian extension of the Eucclhyd scheme to the third-order in his PhD Thesis [Vil12; VMA14].

Around 2009, I had a lot of “philosophical” discussions with P. Hoch on the one hand and E. Labourasse on the other hand, trying to answer the question: how to find a proper way to extend finite-volume schemes in semi-Lagrangian coordinates to curvilinear cells?

Actually, both of them proposed solutions to deal with curvilinear cells in quite similar ways. P. Hoch and his coworkers propose a way to deal with cells with boundaries parameterized as conical curves in [Hoc+11; Ber+12]. On the other hand, A. Claisse, B. Després and E. Labourasse manage to deal with circular edges when they define a treatment of exceptional points for semi-Lagrangian cell-centered schemes [Cla+12].
<figure><img src="image_9.png" /><figcaption>(a) 1D Kidder test. Convergence for reconstructions of various degrees of $\mathbf{u}$ and $\mathbf{p}$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) 1D Kidder test. Convergence for reconstructions of various degrees of conservative variables $\rho$, $\rho\mathbf{u}$ and $\rho E$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) 2D Kidder test. Even if higher-order reconstructions improve accuracy, the numerical integration implied by the $C_{jr}$ prevents high-order of convergence.</figcaption></figure>

**Figure 2.4:** Numerical illustration on the importance of using reconstructions of conservative variables and appropriate quadrature formulas.

<figure><img src="image_12.png" /><figcaption>(a) Deformation of a curvilinear cell: here edges are and remain parabolas.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Notations: a curvilinear cell $j$, its edges $l$ and vertices $r$. The outward unit normal $\mathbf{n}_{jl}$ is a function of $\mathbf{x} \in l$.</figcaption></figure>

**Figure 2.5:** Example of a curvilinear cell and of its motion.

These two approaches have in common the use of formula (2.9) in the case of curvilinear cells: the volume is expressed as a function of some control points $r$ and the deduced $C_{jr}$
satisfy all the properties (2.10)–(2.13) which allow to define a scheme.

Meanwhile, I proposed different approach to deal with curvilinear cells in [9]. The idea is to define an abstract framework to compute edge velocities. The obtained finite-volume scheme family should satisfy conservation and entropy stability properties. Ideally, in the case of polygonal cells one should retrieve known schemes (Glace [Maz07; DM05] and Eucclhyd [Mai+07]), so it should rely on the acoustic solver.

For the sake of simplicity, we restrict here to the 2D case<fnref n="19" /> ($d = 2$). Thus let us first give the general structure of the semi-discrete scheme family. It reads

$$ \forall j \in \mathcal{J}, \quad \left\{ \begin{aligned} \frac{d}{dt} \int_j 1 &= \sum_{l \in \mathcal{L}_j} \int_l \mathbf{u}^* \cdot \mathbf{n}_{jl}, \\ \frac{d}{dt} \int_j \rho &= 0, \\ \frac{d}{dt} \int_j \rho \mathbf{u} &= - \sum_{l \in \mathcal{L}_j} \int_l p_j^* \mathbf{n}_{jl}, \\ \frac{d}{dt} \int_j \rho E &= - \sum_{l \in \mathcal{L}_j} \int_l p_j^* \mathbf{u}^* \cdot \mathbf{n}_{jl}, \end{aligned} \right. \qquad (2.23) $$

where $\mathbf{u}^* : \mathcal{E} \to \mathbb{R}^d$ denotes the velocity of the edges of the mesh and $\forall j, p_j^* : \mathcal{E} \to \mathbb{R}$ is the pressure imposed by the cell $j$ to its edges<fnref n="20" />. One also introduces the unit normal to edge $l$, outgoing from cell $j$: $\forall j, \forall l \in \mathcal{L}_j, \mathbf{n}_{jl} : l \to \mathbb{R}^d$, see Figure 2.5b. At this point, we just assume that $|\mathcal{E}| < +\infty, \mathbf{u}^* \in L^2(\mathcal{E})^d, \forall j \in \mathcal{J}, p_j^* \in L^2(\mathcal{E})$ and that the mesh is regular enough to define $\mathbf{n}_{jl}$ a.e. on $l \in \mathcal{L}$. These are quite weak assumptions.

The abstract structure (2.23) does not impose neither the conservation of momentum nor the conservation of total energy, since one may have $p_j^* \neq p_k^*$ for $j \neq k \in \mathcal{J}_l$ for some edge $l$. Thus as it is the case for nodal solvers (that are mimicked here), we impose a conservation constraint<fnref n="21" />

$$ \forall \mathbf{v} \in L^2(\mathcal{E})^d, \quad \sum_{j \in \mathcal{J}} \sum_{l \in \mathcal{L}_j} \int_l p_j^* \mathbf{v} \cdot \mathbf{n}_{jl} = 0. \qquad (2.24) $$

Thus, it is easy to check that under the constraint (2.24), the abstract scheme (2.23) is conservative in volume, mass, momentum and total energy.

A simple and successful approach to build entropic semi-Lagrangian schemes relies on the use of the acoustic Riemann invariants $dp + \rho c du \cdot n = 0$. This gives a way to link $\mathbf{u}^*$ and $(p_j^*)_{j \in \mathcal{J}}: \forall j \in \mathcal{J}, \forall l \in \mathcal{L}_j$, for almost all $\mathbf{x}$ in $l$,

$$ p_j^*(\mathbf{x}) - p_j(\mathbf{x}) + (\rho c)_j(\mathbf{u}^*(\mathbf{x}) - \mathbf{u}_j(\mathbf{x})) \cdot \mathbf{n}_{jl}(\mathbf{x}) = 0. \qquad (2.25) $$

Actually, if we make the reasonable assumptions that

$$ \forall j \in \mathcal{J}, \forall l \in \mathcal{L}_j, \quad \left| \begin{array}{l} \mathbf{u}_j|_l \in L^2(l)^d, \text{ and} \\ p_j|_l \in L^2(l), \end{array} \right. $$

injecting (2.25) into (2.24) gives the following abstract velocity problem:

$$ \text{find } \mathbf{u}^* \in L^2(\mathcal{E})^d \text{ such that, } \quad \forall \mathbf{v} \in L^2(\mathcal{E})^d, \quad a(\mathbf{u}^*, \mathbf{v}) = \ell(\mathbf{v}), \qquad (2.26) $$

<footnote n="19">¹⁹The 3D extension is straightforward with regard to the abstract scheme structure, but is more technical when building concrete schemes.</footnote>

<footnote n="20">²⁰This implies that $\forall j \in \mathcal{J}, \forall l \in \mathcal{L} \setminus \mathcal{L}_j, p_j^*(\mathbf{x}) = 0$ a.e. on $l$.</footnote>

<footnote n="21">²¹Let us note that (2.24) just imposes weakly $p_j^* = p_k^*$ for $j, k \in \mathcal{J}_l$</footnote>
where the linear forms $a$ and $\ell$ are defined by

$$ 
\begin{array}{l}
a: L^2(\mathcal{E})^d \times L^2(\mathcal{E})^d \to \mathbb{R}, \\
(\mathbf{u}, \mathbf{v}) \mapsto \sum_{j \in \mathcal{J}} \sum_{l \in \mathcal{L}_j} \int_l \mathbf{v}^T A_{jl} \mathbf{u}, \text{ and} \\
\ell: L^2(\mathcal{E})^d \to \mathbb{R}, \\
\mathbf{v} \mapsto \sum_{j \in \mathcal{J}} \sum_{l \in \mathcal{L}_j} \int_l \mathbf{v}^T (p_j \mathbf{n}_{jl} + A_{jl} \mathbf{u}_j),
\end{array}
$$

with $A_{jl} := (\rho c)_j \mathbf{n}_{jl} \otimes \mathbf{n}_{jl}$.

It is easy to show that if $\mathbf{u}^*$ is a solution of (2.26), then the scheme defined by (2.23), (2.24) and (2.25) ensures the growth of physical entropy.

However, one should note that the problem (2.26) is not well-posed: the solution is not unique. The kernel of $a(\cdot, \cdot)$ is the space the tangential velocities to the edges $l$. So, if one only considers the vector space of normal velocities $N_{\mathcal{E}} := \{\mathbf{v} \in L^2(\mathcal{E})^d, \text{s.t. } (I - \mathbf{n} \otimes \mathbf{n}) \mathbf{v} = \mathbf{0}\}$, where the unit normal $\mathbf{n}: \mathcal{E} \to \mathbb{R}^d$ is defined almost everywhere on $\mathcal{E}$, then the problem

$$ \text{find } \mathbf{u} \in N_{\mathcal{E}}, \quad \text{s.t.} \quad \forall \mathbf{v} \in N_{\mathcal{E}}, \quad a(\mathbf{u}, \mathbf{v}) = \ell(\mathbf{v}), $$

admits a unique solution$^{22}$. At this point, the ill-posedness of (2.26) is not an issue in itself since the tangential component of $\mathbf{u}^*$ is not taken into account in (2.25) and (2.23).

In order to define a concrete scheme (i.e. that can be implemented) it remains to define two things, the way that the edges can transform and a vector space $W_h \subset L^2(\mathcal{E})^d$ of finite dimension$^{23}$. One can also use quadrature formulas in order to compute the integrals on the edges.

For instance, setting that the edges must remain straight lines, setting that the velocity field is piecewise linear by edge and continuous at nodes, and finally using a trapezium formula to compute the integrals, one retrieves the Euclyhyd scheme$^{24}$. This confirms that the abstract framework defined in [9] generalizes in some way known methods.

With P. Hoch and E. Labourasse, we have been working together on the derivation of very high-order schemes based on the abstract formulation (2.23)–(2.25). If we have not yet published on the subject, since we still need to find proper ways to deal with limitation, we presented the work in various conferences (see Figure 2.6). I. Marmajou has joined us recently and it has become again an active topic of research for us.

### A conservative slide-line method for compressible flows

In the end of 2011, we were discussing with E. Labourasse of internship subjects and we rapidly figured out that the abstract scheme (2.23), (2.24) and (2.25) was a good candidate to define a slide-line method for cell-centered schemes: for both semi-Lagrangian and indirect ALE methods. The sliding of fluids onto elastic-plastic solids is a topic of industrial interest and a lot of research papers have been devoted to it. On can refer to the review of N. G. Bourago and V. N. Kukudzhanov [BK05] that addresses a variety of different works (more than 600 references).

<footnote>$^{22}$It is a direct consequence of the Lax-Milgram Lemma, since $a(\cdot, \cdot)$ is a bilinear continuous and coercive form on $N_{\mathcal{E}}$ and $\ell(\cdot)$ is linear continuous on $N_{\mathcal{E}}$.</footnote>

<footnote>$^{23}$The choices must be made in order to get a well definite scheme, such that one can compute the edge velocity. It is also at this level that one introduces the continuity of the velocity, which is required to displace the mesh.</footnote>

<footnote>$^{24}$One can also retrieve the Glace[DM05; Maz07] scheme, but it is not that natural. The CAVEAT scheme is also derived straightfully (one still obtains normal velocities at the edges and similarly to [Add+90], the node velocities remain to be determined).</footnote>
<figure><img src="image_14.png" /><figcaption>Figure 2.6: Third-order preliminary results obtained with P. Hoch and E. Labourasse. We use a $\mathbb{P}^2$-isoparametric discretization to compute $\mathbf{u}^*$. The loss of coercivity related to this discretization and its treatment are not discussed in this document.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 2.7: The fluids in $\Omega_1$ and $\Omega_2$ are sliding on the interface $\Gamma$.</figcaption></figure>

Thus we proposed an internship subject on slide-line treatment. We decided to treat the case of a fluid-fluid sliding which is more challenging (as one can guess, the sliding of two fluids may be less stable). Moreover, adapting the method to the elastic-plastic case should be straightforward following [KD10; Mai+13].

After some preliminary results obtained during the internship, we decided to study more deeply the method we proposed. This work, done with S. Bertoluzza and E. Labourasse, has been published in [2].

So, let us consider two connected open sets $\Omega_1$ and $\Omega_2$ (s.t. $\Omega_1 \cap \Omega_2 = \emptyset$). These define two compressible and inviscid fluid domains, modeled by the Euler equations. In other words, one has

$$ \forall i \in \{1, 2\}, \text{ a.e. in } \Omega_i, \quad \left. \begin{array}{l} \partial_t \rho_i + \nabla \cdot (\rho_i \mathbf{u}_i) = 0, \\ \partial_t (\rho_i \mathbf{u}_i) + \nabla \cdot (\rho_i \mathbf{u}_i \otimes \mathbf{u}_i) + \nabla p_i = \mathbf{0}, \\ \partial_t (\rho_i E_i) + \nabla \cdot (\rho_i E_i \mathbf{u}_i) + \nabla \cdot (p_i \mathbf{u}_i) = 0, \end{array} \right\} (2.27) $$

where the notations are the same as in (2.5) and the index $i \in \{1, 2\}$ denotes the fluid $i$. Let us now define $\Gamma := \partial\Omega_1 \cap \partial\Omega_2$, the sliding interface that is assumed regular. Since we consider a perfect sliding, one has

$$ (\mathbf{u}_1 - \mathbf{u}_2) \cdot \mathbf{n}_{\Gamma} = 0, \text{ a.e. on } \Gamma, \qquad (2.28) $$
where $\mathbf{n}_\Gamma$ designates the unit normal to $\Gamma$, see Figure 2.7. In its weak form

$$ \forall \mu \in L^2(\Gamma), \quad \int_{\Gamma} \mu(\mathbf{u}_1 - \mathbf{u}_2) \cdot \mathbf{n}_{\Gamma} = 0, \qquad (2.29) $$

it is a classical formulation that is clearly linked to [BMP94; Ben99; BM92], which studied numerical methods to ensure weak continuity of quantities for non-matching grids in the context of domain decomposition methods. Since the grid velocity problem (2.26) defined in the previous paragraph is the variational formulation of a minimization problem, it seemed a natural starting point to build a slide line method.

Thus, our approach consists in using the abstract semi-Lagrangian framework called in the previous paragraph to define an abstract sliding method. One uses the same structure (2.23) for each fluid. So, we introduce the forms

$$ \forall i \in \{1,2\}, \quad a_i : L^2(\mathcal{E}_i)^d \times L^2(\mathcal{E}_i)^d \to \mathbb{R}, \quad (\mathbf{u}, \mathbf{v}) \mapsto \sum_{j \in \mathcal{J}^i} \sum_{l \in \mathcal{L}_j} \int_l \mathbf{v}^T A_{jl} \mathbf{u}, \text{ and} $$

$$ \forall i \in \{1,2\}, \quad \ell_i : L^2(\mathcal{E}_i)^d \to \mathbb{R}, \quad \mathbf{v} \mapsto \sum_{j \in \mathcal{J}^i} \sum_{l \in \mathcal{L}_j} \int_l \mathbf{v}^T (p_j \mathbf{n}_{jl} + A_{jl} \mathbf{u}_j), $$

with $A_{jl} := (\rho c)_j \mathbf{n}_{jl} \otimes \mathbf{n}_{jl}$, and the edge pressures $p_j^*$ are defined, following (2.25), by

$$ \forall i \in \{1,2\}, \forall j \in \mathcal{J}^i, \quad p_j^*(\mathbf{x}) - p_j(\mathbf{x}) + (\rho c)_j (\mathbf{u}^*(\mathbf{x}) - \mathbf{u}_j(\mathbf{x})) \cdot \mathbf{n}_{jl}(\mathbf{x}) = 0, \text{ for almost all } \mathbf{x} \text{ in } l. $$

Introducing the space of perfectly sliding grid velocities

$$ \mathcal{C} := \left\{ (\mathbf{v}_1, \mathbf{v}_2) \in L^2(\mathcal{E}_1)^d \times L^2(\mathcal{E}_2)^d \text{ s.t. } \forall \mu \in L^2(\Gamma), \int_{\Gamma} \mu(\mathbf{v}_1 - \mathbf{v}_2) \cdot \mathbf{n}_{\Gamma} = 0 \right\}, $$

the abstract sliding velocity problem can be written as, find $(\mathbf{u}_1^*, \mathbf{u}_2^*) \in \mathcal{C}$ such that

$$ \forall (\mathbf{v}_1, \mathbf{v}_2) \in \mathcal{C}, \quad \sum_{i=1}^{2} a_i(\mathbf{u}_i^*, \mathbf{v}_i) = \sum_{i=1}^{2} \ell_i(\mathbf{v}_i). \qquad (2.30) $$

Actually, if $(\mathbf{u}_1^*, \mathbf{u}_2^*) \in \mathcal{C}$ is a solution of (2.30), it is easy to prove that the abstract sliding scheme is conservative in volume, mass, momentum and total energy. Moreover it is entropy stable for piecewise constant data. However, it is the case because the geometry of the slide line is well defined as it corresponds exactly to $\partial\Omega_1 \cap \partial\Omega_2$. As it is depicted on Figure 2.8, using polygonal grids, generally leads to non-matching interfaces, thus the definition of $\Gamma$ itself is part of the problem. So, assuming in the case of polygonal grids, that a mesh of $\Gamma$ has been provided (see Figure 2.8b), we define a concrete $\mathbb{P}^1$ scheme structure as

$$ \forall i \in \{1,2\}, \forall j \in \mathcal{J}^i, \left. \begin{array}{l} V'_j = \sum_{l \in \mathcal{L}_j} \int_l \mathbf{u}_{ih}^* \cdot \mathbf{n}_{jl}, \\ M'_j = 0, \\ M_j \mathbf{u}'_j = - \sum_{l \in \mathcal{L}_j \setminus \Gamma_i} \int_l p_{ih,j}^* \mathbf{n}_{jl} - \sum_{l \in \mathcal{L}_j \cap \Gamma_i} \int_l p_{ih,j}^* \mathbf{n}_{\Gamma,i}, \\ M_j E'_j = - \sum_{l \in \mathcal{L}_j \setminus \Gamma_i} \int_l p_{ih,j}^* \mathbf{u}_{ih}^* \cdot \mathbf{n}_{jl} - \sum_{l \in \mathcal{L}_j \cap \Gamma_i} \int_l p_{ih,j}^* \mathbf{u}_{ih}^* \cdot \mathbf{n}_{\Gamma,i} \end{array} \right\} \qquad (2.31) $$
<figure><img src="image_16.png" /><figcaption>(a) Example of curved meshes which match perfectly the sliding interface.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Example of polygonal meshes. The slide line itself is an unknown of the problem.</figcaption></figure>

Figure 2.8: Comparison of the ideal mesh configuration in the case of the abstract scheme and the practical configuration occurring with polygonal meshes.

where $\Gamma_i \subset \mathcal{L}^i$ is the subset of edges that are sliding on $\Gamma$. One uses a bijective mapping $T_{i,\Gamma}$ of the edges of $\Gamma_i$ to $\Gamma$ in order to evaluate the unit normal $\mathbf{n}_{\Gamma,i}$ to $\Gamma$, outgoing from $\Omega_i$. Observe that in general it does not coincide with the normal to the considered edge. To ensure conservation, one then imposes the following constraint on all $p_{ih,j}^*$'s

$$ \forall (\mathbf{v}_{1h}, \mathbf{v}_{2h}) \in \mathbb{P}^1(\mathcal{E}_1)^d \times \mathbb{P}^1(\mathcal{E}_2)^d, \text{ s.t. } (\mathbf{v}_{1h} - \mathbf{v}_{2h}) \cdot \mathbf{n}_{\Gamma} \text{ on } \Gamma, \\ \sum_{i=1}^{2} \sum_{j \in \mathcal{J}^i} \sum_{l \in \mathcal{L}_j \setminus \Gamma_i} \int_l p_{ih,j}^* \mathbf{v}_{ih} \cdot \mathbf{n}_{jl} + \sum_{i=1}^{2} \sum_{j \in \mathcal{J}^i} \sum_{l \in \mathcal{L}_j \cap \Gamma_i} \int_l p_{ih,j}^* \mathbf{v}_{ih} \cdot \mathbf{n}_{\Gamma,i} = 0. \quad (2.32) $$

**Remark 1.** Let us emphasize that in the previous statement, the term $(\mathbf{v}_{h1} - \mathbf{v}_{h2}) \cdot \mathbf{n}_{\Gamma}$ is written in an improper way since the functions are not defined on the same domain: $\mathbf{v}_{hi}$ are defined on $\mathcal{E}_i$, so they are defined on $\Gamma_i$ and $\mathbf{n}_{\Gamma}$ is defined on $\Gamma$. Actually, the correct way to write these terms is

$$ (\mathbf{v}_{h1} \circ T_{\Gamma,1} - \mathbf{v}_{h2} \circ T_{\Gamma,2}) \cdot \mathbf{n}_{\Gamma} \text{ on } \Gamma, $$

where $\forall i \in \{1,2\}, T_{\Gamma,i} := T_{i,\Gamma}^{-1}$ are the chosen bijective mappings of $\Gamma$ to $\Gamma_i$.

In the following we will use this improper notation for the sake of simplicity.

So, the scheme structure (2.31)–(2.32) defines a family of conservative schemes in mass, momentum and total energy. It remains though to compute the fluxes, and again we use the acoustic Riemann invariants (here in weak form for commodity):

$$ \forall i \in \{1,2\}, \forall j \in \mathcal{J}^i, \forall \mathbf{v}_{ih} \in \mathbb{P}^1(\mathcal{E}_i)^d, \\ \int_l \mathbf{v}_{ih}^T p_{ih,j}^* \mathbf{n}_{jl} - \int_l \mathbf{v}_{ih}^T p_j \mathbf{n}_{jl} + \int_l \mathbf{v}_{ih}^T A_{jl} \mathbf{u}_{ih}^* - \int_l \mathbf{v}_{ih}^T A_{jl} \mathbf{u}_j = 0, \quad (2.33) $$

$$ \text{and} \quad \forall l \in \mathcal{L}_j \cap \Gamma_i, \quad \int_l \mathbf{v}_{ih}^T p_{ih,j}^* \mathbf{n}_{\Gamma,i} - \int_l \mathbf{v}_{ih}^T p_j \mathbf{n}_{\Gamma,i} + \int_l \mathbf{v}_{ih}^T A_{\Gamma,j} \mathbf{u}_{ih}^* - \int_l \mathbf{v}_{ih}^T A_{\Gamma,j} \mathbf{u}_j = 0, \quad (2.34) $$

where $A_{jl} := (\rho c)_j \mathbf{n}_{jl} \otimes \mathbf{n}_{jl}$ and $A_{\Gamma,j} := (\rho c)_j \mathbf{n}_{\Gamma} \otimes \mathbf{n}_{\Gamma}$.

Now injecting (2.33) and (2.34) into (2.32), one gets the following forms

$$ \forall i \in \{1,2\}, \quad a_i^h : \mathbb{P}^1(\mathcal{E}_i)^d \times \mathbb{P}^1(\mathcal{E}_i)^d \to \mathbb{R}, \\ (\mathbf{u}_{ih}, \mathbf{v}_{ih}) \mapsto \sum_{j \in \mathcal{J}^i} \sum_{l \in \mathcal{L}_j \setminus \Gamma_i} \int_l \mathbf{v}_{ih}^T A_{jl} \mathbf{u}_{ih} + \sum_{j \in \mathcal{J}^i} \sum_{l \in \mathcal{L}_j \cap \Gamma_i} \int_l \mathbf{v}_{ih}^T A_{\Gamma,j} \mathbf{u}_{ih}, $$
and

$$ \forall i \in \{1,2\}, \quad \ell_i^h : \mathbb{P}^1(\mathcal{E}_i)^d \to \mathbb{R}, \\ \mathbf{v}_{ih} \mapsto \sum_{j \in \mathcal{J}^i} \sum_{l \in \mathcal{L}_j \setminus \Gamma_i} \int_l \mathbf{v}_{ih}^T (p_j \mathbf{n}_{jl} + A_{jl} \mathbf{u}_j) \\ + \sum_{j \in \mathcal{J}^i} \sum_{l \in \mathcal{L}_j \setminus \Gamma_i} \int_l \mathbf{v}_{ih}^T (p_j \mathbf{n}_{\Gamma,i} + A_{\Gamma,j} \mathbf{u}_j). $$

Finally, in order to use Lagrange multipliers to take into account the perfect sliding constraint, we introduce the forms

$$ \forall i \in \{1,2\}, \quad b_i^h : \mathbb{P}^1(\mathcal{E}_i)^d \times \mathbb{P}^0(\Gamma) \to \mathbb{R}, \\ (\mathbf{v}_{ih}, \mu_h) \mapsto \sum_{l \in \Gamma_i} \int_l \mathbf{v}_{ih} \cdot \mathbf{n}_i^\Gamma \, \mu_h. $$

Introducing $V_h := \mathbb{P}^1(\mathcal{E}_1)^d \times \mathbb{P}^1(\mathcal{E}_2)^d \times \mathbb{P}^0(\Gamma)$, the discrete velocity problem reads$^{25}$,

find $(\mathbf{u}_{1h}^*, \mathbf{u}_{2h}^*, \lambda_h) \in V_h$ such that

$$ \forall (\mathbf{v}_{1h}, \mathbf{v}_{2h}, \mu_h) \in V_h, \quad \left| \begin{array}{l} \sum_{i=1}^{2} a_i^h (\mathbf{u}_{ih}^*, \mathbf{v}_{ih}) + \sum_{i=1}^{2} b_i^h (\lambda_h, \mathbf{v}_{ih}) = \sum_{i=1}^{2} \ell_i^h (\mathbf{v}_{ih}), \\ \sum_{i=1}^{2} b_i^h (\mathbf{u}_{ih}^*, \mu_h) = 0. \end{array} \right. \qquad (2.35) $$

In [2], we show that, as soon as the mesh that defines $\Gamma$ is locally coarser than one of the two meshes $\Gamma_1$ or $\Gamma_2$ then, the saddle point problem admits a unique solution that satisfies

$$ \sum_{i=1}^{2} \| \mathbf{u}_{ih}^* \|_{0, \mathcal{E}_i} + \| \lambda_h \|_{0, \Gamma} \lesssim \sum_{i=1}^{2} \| \ell_i^h \|_{L^2(\mathcal{E}_i)}. \qquad (2.36) $$

The method is proved to be conservative in mass, momentum and total energy. The volume conservation is lost along the slide line. This is not really a surprise since the geometry of the interface is not clearly defined. It implies that the entropy production in the cells along the interface is positive if and only if

$$ p_j \sum_{l \in \mathcal{L}_j \cap \Gamma_i} \int_l (\mathbf{u}_{ih}^* - \mathbf{u}_j) \cdot (\mathbf{n}_{\Gamma,i} - \mathbf{n}_{jl}) \\ \leq \sum_{l \in \mathcal{L}_j \cap \Gamma_i} \int_l (\mathbf{u}_j - \mathbf{u}_{ih}^*)^T A_{jl} (\mathbf{u}_j - \mathbf{u}_{ih}^*) + \sum_{l \in \mathcal{L}_j \cap \Gamma_i} \int_l (\mathbf{u}_j - \mathbf{u}_{ih}^*)^T A_{\Gamma,i} (\mathbf{u}_j - \mathbf{u}_{ih}^*). $$

This defect in entropy production (when the inequality is not satisfied) is localized to the sliding boundary thus should tend to zero with the mesh size, allowing convergence to entropic solutions. In order to increase the stability of the whole scheme, we added subzone entropy [DL12] at the vicinity of slide lines. This is not described here$^{26}$.

Let us emphasize that, this scheme and the scheme proposed in [CDL14] are the only approaches in our knowledge that ensure the conservation of mass, momentum and total energy for sliding in semi-Lagrangian coordinates. In other recent works such

<footnote>$^{25}$Observe that in [2], we choose a $\mathbb{P}^0$ discretization for the Lagrange multipliers living on $\Gamma$ but other choices may be interesting.</footnote>

<footnote>$^{26}$One should however note that increasing stability at the slide line is a classical recipe. See for instance [Car09; Kuc+12].</footnote>
as [Kuc+12; Mor+13] the defect of total energy conservation is often used as an *a posteriori* indicator to assess the quality of the numerical solution.

This section is concluded by a numerical example. We use the test proposed by E. J. Caramana in [Car09]. The computational domain is $\Omega = ]0, 1[ \times ]0, \frac{1}{2}[$. Initially, the fluid state is defined by

$$ \mathbf{u} = 0, \quad \rho = \begin{cases} 1 & \text{if } y < \frac{1}{4}, \\ 10 & \text{elsewhere,} \end{cases} \quad p = \begin{cases} 20 & \text{in } ]0.95, 1[ \times ]0, \frac{1}{4}[, \\ \frac{2}{3} \times 10^{-8} & \text{elsewhere.} \end{cases} $$

The fluid follows a perfect gas law with an adiabatic constant $\gamma = \frac{5}{3}$. Symmetric boundary conditions are imposed on the whole $\partial\Omega$.

This initial data generates a shock wave due to the high pressure zone that travels in the domain. Combined with the density jump at the interface, this leads to sliding along this discontinuity line.

<figure><img src="image_18.png" /><figcaption>(a) Initial 100 × 50 mesh.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) Pressure field — Time $t = 0.3$. Mesh 200 × 100. Lagrangian simulation with stabilization.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(c) Pressure field — Time $t = 0.285$. Mesh 200 × 100. Lagrangian simulation without any stabilization.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(d) Pressure field — Time $t = 0.3$. Mesh 200 × 100. ALE simulation without stabilization.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(e) Density in each domain — Time $t = 0.3$. Mesh 200 × 100. Lagrangian simulation with stabilization.</figcaption></figure>

<caption>Figure 2.9: Caramana test case</caption>

The results are presented on Figure 2.9 where we compare various strategies: with or without subzone entropy stabilization and where we show the efficiency of a simple
ALE strategy which consists in rezoning the grids so that the interfaces are matching at the beginning of each time step.

### Mixture of two compressible flows coupled with friction

In 2014, with E. Labourasse, we proposed to treat the mixture of compressible flows using a multidimensional asymptotic preserving scheme as an internship subject. The subject itself is an extension to the PhD work of C. Enaux [Ena07]<fnref n="27" /> on unstructured grids in dimension 2. Following this work, we chose the Scannapieco-Cheng model [SC02].

We published this work [14] with E. Labourasse and G. Morel. We consider a simplified version of the Scannapieco-Cheng model for the mixture of two fluids $f_1$ and $f_2$. It reads in semi-Lagrangian coordinates

$$ \forall \alpha \in \{f_1, f_2\}, \quad \begin{cases} \rho^\alpha D_t^\alpha \tau^\alpha = \nabla \cdot \mathbf{u}^\alpha, \\ \rho^\alpha D_t^\alpha \mathbf{u}^\alpha = -\nabla p^\alpha - \nu \rho \delta \mathbf{u}^\alpha, \\ \rho^\alpha D_t^\alpha E^\alpha = -\nabla \cdot (p^\alpha \mathbf{u}^\alpha) - \nu \rho \delta \mathbf{u}^\alpha \cdot \bar{\mathbf{u}}, \end{cases} \quad (2.37) $$

where $\rho^\alpha$, $\mathbf{u}^\alpha$ and $E^\alpha$ respectively denote the mass density, the velocity and the total energy density of fluid $\alpha$. Also, $\tau^\alpha = \frac{1}{\rho^\alpha}$ denotes the specific volume. The pressure $p^\alpha$ satisfies the equation of state $p^\alpha := p^\alpha(\rho^\alpha, e^\alpha)$, where $e^\alpha$, the internal energy density, is defined by $e^\alpha := E^\alpha - \frac{1}{2}\|\mathbf{u}^\alpha\|^2$. The total density $\rho$ and the mean velocity $\bar{\mathbf{u}}$ are defined as $\rho := \rho^\alpha + \rho^\beta$ and $\rho\bar{\mathbf{u}} := \rho^\alpha \mathbf{u}^\alpha + \rho^\beta \mathbf{u}^\beta$. The term $\delta\mathbf{u}^\alpha$ is the velocity difference, the $\delta(\cdot)^\alpha$ operator being defined by $\delta\phi^\alpha = -\delta\phi^\beta = \phi^\alpha - \phi^\beta$. Finally, $\nu$ is the friction parameter. Also, remark that the Lagrangian derivative $D_t^\alpha := \partial_t + \mathbf{u}^\alpha \cdot \nabla$, is obviously not the same for each fluid.

The model (2.37) is conservative in volume, mass, in the sum of momenta and in the sum of total energies of the two fluids. The Gibbs formulas written for each fluid yield the following entropy inequalities

$$ \forall \alpha, \beta \in \{f_1, f_2\}, \text{ s.t. } \alpha \neq \beta, \quad T^\alpha D_t^\alpha \eta^\alpha \geq \nu \frac{\tau^\alpha}{\tau^\beta} \delta \mathbf{u}^\alpha \cdot \delta \mathbf{u}^\alpha \geq 0. \quad (2.38) $$

As it as been established in [Ena07], by means of Hilbert expansions, in the limit $\nu \to \infty$, the model (2.37) behaves like

$$ \rho D_t \mathbf{u} = -\nabla (p^\alpha + p^\beta), \quad (2.39) $$

while, for each fluid $\alpha \in \{f_1, f_2\}$, $\beta$ denoting the other one, one has

$$ \begin{aligned} \rho^\alpha D_t^\alpha \tau^\alpha &= \nabla \cdot \mathbf{u}, \\ \rho^\alpha D_t^\alpha E^\alpha &= -\frac{\rho^\alpha}{\rho} \mathbf{u} \cdot \nabla (p^\alpha + p^\beta) - p^\alpha \nabla \cdot \mathbf{u}, \end{aligned} \quad (2.40) $$

where $\mathbf{u}$ is the same velocity for both fluids, and thus the Lagrangian derivative is also the same.

The idea of our work [14] is to provide an indirect ALE method that consists in two phases (see Figure 2.10). In the first phase, starting from a common grid, each fluid evolves in a semi-Lagrangian way (the grid displacement is different for each fluid) and then the fluids are remapped onto a new common grid, allowing to proceed with the next time step.

If the remap phase is quite standard in the finite-volume context, the semi-Lagrangian cell-centered scheme we propose is the novelty.

<footnote n="27">²⁷In this work, dimensions greater than 1 are achieved by means of directional splitting.</footnote>
<figure><img src="image_23.png" /><figcaption>Figure 2.10: **Left:** at time $t = t^n$, both fluids share the same mesh. **Middle:** at the end of the Lagrangian phase, one gets two different meshes, one for each fluid. **Right:** meshes are displaced so that they coincide. Solution is remapped and a new timestep can be performed.</figcaption></figure>

Let us introduce some notations. Let $\rho_r^\alpha := \frac{1}{\#\mathcal{J}_r} \sum_{j \in \mathcal{J}_r} \rho_j^\alpha$ and $\rho_r := \rho_r^\alpha + \rho_r^\beta$. Also, we set the following mean velocities $\bar{\mathbf{u}}_r := \frac{\rho_r^\alpha \mathbf{u}_j^\alpha + \rho_r^\beta \mathbf{u}_j^\beta}{\rho_r^\alpha + \rho_r^\beta}$ and $\bar{\mathbf{u}}_{jr} := \frac{\rho_r^\alpha \mathbf{u}_j^\alpha + \rho_r^\beta \mathbf{u}_j^\beta}{\rho_r^\alpha + \rho_r^\beta}$. The $A_{jr}^\alpha$ matrices are defined by (2.20) or (2.21), and finally the $B_{jr}$ matrices are symmetric positive definite and satisfy $\sum_{r \in \mathcal{R}_j} B_{jr} = V_j I$. In the paper, we discuss and test various choices of $B_{jr}$ that are not recalled here.

Let $\alpha \in \{f_1, f_2\}$ denote one of the two fluids and $\beta$ the other one, we define the scheme,

$$ \left. \begin{array}{l} m_j^\alpha d_t \tau_j^\alpha = \sum_r \mathbf{C}_{jr} \cdot \mathbf{u}_r^\alpha, \\ d_t m_j^\alpha = 0, \\ m_j^\alpha d_t \mathbf{u}_j^\alpha = -\sum_r \mathbf{F}_{jr}^\alpha - \sum_r v \rho_r B_{jr} \delta \mathbf{u}_j^\alpha, \\ m_j^\alpha d_t E_j^\alpha = -\sum_r \mathbf{F}_{jr}^\alpha \cdot \mathbf{u}_r^\alpha - \sum_r v \rho_r \bar{\mathbf{u}}_r^T B_{jr} \delta \mathbf{u}_r^\alpha + \sum_r v \rho_r \bar{\mathbf{u}}_{jr}^T B_{jr} (\delta \mathbf{u}_r^\alpha - \delta \mathbf{u}_j^\alpha), \end{array} \right\} \forall j \in \mathcal{J}, \quad (2.41) $$

where the fluxes are given by

$$ \forall j \in \mathcal{J}, \forall r \in \mathcal{R}_j, \quad \left| \begin{array}{l} \mathbf{F}_{jr}^\alpha = \mathbf{C}_{jr} p_j^\alpha + A_{jr}^\alpha (\mathbf{u}_j^\alpha - \mathbf{u}_r^\alpha) - v \rho_r B_{jr} \delta \mathbf{u}_r^\alpha, \quad \text{and} \\ \sum_j \mathbf{F}_{jr}^\alpha = \mathbf{0}. \end{array} \right. \quad (2.42) $$

Performing an Hilbert expansion, in the limit $v \to +\infty$, the semi-discrete scheme behaves like the following one: $\forall \alpha, \beta \in \{f_1, f_2\}$, with $\alpha \neq \beta$, $\forall j \in \mathcal{J}$,

$$ \left| \begin{array}{l} (m_j^\alpha + m_j^\beta) d_t \mathbf{u}_j = -\sum_r \mathbf{F}_{jr}^\alpha - \sum_r \mathbf{F}_{jr}^\beta, \\ d_t V_j = m_j^\alpha d_t \tau_j^\alpha = \sum_r \mathbf{C}_{jr} \cdot \mathbf{u}_r, \\ d_t m_j^\alpha = 0, \\ m_j^\alpha d_t E_j^\alpha = -\sum_r \mathbf{C}_{jr} p_j^\alpha \cdot \mathbf{u}_r + \sum_r \mathbf{u}_r^T A_{jr}^\alpha (\mathbf{u}_r - \mathbf{u}_j) - \frac{\rho_j^\alpha \rho_j^\beta}{\rho_j} \sum_r \mathbf{u}_j^T \delta \left( \frac{A_{jr}}{\rho_j} \right)^\alpha (\mathbf{u}_r - \mathbf{u}_j), \end{array} \right. \quad (2.43) $$

where $\mathbf{u}_j = \mathbf{u}_j^\alpha = \mathbf{u}_j^\beta$, and where the nodal velocities $\mathbf{u}_r = \mathbf{u}_r^\alpha = \mathbf{u}_r^\beta$ satisfy

$$ \forall j \in \mathcal{J}, \forall r \in \mathcal{R}_j, \quad \left| \begin{array}{l} \mathbf{F}_{jr}^\alpha + \mathbf{F}_{jr}^\beta = \mathbf{C}_{jr} (p_j^\alpha + p_j^\beta) + (A_{jr}^\alpha + A_{jr}^\beta) (\mathbf{u}_j - \mathbf{u}_r), \\ \sum_j \mathbf{F}_{jr}^\alpha = \mathbf{0}. \end{array} \right. \quad (2.44) $$

In [14], we show that the asymptotic scheme is consistent with the asymptotic model. We analyze the fully discrete scheme associated to (2.41)–(2.42) and prove that it is conservative in volume, mass and in the sum of momenta and of total energies. We also
prove that it is entropy stable (the entropy production is lower bounded independently of the value of $v$) if the terms $\bar{u}_{jr}$ and $\delta \mathbf{u}_j^\alpha$ are discretized implicitly. The scheme (2.41)–(2.42) is thus asymptotic preserving$^{28}$.

We conclude this section, by illustrating the effects of the value of $v$ on a Rayleigh-Taylor like instability, see Figure 2.11.

<figure><img src="image_24.png" /><figcaption>Figure 2.11: 80 × 224 mesh. Time $t = 0.7$. Mass fraction of fluid $\alpha$. Influence of the friction parameter. **Left:** $v = 100$. **Middle:** $v = 1000$. **Right:** $v = 10^6$.</figcaption></figure>

### Implicit semi-Lagrangian schemes for compressible flows

In 2018, I was involved in a new topic: fluid-structure interaction in the special case of shocks dynamics. Actually, this was already a field of research in my team in the Eulerian framework (in the HERA code [Jou05]). In this Cartesian AMR code, the treatment is monolithic: fluids are treated as a particular case of elastic-plastic solids. This presents the advantage of simplifying the numerics since the same scheme is more or less used in the whole computational domain. The main defect of this strategy is that it can be very expensive when dealing with thin structures since they may require very fine meshes to capture correctly the geometry. In the case of explicit schemes these very small cells lead to very small time steps too, increasing even more the cost of the calculation.

Thus, various strategies are studied in my team to improve the efficiency of this kind of calculations.

Concerning semi-Lagrangian or indirect ALE methods, starting from the seminal work of M. L. Wilkins [Wil64] various extensions to elastic-plastic flows have been studied in the past years. More recently, in the finite-volume framework, one can refer to [Klu08; KD10; Mai+13] for elastic-plastic extensions. So, implementing an ALE mono-

<footnote>$^{28}$Actually, we also proved that the scheme is consistent with a diffusion equation of the concentration that is obtained by means of first-order Hilbert expansion.</footnote>
lithic solver is not a difficult task in this context. However, it suffers the same kind of
flaws with regard to the computational cost when dealing with thin structures.

Discussing with B. Després, we concluded that if we could define a semi-Lagrangian implicit solver<fnref n="29" /> to treat the thin structures and their vicinity, it could reduce the cost and allow greater time steps. So we wrote a PhD Thesis subject on this topic. This is the PhD Thesis subject of A. Plessier [Ple23], I co-supervise this work with B. Després.

As a starting point, we decided to focus on the definition of an implicit finite-volume scheme to treat compressible flows in semi-Lagrangian coordinates. To do so, we took inspiration in the work of C. Chalons, F. Coquel and C. Marmignon [CCM10] where they define an implicit scheme to solve Euler equations using a predictor-corrector approach. The prediction step consists in the resolution of the isentropic Euler equations to overcome the non-linearities contained in the total energy conservation equation. The correction step allows to retrieve the total energy conservation. In [CCM10], the Authors use a relaxation scheme to get rid of the non-linearities induced by the equation of state in the momentum equation.

In dimension 1, we proposed a method (published in [21]) that takes advantage of the semi-Lagrangian formulation to deal with the isentropic prediction step. It benefits from the fact that one does not need to handle the non-linearities due to the transport terms. Actually, the method can be extended to dimension 2 and the notations in this case are better (the rearrangement of some of the terms gives a simpler understanding of the method and provides some generalizations). Thus in this manuscript, I adopt the 2D notations.

The prediction-step is

$$
\forall j \in \mathcal{J}, \quad
\left\{
\begin{aligned}
\bar{\tau}_j &= \tau_j^n + \frac{\Delta t}{m_j} \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \cdot \bar{\mathbf{u}}_r, \\
\bar{\mathbf{u}}_j &= \mathbf{u}_j^n - \frac{\Delta t}{m_j} \sum_{r \in \mathcal{R}_j} \bar{\mathbf{F}}_{jr}, \\
\bar{\eta}_j &= \eta_j^n,
\end{aligned}
\right.
\qquad (2.45)
$$

where the implicit fluxes are given by

$$
\begin{align}
\forall j \in \mathcal{J}, \forall r \in \mathcal{R}_j, \quad \bar{\mathbf{F}}_{jr} &= \mathbf{C}_{jr} \bar{\mathbf{p}}_j + A_{jr} (\bar{\mathbf{u}}_j - \bar{\mathbf{u}}_r), \nonumber \\
\text{and } \forall r \in \mathcal{R}, \quad \sum_{j \in \mathcal{J}_r} \bar{\mathbf{F}}_{jr} &= \mathbf{0}. \tag{2.46}
\end{align}
$$

Here the notation $\bar{\phi}$ indicates the implicit treatment of the term $\phi$. If at this point, the time discretization of the $C_{jr}$ vectors is not precised, one should observe that in dimension 1, these are the constants 1 or -1. The explicit $A_{jr} = A_{jr}^n$ matrices are symmetric non-negative and such that $A_r := \sum_{j \in \mathcal{J}_r} A_{jr}$ are invertible (one can consider for instance the matrices defined in (2.20) or (2.21)). Let us finally remark that during this phase, writing the pressure $p$ as a function of $\tau$ and $\eta$, since $\bar{\eta}_j = \eta_j^n$, one has

$$
\forall j \in \mathcal{J}, \quad \bar{p}_j = p(\bar{\tau}_j, \eta_j^n) = p_{\eta_j^n}(\bar{\tau}_j). \qquad (2.47)
$$

<footnote n="29"><sup>29</sup>For the kind of applications we have in mind, we believe that an accurate calculation of the acoustic waves inside the thin structures is not relevant.</footnote>
After this step, one knows the implicit fluxes and the correction step consists simply in

$$ \forall j \in \mathcal{J}, \left. \begin{array}{l} \tau_j^{n+1} = \tau_j^n + \frac{\Delta t}{m_j} \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr} \cdot \bar{\mathbf{u}}_r, \\ \mathbf{u}_j^{n+1} = \mathbf{u}_j^n - \frac{\Delta t}{m_j} \sum_{r \in \mathcal{R}_j} \bar{\mathbf{F}}_{jr}, \\ E_j^{n+1} = E_j^n - \frac{\Delta t}{m_j} \sum_{r \in \mathcal{R}_j} \bar{\mathbf{F}}_{jr} \cdot \bar{\mathbf{u}}_r, \end{array} \right. \qquad (2.48) $$

where, one may have noticed that $\tau_j^{n+1} = \bar{\tau}_j$ and $\mathbf{u}_j^{n+1} = \bar{\mathbf{u}}_j$.

It is easy to show that if $(\bar{\tau}_j, \bar{\mathbf{u}}_j, \bar{\eta}_j)_{j \in \mathcal{J}}$ is a solution of the predictor step (2.45)–(2.46), then the predictor-corrector scheme is conservative in volume, mass, momentum and total energy. Moreover, in dimension 1, we prove that the scheme is entropy stable<fnref n="30" /> for any $\Delta t \ge 0$.

Thus, it remains to show that the scheme is well-defined. Observing that, provided that $\frac{\partial}{\partial \tau} p|_{\eta} < 0$, in other words that the problem is hyperbolic, one can express $\forall j \in \mathcal{J}, \bar{\tau}_j = \tau_{\eta_j}^n(\bar{p}_j)$ and then, injecting the fluxes (2.46) into (2.45), one gets

$$ \forall j \in \mathcal{J}, \left. \begin{array}{l} \frac{m_j}{\Delta t} (\tau_{\eta_j^n}^n(\bar{p}_j) - \tau_j^n) - \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr}^T A_r^{-1} \sum_{i \in \mathcal{J}_r} \mathbf{C}_{ir} \bar{p}_i = \sum_{r \in \mathcal{R}_j} \mathbf{C}_{jr}^T A_r^{-1} \sum_{i \in \mathcal{J}_r} A_{ir} \bar{\mathbf{u}}_i, \\ \frac{m_j}{\Delta t} (\bar{\mathbf{u}}_j - \mathbf{u}_j^n) - \sum_{r \in \mathcal{R}_j} A_{jr} A_r^{-1} \sum_{i \in \mathcal{J}_r} A_{ir} \bar{\mathbf{u}}_i + \sum_{r \in \mathcal{R}_j} A_{jr} \bar{\mathbf{u}}_j = \sum_{r \in \mathcal{R}_j} A_{jr} A_r^{-1} \sum_{i \in \mathcal{J}_r} \mathbf{C}_{jr} \bar{p}_i. \end{array} \right. \qquad (2.49) $$

In the case of a perfect gas law, one can show that finding a solution $(p_j, \mathbf{u}_j)_{j \in \mathcal{J}}$ of (2.49) such that $\forall j \in \mathcal{J}, p_j > 0$, can be rewritten as

$$ \text{find } U = \begin{pmatrix} (-p_j)_{j \in \mathcal{J}} \\ (\mathbf{u}_j)_{j \in \mathcal{J}} \end{pmatrix} \in \mathcal{D}, \quad \text{such that} \quad \nabla J(U) = AU, \qquad (2.50) $$

where $\nabla J(U)$ and $AU$ are defined respectively by the left and right hand sides of (2.49), and where $N = \#\mathcal{J}$. The domain $\mathcal{D}$ is defined by $\mathcal{D} := ] -\infty, 0[ \times \mathbb{R}^{dN}$.

If the matrix $A$ is skew-symmetric and if $J$ is a strictly convex function, and under some additional hypothesis that are not recalled here, we prove in [21] that (2.50) admits a unique solution in $\mathcal{D}$. In the case of perfect gases, one can apply this abstract result to conclude that the scheme is well-defined.

To sum up, in dimension 1 and in the case of a perfect gas law, the implicit semi-Lagrangian scheme (2.45),(2.46),(2.48) is well-defined, consistent, conservative and unconditionally stable: $\rho$ and $\epsilon$ remain positive and the entropy increases for any $\Delta t \ge 0$.

This can be extended to more general equations of state<fnref n="31" />. In [21], we show how to adapt it to the case of stiffened gases.

In practice, the non-linear problem (2.49) is solved using a Newton method that converges to the machine precision in just a few iterations (from 2 to 6 iterations depending on the size of $\Delta t$). Numerical tests show that CFL numbers of a few hundreds can be

<footnote n="30">³⁰Actually the predictor scheme satisfies the following inequality $\bar{E}_j \le E_j - \frac{\Delta t}{m_j} \sum_{r \in \mathcal{R}_j} \bar{\mathbf{F}}_{jr} \cdot \bar{\mathbf{u}}_r$. It is consistent with $\rho D_t E + \nabla \cdot p \mathbf{u} \le 0$, and one can verify that $E$ is a mathematical entropy for the isentropic Euler equations, see [CCM10].</footnote>

<footnote n="31">³¹As soon as $p_{\eta_j^n} : \tau \mapsto p_{\eta_j^n}(\tau)$ is a continuous and strictly convex function, the proof that the scheme is well-defined requires only a few adjustments. If $p_{\eta_j^n}$ is not strictly convex, but strictly decreasing, the same kind of results should be obtained.</footnote>
used without any stability problem, which illustrates the unconditional stability of the scheme. As expected, using very large CFL numbers deteriorates the quality of the approximation but, surprisingly, we observe on various tests that contact discontinuities remain precisely located whatever the CFL number is. In [21], we give a theoretical explanation in the case of a 1D Riemann problem. This is actually a very important feature. Since we aim at solving fluid-structure interaction, the location of the interface is essential. Finally, in [21], we also explain how to couple implicit and explicit regions in the same calculation. The obtained implicit-explicit scheme is conservative, consistent and entropy stable.

The ongoing work focuses on the extension to the dimension 2 of the method. As mentioned earlier, in the definition of the predictor scheme (2.45)–(2.46), the time discretization of the $C_{jr}$ vectors is not set$^{32}$. It is important to notice that the predictor scheme is well-defined whatever the choice of the $C_{jr}$ vectors is. However, in dimension 2 or 3, choosing an explicit value $C_{jr} = C_{jr}^n$ defines a scheme that is not entropy stable in general. This a consequence of the time discretization error in the GCL. Actually, considering explicit $C_{jr}$ vectors one obtains

$$ \rho_j^{n+1} = \frac{m_j}{V_j^{n+1}} \neq \frac{1}{\tau_j^{n+1}}, $$

which breaks the entropy stability analysis. Nevertheless, in dimension 2, the GCL error is a linear function of $\Delta t$. Thus, since

$$ C_{jr} := \frac{1}{2} (C_{jr}^n + C_{jr}^{n+1}) \quad \Rightarrow \quad \rho_j^{n+1} = \frac{m_j}{V_j^{n+1}} = \frac{1}{\tau_j^{n+1}}, \quad (2.51) $$

the entropy stability is obtained also in 2D$^{33}$. From the theoretical point of view, it remains to show that the scheme with the new non-linearity imposed by (2.51) remains well-defined. The first numerical tests are encouraging.

## 2.3 Mesh adaptation for semi-Lagrangian compressible flows

As stated previously, when I joined the Troll [Lef+18] team, my first mission was to add and evaluate a layer-based mesh adaptation technique. To be able to define layers, the code requires quasi-logical meshes and the idea is just to collapse layers that are considered too small, and to split the ones that are judged too large. The criteria are purely geometric ones. The method is quite simple$^{34}$ but for quasi-1D flows it is very efficient. The method was not published by itself, but we provide an illustration of its good behavior in [3], see Figure 2.12.

Actually, during my PhD Thesis, I shared the office of N. Dicesare [Di 00] and then of C. Dobrzynski [Dob05]. In both of their PhD Thesis, they used metric-based conforming AMR as a technique to reduce the cost of their simulations while maintaining accuracy. Thus, I had the feeling that it could be interesting to use this kind of approach in the semi-Lagrangian context.

Obviously, adapting the mesh (by changing its connectivity) is not a new idea in the semi-Lagrangian community, and it has been investigated a lot in the 90s (see [FCT85] for instance). In my view, it was not as successful or as popular compared to classic ALE methods for two reasons. First, changing the mesh connectivity is more difficult

<footnote>$^{32}$In dimension 1, there is no choice, the value of $C_{jr}$ is $\pm 1$ according to the node r location in the cell j</footnote>

<footnote>$^{33}$In 3D, the GCL error is a quadratic function of $\Delta t$, thus another approach must be used to ensure entropy stability.</footnote>

<footnote>$^{34}$The main difficulty resided in its implementation in a massively parallel and multi-physics code.</footnote>
<figure><img src="image_25.png" /><figcaption>(a) Meshes at final time $t = 0.2$. **Left:** purely semi-Lagrangian simulation. **Right:** with layer refinement.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(b) Radial cut of the density at final time $t = 0.2$. The simulation using layer refinement is closer to the reference.</figcaption></figure>

Figure 2.12: Illustration of the effect of the dynamic layer refinement on a converging Sod shock tube in 3D.

and more expensive than simply displacing its nodes. Second, providing a conservative remapping using staggered semi-Lagrangian schemes is more difficult for this kind of methods.

However, due to the improvements of calculators, to the progresses made in remeshing techniques (for both triangular and tetrahedral meshes), and to the emergence of finite-volume semi-Lagrangian schemes [DM05; Maz07; Mai+07], it appeared that it was worth looking again at conforming<fnref n="35" /> AMR-like methods in the semi-Lagrangian context.

Also, discussing with C. Dobrzynski at the beginning of her PhD Thesis [Dob05], I felt that the remeshing method that they proposed with P. Frey was a good candidate to deal with mesh adaptation in parallel: they used local mesh modifications in 3D, to achieve adaptation (i.e. not a global remeshing approach).

Thus, as an extension of my previous work, I began to investigate in this direction. At first, we proposed a CEMRACS subject with B. Després [1] on a related topic. The idea was to evaluate the ability of a simple local mesh adaptation strategy to deal with the free fall of a droplet in the air. This work was done in collaboration with É. Bernard, E. Deriaz, B. Després, K. Jurkova and F. Lagoutière.

Meanwhile, I began the development of a metric based AMR-like method in 2D. In this work, published in [10], when dealing with multi-material flows, the mesh adaptation procedure is constrained to preserve pure cells (i.e. each cell of the mesh contains a single material all along the calculation). The problem with this approach is that this remeshing constraint can produce very small cells at the material interfaces. For explicit methods, it may result in arbitrary small time steps than can prevent calculation to reach the final time. So, with I. Marmajou, we improved the method, allowing the creation of multi-material cells and extending the remapping to the second-order. This is the subject of [16].

It is worth noting that when I started working on [10], the study of two other approaches dealing with remeshing by means of connectivity changes also began. Thus at that time, I had discussions on this topic with R. Loubère and with P. Hoch. Actually, R. Loubère and his co-authors proposed the ReALE<fnref n="36" /> method [Lou+10] which in short consists in maintaining a Voronoï grid in a quasi-Lagrangian way all along the calculation. In his work, P. Hoch [Hoc12] proposed a very ambitious method based on polygonal mesh adaptation. This is a very interesting approach since it gives a lot of freedom to the

<footnote n="35">³⁵In this manuscript I focus on conforming methods and I will not discuss the work of R.W. Anderson, N.S. Elliott and R.B. Pember [AEP04], where they propose an ALE method combined with patch-based AMR.</footnote>

<footnote n="36">³⁶This is an acronym for a Reconnection based ALE method.</footnote>
mesh modifications, the difficulty is that it does not relies (as the method we proposed and as the ReALE method) on decades of researches on mesh adaptation. A lot needs to be invented.

Let us describe in a few words the method we proposed in [10; 16]. First of all, fol-
lowing the seminal work of [Bor+97], the mesh adaptation strategy that we use is based
on the construction of unit meshes in a non-euclidean metric space. In other words, the
idea is to build a mesh for which all edges have a length of 1 in a prescribed metric space.
The metric itself is defined in such a way that it controls the adaptation. Thus, one of
the ingredients is the definition of this metric field. Since no a posteriori error estimate
has been yet defined for the approximation of solutions of Euler equations, we use a ge-
ometric error estimate. The idea for such an estimator is simply to put more cells where
the curvature of the solution is high. Actually, the curvature of the graph of a function
and the metric field that drives the mesh refinement is quite straightforward. Assuming
for instance that one wants to adapt the mesh to some scalar field φ, one estimates its
Hessian matrix field ∇²φ = P⁶Λ⁶P⁶⁻¹, where Λ is the diagonal matrix of its eigenvalues
(λᵢ)₁≤i≤d. Then the metric field associated to φ is simply

$$
M^{\phi} := P^{\phi} |\Lambda^{\phi}| P^{\phi^{-1}}.
$$

This allows to define the length<fnref n="37" />

$$
\forall \mathbf{a}, \mathbf{b} \in \mathbb{R}^d, \quad l^\phi(\mathbf{a}, \mathbf{b}) := \int_0^1 (\mathbf{b} - \mathbf{a}, M_{\mathbf{a}+t(\mathbf{b}-\mathbf{a})}^\phi (\mathbf{b}-\mathbf{a}))^{\frac{1}{2}} dt.
$$

This is the formula that we use to check if an edge of the mesh is too long or too short.
Thus, a mesh will be said to be a unit mesh in the metric M^φ if, for each of its edges e,
one has

$$
\sqrt{\frac{1}{2}} \le l^{\phi}(\mathbf{x}_{1}^{e}, \mathbf{x}_{2}^{e}) \le \sqrt{2},
$$

where $\mathbf{x}_1^e$ and $\mathbf{x}_2^e$ are the coordinates of the extremities of $e$. In practice, the formula is
altered to take into account the fact that $l^\phi(\cdot, \cdot)$ is not a distance, see [10; 16] for details.

Actually, this metric formalism is very flexible. For instance, defining an isotropic
metric just consists in replacing the previous definition of M$^{\Phi}$ by

$$
M^{\phi} = \max_{1 \le i \le d} |\lambda_i| I.
$$

It is also easy to take multiple adaptation criteria into account [Bor+97; AF03].

In [10; 16], we define an isotropic mesh adaptation technique based on several physi-
cal criteria (e.g. ρ, u, ε, ...).

As said previously, following [Dob05], a key ingredient in our mesh adaptation method
is the use of three simple local mesh modification patterns:

* a mesh quality pattern: *edge swapping* (see Figure 2.13a),

• a mesh refinement pattern: edge splitting (see Figure 2.13b), and

* a mesh coarsening pattern: *edge collapsing* (see Figure 2.13c).

$$
d_{M^{\phi}}(\mathbf{a}, \mathbf{b}) = \inf_{\gamma} \int_{0}^{1} \left( \gamma'(t), M_{\gamma(t)}^{\phi} \gamma'(t) \right)^{1/2} dt,
$$

where $\gamma$ is a $C^1$-path joining $\mathbf{a}$ to $\mathbf{b}$ ($\gamma(0) = \mathbf{a}$ and $\gamma(1) = \mathbf{b}$).

<footnote n="37"><sup>37</sup>Observe that as soon as M<sup>Φ</sup> is not uniform, l<sup>Φ</sup> is not a distance: since geodesics are not straight lines in that case, l<sup>Φ</sup> does not satisfy triangle inequality. In a metric field, the natural distance is the Riemannian distance defined by</footnote>
<figure><img src="image_27.png" /><figcaption>Figure 2.13: The three mesh modification patterns.</figcaption></figure>

There are two reasons why the use of local remeshing is important in our work.

The first one is a direct consequence of the use of a semi-Lagrangian scheme. Semi-Lagrangian methods tend to produce a somehow adapted mesh: the mesh gets finer in the vicinity of shocks and it just follows the flow at contact discontinuities. Thus one can assume that if a mesh $M^n$ is adapted to some metric field $M^n$ at time $t^n$, then producing an adapted mesh $M^{n+1}$ to the metric field $M^{n+1}$ at time $t^{n+1}$ should require only a few modifications of $M^n$. In other words, one can hope that the metric $M^{n+1}$ is close in some sense to the metric $M^n$ transported by the flow. This is what we observe in practice: only a few cells are changed from one time step to the next one. To fix ideas in our numerical experiments, only a few percent of the cells are changed by the adaptation procedure<fnref n="38" />. This presents the advantage to keep the method cheap<fnref n="39" /> and it produces low numerical dissipation.

The second advantage of local mesh adaptation resides in its implementation itself. Indeed, it allows to define a parallel method. Even more, in our implementation, we ensure that the numerical results, with regard to the semi-Lagrangian flow and mesh adaptation, are exactly the same (bit-to-bit) whichever the number of MPI processes is. This is an important feature in an industrial context.

For details, with regard to our mesh adaptation strategy in the context of the multi-material flows, one should consult [16]. These aspects are not recalled here.

Let us now describe how data are remapped between meshes. Again, it relies strongly on the use of local mesh modification patterns. As illustrated on Figure 2.13 each mesh change is a simple modification of a small region of the mesh. A set of triangles is replaced by a new one that keeps unchanged the geometry of the computational domain<fnref n="40" />. So, the remapping is performed between each pair of old and new cavities, it is thus a local operation. It is cheap since it is performed only where it is needed. The remapped quantities are the conservative ones ($\rho$, $\rho u$ and $\rho E$ in the case of gas dynamics for in-

<footnote n="38">³⁸It clearly depends on the refinement: for coarse meshes it can be around 5% of the cells and around 0.1% on finer grids.</footnote>

<footnote n="39">³⁹Adapting the mesh at each iteration costs 3 to 5 times the cost of the second-order semi-Lagrangian scheme itself.</footnote>

<footnote n="40">⁴⁰In [1; 10; 16] the mesh boundaries are straight. This is not a limitation since one can embed the computational domain $\Omega$ in a box $D$, filling $D \setminus \Omega$ with a fictitious fluid for instance.</footnote>
stance). In [16], we show that the first-order remapping method is conservative and satisfies Maximum Principles<fnref n="41" /> for $ρ, (u_i)_{1\le i\le d}$ and $ϵ = E - \frac{1}{2}\|u\|^2$. Following [HL14], it allows to define a second-order accurate conservative remapping method that is stable: we use the APITALI method (see [BHS20]) to ensure Maximum Principles.

To conclude this section, we precise that in [10] and [16], metric evaluation, local mesh modifications and remapping are intertwined<fnref n="42" />. This allows the convergence of the adaptation loop. One should refer to these two papers for more practical information concerning the method.

We finally illustrate the method’s behavior using two test cases.

The first test consists in the interaction of a planar shock wave in 2D, initially located in the air, with a set of bubbles of another material. Two materials are considered for the bubbles. In the first case, see Figures 2.14b and 2.14c, the bubbles are made of a refrigerant material that behaves like a perfect gas with $γ = 1.25$. In the second case, see Figures 2.14d and 2.14e, one considers helium bubbles which also follow a perfect gas law with $γ = 1.67$. This test has been proposed in [Ban+07]. On Figure 2.14 we compare our solutions for both cases with the ones they obtain at final time $t = 6$. The results are quite similar.

The second test has been designed with I. Marmajou. It is set to assess the robustness and the flexibility of our method. The test consists in two rigid bars rotating in a fluid domain. The fluid itself is made of three constituents of various densities and initially at rest. Dealing with Euler equations it is natural to impose perfect sliding at the bars boundaries. The results we obtain are reproduced on Figure 2.15, the dynamics of the flow is represented by four snapshots at times $\pi/2, \pi, 3\pi/2$ and $2\pi$.

<footnote n="41">⁴¹More precisely, in [16], we establish a Maximum Principle for all the remapped conservative quantities and the deduced specific quantities.</footnote>

<footnote n="42">⁴²After each atomic mesh modification, the numerical solution is remapped locally and the metric field is updated accordingly. One then performs the next mesh modification if needed.</footnote>
<figure><img src="image_28.png" /><figcaption>(a) Initial geometry for the shock interaction with cylindrical bubbles.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(b) Our solution for refrigerant bubbles.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(c) Reference for refrigerant bubbles.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(d) Our solution for helium bubbles.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(e) Reference for helium bubbles.</figcaption></figure>

<caption>Figure 2.14: Shock interaction with cylindrical bubbles.</caption>
<figure><img src="image_33.png" /><figcaption>(a) Initial configuration. The fluid is at rest: $\mathbf{u} = \mathbf{0}$ and $p = 1$. Each fluid follows a perfect gas law ($\gamma = 1.4$).</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(b) $t = \pi/2$.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>(c) $t = \pi$.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(d) $t = 3\pi/2$.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(e) $t = 2\pi$.</figcaption></figure>

Figure 2.15: The mixer test. The geometry is defined on Figure 2.15a. Figures 2.15b to 2.15e show the evolution of the mesh according to time. Meshes are colored by the fraction presence of the middle material.
Conclusion and perspectives

In this manuscript, I presented an overview of my works. Although the two main Chap-
ters are quite different (by their sizes and by the topics they cover), they reflect quite well
the influences that drove my researches.

The most representative example is probably the way I proposed to deal with curvi-
linear cells [9], which are required to define very high-order semi-Lagrangian schemes
for compressible flows. The influence of my finite-element background (coming from the
elliptic community) is quite obvious and defines an original way to compute the fluxes
with regard to hyperbolic problems. If we already applied successfully this framework
to define a conservative slide-line method in [2], and if we already presented first high-
order results in conferences, this work is not finished yet. Following our first encouraging
results, I will soon invest more in this topic with P. Hoch, E. Labourasse and I. Marmajou.
This should lead to new researches: indirect ALE is probably the more natural extension,
treatments of elastic-plastic flows is most likely straightforward with this framework and
defining finite-volume schemes to treat diffusion on curvilinear meshes could also be
studied. In my view, this class of schemes offers a large variety of interesting extensions.

Another research direction is the continuation of the mesh adaptation technique that takes advantage of the semi-Lagrangian solver. As we observed, this approach reduces the number of mesh modifications during adaptation. The next natural work is the extension to 3D. The main difficulty (not considering the implementation) is related to the quality pattern (see Figure 2.13a). Indeed, considering an edge *e*, the number of local remeshing possibilities is given by the Catalan number $\frac{(2n-2)!}{n!(n-1)!}$, where *n* is the number of tetrahedrons connected to *e*. With P. Hoch and I. Marmajou, we shall investigate on how to overcome this difficulty.

Even if it was not covered in this manuscript, with G. Carré and E. Labourasse, we have supervised two internship subjects dedicated to the approximation of compressible Navier-Stokes equations in semi-Lagrangian coordinates. In the first internship, we supervised the work of F. Chopot [Cho16] who studied an Euler/Navier-Stokes coupling<sup>43</sup> in 1D, as a simplified model for compactly supported turbulence. This work was then extended during the internship of J. Patela [Pat20] to dimensions 2 and 3. On the way, we defined a finite-volume scheme for linear elasticity. A publication is in preparation.

In another internship, that I co-supervised with C. Buet, we studied with V. Four-
net [Fou21] an explicit nodal finite-volume scheme for the approximation of the P<sub>N</sub> par-
ticles transport model in 2D. This is an extension for N > 1 of [BDF12], where it is shown
that this kind of scheme preserves the asymptotic on unstructured grids while it is not
the case for classical face-based solvers. A publication is under preparation. Moreover,

<footnote n="43"><sup>43</sup>Molecular diffusion and thermal conduction could be defined on bounded subdomains.</footnote>
an extension of this work, based on the framework I defined in [9], is the subject of a new
internship we proposed with C. Buet.

Also, the results obtained by A. Plessier in her PhD Thesis [Ple23], that I co-supervise with B. Després, are very promising: the implicit scheme for semi-Lagrangian compressible flows is unconditionally stable and very accurate with regard to contact discontinuities [21]. Thus many directions of research are already envisaged: the extension to 3D, the study of a second-order version of the scheme, the definition of a similar scheme in Eulerian coordinates or in ALE, and more obviously the treatment of elastic-plastic flows which will allow us to define an efficient monolithic fluid-structure conservative and entropic solver. In this direction we already have defined with B. Després and E. Labourasse a few internship subjects and a PhD subject is in preparation.

Finally, dealing also with fluid-structure interaction and elastic-plastic flows, I am co-supervising, with R. Abgrall and E. Labourasse, the PhD Thesis of A. Drouard that just began in November 2022. The idea is, following [AT20; AT22], to define and to analyze kinetic schemes to treat hyperelasticity. The aim of this work is to obtain an efficient method that works for large CFL numbers.
# Bibliography

## Publications

[1] É. Bernard, S. Del Pino, E. Deriaz, B. Després, K. Jurkova, and F. Lagoutière. "Lagrangian method enhanced with edge swapping for the free fall and contact problem". In: *ESAIM: Proceedings*. Vol. 24. EDP Sciences. 2008, pp. 46-59 (cit. on pp. 35, 37).

[2] S. Bertoluzza, S. Del Pino, and E. Labourasse. "A conservative slide line method for cell-centered semi-Lagrangian and ALE schemes in 2D". In: *ESAIM: Mathematical Modelling and Numerical Analysis* 50.1 (2016), pp. 187-214 (cit. on pp. 24, 27, 41).

[3] G. Carré, S. Del Pino, B. Després, and E. Labourasse. "A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension". In: *Journal of Computational Physics* 228.14 (2009), pp. 5160-5183 (cit. on pp. 16-19, 34).

[4] G. Carré, S. Del Pino, K. Pichon Gostaf, E. Labourasse, and A. V. Shapeev. "Polynomial Least-Squares reconstruction for semi-Lagrangian cell-centered hydrodynamical schemes". In: *ESAIM: Proceedings*. Vol. 28. EDP Sciences. 2009, pp. 100-116 (cit. on pp. 19, 20).

[5] T. Chacón Rebollo, S. Del Pino, and D. Yakoubi. "An iterative procedure to solve a coupled two-fluids turbulence model". In: *ESAIM: Mathematical Modelling and Numerical Analysis* 44.4 (2010), pp. 693-713 (cit. on p. 6).

[6] S. Del Pino. "3D computing using virtual reality data". In: *ESAIM: Proceedings*. Vol. 10. EDP Sciences. 2001, pp. 267-275 (cit. on p. 5).

[7] S. Del Pino. "Une méthode d'éléments finis pour la résolution d'EDP dans des domaines décrits par géométrie constructive". PhD thesis. Laboratoire Jacques-Louis Lions, Université Pierre et Marie Curie, Dec. 2002 (cit. on p. 4).

[8] S. Del Pino. "A hierarchical and view dependent visualization algorithm for tree based AMR data in 2D or 3D". In: *Proceedings of the 5th Eurographics conference on Parallel Graphics and Visualization*. Eurographics Association. 2004, pp. 49-58 (cit. on p. 1).

[9] S. Del Pino. "A curvilinear finite-volume method to solve compressible gas dynamics in semi-Lagrangian coordinates". In: *Comptes Rendus Mathématique* 348.17-18 (2010), pp. 1027-1032 (cit. on pp. 1, 22, 23, 41, 42).

[10] S. Del Pino. "Metric-based mesh adaptation for 2D Lagrangian compressible flows". In: *Journal of Computational Physics* 230.5 (2011), pp. 1793-1821 (cit. on pp. 35-38).
[11] S. Del Pino, B. Després, P. Havé, H. Jourdren, and P.-F. Piserchia. "3D finite volume simulation of acoustic waves in the earth atmosphere". In: *Computers & Fluids* 38.4 (2009), pp. 765-777 (cit. on p. 11).

[12] S. Del Pino, E. Heikkola, O. Pironneau, and J. Toivanen. "A finite element method for virtual reality data". In: *Comptes Rendus de l'Académie des Sciences-Series I-Mathematics* 330.12 (2000), pp. 1107-1111 (cit. on p. 1).

[13] S. Del Pino and H. Jourdren. "Arbitrary high-order schemes for the linear advection and wave equations: application to hydrodynamics and aeroacoustics". In: *Comptes Rendus Mathematique* 342.6 (2006), pp. 441-446 (cit. on pp. 9, 11).

[14] S. Del Pino, E. Labourasse, and G. Morel. "An asymptotic preserving multidimensional ALE method for a system of two compressible flows coupled with friction". In: *Journal of Computational Physics* 363 (2018), pp. 268-301 (cit. on pp. 29, 30).

[15] S. Del Pino, J.-L. Lions, and O. Pironneau. "A-priori Domain Decomposition of PDE Systems and Applications". In: *Mathematical Modeling and Numerical Simulation in Continuum Mechanics*. Springer, 2002, pp. 125-135 (cit. on p. 1).

[16] S. Del Pino and I. Marmajou. "Triangular metric-based mesh adaptation for compressible multi-material flows in semi-Lagrangian coordinates". In: *Journal of Computational Physics* 478 (Apr. 2023) (cit. on pp. 35-38).

[17] S. Del Pino and B. Maury. "2D/3D turbine simulations with FreeFEM". In: *Numerical Analysis and Scientific Computing for PDEs and their challenging applications*. Ed. by J. Haataja, R. Stenberg, J. Periaux, P. Raback, and P. Neittaanmaki. 2007 (cit. on p. 5).

[18] S. Del Pino and O. Pironneau. "A fictitious domain based general PDE solver". In: *Numerical methods for scientific computing variational problems and applications, Barcelona* (2003) (cit. on pp. 4, 5).

[19] S. Del Pino and O. Pironneau. "Asymptotic analysis and layer decomposition for the couplex exercise". In: *Computational Geosciences* 8.2 (2004), pp. 149-162 (cit. on p. 5).

[20] S. Del Pino, U. Razafison, and D. Yakoubi. "Une borne inférieure pour la constante de la condition inf-sup sur l'opérateur de divergence". In: *Comptes Rendus Mathematique* 346.9-10 (2008), pp. 533-538 (cit. on p. 6).

[21] A. Plessier, S. Del Pino, and B. Després. "Implicit discretization of Lagrangian gas dynamics". In: *ESAIM: Mathematical Modelling and Numerical Analysis* 57 (2023), pp. 717-743 (cit. on pp. 32-34, 42).

## References

[Abg94] R. Abgrall. "On essentially non-oscillatory schemes on unstructured meshes: analysis and implementation". In: *J. Comput. Phys.* (1994) (cit. on p. 20).

[AT20] R. Abgrall and D. Torlo. "High order asymptotic preserving deferred correction implicit-explicit schemes for kinetic models". In: *SIAM Journal on Scientific Computing* 42.3 (2020), B816-B845 (cit. on p. 42).

[AT22] R. Abgrall and D. Torlo. "Some preliminary results on a high order asymptotic preserving computationally explicit kinetic scheme". In: *Communications in Mathematical Sciences* 2 (2022), pp. 297-326 (cit. on p. 42).
[Add+90] F. L. Addessio, J. R. Baumgardner, J. K. Dukowicz, N. L. Johnson, B. A. Kashiwa, R. M. Rauenzahn, and C. Zemach. *CAVEAT: A Computer Code for Fluid Dynamics Problems with Large Distortion and Internal Slip*. Tech. rep. Los Alamos National Laboratory LA-10613, 1990 (cit. on pp. 15, 23).

[AF03] F. Alauzet and P.-J. Frey. *Estimation d'erreur géométrique et métriques anisotropes pour l'adaptation de maillage. Partie I : aspects théoriques*. Tech. rep. 4753. INRIA, 2003 (cit. on p. 36).

[AEP04] R.W. Anderson, N.S. Elliott, and R.B. Pember. "An arbitrary Lagrangian-Eulerian method with adaptive mesh refinement for the solution of the Euler equations". In: *J. Comput. Phys.* 199.2 (2004), pp. 598-617 (cit. on p. 35).

[Att+95] K. Attenborough, S. Taherzadeh, H. E. Bass, X. Di, R. Raspet, G.R. Becker, A. Güdesen, A. Chrestman, G. A Daigle, A. L'Espérance, et al. "Benchmark cases for outdoor sound propagation models". In: *The Journal of the Acoustical Society of America* 97.1 (1995), pp. 173-191 (cit. on p. 11).

[Bab73] I. Babuška. "The finite element method with penalty". In: *Math. Comp.* 27 (1973), pp. 221-228 (cit. on p. 4).

[Ban+07] J.W. Banks, D.W. Scwendeman, A.K. Kapila, and W.D. Henshaw. "A high-resolution method for compressible multi-marterial flow on overlapping grids". In: *J. Comput. Phys.* 223 (2007), pp. 262-297 (cit. on p. 38).

[Ben99] F. Ben Belgacem. "The mortar finite element method with Lagrange multipliers". In: *Numer. Math.* 84 (1999), pp. 173-199 (cit. on p. 25).

[Ber+12] A. Bernard-Champmartin, E. Deriaz, P. Hoch, G. Samba, and M. Schaefer. "Extension of centered hydrodynamical schemes to unstructured deforming conical meshes: the case of circles". In: *ESAIM: Proc.* 38 (2012), pp. 135-162 (cit. on p. 20).

[BHS20] A. Bernard-Champmartin, P. Hoch, and N. Seguin. *Stabilité locale et montée en ordre pour la reconstruction de quantités volumes finis sur maillages coniques non-structurés en dimension 2*. Research Report. Mar. 2020. URL: https://hal.science/hal-02497832 (cit. on pp. 15, 38).

[BM97] C. Bernardi and Y. Maday. *Spectral methods*. North-Holland, Amsterdam, 1997 (cit. on p. 5).

[BMP94] C. Bernardi, Y. Maday, and A. T. Patera. "A new nonconforming approach to domain decomposition: the mortar element method". In: *Nonlinear Partial Differential Equations and Their Applications* (1994). (H. Brezis and J. L. Lions eds.), Pitman, New York, pp. 13-51 (cit. on p. 25).

[BIM11] S. Bertoluzza, M. Ismail, and B. Maury. "Analysis of the fully discrete fat boundary method". In: *Numer. Math.* 118 (2011), pp. 49-77 (cit. on p. 4).

[Bor+97] H. Borouchaki, P.-L. George, F. Hecht, P. Laug, and E. Saltel. "Delaunay mesh generation governed by metric specifications". In: *Finite Elements in Analysis and Design* 25 (1997), pp. 61-83 (cit. on p. 36).

[BK05] N. G. Bourago and V. N. Kukudzhanov. *A review of contact algorithms. The Institute for problems in mechanics of RAS, Izv. RAN, MTT, No. 1, pp. 45-87*. translation into english. 2005 (cit. on p. 23).

[BM92] F. Brezzi and L. D. Marini. "Macro hybrid elements and domain decomposition methods". In: *Optimisation et Contrôle, Meeting in honour of J. Céa*. Ed. by J. Desideri et al. (CÉPADUÈS-Edition, Toulouse, 1993), Apr. 1992, p. 89 (cit. on p. 25).
[BDF12] C. Buet, B. Després, and E. Franck. "Design of asymptotic preserving finite volume schemes for the hyperbolic heat equation on unstructured meshes". In: *Numer. Math.* 122 (2012), pp. 227–278 (cit. on p. 41).

[Bur90] D. Burton. *Exact conservation of energy and momentum in staggered-grid hydro-dynamics with arbitrary connectivity*. Tech. rep. UCRL-JC-105926. Lawrence Livermore National Laboratory, 1990 (cit. on p. 15).

[Car09] E. J. Caramana. "The implementation of slide lines as a combined force and velocity boundary condition". In: *J. Comput. Phys.* 228 (2009), pp. 3911–3916 (cit. on pp. 27, 28).

[Car+98] E. J. Caramana, D. E. Burton, M. J. Shashkov, and P. P. Whalen. "The construction of compatible hydrodynamics algorithms utilizing conservation of total energy". In: *J. Comput. Phys.* 146 (1998), pp. 227–262 (cit. on pp. 15, 17).

[CY18] T. Chacón Rebollo and D. Yakoubi. "A three-dimensional model for two coupled turbulent fluids: numerical analysis of a finite element approximation". In: *IMA Journal of Numerical Analysis* 38.4 (2018), pp. 1927–1958 (cit. on p. 8).

[CCM10] C. Chalons, F. Coquel, and C. Marmignon. "Time-implicit approximation of the multipressure gas dynamics equations in several space dimensions". In: *SIAM* 48 (2010), pp. 1678–1706 (cit. on pp. 32, 33).

[CS07] J. Cheng and C.-W. Shu. "A high order ENO conservative Lagrangian type scheme for the compressible Euler equations". In: *J. Comput. Phys.* 227 (2007), pp. 1567–1596 (cit. on p. 19).

[Cho16] F. Chopot. "Couplage Euler/Navier-Stokes en 1D". MA thesis. Université de Nantes, 2016 (cit. on p. 41).

[CDL14] G. Clair, B. Després, and E. Labourasse. "A one-mesh method for the cell-centered discretization of sliding". In: *Comp. Meth. Appl. Mech. Engrg.* 269 (2014), pp. 315–333 (cit. on p. 27).

[Cla+12] A. Claisse, B. Després, E. Labourasse, and F. Ledoux. "A new exceptional points method with application to cell-centered Lagrangian schemes and curved meshes". In: *J. Comp. Phys.* 231 (2012), pp. 4324–4354 (cit. on p. 20).

[CT65] J. Cooley and J. Tukey. "An algorithm for the machine calculation of complex Fourier series". In: *Math. Comp.* 19 (1965), pp. 297–301 (cit. on p. 3).

[DT04] V. Daru and C. Tenaud. "High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations". In: *J. Comput. Phys.* 193.2 (2004), pp. 563–594. ISSN: 0021-9991 (cit. on pp. 9, 10).

[Des08] B. Després. "Stability of high order finite volume schemes for the 1D transport equation". In: *Finite Volumes for Complex Applications*. Vol. V. London: ISTE, 2008, pp. 337–342 (cit. on p. 10).

[DL12] B. Després and E. Labourasse. "Stabilization of cell-centered compressible Lagrangian methods using subzonal entropy". In: *J. Comp. Phys.* 231.20 (2012), pp. 6559–6595 (cit. on p. 27).

[DM05] B. Després and C. Mazeran. "Lagrangian Gas Dynamics in Two Dimensions and Lagrangian systems". In: *Arch. Rational Mech. Anal.* (2005) (cit. on pp. 13–16, 18, 22, 23, 35).

[Di 00] N. Di Césaré. "Outils pour l'optimisation de forme et le contrôle optimal : application à la mécanique des fluides". PhD thesis. Université Pierre et Marie Curie, 2000 (cit. on p. 34).
[DKR12] V. A. Dobrev, T. V. Kolev, and R. N. Rieben. "High-order curvilinear finite element methods for Lagrangian hydrodynamics". In: *SIAM Journal on Scientific Computing* 34.5 (2012), B606–B641 (cit. on pp. 15, 20).

[Dob05] C. Dobrzynski. "Adaptation de Maillage anisotrope 3D et application à l'aéro-thermique des bâtiments". PhD thesis. Université Pierre et Marie Curie — Paris VI, Nov. 2005 (cit. on pp. 34–36).

[Dub+10] F. Duboc, C. Enaux, S. Jaouen, H. Jourdren, and M. Wolff. "High-order dimensionally split Lagrange-remap schemes for compressible hydrodynamics". In: *Comptes Rendus Mathematique* 348.1 (2010), pp. 105–110 (cit. on p. 11).

[Ena07] C. Enaux. "Analyse mathématique et numérique d'un modèle multifluide multivitesse pour l'interpénétration de fluides miscibles". PhD thesis. École Centrale de Paris, 2007 (cit. on p. 29).

[FGG01] C. Farhat, P. Geuzaine, and C. Grandmont. "The Discrete Geometric Conservation Law and the Nonlinear Stability of ALE Schemes for the Solution of Flow Problems on Moving Grids". In: *Journal of Computational Physics* 174.2 (2001), pp. 669–694. ISSN: 0021-9991 (cit. on p. 16).

[Fou21] V. Fournet. "Construction, analyse et mise en œuvre d'un schéma Volumes Finis nodal pour le modèle $P_N$ dans le cadre du transport de particules". MA thesis. Sorbonne Université, 2021 (cit. on p. 41).

[FCT85] M. Fritts, W.P. Crowley, and H. Trease. "The Free-Lagrange Method". In: *Lecture Note in Physics*. Vol. 238. Springer-Verlag, Mar. 1985 (cit. on p. 34).

[GBM16] G. Georges, J. Breil, and P.-H. Maire. "A 3D GCL compatible cell-centered Lagrangian scheme for solving gas dynamics equations". In: *J. Comput. Phys.* 305 (2016), pp. 921–941 (cit. on pp. 16, 18).

[GR09] C. Geuzaine and J.-F. Remacle. "Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities". In: *International Journal for Numerical Methods in Engineering* 79.11 (2009), pp. 1309–1331 (cit. on p. 4).

[GG95] V. Girault and R. Glowinski. "Error analysis of a fictitious domain method applied to a Dirichlet problem". In: *Japan J. Indust. Appl. Math.* 12.487 (1995) (cit. on p. 3).

[GK98] R. Glowinski and Y. Kuznetsov. "On the solution of the Dirichlet problem for linear elliptic operators by a distributed Lagrange multiplier method". In: *Comptes Rendus de l'Académie des Sciences - Series I - Mathematics* 327.7 (1998), pp. 693–698. ISSN: 0764-4442 (cit. on p. 4).

[GPP94] R. Glowinski, T. Pan, and J. Periaux. "A fictitious domain method for Dirichlet problem and applications". In: *Comput. Methods Appl. Mech. Engrg.* 111 (1994), pp. 283–303 (cit. on p. 4).

[Goa60] W. B. Goad. *A Numerical Method for Two-Dimensional Unsteady Fluid Flow*. Tech. rep. LAMS 2365. Los Alamos National Laboratry, 1960 (cit. on p. 15).

[God59] S. K. Godunov. "A difference scheme for numerical computation of discountinuous solution of hydrodynamics equations". In: *Math. Sib.* 47 (1959) (cit. on p. 14).

[HSV12] A. Hannukainen, R. Stenberg, and M. Vohralík. "A unified framework for a posteriori error estimation for the Stokes problem". In: *Numer. Math.* (2012) (cit. on p. 6).
[Hei+98] E. Heikkola, Y. Kuznetsov, P. Neittaanmäki, and J. Toivanen. "Fictitious Domain Methods for the Numerical Solution of Two-Dimensional Scattering Problems". In: *J. Comput. Phys.* 145.1 (1998), pp. 89–109. ISSN: 0021-9991 (cit. on p. 4).

[HJJ09] O. Heuzé, S. Jaouen, and H. Jourdren. "Dissipative issue of high-order shock capturing schemes with non-convex equations of state". In: *J. Comput. Phys.* 228.3 (2009), pp. 833–860. ISSN: 0021-9991 (cit. on p. 11).

[Hoc12] P. Hoch. "A Discontinuous ALE formulation (DiscALE) for the modeling of polygonal mesh (r,h)-adaptation in finite volume context". preprint. 2012 (cit. on p. 35).

[HL14] P. Hoch and E. Labourasse. "A frame invariant and maximum principle enforcing second-order extension for cell-centered ALE schemes based on local convex hull preservation". In: *Int. J. Numer. Meth. Fluids* 76.12 (2014), pp. 1043–1063 (cit. on p. 38).

[Hoc+11] P. Hoch, P. Navaro, B. Boutin, and E. Deriaz. "Extension of ALE methodology to unstructured conical meshes". In: *ESAIM: Proc.* 32 (2011), pp. 31–55 (cit. on p. 20).

[Jou05] H. Jourdren. "HERA: A hydrodynamic AMR platform for multi-physics simulations". In: *Adaptive Mesh Refinement-Theory and Applications* (2005), pp. 283–294 (cit. on pp. 11, 31).

[Kid74] R. E. Kidder. "Theory of homogeneous isentropic compression and its application to laser fusion". In: *Nuclear Fusion* 14.1 (Jan. 1974), pp. 53–60 (cit. on p. 20).

[Klu08] G. Kluth. "Analyse mathématique et numérique de systèmes hyperélastiques et introduction à la plasticité". PhD thesis. Université Pierre et Marie Curie, 2008 (cit. on pp. 18, 31).

[KD10] G. Kluth and B. Després. "Discretization of hyperelasticity on unstructured mesh with a cell-centered Lagrangian scheme". In: *J. Comp. Phys.* (2010) (cit. on pp. 18, 24, 31).

[Kuc+12] M. Kucharik, R. Loubère, L. Bednárik, and R. Liska. "Enhancement of Lagrangian slide lines as a combined force and velocity boundary condition". In: *Comput. Fluids* (2012) (cit. on pp. 27, 28).

[Lah06] J.-L. Lahaie. "The Tera-10 System: Implementing the Number 1 Supercomputer in Europe". In: *Proceedings of the 2006 ACM/IEEE Conference on Supercomputing. SC '06*. Tampa, Florida: Association for Computing Machinery, 2006, 292–es. ISBN: 0769527000. DOI: 10.1145/1188455.1188760. URL: https://doi.org/10.1145/1188455.1188760 (cit. on p. 13).

[LW60] P.D. Lax and B. Wendroff. "Systems of conservation laws". In: *Commun. Pure Appl Math.* 13.2 (1960), pp. 217–237 (cit. on p. 9).

[Lef+18] E. Lefebvre, S. Bernard, C. Esnault, P. Gauthier, A. Grisollet, P. Hoch, L. Jacquet, G. Kluth, S. Laffite, S. Liberatore, I. Marmajou, P.-E. Masson-Laborde, O. Morice, and J.-L. Willien. "Development and validation of the TROLL radiation-hydrodynamics code for 3D hohlraum calculations". In: *Nuclear Fusion* 59.3 (2018), p. 032010 (cit. on pp. 13, 19, 34).

[Leh87] D. Lehr. *Misty picture event, test execution report*. Tech. rep. ADA283521, 1987 (cit. on p. 12).
[Leo91] B.P. Leonard. "The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection". In: *Computer Methods in Applied Mechanics and Engineering* 88.1 (1991), pp. 17–74. ISSN: 0045-7825 (cit. on p. 10).

[Lew97] R. Lewandowski. *Analyse Mathématique et Océanographie*. Recherches en Mathématiques Appliquées. Masson, 1997 (cit. on p. 6).

[LTW93] J.-L. Lions, R. Temam, and S. Wang. "Models for the coupled atmosphere and ocean". In: *Comput. Mech. Adv.* 1(1).120 (1993). (CAO I,II) (cit. on p. 6).

[LCF16] A. Llor, A. Claisse, and C. Fochesato. "Energy preservation and entropy in Lagrangian space- and time-staggered hydrodynamic schemes". In: *J. Comput. Phys.* 309 (2016), pp. 324–349 (cit. on p. 15).

[LC87] W. E. Lorensen and H. E. Cline. "Marching Cubes: A High Resolution 3D Surface Construction Algorithm". In: *Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques*. SIGGRAPH '87. New York, NY, USA: Association for Computing Machinery, 1987, pp. 163–169. ISBN: 0897912276 (cit. on p. 4).

[Lou+10] R. Loubère, P.-H. Maire, M. Shashkov, J. Breil, and S. Galera. "ReALE: A Reconnection-based Arbitrary-Lagrangian-Eulerian Method". In: *J. Comput. Phys.* (2010) (cit. on p. 35).

[LS05] R. Loubère and M. J. Shashkov. "A subcell remapping method on staggered polygonal grids for arbitrary-Lagrangian-Eulerian methods". In: *J. Comput. Phys.* 209 (2005), pp. 105–138 (cit. on p. 15).

[LF11] G. Luttwak and J. Falcovitz. "Slope limiting for vectors: a novel vector limiting algorithm". In: *Int. J. Numer. Meth. Fluids* 65 (2011), pp. 1365–1375 (cit. on pp. 19, 20).

[Mai+13] P. H. Maire, R. Abgrall, J. Breil, R. Loubère, and B. Rebourcet. "A nominally second-order cell-centered Lagrangian scheme for simulating elastic-plastic flows on two-dimensional unstructured grids". In: *J. Comput. Phys.* 235 (2013), pp. 626–665 (cit. on pp. 24, 31).

[Mai11] P.-H. Maire. *Contribution to the numerical modeling of inertial confinement fusion*. CEA-R-6260, Mar. 2011 (cit. on p. 18).

[Mai+07] P.-H. Maire, R. Abgrall, J. Breil, and J. Ovadia. "A cell-centered Lagrangian scheme for two-dimensional compressible flow problems". In: *SIAM J. Sci. Comput.* 29.4 (2007), pp. 1781–1824 (cit. on pp. 16, 18, 22, 35).

[MN08] P.-H. Maire and B. Nkonga. "Multi-scale Godunov-type method for cell-centered discrete Lagrangian hydrodynamics". In: *J. Comput. Phys.* 228 (2008), pp. 799–821 (cit. on pp. 16, 18).

[Mau01] B. Maury. "A Fat Boundary Method for the Poisson Problem in a Domain with Holes". In: *Journal of Scientific Computing* 16 (2001), pp. 319–339 (cit. on p. 4).

[Mau08] B. Maury. "Numerical Analysis of a Finite Element/Volume Penalty Method". In: *Partial Differential Equations: Modeling and Numerical Simulation*. Ed. by R. Glowinski and P. Neittaanmäki. Dordrecht: Springer Netherlands, 2008, pp. 167–185 (cit. on p. 4).

[Maz07] C. Mazeran. "Sur la structure mathématique et l'approximation numérique de l'hydrodynamique Lagrangienne bidimensionelle". PhD thesis. Université Bordeaux I, 2007 (cit. on pp. 13–16, 22, 23, 35).
[Mor+13] N. R. Morgan, M. A. Kenamond, D. E. Burton, T. C. Carney, and D. J. Ingraham. "An approach for treating contact surfaces in Lagrangian cell-centered hydrodynamics". In: *J. Comput. Phys.* 250 (2013), pp. 527–554 (cit. on p. 28).

[NR50] J. von Neumann and R. D. Richtmyer. "A method for the calculation of hydrodynamics shocks". In: *J. appl. Phys.* 21 (1950), pp. 232–237 (cit. on p. 14).

[Oll96a] C. F. Ollivier-Gooch. *High-order ENO schemes for unstructured meshes based on least-squares reconstruction*. Tech. rep. Mathematics and Computer Science Division, Argonne National Laboratory, 1996 (cit. on p. 20).

[Oll96b] C. F. Ollivier-Gooch. "Quasi-ENO schemes for unstructured meshes based on unlimited data-dependent least-squares reconstruction". In: *J. Comput. Phys.* (1996) (cit. on p. 20).

[Pat20] J. Patela. "Opérateur de diffusion vectoriel pour l'hydrodynamique lagrangienne turbulente sur maillage déformé". MA thesis. Université de Reims Champagne Ardenne, 2020 (cit. on p. 41).

[Pes72] C. Peskin. "Flow patterns around heart valves: A numerical method". In: *J. Comput. Phys.* 10.2 (1972), pp. 252–271. ISSN: 0021-9991 (cit. on p. 3).

[PLT92] O. Pironneau, J. Liou, and T. Tezduyar. "Characteristic-Galerkin and Galerkin Least Squares Space-Time Formulations for the Advection-Diffusion Equation with Time Dependent Domains". In: *Computer Method in Applied Mechanics and Engineering* 100 (1992), pp. 117–141 (cit. on p. 5).

[Ple23] A. Plessier. "Implicit semi-Lagrangian schemes for gas dynamics". PhD thesis. LJLL (UMR 7598) – Laboratoire Jacques-Louis Lions, 2023 (cit. on pp. 32, 42).

[PO01] A. Prosperetti and H.N. Oguz. "Physalis: A New o(N) Method for the Numerical Simulation of Disperse Systems: Potential Flow of Spheres". In: *J. Comput. Phys.* 167.1 (2001), pp. 196–216. ISSN: 0021-9991 (cit. on p. 4).

[RT99] T. Rossi and J. Toivanen. "A Parallel Fast Direct Solver for Block Tridiagonal Systems with Separable Matrices of Arbitrary Dimension". In: *SIAM J. Sci. Comput.* 20.5 (1999) (cit. on p. 3).

[SC02] A. J. Scannapieco and B. Cheng. "A multifluid interpenetration mix model". In: *Physics Letters A* 299 (2002), pp. 49–64 (cit. on p. 29).

[Sch97] J. Schöberl. "NETGEN An advancing front 2D/3D-mesh generator based on abstract rules". In: *Computing and Visualization in Science* 1 (1997), pp. 41–52 (cit. on p. 4).

[Si00] H. Si. "Tetgen: a quality tetrahedral mesh generator and a 3d Delaunay triangler. In: *Weierstrass Institute for Applied Analysis and Stochastics* (2000), pp. 150–171 (cit. on p. 4).

[Ste91] J. L. Steger. "The Chimera method of flow simulation". In: *Workshop on applied CFD*. Univ. of Tennessee Space Institute. 1991 (cit. on p. 1).

[Str68] G. Strang. "On the Construction and Comparison of Difference Schemes". In: *SIAM Journal on Numerical Analysis* 5.3 (1968), pp. 506–517 (cit. on p. 11).

[Swa77] P. N. Swarztrauber. "The methods of cyclic reduction, Fourier analysis and the FACR algorithm for the discrete solution of Poisson's equation on a rectangle". In: *Siam Review* 19.3 (1977), pp. 490–501 (cit. on p. 3).

[TT61] J. Trulio and K. Trigger. *Numerical solution of the one-dimensional Lagrangian hydrodynamic equations*. Tech. rep. UCRL-6267. Lawrence Radiation Laboratory, 1961 (cit. on p. 15).
[Vil12] F. Vilar. "Utilisation des méthodes de Galerkin discontinues pour la résolution de l'hydrodynamique Lagrangienne bi-dimensionnelle". PhD thesis. Université Bordeaux I, 2012 (cit. on p. 20).

[VMA14] F. Vilar, P.-H. Maire, and R. Abgrall. "A discontinuous Galerkin discretization for solving the two-dimensional gas dynamics equations written under total Lagrangian formulation on general unstructured grids". In: *J. Comput. Phys.* 276 (2014), pp. 188–234 (cit. on p. 20).

[Wag87] D. H. Wagner. "Equivalence of the Euler and Lagrangian equations of gas dynamics for weak solutions". In: *J. Diff. Eq.* 68 (1987), pp. 118–136 (cit. on p. 16).

[Wag96] D. H. Wagner. "Conservation laws, coordinate transformations, and differential forms". In: *Proceedings of the Fifth International Conference on Hyperbolic Problems Theory, Numerics, and Applications* (J. Glimm, MJ Graham, JW Grove, and BJ Plohr, eds.), *World Scientific Publishers, Singapore*. 1996, pp. 471–477 (cit. on p. 16).

[Wha96] P. P. Whalen. "Algebraic limitations on two-dimensional hydrodynamics simulations". In: *J. Comput. Phys.* 124 (1996), pp. 46–54 (cit. on p. 17).

[Wil64] M. L. Wilkins. "Calculation of Elastic-Plastic Flow". In: *Methods in Computational Physics*. Vol. 3. Academic Press, 1964, pp. 211–263 (cit. on pp. 15, 31).

[Yak07] D. Yakoubi. "Analyse et mise en œuvre de nouveaux algorithmes en méthodes spectrales". PhD thesis. Laboratoire Jacques-Louis Lions, Université Pierre et Marie Curie, Dec. 2007 (cit. on pp. 5, 6).