# Stable and Conservative Numerical Schemes for Hyperbolic Problems Using Non Convex Polygonal Meshes

Elise Le Mélédo

► To cite this version:

Elise Le Mélédo. Stable and Conservative Numerical Schemes for Hyperbolic Problems Using Non Convex Polygonal Meshes. Numerical Analysis [math.NA]. University of Zuerich, 2021. English.
NNT: tel-03917082

HAL Id: tel-03917082

https://theses.hal.science/tel-03917082

Submitted on 8 Jan 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.

<figure><img src="image_1.png" /></figure>
STABLE AND CONSERVATIVE NUMERICAL SCHEMES
FOR HYPERBOLIC PROBLEMS USING NON-CONVEX
POLYGONAL MESHES

Dissertation

zur

Erlangung der naturwissenschaftlichen Doktorwürde (Dr. sc. nat.)

vorgelegt der

Mathematisch-naturwissenschaftlichen Fakultät

der

Universität Zürich

von

Élise Le Mélédo

aus Frankreich

Promotionskommission

Prof. Dr. Rémi Abgrall (Vorsitz und Leiter der Dissertation)

Prof. Dr. Stefan Sauter

Prof. Dr. Jean Bertoin

Zürich, 2022
<figure><img src="image_2.png" /></figure>
To all who wander
through space
while being
mindful of time
Abstract

Hyperbolic systems of equations are often encountered in problems mod-
elling natural phenomenon, and are known to admit solutions developing sev-
eral types of discontinuities. The design of suitable numerical schemes is
therefore delicate, as they should to be able to handle the emergence of dis-
continuities in the solution, represent them at the discrete level, and grasp
their time-evolving spatial location while ensuring to the approximate solu-
tion a physical meaning. On the top of those already delicate considerations,
another difficulty emerges when considering the particular subclass of hyper-
bolic problems that are further conservative, that is whose time-space dy-
namic preserves physical quantities depending on the state variables. The
numerical scheme has therefore to ensure the associated conservation prin-
ciples at the discrete level, which appears to be often contradictory with the
retrieval of the discontinuity locations. Hence, although those concerns have
been wisely investigated, even recent numerical solvers are lacking precision,
favouring the scheme stability over the spatial accuracy and physical rele-
vance of the solution.

With the aim to draw a step further in the resolution of those issues, we
improved in this thesis existing high order schemes by designing correction
procedures that either ensure the physical relevance of the solution or recover
the discrete conservation principles. We particularly focused on two distinct
aspects of the problem, that are the entropy stability of Flux Reconstruction
schemes and the spatial resolution of interfaces in two-phase flows that pre-
serve the discrete mass conservation principle, both in the two-dimensional
case. Hence considering exclusively schemes that are based on unstructured
grids, we further endeavoured to widen the range of usable meshes to hybrid
polygonal ones, thus resolving conformity issues in the meshing of compli-
cated domains and paving the way towards mesh adaptation techniques.

More precisely, we first investigate in this thesis the entropy stability of
Flux Reconstruction schemes, for which even though an extensive literature
investigating the accuracy, robustness and linear stability exists, considera-
tions with respect to their entropy conservation are still lacking.
As this missing stone appears to be a crucial feature for ensuring the physical relevance of a discrete weak approximate solution to an hyperbolic problem, we explore in this direction by developing a Flux Reconstruction scheme that is entropy conservative and that can be constructed on any hybrid polygonal mesh. To this aim, we reformulate an existing Flux Reconstruction scheme into the Residual Distribution framework so as to apply there known entropy results. A set of conditions on the correction function that are sufficient for the original Flux Reconstruction scheme to be entropy conservative is then naturally derived, yielding to an entropy conservative Flux Reconstruction that is and geometrically flexible provided the existence of correction functions satisfying the obtained constraints for the considered cell geometry.

Hence, we develop in a second part of the thesis a theoretical framework in which the correction function can be defined and characterised upon the previously obtained constraints. In practice, we bridge the geometrical flexibility of the Virtual Elements framework with the boundary variational structure of the Raviart – Thomas spaces, yielding a wide class of $H(\text{div})$-conforming discretisations that can be defined on any polytope, and whose face-wise discretisation setting is independent of both the cell geometry and inner discretisation. Extending the already more general $H(\text{div})$-conforming VEM setting, their definitions are in addition easily amenable to the properties the approximated quantities are wished to fulfil, thus being in particular suitable for discretising the previously investigated correction functions.

In the last part of this thesis, we direct our investigations towards the motion of interfaces in incompressible two-phase flows, whose approximation raises tremendous interesting numerical issues. Indeed, while the total volume of each phase remains constant in time, the interface may develop heavy geometrical distortions and change topology across time, so that a small error in the interface location will have dreadful repercussions on the approximate state dynamic. The numerical schemes tacking two-phase flows problems should therefore provide an accurate description of the interface and its motion while still ensuring the discrete conservation principle, thus allowing to select the proper phase properties on each subdomain and provide a conservative discretisation of the system dynamic.

However, at the discrete level, the conservation property often contradicts with the determination of the interface spatial location. Hence, hybrid corrective method such as the THINC-LS have been developed, initially determining the interface location by advecting a level set before correcting it in each cell upon conservation constraints. While this approach indeed retrieves both a discrete conservation principle and an approximate location for the interface, this achievement is reached to the price of the interface global continuity. We have improved the correction procedure by introducing a new definition of a correction function, allowing to enforce the conservation constraint in any control volume, possibly non-convex, while preserving the continuity of each connected component of the interface.
Acknowledgments

First of all, I would like to express all my gratitude to my supervisor, Professor Rémi Abgrall, for his continuous support throughout all the project, his humanity, his enthusiasm in sharing his knowledge, his trust, and for extensive freedom and research latitude he has conferred to me. All in all, it allowed me to broaden my knowledge and thought me how to juggle between different perspectives and perceptions of the investigated mathematical concepts. I would also like to thank him for giving me the opportunity to participate in conferences and plan short stays at other laboratories. There, I would like to thank the Professor Xiao for having being keen to have me as a visiting student in his laboratory should the pandemic circumstances not had compromised the stay. I would like to thank the Professors Eric Darrigrand and Christophe Chalons that made me discover the field of numerical approximation to hyperbolic problems, and to the Professors Philippo Santambrogio and Tohru Yagi for their recommendation.

I would also like to thank Philipp Öffner for his friendship, support and the pleasant collaborative work. On similar lines, I would like to thank Barbara Re, Hendrik Ranocha and Davide Torlo for the exchange of ideas and passion. More generally, I would like to thank all the members of the working group for having created a wonderful pleasant and collaborative working environment, especially Jianfang Lin, for her tremendous and communicative energy, Paola Bacigaluppi for her wise advises and friendship, Wasilj Barsukov for his general curiosity that made any coffee break instructive, Maria Han Veiga for her consistently joyful mood, Roxana Curtici for her friendship, and Lorenzo Micalizzi for his dynamism and motivational attitude.

On a larger scale, I would like to thank the University of Zurich and the Mathematical institute for creating an uncomplicated and easy-going work environment, erasing the walls between the natural working groups. Complementarily, I would also like to thank all the colleagues who contributed to the general warm and caring atmosphere of the department, in particular Severin, Simon, Andres, Genta, Benjamin, Niklas, Artem, Jacopo, Francesco, Simone, Raul, Karan, Céline, and Marco.

The wonderful IT team is also to be warmly thanked for, especially Carsten, Dario and Benj, for their support, promptness, and any-problem-solving solutions.

Lastly, I would like to thank the Swiss National Foundation for granting the project No. 175784, "Solving advection dominated problems with high order schemes with polygonal meshes: application to compressible and incompressible flow problems", support without which this adventure could never have been possible, and the referees for their careful review of the outcome.

On a more minor note, I would like to thank Tchaikovsky, Sibelius and Grieg for their support in the writing and coding process, as well as Mahler, Bruckner and Shostakovich through difficult times.
<figure><img src="image_3.png" /></figure>
# Introduction

The approximation of solutions to mathematical problems plays a significant role in the visualisation of their features while providing insights to the very nature of the problem itself. In turn, a numerical method that has been designed for taking into account specific features of a given dynamics can then be applied to a large range of other problems where similar behaviours are encountered, thereby providing a general approximation framework. A typical class of problems whose formulation and solutions share common features is the class of problems modelling conservative physical behaviours, the involved principles yielding dynamical systems that have a comparable structure, as pointed out by Poincaré back in 1890 [1].

## 1.1 Problem of interest and general difficulties

The specific structure that we are interested in in this thesis is the class of time-dependent problems governed by a set of conservative equations, whose system is of hyperbolic nature and dynamics preserves the total quantity of each variable in time. Those problems usually deriving directly from physical principles, their numerical approximation find natural applications in a wide spectrum of areas, ranging from engineering and environmental modelling, as e.g. for alloy development [2], weather prediction [3], to finding the origin of a river pollution [4] or even to determine a best wild-fire-fighting strategy [5]. Those applications requiring accurate and reliable approximate solutions to the considered problems, there is a strong demand for accurate solvers. However, the intrinsic nature of the problems makes their development tremendously complicated.

Indeed, any approximation relies on a discretisation of the problem, where the solution is made available only at specific points in the time-space. The dynamic has therefore to be represented accurately from a finite set of so-called degrees of freedom values while satisfying a discrete analogue of the involved physical principles and guaranteeing the relevance of the solution. In par-
ticular, the designed numerical scheme should ensure that the conservation
principle holds at the discrete level for each conserved variable. Otherwise,
the quality of the solution would deteriorate over time, and any convergence
to the solution of the continuous problem would ultimately be lost. For the
solution to be reliable in the long run, the scheme must also be stable, that
is, robust to perturbations that may result from the initial data or from com-
putational errors. However, ensuring stability in the context of hyperbolic
equations is extremely delicate as discontinuities may appear in the solution,
even for smooth initial conditions.

The treatment of these discontinuities is equally of major importance. In-
deed, as they are moving at a finite speed in the spatial domain, their loca-
tion impact significantly the dynamics, making their accurate determination
crucial for the numerical solution to be representative of the comprehensive
dynamics behaviour. However, defining an accurate representation is very
delicate as those discontinuities are possibly of different natures, being either
preserved or weakened in time. Furthermore, as a consequence of the pres-
ence of discontinuities, the global regularity of the sought solution has to be
weakened, leading to a lack of uniqueness. Hence, at the discrete level the
numerical scheme must ensure that the approximate solutions converge to
a physically meaningful one. We therefore consider here an entropy-based
selection criterion by asking the numerical scheme to be entropy stable.

Lastly, it should also be mentioned that in practical applications, the study
domain can have a complex geometry. It is therefore in the interest of numer-
ical schemes to be geometrically flexible, so as to approximate the solution
with the same accuracy all over the computational domain. Moreover, with
the recent growth in the size and complexity of the investigated problems,
there is a strong demand to converge quickly to the solution, i.e. to obtain a
reliable solution without much spatial discretisation effort.

However, the current numerical methods that are dedicated to the approxi-
mation of solutions to conservative hyperbolic problems either severely lack
precision in determining the physically relevant spatial location of the dis-
continuities, or fail to ensure the conservation and stability principles, both
features being crucial to render a physically relevant approximated solution.
**1.2 Scope of the presented work**

In this thesis, we focus on designing high order schemes that provide a local
conservation principle to the retrieved solution while ensuring its homoge-
neous accuracy across the spatial domain and its physical relevance.

**Background choices** To achieve conservation down to the discrete level, we focus on developing gridded numerical methods. Indeed, there, the computational domain is covered with a finite number of cells, small tiles defining control volumes on which the solution's variables can be quantified, preserved, and where a local conservation can be directly enforced. Furthermore, as we aim the developed schemes to be robust with respect to the computational domain geometry and to be flexible with respect to its discretisation, we consider the particular subclass of hybrid unstructured meshes where the cells can be of various shapes and their layout organised in an irregular pattern.

In this regard, the use of non-convex polygonal meshes appears to particu-
larly relevant, as they allow more flexibility in the geometry description, more
robustness towards distortions, and facilitate the development of geometry
dependent techniques such as mesh adaptation techniques by automatically
including hanging nodes and facilitating local mesh coarsening.

Although this class of schemes has been widely investigated in the context
of elliptic problems [6, 7, 8, 9, 10], their use in the context of hyperbolic sys-
tems is relatively recent [11, 12, 13]. Indeed, as for any unstructured mesh,
there is a trade off between the order of the there constructed scheme, in-
dicating by how much the accuracy of the discrete solution improves when
the mesh size is decreased, and the ability of retrieving its possible disconti-
nuities while guaranteeing stability of the approximation [14]. In particular,
the presence of discontinuities yields high order schemes to require specific
treatments as filtering or limiting techniques to ensure stability [15, 16, 17],
and even to ensure that the obtained discrete solution has a physical meaning
[18, 19]. While those treatments have been widely investigated in simplicial
and quads meshes, the range of existing techniques for cells having a higher
number of faces is still limited, their design being much more involved both
geometrically and variationally.

**Considered scope of this thesis** We therefore aim here at improving the quality and robustness of the approximated solution to conservative hyperbolic problems by developing schemes that are both stable and entropy stable even across discontinuities, without sacrificing the conservation. Those discontinuities arising either in the field of a single variable or in the interaction of physical variables, there creating a *so-called* interface, our developments have taken place on two distinct levels.
## 1.3 Thesis contribution

Considering first the case where a discontinuity arises in the field of a single variable, we focused on the scheme stability and physical relevance of the retrieved solution, developing a Flux Reconstruction scheme that can be used on any arbitrary hybrid unstructured mesh whose cells are not necessarily convex and an associated theoretical framework to ascertain its well-definition.

We then focused on discontinuities arising in the interaction of physical variables by improving a correction procedure for mass-conservative in two-phase flows. More precisely, we were able to locate continuously the mixture's interface in two-phase flows while preserving the local mass quantities.

**Development of an entropy stable Flux Reconstruction scheme**

In the first part of the project, we achieved entropy stability by recasting a Flux Reconstruction scheme into the Residual Distribution framework. To this aim, we there used the conservation and stability results of the Residual Distributions schemes to develop an entropy stable Flux Reconstruction scheme in its Residual Distribution formulation, before deriving a determination system on the correction function allowing to convert back the expression of the obtained entropy stable scheme in its original Flux Reconstruction framework.

In order to ascertain the existence of such a correction function, we were led to develop a theoretical framework within which there exists at least one function solving the previously developed determination system. Namely, we constructed a class of $H(\text{div})$-conforming elements that can be set up on any non-convex polytopal cell and that benefit from the properties of the Raviart–Thomas elements on the boundary, using the Virtual Elements framework as a core. We thence obtained a very flexible class of elements that can be used in any dimension, on any possibly non-convex polytope, and whose variational setting can be easily adapted to a large variety of specific needs.

**Improvement of an interface capturing scheme**

In the second part of the project, we restricted ourselves to the two-dimensional case and considered incompressible two-phase flow problems. We particularly focused on retrieving the time-space dependent location of the mixture's interface so that both the conservation of the partial masses and the global continuity of the interface are obtained. To this aim, we focused on the class of level set based schemes and developed a continuous and conservative interface correction procedure by improving the already conservative THINC-LS technique [20].

The THINC-LS being based on the characterisation of a corrective constant, we changed its definition to a space-dependent correction function that preserves or enforces the continuity of the corrected interface location. Paying a particular attention to geometrical issues, we could design a correction
function that can be set up on any non-convex polygonal cell and whose support stays local to the control volume the interface location is corrected on. Hence, the improved correction procedure now allows to correct the initial interface location that has been achieved through any gridded scheme, enforcing the conservation of the local masses and preserving the global continuity of the interface. In addition, in case of an initially discontinuous interface, the global continuity can also be enforced up to a straightforward extension of the procedure.

## Related publications

* *On the connection between residual distribution schemes and flux reconstruction.* Rémi Abgrall, Élise Le Mélédo, Philipp Öffner.
*Arxiv preprint*, 2018, https://arxiv.org/abs/1807.01261.
*(Section 2.5 of Chapter 2)*

* *A class of finite dimensional spaces and $H(\text{div})$-conformal elements on general polytopes.* Rémi Abgrall, Élise Le Mélédo, Philipp Öffner.
*Technical report*, 2019, https://arxiv.org/abs/1907.08678.
*(Sections 3.2 and 3.3 of Chapter 3)*

* *General polytopal $H(\text{div})$-conformal finite elements and their discretisation spaces.* Rémi Abgrall, Élise Le Mélédo, Philipp Öffner.
*ESIAM: M2An*, 2021, https://doi.org/10.1051/m2an/2020048.
*(Section 3.3 of Chapter 3)*

* *A corrective conservative and continuity preserving Level Set approach for two phase flows.* Rémi Abgrall, Élise Le Mélédo.
*In preparation*
*(Section 4.3 of Chapter 4)*

## Related public implementation

* **LESARD: Level Set Algorithm for Residual Distribution**, open source and user-friendly code solving two-dimensional partial differential equation problems. Élise Le Mélédo
  Repository: https://git.math.uzh.ch/elemel/lesard
  Documentation: https://www.math.uzh.ch/pages/lesard/
## 1.4 Outline

All in all, we focus on developing stable, entropy stable, conservative, high order and geometrically flexible schemes for conservative hyperbolic systems. We start in Chapter 2 by introducing the principal features of hyperbolic problems before detailing a construction of an entropy stable Flux Reconstruction scheme. We then present in Chapter 3 the theoretical foundations required by the previously presented numerical scheme, and end in Chapter 4 by detailing the improvement of the THINC-LS correction procedure ensuring the space-time dependent resolution of a continuous and conservative interface location in two-phase flows.

## Thesis overview

<table><tbody><tr><td><b>1 Introduction</b></td><td><b>3</b></td></tr><tr><td>1.1 Problem of interest and general difficulties</td><td>3</td></tr><tr><td>1.2 Scope of the presented work</td><td>5</td></tr><tr><td>1.3 Thesis contribution</td><td>6</td></tr><tr><td>1.4 Outline</td><td>8</td></tr><tr><td><b>2 Preliminaries</b></td><td><b>9</b></td></tr><tr><td>2.1 Considered type of problems</td><td>11</td></tr><tr><td>2.2 A brief algebraic toolbox</td><td>39</td></tr><tr><td>2.3 A brief numerical analysis toolbox</td><td>50</td></tr><tr><td>2.4 Some classical numerical schemes</td><td>58</td></tr><tr><td>2.5 An entropy conservative scheme</td><td>83</td></tr><tr><td><b>3 H(div)-conforming elements on general polytopal meshes</b></td><td><b>97</b></td></tr><tr><td>3.1 Conforming discretisation strategies</td><td>99</td></tr><tr><td>3.2 Description of two conforming elements</td><td>107</td></tr><tr><td>3.3 A framework for arbitrary polytopes</td><td>121</td></tr><tr><td><b>4 A conservative and continuity preserving interface correction</b></td><td><b>161</b></td></tr><tr><td>4.1 Problem statement, modelling, and approximation strategy</td><td>163</td></tr><tr><td>4.2 Description of three fundamental methods</td><td>177</td></tr><tr><td>4.3 A continuity preserving procedure</td><td>184</td></tr><tr><td><b>Bibliography</b></td><td><b>iii</b></td></tr></tbody></table>
# Preliminaries

In order to ease the discussions of the two following chapters, let us present the fundamental notions involved in the description of hyperbolic problems and their approximation strategies. Thereby, we derive in Section 2.1 their general form and state their common properties, before presenting in Section 2.2 algebraic notions that are required for their approximation. Lastly, we summarise in Section 2.3 the basic principles of numerical analysis, and overview the Residual Distribution framework that is underlying in Chapter 4, and for which we developed specific theoretical grounds in Chapter 3.
# Chapter contents

<table><tbody><tr><td>2.1</td><td>Considered type of problems</td><td>11</td></tr><tr><td>2.1.1</td><td>Conservation principle</td><td>12</td></tr><tr><td>2.1.2</td><td>A founding example</td><td>20</td></tr><tr><td>2.1.3</td><td>Hyperbolicity and discontinuity propagation</td><td>27</td></tr><tr><td>2.1.4</td><td>Weak solutions</td><td>28</td></tr><tr><td>2.1.5</td><td>Solution profiles</td><td>29</td></tr><tr><td>2.1.6</td><td>Entropy solutions</td><td>36</td></tr><tr><td>2.2</td><td>A brief algebraic toolbox</td><td>40</td></tr><tr><td>2.2.1</td><td>Multivariate polynomial spaces</td><td>40</td></tr><tr><td>2.2.2</td><td>Polynomial representation and degrees of freedom</td><td>42</td></tr><tr><td>2.2.3</td><td>Basis tuning on degrees of freedom</td><td>45</td></tr><tr><td>2.2.4</td><td>Basis functions adjustment to geometry distortion</td><td>47</td></tr><tr><td>2.3</td><td>A brief numerical analysis toolbox</td><td>51</td></tr><tr><td>2.3.1</td><td>Basic concepts of discretisation</td><td>51</td></tr><tr><td>2.3.2</td><td>Essential properties of a numerical scheme</td><td>56</td></tr><tr><td>2.4</td><td>Some classical numerical schemes</td><td>59</td></tr><tr><td>2.4.1</td><td>Discontinuous Galerkin methods</td><td>59</td></tr><tr><td>2.4.2</td><td>Flux reconstruction schemes</td><td>61</td></tr><tr><td>2.4.3</td><td>Residual distribution framework</td><td>77</td></tr><tr><td>2.4.4</td><td>Methods strengths and drawbacks</td><td>84</td></tr><tr><td>2.5</td><td>An entropy conservative scheme</td><td>86</td></tr><tr><td>2.5.1</td><td>Reformulation of a Flux Reconstruction scheme</td><td>86</td></tr><tr><td>2.5.2</td><td>Discrete conservation of the recast scheme</td><td>88</td></tr><tr><td>2.5.3</td><td>Entropy conservation of the recast scheme</td><td>89</td></tr><tr><td>2.5.4</td><td>Flux Reconstruction nature of the corrected recast scheme</td><td>93</td></tr><tr><td>2.5.5</td><td>Summary of the construction</td><td>96</td></tr><tr><td>2.5.6</td><td>Conclusion and perspectives</td><td>97</td></tr></tbody></table>
## 2.1 Considered type of problems

The type of problems investigated in this thesis generally describe a physical phenomenon whose mathematical model consists in a set of so-called *conservation laws*. In this scope, the solution is sought on a given domain $\Omega \subset \mathbb{R}^d$, $d \in \mathbb{N}$ as a time-space dependent *state vector* $Y: \mathbb{R}_+ \times \Omega \to \mathbb{R}^q$, $q \in \mathbb{N}$ whose dynamic is ruled by the following system.

$$ \frac{\partial Y}{\partial t} + \nabla \cdot \mathcal{F}(Y) = 0 \qquad (2.1a) $$

$$ Y|_{\partial \Omega} = g \qquad (2.1b) $$

$$ Y(0, \cdot) \in \mathcal{D} \subset \mathbb{R}^q, \text{ given} \qquad (2.1c) $$

Each *state variable* contained in the vector $Y$ describes an aspect of the solution whose expression typically involves physical quantities, as for example the height, density, pressure, or direction-wise velocities. Modelled by (2.1a), the dynamic of the state variables is interdependent, formulated through a *flux vector* $\mathcal{F} = (\mathcal{F}_1, \dots, \mathcal{F}_d)^T$ that contains $d$ functions $\mathcal{F}_i: \mathcal{D} \to \mathbb{R}^q$, $i \in [1, d]$ whose definition domain $\mathcal{D}$ corresponds to the subset of $\mathbb{R}^q$ where the value of the state variable have a physical meaning. Each of the flux functions belongs to $C^1(\Omega)$ and contributes to the global dynamic through its variations with respect to the state along the $i^{th}$ spatial direction. This general model is completed by *boundary conditions* specifying the impact of the domain's boundary on the dynamic, here exemplified as Dirichlet with for second member a function $g: \partial \mathbb{R}_+ \times \Omega \to \mathbb{R}^q$, and possibly closure equations.

<figure><img src="image_4.png" /><figcaption>**Figure 2.1.** Example of a spatial (left) and variational (right) domain for $d = 2$, $q = 2$, and $\mathcal{D} = \mathbb{R}_+ \times \mathbb{R}$.</figcaption></figure>

Driven by such a dynamic, each scalar state variable sees its total amount conserved in time, up to possible inflows or outflows through the domain's boundary. The variable $Y$ is hence referred to as the vector of *conservative variables*, and the system (2.1a) is said to be written in a *conservative form*. The classification of (2.1a) as a *system of conservation laws* is therefore natural, each scalar flux function $\mathcal{F}_i$ being in general derived from physical principles where the controlled variable is subject to a *conservation principle*.
### 2.1.1 Conservation principle

As the name suggests, a conservation principle rules the time evolution of a variable representing a physical state that participates to the description of a physical phenomenon by preserving in time a measurable quantity that depends on it. Hence at the equilibrium, this measurable quantity may well take into account any interaction between several of the other state variables and involve source terms (see an illustration in Figure 2.2).

<figure><img src="image_5.png" /><figcaption>**Figure 2.2.** Diagram of external contributions to a given system that can be involved in the definition of a conservation principle.</figcaption></figure>

We derive here a general expression of the conservation principle in the continuous setting, highlighting the essence of the system (2.1a)'s form. The obtained formulation will also allow us to derive naturally discrete analogues in the later Sections 2.5 and 4.3, in accordance with their corresponding discretisation framework. There, the preservation of the discrete conservation principle by the designed numerical scheme is of crucial importance for the retrieved solution to be reliable, and will therefore partly drive the approximation choices.

<figure><img src="image_6.png" /><figcaption>**Figure 2.3.** Illustration of the continuous conservation principle in one dimension, representing the different quantity displacements and transfer that can be involved.</figcaption></figure>

The conservation principle being by nature closely linked to the notion of motion, we first derive the Reynolds transport theorem, from which any conservation principle follows when contextualised in the physical phenomenon of interest. The general form of the system (2.1a) is then naturally retrieved.
### 2.1.1.1 A representation of motion

Key point in the derivation of the Reynolds transport theorem, let us introduce the notion of motion by investigating how a given quantity $q$ is moved through time in the domain $\Omega$, that is, how its trajectory can be described.

**Trajectory representation** Assuming first that the scalar quantity $q$ is initially supported on a single point $x_0 \in \Omega$ and is passive, one can see its displacement as its transport through an underlying a time-space dependent velocity field

$$v: \mathbb{R}_+ \times \Omega \to \mathbb{R}^d \\ (t, x) \mapsto v(t, x). \tag{2.2}$$

Its time-dependent position can then be naturally retrieved through the advection of its initial location $x_0$. Thus, assuming further that $v \in C^1(\mathbb{R}_+ \times \Omega)$, and considering without loss of generality that the initial time associated to the initial location $x_0$ reads $t_0 = 0$, the time-dependent support of the quantity $q$ is always reduced to a single point, and one can retrieve its trajectory $t \mapsto \xi(t, x_0)$ as the solution to the Cauchy problem

$$\left\{ \begin{array}{l} \frac{d\xi}{dt} = v(t, \xi) \\ \xi(0, x_0) = x_0, \end{array} \right. \tag{2.3}$$

The trajectory of the quantity $q$ in the time-space therefore corresponds to the characteristic of the system (2.3) that is issued from the foot $x_0$.

<figure><img src="image_7.png" /><figcaption>**Figure 2.4.** Displacement of a scalar quantity and distortion of control volumes.</figcaption></figure>

**Subdomains distortion** Considering any infinite dimensional collection of initial points $\{x_0\}$ forming a smooth compact subdomain $\omega_0 \subset \Omega$ of non-null measure, one can define the time-dependent distorted subdomain that is associated to the dynamic (2.3) by

$$\omega(t) = \{x \in \mathbb{R}^d, \exists x_0 \in \omega_0, x = \xi(t, x_0)\},$$

defining the image of $\omega_0$ through the flow $v$ at the time $t$. One can therefore define the applications

$$\xi_t: \omega_0 \to \omega(t) \\ x \mapsto \xi(t, x), \tag{2.4}$$

representing the displacement of any scalar quantity $q$ from its initial location $x_0 \in \omega_0$ to its position at the time $t \in \mathbb{R}_+$ in the subspace $\omega(t)$.
**Bijectivity of the time-dependent position mapping** The trajectory being so parametrised, one can particularly notice that any application (2.4) is bijective for any $t$ at least in a neighbourhood of $t_0$, the velocity field being here assumed smooth.

**Proposition 2.1.1** *Local bijectivity of the position mapping*

There exists a time $t_{\max} \in \mathbb{R}_+ \cup \{+\infty\}$ such that for any time $t \in [0, t_{\max})$, the application

$$ \xi_t: \omega_0 \to \omega(t) \qquad (2.4) \\ x \mapsto \xi(t, x), $$

is a bijective mapping.

*Proof.* Let us define the Jacobian of the application (2.4) as $J_{\xi}: t \mapsto |\nabla J_{\xi t}|$, where $|\cdot|$ denotes the determinant. One can then investigate the bijectivity of (2.4) by checking for increasing times $t$ grows whether the matrix $J_{\xi}(t)$ is invertible or not.

There, providing directly an existence statement, we can observe that at the initial time $t_0 = 0$, the relation $\xi(0, x_0) = x_0$ holds by definition of (2.4). Hence ensuring that the matrix $J_{\xi}(0)$ is the identity, the invertibility of the application $x \mapsto \xi(0, x)$ follows. Thus, as in parallel the vector field $v$ transfers its assumed smoothness to $t \mapsto J_{\xi}(t)$ through (2.4), there exists necessarily a neighbourhood $[0, t_{\max})$, $t_{\max} \in \mathbb{R}_+$ of $t_0 = 0$ where for any $0 \le t < t_{\max}$, the Jacobian $J_{\xi t}$ is invertible, concluding the proof by ensuring there the bijectivity of the application (2.4).

Conferring as a by-product a $C^1(\omega_0)$-regularity to the application (2.4) in the local interval $[0, t_{\max})$ by the inverse mapping theorem, Proposition 2.1.1 allows us to derive a relation on the operator $J_{\xi}$ linking directly the time derivative of the trajectory's Jacobian with the underlying velocity field. Indeed, defining the operator

$$ J_{\xi}: [0, t_{\max}) \to C^1(\omega_0) \\ t \mapsto J_{\xi}(t) := |\nabla J_{\xi t}|, $$

one can express its time derivative

$$ \frac{dJ_{\xi}}{dt} = \frac{d}{dt} |\nabla J_{\xi t}| \qquad (2.5) $$

by plugging the definition of the trajectory (2.4). Simply yielding

$$ \frac{dJ_{\xi}}{dt} = \nabla \cdot (v(t, \xi(t, x))) J_{\xi}, \qquad (2.6) $$

this relation will become the cornerstone of the Reynolds transport theorem, allowing us to map the Eulerian and Lagrangian coordinate systems.
### 2.1.1.2 Reynolds transport theorem

Now that we could describe the motion of a scalar quantity $q$ taking values in a subset $I \subset \mathbb{R}$ and supported on a single point whose location is changing in time, let us describe the general motion of the total amount of pointwise quantities that are being transported in a time $t \in \mathbb{R}_+$ from an initial subspace $\omega_0$ to the time-dependent subspace $\omega(t)$. In other words, we investigate the time-evolution of the quantity

$$Q = \int_{\omega(t)} q(t, x) \, dx$$

by reformulating $\frac{\partial Q}{\partial t}$ in terms of the underlying velocity field $v \in C^1(\mathbb{R}_+ \times \Omega)$ and of the infinitesimal time evolution of the pointwise quantity $q: (t, x) \mapsto q(t, x)$ by using the Reynolds transport theorem.

**Theorem 2.1.1** *Reynolds transport theorem*

Let us consider a time interval $[t_0, t_{\max}) \subset \mathbb{R}_+$, a smooth initial subdomain $\omega_0 \subset \Omega$, and a scalar quantity $q: \mathbb{R}_+ \times \Omega \to I$ taking values in $I \subset \mathbb{R}$ that is transported in time across the domain $\Omega$ by an underlying velocity field $v \in C^1(\mathbb{R}_+ \times \Omega)$. Denoting then by $\omega(t)$ the image of $\omega_0$ through the flow $v$ at the time $t$ and describing the trajectory of $q$ from its foot $x \in \omega_0$ by the application $(t, x) \mapsto \xi(t, x)$ that matches a characteristic of the system (2.3), one can derive the relation

$$\frac{d}{dt} \left( \int_{\omega(t)} q(t, x) \, dx \right) = \int_{\omega(t)} \frac{\partial q}{\partial t}(t, x) \, dx + \int_{\partial \omega(t)} q(t, x) (v(t, x) \cdot n) \, d\gamma \quad (2.7)$$

for any $t \in [t_0, t_{\max})$, $t_{\max}$ being the maximum time for which the characteristics (2.4) are bijective, and where $\gamma$ represents the surface measure.

<figure><img src="image_8.png" /><figcaption>**Figure 2.5.** Illustration of the Reynolds theorem.</figcaption></figure>
**Proof of Theorem 2.1.1**

The relation (2.7) is a natural consequence of the expansion of the term $\frac{d}{dt}(\int_{\omega(t)} q(t, x) dx)$ by means of the relation (2.6). Indeed, as we know by assumption that the initial subdomain $\omega_0$ is smooth and that the velocity field $v$ belongs to $C^1(\Omega)$, it comes from the above developments that the subdomain $\omega(t)$ is equally smooth for any $t \in [t_0, t_{\max})$. We can therefore express the term $\frac{d}{dt}(\int_{\omega(t)} q(t, x) dx)$ in the coordinate system of the initial subdomain $\omega_0$ by writing

$$ \frac{d}{dt} \left( \int_{\omega(t)} q(t, x) \, dx \right) = \frac{d}{dt} \left( \int_{\omega_0} q(t, \xi) J_{\xi}(t, \xi) \, d\xi \right), \quad (2.8) $$

where we dropped the dependency of $\xi$ on $t$ and $x$ for the sake of legibility. Considering here only the times $t_0 \le t < t_{\max}$, the quantity $(t, x) \mapsto q(t, x)$ is smooth enough and one can further apply the Leibnitz rule. Reading

$$ \frac{d}{dt} \left( \int_{\omega(t)} q(t, x) \, dx \right) = \int_{\omega_0} \frac{d}{dt} \left( q(t, \xi) J_{\xi}(t, \xi) \right) \, d\xi, \quad (2.9) $$

one can then develop the integrand to obtain

$$ \frac{d}{dt} \left( \int_{\omega(t)} q(t, x) \, dx \right) = \int_{\omega_0} \left( \frac{\partial q}{\partial t}(t, \xi) J_{\xi}(t, \xi) + q(t, \xi) \frac{\partial J_{\xi}}{\partial t}(t, \xi) \right) \, d\xi, \quad (2.10) $$

reducing to

$$ \frac{d}{dt} \left( \int_{\omega(t)} q(t, x) \, dx \right) = \int_{\omega_0} J_{\xi}(t, \xi) \left( \frac{\partial q}{\partial t}(t, \xi) + q(t, \xi) \nabla \cdot v(t, \xi) \right) \, d\xi. \quad (2.11) $$

Using the bijectivity of $\xi_t$ in $[t_0, t_{\max})$ to transfer the expression back in the original coordinate system, we then obtain

$$ \frac{d}{dt} \left( \int_{\omega(t)} q(t, x) \, dx \right) = \int_{\omega(t)} \left( \frac{Dq}{Dt}(t, x) + q(t, x) \nabla \cdot v(t, x) \right) \, dx, \quad (2.12) $$

where $q \mapsto \frac{Dq}{Dt}$ represents here the Lagrangian derivative of $q$ with respect to the time variable. Recalling then that $\frac{Dq}{Dt} = \frac{\partial q}{\partial t} + v \nabla \cdot q$, the above term further reduces to the relation

$$ \frac{d}{dt} \left( \int_{\omega(t)} q(t, x) \, dx \right) = \int_{\omega(t)} \frac{\partial q}{\partial t}(t, x) \, dx + \int_{\partial\omega(t)} (q(t, x) v(t, x) \cdot n) \, d\gamma, \quad (2.13) $$

by applying the divergence theorem, concluding the proof.

Investigating then the evolution of the total quantity $Q$ over time, the application of this theorem in the context of an investigated physical phenomenon directly yields to the expression of a conservation principle.
### 2.1.1.3 Conservation principle

As mentioned at the beginning of Section 2.1.1, any conservation principle rules the time evolution of a quantity representing a physical state that participates to the description of a physical phenomenon by preserving in time a measurable quantity that depends on it. In practice, this quantity measure is supported on an arbitrary compact subset $\omega \subset \Omega$ called control volume, and takes into account the possible source terms and possible physical interactions at the control volume's boundaries.

Hence woven in the considered physical context, asking its value to be preserved in time contributes to defining the local dynamic of the state variables. Furthermore, given a discretisation framework, the value of its discrete analogue contributes to the assessment of the approximate solution reliability. We therefore derive here a general definition, driving the later numerical approximation choices and highlighting the essence of the system (2.1a)'s form.

**Continuous conservation principle** In order to derive a general expression of a continuous conservation principle, we first investigate the rate of change in the amount of a quantity $q: [0, t_{\max}) \times \Omega \to \mathbb{R}^{\tilde{q}}$, $\tilde{q} \ge 1$, over a given compact control volume $\omega \subset \Omega$ through time until $t_{\max} \in \mathbb{R}_+ \setminus \{0\}$. To this aim, we start by observing that in the context of Section 2.1, the quantity $q$ can be expressed in terms of the state variables $(Y_1, \dots, Y_q) \in \mathcal{D}$ and is therefore necessarily integrable under the assumptions of Sections 2.1.1.1 and 2.1.1.2. Hence, its total amount over $\omega$ can be quantified, and its evolution over time investigated through the variations of the application

$$Q: [0, t_{\max}) \to \mathbb{R}^{\tilde{q}} \\ t \mapsto \int_{\omega} q \, dx.$$

In particular, when considering the dynamics of the investigated physical phenomenon on $\omega$ as an isolated system, the total quantity $Q \in \mathbb{R}^{\tilde{q}}$ of $q$ over $\omega$ is preserved in time if and only if

$$\frac{dQ}{dt} = 0, \qquad (2.14)$$

regardless the possible changes in the spatial distribution of $q$ over $\omega$. Going beyond the isolated system and considering $\omega$ as part of the whole domain $\Omega$ and seeing the local dynamic ruling the change in the spatial distribution of the quantity $q$ over $\omega$ as a consequence of the global dynamic modelled on $\Omega$, one has further to take into account the impact of neighbouring behaviours across $\partial\omega$ and possible source terms supported on $\omega$. Hence, contextualising the local preservation of the total quantity amount (2.14) in the physical dynamic investigated globally over $\Omega$ yields the quantity's amount over $\omega$ to fulfil
the local relation

$$ \frac{dQ}{dt} = -S_I - S_B, \qquad (2.15) $$

where $S_I$ represents the total contribution of a given source term $S \in C^1([0, t_{\max}) \times \omega, \mathbb{R}^{\bar{q}})$, i.e.

$$ S_I = \int_{\omega} S(t, x) \, dx, $$

and the term $S_B$ represents the total contribution of neighbouring interactions. There, as here we consider physical systems of equations that model physical behaviours, those interactions comes from principles that are invariant by translation and rotation. The boundary term $S_B$ can be therefore be expressed as

$$ S_B = \int_{\partial\omega} \Theta(t, x, Y) \cdot n \, d\gamma $$

where $\gamma$ represents the boundary measure and for a term $\Theta \in (C^\infty([0, t_{\max}) \times \partial\omega \times \mathcal{D}, \mathbb{R}^{\bar{q}}))^d$ describing the time-space dependent quantity transfer through the boundary. All in all, the local relation (2.15) therefore reads

$$ \frac{d}{dt} \left( \int_{\omega} q(t, x) \, dx \right) = - \int_{\omega} S(t, x) \, dx - \int_{\partial\omega} \Theta(t, x, Y) \cdot n \, d\gamma, \qquad (2.16) $$

from which we can apply the Reynolds transport theorem on the left-hand side, the control volume $\omega$ being arbitrary. Namely, it comes

$$ \frac{d}{dt} \left( \int_{\omega} q(t, x) \, dx \right) = \int_{\omega} \frac{\partial q}{\partial t}(t, x) \, dx + \int_{\partial\omega} (q(t, x) v(t, x) \cdot n) \, d\gamma. $$

Hence, up to assuming $\omega$ to be of a regular shape and $qv \in C^\infty(R_+ \times \Omega)$, on can then the Gauss theorem to obtain the formulation

$$ \frac{d}{dt} \left( \int_{\omega} q(t, x) \, dx \right) = \int_{\omega} \frac{\partial q}{\partial t}(t, x) \, dx + \int_{\omega} \nabla \cdot (q \otimes v) \, dx, \qquad (2.17) $$

where as $v$ belongs to $\mathbb{R}^d$, and where $\otimes$ denotes the tensor product that reduces to $q \otimes v: (q, v) \mapsto qv$ when $q: [0, t_{\max}) \times \Omega \to \mathbb{R}$ and to the projection operator $q \otimes v \cdot n = (b \cdot n) a$ when $q: [0, t_{\max}) \times \Omega \to \mathbb{R}^d$. Hence, combining the relation (2.16) with the reformulation (2.17), we retrieve

$$ \int_{\omega} \frac{\partial q}{\partial t}(t, x) \, dx + \int_{\omega} \nabla \cdot (q \otimes v) \, dx = - \int_{\omega} S(t, x) \, dx - \int_{\partial\omega} \Theta(t, x, Y) \cdot n \, d\gamma. $$

Lastly reformulating the boundary term as

$$ \int_{\partial\omega} \Theta \cdot n \, d\gamma = \int_{\omega} \nabla \cdot \Theta \, dx $$

for a given $s: [0, t_{\max}) \times \omega \times \mathcal{D} \to (\mathbb{R}^{\bar{q}})^d$, we can write

$$ \int_{\omega} \frac{\partial q}{\partial t}(t, x) \, dx + \int_{\omega} \nabla \cdot (q \otimes v + \Theta) \, dx = - \int_{\omega} S(t, x) \, dx, $$

providing the general expression of the *local conservation*, and allowing us to give a definition of the general *continuous conservation principle*.
**Definition 2.1.1** *Continuous conservation principle*

Let us consider the time-space behaviour of a scalar quantity $q \in C^1([0, t_{\max}) \times \Omega, \mathbb{R}^{\tilde{q}})$ whose time-dependent spatial distribution is ruled by its initial state and a velocity vector $v \in C^1(\mathbb{R}_+ \times \Omega)$. There, if the quantity $q$ fulfils at any time $t \in [0, t_{\max})$ the relation

$$ \int_{\omega} \frac{\partial q}{\partial t} dx + \int_{\omega} \nabla \cdot (v \otimes q + \Theta) dx = - \int_{\omega} S(t, x) dx, \quad (2.18) $$

for given functions $\Theta \in (C^1([0, t_{\max}) \times \mathcal{D}))^d$ and $S \in C^1([0, t_{\max}) \times \omega, \mathbb{R}^{\tilde{q}})$, and for any control volume $\omega \subset \Omega$, then the quantity $q$ is said to be conserved. In addition, if $q \in Y$, $Y$ being the state vector describing the dynamic of interest, the quantity $q$ is said to be a *conservative variable*, and the term $\nabla \cdot \Theta$ represents the pointwise contribution of the modelled physics to the investigated quantity.

While we derived in Definition 2.20 a rather general expression of the conservation principle, the source term $S$ will appear to be identically vanishing in all the applications considered in this thesis. When referring to the conservation principle, we will then understand it in its restricted form

$$ \int_{\omega} \frac{\partial q}{\partial t}(t, x) dx + \int_{\omega} \nabla \cdot (q \otimes v + \Theta) dx = 0, \quad (2.19) $$

from which continuity equation immediately follow.

**Continuity equation** The control volume $\omega$ being arbitrary, the relation (2.19) is equivalent to asking the quantity $q$ to verify the Eulerian formulation

$$ \frac{\partial q}{\partial t} + \nabla \cdot (v \otimes q + \Theta) = 0. \quad (2.20) $$

of the conservation principle. Consequently, the conservation relation (2.19) can be directly enforced by driving the dynamic of $q$ through the differential form (2.20), in which the functional term $v \otimes q + \Theta$ is referred to as the *flux*.

**Systems of conservation laws** When a physical phenomenon is modelled by a collection of state variables that are each governed by an equation of the type (2.20), the comprehensive model can be written as a dynamical system of the form (2.1a). Hence, the time-space behaviour of each variable is entirely ruled by the its interaction with the other conservative variables, whose nature is driven from physical principle and amount is determined in accordance with the coexisting conservation constraints.

The general system structure being so predetermined, common behaviours can be observed in solutions, particularly relating to the regularity of the conservative variables. Let us provide an overview of those main features by considering the fundamental example of the Euler equations of gas dynamics.
## 2.1.2 A founding example

A fundamental example of problems having for structure (2.1) is given by the system of Euler equations of gas dynamics, whose study highlights the typical properties and behaviours that are encountered in their solutions. By definition, the set of Euler equations governs the state of a fluid modelled by its density $ρ$, pressure $p$ and velocity $v$ accordingly to the three physical conservation principles of mass, momentum and total energy $E$, neglecting in particular heat transfer phenomenon. We start here by deriving the system of Euler equations, applying in Section 2.1.2.1 the results of Section 2.1.1.3 to those three physical laws. We then present its classical formulation in Section 2.1.2.2 and derive its main features in Section 2.1.2.3, motivating the more general later Sections 2.1.3 to 2.1.6.

### 2.1.2.1 Derivation of three fundamental conservation laws

Let us here restrict ourselves to the scope of the Euler equations of gas dynamic and derive the expressions of the three conservation law of mass, momentum and energy when applied to a mixture consisting of fluids or gases. More precisely, we represent a given mixture by the time-space dependent values of $ρ$, $ρv$, and $ρE$, related to its density $ρ$, velocity $v$, energy $E$, and pressure $p$. Considered as state variables, those three quantities are defined over a domain $Ω ⊂ ℝ^d$, $d ∈ N$ and have a physical meaning whenever their values are contained in $ℝ_+$, $ℝ^d$ and $ℝ$, respectively.

**Mass conservation** Let us first derive the expression of the mass conservation principle, whose literal formulation reads as follows.

> Any system that is closed to any transfer of matter or energy sees its quantity of mass remaining constant over time.

Up to assuming that the investigated mixture dynamic can be represented by a closed system, translating into the absence of mass source term strictly within the domain $Ω$, one can directly translate this physical rule in a continuity equation of the form (2.20). Indeed, knowing that the total mass of the investigated quantity can be represented through the distribution of its density over the domain $Ω$ by the time-dependent expression $M := ∫_Ω ρ dx$, the total variation of mass can then be investigated by applying the Reynolds transport theorem, yielding

$$ \frac{dM}{dt} = \int_{\omega} \frac{\partial \rho}{\partial t} dx + \int_{\partial \omega} (\rho v) \cdot n dx. \qquad (2.21) $$

for any $ω ⊂ Ω$ compact control volume. Applying there the divergence theorem,
we retrieve that

$$\int_{\partial\omega} (\rho v) \cdot n d\gamma = \int_{\Omega} \nabla \cdot (\rho v) dx,$$

for $\gamma$ the surface measure, and as the quantity's mass being conserved through time if and only if $\frac{dM}{dt} = 0$, one can retrieve the mass conservation principle

$$\int_{\omega} \frac{\partial \rho}{\partial t} dx + \int_{\omega} \nabla \cdot (\rho v) dx = 0. \qquad (2.22)$$

As discussed in Section 2.1.1.3, one can then derive the Eulerian formulation

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho v) = 0, \qquad (2.23)$$

providing a first continuity equation contributing to the description of the considered mixture's dynamic.

**Note** Considering $\omega = \Omega$ makes the term $\int_{\partial\Omega} (\rho v) \cdot n dx$ emerge, highlighting the contribution of the net mass income from the boundary to the dynamic and vanishing whenever periodic or wall boundary conditions are considered. ▲

**Conservation of momentum** The expression of the momentum conservation principle is also known as the Newton's second law of motion and can be literal expressed as follows.

The rate of change of a system momentum equals the net sum of all forces acting on it.

Therefore, up to assuming that there exists only pressure forces acting on the mixtures, one can derive the momentum conservation principle in the form of (2.20) similarly as for the mass conservation principle. Indeed, the momentum of the system being defined over any compact control volume $\omega \subset \Omega$ by $\int_{\omega} \rho v dx$, the total contribution of the pressure forces to its local dynamic on $\omega$ can be expressed as a boundary term. Namely, remembering that the velocity vector $v$ belongs to $\mathbb{R}^d$, the conservation principle translates to

$$\frac{d}{dt} \left( \int_{\omega} \rho v \, dx \right) = - \int_{\partial\omega} p n \, d\gamma,$$

guaranteeing an equilibrium relation to the momentum. One can then apply again the Reynolds transport theorem to retrieve the relation

$$\int_{\omega} \frac{\partial}{\partial t} (\rho v) dx + \int_{\partial\omega} \rho v (v \cdot n) d\gamma + \int_{\partial\omega} p n d\gamma = 0, \qquad (2.24)$$

that after having applied the divergence theorem reduces to

$$\int_{\omega} \frac{\partial}{\partial t} (\rho v) dx + \int_{\omega} \nabla \cdot (\rho v \otimes v + p) dx = 0. \qquad (2.25)$$
The control volume $\omega$ being arbitrary and the velocity $v$ being a vector of dimension $d$, this last equation is then equivalent to the Eulerian formulation

$$ \frac{\partial}{\partial t}(\rho v) + \nabla \cdot (\rho v \otimes v + p \operatorname{Id}) = 0, \qquad (2.26) $$

where the term $\operatorname{Id}: \mathbb{R}^d \to \mathbb{R}^d$ denotes the identity. Hence, (2.26) defines a second continuous equation contributing to the description of the considered mixture's dynamic.

**Conservation of energy** A third physical principle involved in the considered description of the mixture dynamic relates to the preservation in time of the total amount of energy that the system has to self-maintain and evolve.

The total energy of any closed system remains constant.

In the considered setting, only two types of energies are involved; the *internal energy* $e$ of the system, used for self-maintenance, and the energy used by the system to set the mixture in motion, called *kinetic* and given by $\kappa = \frac{M}{2}v^2$ where $M$ represents the total system's mass. Hence, the total energy $E$ of the system simply writes

$$ E = e + \kappa = e + \frac{M}{2} \|v\|_2^2, $$

which, when considered over any compact control volume $\omega \subset \Omega$, can be expressed in terms of the mixture's density, i.e.

$$ \int_{\omega} \rho E \, dx = \int_{\omega} \rho e \, dx + \int_{\omega} \frac{\rho}{2} v \cdot v \, dx. $$

Considering then the energy contribution through the control volume boundary, the conservation relation then reads

$$ \frac{d}{dt} \left( \int_{\omega} \rho E \, dx \right) = - \int_{\partial \omega} p v \cdot n \, d\gamma, \qquad (2.27) $$

and by using once again the Reynolds transport theorem we then retrieve

$$ \int_{\omega} \frac{\partial(\rho E)}{\partial t} \, dx + \int_{\partial\omega} (\rho E) v \cdot n \, d\gamma = - \int_{\partial\omega} p v \cdot n \, d\gamma. \qquad (2.28) $$

As this relation holds for any compact control volume $\omega \subset \Omega$, we can derive similarly as in Section 2.1.1.3 the Eulerian formulation

$$ \frac{\partial(\rho E)}{\partial t} + \nabla \cdot ((\rho E + p) v) = 0, \qquad (2.29) $$

providing the last continuity equation defining the system of Euler equations.
### 2.1.2.2 The system of Euler equations

Having derived in Section 2.1.2.1 the fundamental continuity equations resulting from the three conservation law of mass, momentum and energy, let us now combine their dynamics to present the system of Euler equations.

**System dynamics** Recalling that the mixture is here modelled by its density $\rho \in \mathbb{R}_+$, pressure $p \in \mathbb{R}$, velocity $v := (v_1, \dots, v_d)^T \in \mathbb{R}^d$, total energy $E \in \mathbb{R}$, and evolves in time accordingly to the three continuum equations (2.23), (2.26) and (2.29), its comprehensive dynamics then reduces to

$$ 
\left\{
\begin{aligned}
& \frac{\partial}{\partial t}(\rho) + \sum_{i=1}^{d} \frac{\partial}{\partial x_i}(\rho v_i) = 0 \\
& \frac{\partial}{\partial t}(\rho u_i) + \sum_{i=1}^{d} \frac{\partial}{\partial x_i}(\rho v_i v_j + p \delta_{ij}) = 0 \\
& \frac{\partial}{\partial t}(\rho E) + \sum_{i=1}^{d} \frac{\partial}{\partial x_i}(\rho E + p v_i) = 0
\end{aligned}
\right.
\qquad (2.30)
$$

referred to as the *Euler equations*. Following [21], those equations can be reformulated in the form (2.1a) by defining the state vector of conservative variables as $Y = (\rho, \rho v, \rho E)$, where $v \in \mathbb{R}^d$ and where $\mathcal{F} = (\mathcal{F}_1(Y), \dots, \mathcal{F}_d(Y))^T$ is the flux vector containing the $d$ flux functions

$$ \mathcal{F}_i: \mathcal{D} \to \mathbb{R}^q \\ Y \mapsto \begin{pmatrix} \rho v_i \\ u_i v + p e_i \\ \rho v_i (E + p) \end{pmatrix}, \quad i \in [1, d], \qquad (2.31) $$

with $q = d+2$, and where $e_i$, $i \in [1, d]$ represents the $i^{th}$ vector of the canonical basis of $\mathbb{R}^d$, the term $v_i v + p e_i$ alone taking values in $\mathbb{R}^d$. Furthermore, please note that as the domain $\mathcal{D} \subset \mathbb{R}^q$ corresponds to the subspace of $\mathbb{R}^q$ where the value of the state variable have a physical meaning, reducing here to $\mathcal{D} = \mathbb{R}_+ \times \mathbb{R}^d \times \mathbb{R}$, the pressure to stay positive and the velocity finite.

**System closure** The system (2.30) is however not self-sufficient to describe completely the mixture behaviour, and has thus to be closed by linking the pressure $p$ with the energy through an Equation of State that takes into account the fluid properties. This relationship can be modelled e.g. by one of the three assumptions on the fluid nature reported in the below table, where $e$ denotes the specific internal energy of the fluid given by $E = \rho e + \frac{1}{2}\rho \|v\|_2^2$, $\gamma$ its specific heat ratio, $p_\infty$ its stiffness parameter, and $\rho_0$ its reference density.

<table><thead><tr><th>Ideal gas</th><th>Stiff gas</th><th>Mie-Grünesien</th></tr></thead><tbody><tr><td>ρe = <sup>p</sup>&frasl;<sub>γ-1</sub></td><td>ρe = <sup>p+γp<sub>∞</sub></sup>&frasl;<sub>γ-1</sub></td><td>p = p<sub>∞</sub>(ρ) + γ(ρe - ρ<sub>0</sub>(ρ))</td></tr></tbody></table>
### 2.1.2.3 Properties of the system

This system of Euler Equations is particularly representative of the problems of type (2.1) as it contains all the major properties common to such systems. Let us explore some of them in the following paragraph, particularly emphasising on its symmetrical and hyperbolic nature, and consequently on the developments of pathological behaviours in the solutions.

**A hidden conservation law: the system's entropy** The structure of the system hides another intrinsic conservation law, ensuring that the energy used by the system to evolve the solution matches exactly the one that is available, i.e. not gaining extra energy nor consuming its internal energy. Let us derive its actual expression by denoting $\iota := \rho e$ the internal energy and deriving the relation $E = \rho e + \frac{1}{2}\rho \|v\|_2^2$ with respect to the time, i.e.

$$ \frac{\partial \iota}{\partial t} = \frac{\|v\|^2}{2} \frac{\partial \rho}{\partial t} - v \cdot \frac{\partial}{\partial t} (\rho v) + \frac{\partial E}{\partial t}. $$

We can here observe that the relation decomposes on the vector $\tilde{Y} := (\rho, v, E)^T$. Thus, up to a straightforward change of variables in (2.31), we can multiply the set of equations by the vector $U = (\frac{\|v\|^2}{2}, -v, 1)$ and obtain the relation

$$ \left( \frac{\|v\|_2^2}{2}, -v, 1 \right) \left( \frac{\partial \tilde{Y}}{\partial t} \right) + \sum_{i=1}^{d} \left( \frac{\|v\|_2^2}{2}, -v, 1 \right) \left( \mathcal{F}_i(\tilde{Y}) \right) = 0, \quad (2.32) $$

reducing after algebraic manipulations to

$$ \frac{\partial \iota}{\partial t} + v \cdot \nabla \iota + (p + \iota) \nabla \cdot v = \frac{\partial e}{\partial t} + v \cdot \nabla e + p \nabla \cdot v = 0, $$

where we have used that $\iota = \rho e$. Introducing then the Gibbs relation

$$ T \frac{\partial s}{\partial t} = \frac{\partial e}{\partial t} - \frac{p}{\rho^2} \frac{\partial \rho}{\partial t} $$

where $s$ is the specific entropy depending only on the variable $t$ and $T$ is the absolute temperature, we finally obtain

$$ \begin{aligned} T \left( \frac{\partial s}{\partial t} + v \cdot \nabla s \right) &= \frac{\partial e}{\partial t} + v \cdot \nabla e - \frac{p}{\rho} \left( \frac{\partial \rho}{\partial t} + v \cdot \nabla \rho \right) \\ &= -\frac{p}{\rho} \left( \rho \nabla \cdot v + \frac{\partial \rho}{\partial t} + v \cdot \nabla \rho \right) = 0. \end{aligned} $$

Thus, knowing from thermodynamic considerations that $T \in \mathbb{R}_+$ and defining the function

$$ \begin{aligned} S: \mathcal{D} &\to \mathbb{R} \\ Y &\mapsto -\rho s(v), \end{aligned} $$

it comes the relation

$$ \frac{\partial S}{\partial t} + \nabla \cdot (S v) = 0. \quad (2.33) $$
The function $S$ being a convex function of $u$, it can be seen as an entropy flux. Hence, the relation (2.33) defines an entropy conservation law. Furthermore, we can now go back to the relation (2.32) and derive the equivalence

$$ U \frac{\partial \tilde{Y}}{\partial t} + \sum_{i=1}^{d} U \frac{\partial \mathcal{F}_i}{\partial x_i}(\tilde{Y}) \iff \frac{\partial S}{\partial t} + \sum_{i=1}^{d} \frac{\partial}{\partial x_i} (v_i S) = 0. \quad (2.34) $$

Furthermore, as the Euler equations are invariant by rotation we can define an entropy flux by the vector valued function $\mathcal{G}: Y \to S u$ and observe that

$$ U = (\nabla_{\tilde{Y}} S)^T \quad \text{and} \quad \nabla_{\tilde{Y}} \mathcal{G}_i = U \nabla_{\tilde{Y}} \mathcal{F}_i \quad \forall i \in [1, d], \quad (2.35) $$

relation defining the couple $(S, \mathcal{G})$ as an entropy function and flux pair. Furthermore, besides providing an characterisation of the system's energy, this last relation also allows us to observe the symmetry of the system (2.30).

**System symmetry** Let us observe the symmetry of the system (2.30). We there start by remembering that the functions $S$ and $\mathcal{F}$ are assumed to be smooth, and that the vector $U$ is a smooth bijective function of the vector variable $Y$. Thus, assuming for now that the solution $Y$ is smooth, we can apply the chain rule and obtain from (2.35)

$$ \left\{ \begin{aligned} U \frac{\partial \tilde{Y}}{\partial t} &= \frac{\partial}{\partial t} (U \tilde{Y}) - \left(\frac{\partial U}{\partial t}\right) \tilde{Y} = \frac{\partial S}{\partial t} \\ U \frac{\partial \mathcal{F}}{\partial t} &= \frac{\partial}{\partial t} (U \mathcal{F}) - \left(\frac{\partial V}{\partial t}\right) \mathcal{F} = \frac{\partial \mathcal{G}}{\partial t}. \end{aligned} \right. \quad (2.36) $$

Hence, there exists $\vartheta: \mathcal{D} \to \mathbb{R}$ function of $\tilde{Y}$, and a vector $\Upsilon = (\Upsilon_1, \dots, \Upsilon_d)$ of $d$ functions $\Upsilon_i: \mathcal{D} \to \mathbb{R}^q$ such that $S = U \tilde{Y} - \vartheta$ and $\mathcal{G} = U \mathcal{F} - \Upsilon$. Plugging then their expression in (2.35), it comes

$$ \left\{ \begin{aligned} U^T &= \frac{\partial S}{\partial U^T} = U^T + \frac{\partial U}{\partial U} \tilde{Y} - \left(\frac{\partial \vartheta}{\partial U}\right)^T \\ \frac{\partial \mathcal{G}}{\partial U^T} &= U \left(\frac{\partial \mathcal{F}}{\partial U}\right)^T = U \left(\frac{\partial \mathcal{F}}{\partial U}\right)^T + \mathcal{F} - \frac{\partial \Upsilon}{\partial V}. \end{aligned} \right. \quad (2.37) $$

Knowing that $\frac{\partial U}{\partial U} \tilde{Y} = \tilde{Y}$, it thus comes by identification that

$$ \tilde{Y} = \left(\frac{\partial \vartheta}{\partial U}\right)^T \quad \text{and} \quad \mathcal{F} = \left(\frac{\partial \Upsilon}{\partial U}\right)^T. $$

The functions $\vartheta$ and $\Upsilon$ are therefore potentials of the solution and the system flux, respectively. Furthermore, considering the Jacobian of the flux with respect to the vector variable, we observe

$$ \left(\frac{\partial \mathcal{F}}{\partial U}\right)^T = \left(\frac{\partial^2 \Upsilon}{\partial U^2}\right)^T, \quad (2.38) $$

showing the symmetrical structure. There is then only left to show that this symmetry transfers to the system.
To this aim, we start by denoting $A_0$ the Hessian of the entropy function $S$ with respect to the state vector $\tilde{Y}$, and observe that since $U^T = \nabla_{\tilde{Y}} S$, $A_0 = (\frac{\partial U}{\partial \tilde{Y}})$. Using again the chain rule and transposing (2.38), we there obtain

$$ \frac{\partial \mathcal{F}}{\partial U} = \frac{\partial \mathcal{F}}{\partial \tilde{Y}} \frac{\partial \tilde{Y}}{\partial U} = \frac{\partial \mathcal{F}}{\partial \tilde{Y}} A_0^{-1} = \frac{\partial^2 \Upsilon}{\partial U^2}. $$

Defining then a set $A = (A_1, \dots, A_d)$ of $d$ matrices $A_i \in \mathcal{M}^{q \times q}(\mathbb{R})$, $i \in [1, d]$, we can consider any $\alpha = (\alpha_1, \dots, \alpha_d)^T$ vector of $\mathbb{R}^d$ and construct the matrix

$$ \partial \mathcal{F}_{\alpha} := \sum_{i=1}^{d} \alpha_i \frac{\partial \mathcal{F}}{\partial \tilde{Y}} = \sum_{i=1}^{d} \alpha_i A_0 \left( \frac{\partial^2 \Upsilon_i}{\partial U^2} \right)^T = A_0 \sum_{i=1}^{d} \alpha_i \left( \frac{\partial^2 \Upsilon_i}{\partial U^2} \right)^T. \quad (2.39) $$

This last matrix being symmetric for any $\alpha \in \mathbb{R}^d$, the system (2.1a) itself is, yielding immediately its hyperbolicity.

**System hyperbolicity** A straightforward consequence of the symmetric nature of the system of Euler equations is its hyperbolicity. To observe it, one can simply use the relation (2.39) combined with the fact that $A_0 > 0$ to derive

$$ A_0^{1/2} (\partial \mathcal{F}_\alpha) A_0^{-1/2} = A_0^{1/2} \left( \sum_{i=1}^{d} \alpha_i \left( \frac{\partial^2 \Upsilon_i}{\partial U^2} \right)^T \right) A_0^{1/2} := \partial \Upsilon_\alpha, \quad (2.40) $$

from where it follows that the matrix $\partial \mathcal{F}_\alpha$ is diagonalisable in $\mathbb{R}$ as a consequence of the following lemma.

### Lemma 2.1.1 Orthogonality of the eigenvectors

For any $\alpha$ given vector of $\mathbb{R}^d$, the set of eigenvectors $\{r_1, \dots, r_d\}$ associated to the matrix $\partial \mathcal{F}_\alpha$ is orthogonal for the operator

$$ \langle \cdot, \cdot \rangle_{A_0} : \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R} \\ (\tilde{\eta}, \eta) \mapsto \tilde{\eta}^T A_0 \eta $$

**Proof.** The result is immediate by combining the relations (2.39) and (2.40). Indeed, since the matrix $\partial \Upsilon_\alpha$ is symmetric, it admits a set of orthonormal eigenvectors for the standard Euclidean scalar product, that we will denote $\{\eta_1, \dots, \eta_n\}$. Denoting then $\{\lambda_j\}_{j \in [1, d]}$ the set of eigenvalues of the matrix $\partial \Upsilon_\alpha$, one can write the relation $\partial \Upsilon_\alpha \eta_j = \lambda_j \eta_j$ up to a possible permutation of the eigenvalues indices. It then comes by (2.40) that

$$ A_0^{1/2} \left( \sum_{i=1}^{d} \alpha_i \left( \frac{\partial^2 \Upsilon_i}{\partial U^2} \right)^T \right) A_0^{1/2} \eta_j = \lambda_j \eta_j = A_0^{-1/2} (\partial \mathcal{F}_\alpha) A_0^{1/2} \eta_j, $$

yielding finally $(\partial \mathcal{F}_\alpha) A_0^{1/2} \eta_j = \lambda_j A_0^{1/2} \eta_j$, and concluding the proof.

As a consequence, the set of eigenvectors $\{r_i\}_{i \in [1, d]}$ of $(\partial \mathcal{F}_\alpha)$ is of the form $\{A_0^{1/2} \eta_i\}_{i \in [1, d]}$ where $\{\eta_j\}_{j \in [1, d]}$ is an orthogonal set of eigenvectors of $\partial \Upsilon_\alpha$. Thus, $\eta_i^T \eta_j = (A_0^{-1/2} r_j)^T (A_0^{-1/2} r_i) = r_j^T A_0^{-1} r_i = \langle r_j, r_i \rangle_{A_0^{-1}} = \delta_i^j$, ensuring that the matrix $\partial \mathcal{F}_\alpha$ is diagonalisable for any $\alpha \in \mathbb{R}^d$ and defining the hyperbolic nature of the system of Euler equations.
While the problem based on the system of Euler equations investigated in Section 2.1.2 remains an example among the class of problems that can be formulated as (2.1), its hyperbolic nature is a crucial feature that is common to all the problems considered in this thesis. We therefore now explore on a more general level the implications that the hyperbolic nature of a general system has on its solutions.

### 2.1.3 Hyperbolicity and discontinuity propagation

As seen in the above Section 2.1.2, the hyperbolicity of the system (2.1a) usually comes from its symmetry, being itself a consequence of the physical principles considered for modelling the dynamics of the investigated phenomenon. Stating in essence that the Jacobian of the system is diagonalisable, a possible definition reads as follows.

**Definition 2.1.2 Hyperbolicity of a conservative system**

Let $A_i(Y) \in \mathbb{R}^{q \times q}$ be the Jacobian matrix associated to the flux function $F_i: \mathbb{R}^q \to \mathbb{R}^q, i \in [1, d]$ evaluated on a given state $Y \in \mathbb{R}^q$. Then, if for any $\alpha = (\alpha_1, \dots, \alpha_d) \in \mathbb{R}^d$ and any state $Y \in \mathbb{R}^q$ the matrix $\sum_{i=1}^d \alpha_i A_i(Y)$ has real eigenvalues whose associated eigenvectors span $\mathbb{R}^q$, the system (2.1a) is said to be hyperbolic.

While the hyperbolicity of the system (2.1a) is a necessary condition for the problem (2.1a) to be well posed in the sense of Hadamard, it is in general far to be a sufficient condition [22, 23, 24, 25]. Furthermore, any physically relevant problem formulated in the form (2.1a) will necessarily have solutions developing a finite set of discontinuities, even for smooth initial data [26, 27].

**Example**

Let us observe the development of a discontinuity in the one-dimensional case by investigating the behaviour of the solution $Y: \mathbb{R}_+ \times \mathbb{R} \to \mathbb{R}$ to the problem

$$ \left\{ \begin{array}{l} \frac{\partial Y}{\partial t} + \frac{\partial F}{\partial x}(Y) = 0 \\ Y(0, \cdot) \in C^1(\Omega) \text{ given,} \end{array} \right. \qquad (2.41) $$

where $F: \mathbb{R} \to \mathbb{R}$ represents a scalar flux. The initial condition being smooth, we can denote by $A = \frac{dF}{dx}$ the flux's Jacobian and start by deriving the non-conservative form of the equation (2.41), that is, $\frac{\partial Y}{\partial t} + A(Y) \frac{\partial Y}{\partial x} = 0$. One can then look for a solution by considering the associated characteristic curves $\{t \mapsto x_\xi(t)\}_{\xi \in \mathbb{R}}$, solutions to the ordinary differential equations

$$ \left\{ \begin{array}{l} \frac{dx_{\xi}}{dt} = A(Y(x_{\xi}(t), t)) = 0 \\ x_{\xi}(0) = \xi \in \mathbb{R}. \end{array} \right. $$
Each solution being constant along the characteristic, we can write

$$ \frac{dY(x_{\xi}(t), t)}{dt} = \frac{\partial Y}{\partial t}(x_{\xi}(t), t) + A(Y(x_{\xi}(t), t)) \frac{\partial Y}{\partial x}(x_{\xi}(t), t) = 0, $$

providing the time-dependent solution for each initial point in $\mathbb{R}$ by tracing back the relation $x_{\xi}(t) = \xi + t A(Y(0, \xi))$. Assuming further that $\xi \mapsto A(Y(0, \xi))$ is decreasing on some interval $I \subset \mathbb{R}$, there exists two characteristic foots $r, s \in \mathbb{R}$ for which $A(Y(0, \xi_r(t))) < A(Y(0, \xi_s(t)))$ from a time $t < \infty$ on. The characteristics would then intersect, breaking the hitherto smoothness and defining two different states for a same point in the time-space whenever $Y(0, \xi_r) \neq Y(0, \xi_s)$.

<figure><img src="image_9.png" /><figcaption>**Figure 2.6.** Emergence of discontinuities (red circles).</figcaption></figure>

As each discontinuity will be either preserved or smoothened in time from their appearance on, the problem (2.1) would not admit classical solutions at any time, being sought under strong regularity conditions. We therefore need a more general concept of solutions, that is, weak solutions.

### 2.1.4 Weak solutions

The concept of weak solutions relies on reducing the regularity requirement on the sought solutions while still providing a meaning for the partial differential equations that rule them. Especially, in the context of hyperbolic systems of conservation laws, it is known a-priori that the locations in $\Omega$ where the solution $Y$ is discontinuous form a subspace of null measure. One can thus safely look for solutions living in $Y \in (L_{\text{loc}}^{\infty}(\mathbb{R}_+ \times \Omega))^q, q \in \mathbb{N}$ that solve (2.1a) in the sense of distributions.

**Definition 2.1.3** *Weak solutions*

Let us consider the problem (2.1a) for a given initial condition $Y(0, \cdot) \in L_{\text{loc}}^{\infty}(\Omega, \mathbb{R}^q)$. There, any function $Y \in L_{\text{loc}}^{\infty}(\mathbb{R}_+ \times \Omega, \mathbb{R}^q)$ satisfying

$$ \int_{\mathbb{R}_+} \int_{\Omega} Y \frac{\partial \varphi}{\partial t} + \mathcal{F}(Y) \cdot \nabla \varphi \, dx \, dt + \int_{\Omega} \varphi(0, x) Y(0, x) \, dx = 0 \quad (2.42) $$

for any $\varphi \in C_0^{\infty}(\mathbb{R}_+ \times \Omega, \mathbb{R})$ provides a weak solution to the problem (2.1a).

The functions $\varphi$ are called test functions, whose compact support guarantees the convergence of the integral and makes the interpretation of the system (2.1a) in the sense of distributions possible. Please also note that even though the test functions are here simply presented as scalar functions, alternative vectorial formulations can equally be derived (see e.g. Section 3.1.3).
With such a definition, the solution's derivative are not defined at specific points in $\Omega$ that either vanish later in time or are transported as such along given characteristics. Hence, more than considering the system (2.1a) in the sense of distributions, all the relations (2.19), (2.20) and (2.41) are also holding only almost everywhere, and the held discussions have hence to be understood in the distribution sense as well. Lastly, requiring regularity assumptions, the intuitive developments of Sections 2.1.2.1 and 2.1.2.3 only hold on smooth parts of the solution, and hold globally in a weaker sense, see e.g. [28] for extensive details.

Even though the existence of weak solutions has been proven [29, 27], uniqueness can in general not be guaranteed, and the exact shape of any solution is generally unknown. However, the impact of the created discontinuity on the later solution's development can be inferred by analysing the eigenstructure of the flux Jacobian upon the emergence of discontinuities.

### 2.1.5 Solution profiles

Now that we saw that there exists pathological behaviours and that a new concept of solutions could have been defined, let us describe the main solution profiles that can be encountered. We start in the one-dimensional case where extensive results are available, before opening on the multi-dimensional case and its subtleties. In any case, note that as a consequence of the regularity weakening introduced in Section 2.1.4, the upcoming equations (2.43)–(2.45) and (2.52) have to be understood in the distribution sense.

#### 2.1.5.1 Solution profile in one spatial dimension

Let us first investigate the possible types of developments for a discontinuity present in the solution to a one spatial dimensional problem. To this aim, we consider again the one-dimensional problem (2.41) in its non-conservative formulation and investigate the behaviour of the solution from the time at which the discontinuity is created. Assuming without loss of generality that the discontinuity appears at $t = 0$ and is located at $x = 0$, we can resolve the behaviour of the solution around the discontinuity by studying the solution to the well-known Riemann problem

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial Y}{\partial t} + \frac{\partial F}{\partial x}(Y) = 0 \\[1.5em] Y(0, x) = \begin{cases} Y_L & \text{if } x < 0 \\ Y_R & \text{otherwise} \end{cases} \end{array} \right. \qquad (2.43) $$

on $\mathbb{R}_+ \times \mathbb{R}$, where $Y_L, Y_R \in \mathbb{R}^q$ are the values of the discontinuity coming from the left and right sides, respectively. More precisely, the impact of the initial
discontinuity on the solution to the problem (2.1a) can be inferred by defining an immediate vicinity $I \subset \mathbb{R}$, $0 \in I$ and analysing there the behaviour of the solution to (2.43), the vicinity interval possibly widening in time according to the discontinuity's development.

**Example**

In the previous example, the discontinuity would appear at $t_b$ and be located at $\xi_r(t_b) = \xi_s(t_b)$. The left and right side values of the discontinuity would thus respectively read $Y_L = Y(0, \xi_r)$ and $Y_R = Y(0, \xi_s)$, obtained by tracing back their respective characteristics.

The widening of the interval of interest depends on the development of the discontinuity, whose nature and properties are intricate with the eigenstructure of the Jacobian matrix $A = \frac{\partial F}{\partial x} \in \mathbb{R}^{q \times q}$, composed from $\lambda_k$, $k \in [1, q]$ non-zero eigenvalues and $\{r_k\}_{k \in [1, q]}$ linearly independent eigenvectors. In particular, it is well known (see e.g. [30, 31, 32]) that the initial discontinuity splits into $q-1$ simply connected intermediate states that are separated from each others by so-called waves (see Figure 2.7). Being possibly of different natures, those waves are travelling at the speed of $\lambda_k(Y)$. They define further characteristic fields, the nature of which being qualified either as linearly degenerated or genuinely non linear, depending on whether the term $\nabla\lambda_k(\tilde{Y}(\xi)) \cdot r_k(\tilde{Y}(\xi))$ vanishes or not, respectively.

The field $Y$ being discontinuous at $(x, t) = (0, 0)$, the $k^{th}$ wave is here asked to travel simultaneously at possibly different speeds $\lambda_k(Y)$, yielding either compressive or expansive scenarios, as depicted in Figure 2.7. Furthermore, even though the waves propagation pattern is known, neither the intermediate state values nor the specific wave's profile are. Thus, solving (2.43) reduces to finding simultaneously a set of intermediate states $Y_k$, $k \in [1, q-1]$, and a function $Y: \mathbb{R} \to \mathbb{R}^q$ that verifies (2.43) and whose expression specifies the nature of each wave connecting the left and right intermediate states.

The problem being therefore global, let us first investigate the possible waves natures before defining admissible global solutions.

<figure><img src="image_10.png" /><figcaption>(a) Discontinuity decomposition</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Transition layers at a time $t > 0$</figcaption></figure>

**Figure 2.7.** Discontinuous solution splitting into several constant states (a) linked by (grey) transition layers (b).
**General nature of some solutions** In order to characterise the possible types of waves that would provide a local solution to (2.43) and describe the change from one intermediate state to another, let us consider a generic $k^{\text{th}}$ wave, $k \in [1, q]$, and consider the problem

$$ 
\left\{
\begin{array}{l}
\displaystyle \frac{\partial Y}{\partial t} + \frac{\partial \mathcal{F}}{\partial x}(Y) = 0 \\[1.5em]
Y(0, x) = \begin{cases}
Y_{k-1} & \text{if } x < t \lambda_k(Y_{k-1}) \\
Y_k & \text{if } x > t \lambda_k(Y_k)
\end{cases}
\end{array}
\right.
\qquad (2.44)
$$

where $Y_0 := Y_L$ and $Y_q := Y_R$ represent the left and right original states. We can here first observe that the system (2.44) is quasi-linear, and therefore look for self-similar solutions that can be either non-constant but regular or discontinuous piecewise constant.

**Continuity recovery and rarefaction waves** Let us first investigate if there exists non-constant self-similar regular solutions by plugging their general expression directly in (2.44). Namely, we define a function $\tilde{Y}: \mathbb{R} \to \mathbb{R}^q$ that verifies $Y(x, t) = \tilde{Y}(\frac{x}{t})$ at any point $x \in \mathbb{R}$ and for any time $t \in \mathbb{R}_+$ and set $\xi := \frac{x}{t}$. The function $\tilde{Y}$ being thus self-similar solution to (2.44) can then be given by solving

$$ 
\left\{
\begin{array}{l}
\displaystyle \frac{d\tilde{Y}}{d\xi}(\xi) = r_k(\tilde{Y}(\xi)) \\[1.5em]
\lambda_k(\tilde{Y}(\xi)) = \xi,
\end{array}
\right.
\qquad (2.45)
$$

for some $k^{\text{th}}$ eigenvalue $\lambda_k$ and its associated eigenvector $r_k$ of $A = \frac{\partial \mathcal{F}}{\partial x} \in \mathbb{R}^{q \times q}$, both depending on the state $Y$. Deriving the last equation with respect to $\xi$, one can first observe that the relation $\nabla \lambda_k(\tilde{Y}(\xi)) \cdot r_k(\tilde{Y}(\xi)) = 1$ is a necessary condition for a solution to exist, requiring in particular the $k^{\text{th}}$ characteristic field to be genuinely non-linear. When this condition is fulfilled, the regularity of $\mathcal{F}$ ensures the existence and uniqueness of a solution $\tilde{Y}: \mathbb{R} \to \mathbb{R}^q$ to (2.44) for any values of the intermediate states $Y_{k-1}$ and $Y_k$. The $k^{\text{th}}$ wave then takes the form

$$ Y(x, t) = \begin{cases} Y_{k-1} & \text{if } \xi \le \lambda_k(Y_{k-1}), \\ \tilde{Y}(\frac{x}{t}) & \text{if } \lambda_k(Y_{k-1}) \le \xi \le \lambda_k(Y_k) \\ Y_k & \text{if } \xi \ge \lambda_k(Y_k), \end{cases} $$

recovering the continuity of the solution to (2.43) for any time $t > 0$ by connecting smoothly the two intermediate states. Such solutions are referred to as *rarefaction waves*, an illustration of which has been provided in Figure 2.7b while exemplifying the solution behaviour in the first and third transition layer.
**Discontinuity propagation** Let us now ask ourselves if the discontinuity can be preserved through time, and look for discontinuous piecewise constants self-similar solutions that connect two intermediate states $Y_{k-1}$ and $Y_k$. There, defining functions of the type

$$ Y(t, x) = \begin{cases} Y_{k-1} & \text{if } x < c_k t, \\ Y_k & \text{if } x > c_k t, \end{cases} \qquad (2.46) $$

for a coefficient $c_k \in \mathbb{R}$ representing the wave's propagation speed to be determined, we notice that any such function is automatically a solution to (2.44). Hence, there exists an infinite amount of solutions propagating the discontinuity without altering it.

**A multitude of solutions** As we investigated the admissible wave types, we observed that there is no guarantee of their uniqueness even for given intermediate states. Hence, the uniqueness of the solution cannot be guaranteed for the problem (2.43), and selection criteria are therefore required to select the physically relevant ones among their possibly infinite admissible set. Let us first derive specific criterion for selecting locally the relevant discontinuous wave, before deriving a global criterion allowing to select the relevant intermediate states.

### 2.1.5.2 Relevant solutions in one spatial dimension

Let us first focus on selecting physical waves of the type (2.46). To this aim, we start by remembering that as being part of the global solution to (2.43) and more largely to (2.1a), a wave is physically relevant only if its expression fulfils the conservation principle (2.20) on its definition domain, that is, across the discontinuity. It therefore reduces to ask for the propagation speed $c_k \in \mathbb{R}$ to guarantee the conservation principles (2.19) associated to the system (2.1a) at almost anytime and for any compact subdomain containing the discontinuity. Writing

$$ \int_{\omega} \int_{0}^{T} \frac{\partial Y}{\partial t} + \frac{\partial \mathcal{F}}{\partial x}(Y) \, dx \, dt = 0 $$

for a given finite time $T \gg 1$ and any compact support $\omega \subset \mathbb{R}$ containing the discontinuity, it reduces after straightforward computations to

$$ -c_k (Y_k - Y_{k-1}) + (\mathcal{F}(Y_k) - \mathcal{F}(Y_{k-1})) = 0, \qquad (2.47) $$

the expression (2.47) being known as the *Rankine-Hugoniot relations*. Seen in a different light, those relations can also be referred to as the *jump conditions*, relating the two states of each scalar variable to a same discontinuity's propagation speed.
**Contact discontinuities**
Coupled with those last relations, the Riemann problem (2.43) admits a unique discontinuous solution whenever the associated $k^{th}$ characteristic field is linearly degenerated. Indeed, applying (2.47) to a field where $\nabla\lambda_k(Y) \cdot r_k(Y) = 0$ necessarily reduces to $c_k = \lambda_k(Y_{k-1}) \equiv \lambda_k(Y_k)$, defining a wave that propagates through the solution at the same velocity as the two characteristics emanating from the left and right intermediate states. The discontinuity may therefore not appear in all the conservative variables of the problem (2.1a), especially when they directly relate to the internal velocity of the investigated quantity. Hence providing an additional relation between the intermediate states, this type of wave is referred to as a contact discontinuity and is the only one encountered when the flux $\mathcal{F}$ is linear.

**Shock discontinuities**
In the case of genuinely non-linear characteristic however, the Rankine-Hugoniot relations do not suffice anymore to ensure uniqueness as the non linearity prevents a straightforward relationship between the waves and characteristic velocities. The discontinuity is then referred to as a shock, and one requires a further selection criterion.

<figure><img src="image_12.png" /><figcaption>**Figure 2.8.** Possible solution profiles around shock discontinuities (thick line). A physical solution sees its characteristics (thin lines) merging into the shock.</figcaption></figure>

We therefore ask the solution to be physically realistic, namely that existing quantities should rather merge into the shock than being spontaneously created along it (see Figure 2.8 for a characteristic’s perspective). Referred to as the *Lax shock conditions*, this principle intuitively writes

$$ \left\{ \begin{array}{l} \lambda_k(Y_k) < c < \lambda_{k+1}(Y_k) \\ \lambda_{k+1}(Y_{k-1}) < c < \lambda_k(Y_{k-1}) \end{array} \right. \qquad (2.48) $$

for convex fluxes, and has been generalised by Oleinik [33] to

$$ \left\{ \begin{array}{l} \mathcal{F}(Y_k) - \mathcal{F}(Y_{k+1}) < c_k (Y_k - Y_{k+1}) \\ c_k (Y_k - Y_{k-1}) < \mathcal{F}(Y_k) - \mathcal{F}(Y_{k-1}). \end{array} \right. \qquad (2.49) $$

This selection criterion so defined, the problem (2.44) coupled with the two Rankine-Hugoniot (2.47) and Oleinik (2.49) relations admits a unique discontinuous solution [22, 34], whose discontinuity is fully preserved through time and travels at a shock speed $c_k$.
**Uniqueness of the global solution** So far, analysing the structure of the characteristic fields allowed us to derive the possible natures of each wave, and asking for the physical relevance of the transition profile between two given intermediate states made the selection a single wave possible in case of contacts and shocks. As a consequence, a unique solution to (2.1a) can already be defined in the linear case, the global solution being a piecewise constant function where the states are separated by contact waves whose velocities are the solutions to a system of Rankine-Hugoniot relations.

However, in the non-linear case those two conditions alone do not necessarily guarantee the physical relevance of the global solutions nor a uniqueness property. Indeed, the values of the intermediate states could not be uniquely determined but only interrelated at contacts, and there is still an indetermination of the wave type itself for genuinely non-linear characteristic fields (see Figure 2.9 for an illustration). If the existence of a physical concept ensuring uniqueness is still an active debate, it is nonetheless possible to ask for a last condition that ensures the globally solutions to (2.1a) to be physically admissible. Typical criterion relates to the notions of viscosity [35, 36], entropy [33, 30], and entropy rate [37, 38]. We choose here to focus on preventing unrealistic behaviour by including considerations about the energy dissipation, and thus consider in the following entropy solutions, whose definition will be specified in Section 2.1.6.

<figure><img src="image_13.png" /><figcaption>**Figure 2.9.** An example of possible solutions that may fulfil an entropy criterion.</figcaption></figure>

**Note** In the scalar one-dimensional case, the Rankine-Hugoniot and Oleinik conditions are enough to single out the physically relevant solution. Indeed, as there is no intermediate state the global admissibility of the solution is immediate given the admissibility of the only wave itself. And since $q = 1$, the Oleinik relations reduces to simply choosing a shock wave when $Y_L > Y_R$ (the state is compressed) and a rarefaction otherwise (the state is continuously relaxed). In case of a shock, the Rankine-Hugoniot then naturally determines its speed given the left and right states, while in relaxation case, the $C^1$-regularity of $F$ grants uniqueness to the problem (2.44).
**Example**

Let us consider the Riemann Problem for the one-dimensional space Euler equations that has been previously detailed in Section 2.1.2.2, where we consider $Y_L > Y_R$. There, it is well-known that the three eigenvalues of the flux Jacobian read $\lambda_1 = v - c$, $\lambda_2 = v - c$, and $\lambda_3 = v + c$, where $c = \sqrt{\frac{\gamma_p}{\rho}}$. Moreover, there exists a unique solution that is further entropic, for which the characteristic fields associated to $\lambda_1, \lambda_2, \lambda_3$ are respectively a relaxation wave, a contact wave, and a shock wave.

<figure><img src="image_14.png" /><figcaption>**Figure 2.10.** Density profile of the entropy weak solution to the one-dimensional Euler equations. Dashed line: initial condition. Solid line: after the discontinuity split.</figcaption></figure>

### 2.1.5.3 Solutions to higher dimensional problems

The above qualitative observations transfer to the higher spatial dimensional case, in which the lack of uniqueness for the solutions to the problem (2.1) is even more blatant. However, if the development of discontinuities still occurs and their classification as shocks, contacts or relaxations is still relevant, it is in general not possible to investigate *a-priori* the nature of the characteristic fields. Indeed, the divergence operator combines the effects of the direction- wise fluxes to the solution, making the interplay hard to analyse from the eigenvalues of the individual matrices of the non-conservative form. If some approaches tried to determine a suitable analytical framework by considering directional behaviours [39, 40], there is still no consensual theory. However, some admissibility criterion to retrieve physically relevant waves could have been generalised. In particular, when considering a $k^{th}$ characteristic field for which the associated wave is discontinuous, the Rankine-Hugoniot relations become

$$c_k[Y] = \sum_{i=1}^{d} n_i[\mathcal{F}_i(Y)], \qquad (2.50)$$

where $n = (n_1, \dots, n_d) \in \mathbb{R}^d$ represents the unit normal to the discontinuity's characteristic and $[\cdot]$ denotes the jump across it. Complementarily, the extension of the Lax shock condition to higher spatial dimensions directly leads to the more general notion of entropy solutions.
## 2.1.6 Entropy solutions

As seen in the above sections, there is generally no uniqueness of the solutions to the problem (2.1a), and additional criterion are required to extract the subset of physically relevant ones. If the conservation principle is obtained by asking the Rankine-Hugoniot relations to hold across any discontinuity, those relations are not enough for preventing the selected solutions to show unphysical behaviours (see the last paragraph of Section 2.1.3).

**A general notion of entropy** A popular way to exclude physically unreal-istic behaviours from the set of conservative solutions is to pay attention to the amount of internal energy that is available to the system for performing motions or transformations, quantity defined as the system's entropy. There,

considering physically relevant behaviours reduces to asking the system to consume no more energy during its evolution than its entropy state allows. The system's behaviour being dependent on the solution state, this last requirement finally reduces to ensuring that the internal energy of the solution itself is never decreasing, providing a selection criteria at the solution level.

However, the representation of the available internal energy for a system to evolve is not necessarily unique, and neither necessarily

existing. When it does, which is the case for most systems that are derived from physical principles, it relates to so-called entropy fluxes that describe the use of the energy during the system's transformation.

<figure><img src="image_15.png" /><figcaption>**Figure 2.11.** Schematic description of the system entropy.</figcaption></figure>

**Definition 2.1.4** *Entropy relations*

A convex function $S: \Omega \to \mathbb{R}$ is an entropy function associated to the problem (2.1a) if there exists $d$ functions $G_j: \Omega \to \Omega$ functions, called entropy fluxes verifying

$$ \nabla_Y G_i(Y) = (\nabla_Y S(Y))^T A_i(Y), \quad \forall i \in [1, d], $$

where the term $A_i$ denotes the Jacobian of the flux function $F_i, i \in [1, d]$.

**Note** As it is easier to deal with convex functions, the mathematical entropy is defined in Definition 2.1.4 as the opposite of the physical entropy, without changing its fundamental meaning. Thus, being convex, the Lax and Oleinik conditions derived in the previous Section 2.1.5.2 can be seen as resulting from a particular choice of entropy functions.

<figure><img src="image_16.png" /></figure>
**Entropy solutions and discontinuities** So defined, one can observe that there exists as many entropy conservation laws as admissible entropy pair $(S, G)$. Indeed, for any entropy pair $(S, G)$ given by Definition 2.1.4, multiplying the system (2.1a) by the entropy functions' state-gradient $\nabla_Y S(Y)$ yields

$$ 
\begin{aligned}
& \nabla_Y S(Y) \frac{\partial Y}{\partial t} + \sum_{i=1}^{d} (\nabla_Y S(Y))^T A_i(Y) \frac{\partial Y}{\partial x_i} = 0 \\
\Leftrightarrow & \frac{\partial S(Y)}{\partial t} + \sum_{i=1}^{d} \nabla_Y G(Y) \frac{\partial Y}{\partial x_i} = 0 \\
\Leftrightarrow & \frac{\partial S(Y)}{\partial t} + \sum_{i=1}^{d} \frac{\partial G_i}{\partial x_i}(Y) = 0,
\end{aligned} 
\quad (2.51) $$

holding wherever the solution is smooth. Any smooth solution to (2.1a) is therefore naturally entropic and physically relevant. On non-smooth parts of the solution however, that is in case of contacts or shocks, the relation (2.51) cannot be derived as such, the chain rule not applying. Then, one can show that across a contact the relation also holds in a distribution sense, provided that the Rankine-Hugoniot relations (2.50) are satisfied. It is however not the case across shocks, where the second law of thermodynamics tells us that even if the internal energy of the system available for it to evolve the solution should still not decrease, it may be dissipative. The entropy constraint associated to a given entropy pair $(S, G)$ is therefore relaxed to asking

$$ \frac{\partial S(Y)}{\partial t} + \sum_{i=1}^{d} \frac{\partial G_i}{\partial x_i}(Y) \le 0. \quad (2.52) $$

to hold in a distribution sense. The selection of a physically relevant solution is then solely driven by the relation (2.52), being only active across shocks.

**Weak entropy solutions** As seen in the above paragraph, physically admissible solutions should match with the time-space dynamic as dictated by the problem (2.1) while fulfilling the entropy relations (2.52) in the sense of distributions. Hence, all considerations taken into account, any physically admissible solution should enters in the scope of weak entropy solutions.

**Definition 2.1.5** *Weak entropy solutions*

A function $Y \in C^\infty(\mathbb{R}_+ \times \Omega)$ is said to be a weak entropy solution if it provides a solution to the problem (2.42), and that for any entropy pair $(S, G)$ given in Definition 2.1.4 it holds

$$ \int_{\mathbb{R}_+} \int_{\Omega} S \frac{\partial \varphi}{\partial t} + G(Y) \cdot \nabla \varphi \, dx \, dt + \int_{\Omega} \varphi(0, x) S(0, x) \, dx \le 0. \quad (2.53) $$
In general, not any weak solution given by Definition 2.1.4 fulfils the relation (2.53) for all entropies of the system (2.1a). However, it is well-known that any weak solution that can be written as a vanishing viscosity limit associated to the strong form of the problem (2.1a) is entropic.

**Theorem 2.1.1** *Vanishing viscosity limit*

Let us define the sequence of smooth functions $\{Y^\epsilon \in C^\infty(\mathbb{R}_+ \times \Omega)\}_{\epsilon>0}$ that satisfy the parabolic problem

$$ \frac{\partial Y^\epsilon}{\partial t} + \nabla \cdot \mathcal{F}(Y^\epsilon) = \epsilon \Delta Y. \qquad (2.54) $$

If

* the sequence $\{Y^\epsilon\}_{\epsilon>0}$ admits a limit as $\epsilon \to 0$,

* there exists a constant $C \in \mathbb{R}$ such that $\sup_{\epsilon>0} \|Y^\epsilon\|_{L^\infty(\mathbb{R}_+ \times \Omega, \mathbb{R}^q)} < C$,

then the function $Y = \lim_{\epsilon \to 0} Y^\epsilon$ belongs to $L_{\text{loc}}^\infty(\mathbb{R}_+ \times \Omega)$, is a weak solution to the problem (2.42), and satisfies the relations

$$ \frac{\partial S(Y)}{\partial t} + \sum_{i=1}^{d} \frac{\partial \mathcal{G}}{\partial x_i}(Y) \le 0 \qquad (2.55) $$

in the sense of distributions for any entropy pair $(S, \mathcal{G})$ associated to (2.1a).

**Proof.** The inequality comes by the convexity of the entropy function, obtained when multiplying (2.54) by the gradient of $S$. The existence to viscosity solution is only due to the parabolic nature of (2.54) and can be found in e.g. [24]. Classical, the convergence proof can be found in e.g. [21].

However, even though this relationship provides a powerful selection criterion, it does not apply to any problem of the form (2.1a). Indeed, in the nonlinear case, due to possible oscillations the norm $\sup_{\epsilon>0} \|Y^\epsilon\|_{L^\infty(\mathbb{R}_+ \times \Omega, \mathbb{R}^q)}$ may not be bounded even if all the terms $\|Y^\epsilon\|_{L^\infty(\mathbb{R}_+ \times \Omega, \mathbb{R}^q)}$ are. Hence, there is no guarantee here that the viscous solutions converge to an entropy solution.

**A note about uniqueness of entropy solutions** While asking a weak solution to be entropic in the sense of Definition 2.1.5 ensure its physical relevance with respect to the associated entropy pairs $(S, \mathcal{G})$, an entropy solution is still far from being unique [41]. Furthermore, if the conservative (2.50) and entropy (2.51) relations allow to select physically relevant solutions, they do not enforce smooth solutions nor smooth the discontinuous ones out. Indeed, the discontinuities are intrinsic features of the solutions and not artefacts.
## 2.2 A brief algebraic toolbox

Having defined the essential properties that are common to the investigated problems, let us now introduce some algebraic notions that are of importance in their approximation. We start in Section 2.2.1 by detailing the properties of multivariate polynomial spaces, before presenting the concept of degrees of freedom in Section 2.2.2. We then describe their relationship with the notion of basis functions in Section 2.2.3, which we address some delicate geometrical concern in Section 2.2.4.

### 2.2.1 Multivariate polynomial spaces

Finding a discrete approximation of the solution to the problem (2.1a) involves in particular the definition of finite dimensional spaces where the approximate solution can be represented. Hence, the development of suitable approximation spaces is usually based on polynomial ones, whose dimension are known and whose definitions are easy to handle. In particular, the representation of a solution living on a multi-dimensional domain $K \subset \mathbb{R}^d, d \in \mathbb{N}$ involves an interplay between the components of the position variables $x = (x_1, \dots, x_d)$, making necessary the use of multivariate spaces.

**Finite dimensional multivariate polynomial spaces** By nature, any finite dimensional multivariate polynomial space is necessarily defined as a restriction of the more general space $\mathbb{P}(K)$ of all $d$-variate polynomials having the polytope $K$ for definition domain. The selection of a finite dimensional subspace of discretisation $V$ is therefore usually based on the polynomial degree, further providing a simple indicator on the discretisation fineness. The definition of the degree not being unique, there exists several classical definitions of finite dimensional polynomial spaces.

**Note** As any polynomial $p \in V$ can be written as a finite sum of monomial terms $p_l, l \in [1, m]$, $m \in \mathbb{N}$, whose degree is then always defined by $\deg(p) = \max(\{\deg(p_l)\}_{l \in [1, m]})$, we only discuss here the degree on monomials, the definition immediately transferring to the polynomial case.

**Polynomial degrees and finite dimensional subspaces** Let us derive three common classes of finite dimensional polynomial subspaces by analysing the natural definitions for the degree of the multivariate monomial

$$p: x \mapsto \prod_{i=1}^{d} x_i^{\alpha_i} \in \mathbb{P}(K).$$

There, a first natural choice of consists in considering the monomial as a single entity and interpreting the degree as a mixed growing rate. The degree
then reads $\deg_1(p) = \sum_{i=1}^d \alpha_i$, leading to the space

$$ \mathbb{P}_k(K) = \{q \in \mathbb{P}(K), \deg_1(q) \le k\}. $$

Another natural choice would rather be to consider the degree as representing the maximal unidirectional growth rate of the monomial, regardless the variable it applies to. Its definition would then read $\deg_2(p) = \max(\{\alpha_i\}_{i \in [1, d]})$, yielding to defining the finite dimensional subspace

$$ \mathbb{Q}_k(K) = \{q \in \mathbb{P}(K), \deg_2(q) \le k\}. $$

There, keeping the same definition of the degree, one can be more distinctive by prescribing the polynomial behaviour coordinate-wise, that is, imposing a maximal growth rate for each variable. The associated finite dimensional polynomial space thus directly writes

$$ \mathbb{P}_{k_1, \dots, k_d}(K) = \left\{ q \in \mathbb{P}(K), q = \sum_{l=1}^{m} \left( \beta_l \prod_{i=1}^{d} x_i^{\alpha_{il}} \right), \text{ s.t. } \begin{array}{l} \forall l \in [1, m], \\ \forall i \in [1, d], \\ \max_{l \in [1, m]} \{\alpha_{il}\} \le k_i \end{array} \right\}. $$

where $m \in \mathbb{N}$ represents the number of monomials composing the polynomial $q \in \mathbb{P}_k(K)$. While this definition might be delicate to manipulate in practice, requiring an *a-priori* knowledge of the suitable discretisation structure to prescribe the differentiated behaviours appropriately, the obtained space is very powerful for discretising an anisotropic problem, or when the quality of the discretisation projections is wished to be detailed variable-wise.

### Finite dimensional subspaces dimensions and embeddings

The main classical definitions of finite dimensional polynomial subspaces being set, let us now express their dimension with respect to the value of the selected polynomial degree, property that will be of crucial importance in the developments of Chapter 3. Using algebraic results [42] it directly comes

$$ \dim \mathbb{P}_k(K) = \frac{(k+d)!}{k!d!}, \quad \dim \mathbb{Q}_k(K) = (k+1)^d, \quad \text{and} \quad \dim \mathbb{P}_{k_1, \dots, k_d} = \prod_{l=1}^{d} (k_l + 1). $$

Furthermore, observing the structure and dimensions of those spaces, one can naturally obtain the embedding relations

$$ \mathbb{P}_{k_1, \dots, k_d}(K) \subset \mathbb{Q}_{\max\{k_i, i \in [1, d]\}}(K) \quad \text{and} \quad \mathbb{P}_{k, \dots, k} = \mathbb{Q}_k(K), $$

of importance when deriving projectors on polynomial discretisation spaces.

### Choices of polynomial discretisation spaces

As a last note, let us point out that while the discretization objective guides the choice of polynomial space structure, the geometry also plays a role. Indeed, the geometry determines the support of the polynomial functions, whose characterisation nature is also possibly related to it. The dimension of the space must therefore often be in accordance with the geometry of $K$, typically yielding to choosing e.g. spaces based on $\mathbb{P}_k(K)$ when $K$ is a simplicial form and on $\mathbb{Q}_k(K)$ when $K$ is a polytope with a higher number of faces.
## 2.2.2 Polynomial representation and degrees of freedom

Once a functional is discretised in a finite dimensional space, it is necessary to define a way it will be represented, and thus controlled by the numerical scheme. This representation can be either implicit through the definition of quantifiers, or explicit by using a set of basis functions. Either way, the discrete quantity will be controlled by the numerical scheme from a finite set of degrees of freedom values.

**Discrete representation and degrees of freedom** Let us consider some finite dimensional space $V$. Seeing it as a discretisation space, we can approximate any quantity of interest by projecting it onto $V$ and decompose it on a set of basis functions. There, one can see the coefficients of the decomposition as freedom in the definition of the function lying in the space $V$, freedom being expressed in terms of the behaviours described by each of the basis functions.

**Example**

Let us consider the finite dimensional space $P_2(\mathbb{R})$. There, any projected quantity is a polynomial of degree at maximum two and can be decomposed on any basis of $P_2(\mathbb{R})$. For example, the polynomial $q: x \to 2x^2+3x+2$ can be decomposed on the canonical basis $\{x^2, x, 1\}$ and represented by its corresponding values 2, 3 and 2. Each of the monomials $x^2$, $x$ and 1 then represents a type of behaviour that one can encounter in any function living in $V$, their pondered completely describing any discrete function in $V$. Therefore, not being unique the choice of the basis function emphasises a set of typical behaviours that are present in any function of $V$ and whose weighted combination is enough to entirely describe any of them.

Another way to represent the quantity of interest is to fix those freedom not by projecting on the basis functions, but by specifying values of some quantifiers that are characterising the functions living in the discretisation space. The definition of a quantifier is then called a *degree of freedom*.

**Definition 2.2.1 Degrees of freedom**

A degree of freedom $\sigma: V \to \mathbb{R}$ is a linear form supported on $K$ and acting on $V$ whole value characterises some common feature to any function of $V$.

Provided that they define linear forms and that their definition is supported on $K$, the shape of the degrees of freedom is up to one's taste, its value then settling the function's freedom along the obtained property. Thus, to entirely determine any function of the discretisation space, one needs as many degrees of freedom as the dimension of the space, each of them emphasising a type of behaviour that can be found in any function of $V$ and being independent from the others. This set of degrees of freedom thus lives in the space that is dual to $V$, and the collection of their values entirely determines any function of $V$.
### Example

(Continued) Knowing that the polynomial $q: x \mapsto 2x^2 + 3x + 2$ belongs to $\mathbb{P}_2(\mathbb{R})$, it can be entirely determined by either of the following relations.

$$ (a) \begin{cases} q(0) = 2 \\ q(1) = 7 \\ q(-1) = 1 \end{cases} \quad (b) \begin{cases} q(0) = 2 \\ q(1) = 7 \\ \int_0^1 x q \, dx = 3 \end{cases} \quad (c) \begin{cases} \int_0^1 q \, dx = \frac{25}{6} \\ \int_0^1 x q \, dx = \frac{5}{2} \\ \int_0^1 x^2 q \, dx = \frac{115}{60} \end{cases} $$

There, the relations (a) are characterising the polynomial through *point values*, while the configuration (c) involves only *moments*. One can also mix the characterisation types and retrieve a *mixed* configuration as shown in the configuration (b). The degrees of freedom are then the linear forms that to $q$ associate the left hand side of those relations.

To see this how this characterisation acts on the discretisation space $V$, we can derive the linear system corresponding to a generic function of $\mathbb{P}_2(\mathbb{R})$ that is quantified by any set of degrees of freedom (a), (b) or (c). There, given any basis $\mathcal{B} = \{\varphi_1, \varphi_2, \varphi_3\}$ of $\mathbb{P}_2(\mathbb{R})$, each of the relations (a), (b) or (c) provides the values of the three weights coefficients associated to the basis $\mathcal{B}$. For example, we can look for a function $q$ of $\mathbb{P}_2(\mathbb{R})$ under the form $q: x \mapsto ax^2 + bx + c$, using the canonical basis for $\mathcal{B}$ and real coefficients $a, b$ and $c$. Seeing the configuration (b) as a set of degrees of freedom, we can derive the following determination system.

$$ \left\{ \begin{array}{l} a \times 0^2 + b \times 0 + c = 2 \\ a \times 1^2 + b \times 1 + c = 7 \\ \int_0^1 x (ax^2 + bx + c) \, dx = 3 \end{array} \right. $$

By linearity of the integral, it reduces to

$$ \begin{pmatrix} 0 & 0 & 1 \\ 1 & 1 & 1 \\ \frac{1}{3} & \frac{1}{2} & 1 \end{pmatrix} \begin{pmatrix} a \\ b \\ c \end{pmatrix} = \begin{pmatrix} 2 \\ 7 \\ 3 \end{pmatrix}, $$

from which the coefficients $a = 2$, $b = 3$ and $c = 2$ follow. Note also that at this point there is no particular relation between the chosen basis $\mathcal{B}$ and the used linear forms.

Lastly, one can observe that the rows of the above matrix, expressing the shape of the degrees of freedom applied to each basis function, are linearly independent. Thus, as the quantifiers act on $V$ and as we have as many rows as its dimension, the matrix is fully ranked, and invertible. The set of relations (b) then characterises uniquely a function of $\mathbb{P}_2(\mathbb{R})$ for any right-hand-side values.
The property that among a set of degrees of freedom each of them characterises a unique type of behaviour present in any function living in V translates into the linear independence of the equations obtained when the degrees of freedom are applied to any function of V. Thus, when there are as many relations as the dimension of the space V, the given set of degrees of freedom is *unisolvent* in V.

**Theorem 2.2.1** *Unisolvence of a set of degrees of freedom*

When a set of degrees of freedom {$\sigma_i$}$_{i \in [1, \dim V]}$ completely determines any element living in V, i.e. when for any given basis $\mathcal{B} := \{\varphi_j\}_{j \in [1, \dim V]}$ the matrix $(\sigma_i(\varphi_j))_{i,j}$ is invertible, the set of degrees of freedom is said to be *unisolvent* for the space V.

The definition of the degrees of freedom emphasises some properties or behaviours of functions one may wish to analyse. In particular, for any set of degrees of freedom there exists a basis of V that is dual to its own, i.e. writing {$\sigma_i$}$_{i \in [1, V]}$ a generic set of degrees of freedom, there exists a couple ($\{\sigma_i\}_i$, $\{\varphi_j\}_j$), $(i, j) \in [1, \dim V]^2$ such that

$$ \forall (i, j) \in [1, V]^2, \quad \sigma_i(\varphi_j) = \delta_{ij}. \qquad (2.56) $$

where {$\varphi_j$}$_{i \in [1, V]}$ is the dual basis representing all the behaviours the functions of V are decomposed on, while the degrees of freedom {$\sigma_i$}$_{i \in [1, V]}$ give a direct access the corresponding quantified properties.

**Example**

The following degrees of freedom are dual to the canonical basis of $\mathbb{P}_2(\mathbb{R})$.

$$ \left\{ \begin{array}{l} \sigma_1: q \mapsto q(0) \\ \sigma_2: q \mapsto \frac{q(2) - q(-1)}{3} \\ \sigma_3: q \mapsto q(1) - q(0) - \frac{q(2) - q(-1)}{3} \end{array} \right. $$

Indeed, denoting the functions of the canonical basis by $\varphi_1 := x \mapsto 1$, $\varphi_2 := x \mapsto x$, and $\varphi_2 := x \mapsto x^2$, it comes

$$ \{\sigma_i\}_{i \in [1, V]} = \begin{pmatrix} 1 & 0 & 0 \\ \frac{1-1}{3} & \frac{2-(-1)}{3} & 0 \\ (1^2 - 1^2 - \frac{1^2-1^2}{3}) & (1 - 0 - \frac{1^2-(-1)^2}{3}) & (1 - 0 - \frac{1^2-(-1)^2}{3}) \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} $$

One can further notice that computing the value of the degrees of freedom for a specific function provides the coefficients of its projection onto the dual basis. Note however that counter-intuitively, sets of moment based degrees of freedom (as e.g. the configuration (c) of the previous example) are not necessarily dual to their projection basis.
Let us also point out that even though this dual basis exists, its knowledge is not necessary to express the discretised quantities. Indeed, as shown above the degrees of freedom values are enough to characterise any function living in $V$ on their own. Thus, the knowledge of the values taken by each of the linear form defining a unisolvent set of degrees of freedom in $V$ is enough to represent uniquely a function of $V$. Furthermore, the same way it is possible to define several basis for the space $V$, it is possible to design several sets of set of degrees of freedom, and therefore to have several representations of the same quantity. As those two sets will be dual to two different basis of $V$, we can link them through a transfer matrix (see Section 2.2.3 for details).

In practice, the two common types of degrees of freedom are so-called *point-wise values* and *moments*. The formers consist in evaluating the discretised quantities at a precise point in $K$, while the latter is averaging the quantity tested against some function characterising a target behaviour on a subset of $K$. In particular, when using moment based degrees of freedom one can test functions living in $V$ against some subspace contained in $V$. However, testing against a full subspace is computationally unrealistic. Therefore, practically one only computes the degrees of freedom from a basis of the subspace.

### 2.2.3 Basis tuning on degrees of freedom

Two different sets of basis functions that generate a same vector space do not necessarily share the same properties. However, it is possible to find an bijective mapping between each of the sets. As a consequence, one can design a specific basis from any generic basis and define the set of degrees of freedom that express some wished properties through a duality relationship.

While aiming for a customised characterisation of functions living in $V$, one should pay attention so that the desired properties are compatible with the considered space. Typically, they should not contradict each other nor be equivalent, and for each property there should exist functions that belong to the vector space fulfilling it. Furthermore, the set of relations defining the properties should be of duality type, and the set of the corresponding quantification functions should have the dimension of the vector space. We assume in this paragraph that all the above requirements hold, and detail the mapping.

<figure><img src="image_17.png" /><figcaption>**Figure 2.12.** Duality between a set of degrees of freedom and its associated basis.</figcaption></figure>
Considering the vector space $V$ endowed with a basis $\{\phi_j\}_{j \in [1, \dim V]}$ and a set of real valued linear applications $\{\sigma_i\}_{i \in [1, \dim V]}$ acting on $V$, let us derive a mapping that constructs the basis $\{\varphi_j\}_{j \in [1, \dim V]}$ as the one that is dual to $\{\sigma_i\}_i$. First of all, as $\{\phi_j\}_j$ and $\{\varphi_j\}_j$ span the same vector space, we can write

$$ \varphi_j = \sum_{m=1}^{\dim V} \alpha_{jm} \phi_m \qquad (2.57) $$

for some coefficients $\{\alpha_{jm}\}_{(j,m) \in [1, \dim V]^2}$ that have to be found. Using then the degrees of freedom and the definitions of their related basis $\varphi$, the duality relations (2.56) expand to

$$ \sigma_i \left( \sum_{m=1}^{\dim V} \alpha_{jm} \phi_m \right) = \delta_{ij} $$

for any $i, j \in [1, \dim V]$. Then, by the linearity of the degrees of freedom we obtain

$$ \sum_{m=1}^{\dim V} \alpha_{jm} \sigma_i(\phi_m) = \delta_{ij} $$

for any $i, j \in [1, \dim V]$, relations that when rewritten in a matrix form read

$$ \Lambda A^T = I_N \qquad (2.58) $$

for $(\Lambda)_{ij} = \sigma_i(\phi_j)$ and $(A)_{jm} = \alpha_{jm}$. Assuming further that the set $\{\sigma_i\}_i$ is unisolvent for $V$, the matrix $\Lambda$ is then invertible and the coefficients $\alpha_{im}$ are simply obtained by computing $A = (\Lambda^{-1})^T$. The basis functions $\{\varphi_j\}_j$ then simply read through the relations (2.57). Furthermore, the partial derivatives being linear operator one can retrieve analogously the relations

$$ \left\{ \begin{array}{l} D_x \phi_{ij} = \sigma_i(\partial_x \varphi_j) \\ (D_x \varphi)_{ij} = (\sigma_i(\partial_x \varphi_j))_{ij} \\ D_x \phi A^T = D_x \phi \Lambda^{-1} = D_x, \varphi \end{array} \right. $$

for any $i, j \in [1, \dim V]$ and where $\varphi = (\varphi_1, \dots, \varphi_{\dim V})^T$. A comprehensive bijective mapping that turns any basis of $V$ to the one corresponding to set of degrees of freedom unisolvent for $V$ thus reads as follows.

### Theorem 2.2.2 Tuning of basis functions

Let $V$ be vector space endowed with a basis $\{\phi_j\}_{j \in [1, \dim V]}$ and define a set of degrees of freedom $\{\sigma_i\}_{i \in [1, \dim V]}$ that is unisolvent for $V$. The basis $\{\varphi_j\}_j$ that is dual to $\{\sigma_i\}_i$ can then be constructed from the relations

$$ \varphi_j(x) = \sum_{m=1}^{N} \alpha_{jm} \phi_m(x) \quad \text{and} \quad \partial_x \varphi_j(x) = \sum_{m=1}^{N} \alpha_{jm} \partial_x \phi_m(x), \qquad (2.59) $$

where the coefficients are given by $(\alpha_{jm})_{jm} = ((\sigma_i(\phi_j))_{ij})^{-T}$.
To conclude this section, let us mention that the choice of the degrees of freedom {$σ_i$}$_i$ matches the choice of the properties one wants to instill to the basis functions. Indeed, applying them to the set of their dual basis functions {$φ_j$}$_j$ will impose by duality the relations (2.56), defining properties for the basis functions. Therefore, the relations (2.59) will automatically tune the basis functions towards all the desired ($\dim V$)$^2$ properties.

**Example**

* A moment based degree of freedom by the application $\sigma: q \mapsto \frac{1}{|K|} \int_K q \, dx$ quantifies $q$ by its mean value on the cell $K$, and gives $\dim V$ properties to the set of dual basis functions. Indeed, for some basis function $\varphi_i$, $i \in [1, \dim V]$, it then holds $\frac{1}{|K|} \int_K \varphi_i \, dx = 1$, and for any $j \in [1, \dim V]$, $j \neq i$, it reduces to $\frac{1}{|K|} \int_K \varphi_j \, dx = 0$.

* The same observation also holds for pointwise degrees of freedom. There, assuming that one would like to prescribe the midpoint value to be zero for one specific basis function and one for all the others, one can simply set $\sigma: q \mapsto 1 - q(\frac{1}{2})$ as a quantification function.

### 2.2.4 Basis functions adjustment to geometry distortion

When designing a mesh-based numerical scheme, it is often useful to build a local discretisation setting $(K, V, \{\sigma\}_{i \in [1, \dim V]})$ on a reference geometry $\tilde{K}$ and map it onto the mesh cells $K$. However, by simply choosing any bijective geometrical mapping that would automatically describe the distortion between the two geometries and modify accordingly the definition of the space of discretisation $V$ through an admissible change of variables, the nature of the chosen characterisation $\{\sigma\}_{i \in [1, \dim V]}$ may be substantially modified.

Typically, a function of $V(\tilde{K})$ and its mapping in $V(K)$ may be characterised by different values of their respective degrees of freedom whenever the mapping consists in the simple change of variable. Indeed, the quantities estimated from the function by the degrees of freedom may not be preserved through the change of variables, especially when they are defined as moments or when they involve the function's normal components. In the former case, the definition of the degree of freedom would be sensitive to the change in the function shape and to the distortion of their support, while in the latter there would be sensitivity to the orientation of the cell. As a consequence, the characterisation nature defined through the degrees of freedom and thus the properties instilled to their dual basis functions would become cell-dependent. One has then to define a transformation between the elements of $V(\tilde{K})$ and $V(K)$ of which the chosen set of degrees of freedom are invariants.
As in our case we will focus on developing $H(\mathrm{div})$-conforming elements, we are specifically interested in preserving the conformity, that is preserving the continuity of normal components $q \mapsto q \cdot n$ and the values of the boundary flux integrals $q \mapsto \int_F q \cdot n d\gamma(x)$, where $\gamma$ is the path skimming a face $f$ of an polytope $K$, and $n$ the normal to the face $f$. A mapping that obtains those properties is e.g. the Piola transform, that is constructed out of any orientation-preserving bijective geometrical mapping.

**General definition of the Piola transform** Let us consider any reference shape $\tilde{K} \subset \mathbb{R}^d$ and a given cell $K \subset \mathbb{R}^d$. Let us also define a bijective geometrical mapping $F: \tilde{K} \to K$ that maps the geometry $\tilde{K}$ on $K$, denote its Jacobian by $J_F$ and the absolute value of its determinant by $|J_F| = |\det(J_F)|$. Then, provided that $V \subset L^2(K)$, one can construct a Piola transform as follows.

**Definition 2.2.2** *Piola transform*

Given any $\tilde{K}$, $K$ geometrical shapes contained in $\mathbb{R}^d$, functional subspaces $V(\tilde{K}) \subset L^2(\tilde{K})$, $V(K) \subset L^2(K)$, and a bijective mapping $F: \tilde{K} \to K$, the *Piola transform*

$$ \mathcal{P}: V(\tilde{K}) \to V(K) \\ \tilde{\phi} \mapsto \phi := \frac{1}{|J_F|} J_F \phi. \qquad (2.60) $$

is bijective, $H(\mathrm{div})$-conforming, orientation, and regularity preserving.

So defined, the same map also transfers to the gradient, where it particularly holds

$$ \nabla_x \phi = \frac{1}{|J_F|} J_F (\nabla_x \tilde{\phi}) J_F^{-1} \quad \text{and} \quad \nabla_x \cdot \phi = \frac{1}{|J_F|} \nabla_x \cdot \tilde{\phi}, $$

where $x$ represents the coordinate system of the cell $K$ and $x$ the one of the reference shape $\tilde{K}$. In addition, denoting by $\gamma$ and $\tilde{\gamma}$ the respective boundary measures, we have [43] for any function $\phi \in H(\mathrm{div}, K)$ and any $\tilde{v} \in H^1(K)$ that

$$ \left\{ \begin{aligned} \int_K \phi_i \cdot \phi_j \,dx &= \int_{\tilde{K}} \frac{1}{|J_F|} J_F \tilde{\phi}_i \cdot J_F \tilde{\phi}_j \,dx \\ \int_{\partial K} \phi \cdot n v \,d\gamma &= \int_{\partial \tilde{K}} \tilde{\phi} \cdot n \tilde{v} \,d\tilde{\gamma}, \\ \int_K \phi \cdot \nabla_x v \,dx &= \int_{\tilde{K}} \tilde{\phi} \cdot \nabla_x \tilde{v} \,dx \\ \int_K (\nabla_x \cdot \phi) v \,dx &= \int_{\tilde{K}} (\nabla_x \cdot \tilde{\phi}) \tilde{v} \,dx, \end{aligned} \right. $$

where $\tilde{v} = v \circ F^{-1}$ and where the relation $dx = |J_F|dx$ has been used. As a consequence, the degrees of freedom and basis functions can be naturally defined through the reference element $\tilde{K}$ on any distorted element $K$.
## Orientation preserving map for simplicial and quadrilateral elements

Let us now show how the Piola transformation is defined in practice by building an orientation-preserving transformation $F$ in the two simple case of the two-dimensional simplicial and quadrilateral geometries. In the case of simplicial geometries, such a transformation can be simply obtained by writing

$$F: \tilde{K} \to K \qquad \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \mapsto \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} x_{1,1} + (x_{1,2} - x_{1,1})x_1 + (x_{1,3} - x_{1,1})x_2 \\ x_{2,1} + (x_{2,2} - x_{2,1})x_1 + (x_{2,3} - x_{2,1})x_2 \end{pmatrix},$$

where the term $x_{i,j}$ represents the $i^{th}$ coordinate of the $j^{th}$ vertex of the cell $K$, and $i \in [1, 2]$, $j \in [1, n]$. Its Jacobian matrix then reads

$$J_F = \begin{pmatrix} \frac{\partial F_{x_1}}{\partial x_1} & \frac{\partial F_{x_1}}{\partial x_2} \\ \frac{\partial F_{x_2}}{\partial x_1} & \frac{\partial F_{x_2}}{\partial x_2} \end{pmatrix} = \begin{pmatrix} x_{1,2} - x_{1,1} & x_{1,3} - x_{1,1} \\ x_{2,2} - x_{2,1} & x_{2,3} - x_{2,1} \end{pmatrix},$$

and the definition of the Piola transform follows from (2.60). For quadrilateral shapes, the mapping further involves the term $x_1x_2$, reading

$$F: \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \mapsto \begin{pmatrix} x_2 \\ x_1 \end{pmatrix} = \begin{pmatrix} a_1 + b_1 x_1 + c_1 x_2 + d_1 x_1 x_2 \\ a_2 + b_2 x_1 + c_2 x_2 + d_2 x_1 x_2 \end{pmatrix},$$

where the coefficients $\{a_i, b_i, c_i, d_i\}_{i=1,2}$ are determined by the relations

$$A \begin{pmatrix} a_1 \\ b_1 \\ c_1 \\ d_1 \end{pmatrix} = \begin{pmatrix} 0 \\ 1 \\ 1 \\ 0 \end{pmatrix} \quad \text{and} \quad A \begin{pmatrix} a_2 \\ b_2 \\ c_2 \\ d_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 1 \\ 1 \end{pmatrix}, \quad A = \begin{pmatrix} 1 & x_1 & y_1 & x_1 y_1 \\ 1 & x_2 & y_2 & x_2 y_2 \\ 1 & x_3 & y_3 & x_3 y_3 \\ 1 & x_4 & y_4 & x_4 y_4 \end{pmatrix}.$$

Deriving then straightforwardly the Jacobian ends the definition of the Piola transform for quadrilateral elements.

<figure><img src="image_18.png" /><figcaption>**Figure 2.13.** Orientation-preserving mapping in the simplicial (grey area) and quadrangle (grey and white areas) cases.</figcaption></figure>
**Orientation preserving map for non-convex polytopal elements** The definition of an orientation preserving map for general non-necessarily convex polytopal elements is however much more delicate. Indeed, in addition to finding a mapping that is valid for any number of faces, one has to determine whether the normals orientation should be preserved or flipped, as depicted in Figure 2.14. This choice is highly dependent on the position of the element $K$ with respect to the axes and on the local convexities of the reference geometry $\tilde{K}$. One thus needs to redefine the concept of orientation for a non-convex element, and determine a transformation that depends on the cell shape, at least locally. Being a research topic in itself, we rather simply consider here mappings made of subtesselations of the polytope into simplicial shapes.

**Figure 2.14.** Difficulty in defining the notion of orientation when there is a change in the element local convexity.

<figure><img src="image_19.png" /><figcaption></figcaption></figure>
## 2.3 A brief numerical analysis toolbox

The investigated problem being described and the algebraical notions used in approximation theory being defined, we can now present the core principles of numerical analysis that allow the approximation of partial differential equations under a controlled error. We start in Section 2.3.1 by describing the general idea of a discretisation, and detail in Section 2.3.2 the fundamental features that a numerical scheme should have for the approximate solution to be trustworthy. In the last Section 2.4, we present an example of a scheme satisfying all the fundamental features, detailing in the process the Residual Distribution and Flux Reconstruction frameworks that define the scope of our subsequent work.

### 2.3.1 Basic concepts of discretisation

The numerical approximation of solutions to hyperbolic problems of the form (2.1a) relies on discretising the associated partial differential operator that is acting on the two levels that are the time evolution and the state-dependent system dynamics. The state usually not being constant over the domain $\Omega$, any approximation depends on a discretisation of the entire time space $\mathbb{R}_+ \times \Omega$. Furthermore, these two discretisations should be exploited by the approximation scheme in a way that ensures the resulting solution to be physically relevant in the sense of Definition 2.1.5, that is, being conservative and accurately representing the types and time-dependent locations of discontinuities.

#### 2.3.1.1 Discretisation approaches

In numerical analysis, any approximate solution to the problem (2.1a) is given for a bounded time interval $[0, T] \subset \mathbb{R}_+$, $T \in \mathbb{R}_+$ and is always finite dimensional. More precisely, the discrete solution lives in a variational space $W^h$ of finite dimension $|W^h|$ and is represented by the values $Y^h = \{Y_\sigma\}_{\sigma \in [1, |W^h|]}$ of a finite set of degrees of freedom $\{\sigma\}_{\sigma \in [1, |W^h|]}$ from which a global solution approximation $Y^h: \Omega \to W$ can be reconstructed. The value of those degrees of freedom are obtained and updated through numerical schemes, whose construction can follow or combine several flux discretisation approaches.

<figure><img src="image_20.png" /><figcaption>**Figure 2.15.** Illustration of the characterisation of a discrete solution living in a finite dimensional variational space.</figcaption></figure>
**Strategies to handle the spatio-temporal interplay** Constructing those numerical schemes reduces to building a discrete partial differential operator that approximates both the time derivative and the spatial divergence.

To this aim, a first approach treats simultaneously the two discretisation aspects by considering the time as an extra variable to be included in the divergence operator. Entering the class of *time-space* methods, any discretisation is then built on an augmented flux and is therefore performed on a $d + 1$ dimensional space. If this type of approximation is comprehensive and suggests a global robustness of the solution, it leads however to the solving of large algebraic systems. Furthermore, the time discretisation being enclosed in a global setting that has been entirely prescribed, the set of points in the time-space on which the solution is available is pre-determined. Hence, there cannot be any time step adaptation with respect to possibly rapid changes in the solution, behaviour that is commonly encountered in the scope of hyperbolic problems due to the development of discontinuities and presence of heavy distortions. Despite promising tentatives [44, 45], this approach therefore seems unsuitable to our case.

A second approach, considered here, uses independent discretisations for the time and for the dynamic of the system. The solution there lives in a variational space that is built on an approximation of the domain $\Omega$ and is represented through a set of time-dependent degrees of freedom values that are updated every *time step*. The coupling between the temporal and spatial behaviour being then only active a specific points of the time space, this formulation already represents a first approximation of the problem (2.1a). Hence, such formulations are referred to as *semi-discrete*, and read

$$ \frac{dY^h}{dt} + D_{\mathcal{F}}^h(Y^h) = 0, \qquad (2.61) $$

where $D_{\mathcal{F}}^h: W \to \mathbb{R}^q$ is a given discrete operator that approximates the flux's divergence and encompasses the boundary conditions. Considering this type of formulations, it is then possible to define the expression of $D_{\mathcal{F}}^h$ independently from the approximation of the plain temporal derivative, their conjoint use providing anyhow a comprehensive numerical scheme.

**Strategies for approximating the differential operator** Like in the approximation of the time-dependent problem, the discretisation of the differential operator $D_{\mathcal{F}}^h$ can also be investigated from different perspectives. The two main approaches define the classes of collocation and grid-based methods.

Collocation methods [46, 47] structure the solution from the definition of basis function themselves and are particularly popular for their ability to handle multivariate interpolation over discontinuous fields. Further benefiting from a remarkable accuracy, they however require an extreme carefulness in the long run due to the high conditioning of their associated system.
More robust and achieving more easily accurate behaviours in long term simulations, gridded methods rely on a predefined spatial structure of the discrete solution. The connectivity between discrete points playing then a key role in the approximation of the flux's divergence, the use of Cartesian grids has been shown to be very efficient [48, 49]. However, they fail in practical applications where complex geometries are involved, as they are usually not boundary conforming and therefore loose accuracy if no local and specialized adaptation is made. We therefore choose here to focus on methods based on unstructured grids, conferring more flexibility towards the computational domain's geometry and allowing to define customisable spaces of discretisation.

### 2.3.1.2 Unstructured spatial discretisation

Defining a numerical scheme that approximates the system dynamics on an unstructured grid relies on three discretisations, one geometrical structuring the spatial computational domain, one variational representing the discrete solution, and one approximating the differential operator itself.

**Geometrical discretisation** Discretising the problem geometry by means of an unstructured grid comes down to defining a covering of the domain $\Omega$ by a finite number of cells $K \subset \Omega$ whose shape can be *curvilinear* or *polytopal*. The collection of all the cells, denoted by $\tau^h = \{K\}_K$, then forms a *mesh* whose frame determines the connectivity between the vertices of all the cells. Belonging to this frame, each generic face $f$ of a cell's boundary $\partial K$ defines either an interface with a neighbouring cell $K^-$ or approximates locally the domain's boundary $\partial \Omega$. We then denote by $\partial\tau^h$ the set of faces that approximates the domain boundary, and by $\Gamma$ a generic boundary face belonging to $\partial\tau^h$. All in all, the space covered by the mesh $\tau^h$ defines the *computational domain* and will be denoted $\Omega$.

<figure><img src="image_21.png" /><figcaption>**Figure 2.16.** Example of a boundary-conforming hybrid polygonal mesh providing an approximation of $\Omega$ and a structure on which to define local variational spaces.</figcaption></figure>
If any covering defines a mesh, one should still pay attention to its *quality*, that is, the ability of the mesh to support the discretised solution so that the approximation error due to the tessellation is homogeneous throughout the computational domain Ω. This mesh quality depends on the overall layout of the cells vertices, particularly impacting the homogeneity of the cells' size and aspect ratio. Given a quality mesh, one can then define the *mesh size h* as the maximum radius of each of the cells' circumscribed sphere. Lastly, let us point out that a mesh containing *hanging nodes*, that is, containing cells for which there exists a vertex that intersects with a face of another cell, are called *non-conformal* and are usually mishandled by numerical schemes. Therefore, it may be interesting to consider *hybrid meshes*, composed of cells that do not necessarily share the same topology.

<figure><img src="image_22.png" /><figcaption>**Figure 2.17.** Example of a pathological hybrid mesh containing quads and simplices, showing hanging nodes and heavy distortions. Note that the exact same mesh seen as polygonal one would not have any hanging nodes, being then absorbed in the cell geometries. Namely, the pictured polygonal mesh would containing simplicial, quadrilateral, and pentagonal cell shapes (the two vertically oriented cells).</figcaption></figure>

In this thesis, we are particularly interested in using non-convex poly-
topal hybrid meshes, allowing flexibility in both the geometry approximation
and in the scheme design. Indeed, compared to the widely used simplicial
or quadrilateral meshes, polytopal meshes allow an automatic inclusion of
hanging nodes and are very efficient in mesh refinement. Allowing further
non-convexity then alleviates coarsening issues and confers a robustness to
distortions, making them suitable for developing mesh-adaptation techniques.
However, the definition of variational spaces that can represent the discrete
solution over any polytopal shape with a same representation quality is ex-
tremely delicate, as it will be observed in Chapter 3.
**Variational discretisation** After defining the mesh, one can specify on each cell K a finite dimensional variational space W(K) containing the local approximate solution Y^K. The approximation nature is then defined in W(K) through the choice of a set of degrees of freedom Σ_{W(K)} := {σ}σ∈[1,|W(K)|] whose values Y^K = {Y_σ^K}σ∈Σ_{W(K)} characterises uniquely the discrete solution on the cell K, and from which the approximate solution Y^K can be reconstructed. Indeed, the triplet (K, W(K), Σ_{W(K)}) thus forming an element, one can retrieve the basis functions of W(K) that are dual to the set Σ_{W(K)} and provide a direct local representation Y^K: K → R^q of the approximate solution given the values Y^K. This local representation can then be used in turn to define the global approximate solution reconstruction

$$
\begin{equation} \tag{2.62}
\begin{aligned}
Y^h: \Omega &\longrightarrow \mathbb{R}^q \\
x &\longmapsto Y^K(x), \quad K \text{ s.t. } x \in K,
\end{aligned}
\end{equation}
$$

living in W^h = {W(K), K ∈ τ^h}, but also to design the discrete differential operator, typically allowing flux evaluations at quadrature points that are not necessarily matching the pointwise degrees of freedom, if even any.

A classical choice is to consider W(K) as a finite dimensional polynomial space and to endow it with a nodal basis. Intuitive discretisation setting, the characterisation of the solution there consists in point-values, while its representation is simply given by a polynomial interpolation driven by the shape of the nodal basis functions.

**Approximation of the system dynamics** To finalise the discretisation of the system dynamics, it remains to define the discrete operator $D_F^h$ so that it provides a relevant approximation of the continuous differential operator. Especially, in light of the discussion of Section 2.1, one should derive a discrete expression starting from a conservative formulation. Here again, two main perspectives are possible. One may indeed consider either the strong formulation

$$
\left\{
\begin{array}{l}
\nabla \cdot \mathcal{F}(Y^h) = 0 \\
Y^h|_{\partial\Omega} = g
\end{array}
\right.
\qquad (2.63)
$$

and approximate the flux's divergence by evaluating the flux itself across the cells, or consider the steady analogue of the weak formulation (2.42) split over the mesh $\tau^h$ i.e.

$$
\sum_{K \in \tau^h} \int_K \varphi (\nabla \cdot \mathcal{F}(Y^h)) \, dx + \sum_{\Gamma \in \partial \tau^h} \int_{\Gamma} \varphi g \, dx = 0, \quad \forall \varphi \in C_0^\infty(K), \forall K \in \tau^h \quad (2.64)
$$

approximating then each surface and volume integrals. There, it may also be
suitable to emphasise the impact of the flux behaviour on the boundary by
reformulating the inflow boundary conditions, yielding

$$ \sum_{K \in \tau^h} \int_K \varphi(\nabla \cdot \mathcal{F}(Y^h)) \, dx + \sum_{\Gamma \in \partial \tau^h} \int_{\Gamma} \varphi(\mathcal{F}(Y^h) \cdot n - F_b(Y^h, g)) \, dx = 0 \quad (2.65) $$

for any $K \in \tau^h$, any $\varphi \in C_0^\infty(K)$, and where $F_b$ is a flux consistent with the relation $\frac{d\mathcal{F}}{dY} \cdot n < 0$. The discrete operator defined, one can then investigate its properties and the quality of the approximation it provides given a mesh size.

**Note** In this thesis, we are considering a semi-discrete discretisation and only focus on defining the operator $D_F^h$. We therefore drop in below developments the time considerations whenever they are not impacting the spatial dynamic. When they do, the related statements have to be understood in light of the temporal approximation, chosen here as the already developed Deferred Correction framework [50].

## 2.3.2 Essential properties of a numerical scheme

Let us now present the essential features that a numerical scheme should possess for the approximate solution to be reliable. We start by stating the properties allowing trust in the discrete solution, before deriving estimates on the scheme performance.

**Reliability of the approximate solution** We first consider the essential properties of the numerical scheme that is defined by the expression of the discrete differential operator $D_F^h$ and whose application to any boundary data $g$ gives the values of $Y^K$ on each element $(K, W(K), \{\sigma\})$, $K \in \tau^h$, consequently providing the approximate solution reconstruction $Y^h$.

A first natural property to ask when aiming for the relevance of the approximate solution is the convergence, expressing that the approximate solution tends to the exact one as the mesh is refined.

<figure><img src="image_23.png" /><figcaption>Figure 2.18. Illustration of the convergence property.</figcaption></figure>

### Definition 2.3.1 Convergence

A numerical scheme is said to be convergent when the error

$$ \lim_{h \to 0} \|Y - Y_h\| \to 0 $$

for a given $\|\cdot\|$ that is suited to estimate the deviation of the quantities expressed in $W^h$ from a possibly discontinuous reference one.

A second property that is necessary for the solution to be reliable is called consistency and assesses the ability of the scheme expression to approximate the continuous differential operator.
**Definition 2.3.2** *Consistency*

Given the steady analogue of the continuous problem (2.1a) and its solution by $Y$, a numerical scheme is said to be consistent if the discrete differential operator $D_{\mathcal{F}}^h$ converges towards the continuous operator $\nabla \cdot \mathcal{F}$ when $h \to 0$ in the norm the discretisation space of the operator is endowed with.

Finally, the numerical scheme must also be robust, that is, not amplifying perturbations to a greater extent than the continuous problem would.

<figure><img src="image_24.png" /><figcaption>**Figure 2.19.** Illustration of the stability property.</figcaption></figure>

**Definition 2.3.3** *Stability*

A numerical scheme is said to be stable if and only if there exists a unique solution $Y^h$ to the problem $D_{\mathcal{F}}^h(Y^h) = 0$, and that this solution depends continuously on the changes in the given data, here the boundary conditions.

Whenever a numerical scheme fulfils those three properties, the approximate solution is *a-priori* trustworthy. Therefore, seen in the context of hyperbolic systems, any convergent, consistent and stable operator $D_{\mathcal{F}}^h$ would provide a reliable representation of the state dynamics to the temporal scheme tackling the plain temporal derivative in (2.61). Furthermore, when the comprehensive time-space scheme (2.61) tackling the problem (2.1a) is also consistent and convergent, it has been shown [51] that the approximate solution converges to a weak solution of the system (2.1a) as the mesh size tend to zero, statement known as the Lax-Wendroff theorem.

**Theorem 2.3.1** *Lax-Wendroff convergence theorem*

Let $Y(0, \cdot) \in L^\infty(\Omega)$ be any bounded initial condition, $Y^h: \mathbb{R}_+ \times \Omega \to \mathbb{R}^q$ be an approximate solution to the problem (2.1a) obtained through a consistent and conservative scheme, and $Y \in L^2(\mathbb{R}_+ \times \Omega)$ be a square integrable function. Then, whenever

* $Y^h$ is uniformly bounded for any reasonable mesh size, i.e.

$\exists h_{\max} \in \mathbb{R}_+^* \text{ s.t. } \forall h \in (0, h_{\max}), \|Y^h\|_{L^\infty(\mathbb{R}_+ \times \Omega)} \le C \in \mathbb{R},$

* the solution converges almost everywhere to the exact solution, i.e.

$\lim_{h \to 0} \|Y - Y^h\|_{L_{\text{loc}}^2(\mathbb{R}_+ \times \Omega)} = 0,$

the function $Y$ is a weak solution of the problem (2.1a).

However, while the theorem ensures the obtained solution to be a weak solution of the problem (2.1), it does not guarantee its physical relevance.
**Physical relevance of the approximate solution** As seen in Section 2.1, the lack of uniqueness of weak solutions to the problem (2.1a) can be compensated by singling out the physically relevant solutions. This criterion transfers to the discrete level, where we ask the numerical scheme to fulfil the discrete analogue on the entropy relation (2.53) for a given entropy functions and flux pair.

**Definition 2.3.4** *Entropy stable and conservative schemes*

A numerical scheme is said to be entropy stable with respect to the entropy pair $(S, G)$ if for the discrete approximate solution $Y^h$ and its reconstruction $Y^h: \mathbb{R}_+ \times \Omega \to \mathbb{R}^q$ it holds

$$ \frac{d}{dt}(S(Y^h)) + D_{G_n}^h(Y^h) \le 0, \qquad (2.66) $$

where $G_n$ is a numerical entropy flux consistent with $G$ and $D_{G_n}^h$ is the discrete flux divergence operator $D_F^h$ defined on the numerical entropy flux $G_n$. The scheme is further said to be *entropy conservative* when the equality holds.

Once a scheme is stable, consistent, conservative and entropy conservative, we are ensured that the approximate solution converges to a weak solution of the system (2.1a) that is physically relevant in the sense of the considered entropy pair. One can however wonder about the performance of a scheme, that is, how much effort should one put in the geometrical meshing to obtain a meaningful solution.

**Performance of a numerical scheme** The performance of a scheme can be assessed by estimating the overall error of the approximate solution given the size of the mesh it has been computed from. The estimate retrieved from the scheme's definition indeed not only indicates how fast the solution will converge to a weak solution of the system (2.1a), but also provides the expected quality of the approximate solution for a given mesh size.

**Definition 2.3.5** *Convergence order and method accuracy*

Provided a smooth initial condition $Y(0, \cdot)$, the scheme's order of accuracy is given by the highest value of $p$ such that

$$ \|Y - Y_h\| \le C h^p, \quad C \in \mathbb{R}, \qquad (2.67) $$

for a given norm $\|\cdot\|$ that is suited to estimate the deviation of the quantities expressed in $W^h$ from a possibly discontinuous reference one. The constant $C$ is here independent of the mesh size and of the approximate solution, but may depend on the initial and boundary conditions, and even on the exact solution itself. The relation (2.67) is therefore an *a-priori estimate*.
## 2.4 Some classical numerical schemes

With those theoretical principles in mind, let us now practically construct a numerical scheme that tackles two-dimensional problems of the type (2.1a) while being entropy stable, conservative, and applicable to on any polygonal mesh. In this perspective, we consider the formulation (2.61) in which the plain temporal derivative is approximated by a Deferred Correction approach [50], and focus only on deriving a discrete expression for the flux divergence operator $D_F^h$. We there choose to employ the entropy conservative Flux Reconstruction scheme developed in [52], repeating here its construction to highlight its range of application and deriving conditions under which the scheme can be applied on any unstructured polygonal mesh. This scheme being built by exploiting the existing links between the Flux Reconstruction schemes, the Discontinuous Galerkin schemes, and the Residual Distribution framework, we start by describing those three concepts in their respective Sections 2.4.1, 2.4.2, and 2.4.3, before deriving the entropy conservative Flux Reconstruction scheme in Section 2.5.

### 2.4.1 Discontinuous Galerkin methods

Introduced by Reed and Hill in [53] and later developed by Cockburn and Shu [54], the Discontinuous Galerkin method is very popular in applications involving fluid dynamics or wave propagation phenomenon, as e.g. meteorology, oceanography, and seismology. This popularity is mainly due to the flexibility conferred to the design of both the flux discretisation and the solution representation, in particular by not requiring their overall continuity.

**Discretisation setting** Freed from the global continuity requirement, the approximate solution is then considered with a global $L^2(\Omega)$ regularity, living in the finite dimensional *trial* subspace

$$W^h = \{ w \in L^2(\Omega, \mathbb{R}^q) \text{ s.t. } \forall K \in \tau^h, w|_K \in W(K) := (\mathbb{P}_K(K))^q \}.$$

The system dynamics is then approximated starting from the conservative variational formulation (2.65), where the *test space* is approximated by

$$V^h = \{ v \in L^2(\Omega) \text{ s.t. } \forall K \in \tau^h, v|_K \in V(K) := \mathbb{P}_K(K) \},$$

within which all test functions have a compact support and are locally regular. Hence, one can define local sets of degrees of freedom $\Sigma_{V(K)} := \{\Sigma_{V(K)}, K \in \tau^h\}$ that characterise the functions living in a local space $V(K)$, and approximate the system dynamics (2.65) on each element $(K, V(K), \Sigma_{V(K)})$.
**Formulation of the approximated dynamic** Given the above discretisation choices, the approximation of the formulation (2.65) reduces to asking

$$ \sum_{K \in \tau^h} \int_K \varphi (\nabla \cdot \mathcal{F}(Y^K)) dx + \sum_{\Gamma \in \partial \tau^h} \int_{\Gamma} \varphi (\mathcal{F}(Y^K) \cdot n - F_b(Y^K, g)) dx = 0 $$

only for any function $\varphi$ belonging to $V^h$. Integrating then by parts, it comes

$$ - \sum_{K \in \tau^h} \left( \int_K \nabla \varphi \cdot \mathcal{F}(Y^K) dx + \int_{\partial K} \varphi \mathcal{F}(Y^K) \cdot n dx \right) + \sum_{\Gamma \in \tau^h} \int_{\Gamma} \varphi (\mathcal{F}(Y^{\Gamma}) \cdot n - F_b(Y^{\Gamma}, g)) dx = 0, $$

making the approximation of the differential operator $\mathcal{D}_{\mathcal{F}}^h$ entirely local to one cell, the dynamic being entirely represented by the sole flux's expression. Thus, to complete the approximation one has just left to define a discrete numerical flux $F_n$ that is consistent with the continuous one and that considers the impact of the neighbouring behaviour to the local dynamic through the element boundary. The final formulation thus reads

$$ - \sum_{K \in \tau^h} \left( \int_K \nabla \varphi \cdot \mathcal{F}(Y^K) dx + \int_{\partial K} \varphi F_n(Y^K, Y^{K-}) dx \right) + \sum_{\Gamma \in \tau^h} \int_{\Gamma} \varphi (\mathcal{F}(Y^{\Gamma}) \cdot n - F_b(Y^{\Gamma}, g)) dx = 0 \quad (2.68) $$

for any $\varphi \in V^h$ and any $K \in \tau^h$, the actual definition of the scheme being then given by the choice of the two numerical fluxes' expressions. Furthermore, while the boundary flux $F_b$ has to be consistent with $\frac{d\mathcal{F}}{dY} \cdot n < 0$, the expression of $F_n$ shall be consistent with $\mathcal{F}$ while ensuring the scheme's conservation. One therefore asks

$$ \left\{ \begin{array}{l} F_n(Y^K, Y^K) = \mathcal{F}(Y^K) \cdot n \\ F_n(Y^K, Y^{K-}) = -F_{-n}(Y^{K-}, Y^K), \end{array} \right. $$

under which condition a typical choice is given by

$$ F_n(a, b) = \frac{1}{2} (\mathcal{F}(a) \cdot n + \mathcal{F}(b) \cdot n) - \Delta(\mathcal{F}, a, b, n), \quad (2.69) $$

where $\Delta$ represent a numerical dissipation operator, see e.g. [21] for possible choices. The scheme so defined, one can obtain its practical expression by decomposing any function $\varphi$ on the basis $\{\varphi_{\sigma}\}_{\sigma \in [1, |V(K)|]}$ of the space $V(K)$ that is dual to the set of degrees of freedom $\Sigma_{V(K)}$. A direct representation of the approximate solution can equally be obtained by endowing the space $W(K)$ with a set of degrees of freedom $\Sigma_{W(K)} := \{\sigma\}_{\sigma \in [1, |W(K)|]}$ that are compactly supported in $K$. The above formulations are then made entirely explicit, and the discrete solution is entirely characterised in each cell $K$ by its local degrees of freedom values. Endowing further each local space $W(K)$ with the basis that is there dual to $\Sigma_{V(K)}$ then provides a global reconstruction of the solution as a piecewise polynomial that is possibly discontinuous across the cells faces.
## 2.4.2 Flux reconstruction schemes

Introduced in particular by Huynh [19] and extended by Jameson [55], the Flux Reconstruction approach mainly differs from the Galerkin methods by relying entirely on the estimation of pointwise quantities rather than on the approximation of local integrals. It also enforces the continuity of the discrete flux across the element boundaries with a correction function, making its definition very flexible. Furthermore entirely discrete and naturally conservative, those methods require however specific treatments to be stable.

### 2.4.2.1 Discretisation framework

Easing the technique presentation, we focus on the two spatial dimensional case ($d = 2$) in which the geometrical discretisation provides a conformal mesh $\tau^h$ of the domain $\Omega \subset \mathbb{R}^2$. There constructed, any Flux Reconstruction technique relies on the mesh to structure the discrete solution and to define a consistent approximation of the dynamics based on the strong conservative formulation (2.63), requiring in particular a direct representation of both the solution and the flux.

**Solution discretisation** The solution is simply represented in each cell $K \in \tau^h$ by the multivariate polynomial $Y^K: K \to \mathbb{R}^q$ of degree $k \in \mathbb{N}$ that interpolates given point-values accordingly to the associated nodal basis of $\mathbb{P}_k(K)$. The characterisation of the approximate solution is therefore immediately defined by the dual set of $N_Y := \frac{1}{2}(k+1)(k+2)$ degrees of freedom $\Sigma_Y = \{\mathcal{G}: q \mapsto q(x_m), x_m \in \mathcal{X}_Y\}$, where $\mathcal{X}_Y = \{x_m \in K\}_{m \in [1, N_Y]}$ represents a set of disjoint evaluation points whose associated degrees of freedom values $Y^K := \{Y_{\mathcal{G}}^K\}_{\mathcal{G} \in \Sigma_Y}$ determine the discrete solution in $(\mathbb{P}_k(K))^q$. Hence, denoting by $\{l_{\mathcal{G}}: K \to \mathbb{R}\}_{\mathcal{G} \in \Sigma_Y}$ the set of Lagrangian functions that is built from the points $\mathcal{X}_Y$, the local approximation of the solution naturally reads

$$Y^K(x) = \sum_{\mathcal{G} \in \Sigma_Y} Y_{\mathcal{G}}^K l_{\mathcal{G}}(x). \qquad (2.70)$$

Each local representation being compactly supported on $K$, the approximate solution is naturally defined over $\Omega$ as the piecewise polynomial given straightforwardly through (2.62). However, as the interpolation nature generally prevents the reconstruction to be continuous on the entire computational domain regardless the spatial distribution of the points $\mathcal{X}_Y$ in each element $K$, the approximate solution only benefits from a global $L^2(\Omega)$ regularity, hence living in the space

$$W^h = \{w \in L^2(\Omega, \mathbb{R}^q) \text{ s.t. } \forall K \in \tau^h, w|_K \in W(K) := (\mathbb{P}_k(K))^q\}.$$
**Flux discretisation** The flux $\mathcal{F} = (\mathcal{F}_1, \mathcal{F}_2)^T$ is discretised similarly, each approximate flux function living in the finite dimensional variational space

$$V^h = \{v \in L^2(\Omega) \quad \text{s.t.} \quad \forall K \in \tau^h, v|_K \in V(K)\}$$

that is built on given local discretisation spaces $V(K)$. The global approximation of the flux is therefore given by

$$F^h: \Omega \longrightarrow \mathbb{R}^q \\ x \longmapsto F^K(x), \quad K \text{ s.t. } x \in K,$$

where the term $F^K: K \to \mathbb{R}^q \times \mathbb{R}^q$ denotes a local *flux reconstruction* that is compactly supported on $K$. Defined in a corrective approach, this reconstruction simply writes

$$F^K := F^{K,D} + F^{K,C} = \begin{pmatrix} F_1^{K,D} \\ F_2^{K,D} \end{pmatrix} + \begin{pmatrix} F_1^{K,C} \\ F_2^{K,C} \end{pmatrix}, \quad (2.71)$$

where $F^{K,D} = (F_1^{K,D}, F_2^{K,D})^T$ denotes a direct approximation of the local flux functions and $F^{K,C} = (F_1^{K,C}, F_2^{K,C})^T$ its correction terms whose expressions complete the definition of the numerical scheme. Hence, each local discretisation space can be written as

$$V(K) = V^D(K) + V^C(K)$$

for $V^D(K)$ and $V^C(K)$ two finite dimensional variational spaces in which any initial flux approximation and correction can be expressed, respectively. While this correction procedure appears to be a natural approach, those two spaces are however not necessarily in direct sum, making the design of the correction function's characterisation subject to careful considerations.

### ▷ Initial flux approximation.
In classical Flux Reconstruction approaches, which the approach [56] that we follow here is part of, the initial approximation of each flux function $F_i^{K,D}$, $i \in [1, d]$, is typically sought as for the solution, i.e. considered in the space

$$V^D(K) = (\mathbb{P}_k(K))^d$$

and defined as the interpolation of given point-values by means of the associated nodal basis of $\mathbb{P}_k(K)$. Namely, denoting by $N_F$ the dimension of $\mathbb{P}_k(K)$ and considering a same set of disjoint points $\varkappa_F = \{x_j \in K\}_{j \in [1, N_F]}$ where all the discrete flux functions are given by their respective values $F_i^{K,D} := \{F_{i,\sigma}^{K,D}\}_\sigma$ in $(\mathbb{P}_k(K))^q$, one considers the set of pointwise degrees of freedom $\Sigma_F := \{\sigma: q \mapsto q(x_j), x_j \in \varkappa_F\}$ whose vectorial values characterise the approximate flux functions $F_i^K$, $i \in [1, 2]$. The initial flux representation then writes

$$F_i^{K,D}(x) = \sum_{\sigma \in \Sigma_F} F_{i,\sigma}^K l_{\sigma}(x), \quad i \in [1, 2], \quad (2.72)$$

for $\{l_\sigma: K \to \mathbb{R}\}_{\sigma \in \Sigma_F}$ the set of Lagrangian functions built from the points $\varkappa_F$. Straightforward, this representation is however not continuous across the cells, leading to the first approximation $F^{K,D}$ to be called *discontinuous flux*.
▷ **Correction of the initially approximate flux.** The correction term $F^{K,C}$ is then defined so as to enforce the cross-cell continuity, and is hence usually thought as being characterised on the boundary through quantifiers that are directly supported on a given face alone, thus easily and locally connecting the neighbouring cells. The overall boundary characterisation, supported in practice by a discrete stencil that we will denote $\varkappa_C$, may however not be enough to characterise the correction function entirely on $K$. Indeed, its expression may benefit from an extra freedom to be tuned so that the obtained function shape propagates the correction amount within the cell under given inner constraints.

Hence intricate, the design of a suitable general shape for the correction function and the choice of relevant determination constraints are crucial to the well-definition of the comprehensive numerical scheme. However, even though the determination constraints can be designed directly upon the properties one wants to instil to the numerical scheme, as shown in the second next paragraph of this section, the design of the general shape of the correction function requires to construct a suitable discretisation space $V^K(K)$ and to define a characterisation method through the choice of a set of degrees of freedom $\Sigma_C$. Requiring careful considerations and not straightforward, those last two aspects will be the object of the entire Section 2.5, in turn motivating the entire Chapter 3.

▷ **Validity of the flux discretisation.** While the representation of the comprehensive corrected flux approximation (2.71) may seem intuitive, its characterisation method has to be carefully considered. Indeed, as the spaces $V^K(K)$ and $V^K(K)$ are not necessarily in direct sum, there may be an interplay between the two characterisation methods when considering the comprehensive determination of any corrected approximate flux functions in $V(K)$. In other words, the set of degrees of freedom $\Sigma_F + \Sigma_C$ may not coincide with the set of degrees of freedom $\Sigma_{V(K)}$ characterising the comprehensive corrected approximate flux functions (2.71) in $V(K)$.

As a consequence, one has to take care that the chosen characterisation $\Sigma_C$ determining the correction function in $V^K(K)$ combined with the characterisation $\Sigma_F$ of $V^K(K)$ ensure that the associated comprehensive characterisation $\Sigma_{V(K)}$ of the corrected approximate flux functions (2.71) is unisolvent for $V(K)$. If it may *a-priori* be delicate, it is however enough to ensure that there exists subsets $\Sigma_C \subset \Sigma_C$ and $\Sigma_F \subset \Sigma_F$ that both characterise the subset $V^K(K) \cap V^D(K)$, possibly in a different way, as shown in the following proposition.
**Proposition 2.4.1** *Sufficient conditions for unisolvence*

If there exists subsets $\Sigma_C \subset \Sigma_C$ and $\Sigma_F \subset \Sigma_F$ that both characterise the subset $V^C(K) \cap V^D(K)$, even in a different way, then

* the comprehensive characterisation $\Sigma_{V(K)}$ that is associated to the independent characterisations $\Sigma_F$ and $\Sigma_C$ of $V^{D,K}$ and $V^{C,K}$, respectively, is unisolvent in $V(K)$.

* the corrected approximate function $F^K$ given by (2.71) is uniquely determined in $V(K)$

*Proof.* Let us denote by $\{l_\sigma\}_{\sigma \in [1, N_F]}$ the set of basis functions that is dual to $\Sigma_F$ in $V^D(K)$ and by $\{\psi_j\}_{j \in [1, N_C]}$ the set of basis functions that is dual to $\Sigma_C$ in $V^C(K)$. Then, one can write any function $F^K \in V(K)$ in the form

$$F^K = \sum_{\sigma \in \Sigma_F} F_\sigma^K l_\sigma + \sum_{j \in [1, N_C]} \alpha_j \psi_j.$$

Expressing then explicitly the parts of the function that belong to $V^C(K) \cap V^D(K)$, one retrieves

$$F^K = \sum_{\substack{\sigma \in \Sigma_F, \\ \sigma \notin \Sigma_{F \cap}}} F_\sigma^K l_\sigma + \left( \sum_{\sigma \in \Sigma_{F \cap}} F_\sigma^K l_\sigma + \sum_{\substack{j \in [1, N_C], \\ \sigma_j \in \Sigma_{C \cap}}} \alpha_j \psi_j \right) + \sum_{\substack{j \in [1, N_C], \\ \sigma_j \in \Sigma_{C \cap}}} \alpha_j \psi_j.$$

There, one can then consider $\Lambda$ the matrix that transfers the set of basis functions that are dual to $\Sigma_C$ in $V^C(K) \cap V^D(K)$ to the set of basis functions that are dual to $\Sigma_F$ in $V^C(K) \cap V^D(K)$. Denoting then $(\Lambda \alpha)_\sigma$ the value of the correction function that is associated to the function $l_\sigma$ in its decomposition on the new basis representation, one can write

$$F^K = \sum_{\substack{\sigma \in \Sigma_F, \\ \sigma \notin \Sigma_{F \cap}}} F_\sigma^K l_\sigma + \sum_{\sigma \in \Sigma_{F \cap}} (F_\sigma^K + (\Lambda \alpha)_\sigma) l_\sigma + \sum_{\substack{j \in [1, N_C], \\ \sigma_j \in \Sigma_{C \cap}}} \alpha_j \psi_j \\ = \sum_{\sigma \in \Sigma_{V(K)}} \beta_\sigma \Psi_\sigma,$$

for the basis of $\Psi = \{\Psi_\sigma\}_{\sigma \in \Sigma_{V(K)}}$ of $V(K)$ and associated coefficients $\{\beta_\sigma\}_{\sigma \in \Sigma_{V(K)}}$ obtained through a straightforward identification between the two expressions. Hence, as by construction the basis functions of the set $\{l_\sigma\}_{\sigma \in F} \cup \{\psi_j\}_{j \in [1, N_C], \sigma_j \in \Sigma_C}$ are linearly independent and of dimension $\dim V(K)$, the function $F^K$ is uniquely defined in $V(K)$, and the characterisation of the correction function viable in the more general context of the comprehensive discretisation.

Note however that the overall characterisation is not necessarily nodal, and that the identified basis $\Psi_\sigma$ does not necessarily correspond to the dual basis of the comprehensive characterisation $\Sigma_{V(K)}$ of $V(K)$. Indeed, even though the value of $\alpha$ is still to be considered at point values, its transfer through $\Lambda$ takes away its spatial connection, the initial point-value $(x_j, F_\sigma^K)$ being simply modulated from non-point-value information.
**Dynamic evolution** To update the state variables, the Flux Reconstruction technique uses the strong formulation (2.63) in which the solution approximation and the flux functions are directly employed. Hence, the discrete solution $\Upsilon^K$ to the problem (2.1) sees its state values $\{\Upsilon_{\mathcal{S}}^K\}_{\mathcal{S} \in \Sigma_{\gamma}}$ updated independently at their corresponding point in $\varkappa_{\gamma}$ by a straightforward evaluation of the approximated pointwise dynamic expressed in (2.72) and (2.70). Any updated value $\Upsilon_{\mathcal{S}}^K, \mathcal{S} \in \Sigma_{\gamma}$ is therefore practically obtained through the relation

$$ \frac{d\Upsilon_{\mathcal{S}}^K}{dt} = -(\nabla \cdot F^K)(x_m) \\ = -(\nabla \cdot F^{K,D})(x_m) - (\nabla \cdot F^{K,C})(x_m) \\ = -\left(\sum_{\sigma \in \Sigma_{\mathcal{F}}} F_{1,\sigma}^K \frac{\partial l_{\sigma}}{\partial x_1}(x_m) + \sum_{\sigma \in \Sigma_{\mathcal{F}}} F_{2,\sigma}^K \frac{\partial l_{\sigma}}{\partial x_2}(x_m)\right) - (\nabla \cdot F^{K,C})(x_m), \quad (2.73) $$

where the specific point $x_m \in \varkappa_{\gamma}$ is the one supporting the definition of the degree of freedom $\mathcal{S}$ in its pointwise expression, and where the term $(\nabla \cdot F^{K,C})(x_m)$ represents the divergence of the correction function evaluated at the point $x_m \in \varkappa_{\gamma}$, to be further discretised accordingly to the definition of the correction function. The choice of the correction function shape along with the selection of the flux and discrete solution points $\varkappa_{\gamma}$ therefore jointly define the flux reconstruction scheme.

<figure><img src="image_25.png" /><figcaption>**Figure 2.20.** Update procedure for general Flux Reconstruction schemes. Left: Example of a discrete point layout. Right: Update procedure in the one-dimensional case (analogous in higher dimensions).</figcaption></figure>

As those discretisation choices are not bound a-priori to any constraint other than requiring the existence of a pointwise reconstruction and ensuring the cross-cell continuity of the flux, the design of flux reconstruction schemes is extremely flexible and encompasses most of the classical schemes as for example Discontinuous Galerkin [57] and Spectral Element Methods [58] in their fully discretised formulations. However, the schemes properties only rely on the articulation of the discrete flux and solution points with the definition of the correction function.
Each approximation design is therefore extremely delicate, in particular when tackling hyperbolic problems of the form (2.1). Indeed, there, although the conservation properties can be readily obtained by imposing the continuity of the flux, the scheme long-term stability and the convergence of the approximate solution towards a physically relevant one remain to be ensured from the definition of the correction function alone. Focusing on the former aspect, a flux reconstruction scheme guaranteeing linear stability has already been proposed in [56], and non-linearity aspects have been investigated in [59]. We choose to focus here on the latter aspect and detail a choice of correction function and discretisation points $\kappa_F$ that preserves the consistency, enforces the conservation, and guarantees entropy stability regardless the chosen points $\kappa_Y$ where to define the discrete solution.

### 2.4.2.2 A possible definition of a correction function

For the sake of the example and in view of later use, let us now derive step by step the class of flux reconstruction schemes introduced in [56]. To this aim, we place the general flux reconstruction framework presented in the above paragraph in perspective with the considerations of Section 2.3.2, and determine conditions on the layout of the discrete point sets $\kappa_Y, \kappa_F, \kappa_C$ and on the expression of the correction function for the resulting numerical scheme to be consistent and conservative.

**Choice of the discretisation points** Let us start by choosing sets of discretisation points $\kappa_Y, \kappa_F$, and $\kappa_C$ whose spatial distributions are respectively appropriate for representing the approximate solution in accordance with its multivariate polynomial nature, for allowing an accurate evaluation of the initial flux approximation, and for providing a control of the correction that is coherent with the overall procedure and enforces extra properties to the numerical scheme. Primarily focusing on the accuracy of the initially approximate flux evaluation, we can first obtain a desirable relationship between the flux and the solution discretisation points. Indeed, as any flux reconstruction technique entirely relies on pointwise values and polynomial reconstructions, the sole flux evaluation at one of its discretisation points is generally subject to an approximation error whenever the point location does not coincide with any of the solution points, the solution value being then approximated by its reconstruction process. It is therefore suitable to represent the discrete flux with a degree less or equal than for the solution, and set $\kappa_F \subset \kappa_Y$.

The actual distribution of the flux nodes, and therefore of part of the solution nodes, is *a-priori* free in each cell $K$. Typically, considering them on the boundary yields purely interfacial approaches as [60], while considering them within the cell makes the flux reconstruction scheme closer to classical
Discontinuous Galerkin approaches [61]. Adopting here the latter perspective, we consider sets $\kappa_F$ and $\kappa_Y$ where all the points lie strictly within the cell $K$, enhancing in particular the robustness of the solution within the cell and the flexibility of the corrected flux discretisation at the boundary.

The flexibility of the correction process is guaranteed along with its relevance by choosing a compactly supported correction function whose behaviour is in particular characterised from boundary constraints. Indeed, as evoked before, the cross-cell continuity of the corrected approximate flux can be easily obtained by coupling the behaviour of the two adjacent cells at the boundary and by propagating the information towards the nodes of the initial flux approximation. The definition of the correction function is then very flexible, its shape being in particular constrained by continuity and consistency relations at given points on the boundary. Thus, choosing to constrain the boundary behaviour of the correction function from a stencil $\kappa_C = \{x_j \in \partial K\}_j$ by comparing the shape of the initial flux approximation with the continuous one will naturally enforce the cross-cell continuity of the corrected flux, at least in the sense of the designed quantifiers. Furthermore, as then $\kappa_C \neq \kappa_F$, the correction function would be meaningful, not reducing to an ad-hoc zero interpolant. The stencil characterising the boundary behaviour of the correction function being selected, let us now precise the boundary characterisation type, points number, and points distribution while deriving a possible correction function shape. Possible inner characterisations will be described later in Section 2.5.

**Choice of the nature of the correction function** Let us first select a nature of the correction function that is compatible with the original discretisation performed on each cell $K$. We thus consider the discrete expression (2.73) and recall that both the initially approximate flux $F^{K,D}$ and the local approximate solution $Y^K$ are polynomial of degree $k$, and observe that it appears suitable to ask the divergence of the correction function $F^{K,C}$ to belong to $\mathbb{P}_k(K)$. It is however not necessary, the correction being propagated within the cell only to the specific flux points and the shape of the correction function not impacting the solution representation.

A-contrario, we do require the correction function to have a normal component belonging to $\mathbb{P}_k(f)$ on each of the cell faces $f \in \partial K$ for the correction procedure to act on the same discretisation space as the one in which the initially approximate flux is discretised in. It thus reduces to look for a correction function fulfilling

$$ (F^{K,C} \cdot n)|_{\partial K} \in \{w \in L^2(\partial K), \forall f \in \partial K, w|_f \in \mathbb{P}_k(f)\}, \qquad (2.74) $$

and optionally

$$ \nabla \cdot F^{K,C} \in \mathbb{P}_k(K). \qquad (2.75) $$
The normal component of the correction function thus being a polynomial on the boundary, the point based approach of the flux reconstruction technique can be preserved by designing a characterisation based exclusively on disconnected point values, not on moments employing the stencil $\varkappa_C$ in a quadrature scheme. The stencil $\varkappa_C$ is then used to define point-values degrees of freedom entirely characterising the correction function on the boundary. According to the property (2.74) of the correction function, we then distribute $k+1$ points discretisation points per face, yielding $N_C := |\varkappa_C| = n(k+1)$. The type of characterisation being now defined, let us select the constraints expressions.

**Obtention of the approximate flux continuity** The point-values constraints are then defined so as to erase the gaps between the initially approximate flux $F^{K,D}$ and the continuous flux at the associated boundary points by evaluating their difference, that is, by prescribing for any $x_j \in \varkappa_C$

$$F^{K,C}(x_j) = (\mathcal{F}(Y_j^h) \cdot n - F^{K,D}(Y_j^h)).$$

However, if the reconstruction of $F^{K,D}$ grants a simply defined flux value at any point $x_j \in \partial K$, this is not the case for the continuous flux evaluation, the solution at the boundary being itself possibly discontinuous. One is therefore led to define a numerical flux that is consistent with the continuous one and couples the normal flux information coming from either side of the boundary. In practice, the two different values obtained for $Y^h(x_j)$ for any $x_j \in \varkappa_C$ during its reconstruction, i.e. $Y_j^K := Y^K(x_j)$ for the value of $Y^h$ computed in the current element and $Y_j^{K-} := Y^{K-}(x_j)$ for the value computed using the information from the neighbouring element, are considered simultaneously to define a consistent numerical flux. Generically reading

$$\mathcal{F}^I: \quad \mathcal{D} \times \mathcal{D} \quad \to \mathbb{R}^q \\ (Y^K, Y^{K-}) \mapsto \mathcal{F}_n(Y^K, Y^{K-}),$$

it can typically be defined analogously as in Discontinuous Galerkin approach (see the equation (2.69)) by the expression

$$\mathcal{F}_n(Y^K, Y^{K-}) := \frac{1}{2} (\mathcal{F}(Y^K) \cdot n + \mathcal{F}(Y^{K-}) \cdot n) - \Delta(\mathcal{F}, n, Y^K, Y^{K-}),$$

where $\Delta$ represents a numerical dissipation operator (see e.g. [21] for more details). The discontinuities of the flux are then compensated at each point of $\varkappa_C$ by applying the approximated point value constraint

$$F^{K,C}(Y_j^K) = (\mathcal{F}^I(Y_j^K, Y_j^{K-}) \cdot n - F^{K,D}(Y_j^h)) \quad \forall x_j \in \varkappa_C, \qquad (2.76)$$

completely determining as a by-product the polynomial shape of the correction function on each face $f \in \partial K$.
**Propagation of the boundary correction within the cell** Once characterised from pointwise constraints, the correction amount is propagated towards the interior of the cell via space-dependent lifting functions that are defined everywhere on $K$, thus reaching out the points where the solution and initial flux approximation are discretised. Following here the approach of [56], the correction function is therefore expressed over the entire cell $K$ as

$$F^{K,C} = \sum_{j \in [1, N_C]} \left( \mathcal{F}^I(Y_j^K, Y_j^{K-}) \cdot n - F^{K,D}(Y_j^h) \cdot n \right) h_j(x) \quad (2.77)$$

where each term $h_j$ represents the lifting function that is associated with the amount corrected at the point $x_j \in \mathcal{X}_C$ and whose shape will be progressively derived. To this aim, we first notice that the term

$$\alpha_j := \left( \mathcal{F}^I(Y_j^K, Y_j^{K-}) \cdot n - F^{K,D}(Y_j^K) \cdot n \right)$$

is a constant, making the expression of the correction function entirely dependent on the shape of the lifting functions $h_j$. Hence, as the definition (2.77) is constructed as a linear combination of all the functions $h_j, j \in [1, N_C]$, each of them has in particular to fulfil the property (2.74). Their expression being further directly involved in the update process (2.73) under the divergence operator, we then define all the lifting function as being vectorial, i.e.

$$h_j: K \to \mathbb{R}^d \\ x \mapsto h_j(x), \quad \forall j \in [1, N_C]$$

The same information $\alpha_j$ is therefore propagated by (2.77) similarly in all spatial directions accordingly to the variations of the constrained function $h_j$, making the actual expression of the lifting functions $h_j, j \in [1, N_C]$ play an extensive role both in the validity of the correction process and in the way the correction amount echoes within the cell and on the boundary between the control points. In particular, while the vectorial expression of each lifting function is a-priori free, its boundary characterisation should preserve the continuity constraints and ensure the scheme conservativeness. Its inner characterisation can then be used later to enhance stability properties.

**Preservation of the cross-cell continuity and consistency properties**

Let us provide a first characterisation of any single lifting function $h_j, j \in [1, N_C]$ by guaranteeing the preservation of the cross-cell continuity obtained in (2.76). We there simply notice that given the definition (2.77), a simple way to ensure it is to prevent any interplay between the the individual lifting functions at the discrete correction points $\mathcal{X}_F$. We therefore require the normal component of the lifting function $h_j$ to vanish at every other point of $\mathcal{X}_C$ than $x_j$ itself. In other words, we ask each lifting function to fulfil

$$h_j(x_l) \cdot n|_{\partial K} = \delta_{jl} \quad \forall x_k \in \mathcal{X}_C, \forall l \in [1, N_C]. \quad (2.78)$$
One can here notice that as the normal component of each lifting function restricted to the boundary is by (2.74) a polynomial of order less or equal than $k$, and as the set $\varkappa_C$ distributes $k+1$ nodes on each face $f$ of the boundary $\partial K$, the normal component of any $h_j$ is already entirely determined by the condition (2.78) alone. Furthermore, the consistency of the reconstructed flux at the correction points is naturally ensured, the expression (2.77) guaranteeing by definition that $F^K(Y_j^K) \cdot n = \mathcal{F}^I(Y_j^K, Y_j^{K-}) \cdot n$.

**Obtaining local conservation** Let us now derive a further condition on the lifting functions so that the conservation of the scheme is ensured at the discrete level. To this aim, we investigate the conditions under which the relation (2.20) holds at the discrete level for any control volume $K \in \tau^h$. Reducing here to asking

$$ \int_K \nabla \cdot (F^{K,D}(Y^K) + \sum_{j \in [1, N_C]} \alpha_j h_j) d\gamma = 0, \quad (2.79) $$

we start by applying the Gauss divergence theorem to retrieve

$$ \int_{\partial K} (F^{K,D}(Y^K) \cdot n + \sum_{j \in [1, N_C]} \alpha_j h_j \cdot n) d\gamma = 0. \quad (2.80) $$

Making then the term $\alpha_j$ explicit, we can derive

$$ \int_{\partial K} \left( F^{K,D} \cdot n + \sum_{j \in [1, N_C]} \left( \mathcal{F}^I(Y_j^K, Y_j^{K-}) \cdot n - F^{K,D}(Y_j^K) \cdot n \right) h_j \cdot n \right) d\gamma = 0, \quad (2.81) $$

whose satisfaction for any cell $K \in \tau^h$ yields a natural assumption on the shape of the lifting functions. Indeed, for the discrete conservation principle (2.81) to hold without having to impose global constraints, it requires the corrected flux on a given face to be independent from its behaviour on any other one. We therefore ask each function $h_j$, $j \in [1, N_C]$ to have a boundary support restricted to the face on which its associated point $x_j \in \varkappa_C$ lies, i.e.

$$ \forall j \in [1, N_C], \quad \forall f \in \partial K, \quad \text{s.t. } x_j \notin f, \quad (h_j \cdot n)|_f = 0. \quad (2.82) $$

**Note** Such a condition means in watermark that each lifting function propagates the correction amount obtained from a given face towards the interior of the cell without impacting any other face. Compatible with the relation (2.78), it is hence a very natural assumption. ▲

The boundary behaviour of each lifting function being compactly supported on its associated face, the condition (2.81) then reduces to asking

$$ \int_f \left( F^{K,D} \cdot n + \sum_{\substack{j \in [1, N_C], \\ x_j \in \varkappa_C \cap f}} \left( \mathcal{F}^I(Y_j^K, Y_j^{K-}) \cdot n - F^{K,D}(Y_j^K) \cdot n \right) h_j \cdot n \right) d\gamma = 0. \quad (2.83) $$

for any face $f \in \partial K$ of any cell $K \in \tau^h$. The interplay between the weighted lifting functions makes the expression (2.83) too involved to derive a general
shape of the functions $h_j$ that would enforce the conservation property on top of the other conditions. However, as that we are here focusing only on the discrete level, the relation (2.83) has in practice to hold only at given quadrature points $\varkappa_G$. It is therefore enough to ask

$$ \sum_{g \in [1, |\varkappa_G|]} \omega_g \left( F^{K,D} \cdot n + \sum_{\substack{j \in [1, N_C], \\ x_j \in \varkappa_C \cap f}} \left( \mathcal{F}^I(Y_j^K, Y_j^{K-}) \cdot n - F^{K,D}(Y_j^K) \cdot n \right) h_j \cdot n \right) = 0 \quad (2.84) $$

for any quadrature rule made of the points $\varkappa_G$ and their associated weights $\{\omega_g\}_{g \in [1, |\varkappa_G|]}$, holding in particular when asking the relation

$$ F^{K,D}(Y_g^K) \cdot n + \sum_{\substack{j \in [1, N_C], \\ x_j \in \varkappa_C \cap f}} \left( \mathcal{F}^I(Y_j^K, Y_j^{K-}) - F^{K,D}(Y_j^K) \cdot n \right) h_j(x_g) \cdot n = \mathcal{F}^I(Y_g^K, Y_g^{K-}). $$

to be fulfilled for any $g \in [1, |\varkappa_G|]$. Generally strong, this assumption naturally holds when simply choosing the quadrature points as the correction nodes, i.e.

$$ \varkappa_G := \varkappa_C. \qquad (2.85) $$

Indeed, one can then apply the relations (2.84) to the nodes $\varkappa_C$, that reduces by virtue of (2.78) to

$$ F^{K,D}(Y_j^K) \cdot n + \left( \mathcal{F}^I(Y_j^K, Y_j^{K-}) - F^{K,D}(Y_j^K) \cdot n \right) = \mathcal{F}^I(Y_j^K, Y_j^{K-}), $$

holding by for any $j \in [1, N_C]$. The discrete conservation is therefore achieved in the sense of all quadrature rules having for stencil (2.85).

### Summary of the requirements on the correction function design
So far, we could set a general structure of the correction function $F^{K,C}$ and derive conditions over the lifting functions $h_j$, $j \in [1, N_C]$ so that the resulting numerical is consistent and conservative. However, no proper shape for the lifting function could yet be derived. To this aim, let us summarise the admissibility conditions applying on the lifting functions $h_j$, $j \in [1, N_C]$.

#### Admissibility conditions 2.4.1 Admissibility criterion on lifting functions

Let $h_j$, $j \in [1, N_C]$ be a lifting function associated to the discretisation point $x_j \in \varkappa_C$. Its definition can contribute to the expression (2.77) of the correction function $F^{K,C}$ whenever

$$ (h_j^{K,C} \cdot n)|_{\partial K} \in \{w \in L^2(\partial K), \forall f \in \partial K, w|_f \in \mathbb{P}_k(f)\}, \quad (2.74) $$

$$ h_j(x_l) \cdot n|_{\partial K} = \delta_{jl}, \quad \forall j \in [1, N_C], \forall x_l \in \varkappa_C, \quad (2.78) $$

$$ (h_j \cdot n)|_f = 0, \quad \forall j \in [1, N_C], \forall f \in \partial K, \text{s.t. } x_j \notin f. \quad (2.82) $$

We can here observe that no constraint has been yet formulated within the cell, allowing us to define freely the lifting functions as members of a finite dimensional variational space in which the properties (2.74), (2.78), and (2.82) are either obtained naturally or can be enforced.
**A framework for admissible correction functions** When the cells $K$ all have a simplicial shape, each lifting function can be taken in the $(N_C-1)^{th}$ order Raviart – Thomas space, that is as $RT_k(K) := (\mathbb{P}_k(K))^d + x \mathbb{P}_k(K)$ for $k := N_C - 1$ (see Section 3.2.1 for a more detailed description).

Indeed, there, the condition (2.74) is obtained by construction, and the nature of the Raviart – Thomas space allows to consider this sole polynomial characterisation as being part of the global characterisation of any function living in the $(N_C-1)^{th}$ order Raviart – Thomas space. Thus, given a basis function of Raviart – Thomas, one can initially define

<figure><img src="image_26.png" /><figcaption>**Figure 2.21.** Choice of a set of lifting functions.</figcaption></figure>

the set of functions $h_j$ as the subset of basis functions that characterise the boundary behaviour of the normal components, requiring it to form the nodal basis of $RT_k(K)|_{\partial K}$ that is associated to the discretisation points $\varkappa_C$. Namely, for all $j \in [1, N_C]$ the lifting functions thus read

$$h_j: K \to \mathbb{R}^d \\ x \mapsto (x + n_f) l_{f,j}(x), \quad f \in \partial K \text{ s.t. } x_j \in f, \qquad (2.86)$$

where $l_{f,j}$ denotes here the Lagrangian function that corresponds to the node $x_j$ among the set that is associated to the points $\{x_j \in \varkappa_C, x_j \in f\}$, and for which the detail of the construction can be found in e.g. [62] or [63]. Indeed providing a basis of $RT_k(K)|_{\partial K}$, this choice however does not correspond to any subset of the classical Raviart – Thomas elements basis, being entirely modal. If the discretisation space is chosen as being the Raviart – Thomas one, the considered element are thus differ from the original setting.

This choice yields to define correction functions (2.77) that also lives in the $(N_C-1)^{th}$ order Raviart – Thomas space by linearity of their construction with respect to the lifting functions. Hence, the finite dimensional space that is associated to the discrete correction functions reads

$$V^C(K) := (\mathbb{P}_k(K))^d + x \mathbb{P}_{[k]}(K),$$

from which the comprehensive flux discretisation framework can be inferred. Considering conjointly the two variational spaces $V^D(K)$ and $V^C(K)$ in which parts of the corrected approximate flux function can be expressed, one can observe that the comprehensive flux discretisation takes place in the finite dimensional polynomial space

$$V(K) = V^D(K) + V^C(K) = (\mathbb{P}_k(K))^d \oplus x \mathbb{P}_{[k]}(K),$$

where $\zeta = \{\zeta_i, i \in [1, d]\}$ is the set of all permutations of $(1, \dots, d)$, see Section 3.2.1.1 of Chapter 3 for more details, and that the overall characterisation is unisolvent in $V(K)$.
**Proposition 2.4.2** *Unisolvence of the flux characterisation in $V(K)$*

Given the initial flux approximation (2.72) and the correction functions (2.77) constructing from the lifting functions (2.86), the overall characterisation $\Sigma_{V(K)}$ is unisolvent in $V(K)$.

**Proof.** As in this specific case we can derive the relation $V^D(K) \subset V^C(K) = V(K)$, the above statement follows immediately from Proposition 2.4.1. Indeed, the space $V(K)$ being here a polynomial, any subset $\Sigma_{C^\cap}$ of $\Sigma_C$ that contains $\frac{(k)(k+1)}{2}$ degrees of freedom would characterise any polynomial of degree $k$ defined on $K$. As the dimension of $\Sigma_C$ is by construction bigger or equal to $\frac{(k)(k+1)}{2}$, the hypothesis of Proposition 2.4.1 apply and the unisolvence of the comprehensive flux discretisation associated to the Flux Reconstruction scheme (2.88) follows. ■

Once the boundary setting is defined, it is then possible whenever $N_C > 1$ to refine the definition of each lifting function. Indeed, in this case, the dimension of the Raviart–Thomas space is greater than the dimension of its restriction on the boundary. The description of the Raviart–Thomas space is then completed by other basis function, whose normal component vanish on all boundaries, the polynomial representation of the any function’s normal component being already complete by using the set of nodal basis functions of $RT_{N_C-1}(K)|\partial K$. One can therefore modify each of the predefined correction functions $h_j, j \in [1, N_C]$ by adding a linear combination of those complementary basis functions without destroying the properties required for fulfilling the Admissibility conditions 2.4.1. The characterisation of this linear combination can therefore be done upon constraints ensuring stability and physical relevance of the solution, as we will see in Section 2.5.

<figure><img src="image_27.png" /><figcaption>**Figure 2.22.** Necessity of a flexible setting that guarantees the lifting functions a same discretisation framework for any polygonal shape.</figcaption></figure>

However, up to the author knowledge the definition of the Raviart–Thomas space is limited to simplicial cells, and can be straightforwardly extended only to quadrilateral cells. A new setting is therefore required for allowing the use of general polygonal shapes, and will be the object of the entire Chapter 3.
### 2.4.2.3 Solution update given the reconstructed flux

Given the definition (2.77) of the correction function built on functions fulfilling the Admissibility conditions 2.4.1, the flux reconstruction scheme as presented in [56] is now conservative and consistent. We can therefore derive the proper expression of the update procedure initiated in (2.73), by providing a precise definition of the operator $D^h$. It comes

$$ \frac{dY_{\mathcal{G}}^K}{dt} = -(\nabla \cdot F^{K,D})(x_m) - (\nabla \cdot F^{K,C})(x_m) \quad (2.87) $$

$$ = -\nabla \cdot (F^{K,D})(x_m) - \sum_{j \in [1, N_C]} (\alpha_j \nabla \cdot h_j(x_m)) \quad (2.88) $$

that in practice reduces to (2.73)

$$ \frac{dY_{\mathcal{G}}^K}{dt} = -\left( \sum_{\sigma \in \Sigma_{F_1}} F_{1, \sigma}^K \frac{\partial l_{\sigma}}{\partial x_1}(x_m) + \sum_{\sigma \in \Sigma_{F_2}} F_{2, \sigma}^K \frac{\partial l_{\sigma}}{\partial x_2}(x_m) \right) - (\nabla \cdot F^{K,C})(x_m) \\ = -\left( \sum_{\sigma \in \Sigma_{F_1}} F_{1, \sigma}^K \frac{\partial l_{\sigma}}{\partial x_1}(x_m) + \sum_{\sigma \in \Sigma_{F_2}} F_{2, \sigma}^K \frac{\partial l_{\sigma}}{\partial x_2}(x_m) \right) - \sum_{j \in [1, N_C]} \alpha_j \nabla \cdot h_j(x_m), $$

ending the definition of the numerical scheme and allowing us to design further constraints on the lifting functions to enforce stability properties.

<figure><img src="image_28.png" /><figcaption>Figure 2.23. Update procedure of the Flux Reconstruction scheme built upon the chosen correction functions and discretisation points.</figcaption></figure>

Let us now derive a possible definition of constraints to apply on the lifting functions for the resulting Flux Reconstruction scheme to be entropy conservative and to provide a physically relevant solution. To this aim, we follow the construction of [52] and rely on a parallel between the Flux Reconstruction and Residual Distribution frameworks, deriving an entropy conservative scheme in the Residual Distribution framework by introducing a correction term. This term is finally exploited to provide conditions over the lifting functions for the Flux Reconstruction scheme to be entropy conservative.
### 2.4.3 Residual distribution framework

Very general, the Residual Distribution approach tackles the problem (2.1) by constructing a discrete dual formulation. Namely, instead of expressing the solution from a finite dimensional set of basis functions and update it by approximating the dynamics locally on control volumes that match with the mesh cells, Residual Distribution schemes exploit the definition of the degrees of freedom and act on the mesh variational dual.

**General concept** The definition of a Residual Distribution scheme classically starts by selecting a finite dimensional variational space $V^h$ in which the operator $\mathcal{D}_F^h$ will be discretised, a trial space $W^h$ where to express the discretised solution, and a set of degrees of freedom $\Sigma^h$ having for dimension $|W^h|$ and whose values characterise any function of the trial space.

The structure of those spaces are chosen piecewise, that is, there exists a conformal mesh $\tau^h$ and finite dimensional variational spaces $V$ and $W$ so that $V^h := \{V(K), K \in \tau^h\}$, $W^h := \{W(K), K \in \tau^h\}$. Hence, one can define local sets of degrees of freedom $\Sigma_{V(K)}$ describing the local dynamic in the element $(K, V(K), \Sigma_{V(K)})$, and obtain a local description of the solution for each cell $K$ by extracting the set of local degrees of freedom $\Sigma_{W(K)}$ out of $\Sigma^h$ that characterise locally the solution in the element $(K, W(K), \Sigma_{W(K)})$. Furthermore, the treatment of the boundary conditions may require specific treatments impacting the value of the solution degrees of freedom that are entirely supported on the boundary. We therefore denote for the sake of later use $\Sigma_{W(\Gamma)} := \{\mathcal{G} \in \Sigma_{W(K)}|_{\Gamma}, K \text{ s.t. } \Gamma \in K, K \in \tau^h\}$ the set of degrees of freedom whose support are contained on the boundary face $\Gamma \in \partial\tau^h$.

**Note** Whenever the setting $W^h$ provides a continuous solution over $\Omega$, $\Sigma_{W(K)} \cap \Sigma_{W(K^{-})} \neq \emptyset$. However, any discontinuous setting would never share some degrees of freedom between elements, as even though their expressions may be identical, the restriction to the boundary of the space they are affiliated with would be different, and so might be their associated values. ▲

The description of the solution and dynamics so defined, the total impact of the local fluxes on the approximate dynamics is then evaluated for each primal element, contributing to the values of the degrees of freedom characterising the functions living in $W(K)$ accordingly to a given weighting (see an illustration in Figure 2.24). Both the choice of the weighting proportions and of the initial dynamics approximation therefore define the numerical scheme.

**Framework definition** Let us now derive the general expression of a Residual Distribution scheme by considering the formulation (2.63) and a conformal mesh $\tau^h$ that approximates the domain $\Omega$. There, the first step of the
approximation procedure is to consider the impact of the relation (2.63) over the entire computational domain, that is, to ask the total contribution of the dynamic to the state evolution to preserve all the conservation principles acting on $\Omega$, namely

$$ \int_{\tau^h} \nabla \cdot \mathcal{F}(Y^h) \, dx - \int_{\partial \tau^h} g \, dx = 0. \qquad (2.89) $$

Using then the piecewise structure of the approximate solution $Y^h$ and fore-seeing a local piecewise dynamics approximation, we can split the relation over each cell and write

$$ \sum_{K \in \tau^h} \int_K \nabla \cdot \mathcal{F}(Y^h) \, dx - \sum_{\Gamma \in \partial \tau^h} \int_{\Gamma} \mathcal{F}(Y^h) \cdot n - F_b(Y^h, g) \, dx = 0, \qquad (2.90) $$

where we reformulated the boundary conditions in the same manner as when deriving the relation (2.65). Thereon, given that each local solution has a compact support on $W(K)$, we can consider the total contribution of the local fluxes to the global dynamic by defining the *total residuals*, the set of fluctuation terms

$$ \Phi^K(Y^h) = \int_K \nabla \cdot \mathcal{F}(Y^h) \, dx \qquad (2.91a) $$

$$ \Phi^{\Gamma}(Y^h) = - \int_{\Gamma} \mathcal{F}(Y^h) \cdot n + F_b(Y^h, g) \, dx, \qquad (2.91b) $$

associated to the element $(K, W(K), \Sigma_{W(K)})$, $K \in \tau^h$. The relation (2.90) then reduces to

$$ \sum_{K \in \tau^h} \Phi^K(Y^h) + \sum_{\Gamma \in \partial \tau^h} \Phi^{\Gamma}(Y^h) = 0, \qquad (2.92) $$

expression from which a dual formulation of the problem (2.89) can be derived by having fractions of these residuals contribute to the values of the degrees of freedom supported on the cell. To this aim, we define for any cell $K \in \tau^h$ a set of coefficients $\beta_{\mathcal{S}}^K, \mathcal{S} \in \Sigma_{W(K)}$ and $\beta_{\mathcal{S}}^{\Gamma}, \mathcal{S} \in \Sigma_{W(K)|_{\Gamma}, K \text{ s.t. } \Gamma \in K}$ so that

$$ \sum_{\mathcal{S} \in \Sigma_W(K)} \beta_{\mathcal{S}}^K = 1 \quad \text{and} \quad \sum_{\mathcal{S} \in \Sigma_W(\Gamma)} \beta_{\mathcal{S}}^{\Gamma} = 1 \qquad (2.93) $$

and construct the associated partial residuals

$$ \Phi_{\mathcal{S}}^{K}(Y^{h}) = \beta_{\mathcal{S}}^{K} \Phi^{K}(Y^{h}) \quad \forall \mathcal{S} \in \Sigma_{W(K)} \qquad (2.94a) $$

$$ \Phi_{\mathcal{S}}^{\Gamma}(Y^h) = \beta_{\mathcal{S}}^{\Gamma} \Phi^{\Gamma}(Y^h) \quad \forall \mathcal{S} \in \Sigma_{W(\Gamma)}, \qquad (2.94b) $$

Since then by the (2.93) we obtain the conservation property

$$ \sum_{\mathcal{S} \in \Sigma_{W(K)}} \beta_{\mathcal{S}}^{K} \Phi^{K}(Y^{h}) = \Phi^{K}(Y^{h}) \quad \text{and} \quad \sum_{\mathcal{S} \in \Sigma_{W(\Gamma)}} \beta_{\mathcal{S}}^{\Gamma} \Phi^{\Gamma}(Y^{h}) = \Phi^{\Gamma}(Y^{h}), \qquad (2.95) $$

the relation (2.92) becomes

$$ \sum_{K \in \tau^h} \left( \sum_{\mathcal{S} \in \Sigma_{W(K)}} \Phi_{\mathcal{S}}^{K}(Y^h) \right) + \sum_{\Gamma \in \partial \tau^h} \left( \sum_{\mathcal{S} \in \Sigma_{W(\Gamma)}} \Phi_{\mathcal{S}}^{\Gamma}(Y^h) \right) = 0. \qquad (2.96) $$