QUEEN M Jens Mueller

Président O F London

BCE KARMIEL, Invité(e) M. CORENTIN M Amir Adler

Lapeyre

Membre M Cerfacs

Ignacio José

Onera Cardesa

Membre Toulouse

Andrea Mme

Universitat Beck

Membre M Stuttgart

Toulouse Thierry Poinsot

Membre M Inp

Dantas Luciano Drozda

Martins

M Thierry

Poinsot M Corentin

Lapeyre M Tapan

Kumar Sengupta

Iit Kanpur

Keywords: Computational Fluid Dynamics, Local Transfer function Analysis, Data-driven numerical schemes, Machine Learning

The numerical simulation of fluid flows has become an essential part of the virtual prototyping in the aerospace industry. It is made possible through Computational Fluid Dynamics (CFD) tools. These tools rely on algebraic relations between the values stored at mesh points called numerical schemes, which must fulfill specific criteria to lead to reliable simulations. Namely, numerical schemes must be accurate without sacrificing stability (among other conditions). Traditionally, accuracy and stability conditions have been derived for problems governed by linear equations on regular meshes. However, fluid motion equations are non-linear, and meshes used in industrial CFD can be irregular (especially when dealing with complex geometries). This thesis proposes a novel framework to analyze the accuracy and stability of numerical schemes for non-linear systems of equations and irregular meshes. This framework also establishes conditions for optimizing numerical schemes locally in time and space. It is named Local Transfer function Analysis (LTA). Under the light of LTA, each mesh element acts as an impedance block that resists to the solution propagation over time and space. The optimization of numerical schemes becomes an impedance-matching problem. It is solved by minimizing the value of an objective function. The LTA objective function measures the distance between the dynamics predicted by the numerical scheme and some reference dynamics. Using reference data in this optimization process turns LTA into a tool that generates accurate and stable data-driven numerical schemes. More recently, Machine Learning (ML) has emerged as a field dedicated to extracting knowledge from data. This thesis then proposes employing ML architectures to find optimal values for the parameters of a numerical scheme in the sense of the LTA objective function. The method is applied first on 1D problems modeled by the Convection/Burgers' equations and finally extended to 2D applications governed by the Convection/Euler equations.

grad x F Gradient of a scalar-valued function F with respect to x

Symbols x

Tensor of order 1 (i.e., a vector)

x Tensor of order greater than unity (1962). From [START_REF] Rosenblatt | The perceptron: A probabilistic model for information storage and organization in the brain[END_REF]. = exp(3x) at x = 0.2 using the approximation (1.12). . . 1.11 Computational graph and forward primal evaluation trace for a program evaluating f (x 1 , x 2 ) := sin(x 1 + x 2 ), where x 1 , x 2 ∈ R. 1.12 (Top) Schematic representation of the sliding ladder problem.

List of Figures

(Middle and bottom) Computational graph, forward primal, and forward tangent evaluation traces for a program evaluating the angle y 1 (θ sep in Eq. (1.17)) at which the falling ladder loses contact with the vertical wall and the speed y 2 (s f in Eq. (1.18)) at which it hits the floor. . . . . . . . . . . . . . .

2.8

Discretization coefficients along the x-direction generated by a CNN-PAL model for the vortex convection analyzed in experiment n.1 (depicted in Tbl. 2.2). The "UNTRAINED" columns refer to the coefficients returned by the CNN-PAL model before training, while the "TRAINED" columns relate to the coefficients returned by the same CNN-PAL model after training. The coefficients are assigned to the neighborhood {(i-2, j), (i-1, j), (i, j), (i+1, j), (i+2, j)} around each mesh point (i, j), composing a centered stencil of size N = 5 (see Tbl.

Introduction

In the aerospace industry, engineering systems design must consider the effects of fluid motion around airfoils or inside nozzles, among other components. Nevertheless, real-world prototyping of these systems can be excessively expensive. Their numerical simulation through Computational Fluid Dynamics (CFD) has become essential to explore design spaces and speed up development [START_REF] Hirsch | Numerical Computation of Internal and External Flows[END_REF]. First frameworks such as Boeing's PAN AIR [START_REF] Magnus | PAN AIR: A computer program for predicting subsonic or supersonic linear potential flows about arbitrary configurations using a higher order panel method[END_REF] employed simple potential flow models, where the velocity field is assumed to be irrotational. Over the years, the increasing availability of computational power has allowed manufacturers to employ more detailed mathematical models to improve predictions. The most accurate model that describes the motion of viscous fluid substances is the set of Navier-Stokes (NS) equations. In the motion of viscous fluids, velocity and pressure fluctuations of chaotic nature may arise whenever inertial forces are strong enough compared to viscous forces. These fluctuations constitute a phenomenon called turbulence. [START_REF] Kolmogorov | The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds' Numbers[END_REF] showed that turbulent flows are characterized by an energy transfer from their largest to their smallest eddies. Spatial discretization of the NS equations down to the length scale of the flow corresponding to the smallest eddies remains prohibitively expensive in industrial configurations [START_REF] Sagaut | Large eddy simulation for incompressible flows: an introduction. 3rd ed. Scientific computation[END_REF]. State-ofthe-art simulation frameworks still consider modified forms of the underlying Partial Differential Equations (PDEs). In combustion and aero-acoustics, for instance, a filtered version of the NS equations is solved instead, such as in Large-Eddy Simulations (LES) [START_REF] Poinsot | Theoretical and numerical combustion[END_REF]Wagner et al., 2007). In LES, the effects of eddies of size under a given threshold are modeled using subgrid closure rules. The first LES is credited to [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers[END_REF] for the case of a simple plane Poiseuille flow. More recently, Pérez [START_REF] Arroyo | Towards the Large-Eddy Simulation of a full engine: Integration of a 360 azimuthal degrees fan, compressor and combustion chamber. Part I: Methodology and initialisation[END_REF] performed the first high-fidelity LES of a full aircraft engine, thanks to the evolution of computational resources.

INTRODUCTION

The increase in available computational power has also contributed to the emergence of large-scale applications in an apparently uncorrelated field: Machine Learning (ML). It aims at extracting knowledge from data to improve task automation [START_REF] Jordan | Machine learning: Trends, perspectives, and prospects[END_REF]. An ML architecture is a parametrized function to be adjusted in order to execute a given task better. The process of using data to tune the parameters of an ML model is called training. The widespread availability of Graphics Processing Units (GPUs) enabled the training of ML models at a much faster pace than with previous devices [START_REF] Chellapilla | High Performance Convolutional Neural Networks for Document Processing[END_REF][START_REF] Steinkraus | Using GPUs for machine learning algorithms[END_REF]. As a consequence, during the past two decades, ML applications in a myriad of technology areas have been developed. Some examples are the automation of speech recognition [START_REF] Nassif | Speech Recognition Using Deep Neural Networks: A Systematic Review[END_REF], image segmentation [START_REF] Minaee | Image Segmentation Using Deep Learning: A Survey[END_REF], and autonomous driving [START_REF] Kiran | Deep Reinforcement Learning for Autonomous Driving: A Survey[END_REF]. Recently, the ML-based chatbot GPT-4 exhibited human-level performance on various professional and academic benchmarks (OpenAI, 2023). The ability of an ML model to solve problems across different contexts or disciplines is called Artificial General Intelligence (AGI) [START_REF] Goertzel | Artificial General Intelligence: Concept, State of the Art, and Future Prospects[END_REF]. The capacities of GPT-4 in handling tasks that span mathematics, vision, coding, law, psychology, medicine, and more have sparked a debate on whether it could be considered an early (yet incomplete) version of an AGI system [START_REF] Bubeck | Sparks of Artificial General Intelligence: Early experiments with GPT-4[END_REF]. It opens promising perspectives for using ML models in all aspects of human life.

Most of the ML applications mentioned above employ (artificial) neural networks: a class of ML architectures inspired by information processing in biological systems [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. They are universal approximators of continuous functions over compact domains [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF]. This fact triggered the scientific community's interest in possible ML applications in fields such as CFD, from turbulence modeling [START_REF] Duraisamy | Turbulence Modeling in the Age of Data[END_REF] to the development of ML-based numerical solvers [START_REF] Karniadakis | Physics-informed machine learning[END_REF]. The former is concerned with using neural networks as surrogate models of the unresolved scales of fluid motion. The latter aims to employ neural networks to resolve the governing equations. Its goal is to improve the accuracy of the numerical predictions at cost levels similar to those of a traditional solver. This thesis will only address the development of ML-based numerical solvers. In this context, [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF] proposed to approximate the fluid flow function of space and time by directly minimizing the solution residual on a sample of observation points, a popular approach known as Physics-Informed Neural Networks (PINNs). This is a powerful method for problems where multiple constraints compete or when parts of the problem are not well known INTRODUCTION (e.g., boundary conditions in hydrology, initial conditions in weather predictions. . . ). However, for well-defined problems, it is less accurate compared to direct resolution [START_REF] Fuks | Limitations of physics informed machine learning for nonlinear two-phase transport in porous media[END_REF]. In the scope of fluid mechanics, P. [START_REF] Sharma | A Review of Physics-Informed Machine Learning in Fluid Mechanics[END_REF] noted that the error levels in the prediction of turbulent flows using PINNs are often orders of magnitude higher than those of laminar flows. This is due to problems associated with non-linearity, highdimensionality, and multi-scale physics. S. [START_REF] Wang | When and why PINNs fail to train: A neural tangent kernel perspective[END_REF] showed that PINNs suffer from spectral bias, a well-known pathology that prevents deep fully-connected neural networks from learning high-frequency functions [START_REF] Rahaman | On the Spectral Bias of Neural Networks[END_REF]. [START_REF] Markidis | The Old and the New: Can Physics-Informed Deep-Learning Replace Traditional Linear Solvers?[END_REF] assessed the potential of PINNs as linear solvers in the case of Poisson's equation (found in incompressible flow solvers, among many other applications in scientific computing). Markidis (2021) also found that PINNs learn the low-frequency components of the solution field more easily than the high-frequency ones, showing that accuracy is a limiting factor for using PINN linear solvers alone, and eventually recommended applying hybrid strategies where PINNs are integrated into traditional linear solvers. By doing this, one could develop ML-based solvers that are competitive with state-of-the-art ones, both performance-and accuracy-wise.

Combining traditional numerical methods with ML to propose a new class of ML-based numerical schemes for CFD is the main focus of this thesis. This project does not aim to rely on an ML model alone to unveil some system dynamics but to use the traditional methods as inductive biases. More precisely, ML models will be used to tune the parameters of numerical schemes. Tuning discretization parameters has been essential to practical numerical algorithms to solve PDEs. They aim to ensure that the numerical solution satisfies certain desired but critical characteristics of the physical problem. These characteristics relate mainly to stability, dispersion, and dissipation numerical errors. [START_REF] Rajpoot | Optimal time advancing dispersion relation preserving schemes[END_REF] and T. K. [START_REF] Sengupta | Space-time discretizing optimal DRP schemes for flow and wave propagation problems[END_REF] proposed optimized compact finite-difference schemes that accurately capture the dispersion, dissipation, and group velocity characteristics of the physical problem through a comprehensive approach considering both temporal and spatial discretizations along with the governing equation. [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF] performed parameter optimization of compact finite-difference schemes to provide improved representation of a range of spatial scales. However, they considered exact time advancement, which leads to an incorrect evaluation of the dispersive error characteristics. [START_REF] Tam | Dispersion-Relation-Preserving Finite Difference Schemes for Computational Acoustics[END_REF] introduced finite-difference schemes whose parameters are optimized to reduce truncation error for first order spatial derivative. Again, INTRODUCTION space and time discretizations were treated independently, which results in an incorrect numerical dispersion relation. Besides, they used a four time level method which introduces two spurious numerical modes apart from the physical mode (as explained in T. K. Sengupta and Dipankar (2004)). In the realm of ML-based schemes, [START_REF] Bar-Sinai | Learning data-driven discretizations for partial differential equations[END_REF] and [START_REF] Kochkov | Machine learning-accelerated computational fluid dynamics[END_REF] employed ML models to predict coefficients of a finite-difference method, leading to accurate numerical solutions on relatively coarse grids. However, they did not establish conditions for the stability of the resulting scheme. Using neural networks, [START_REF] Stevens | Enhancement of shock-capturing methods via machine learning[END_REF] improved a fifth-order shock-capturing Weighted Essentially Non-Oscillatory (WENO) scheme [START_REF] Jiang | Efficient Implementation of Weighted ENO Schemes[END_REF]. Their ML model predicted optimal perturbations to WENO coefficients. Nevertheless, the only constraint applied to the perturbed coefficients was an affine transformation which guarantees that they sum up to one. Hence, the resulting numerical method was only first-order accurate. [START_REF] Kossaczká | Enhanced fifth order WENO shock-capturing schemes with deep learning[END_REF] proposed instead to perturb the smoothness indicators of the WENO method. Their resulting scheme presented less diffusion and less overshoot in shocks than the original WENO scheme while maintaining high order of accuracy in the smooth regions of the flow. [START_REF] Bezgin | A data-driven physics-informed finite-volume scheme for nonclassical undercompressive shocks[END_REF] applied neural networks to adaptively reconstruct nonlinear fluxes in the context of nonclassical shock waves, i.e., shock waves which violate the Lax entropy condition [START_REF] Lefloch | Hyperbolic Systems of Conservation Laws[END_REF].

This thesis concerns developing ML-based numerical schemes on irregular and unstructured grids typical of industrial CFD. First, a novel spectral analysis that can be applied to non-linear systems of equations and irregular meshes is proposed. It is worth noting that the proposed analysis does not consider aliasing errors, which results from evaluation of product terms in a finite resolution grid (T. [START_REF] Sengupta | High Accuracy Computing Methods: Fluid Flows and Wave Phenomena[END_REF]. Such a spectral analysis establishes conditions for optimizing numerical schemes locally in time and space. An ML model is then used to predict optimal values of a parameter of a traditional numerical method used in LES. Finally, applications to problems governed by the Convection/Burgers' equations on 1D irregular meshes and the Convection/Euler equations on 2D irregular triangular meshes demonstrate how the ML-based numerical scheme can outperform a traditional one accuracy-wise.

This manuscript is structured as follows:

Chapter 1 -This chapter introduces the reader to the Machine Learning field. It will be seen that commonly used neural network models like the Multilayer

INTRODUCTION

Perceptron (MLP), the Convolutional Neural Network (CNN), and the Graph Neural Network (GNN) arose from problems outside the scope of scientific discovery but have found many applications therein due to their specific features.

Chapter 2 -This chapter briefly reviews the development of numerical schemes for hyperbolic conservation laws, typical to CFD. Subsequently, the datadriven numerical scheme of [START_REF] Bar-Sinai | Learning data-driven discretizations for partial differential equations[END_REF] and [START_REF] Kochkov | Machine learning-accelerated computational fluid dynamics[END_REF] is presented. Their method employed a CNN model to explore the space of possible finite-difference spatial discretization coefficients obtained via Taylor series expansions. It is applied to the system of compressible Euler equations in 2D cartesian grids. The method does not guarantee stability despite improving accuracy levels on a particular set of problems.

Chapter 3 -A theoretical framework that establishes accuracy and stability conditions for coupled ML-PDE systems is developed in this chapter. A novel spectral analysis based on local spatio-temporal features named the Local Transfer-function Analysis (LTA) is proposed. One-dimensional problems governed by the Convection and Burgers' equations are studied under LTA. The considered numerical scheme belongs to the family of Taylor-Galerkin schemes used in LES.

Chapter 4 -The Two-step Taylor-Galerkin C (TTGC) numerical scheme of [START_REF] Rudgyard | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF] is modified to incorporate a cell-valued parameter provided by a GNN model, which is constrained to respect stability conditions predicted by LTA. The resulting numerical method is termed ML-TTGC and can address some limitations of the original TTGC scheme by successfully damping spurious oscillations on irregular meshes and preserving amplitudes of high wavenumber waves in the linear convection problem.

Chapter 5 -Finally, an extension of the method to higher dimensions is sought in Chapter 5. In particular, an implementation of the TTGC scheme for 2D Convection and Euler problems on unstructured meshes is coupled to a GNN model identical to the one used in the 1D experiments of Chapter 4. The resulting framework confirms ML-TTGC's capacity to damp wiggles on irregular meshes. Moreover, in the case of a double

INTRODUCTION

shear layer problem governed by the Euler equations, ML-TTGC proves to be more stable than the original scheme.

Fig. 1 summarizes the manuscript's structure and highlights the relationships between the chapters.

Figure 1: Manuscript's structure. Arrows denote relationships between chapters. More precisely, an arrow from chapter A to chapter B shows that concepts defined in chapter A will be employed in chapter B.

Chapter 1

Introducing machine learning

This chapter surveys the field of Machine Learning (ML). Sec. 1.1 presents standard neural network architectures like multilayer perceptrons and convolutional neural networks. They have been applied by researchers in many scientific applications, despite originating in the image recognition realm. This section also presents graph neural networks. Unstructured data (typical to industrial CFD applications) are naturally expressed in graphs, such that graph neural networks are employed in the data-driven numerical schemes developed in Chapters 4 and 5. These neural network models are trained using gradient-based methods. Sec. 1.2 introduces the reader to algorithmic differentiation, an efficient way to compute gradients of the objective function for training with respect to the ML model parameters. This chapter ends in Sec. 1.3 by presenting some algorithmic differentiation frameworks.

Machine learning: from checkers to scientific computing

IBM engineer Arthur Samuel is credited for popularizing the term machine learning in his attempts to program computers to learn from data during the 1950s. He believed that learning devices would suppress the need for detailed programming effort in problem-solving, as machine-learned rules replace hard-coded algorithmic steps. The realm of games seemed to him a convenient way to develop the first learning techniques due to their relative Figure 1.1: Arthur Samuel and IBM 700 (1956). From [START_REF] Gladchuk | A history of Machine Learning from 1940s to present days[END_REF].

boundedness in comparison to life problems in general. In 1952, Samuel then created a checkers-playing program that would not build all possible paths until the game's conclusion, but rather evaluates a scoring function estimating winning chances at any given game state [START_REF] Samuel | Some studies in machine learning using the game of checkers[END_REF]. Such a scoring function also depended on a history of game actions, from which the machine is expected to continuously improve while playing. Samuel's methods relate to modern techniques of the so-called reinforcement learning field (Richard S. [START_REF] Sutton | Reinforcement Learning: An Introduction, 2nd Edition[END_REF]. Remarkably, his program was able to play better than himself after some 10 hours of device-playing time. Fig. 1.1 shows him with an IBM 700 and a checkers board.

Later on, the American psychologist Frank [START_REF] Rosenblatt | The perceptron: A probabilistic model for information storage and organization in the brain[END_REF] invented the perceptron, a learning unit inspired from human brain neurons that would trigger the modern deep learning field. Rosenblatt's perceptron was designed to become a machine for pattern recognition (capable of distinguishing geometrical shapes, letters of the alphabet, among others), under the ideas stated in the following words : " [. . . ] our objective has been to discover a physical system, or abstract model, which will be capable of "perceiving" its environment, and learning to recognize those objects or events which it has perceived in the past. However, since it is our purpose to understand the actual mechanisms employed by the brain, rather than simply to construct a new type of computing device, the perceptron models are constrained in their organization and dynamic properties by what is known of the biological nervous system. Rather than attempting to "invent" or "construct" a machine which will calculate such things as similarities or geometrical properties of stimuli, the approach has been to begin with a hypothetical network of idealized neurons, or nerve cells, resembling the brain in its general organization, and then analyze the system mathematically to determine whether or not it possesses "psychological" properties of interest. " (From Rosenblatt (1962), p. 574) Its original model consisted of three layers S, A, R connected in series to compose a perception mechanism [START_REF] Van Der Malsburg | Frank Rosenblatt: Principles of Neurodynamics: Perceptrons and the Theory of Brain Mechanisms[END_REF], as illustrated in Fig. 1.2. The layer S corresponds to a sensorial surface (e.g., eye's retina) that reacts to environment stimuli thanks to its set of sensory units. Their response then fed layer A as a set of association cells responsible for detecting properties that are particular to the object to be recognized. Finally, layer R of recognition units receives the outputs from A and should fire only for a specific type of pattern projected on S (e.g., the letter C shown in Fig. 1.3(a)).

In mathematical terms, the output signal y from a given R-unit can be expressed as [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]

y := h i w i ϕ(x i ) ≡ h w ⊤ ϕ(x) , (1.1)
where x is the input signal to the A-layer. This signal is a vector of numerical values that represent the object to be recognized. Hereafter, any tensor of numerical values which identify a measurable property is referred to as a feature. The input features x are transformed by some non-linear function ϕ (e.g., a sigmoid). Finally, the variable w stands for the (adaptive) weights of the A-layer1 and h is a step function such that

h(z) :=    +1 if z ≥ 0, -1 if z < 0.
Figure 1.2: Schematic representation of the original (photo-)perceptron. Layers S, A and R refer to the sets of sensory, association and recognition units/cells, respectively. In the example above, w 1 , w 2 , w 3 , w 4 are the weights of the A-layer and h is the threshold function accounting for neuron "firing". The threshold function h accounts for neuron "firing" for a specific pattern on screen.

The perceptron learning algorithm can be described as follows: for a given pattern x on screen, the output signal y is computed and, if the pattern is classified incorrectly (for instance, the neuron "fired" while it should not), a change of the weights values w take place. More precisely, the following updating rule is called for each misclassified pattern x in a dataset:

w ←-w + η ϕ(x) t, (1.2)
where t is an identifier for the class of the pattern and η is some step size.

The perceptron was designed for binary classification, so that t ∈ {-1, +1}.

If the dataset the model learns from is linearly separable2 , the perceptron convergence theorem as stated in [START_REF] Rosenblatt | The perceptron: A probabilistic model for information storage and organization in the brain[END_REF] guarantees that the learning algorithm above will find an exact solution in a finite number of steps. However, in practice, it can be difficult to distiguish a nonseparable problem from one that is simply slow to converge. Besides, many solutions to a separable problem can exist (as shown in Fig. 1.4), and a perceptron's converged state will depend on the weights initialization. Besides these limitations, Rosenblatt's perceptron as a neurophysiological model was criticized by [START_REF] White | Principles of Neurodynamics: Perceptrons and the Theory of Brain Mechanisms[END_REF], which stated that it lacked detailing of the temporal characteristics of a neurone firing and that most of the reported experiments dealt with figures that illuminate a large fraction of the retina whereas the device performance in other cases remained uninvestigated. The most impactful review of this work came only later on in 1969, though: [START_REF] Minsky | Perceptrons: An introduction to computational geometry[END_REF] in their famous book Perceptrons show that the original perceptron was able to learn only linearly separable problems. Despite acknowledging that increasing the number of A-layers could help solving non-linear ones, the reflections developed inside this book led to a substancial decline in research funding for neural networks.

The period between 1969 and the early 1980s is usually recognized as an Artificial Intelligence (AI) winter, where the field was surrounded by pessimism and criticism. Renewed interest would come after the works of [START_REF] Hopfield | Neural networks and physical systems with emergent collective computational abilities[END_REF] on associative memory systems. Hopfield networks are models that can reconstruct data after being fed with partial versions of the same data. For example, if the dataset some Hopfield model has learned from is made of photograph samples, the given model exposed to a piece of one of these photographs is capable of reconstructing the original. It was theoretically proven that a Hopfield network can reliably reconstruct up to 0.138×N samples, where N is the number of learning units in the model [START_REF] Stuart Russell | Artificial Intelligence: A Modern Approach. 3rd[END_REF]. Finally, the works of [START_REF] Rumelhart | Learning representations by back-propagating errors[END_REF] come as a major milestone in AI research. They proposed a new architecture called the multilayer perceptron, which is depicted in the next section.

Multilayer perceptrons

As opposed to Rosenblatt's original system, currently named single-layer perceptron3 , multilayer perceptrons (MLPs) have at least three A-R sets of layers. In addition to this, the threshold function4 of each perceptron composing an MLP is arbitrary. Finally, MLPs can either perform (multiclass) classification or regression, while Rosenblatt's perceptron was restricted to binary classification tasks. By solving the XOR problem (to predict outputs of exclusive or logic gates) and by detecting input symmetry (to determine if two strings are symmetric about their center), among other problems (Rumelhart, McClelland, et al., 1988), MLPs have re-ignited interest in neural networks for machine learning. Fig. 1.5 illustrates an MLP in the modern sense. For the sake of rigorousness, the reader will find below a formal definition of the MLP architecture : Definition 1 (Multilayer perceptron (MLP)). An N -layer multilayer perceptron is the composition

d 1 • d 2 • . . . • d N (1.3) of functions d k : R N in k --→ R N out k x --→ σ k (w k x + b k ), (1.4) (k ∈ 1, N , N > 2) (N in k ∈ N, N out k ∈ N)
where σ k is a non-linear activation function,

w k ∈ R N out k × N in
k is a matrix of trainable weights and b k ∈ R N out k is a vector of trainable biases. Each d k is named a dense layer of the MLP. In order to exist, the composition (1.3) must obey the constraint

N in k = N out k-1 (k > 1).
The layer d 1 is called the input layer, d N is named the output layer, and the layers d k with k ∈ 2, N -1 are referred to as the hidden layers of the MLP.

In the particular case where

N in k = N out k = N feat ∀k ∈ 2, N -1 ,
one says that N feat is the number of hidden features of the MLP.

Figure 1.5: Schematic representation of a 3-layer MLP with its sets {w idx_layer, idx_output, idx_input } of weights and {b idx_layer, idx_output } of biases. Connections between nodes of the computational graph store weights that take part in the matrix-vector multiplication explicited in the bottom part of the image. Since all nodes from one layer are linked to all nodes of the subsequent one, they are usually referred to as densely-connected or simply dense layers.

Image recognition, the primary task addressed by Rosenblatt's perceptrons can instead be tackled by MLPs. For instance, [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] showed that a simple 2-layer MLP can reach an accuracy level of 95.3% on a 10-label classification problem to recognize handwritten digits. The famous MNIST (Modified National Institute of Standards and Technology) dataset which contains 70, 000 samples of 28 × 28 pixel images of digits written by American Census Bureau employees and high-school students was used to train the MLPs.

Convolutional neural networks

It is possible to constrain MLPs to solve the MNIST recognition problem more efficiently. These constrained networks were termed convolutional neural networks (CNNs). The idea is to group spatially correlated pixels from the images and process them through shared weights. Fig. 1.7 illustrates the concept with an example. CNNs employ fewer trainable parameters than MLPs to reach similar accuracy levels in an image recognition task. On the one hand, each of the output pixels of a dense layer is a function of all of its input pixels. On the other hand, each of the output pixels of a convolutional layer is a function of a small amount of its input pixels. More precisely, the value of a pixel in the output of a convolutional layer is influenced only by a local neighborhood of pixels from the input. This inherent bias embedded in such CNN architectures towards the local neighborhood of pixels is called the locality inductive bias. Furthermore, the trainable kernel and biases of a convolutional layer are independent of a pixel's position in an image. CNNs transform portions of an image independent of their location in the picture itself5 , a property called spatial translation invariance.

A formal definition of a CNN follows:

Definition 2 (Convolutional neural network (CNN)). An N -layer convolutional neural network is the composition

c 1 • c 2 • . . . • c N (1.5) of functions c k : R D in k --→ R D out k x --→ c k (x)[m, n, q] := σ k   b k [q] + ∆ x k i=-∆ x k ∆ y k j=-∆ y k p w k [i, j, p, q] x[m + i, n + j, p]    , (1.6) (m ∈ 1, N x , n ∈ 1, N y , p ∈ 1, N in k , q ∈ 1, N out k ) (k ∈ 1, N , N > 2) (D in k ≡ N x × N y × N in k ) (D out k ≡ (N x -∆ x k + 1) × (N y -∆ y k + 1) × N out k ) (N in k ∈ N, N out k ∈ N, N x ∈ N, N y ∈ N, ∆ x k ∈ N, ∆ y k ∈ N)
where σ k is a non-linear activation function, N x and N y are the number of pixels composing the image along the x and y directions, respectively,

w k ∈ R (2∆ x k +1) × (2∆ y k +1) × N in k × N out k is a trainable convolutional kernel/filter and b k ∈ R N out k
is a vector of trainable biases. Each c k is named a convolutional layer of the CNN. In order to exist, the composition (1.5) must obey the constraint

N in k = N out k-1 (k > 1).
The layer c 1 is called the input layer, c N is named the output layer, and the layers c k with k ∈ 2, N -1 are referred to as the hidden layers of the CNN.

In the particular case where

N in k = N out k = N feat ∀k ∈ 2, N -1 ,
one says that N feat is the number of hidden features of the CNN.

Figure 1.7: Schematic comparison between a dense layer (top) and a convolutional layer (bottom). A 3 × 3 image with pixel values {x 1 , . . . , x 9 } is given as input and converted to an 8-pixel representation with values {y 1 , . . . , y 8 }. While the dense layer requires all-to-all connections between input and output neurons, the convolutional layer shares a relatively small set of weights (also known as convolutional kernel/filter) among spatially correlated pixels of the input image (highlighted in red). This allows CNNs to perform image recognition tasks with far fewer trainable parameters than corresponding MLPs.

In a scientific context, many applications for CNNs can be found. In the realm of premixed turbulent combustion modelling, for instance, [START_REF] Lapeyre | Training convolutional neural networks to estimate turbulent sub-grid scale reaction rates[END_REF] trained CNNs to estimate sub-grid flame surface density. Their CNN model was able to efficiently extract flame topology and predict subgrid scale wrinkling, outperforming classical algebraic models. In the field of composite materials, [START_REF] Ejaz | Convolutional neural networks for approximating electrical and thermal conductivities of Cu-CNT composites[END_REF] trained CNNs to map images of copper-carbon nanotube networks to corresponding electrical and thermal conductivities. They highlight the runtime savings obtained from running a CNN surrogate model for conductivity as opposed to a convectional numerical simulation. Regarding data-driven discretizations, the main focus of this thesis, [START_REF] Bar-Sinai | Learning data-driven discretizations for partial differential equations[END_REF] and [START_REF] Kochkov | Machine learning-accelerated computational fluid dynamics[END_REF] employed CNNs to approximate spatial derivatives in coarse grids based on downsampled dynamics from high-resolution numerical simulations. In Chapter 2, the reader will find an application of this method. It is restricted to structured grids, though, where instantaneous physical fields can feed a CNN as simple images. In order to overcome this restriction, a neural network architecture that handles grids as graphs can be used. The next section is dedicated to describing this architecture, particularly suitable for industrial CFD where unstructured grids are often used. [START_REF] Battaglia | Relational inductive biases, deep learning, and graph networks[END_REF] presented the CNN properties of spatial translational invariance and locality inductive bias within a more general framework. In such a framework, learning concerns the representation of entities and their relations. In the specific case of CNNs, the pixels of an image are the entities whose relations are synthesized by the locality inductive bias and the spatial translation invariance properties. In the more general case where the relations between entities are arbitrary, they proposed employing graphs for their representation. The set of nodes composing a graph stand for the entities, while its set of edges represent their relations. The resulting ML model destined to learn from the entities/relations representation in a graph was named graph neural network (GNN). In a GNN, the relations between entities are encoded as messages carried by the edges. These messages are exchanged between nodes and transformed across the layers of a GNN, which are termed message-passing layers. Tbl. 1.1 summarizes the entities and their relations for the neural network models defined so far. Finally, Fig. 1.8 depicts the operations performed in a message-passing layer of a GNN. et al. (2018).

Graph neural networks

Figure 1.8: Schematic representation of operations performed in a messagepassing layer. First, update of the edge features (i.e., the messages) takes place by the action of ϕ edge on the current features of each edge, along with the features of the sender and receiver nodes composing it. Last, ϕ node updates node features by acting on the current features of each node along with the updated features of the incoming edges (i.e., the updated incoming messages).

The building blocks ϕ edge and ϕ node are MLPs or CNNs, typically.

A formal definition of a GNN follows:

Definition 3 (Graph neural network (GNN)). An N -layer graph neural network is the composition

g 1 • g 2 • . . . • g N (1.7)
of functions

g k : R D in, node k , R D in, edge k --→ R D out, node k , R D out, edge k (x node , x edge ) --→ g k (x node , x edge ) := (ϕ edge k • ϕ node k )(x node , x edge ), (1.8) (k ∈ 1, N , N > 2) (D in, node k ≡ N node × N in, node k ) (D in, edge k ≡ N edge × N in, edge k ) (D out, node k ≡ N node × N out, node k ) (D out, edge k ≡ N edge × N out, edge k ) (N in, node k ∈ N, N out, node k ∈ N, N in, edge k ∈ N, N out, edge k ∈ N) (N node ∈ N, N edge ∈ N)
where N node and N edge are the number of nodes and edges in the graph, respectively, ϕ edge k is a building block responsible for updating edge features and ϕ node k is a building block responsible for updating node features. More precisely, ϕ edge k and ϕ node k are independent MLPs (as per Def. 1) or CNNs (as per Def. 2). The reader is referred to [START_REF] Battaglia | Relational inductive biases, deep learning, and graph networks[END_REF] for details on the actions of ϕ edge k and ϕ node k . Each g k is named a message-passing layer of the GNN. In order to exist, the composition (1.7) must obey the constraints

N in, node k = N out, node k-1 (k > 1) N in, edge k = N out, edge k-1 (k > 1).
The layer g 1 is called the input layer, g N is named the output layer, and the layers g k with k ∈ 2, N -1 are referred to as the hidden layers of the GNN. In the particular case where N in, node

k = N out, node k = N in, edge k = N out, edge k = N feat ∀k ∈ 2, N -1 ,
one says that N feat is the number of hidden features of the GNN.

GNNs have found successful application to problems dealing with unstructured data. Social networks for instance are naturally represented by graphs: [START_REF] Fan | A Characteristic-Featured Shock Wave Indicator for Conservation Laws Based on Training an Artificial Neuron[END_REF] built social recommender systems based on GNNs. In a scientific context, Deepmind's AlphaFold [START_REF] Jumper | Highly accurate protein structure prediction with AlphaFold[END_REF] employed GNNs to accurately predict protein 3D structure resulting from an amino acid sequence. [START_REF] Belbute-Peres | Combining Differentiable PDE Solvers and Graph Neural Networks for Fluid Flow Prediction[END_REF] accelerate predictions of fluid flows around airfoils by combining GNNs with a differentiable CFD solver (namely, SU2 by [START_REF] Economon | SU2: An Open-Source Suite for Multiphysics Simulation and Design[END_REF]). Unstructured meshes can be naturally assimilated into graphs. This motivates the use of GNNs for developing data-driven numerical schemes on 1D irregular and 2D unstructured grids in Chapters 4 and 5.

The ML models defined so far can learn to perform a given task by adjusting the values of their weights and biases, a process called training. Learning methods can be organized in three main paradigms [START_REF] Jordan | Machine learning: Trends, perspectives, and prospects[END_REF]: supervised, unsupervised and reinforcement learning. In supervised learning, the ML model learns by comparing its predictions to available reference data. In unsupervised learning, no reference data exist, and the model learns only by analyzing its inputs under assumptions about its structural properties (e.g., algebraic, combinatorial or probabilistic). Finally, in reinforcement learning, reference data is absent as well but the inputs contain indications as to whether an action of the ML model is correct or not. The development of data-driven numerical schemes in Chapters 2, 4, and 5 makes use of supervised learning because either the analytical solution dynamics is available for the cases studied or reference data can be obtained from a highfidelity numerical solver. The next section explains the training procedure in the supervised learning paradigm.

The supervised learning paradigm

In supervised learning, the model parameters are optimized by minimizing a cost, which is a measure of the distance between the network prediction to some reference value. Such a cost is a function of the neural network trainable parameters and is called the loss function in the ML community. The training procedure modifies the values of the trainable parameters in order to minimize the mean loss function value across a dataset. The optimization problem could be formulated as

W ⋆ := argmin W 1 N s Ns i=1 L(W, u i , r i ) , (1.9)
where W is the set of trainable parameters of the neural network. evaluates the gradient of the loss function with respect to the neural network parameters in a scheme known as the backpropagation algorithm. The name comes from the fact that the gradient computation starts from the output and propagates back to the inputs of the neural network model. Once this gradient is computed, an update rule is applied to the trainable parameters following the opposite direction of the gradient (similar to Eq. (1.2)). In its simplest form,

W ←-W -η • ∂ ∂W 1 N s Ns i=1 L(W, u i , r i ) , (1.10)
where η is the gradient descent step size, also known as the learning rate of the training procedure. The process is repeated iteratively until specific stopping criteria are met. Given that training sets can reach sizes of the order of thousands (as MNIST) to millions (as ImageNet by [START_REF] Russakovsky | ImageNet Large Scale Visual Recognition Challenge[END_REF] for generic categorical object recognition) of samples, applying the updating rule (1.10) is unfeasible in practice. An approximation of the gradient descent procedure takes place by computing gradients only for a randomly selected batch of samples at a time. Let U b ⊂ U be a batch of training samples of size N b and R b ⊂ R the corresponding batch reference set. The trainable parameters are updated as

W ←-W -η • ∂ ∂W   1 N b N b j=1 L(W, u j , r j )   .
(1.11)

When the batch size N b = 1, meaning that the trainable parameters W are updated one sample at a time, the learning process is said to occur by online training, and mainly relates to systems where data is generated as a function of time (e.g., stock price prediction as shown in [START_REF] Hoi | Online learning: A comprehensive survey[END_REF]).

In this case, the method that uses the updating rule (1.11) to solve the optimization problem (1.9) is called stochastic gradient descent [START_REF] Goodfellow | Deep Learning[END_REF]. When the batch size N b > 1, the method is called mini-batch gradient descent. In mini-batch gradient descent, reducing training time is made possible by distributing the computation of the gradients of the loss function among multiple processors, each of which handling one sample from the batch. Parallelization across the batch dimension is a common feature among modern AD tools in ML frameworks (e.g., TensorFlow and PyTorch, see Sec. 1.3). Random batches U b are sampled without replacement from the training set U until the set is exhausted, which is called an epoch. Typically, a training procedure involves running multiple epochs, and the training is terminated when a specific stopping criterion is met. A widely used stopping criterium compares the mean loss function value of the training set U with that of a different dataset, which is called the validation set (noted V herein). During training, the model is expected to improve performance metrics for both sets U and V. However, after some point, it will keep improving them only for the training set U. Generalization is an ML model's capacity to perform well on data it has never seen before [START_REF] Chollet | Deep learning with Python[END_REF]. The neural network is losing its generalization capacity from the point where its performance for the validation set V starts degrading. The training procedure must stop before that happens. Fig. 1.9 illustrates how the performance metrics of sets U and V may evolve along epochs. Finally, Alg. 1 summarizes the training procedure described throughout this section.

Accurate gradient computation is essential to obtain a reliable neural network model from the training procedure. An efficient technique to accurately evaluate gradients for computer programs is algorithmic differentiation, described in the next section. Figure 1.9: Performance metrics for training and validation sets at the end of each epoch during the training procedure. In (a), the neural network model can still potentially improve its predictions for both sets, in which case it is said to be underfitted. In (b), after a given point, it worsens validation set predictions for a number of epochs, in which case it is said to be overfitting on the training set. The optimal set of trainable parameters for the neural network can be found then just before the model starts overfitting, and is indicated by a star in the figure. for sample ∈ batch do 10: Algorithmic Differentiation (AD) comprises a collection of methods for converting a program that computes numerical values of a function into a program that computes numerical values of the function's derivatives. These derivatives can include first-order derivatives, such as the function's gradient or the Jacobian of a set of constraints, higher-order derivatives like the Hessian matrix multiplied by a directional vector, or even nested derivatives [START_REF] Bartholomew-Biggs | Automatic differentiation of algorithms[END_REF]. In AD, the derivatives are computed with approximately the same accuracy and efficiency as the original function evaluations. This is what distinguishes AD from two other ways of computing derivatives numerically. They are [START_REF] Baydin | Automatic Differentiation in Machine Learning: a Survey[END_REF] • Numerical Differentiation (ND): it consists in evaluating the program at some sample points and using finite difference formulas to obtain approximations of its derivatives. For example, the method dfdx defined as

L ← L + 1 N b L(W,
1 function dfdx(x0 :: Float64) 2 h = 1e-6 3 return ( f(x0 + h) -f(x0 -h) ) / 2h 4 end
uses ND to approximate the derivative of the method f representing a scalar map f : R → R at some real value x 0 . More precisely, dfdx makes use of a second order finite difference approximation:

f ′ (x 0 ) = f (x 0 + h) -f (x 0 -h) 2h + O(h 2 ).
(1.12)

The approximation (1.12) results from the truncation of the Taylor series expansions of f (x 0 + h) and f (x 0 -h) around the value x 0 , and the error associated with it is called the truncation error. Asymptotically, as h → 0, the truncation error diminishes. On the other hand, f(x0 + h)f(x0 -h) becomes closer to machine precision, which introduces the so-called round-off error. It grows as h → 0. The presence of both round-off and truncation errors is one of the main issues behind ND, since their asymptotic behaviors are opposite to each other. Fig. 1.10 illustrates it.

Round-off error dominant

Truncation error dominant

Figure 1.10: Absolute error of the ND-evaluated derivative of the function f (x) := exp(3x) at x = 0.2 using the approximation (1.12).

• Symbolic Differentiation (SD): it makes use of hard-coded differentiation rules for a finite set of operations and functions (e.g., +, -, ×, ÷, log, sin, . . . ). This way, SD does not suffer from the numerical errors present in ND. To illustrate it, one can consider the real-valued function f (x) := exp(3x). In ND, truncation and round-off errors are present for employing truncated Taylor series expansions (see approximation (1.12) and Fig. 1.10). In SD, on the other hand, one can simply code the differentiation rule for exponentials and avoid these sources of error. In SD, the derivative of f is directly evaluated from its analytical expression f ′ (x) = 3 exp(3x). Nevertheless, SD presents some limitations. The program to be differentiated must not contain control flow mechanisms (e.g., loops, conditionals). Moreover, careless application may lead to increasingly long symbolic expressions and performance issues, a phenomenon known as expression swell. To illustrate it, one can consider the function [START_REF] Speelpenning | Compiling fast partial derivatives of functions given by algorithms[END_REF]. The gradient of g is expressed as:

g(x 1 , . . . , x n ) := n i=1 x i , (1.13) with n > 1 (∀i ∈ 1, n , x i ∈ R)
∂g ∂x i =       x 2 x 3 . . . x n x 1 x 3 . . . x n . . . x 1 x 2 . . . x n-1      
.

(1.14)

In a (naive) implementation of SD, one builds the symbolic expression for Eq. (1.14) without considering common subexpressions. In this case, its memory allocation grows exponentially with n.

In order to avoid the numerical errors present in ND and the expression swell issue in (naive) SD, AD in its basic form proposes to break the computer program down into elementary pieces to which SD will be applied and store intermediate numerical results whenever necessary. By interpreting the program as a composition of differentiable pieces, the chain rule of calculus can be applied to derive gradients for the entire chain of operations. [START_REF] Laue | On the Equivalence of Automatic and Symbolic Differentiation[END_REF] states that there is an equivalence of AD and SD when the implementation of SD stores intermediate results as well as pointers to identify common subexpressions.

In practical implementations of AD, track records of operations called evaluation traces (also known as [START_REF] Wengert | A Simple Automatic Derivative Evaluation Program[END_REF] lists) are built. The chain of operations constituting the computer program to be differentiated is referred to as the forward primal evaluation trace (or forward pass). Fig. 1.11 shows the forward primal evaluation trace for a sample program. To compute gradients, two distinct evaluation traces can be built: (i) the forward tangent (Sec. 1.2.1) and (ii) the adjoint (Sec. 1.2.2). Here, the variable naming convention of [START_REF] Griewank | Evaluating Derivatives[END_REF] is used: v i-N with i ∈ 1, N are the N input variables, v j with j ≥ 1 are the so-called intermediate variables and y k with k ∈ 1, M are the M output variables of the evaluation trace.

Figure 1.11: Computational graph and forward primal evaluation trace for a program evaluating f (x 1 , x 2 ) := sin(x 1 + x 2 ), where x 1 , x 2 ∈ R.

Forward mode of AD

Forward mode AD is best suited for the differentiation of programs where the number of output variables is greater than the number of input variables, being mostly used for sensitivity analysis of physical systems. For example, T. W. [START_REF] Hughes | Forward-Mode Differentiation of Maxwell's Equations[END_REF] applied the forward mode of AD on Maxwell's equations to compute the sensitivity of a spatial near-field intensity distribution (defined for each cell of a 2D domain) to the dielectric constant of the medium (a scalar value, typically). This mode evaluates the derivatives

vj, i := ∂ v j ∂ v i-N (j ≥ 1, i ∈ 1, N ) (1.15) ẏk, i := ∂ y k ∂ v i-N (k ∈ 1, M ) (1.16)
of a given forward primal evaluation trace. In other words, the forward mode of AD computes the derivatives of all the intermediate and output variables with respect to each of the input variables. The resulting evaluation trace is called forward tangent. For the sake of simplicity, whenever the number of input variables N = 1, the second index of the derivatives vj, i and ẏk, i is dropped. Hence, they are represented simply as vj and ẏk , respectively.

As an example, one can consider the sliding ladder problem: it consists on placing a ladder as a rigid and uniform rod in a 2D plane that lays over two perpendicular frictionless walls, as illustrated in the top part of Fig. 1.12. [START_REF] Kapranidis | Variations of the Sliding Ladder Problem[END_REF] showed that as the ladder falls under the action of gravity, there will be a separation point at which it loses contact with the vertical wall (wall normal diminishes). The separation angle and the speed at which the rod hits the floor are respectively given by θ sep (α) = arccos 2 3 cos α (1.17)

s f (α) = gℓ 6 cos α - 4 9 cos 3 α , (1.18)
where g, ℓ and α stand for local gravity, ladder length and initial ladder angle with the wall, respectively. For fixed values of g and ℓ, one may be interested in knowing the sensitivity of {θ sep , s f } to the parameter α. Equivalently, it means evaluating how much a perturbation in the value of α affects the values of {θ sep , s f }. This boils down to computing the gradients ∂θ sep /∂α and ∂s f /∂α. They are computed at the end of the forward tangent evaluation trace shown at the bottom of Fig. 1.12. Finally, the number of forward tangent passes equals the number N of input variables, i.e., it is independent of the number M of output variables [START_REF] Baydin | Automatic Differentiation in Machine Learning: a Survey[END_REF].

Figure 1.12: (Top) Schematic representation of the sliding ladder problem. (Middle and bottom) Computational graph, forward primal, and forward tangent evaluation traces for a program evaluating the angle y 1 (θ sep in Eq. (1.17)) at which the falling ladder loses contact with the vertical wall and the speed y 2 (s f in Eq. (1.18)) at which it hits the floor.

Reverse mode of AD

The Finnish mathematician and computer scientist Seppo Linnainmaa introduced the reverse mode of AD for efficiently computing derivatives of composite functions in his master's thesis in 1970 in Finnish [START_REF] Linnainmaa | The representation of the cumulative rounding error of an algorithm as a Taylor expansion of the local rounding errors[END_REF]. The English speaking world was largely unaware of the method until the publication of the English translation (Linnainmaa, 1976) 7 . In this paper, he applied reverse mode of AD to evaluate error propagation in floating point arithmetic. The reverse mode of AD is best suited for differentiating programs where the number of input variables is much higher than the number of output ones. This is the case for neural networks, where one is interested in minimizing a scalar-valued loss function by changes in the values of the trainable weights and biases. The number of trainable parameters reached the outstanding figure of 530 billion for Megatron-Turing NLG 530B, a large-scale generative language model [START_REF] Smith | Using DeepSpeed and Megatron to Train Megatron-Turing NLG 530B, A Large-Scale Generative Language Model[END_REF]. The first application of reverse mode of AD to neural networks training is credited to [START_REF] Werbos | Beyond regression: new tools for prediction and analysis in the behavioral sciences[END_REF]. This mode evaluates the derivatives

v k, j := ∂ y k ∂ v j (k ∈ 1, M , j ≥ 1 -N ) (1.19)
of a given forward primal evaluation trace. In other words, the reverse mode of AD computes the derivatives of each of the output variables with respect to all the intermediate and input variables. The resulting evaluation trace is called adjoint and each v k,j is named a reverse derivative of this trace.

For the sake of simplicity, whenever the number of output variables M = 1, the first index of the reverse derivatives v k, j is dropped. Hence, they are represented simply as v k .

As an example, a dense layer (as stated inside Def. 1) is considered. It outputs a scalar value (a loss function value, typically) and is composed of the set {w 1 , w 2 , w 3 } ∈ R 3 of (trainable) weights and b 1 ∈ R as (trainable) bias. Given {x 1 , x 2 , x 3 } ∈ R 3 as input, the output of such a dense layer can be expressed as

L(w 1 , w 2 , w 3 , b 1 ) = σ(w 1 x 1 + w 2 x 2 + w 3 x 3 + b 1 ), (1.20)
where ∀x ∈ R, σ(x) := max(0, x). The activation function σ is called Rectified Linear Unit (ReLU). It was introduced by [START_REF] Fukushima | Visual Feature Extraction by a Multilayered Network of Analog Threshold Elements[END_REF] in the context of visual feature extraction using neural networks. The function σ is not differentiable everywhere. By definition, its derivative σ ′ is expressed as

σ ′ (x) :=    1 if x > 0 0 if x < 0, (1.21)
and has no defined value at x = 0. However, when using a computer program to execute the adjoint evaluation trace, one needs to define a numerical value for the derivative σ ′ at x = 0. The ML libraries TensorFlow, PyTorch, and JAX (see Sec. [START_REF] Bertoin | Numerical influence of ReLU'(0) on backpropagation[END_REF] reported that the choice of value for σ ′ at x = 0 influences the training results, and that the choice σ ′ (0) = 0 provides the best accuracy levels for neural networks trained via stochastic gradient descent (see Sec. 1.1.4). Fig. 1.13 shows the computation of all reverse derivatives related to the computational graph of the dense layer whose output is expressed in Eq. (1.20).

1.3.1-1.3.2) enforce σ ′ (0) = 0.
Finally, the number of adjoint passes equals the number M of output variables, i.e., it is independent of the number N of input variables [START_REF] Baydin | Automatic Differentiation in Machine Learning: a Survey[END_REF]. This justifies its use in the machine learning field, where the loss function typically returns a scalar value and the number of trainable weights and biases is large.

The reverse mode of AD is more computationally efficient than the forward mode in the development of ML-based numerical schemes, which is the main focus of this thesis. In Sec. 1.3, the reader will find an overview of AD tools with reverse-mode support.

An overview of AD frameworks

This section will familiarize the reader with existing differentiation tools that could be applied to obtain the adjoint of codes composing ML systems and PDE solvers. Among them, TensorFlow and PyTorch focus more on the former, whereas the others better handle generic control flow typical to the latter.

TensorFlow & PyTorch (Python/C++)

TensorFlow is an open source platform for machine learning originally developed by the Google Brain team to conduct research in AI within the organization, and then publicly released under the Apache 2 license [START_REF] Metz | Google Just Open Sourced the Artificial Intelligence Engine at the Heart of Its Online Empire[END_REF]. PyTorch is an open source machine learning framework developed by Meta AI. It was released in 2016 under Modified BSD license (Release alpha-1 release • pytorch/pytorch 2022). They are grouped for employing model subclassing and object-oriented paradigm in ML pipelines. Besides, users need to become familiar with their Domain Specific Languages (DSLs), under which classic layers in the ML community (e.g., dense and convolutional ones) are pre-defined, but custom layers are created at the cost of a rather verbose code dependent on the availability of elemental building blocks by the DSL.

JAX (Python)

JAX is a just-in-time compiler of Python/NumPy-similar code released by the Google Brain team [START_REF] Frostig | Compiling machine learning programs via high-level tracing[END_REF]. Like TensorFlow and PyTorch, JAX introduces a DSL that is intended to be similar to plain Python/NumPy in order to ease more complex control flow differentiation, typical to coupled ML-PDE systems. However, it only handles pure functions, i.e., input arguments cannot be mutated in the scope of the method. This may turn the coding experience cumbersome.

Tapenade (Fortran/C)

Tapenade is a source-to-source AD engine released by [START_REF] Hascoët | The Tapenade Automatic Differentiation tool: Principles, Model, and Specification[END_REF]. Adjoint source code is accessible and can be modified for performance enhancement through parallelization. However, no pre-defined adjoints of typical building blocks of neural networks like dense and convolutional layers are available. This might repel ML users, especially from industry. The development of the tool toward ML applications would ease the integration of existing legacy industrial CFD codes (mainly written in Fortran/C languages) to neural networks and eventually open a new path to high-performance machine learning for scientific applications. The reader can find in Appendix A an example of Tapenade-generated adjoint of a program.

Zygote (Julia)

Zygote is a source-to-source differentiation tool in Julia [START_REF] Bezanson | Julia: A Fresh Approach to Numerical Computing[END_REF] released by Innes (2018a). It is the AD system of the language's main machine learning framework : Flux [START_REF] Innes | Fashionable Modelling with Flux[END_REF][START_REF] Innes | Flux: Elegant Machine Learning with Julia[END_REF]. It handles generic control flow in Julia (apart from array mutation, as explained in Appendix B) and coupling PDE solvers coded in Julia to neural networks is as easy as writing any other Julia code. It is used to differentiate dU_ttgc_reassignment in Sec. B.3.

Enzyme (LLVM Intermediate Representation)

Enzyme is an AD tool for differentiating code in LLVM Intermediate Representation (IR) released by [START_REF] Moses | Instead of Rewriting Foreign Code for Machine Learning, Automatically Synthesize Fast Gradients[END_REF]. LLVM is a compiler and a toolkit for building compilers that transforms source code from a myriad of languages (e.g., Fortran, C++, Swift, Rust, Julia) into an IR that is efficiently converted to machine code (The LLVM Compiler Infrastructure Project 2022). It was designed for hardware portability, as it offers machineindependent primitive types [START_REF] Yegulalp | What is LLVM? The power behind Swift, Rust, Clang, and more[END_REF]. Enzyme is used to differentiate dU_ttgc_mutation in Sec. B.3. This chapter detailed the concepts required to understand the ML components of the data-driven numerical schemes developed in Chapters 2, 4, and 5. The next chapter is dedicated to reviewing concepts of numerical methods for PDEs. They precede the spectral analysis proposed in Chapter 3, which helps to define stability constraints for these data-driven numerical schemes.

Chapter 2 Introducing data-driven numerical schemes

The numerical simulation of fluid flows makes use of numerical schemes, which are algebraic relations between the space, time, and flow quantities represented in a computer. This chapter presents classical results related to the numerical approximation of hyperbolic systems of conservation laws essential to Computational Fluid Dynamics (CFD). Among them, a theorem by the Russian mathematician Sergei K. [START_REF] Godunov | A difference method for the numerical computation of discontinuous solutions of hydrodynamic equations[END_REF] states that highly accurate linear numerical schemes cannot obey a property called monotonicity. This property ensures that the numerical scheme cannot generate new extrema other than those eventually present in the initial condition [START_REF] Hirsch | Numerical Computation of Internal and External Flows[END_REF]. Lack of monotonicity harm the reliability of the numerical simulation since it leads to unphysical values of quantities that must remain positive such as density, turbulent kinetic energy, or combustion mixture fraction. Spurious oscillations may also lead to unphysical values of certain quantities (T. K. [START_REF] Sengupta | Spurious waves in discrete computation of wave phenomena and flow problems[END_REF]. These oscillations are specifically named Gibbs' phenomena when originating from shock waves or sharply discontinuous solutions (T. K. [START_REF] Sengupta | High Accuracy Compact Schemes and Gibbs' Phenomenon[END_REF]). Godunov's theorem and the emergence of spurious oscillations (among them, Gibbs' phenomenon) in numerical simulations have significantly impacted the development of numerical schemes used in CFD. It should be highlighted that while state variables display discontinuity near shocks, the corresponding flux terms are continuous functions. Spectral accuracy of schemes is capable of curing Gibbs' phenomenon, as it is done by Essentially Non-Oscillatory (ENO) and Weighted Essentially Non-Oscillatory (WENO) schemes [START_REF] Zhang | ENO and WENO Schemes[END_REF] applied to the state variables by one-sided differentiation. ENO and WENO schemes perform particular reconstructions of numerical fluxes which are nonlinear. The neural network models introduced in the previous chapter are nonlinear functions that can be employed to design numerical schemes. Sec. 2.3 presents a state-of-the-art method that approximates the derivatives of the numerical fluxes using a constrained Convolutional Neural Network (CNN) (see Def. 2). The constrained CNN is trained to accurately reproduce the convection of isentropic vortices and to evolve a double shear layer initial condition following the Euler equations (an inviscid version of the Navier-Stokes equations) on 2D cartesian grids.

Generalities about the numerical approximation of hyperbolic conservation laws

As they relate to transport phenomena, conservation laws of hyperbolic nature play an important role in Computational Fluid Dynamics (CFD) [START_REF] Hirsch | Numerical Computation of Internal and External Flows[END_REF]. In their differential form (no source terms):

∂U ∂t + div F (U ) = 0, (2.1)
where U is the state variable, F is a convective flux, and t is the temporal coordinate. The first part of this chapter concerns only one-dimensional equations such that

( U , F ) ≡ ( u(x, t), c u(x, t) ) (2.2) ( U , F ) ≡ ( u(x, t), u(x, t) 2 /2 ), (2.3) where c ∈ R and u : (x, t) ∈ R×[0, ∞) → u(x, t) ∈ R. The equivalences (2.
2) and ( 2.3) respectively lead to the following Partial Differential Equations (PDEs):

u t + c u x = 0 (2.4) u t + u u x = 0, (2.5)
where the subscripts (•) t = ∂ ∂t and (•) x = ∂ ∂x denote the partial derivatives with respect to the time and space, respectively. Eq. (2.4) is referred to as the (linear) Convection equation. It models the transport of a quantity u by the constant convection velocity c. The Convection equation is also regarded as a wave propagation equation, where the wave amplitude u has a phase propagating speed equal to c. Eq. (2.5) is called the (inviscid) Burgers' equation, named after [START_REF] Burgers | A Mathematical Model Illustrating the Theory of Turbulence[END_REF]. It appears in studies of gas dynamics [START_REF] Shefter | Quasiperiodic Solutions in Weakly Nonlinear Gas Dynamics. Part I. Numerical Results in the Inviscid Case[END_REF] and of traffic flow [START_REF] Lighthill | On Kinematic Waves. II. A Theory of Traffic Flow on Long Crowded Roads[END_REF]. Its viscous counterpart 1

u t + u u x -ν u xx = 0 (ν > 0), (2.6)
can be viewed as a model for decaying free turbulence [START_REF] Cole | On a Quasi-Linear Parabolic Equation Occurring in Aerodynamics[END_REF].

The problem that consists in finding u := u(x, t) solutions to Eq. (2.4) or to Eq. (2.5) for a given initial condition u 0 := u 0 (x) is called an Initial Value Problem (IVP). Even when initial conditions are smooth, nonlinear hyperbolic conservation laws can develop jump discontinuities that propagate as shock waves [START_REF] Dafermos | Introduction to the Theory of Hyperbolic Conservation Laws[END_REF]. An important IVP studied in the theory of hyperbolic PDEs is the Riemann problem, which contains a discontinuity [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]. It is formally defined as follows:

Definition 4 (1D Riemann problem). Let u L ∈ R, u R ∈ R, and x 0 ∈ R such that u L < u R . The IVP with initial condition u 0 (x) =    u L if x ≤ x 0 , u R if x > x 0 (2.7)
is called a 1D Riemann problem. The variables u L and u R are referred to as the left and right states of the 1D Riemann problem, respectively.

One can also state a version of the problem for periodic domains as follows:

Definition 5 (1D periodic Riemann problem). Let u min ∈ R, u max ∈ R, x 0 ∈ R, x 1 ∈ R, and ℓ ∈ R such that 0 ≤ u min < u max and 0 < x 0 < x 1 < ℓ. The IVP with initial condition u 0 (x) =    u max if x ∈ [x 0 , x 1 ] , u min if x ∈ [0, x 0 ) ∪ (x 1 , ℓ] (2.8)
is called a 1D periodic Riemann problem. (2.4) andthe Burgers' equation (2.5). They are obtained using the method of characteristics described in Chapter I (Sec. 4.1) of [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]. On the one hand, the solution to any IVP following the Convection equation (2.4) simply corresponds to the propagation of the initial condition u 0 at the convection speed c, i.e.,

u(x, t) ≡ u 0 (x -ct) (t > 0).
On the other hand, the solution to the 1D periodic Riemann problem following the Burgers' equation (2.5) is given by

u(x, t) ≡        (x -x 0 )/t, x ∈ (x 0 + u min t, x 0 + u max t] u max , x ∈ [x 0 + u max t, x 1 + ct] u min , elsewhere (t ∈ (0, t m )), (2.9)
where c := (u max + u min )/2 and t m := 2(x 1 -x 0 )/(u max -u min ). The 1D periodic Riemann problem following the Burgers' equation (2.5) admits solutions other than Eq. (2.9). Nevertheless, this is the only physically admissible one for respecting an entropy condition defined by [START_REF] Lax | Shock Waves and Entropy[END_REF]. Eq. (2.9) is called an entropy solution. One can prove the uniqueness of the entropy solution to any IVP following the Burgers' equation (2.5) (see Thm. 5.4 in Chapter I of [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]). However, it remains an open problem in the case of nonlinear systems of hyperbolic conservation laws. The 1D Riemann problem is suitable for testing numerical methods that approximate Eqs. (2.4)-(2.5), thanks to the availability of analytical solutions. Its numerical simulation requires representing the state variable u using a finite number of values. The spatial and temporal simulation domains must also be represented using a finite number of points in space and time, respectively. Such a process is called discretization. A numerical scheme is a function that employs spatial and temporal discretizations to evolve the discrete values of the state variable. A formal definition follows: Definition 6 (Numerical scheme). Let I be the set of points representing a spatial discretization and {t n } n∈ 1,N be the set of N points representing a temporal discretization. A numerical scheme for the PDE (or system of PDEs) E is a function

S E : R m , R m|K i | --→ R m (U n i , {U n k } k∈K i ) --→ U n+1 i , (2.10)
where m ∈ N is the number of equations, K i ⊂ I is a subset of points of the spatial discretization and |K i | is the number of points in K i . The variable U n i represents the discrete value of the state variable U at time level t n and spatial point i ∈ I. Furthermore, S E is called a linear scheme if and only if

∃ c k ∈ R such that U n+1 i := S E (U n i , {U n k } k∈K i ) = k∈K i c k U n k .
Provided that the numerical solution approximates the solution of the considered PDE (or system of PDEs), it is important to quantify how the approximation error levels evolve with the number of points used in the spatial and temporal discretizations. In general, numerical schemes that require fewer discretization points to achieve a given approximation error value are preferred. In the literature, it is common to compare numerical schemes by their truncation error behaviors using the concept of order of accuracy [START_REF] Strikwerda | Finite difference schemes and partial differential equations[END_REF]. It can be defined for equally spaced cartesian grids as follows:

Definition 7 (Order of accuracy on equally spaced cartesian grids). Let U be a smooth solution to Eq. (2.1). A numerical scheme S E is said to be p-th order accurate in space and q-th order accurate in time if and only if

{U n+1 i,exact } i∈I -{S E (U n i , {U n k } k∈K i )} i∈I 2 = O ∆x p+1 + ∆t q+1 , (2.11) (as ∆t → 0, ∆x → 0)
where ∆x and ∆t denote the spatial and temporal discretization parameters, respectively. More precisely, ∆x is the cartesian grid spacing and ∆t = t n+1 -t n ∀n is a fixed time step size. Unless otherwise stated, p and q are the largest positive integers for which Eq. (2.11) holds. Finally, the left-hand side of Eq. (2.11) is called the truncation error of the numerical scheme.

Numerical schemes for hyperbolic conservation laws are also compared in terms of physical admissibility of the numerical solutions they provide. In gas dynamics, for instance, it is essential that numerical values for density remain positive. In reactive flows, mass fraction values must stay in the interval [0, 1]. An important property which helps ensuring that these values remain in a physically admissible interval is that of monotonicity. This property states that no new extrema will be generated by the numerical scheme, aside from those present in the initial condition. Nevertheless, [START_REF] Godunov | A difference method for the numerical computation of discontinuous solutions of hydrodynamic equations[END_REF] showed that one cannot formulate a monotone linear scheme for the 1D Convection equation (2.4) with spatial order of accuracy higher than one:

Theorem 1 (Godunov's order barrier). Let S linear Convection be a linear scheme for the 1D Convection equation (2.4) which is p-th order accurate in space. If p ≥ 2, then S linear Convection cannot be monotone.

Godunov's order barrier Thm. 1 has played an important role in the development of numerical schemes for the 1D Convection equation (2.4) and for the equations of fluid motion in CFD, by extension. The next section depicts numerical schemes of historical relevance.

Numerical analysis of schemes

A numerical scheme is said to be consistent if it tends to the considered PDE (or system of PDEs) as the discretization time and space steps tend to zero. Furthermore, it is said to be stable if the numerical solution remains bounded over time. This section is dedicated to show accuracy properties of consistent numerical schemes of historical importance. In what follows, only 1D equally spaced grids (with mesh spacing ∆x > 0) and uniform time discretizations (with time step ∆t > 0) are considered.

First Order Upwind (FOU) and Lax-Wendroff (LW)

The First Order Upwind (FOU) scheme proposed by [START_REF] Courant | On the solution of nonlinear hyperbolic differential equations by finite differences[END_REF] uses Taylor expansions of the perturbed states u n+1 i and u n i-1 in time and space, respectively. Firstly, one can notice that

u n+1 i = u n i + ∆t u n t | x=x i + O(∆t 2 ) (2.12)
leads to the following expression of the first order derivative in time u n t | x=x i :

u n t | x=x i = u n+1 i -u n i ∆t + O(∆t).
(2.13)

Similarly,

u n i-1 = u n i -∆x u n x | x=x i + O(∆x 2 ) (2.14)
leads to the following expression of the first order derivative in space u n x | x=x i :

u n x | x=x i = u n i -u n i-1 ∆x + O(∆x).
(2.15)

Finally, one obtains the first order accurate in time and space FOU scheme for the 1D Convection equation (2.4) as

S FOU 1D Convection : R, R 2 --→ R u n i , {u n k } k∈{i-1,i} --→ u n+1 i := u n i -N c (u n i -u n i-1 ), (2.16)
where N c := c ∆t/∆x. The FOU scheme can also be used to discretize the 1D inviscid Burgers' equation (2.5). It similarly applies Eq. (2.13) to express the time derivative u n t , and uses a Taylor expansion of the associated flux F ≡ f (u) := u 2 /2 ∀u ∈ R to express its spatial derivative. The FOU scheme for the 1D Burgers' equation (2.5) reads

S FOU 1D Burgers' : R, R 2 --→ R u n i , {u n k } k∈{i-1,i} --→ u n+1 i := u n i - ∆t 2∆x (u n i ) 2 -u n i-1 2 .
(2.17) [START_REF] Lax | Systems of conservation laws[END_REF] employed second-order Taylor expansions to produce a scheme of better accuracy levels when compared to first-order methods like the FOU scheme. In the Lax-Wendroff (LW) scheme, the firstorder time derivative u n t | x=x i is expressed as

u n t | x=x i = u n+1 i -u n i ∆t - ∆t 2 u n tt | x=x i + O(∆t 2 ).
(2.18)

The time derivative u n tt in Eq. (2.18) is replaced by an equivalent spatial derivative using the 1D Convection equation (2.4), i.e.,

u n tt ≡ -c (u n x ) t ≡ -c (u n t ) x ≡ c 2 u n xx implies that u n t | x=x i = u n+1 i -u n i ∆t - c 2 ∆t 2 u n xx | x=x i + O(∆t 2 ). (2.19)
Eq. ( 2.19) leads to the following semi-discrete form of the 1D Convection equation (2.4):

u n+1 i -u n i ∆t - c 2 ∆t 2 u n xx | x=x i + c u n x | x=x i + O(∆t 2 ) = 0. (2.20)
The first and second order spatial derivatives u n x and u n xx in Eq. (2.20) are discretized using the following second-order Taylor expansions:

u n x | x=x i = u n i+1 -u n i-1 2∆x + O(∆x 2 ), u n xx | x=x i = u n i+1 -2u n i + u n i-1 ∆x 2 + O(∆x 2 ).
Finally, the LW scheme for the 1D Convection equation (2.4) is expressed as

S LW 1D Convection : R, R 3 --→ R u n i , {u n k } k∈{i-1,i,i+1} --→ u n+1 i := u n i - N c 2 (u n i+1 -u n i-1 ) + N 2 c 2 (u n i+1 -2u n i + u n i-1 ). (2.21)
In order to derive the LW scheme for the 1D Burgers' equation (2.5), the time derivative u n tt in Eq. (2.18) is replaced by a spatial derivative of the associated flux

F ≡ f (u) := u 2 /2 ∀u ∈ R, i.e., 2 u n tt ≡ -(f n x ) t ≡ -(f n t ) x ≡ - ∂f n ∂u u n t x ≡ ∂f n ∂u f n x x
implies that

u n t | x=x i = u n+1 i -u n i ∆t - ∆t 2 ∂f n ∂u f n x x x=x i + O(∆t 2 ). (2.22)
Eq. (2.22) leads to the semi-discrete form of the 1D Burgers' equation (2.5) below:

u n+1 i -u n i ∆t - ∆t 2 ∂f n ∂u f n x x x=x i + f n x | x=x i + O(∆t 2 ) = 0.
(2.23)

The spatial derivative f n x in Eq. (2.23) is discretized using the following second-order Taylor expansion:

f n x | x=x i = f n i+1 -f n i-1 2∆x + O(∆x 2 ).
(2.24)

The spatial derivative

∂f n ∂u f n x x
in Eq. (2.23) can be discretized as follows [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]:

∂f n ∂u f n x x x=x i = 1 ∆x 2   ∂f n ∂u i+1/2 (f n i+1 -f n i ) - ∂f n ∂u i-1/2 (f n i -f n i-1 )   ,
(2.25)

where ∂f n ∂u i±1/2 := 1 2 ∂f n ∂u i + ∂f n ∂u i±1
is the average value of the flux Jacobian in the neighborhood {i, i ± 1}. Finally, by applying Eqs. (2.24)-(2.25) into Eq. ( 2.23), one can express the LW scheme for the 1D Burgers' equation (2.5) as follows:

S LW 1D Burgers' : R, R 3 --→ R u n i , {u n k } k∈{i-1,i,i+1} --→ u n+1 i := u n i - ∆t 4∆x u n i+1 2 -u n i-1 2 + ∆t 2 8∆x 2 (u n i+1 + u n i ) u n i+1 2 -(u n i ) 2 -(u n i + u n i-1 ) (u n i ) 2 -u n i-1 2 .
(2.26)

The LW scheme enjoys less dissipation than the FOU scheme. However, Gibbs' phenomenon is observed when evolving initial conditions with steep gradients, as shown in Fig. 2.4. It appears as a consequence of Godunov's order barrier Thm. 1. In order to overcome the limitations predicted by the theorem, numerical schemes for the Convection equation (2.4) must be nonlinear. Sec. 2.2.2-2.2.3 introduce nonlinear schemes of historical relevance.

Monotonic Upstream-centered Schemes for Conservation Laws (MUSCL)

Van [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF] introduced the class of Monotonic Upstream-centered Schemes for Conservation Laws (MUSCL), which can achieve spatial second order of accuracy while preventing the emergence of Gibbs' phenomenon in the presence of shocks or sharp discontinuities. Given the semi-discrete form of the hyperbolic conservation law (2.1) below:

u n t | x=x i + L(u n i ) = 0,
where

L(u n i ) := 1 ∆x f i+1/2 (u n i ) -f i-1/2 (u n i ) ≡ 1 ∆x f n i+1/2 -f n i-1/2 ,
MUSCL seek suitable nonlinear approximations for the numerical fluxes f n i±1/2 , defined at the boundaries of a cell surrounding the mesh point i (see Fig. 2.3). The computation of these numerical fluxes depend on values of the state variable extrapolated to the right and left sides of the cell's boundaries, and are denoted as u R, n i±1/2 and u L, n i±1/2 , respectively. The extrapolations u R, n i±1/2 and u L, n i±1/2 are nonlinear functions of the state variable values in the neighborhood {i -

1, i, i + 1, i + 2}. More precisely, u R, n i+1/2 := u n i+1 -ϕ(r n i+1 ) u n i+2 -u n i+1 2 , u L, n i+1/2 := u n i + ϕ(r n i ) u n i+1 -u n i 2 ,
where

r n i := (u n i -u n i-1 )/(u n i+1 -u n i
) is an indicator of the local smoothness of the numerical solution and ϕ is a function of this smoothness indicator that limits the slope of the piecewise linear approximation of the state variable inside the cell. The function ϕ is commonly referred to as flux limiter. Flux limiters help to prevent the emergence of Gibbs' phenomenon in the numerical solution. i±1/2 and u L, n i±1/2 of the state variable at the boundaries of the cell [x i-1/2 , x i+1/2 ], where x i±1/2 := (x i + x i±1 )/2. Fig. 2.4 shows numerical experiments performed with the second-order accurate in space MUSCL of [START_REF] Kurganov | New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations[END_REF], for which

f i±1/2 (u n i ) := 1 2 f u R, n i±1/2 + f u L, n i±1/2 -a i±1/2 u R, n i±1/2 -u L, n i±1/2 , (2.27)
where a i±1/2 is the maximum absolute value of the eigenvalues of the flux Jacobian ∂f n /∂u. The generalized minmod flux limiter of Van Leer (1979) defined as

ϕ θ (r n i ) := max [0, min (θ, θr n i , 0.5(1 + r n i ))] (θ ∈ [1, 2]) (2.28)
is used with θ = 2 (value of the parameter for which the dynamics is least dissipative). To evolve the numerical solution over time, the three-stage, third-order Strong-Stability-Preserving (SSP) explicit Runge-Kutta time discretization method of Gottlieb et al. ( 2001) is used here. It is given by

u (1) i = u n i + ∆t L(u n i ), u (2) i = 3 4 u n i + 1 4 u (1) i + 1 4 ∆t L u (1) i ,
(2.29)

u n+1 i = 1 3 u n i + 2 3 u (2) i + 2 3 ∆t L u (2) i .

Weighted Essentially Non-Oscillatory (WENO)

The class of Weighted Essentially Non-Oscillatory (WENO) schemes was designed to achieve arbitrarily high order of accuracy in smooth regions while also resolving singularities such as discontinuities and sharp gradients in an accurate and non-oscillatory fashion [START_REF] Zhang | ENO and WENO Schemes[END_REF]. The fifth-order accurate in space WENO (WENO5) scheme of [START_REF] Jiang | Efficient Implementation of Weighted ENO Schemes[END_REF] has been used in most applications and is considered here.

Similarly to MUSCL, WENO schemes look for suitable forms of the numerical flux f n i±1/2 at cells' boundaries. They reconstruct the numerical flux at a given boundary as a weighted sum of the flux values at neighboring mesh points. The set of weights is obtained as a function of indicators of the local smoothness of the solution, which can adapt the reconstructed flux to accurately capture steep gradients and shocks. The WENO5 method assumes the final form below:

u n t | x=x i + L(u n i ) = 0,
where

L(u n i ) := 1 ∆x f i+1/2 (u n i ) -f i-1/2 (u n i ) ≡ 1 ∆x f n i+1/2 -f n i-1/2 ,
with the flux decomposition

f n i+1/2 ≡ f +, n i+1/2 + f -, n i+1/2 ,
and the reconstruction

f +, n i+1/2 := w + 0 2 6 f + i-2 - 7 6 f + i-1 + 11 6 f + i + w + 1 - 1 6 f + i-1 + 5 6 f + i + 2 6 f + i+1 + w + 2 2 6 f + i + 5 6 f + i+1 - 1 6 f + i+2 , (2.30) f -, n i+1/2 := w - 2 - 1 6 f - i-1 + 5 6 f - i + 2 6 f - i+1 + w - 1 2 6 f - i + 5 6 f - i+1 - 1 6 f - i+2 + w - 0 11 6 f - i+1 - 7 6 f - i+2 + 2 6 f - i+3 , (2.31)
where the expressions of the weights w

+/- k k∈ 0,2
and of the fluxes f

+/- i+k k∈ -2,3
can be found in R. [START_REF] Wang | Linear Instability of the Fifth-Order WENO Method[END_REF]. Fig. 2.4 illustrates the output of numerical simulations performed with the WENO5 spatial discretization method coupled with the third-order SSP explicit Runge-Kutta time-stepping method (2.29).

In order to verify the spatial order of accuracy (see Def. 7) of the implementations of the numerical schemes depicted in this section, one can generate a so-called grid convergence error plot. It consists in measuring the mean error of the numerical solution over one box turn of a smooth initial condition projected over equally spaced meshes of different resolutions. Fig. 2.5 illustrates grid convergence error plots for the FOU, LW, MUSCL of [START_REF] Kurganov | New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations[END_REF], and WENO5 schemes. It is worth noting that the WENO5 scheme exhibits third-order accurate behavior for a time step size ∆t = O(∆x). Such a behavior arises because the time stepping method (2.29) is only third-order accurate, and the total error ϵ of the numerical solution then reads

ϵ = O ∆t 3 + ∆x 5 ∆t=O(∆x) = O ∆x 3 + ∆x 5 = O ∆x 3 .
The grid convergence error plot for the WENO5 scheme coupled with the time stepping method (2.29) must be produced using a time step ∆t = O ∆x 5/3 to reveal the spatial order of accuracy of the implementation of the WENO5 method.

MUSCL and WENO schemes employ suitable nonlinear functions of the values of the state variable to provide accurate numerical solutions without spurious oscillations. However, they present applicability limitations. On the one hand, MUSCL uses nonlinear flux limiters formulated upon heuristic arguments, preventing a given limiter from performing well for all problems [START_REF] Waterson | Design principles for bounded higher-order convection schemes -a unified approach[END_REF]. On the other hand, WENO schemes of spatial order of accuracy higher than two may use 3 to 10 times more CPU time than a second-order scheme [START_REF] Shu | High Order Weighted Essentially Nonoscillatory Schemes for Convection Dominated Problems[END_REF]. The recent field of Machine Learning (ML) introduced in the previous chapter relies on nonlinear functions that can be used to address these limitations. Nguyen-Fotiadis et al. ( 2022) employed neural networks as flux limiters in the discretization of the viscous Burgers' equation and showed that they outperform 11 classical flux limiters over a range of mesh resolutions. [START_REF] Bezgin | WENO3-NN: A maximum-order three-point data-driven weighted essentially non-oscillatory scheme[END_REF] used neural networks to enhance a third-order WENO scheme. The resulting data-driven numerical method showed better accuracy levels than the more computationally demanding fifth-order WENO5 scheme at critical points. In this context, the next section investigates the ability of the state-of-the-art data-driven spatial discretization method proposed by [START_REF] Bar-Sinai | Learning data-driven discretizations for partial differential equations[END_REF] and [START_REF] Kochkov | Machine learning-accelerated computational fluid dynamics[END_REF] to provide accurate and stable dynamics.

Data-driven spatial discretizations on 2D cartesian grids

The FOU and LW schemes described in Sec. 2.2.1 employ Taylor series expansions to approximate the spatial derivative in the Convection and Burgers' equations (2.4)-(2.5). The derivative at a given point in space is expressed as a weighted sum of the values of the state variable in a neighborhood of such a point. The set of weights that can approximate the spatial derivative at a given order of accuracy is not unique [START_REF] Fornberg | Generation of finite difference formulas on arbitrarily spaced grids[END_REF]. [START_REF] Bar-Sinai | Learning data-driven discretizations for partial differential equations[END_REF] and [START_REF] Kochkov | Machine learning-accelerated computational fluid dynamics[END_REF] trained neural networks to provide suitable values of the weights for a myriad of problems. More precisely, they aimed to recover features from high-resolution simulations on relatively coarse grids with an appropriate choice of values for the weights. High-resolution simulations of turbulent flows that capture the entire range of eddies' wavelengths are prohibitively expensive. In industrial applications, one performs coarsegrained Large-Eddy Simulations (LES), which account for the effects of unresolved small-scale features using a subgrid-scale model. However, subgridscale models are often formulated using ad hoc assumptions [START_REF] Sagaut | Large eddy simulation for incompressible flows: an introduction. 3rd ed. Scientific computation[END_REF]. [START_REF] Kuramoto | Diffusion-Induced Chaos in Reaction Systems[END_REF][START_REF] Sivashinsky | Nonlinear analysis of hydrodynamic instability in laminar flames-I. Derivation of basic equations[END_REF]. tions. This section studies periodic problems governed by the 2D Euler equations (an inviscid version of the compressible Navier-Stokes equations) on cartesian grids. They can be expressed using Eq. (2.1) with the following choice of state variable U and convective flux F :

U :=      ρ ρu ρv E      , F (U ) :=      ρu ρv ρu 2 + p ρuv ρvu ρv 2 + p u(E + p) v(E + p)      , (2.32)
where ρ, p, and E are the fluid's density, pressure, and energy fields, respectively; and u and v are the components of the flow velocity field along the x and y directions, respectively. A presentation of the considered neural network model that approximates the spatial derivatives of the flux F follows. 2) followed by a non-trainable layer called the Polynomial Accuracy Layer (PAL). The density (ρ), the components of the momentum along the x and y directions (ρu and ρv, respectively), and the energy (E) fields of a given initial condition are used as input features of the CNN. Subsequently, the CNN output features are split among eight channels. Each channel corresponds to a spatial derivative to be approximated (one per scalar flux, see Eq. (2.32)). The PAL projects the features of each channel on the space of weights that approximate spatial derivatives at a given enforced order of accuracy. Finally, the third-order SSP Runge-Kutta time-stepping method (2.29) is used to evolve the initial condition over time.

Neural network model

In order to explain how the PAL works, a 1D field f := f (x) is considered, for the sake of simplicity. The PAL computes a set of weights that is used to approximate the ℓ-th derivative f (ℓ) (x 0 ) of the field f for a given number N of points in the neighborhood of x 0 . The Taylor series of f around the neighbor x n := x 0 + h n (h n > 0) can be expressed as follows:

f (x 0 + h n ) = f (x 0 ) + . . . + f (N -1) (x 0 ) h N -1 n /(N -1)! + O(h N n ).
(2.33) By approximating f (ℓ) (x 0 ) as a sum of the neighboring values of f , weighted by a set of coefficients α := (α i ) i∈ 1,N to be determined, one obtains:

f (ℓ) (x 0 ) ≈ N n=1 α n f (x 0 + h n ) Eq. (2.33) ≈ n α n f (x 0 ) + . . . + n α n f (ℓ) (x 0 ) h ℓ n /ℓ! + . . . + n α n f (N -1) (x 0 ) h N -1 n /(N -1)! + O(α 1 h N 1 ).
(2.34)

By comparing left-and right-hand sides of Eq. (2.34), the following set of constraints arise:

n α n h ℓ n /ℓ! = 1, (2.35) n α n h 0 n /0! = 0, . . . n α n h N -1 n /(N -1)! = 0.
Alternatively, α is solution to the linear system

A f α = b f , (2.36)
with

A f :=     h 0 1 . . . h 0 N . . . . . . . . . h N -1 1 . . . h N -1 N     N ×N , b f := ℓ!     δ 0,ℓ . . . δ N -1,ℓ     N ×1
, where δ n-1,ℓ is the Kronecker's delta, ∀ n ∈ 1, N . The matrix A f is full rank and hence the coefficients that solve the linear system (2.36) are unique. The resulting approximation of

f (ℓ) is O(h N -ℓ 1
), since Eq. (2.35) ensures that

O(α 1 h N 1 ) ∼ O(h N -ℓ 1
).

However, unique solutions prevent using the input of the PAL (i..e., the CNN output) to approximate spatial derivatives. To circumvent this, one must enforce lower accuracy in the approximation of the ℓ-th derivative f (ℓ) , which leads to considering fewer terms in the Taylor expansion (2.33). Equivalently, for a desired order of accuracy O(h m 1 ), with m < N -ℓ, one has to solve the linear system

A α = b,
(2.37)

with

A :=     h 0 1 . . . h 0 N . . . . . . . . . h m+ℓ-1 1 . . . h m+ℓ-1 N     (m+ℓ)×N , b := ℓ!     δ 0,ℓ . . . δ m+ℓ-1,ℓ     (m+ℓ)×1
.

In such a case, the linear system is under-determined and hence admits an infinite number of solutions. Here, any solution α can be written as the sum of an arbitrary fixed solution α p and a vector α null from the nullspace of the matrix

A: A α = A α p + A α null = A α p = b.
In this work, α p is the solution of the fully determined linear system (2.36) (which is solution of the under-determined system (2.37) as well) and α null is the projection of the output of the CNN on the nullspace of A (so that A α null = 0). This projection is performed as follows: let A ≡ U ΣV ⊤ be the singular value decomposition of A, where U and V are square matrices whose columns are the so-called left-and right-singular vectors of A, respectively, and Σ is a diagonal matrix of the so-called singular values of A. It is a well-known fact that the last N -rank(A) columns of V span the nullspace of A. The matrix formed by these columns is denoted N A . The projection of the output x of the CNN on the nullspace of A can be expressed as α null := N A x. Such a projection ensures that PALs return discretization coefficients α that approximate spatial derivatives in cartesian grids with a constrained order of accuracy.

Only first-order derivatives of the fluxes will be approximated using a PAL (i.e., ℓ = 1). For equally spaced grids and centered stencils, Eq. (2.37) then reduces to solving

         1 1 . . . 1 -p -p + 1 . . . p (-p) 2 (-p + 1) 2 . . . p 2 . . . . . . . . . . . . (-p) m (-p + 1) m . . . p m          α =          0 1/h 0 . . . 0          , (2.38)
where h is the cell size and p ∈ N is such that N ≡ 2p + 1 is the number of mesh points of the centered stencil.

A series of experiments with CNN-PAL neural network models that employ Eq. (2.38) was conducted. The parameters common to all the models used in these experiments are depicted in Tbl. 2.1 below. The next section is dedicated to show their outcomes.

Numerical experiments

In this section, CNN-PAL neural network models are trained to predict the dynamics of a given initial condition. It is worth noting that a dimensionless version of the Euler equations is considered. It reads

∂U ′ ∂t ′ + div ′ F ′ (U ′ ) = 0, (2.39)
where

U ′ :=      ρ ′ ρ ′ u ′ ρ ′ v ′ E ′      , F ′ (U ′ ) :=      ρ ′ u ′ ρv ′ ρ ′ u ′2 + p ′ ρ ′ u ′ v ′ ρ ′ v ′ u ′ ρ ′ v ′2 + p ′ u ′ (E ′ + p ′ ) v ′ (E ′ + p ′ )      ,
(2.40) and ρ ′ , p ′ and E ′ are the dimensionless fields of density, pressure and energy, respectively; u ′ and v ′ are the dimensionless components of the flow velocity field along the x and y directions; t ′ is the dimensionless temporal coordinate; div ′ is the dimensionless divergence operator. The dimensionless variables relate to their original counterparts by the following expressions:

ρ ′ := ρ/ρ charact , u ′ := u/u charact , v ′ := v/u charact , E ′ := E/(ρ charact u 2 charact ),
(2.41)

p ′ := p/(ρ charact u 2 charact ), t ′ := t/(l charact /u charact ), div ′ := l charact div,
where ρ charact , u charact , and l charact are called the characteristic density, speed and length, respectively. The characteristic quantities are specific to each initial condition and simulation domain. The training procedure using the original variables has proven to be numerically unstable, unlike when using their dimensionless counterparts. This can be explained by the fact that the scales of the original variables are typically very dissimilar (with E ∼ O(10 5 ) whereas ρ ∼ O(1)), which may affect the convergence of gradient-based optimization methods used in training [START_REF] Wan | Influence of feature scaling on convergence of gradient iterative algorithm[END_REF]. The numerical instabilities are related to the exploding gradients problem in ML [START_REF] Bengio | Learning long-term dependencies with gradient descent is difficult[END_REF]. During backpropagation (see Sec. 1.1.4), the successive multiplication of gradient values with relatively large norm may lead to a value which cannot be represented with finite precision floating point arithmetic. The ML literature encourages performing input normalization to ensure the stability of the training procedure (Shalev-Shwartz and Ben-David, 2014). Moreover, it helps to faster achieve a given accuracy level during training [START_REF] Sola | Importance of input data normalization for the application of neural networks to complex industrial problems[END_REF]. This is because input normalization reduces the distance from the initial state of the trainable parameters of the neural network to the optimal state. The trainable parameters are typically initialized to random values within the (-1, 1) interval. When the input normalization compresses the entire search space into a unitary hypercube, the distance that the backpropagation algorithm (see Sec. 1.1.4) needs to traverse in each iteration is reduced. However, input normalization is a non-physical transformation and the normalized input cannot be used by the PAL. Employing dimensionless variables instead of applying input normalization allows to (partially) address the aforementioned scaling issue in a physical way. Finally, one may notice that the algorithm that implements a numerical scheme for the Euler equations (2.1)-(2.32) can be employed without modifications for the discretization of the dimensionless Euler equations (2.39)-(2.40).

In a first experiment (n.1), a CNN-PAL model is trained to transport a single isentropic vortex by a uniform flow. This case is used to test a numerical scheme's capability to preserve vorticity in unsteady inviscid flows (VI1 Vortex transport by uniform flow | HiOCFD5 2022). The initial condition is set by superimposing vortex-related velocity fields over homogeneous fields. The components u 0 and v 0 of fluid's velocity along the x and y directions, as well as the medium temperature field T 0 are given by

u 0 (x, y) := u ∞ 1 -β y -y 0 R exp(-r(x, y) 2 /2) (2.42) v 0 (x, y) := u ∞ β x -x 0 R exp(-r(x, y) 2 /2) (2.43) T 0 (x, y) := T ∞ - u 2 ∞ β 2 2 C p exp(-r(x, y) 2 ) (2.44) ((x, y) ∈ [0, 1] × [0, 1]),
where u ∞ and T ∞ are the initialized homogeneous values of velocity and temperature (with null component of velocity along the y direction), C p is the fluid's heat capacity at constant pressure, β is the vortex strength, R is the vortex radius, and r(x, y) := (x -x 0 ) 2 + (y -y 0 ) 2 /R is the dimensionless distance to the vortex core position (x 0 , y 0 ). The fluid is considered to be an ideal gas. Given the adiabatic nature of the vortex transport, the initial density (ρ 0 ) and pressure (p 0 ) fields read

ρ 0 (x, y) := ρ ∞ (T 0 (x, y) / T ∞ ) 1 1-γgas (2.45) p 0 (x, y) := ρ 0 (x, y) R gas T 0 (x, y), (2.46)
where R gas = 287.15 J • kg -1 • K -1 is the ideal gas constant and ρ ∞ := p ∞ /(R gas T ∞ ) is the initialized homogeneous value of density, defined for a given initialized pressure p ∞ . The value of R gas relates to dry air at 15 • C (or 288.15 K) and 101, 325 P a, for which the adiabatic index value γ gas = 1.4. Finally, the characteristic quantities of the problem are given by (1928), the Courant-Friedrichs-Lewy condition number N c limits the simulation time step size for a given grid spacing in order to correctly capture wave propagation. The international system of units is used to express all values.

ρ charact := ρ ∞ , u charact := u ∞ ≡ M ∞ c ∞ ,
The dynamics of an isentropic vortex initial condition under the Euler equations corresponds to its convection at the speed u ∞ along the x direction. In order to match the analytical vortex dynamics, the baseline CNN-PAL model depicted in Tbl. 2.1 is trained to minimize the distance between its predictions over time and the exact solution. For a given initial condition U 0 , the loss function reads (see Sec. 1.1.4)

L(U 0 ) := 1 N step Nstep k=1 U k CNN-PAL -U k exact 2 2 , (2.47)
where U k CNN-PAL and U k exact are the predicted and analytical solution states at time step k, respectively. The variable N step denotes the length of the forecasting window which the CNN-PAL model is trained to accurately predict. Since CNNs benefit from spatial translation invariance (see Sec. 1.1.3) and the vortex dynamics simply corresponds to the convection of the initial condition, one can expect that even a relatively small value of N step leads to a trained CNN-PAL model that is capable to accurately perform longrange forecasts. Two independent CNN-PAL models are trained to predict N step = 4 and N step = 16 time steps. The outcomes of the training procedures including error levels after one box turn of the initial condition are summarized in Fig. 2.7. Additionally, the reader can find in Fig. 2.8 the spatial discretization coefficients (i.e., α solution to Eq. (2.38)) predicted by the CNN-PAL model before and after the training procedure with N step = 16. Both sets of coefficients approximate the spatial derivatives of the fluxes along the x direction at the enforced order of accuracy m = 2 (see Tbl. 2.1). However, the ones predicted by the CNN-PAL model after the training procedure lead to a less dissipative evolution of the vortex. Finally, Fig. 2.8 also shows that these coefficients better capture the symmetry of the problem than those of the untrained model.

A second experiment (n.2) consisted of training CNN-PAL models for the convection of a set of vortices (instead of a single one as in experiment n.1). They are trained to predict N step = 16 time steps of the initial condition, which led to better long-range forecasting accuracy levels in experiment n.1. was used. The training datasets are composed of vortices whose strength (β), radius (R) and mach number (M ∞ ) were drawn from uniform distributions, as depicted in Tbl. 2.2. To assess the generalization capabilities of the trained CNN-PAL models, one also tracks the evolution of the mean value of the loss function (2.47) for an indepedent set of vortices, not used to train the models. This set is referred to as the validation dataset, which contains 1024 vortices whose strength (β), radius (R) and mach number (M ∞ ) were drawn from uniform distributions in ranges of values identical to those of the training datasets. Fig. 2.9 shows the evolution of the mean value of the loss function (2.47) evaluated for the validation dataset. One can observe that the difference in the sizes of the training datasets plays an important role in generalization. The CNN-PAL model trained using N s = 128 samples (see Fig. 2.9(b)) better predicts N step = 16 time steps of the previously unseen vortices composing the validation dataset than the CNN-PAL model trained using N s = 32 samples (see Fig. 2.9(a)).

In a last experiment (n.3), CNN-PAL models are trained for a different initial condition. Namely, the (periodic) double shear layer. In this case, the components u 0 and v 0 of fluid's velocity along the x and y directions are given by

u 0 (x, y) := 0.5 u ∞ {tanh [20(y -0.25)] -tanh [20(y -0.75)] -0.2} (2.48) v 0 (x, y) := 0.5 u ∞ sin (4πx) (2.49) ((x, y) ∈ [0, 1] × [0, 1]),
where u ∞ is some preset value of velocity. The fluid's density, pressure, and temperature are initialized to the homogeneous values of ρ ∞ , p ∞ , and T ∞ , shown in Tbl. 2.2. The characteristic quantities of the problem are given by

ρ charact := ρ ∞ , u charact := u ∞ ≡ M ∞ c ∞ , and l charact := 1,
where M ∞ is the Mach number and c ∞ := γ gas R gas T ∞ is the speed of sound in the initialized homogeneous medium.

Unlike the time evolution of the isentropic vortex initial conditions used in experiments n.1 and n.2, that of the double shear layer cannot be derived analytically. In this experiment, reference dynamics for the training of the CNN-PAL models is produced using a WENO5 scheme implementation for the Euler equations, combined with the SSP Runge-Kutta third-order timestepping method (2.29). For systems of conservation laws such as the Euler equations, the WENO flux reconstruction procedure can be performed in one of the following ways: component-or characteristic-wise. The componentwise reconstruction is the one described in Sec. 2.2.3. The characteristicwise reconstruction consists in applying the operations of the component-wise reconstruction to the fluxes projected on a so-called characteristic space. The projection is performed by the averaged flux Jacobian of the system, which is recomputed at each time step. This makes characteristic-wise WENO schemes much more computationally demanding than their component-wise counterparts. Nevertheless, they produce better nonoscillatory results for spatial orders of accuracy higher than 3 [START_REF] Shu | High Order Weighted Essentially Nonoscillatory Schemes for Convection Dominated Problems[END_REF]. For this reason, a characteristic-wise WENO5 scheme is used here to generate reference data for training. The reader is referred to the Algorithm 4.8 of [START_REF] Shu | High Order Weighted Essentially Nonoscillatory Schemes for Convection Dominated Problems[END_REF] for details on the characteristic-wise flux reconstruction procedure. Fig. 2.10 shows grid convergence error plots for the characteristic-wise WENO5 scheme implementation considering the (smooth) isentropic vortex initial condition from experiment n.1.

For a given double shear layer initial condition U 0 , the loss function reads only the density field, which is the one with the lowest error levels:

L(U 0 ) := 1 N step Nstep k=1 U k CNN-PAL -U k WENO5 2 2 , ( 2
L(U 0 ) := 1 N step Nstep k=1 ρ k CNN-PAL -ρ k WENO5 2 2 .
(2.51)

One can expect that training CNN-PAL models using the loss function (2.51) will reduce prediction error levels not only for density, but for all fields that compose U , since the time evolution of each field relies on other fields as well via the Euler equations. This is demonstrated by the results shown in Fig. 2.12.

The three experiments described in the previous paragraphs (summarized in Tbl. 2.2) demonstrate the capabilities of a CNN-PAL model to evolve an initial condition in an accurate fashion. However, the enforcement of stability conditions in the design of CNN-PAL models was not addressed by their original authors, namely [START_REF] Bar-Sinai | Learning data-driven discretizations for partial differential equations[END_REF] and [START_REF] Kochkov | Machine learning-accelerated computational fluid dynamics[END_REF]. The latter claimed that training for the prediction of longer time series improved a model's stability. Fig. 2.13 shows the evolution of the energy field in a long-range forecasting of the double shear layer from experiment n.3. The forecasting is performed by trained CNN-PAL models from all three experiments. One can notice that even for the models trained for the prediction of longer time series, the energy field experience a growth after a certain time, which indicates lack of stability of the trained CNN-PAL models. In the next chapter, a novel spectral analysis that helps to design accurate and stable data-driven numerical schemes is proposed. Chapter 3

Proposing a novel spectral analysis

The previous chapter showed that the lack of stability constraints in formulating data-driven numerical schemes could lead to unstable predicted dynamics. This chapter introduces a novel spectral analysis that helps to generate stable and accurate ML-based numerical methods called the Local Transfer function Analysis (LTA). LTA transforms the problem of developing stable and accurate numerical schemes into a constraint impedance-matching problem in the spectral domain. Furthermore, using LTA, one can describe the complete dynamics of the discrete numerical solution including some specific numerical artifacts. The considered numerical method belongs to the Taylor-Galerkin family of finite-element schemes. They are employed in CERFACS' solver AVBP for large-eddy simulations in an industrial context, which is the ultimate application of the methodology proposed in this thesis. For the linear convection equation, LTA predicts the generation of spurious waves at the boundary between two distinct uniformly spaced domains and the presence of a local linear instability at the junction. For the inviscid Burgers' equation, LTA is applied at a shock front to explain the emergence of Gibbs' phenomenon near the vicinity of the shock. LTA also shows the lack of numerical dissipation at high wavenumbers for the Two-step Taylor-Galerkin C (TTGC) scheme, which necessitates adding a diffusion term for stability.

The Taylor-Galerkin family of schemes

The Taylor-Galerkin (TG) schemes were designed to numerically solve convection problems with high numerical accuracy and low numerical dissipation and dispersion errors. The transport of some arbitrary quantity u = u(x, t) over time and space (in 1D) can be modelled using the simple PDE shown in Eq. (3.1):

u t = -cu x , (3.1)
where c is the convective speed and u(x, 0) = u 0 (x) is the initial solution. This convection operator can also be viewed as a one-sided wave equation (wave moving only to the right). The subscripts (•) t = ∂ ∂t and (•) x = ∂ ∂x denote the partial derivatives with respect to the time and space, respectively. For simplicity, spatial periodicity with a period x ∈ [0, 1] is assumed. The first step of the TG formulation involves the approximation of the time derivative in Eq. (3.1) using truncated Taylor series. Two levels of truncation with leading order of O(∆t 2 ) and O(∆t 4 ) are shown in Eq. (3.2)-(3.3), respectively:

u n t = u n+1 -u n ∆t - ∆t 2 u n tt + O(∆t 2 ) (3.2) u n t = u n+1 -u n ∆t - ∆t 2 u n tt - ∆t 2 6 u n ttt - ∆t 3 24 u n tttt + O(∆t 4 ), (3.3)
where the solution at time t = t n is denoted by u n and the solution at a perturbed time t = t n + ∆t is u n+1 . The temporal derivatives in the truncated Taylor-series expression can be transformed into spatial derivatives using Eq. (3.1). The simplest one-step numerical scheme is obtained from Eq. (3.2) as [START_REF] Roig | One-step Taylor-Galerkin methods for convection-diffusion problems[END_REF])

δu n = -(c∆t)u n x + (c∆t) 2 2 u n xx , (3.4)
where δu n := u n+1 -u n . Eq. (3.4) is similar to the Lax-Wendroff (FV-LW) scheme used in a finite-volume context. Since the spatial derivatives are obtained using a weak-formulation in the Galerkin finite-element method, one obtains an additional mass-matrix on the RHS of Eq. (3.4), which distinguishes it from the FV-LW. More precisely, a linear basis function ψ also called the hat function (plotted in Fig. 3.1) is considered to project the continuous operators to the discrete space of the basis function. It is defined as shown in Eq. (3.5), where h i defines the cell size or the width of the 1D finite-element e i with end nodes i and i + 1. Note that the basis function satisfies the sifting property ψ i (x j ) = δ ij , where δ ij is the Kronecker delta function.

ψ i (x) :=      0 for x < x i-1 or x > x i+1 h -1 i-1 (x -x i-1 ) for x i-1 ≤ x < x i 1 -h -1 i (x -x i ) for x i ≤ x < x i+1 (3.5)
Using this basis allows to express the solution u in terms of its nodal values u j defined at node j : u(x) = N j=1 u j ψ j (x). This basis when introduced into the semi-discrete temporal TG formulation and upon taking the inner product (Galerkin projection) with this same basis yields the complete weak formulation of the TG finite-element method. The weak formulation of the scheme (3.4) is shown in Eq. (3.6), where ⟨a, b⟩ = ab dx defines the inner product of two arbitrary functions a = a(x) and b = b(x):

N j=1 ⟨δu n j ψ j , ψ i ⟩ = -(c∆t) N j=1 ⟨u n j (ψ x ) j , ψ i ⟩ - (c∆t) 2 2 N j=1 ⟨u n j (ψ x ) j , (ψ x ) i ⟩.
(3.6)

In discretization (3.6), the last inner product is obtained through integration by parts, eliminating the need for the basis function to be twice-differentiable. According to Sec. 12.15.5 of T. [START_REF] Sengupta | High Accuracy Computing Methods: Fluid Flows and Wave Phenomena[END_REF], this form is considered to be beneficial when the physical solution admits discontinuities. Furthermore, if the governing equation can be cast in self-adjoint form, the linear algebraic equation which results from the discretization using integration by parts involves symmetric matrix. This leads to significantly reduced computational effort by achieving faster convergence. Discretization (3.6) is hereafter referred to as TG-LW for convenience. Assuming a regular mesh (i.e., h i ≡ h), the integrals in Eq. (3.6) can be evaluated and simplified for the linear basis function as follows:

h -1 M ij a j = N j=1 ⟨a j ψ j , ψ i ⟩ = a i+1 + 4a i + a i-1 6h (3.7) K ij a j = N j=1 ⟨a j ψ x,j , ψ i ⟩ = a i+1 -a i-1 2 (3.8) h -1 D ij a j = N j=1 ⟨a j ψ x,j , ψ x,i ⟩ = a i+1 -2a i + a i-1 h , (3.9)
where a i is some arbitrary nodal value at node i of function a; M , K and D are the mass, stiffness and damping matrices. The weak form of TG-LW can be written in terms of these matrices as

M ij δũ n j = -N c K ij u n j + N 2 c 2 D ij u n j , (3.10)
where N c := h -1 (c∆t) is the Courant-Friedrichs-Lewy (CFL) number, a nondimensional quantity that relates to stability conditions for the numerical scheme as shown in Sec. 3.2. [START_REF] Donea | An analysis of time discretization in the finite element solution of hyperbolic problems[END_REF] developed two-step TG schemes in order to increase time discretization accuracy. It starts by noticing that the 4th-order truncation of the time derivative u t from Eq. (3.3) can be rewritten as where δũ n := ũn -u n and the parameters α = 1 3 and ᾱ = 1 12 produce the so-called fourth-order TTG4A scheme. By choosing ᾱ = 1 9 , one obtains a third-order scheme named TTG3. The Galerkin projection using the linear basis function ψ generates the following weak form for TTG3/4A:

u n t = u n+1 -u n ∆t - ∆t 2 ũn tt + O(∆t 4 ), ( 3
M ij δũ n j = -αN c K ij u n j + ᾱN 2 c D ij u n j (3.14) M ij δu n j = -N c K ij u n j + N 2 c 2 D ij ũn j . (3.15)
Colin and Rudgyard (2000) proposed a six-parameter TG scheme (namely, {α, β, θ 1 , θ 2 , ϵ 1 , ϵ 2 }) that gives third-order accuracy with lower overall dissipation (Eq. (3.16)-(3.17)):

δũ n = α∆tu n t + β∆t 2 u n tt (3.16) δu n = ∆t(θ 1 u n t + θ 2 ũn t ) + ∆t 2 (ϵ 1 u n tt + ϵ 2 ũn tt ).
(3.17)

By performing the Fourier stability analysis assuming a linear finite-element basis, they obtained the optimal values for the six-parameter model. The optimal values reduce the scheme to the single parameter Two-step Taylor-Galerkin C (TTGC)-γ scheme for which

α = 1 2 -γ, β = 1 6 , θ 1 = 0, θ 2 = 1, ϵ 1 = γ, ϵ 2 = 0. (3.18)
Note that γ is a free parameter that controls the dissipation at high wavenumbers. The TTGC-γ scheme is third-order accurate on regular meshes. For the parameter setting α = 1 3 , β = ᾱ, θ 1 = 1, θ 2 = 0, ϵ 1 = 0 and ϵ 2 = 1 2 , the TTG4A (ᾱ = 1 12 ) and TTG3 (ᾱ = 1 9 ) schemes of [START_REF] Donea | An analysis of time discretization in the finite element solution of hyperbolic problems[END_REF] can be derived from this same six-parameter family. The weak form of TTGC-γ can be written in terms of the mass, stiffness and damping matrices from Eq. (3.7)-(3.8)-(3.9) as follows:

M ij δũ n j = -αN c K ij u n j + βN 2 c D ij u n j (3.19) M ij δu n j = -N c K ij ũn j + γN 2 c D ij u n j .
(3.20)

Stability, dissipation and dispersion analysis

Recently, [START_REF] Najafiyazdi | Low-dissipation low-dispersion explicit Taylor-Galerkin schemes from the Runge-Kutta kernels[END_REF] derived a low-dispersion and dissipation TTG scheme starting from the TTGC. They show that using multi-staging one can prevent the evaluation of the third-order spatial derivative that is tedious to obtain for Euler and Navier-Stokes equations. In addition, they provide the connection between the multistage Runge-Kutta (RK) schemes and the TG family of schemes. In fact, TTGC is a sub-class of the Turantype RK methods. An important problem not fully addressed by [START_REF] Rudgyard | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF] is the determination of the optimal γ considering a range of N c values. The authors utilize a single global optimal value for γ for a fixed N c . Using a single global value for γ as in [START_REF] Rudgyard | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF] demands additional artificial dissipation when solving non-linear system with shock or jump discontinuities on irregular meshes [START_REF] Roux | Analysis of unsteady reacting flows and impact of chemistry description in Large Eddy Simulations of side-dump ramjet combustors[END_REF]. In this section, spectral properties of TG schemes in general are depicted and a novel spectral analysis based on local spatio-temporal features is proposed in order to theoretically approach the aforementioned problem for the TTGC scheme in particular.

Global Spectral Analysis (GSA)

Traditional spatial and temporal order of error analysis merely gives the leading truncation error in space and time but does not shed light on the stability and spatio-temporal evolution of the numerical solution. The numerical stability of a discretization of a PDE is classically studied using the von Neumann analysis [START_REF] Charney | Numerical Integration of the Barotropic Vorticity Equation[END_REF]. It is based on the decomposition of the numerical error into Fourier series. The numerical scheme is said to be stable if the errors do not grow unboundedly over time.

Nevertheless, the von Neumann analysis presents limitations: it is applicable to spatially periodic problems only; it implicitly assumes no interactions among the Fourier modes of the numerical error; and it assumes that the numerical error and the numerical solution follow the same dynamics [START_REF] Sagaut | Global spectral analysis: Review of numerical methods[END_REF]. T. [START_REF] Sengupta | Analysis of central and upwind compact schemes[END_REF] addressed these limitations by proposing a novel analysis named the Global Spectral Analysis (GSA). GSA provides information on both the stability and spatio-temporal behavior of the numerical solution. This is achieved by transforming the discretized form and the governing equation into the spectral space using a bi-directional Fourier-Laplace transformation:

u(x, t) = Ω K U (k, ω)e i(kx-ωt) dk dω, (3.21)
where k and ω are the spatial wavenumber and temporal frequency, respectively. The convective problem in Eq. (2.4) transforms to the simplified ex-pression ω = kc in the spectral space. This is the physical dispersion relation that relates the phase speed (c) to the wavenumber (k) and circular frequency (ω). In addition to the individual phase speed of a wave it is also possible to define the velocity associated with a group of waves called the group velocity, v g := dω/dk. For linear dispersive wave propagation, the energy propagates at the group velocity (and not at the phase speed), making it much more important than the phase. Rayleigh (1899) was instrumental in establishing the theory of group velocity, as opposed to phase speed. In addition, for non-dispersive linear propagation, the phase speed and group velocity are the same. Due to the discretization procedure, the numerical solution rather follows the so-called numerical dispersion relation ω N = c N k, where c N is the numerical phase speed (T. K. Sengupta and Dipankar, 2004). Since numerical solutions are almost always dispersive, matching of numerical and exact group velocities becomes critical to accurately capture the energy dynamics of the system. The amplification factor (G) is another important quantity that measures the amount of dissipation and that affects the dispersion the solution undergoes within a time interval of ∆t. More precisely, G affects the value of the numerical phase speed c N , which differs from its physical value c. Since the linear convection in Eq. (2.4) is non-dissipative and non-dispersive, its amplification factor is G phy := e -ikc∆t ≡ cos(kc∆t)-i sin(kc∆t), for which |G phy | = 1.

To obtain the numerical amplification, the existence of a continuous initial condition u(x, t n ) at some time level t n is considered. Let U (k, t n ) be the spatial Fourier transform of the function at t n and U (k, t n +∆t) is the solution at t n + ∆t. The functions can be defined using their Fourier transforms as shown below:

u(x, t n ) := K U (k, t n )e ikx dk (3.22) and u(x, t n + ∆t) := K U (k, t n + ∆t)e ikx dk. (3.23)
The discrete sampling of this function at points x i is performed by sifting with the Dirac delta function δ(x-x i ) transforming the samples into a generalized function representation as shown in Eq. (3.24)-(3.25):

u n i = u(x, t n ) * δ(x -x i ) = X u(x, t n )δ(x -x i )dx = X K U (k, t n )e ikx dk δ(x -x i )dx = K X U (k, t n )e ikx δ(x -x i )dx dk = K U (k, t n )e ikx i dk (3.24) u n+1 i = u(x, t n + ∆t) * δ(x -x i ) = K U (k, t n + ∆t)e ikx i dk (3.25)
Note that X in the integral is the spatial domain of the problem and * is the generalized product [START_REF] Kanwal | Generalized Functions[END_REF][START_REF] Lighthill | An Introduction to Fourier Analysis and Generalised Functions[END_REF]. In general, the limits of the integration is -∞ < k < ∞ but in the case of equally spaced spatial samples (h), the Nyquist sampling theorem states that the numerical scheme can represent wavenumbers only within the range -π < p < π, where p := kh. The numerical amplification G num is then defined as the ratio of the spatial Fourier transforms at time t n and t n + ∆t as shown below (T. [START_REF] Sengupta | Error dynamics: Beyond von Neumann analysis[END_REF]:

G num := U (k, t n + ∆t) U (k, t n ) := U n+1 U n (3.26)
Here the numerical phase c N and group v gn velocities are obtained using the real ℜ(G num ) and imaginary ℑ(G num ) parts of the numerical amplification as follows (T. [START_REF] Sengupta | Error dynamics: Beyond von Neumann analysis[END_REF]:

c N c := ϕ c pN c (3.27) v gn c := 1 N c dϕ c dp , (3.28)
where

ϕ c := -arctan - ℑ(G) ℜ(G) , N c := c∆t h , (3.29)
and

dϕ c dp := 1 |G| 2 ℑ(G) dℜ(G) dp -ℜ(G) dℑ(G) dp .
(3.30)

An ideal numerical scheme for the convection equation must have unit value for the above defined ratios, i.e., satisfies the condition

|G num |/|G phy | = 1, c N /c = 1 and v gn /v g = 1.
Due to the presence of finite dissipation and dispersion errors in numerical schemes, these ratios tend to deviate from this unit value. T. [START_REF] Sengupta | Error dynamics: Beyond von Neumann analysis[END_REF] give an exact dynamics of the solution error ϵ := (u num -u) in the numerical solution u num using the above defined ratios using the error propagation equation shown below:

∂ϵ ∂t + c ∂ϵ ∂x = -c 1 - c N c ∂u num ∂x - v gn -c N k   k 0 ik ′ U 0 |G num | t ∆t e ik ′ (x-c N t) dk ′   dk - ln (|G num |) ∆t U 0 |G phy | t ∆t e ik(x-c N t) dk (3.31)
It is interesting to note that the error follows the same dynamics as the exact equation with extra source terms defined by the dispersion and dissipation errors. The main highlight of GSA is that it identifies that the constant phase in Eq. (2.4) does not remain a constant, rather the numerical phase speed (c N ) is found to be a function of wavenumber (k), an attribute causing dispersion of numerical solution (T. K. [START_REF] Sengupta | A new alternating bi-diagonal compact scheme for non-uniform grids[END_REF]. For a detailed comparison of GSA with other spectral analysis the reader is referred to T. K. Sengupta and Dipankar (2004); T. [START_REF] Sengupta | High Accuracy Computing Methods: Fluid Flows and Wave Phenomena[END_REF].

Application of GSA to Taylor-Galerkin schemes

In order to obtain the numerical amplification, the spatial Fourier transform is applied to the weak numerical discretization (generalized function representation for the samples) in Eq. (3.19)-(3.20) for the TTGC-γ scheme over equally spaced spatial samples and one obtains the following expression

Ũ n = U n -iαN c M -1 K U n + βN 2 c M -1 D U n (3.32) U n+1 = U n -iN c M -1 K Ũ n + γN 2 c M -1 D U n , (3.33) Method Amplification factor (G) Real part ℜ(G) Imaginary part ℑ(G)
Pure convection from Eq. (3.1) cos(pN c ) -sin(pN c )

FV-LW (see Roig ( 2007))

1 + 2 -1 N 2 c D -N c K TG-LW from Eq. (3.4) 1 + 2 -1 N 2 c A D -N c A K TTG3/4A (see Donea et al. ( 1987
)) 1 + 2 -1 N 2 c A D(1 + ᾱN 2 c A D) -N c A K(1 + 6 -1 N 2 c A D) TTGC-γ from Eq. (3.35) 1 + γN 2 c A D -α(N c A K) 2 -N c A K 1 + βN 2 c A D
Table 3.1: Amplification factor (real and imaginary part) for indicated numerical method; where A ≡ M -1 .

where the transformed mass, stiffness and damping matrices are given by

M = 3 -1 [2 + cos(p)] , K = sin(p), D = 2 [cos(p) -1] ,
and p = kh, where k is the spatial wavenumber of the solution in the transformed plane. The amplification factor (G num ) for the TTGC-γ numerical discretization is shown in Eq. (3.34)-(3.35):

G = Ũ n U n = 1 -iαN c M -1 K + βN 2 c M -1 D (3.34) G num = U n+1 U n = 1 -iN c M -1 K G + γN 2 c M -1 D (3.35)
Notice that the numerical amplification G num is a complex function which varies with N c and p, i.e., G num = G num (N c , p). Therefore, finite dissipation and dispersion errors in the numerical solution are almost always found.

To compare the numerical and physical amplification the following ratios are considered: (i) damping ratio (|G num |/|G phy |), (ii) phase velocity ratio (c N /c) and (iii) group velocity ratio (v gn /v g ). A summary of the amplification factor for the various TG numerical schemes is provided in Tbl.

Group velocity preservation

Group velocity preservation is critical compared to preserving the phase speed because, in dispersive systems, the group velocity is the true measure of energy propagation [START_REF] Rayleigh | Scientific Papers[END_REF]. The Group Velocity Preserving (GVP) region (grey shaded area in the v gn /v g contours of Fig. 3.2) is defined as the numerical group velocity v gn within the ±5% tolerance of the physical value. Fig. 3.2(ii) shows the GVP region for various numerical schemes. It is worth noting that TG schemes have much larger GVP regions compared to the second order FV-LW scheme, which makes them more suitable for unsteady computations. This feature is attributed to the additional massmatrix term that bestows a compact stencil to the gradient and diffusion terms.

Filtering effect of the mass-matrix

Two features are common to the Taylor-Galerkin schemes considered here (TG-LW, TTG3, TTG4A and TTGC-γ), namely: (i) they have similar GVP regions and (ii) their GVP regions are much larger than that of FV-LW (for relatively low values of the CFL condition number). This is mainly attributed to the addition of mass-matrix which produces a compact stencil for the gradient evaluation similar to the compact schemes in finite-difference. This is the reason one can find a stark contrast in the GVP region between the FV-LW and other TG schemes. Note that the additional mass-matrix term (A ≡ M -1 ) improves significantly the GVP region but at the cost of reduced stability in TG-LW scheme. To illustrate this, the values of A are plotted in Fig. 3.3 for the Nyquist range 0 ≤ kh ≤ π: the inverse massmatrix term A ≡ M -1 compensates the dissipation and dispersion at high wavenumber components by adding anti-diffusion. The role of A in GVP can be confirmed by noticing that the GVP region of the simple one-step TG-LW is quite similar to that of other two-step schemes for relatively low values of the CFL condition number. This is the primary reason approximations to mass-matrix (by lumping) generate unwanted dispersion errors [START_REF] Guermond | A correction technique for the dispersive effects of mass lumping for transport problems[END_REF]. Petrov-Galerkin (PG) schemes [START_REF] Cardle | A modification of the Petrov-Galerkin method for the transient convection-diffusion equation[END_REF] exploit this filtering property of mass-matrix to introduce upwind basis function that adds stabilization by modifying the matrix as M * = I + 6 -1 (1 -β c )D, where β c is a parameter that controls stabilization. For (i) β c = 0 one recovers the exact mass-matrix (M * = M ), (ii) β c = 1 one obtains a block-diagonal form The dissipation and GVP regions for the TG-LW scheme are plotted in Fig. 3.4 for β c ∈ {0, 0.5, 0.85, 1}. The stability of the TG-LW scheme is greatly improved by higher values of β c . At the same time the operational CFL of the scheme can also be increased significantly by varying the parameter. For convection-diffusion equation, [START_REF] Cardle | A modification of the Petrov-Galerkin method for the transient convection-diffusion equation[END_REF] gives a variety of choices for β c based on the type of temporal discretization. [START_REF] Berzins | Modified mass matrices and positivity preservation for hyperbolic and parabolic PDEs[END_REF] suggests to use a non-linear function for β c that preserves the positivity of the mass-matrix for stabilization and shock capturing. Figure 3.6: Global CFL number versus maximum stable γ = γ max ; stable region obtained from GSA indicated in grey shaded area.

Need for tuning parameters and adaptive nature of TTGCγ scheme

An interesting feature is that the TTGC-γ scheme of [START_REF] Rudgyard | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF] can vary its spectral properties at high wavenumbers by varying γ. The variation of the modulus of the amplification factor of the TTGC-γ scheme by varying γ is shown in Fig. 3.5. This is not possible in TTG3/4A that does not have this tunable parameter. The possibility of redesigning the TTGC-γ such that one varies γ based on the CFL and attenuate spurious waves to some extent is raised. But one must be careful to ensure stability because at higher values of γ the scheme becomes unstable for larger N c . The range of maximum allowable γ = γ max can be deduced from the stability analysis and is shown in Fig. 3.6.

Locally tunable TTGC-γ scheme (TTGC-γ L )

Using GSA, it was shown that parameters with local behavior in the numerical scheme are essential for achieving optimal accuracy. In irregular meshes, even if the convection speed c and time step ∆t are kept constant, the local CFL number N c (= h -1 c∆t) of the finite-element will differ from one element to another due to variation in local cell size (h). Note that non-linear convective equations even on regular meshes (uniformly spaced) give rise to variable convective speeds across elements. Therefore, it is clear that for irregular meshes and non-linear convective systems, choosing the right value of γ becomes critical to the successful application of the TTGC scheme. These results motivate the construction of a locally tunable TTGC-γ scheme called TTGC-γ L . This scheme allows local variations in γ to minimize errors due to numerical dispersion/dissipation contrary to fixing a global γ value in the original TTGC-γ formulation. The weak form of the original TTGC-γ scheme is stated as

N j=1 ⟨δũ n j ψ j , ψ i ⟩ = -α(c∆t) N j=1 ⟨u n j (ψ x ) j , ψ i ⟩ -β(c∆t) 2 N j=1 ⟨u n j (ψ x ) j , (ψ x ) i ⟩ (3.36) N j=1 ⟨δu n j ψ j , ψ i ⟩ = -(c∆t) N j=1 ⟨ũ n j (ψ x ) j , ψ i ⟩ -γ(c∆t) 2 N j=1 ⟨u n j (ψ x ) j , (ψ x ) i ⟩, (3.37)
where α = 2 -1 -γ and β = 6 -1 . To build TTGC-γ L , the coefficients of the original TTGC-γ scheme are redefined using elemental γ values, i.e., γ can vary from one element to another. Given that the support of

ψ i is [x i-1 , x i+1 ],
the sums in Eqs. (3.36)-(3.37) are expressed in the scope of TTGC-γ L as

N j=1 ⟨δu n j ψ j , ψ i ⟩ ←-δu n i-1 ⟨ψ i-1 , ψ i ⟩ + δu n i ⟨ψ i , ψ i ⟩ + δu n i+1 ⟨ψ i+1 , ψ i ⟩ (3.38) α N j=1 ⟨u n j (ψ x ) j , ψ i ⟩ ←-u n i-1 ⟨α e i-1 (ψ x ) i-1 , ψ i ⟩ + u n i ⟨α e i-1 (ψ x ) i | e i-1 , ψ i | e i-1 ⟩ + ⟨α e i (ψ x ) i | e i , ψ i | e i ⟩ + u n i+1 ⟨α e i (ψ x ) i+1 , ψ i ⟩ (3.39) γ N j=1 ⟨u n j (ψ x ) j , (ψ x ) i ⟩ ←-u n i-1 ⟨γ e i-1 (ψ x ) i-1 , (ψ x ) i ⟩ + u n i ⟨γ e i-1 (ψ x ) i | e i-1 , (ψ x ) i | e i-1 ⟩ + ⟨γ e i (ψ x ) i | e i , (ψ x ) i | e i ⟩ + u n i+1 ⟨γ e i (ψ x ) i+1 , (ψ x ) i ⟩, (3.40)
where γ e i is the value of γ associated with the finite-element e i (with end nodes i and i + 1) and α e i := 2 -1 -γ e i . Finite element codes on unstructured grids make use of the so-called element sub-domain paradigm, which consists in splitting any integrals into a sum over the sub-domains or elements [START_REF] Löhner | Applied Computational Fluid Dynamics Techniques: An Introduction Based on Finite Element Methods[END_REF]. The restrictions of the sums (3.38)-(3.39)-(3.40) to an element e i are expressed as

  N j=1 ⟨δu n j ψ j , ψ i ⟩   e i ←-δu n i ⟨ψ i | e i , ψ i | e i ⟩ + δu n i+1 ⟨ψ i+1 , ψ i ⟩ (3.41)   α N j=1 ⟨u n j (ψ x ) j , ψ i ⟩   e i ←-u n i ⟨α e i (ψ x ) i | e i , ψ i | e i ⟩ + u n i+1 ⟨α e i (ψ x ) i+1 , ψ i ⟩ (3.42)   γ N j=1 ⟨u n j (ψ x ) j , (ψ x ) i ⟩   e i ←-u n i ⟨γ e i (ψ x ) i | e i , (ψ x ) i | e i ⟩ + u n i+1 ⟨γ e i (ψ x ) i+1 , (ψ x ) i ⟩. (3.43)
Considering the coordinate x ∈ [0, h i ] for the element e i , one has ψ i (x) := 1 -x/h i and ψ i+1 (x) := x/h i . This results in the following expressions for the inner products in the restrictions (3.41)-(3.42)-(3.43):

⟨ψ i | e i , ψ i | e i ⟩ = h i 0 (1 -x/h i ) 2 dx = h i /3 ⟨ψ i+1 , ψ i ⟩ = h i 0 (1 -x/h i )x/h i dx = h i /6 ⟨α e i (ψ x ) i | e i , ψ i | e i ⟩ = h i 0 (1 -x/h i )(-α e i /h i ) dx = -α e i /2 ⟨α e i (ψ x ) i+1 , ψ i ⟩ = h i 0 (1 -x/h i )(α e i /h i ) dx = α e i /2 ⟨γ e i (ψ x ) i | e i , (ψ x ) i | e i ⟩ = h i 0 (-γ e i /h i )(-1/h i ) dx = γ e i /h i ⟨γ e i (ψ x ) i+1 , (ψ x ) i ⟩ = h i 0 (γ e i /h i )(-1/h i ) dx = -γ e i /h i .
Finally, the TTGC-γ L scheme can be represented in a compact and convenient elemental matrix form as shown in Eqs. where the mass M e , stiffness K e , dissipation D e matrices, and the unknown nodal degree of freedom u e are given by

M e = h i 3 1 1 2 1 2 1 , K e = 1 2 -c c -c c , D e = 1 h i -c 2 c 2 c 2 -c 2 and u e = u i u i+1 .
(3.46) Assembly of the local mass, stiffness and damping is performed elementwise to obtain the global ones. The symbol P e represents a permutation matrix that maps the local degrees of freedom into the correct row entries in the global matrix. Inspired from the edge-based FEM schemes of [START_REF] Löhner | Applied Computational Fluid Dynamics Techniques: An Introduction Based on Finite Element Methods[END_REF]; [START_REF] Luo | Edge-based finite element scheme for the Euler equations[END_REF] and control-volume FEM (CV-FEM) of [START_REF] Baliga | A new finite-element formulation for convection-diffusion problems[END_REF] and that of [START_REF] Bochev | A new Control Volume Finite Element Method for the stable and accurate solution of the drift-diffusion equations on general unstructured grids[END_REF]; [START_REF] Bochev | A multiscale control volume finite element method for advection-diffusion equations[END_REF], a splitting of the integrals into a diagonal term K diag e and a dual control-volume edge-flux contribution K cv e to the stiffness matrix is introduced and yields the following:

K e = -c 0 0 c - 1 2 -c -c c c = K diag e - 1 2 K cv e (3.47)
In the split stiffness term, only the off-diagonal edge contribution is weighted using the element local α e value, i.e., the first step of the TTGC-γ L becomes where α = 2 -1 -γ is a global value as used in the original TTGC-γ of [START_REF] Rudgyard | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF]. Note that this split flux formulation is introduced only for the stiffness term in the first sub-step. The main highlight of the TTGC-γ L scheme is that the stabilization comes directly from the problem of adjusting γ e such that the numerical solution matches the physical problem while maintaining the stability. For γ e > 0, the off-diagonal flux residual is weighted by α e = 2 -1 -γ e so it behaves like a local flux limiter for the convective term on the first step of TTGC-γ L . For the second TTGC-γ L step it provides additional damping (since γ e > 0). By adjusting γ e , unresolved waves can be dissipated while excessive dissipation can be reduced in regions that are fully resolvable by the numerical scheme.

Extension to inviscid Burgers' equation

The inviscid Burgers' equation (3.49) in non-conservation form with periodic boundary conditions is considered:

u t + u u x = 0. (3.49)
Applying the weak formulation using Taylor-Galerkin expansion to Eq. (3.49), one obtains the following elemental matrix form of the TTGC-γ L scheme for the inviscid Burgers' equation similarly to the linear convection as where the mass M e , stiffness K e , dissipation D e matrices and the unknown nodal degree of freedom u e are given by

M e = h i 3 1 1 2 1 2 1 , K e = 1 6 -2u i u i+1 + u i -(u i+1 + u i ) 2u i+1 , (3.52) D e = 1 3h i -u 2 i u 2 i+1 u 2 i -u 2 i+1
and u e = u i u i+1 .

(3.53) P e is a permutation matrix that maps the local degrees of freedom into the correct row entries in the global matrix. By specifying the appropriate initial conditions for the state u(x, t = 0) the solution can be marched forward in time. Similar to the convection discretization, the splitting of the stiffness is introduced in the first TTGC-γ L step into a diagonal term K diag e and a dual control-volume edge-flux contribution K cv e for the stiffness matrix that gives

K e = 1 3 -u i 0 0 u i+1 - 1 6 -u i+1 -u i+1 u i u i = K diag e - 1 2 K cv e .
(3.54)

The splitting and overall form of the final stiffness terms in Eq. (3.47) and Eq. (3.54) are quite similar. The elemental α e weights the off-diagonal term, whereas a constant global optimal value α = 2 -1 -γ weights the maindiagonal terms (similar to Eq. (3.48)). In the next section, the analytical framework to optimize the local γ e values that lead to a stable and accurate numerical scheme is presented.

Local Transfer function Analysis (LTA)

GSA gives the global behavior of the numerical solution vis-a-vis the exact one (dissipation, dispersion, energy propagation speed, etc.). However, the analysis of numerical schemes on non-uniform meshes has been restricted to a subset of meshes whose global structure is known (T. K. [START_REF] Sengupta | Analysis of anisotropy of numerical wave solutions by high accuracy finite difference methods[END_REF][START_REF] Sengupta | A new alternating bi-diagonal compact scheme for non-uniform grids[END_REF][START_REF] Sharma | Hybrid sixth order spatial discretization scheme for non-uniform Cartesian grids[END_REF]. This section considers the numerical convection of a wave packet (WP) with a central wavenumber k and phase speed c = 1 on a mesh that has a jump discontinuity in the mesh spacing at some location x 0 is analyzed. The mesh is equally distributed with a spacing of h L and h R to the left and right of x 0 with a jump in size τ = h R h -1 L = 2 at x 0 . Fig. 3.7 displays the numerical solution at various time steps for the TTGC-γ scheme with a constant γ = 0.01 and central frequency of kh L = 0.8. The location of x 0 is indicated using the blue dot and the yellow dots are the mesh nodes. The frequency of this WP is chosen such that the wave is fully resolved (i.e., inside the GVP region) on the left and slightly above the GVP region on the right. This results in the generation of upstream propagating numerical waves that reflect off the junction x 0 toward the left boundary. Non-linear equations admit shock solutions that produce CFL (N c ) jumps across the shock even on regular meshes and result in a similar behavior. Due to the local nature of this problem, GSA in its current form cannot be used for this analysis.

This motivates the extension of GSA by a new approach which will be named Local Transfer function Analysis (LTA). In LTA, the local numerical amplification G represents the transfer function (of the discrete form) that dictates the dynamics of all degrees of freedom. Therefore, one can conceive these local transfer functions (elemental or nodal degree of freedom) as transfer function blocks forming a circuit block diagram as shown in Fig. 3.8. The impedance or the local transfer function can be obtained as follows. The

:= τ -1 h -1 L c ∆t = τ -1 N L c .
In doing so, the mass, stiffness, and damping matrices for each side (left and right elements) can be obtained, which leads to the following form of TTGC:

M L ũn i -M L u n i + α L (c∆t)K L u n i -β(c∆t) 2 D L u n i = -(M R ũn i -M R u n i + α R (c∆t)K R u n i -β(c∆t) 2 D R u n i ) (3.55) M L u n+1 i -M L u n i + (c∆t)K L ũn i -γ L (c∆t) 2 D L u n i = -(M R u n+1 i -M R u n i + (c∆t)K R ũn i -γ R (c∆t) 2 D R u n i ), (3.56)
where

M L u i = h L 3 u i-1 2 + u i , M R u i = τ h L 3 u i+1 2 + u i , K L u i = - u i-1 + u i 2 , K R u i = u i+1 + u i 2 , D L u i = u i-1 -u i h L and D R u i = u i+1 -u i τ h L . (3.57)
Following the generalized function GSA approach, one can substitute the Fourier transform counterpart of the operators to this left/right TTGC form to obtain

M L ( G -1) + α L N c K L -βN 2 c D L = -M R ( G -1) + α R N c K R -βN 2 c D R (3.58) M L (G -1) + N c K L G -γ L N 2 c D L = -τ M R (G -1) + N c K R G -γ R N 2 c D R , (3.59)
where

M L = 1 3 e -ip 2 + 1 , M R = τ 3 e iτ p 2 + 1 , K L = - e -ip + 1 2 , K R = e iτ p + 1 2 , D L = e -ip -1, D R = e iτ p -1 τ , p = kh L and N c = N L c . (3.60)
Rearranging into a new form, TTGC step 1: impedance of G L contributing from the left element and G R from the right element for a given junction i. Therefore, for maximum transmission of the wave power, one has to satisfy the impedance matching condition, i.e., G L = G † R , where (•) † is the complex conjugate. For minimum reflection, the matching condition becomes G L = G R . Alternatively, one can also match the nodal impedances, i.e., G i+1 = G i or G i+1 = G † i , where the nodal impedance

G = M L -α L N c K L + βN 2 c D L M L + M R + M R -α R N c K R + βN 2 c D R M L + M R G = GL + GR (3.61) TTGC step 2: G = M L -N c K L G + γ L N 2 c D L M L + M R + M R -N c K R G + γ R N 2 c D R M L + M R G = G L + G R (3.
G i := G i | L + G i | R
is simply the sum of the left/right elemental impedances. Therefore, one has the freedom to change the local characteristics of the numerical scheme by impedance matching for (i) minimizing wave reflection or (ii) minimizing the dissipation of the wave energy at the junction. Note that perfect matching by tuning the numerical scheme might not be possible over the entire range of N c and p. Additional damping is usually added (especially near discontinuities) in order to filter such problematic waves (idealized as dash-pots in LTA, as illustrated in Fig. 3.8).

Further analysis of the junction problem

The |G| contours of the nodal impedance G i-1 , G i and G i+1 for the propagation of the WP across a junction of mesh discontinuity (described in the beginning of Sec. 3.2.4 and summarized in Fig. 3.7) is plotted in Fig. 3.9(a-c).

An important feature to notice in Fig. 3.9 is the unstable region (shown in red contour lines) of the numerical amplification G i at the junction. This implies that the wave gains energy while crossing the junction. It is clearly evident in the numerical solution at t = 26∆t and t = 38∆t plotted in Fig. 3.7 that the wave gains amplitude from a value u ≈ -0.5 to a value u ≈ -0.75. The reason this local instability does not develop into a global one is due to the stability of the impedance to the left and right of the junction and that the numerical phase speed at the junction is non-zero. As a result, the wave simply moves away from the junction before it can gain sufficient amplitude.

The spectral spread of the WP in terms of the local Nyquist frequency p is indicated by the vertical striped region. The values N L/R c = {0.4, 0.2} of the local CFL numbers on the left and right side of the junction were deliberately chosen in a region of zero numerical dissipation, i.e., |G| = 1 for the entire spectrum of the WP. The nodal impedance G i-1 is within the DRP region (±5%) but, in the case of the nodal impedance G i+1 , a small fraction of the WP falls outside DRP. This leads to a mismatched impedance that produces reflections. [START_REF] Vichnevetsky | Fourier Analysis of Numerical Approximations of Hyperbolic Equations[END_REF] show that the semi-discrete form of the convection equation is also a consistent discretization for the wave equation. Hence, it admits both a downstream numerical solution called the p-wave (physical) and an upstream numerical solution called the q-wave (spurious). They also mention that at discontinuities of the computational domain, spurious q-waves are generated, for example, at computational boundaries. The junction problem can also be viewed as a discontinuity in the computational domain that produces spurious reflected q-waves propagating upstream. At the computational boundary, [START_REF] Vichnevetsky | Fourier Analysis of Numerical Approximations of Hyperbolic Equations[END_REF] employ time-Fourier transforms to derive the reflection ratio ρ(ω) = q(ω)p -1 (ω) in terms of the q-wave and p-wave solution to the semi-discrete form. Using the LTA, one can obtain the reflection ratio of the junction problem as shown below:

ρ i-1 2 = G i-1 -G i G i-1 + G i and ρ i+ 1 2 = G i -G i+1 G i + G i+1 (3.63)
Unlike the problem at the computational boundary, the interior junction problem requires matching the two nodal impedance to the left and right of the junction to avoid spurious reflection. Note that the reflection formula in Eq. (3.63) is quite general and can be applied at the computational boundary to obtain the boundary reflection ratio similar to [START_REF] Vichnevetsky | Fourier Analysis of Numerical Approximations of Hyperbolic Equations[END_REF]. |ρ| gives the magnitude of the reflected wave in terms of the incident wave at the junction between the two nodal impedance and is plotted in Fig. 3.10 for the range of CFL number 0 < N c < 1 and frequency 0 < p < π/2. For the given spectral spread of the WP the LTA predicts a reflected q-wave with amplitude between 1-5% of the incident wave at the junction, which is confirmed by the numerical simulation (see Fig. 3.7). Note that qwaves (reflections) are present almost for the entire Nyquist spectrum of the incident WP at the junction albeit the amplitudes are very low for kh L < 0.5. For an irregular mesh, several such junctions make the reflected and incident wave dynamics highly complex.

LTA and the spatial spectrogram

It is possible to leverage a spatial analogue of the time-frequency spectral analysis, i.e., space-wavenumber (SW) analysis along with LTA to yield some interesting results. The idea is to construct a spectrogram of the function in the spatial domain and analyze the individual spectrogram snapshots using LTA to infer global error behavior. To extract this localized SW, spatially compact Fourier transform (SFT) can be used. SFT is a linear space-Figure 3.11: Space-wavenumber idealization and local neighborhood of nodal and elemental impedance to the continuous initial condition u(x, t n ) leading to the solution at the next step u(x, t n+1 ).

wavenumber (SW) transform that correlates the spatial signal with a family of waveform that are well concentrated in spatial and wavenumber domain.

The localized SW kernel is chosen such that it is centered around a spatial location ξ yielding the following forward SFT, U (x, k) of a function u(x):

U (x, k) = u(ξ)g(ξ -
x)e -ikx dξ.

(3.64)

Here g(ξ, x) is the symmetric, normalized (i.e., ∥g∥ = 1) and real valued window function that is compactly supported in space. One can recover the inverse SFT as shown in Eq. (3.65):

u(x) = U (x, k)g(x -ξ)e ikx dkdξ (3.65)
and the Parseval's identity for the SFT as shown in Eq. (3.66), respectively:

|u(x)| 2 dx = | U (x, k)| 2 dxdk. (3.66)
Let u(x, t) be a spatio-temporally evolving function which is sampled in time at two time instances t n and t n+1 . Fig. 3.11 shows the two snapshots of the solution u(x, t n ) and u(x, t n+1 ) and an underlying spatial discretization for numerically evolving the solution. One can now make use of the generalized product with the delta function to obtain the numerical amplification in the SW domain (similar to the approach shown in Sec. 3.2.1):

U j (k, t n ) = U (x, k, t n )δ(x -x j )dx = u(ξ, t n )g(ξ -x)e -ikx δ(x -x j )dx dξ = u(ξ, t n )g(ξ -x j )e -ikx j dξ, (3.67)
and the generalized product of the delta function with the inverse SFT yields

u j (t n ) = U (x, k, t n )g(x -ξ)e ikx dkdξ δ(x -x j )dx = U j (k, t n )g(x j -ξ)e ikx j dkdξ. (3.68)
The window function g(x j -ξ) must have compact support and it is useful to make its support equal to the Nyquist limit of the underlying discretization (shaded region in Fig. 3.11). [START_REF] Vichnevetsky | Fourier Analysis of Numerical Approximations of Hyperbolic Equations[END_REF] applied the Parseval's identity to connect the l2 norm of solution error in the physical space to the spectral domain. Motivated by his work, the l2 solution error will be expressed in terms of the local impedance model used in LTA. It is considered the squared l2 norm of the error between the numerical and exact solution at time step n + 1 assuming a known continuous initial solution at time step n. The error is denoted by L2(u n+1 ) and defined below:

Error norm and maximum power impedance matching

L2(u n+1 ) := |ϵ(x, t n+1 )| 2 dx = |u(x, t n+1 ) -ū(x, t n+1 )| 2 dx ≈ j h j+1/2 |ϵ n+1 j | 2 = j h j+1/2 |u n+1 j -ūn+1 j | 2 . (3.69)
The index j denotes the nodal position as illustrated in Fig. 3.11. In LTA, the nodal impedance G j at node j is idealized as the sum of the elemental impedances to the left (G L ) and right (G R ) of the node i.e., G j = G L + G R . Now the error εn+1 in the Fourier space is given by

εn+1 = U n+1 - Ū n+1 ≈ j εn+1 j = j U n j (G j -G phy ). (3.70)
By applying Parseval's identity to the error ϵ one can connect the l2 norm of error in the physical and spectral space as

L2(u n+1 ) ≈ j h j+1/2 |ϵ n+1 j | 2 dx = j |ε n+1 j | 2 dk = j | U n j | 2 (G j -G phy )(G † j -G † phy )dk.
(3.71)

By using the polar form of the complex impedance contributed by the left and right finite-element of node j denoted as

G L = |G L |e iϕ L and G R = |G R |e iϕ R
then the nodal impedance can be expressed using the polar form as

|G j |e iϕ j = |G L |e iϕ L + |G R |e iϕ R .
Note that the physical impedance for the convection problem is G phy = |G phy |e iϕ phy = e -ikc∆t . Now the l2 norm of solution error in Eq. (3.71) can be expressed using this polar form of nodal impedance as

L2(u n+1 ) ≈ j | U n j | 2 |G j | 2 + |G phy | 2 dk -2 j | U n j | 2 |G j ||G phy | cos (ϕ j -ϕ phy ) dk (3.72)
and using the polar form of elemental impedance, the l2 norm of solution error can be expressed using the form below: Otherwise, the minimum solution can lead to an unstable scheme. [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF] reasoned that the linear stability of a scheme does not imply its convergence to discontinuous solutions of nonlinear problems. He introduced a stronger stability criterion, namely the uniform boundedness of the total variation of the numerical solution. This condition implied convergence to the weak solution. The 1-norm total-variation of a discrete quantity u is given by T

L2(u n+1 ) ≈ j | U n j | 2 |G L | 2 + |G R | 2 + |G phy | 2 dk + 2 j | U n j | 2 |G L ||G R | cos (ϕ L -ϕ R ) dk -2 j | U n j | 2 |G L ||G phy | cos (ϕ L -ϕ phy ) dk -2 j | U n j | 2 |G R ||G phy | cos (ϕ R -ϕ phy ) dk. ( 3

Total-variation diminishing (TVD) condition

V (u) = N i=1 |u i+1 -u i |.
Harten defines a numerical scheme E h to be total variation diminishing (TVD) if

T V (E h • u) ≤ T V (u) (3.74)
and showed that numerical schemes with TVD property can successfully capture shocks. Motivated by the work of Harten, a connection between the TVD and the impedance matching condition is sought. It is considered the 2-norm total variation (TV2) defined as

T V 2(u n ) = 1 2 du dx 2 dx = 1 2 lim h→0 u n (x + h) -u n (x) h 2 dx ≈ 1 2 i h i+ 1 2 |ϵ n h | 2 = 1 2 i |u n i+1 -u n i | 2 h i+ 1 2 .
(3.75)

Let εn h be the error in the spectral space for the element formed by the nodes i + 1 and i. Using the fact that translation in x by h becomes simple multiplication by e ikh i+1/2 in the spectral space, one obtains

εn h i+1/2 = U n i+1 e ikh i+1/2 -U n i = U n i e ikh i+1/2 -1 (3.76)
Similarly, one can express εn+1

h as εn+1 h i+1/2 = U n+1 i+1 e ikh i+1/2 -U n+1 i = U n i G i+1 e ikh i+1/2 -G i (3.77)
Using Parseval's identity, one can show that the TV2 from Eq. (3.75) in the spectral space is given by

T V 2(u n ) ≈ 1 2 i |ε n h i+1/2 | 2 dk = 1 2 i | U n i | 2 (e ikh i+1/2 -1)(e -ikh i+1/2 -1)dk = i |u n i | 2 - i | U n i | 2 cos(kh i+ 1 2 )dk.
(3.78)

Similarly, one can obtain the spectral equivalent of T V 2(u n+1 ) as

T V 2(u n+1 ) ≈ i |u n+1 i | 2 - 1 2 i | U n i | 2 G i+1 G † i e ikh i+1/2 + G † i+1 G i e -ikh i+1/2 dk (3.79)
For the solution at time n+1 to be bounded, the T V 2 norm of the solution at n+1 has to be bounded by the exact solution T V 2(u n+1 exact ), i.e., T V 2(u n+1 ) ≤ T V 2(u n+1 exact ). Since the physical amplification for the convection equation is |G phy | = 1, the condition simplifies to

T V 2(u n+1 ) ≤ T V 2(u n ).
(3.80)

For a stable numerical scheme, one can require |G i+1 |, |G i | ≤ 1, which will automatically satisfy T V 2 condition if one imposes the stronger element-wise condition:

| U n i | 2 G i+1 G † i e ikh i+1/2 + G † i+1 G i e -ikh i+1/2 -2 cos(kh i+1/2 ) dk ≤ 0. (3.81)
Firstly, a striking similarity of the T V 2 condition in Eq. (3.80) with the one derived by Harten can be found. Secondly, Eq. (3.81) shows that the T V 2 equality is valid strictly element-wise if the impedance satisfies the matching condition G i+1 = G i . This brings out the one-to-one relationship between the TVD and the impedance matching condition. Satisfying the TVD condition is critical to suppress the spurious waves generated near shock because it leads to minimum reflection impedance matching (not maximum power transmission) as shown by the local transfer function analysis. To achieve minimum reflection, TVD schemes typically sacrifice the conservation of total energy in the wave (power loss), which is what one finds in practice. The 2-norm T V 2 is considered (rather than the 1-norm T V ) because it directly relates to the total energy of the wave, which is conserved in pure convection. One arrives at the same expression by applying the 1-norm TVD condition of Harten point-wise and squaring to obtain the point-wise 2-norm TVD.

LTA of the inviscid Burgers' equation

From the LTA perspective, both linear and non-linear convection are equivalent because linearization or averaging [START_REF] Griewank | Efficient accurate numerical treatment of the modified Burgers' equation[END_REF][START_REF] Ramadan | Numerical treatment for the modified burgers equation[END_REF][START_REF] Ucar | Numerical Solutions of the Modified Burgers' Equation by Finite Difference Methods[END_REF] is used to march the solution from one time step to another. This attribute motivates extending LTA to the non-linear inviscid Burgers' equation shown in Sec. 3.2.3.1. To derive the LTA for non-linear Burgers' equation, it is made the key assumption that the linearization is performed with respect to the time level n (marching solution to n+1) and keeps the same throughout the sub-step in time leading to n+1.

For N c < 1, this is a good approximation, since the changes to the linearization at sub-steps are minimal. A sub-stepwise impedance matching can also be performed but this approach is not shown in the current work. Once the linearization is complete, the convection velocity c i = u n i should be maintained a constant from n to n + 1 as per the assumption. Let τ cR | i := c i+1 /c i and τ cL | i := c i-1 /c i , be two ratios which measure the deviation in velocity between the nodes i + 1, i and i -1, i, respectively. Applying the Fourier-Laplace transform to this linearized system and rearranging, one obtains a two step impedance formula that looks exactly the same as the convection case (shown in Eq. (3.61) and Eq. (3.62)); but the stiffness, damping and CFL number N c are modified as

K L = - τ cL (e -ip + 1) 6 , K R = τ cR (e τ ip + 1) 6 , D L = τ 2 cL e -ip -1 3 , D R = τ 2 cR e iτ p -1 3τ and N c = h -1 L c i ∆t. (3.82) 2 2 2 2 (a) (b) (c) (d) u 1 0 x τcL = 1 τcR = 1 τcL = 1 τcR = 0 τcL = 0 τcR = 1

Analysis at the Burgers' shock

A simplified scenario of a Burgers' initial solution containing a shock spread over three nodes i + 1, i and i -1 is considered. The value of u transitions from u = 1 (before shock) to u = 0 (after shock). Besides, mesh distribution of nodes is assumed to be uniform, i.e., τ = 1. A schematic of the problem is illustrated in Fig. 3.12(a). The LTA dissipation (|G|) and phase angle (ϕ) for the three nodal impedances (linearized) shown in Fig. 3.12(b-d) are plotted in Fig. 3.13. Behind the shock, large variations in phase angle in comparison to aft of the shock take place. This indicates that heavier reflections should occur behind the shock compared to the aft of the shock. Fig. 3.14 shows the numerical solution to a moving shock problem with the given initial conditions for the first two steps using TTGC scheme and one can observe exactly the same behavior as predicted by LTA. [START_REF] Book | The Conception, Gestation, Birth, and Infancy of FCT[END_REF] made similar observations that the phase angle variation was the primary cause for the appearance of over/undershoot at shocks.

|G| i contours at the junction i reveals a locally linear unstable region (|G| > 1) indicated by the grey area in Fig. 3.13. Therefore, the spurious reflections will grow and sustain behind the shock. Comparing the numerical results at t = ∆t and t = 2∆t (see Fig. 3.14), one finds the same behavior as predicted by the LTA. For a standing shock problem, the reflected numerical waves have phase speeds close to zero, resulting in piling up of reflected waves at each time step due to this local instability, which continuously feeds energy leading to numerical blowup. Here, simply adjusting γ is not sufficient to remove this local instability and demands additional external To test the improved stability using artificial dissipation, two variants are considered; where an explicit dissipation term is added (i) at the end of every time step and (ii) at the end of every sub-step of the TTGC-γ scheme.

The elemental form of the explicit dissipation term D * e and its left/right local transfer function D * ,L/R are given by

D * e = ν e ∆t h 2 i -1 1 1 -1 , D * ,L = Pe L (e -ip -1) and D * ,R = Pe R (e ip -1),
(3.84) where ν e is the element local value of the coefficient of the artificial viscosity term and Pe

L/R = ν L/R ∆t h 2
are the non-dimensional Péclet number (S. [START_REF] Sengupta | Global spectral analysis of the Lax-Wendroff-central difference scheme applied to Convection-Diffusion equation[END_REF] of the left and right elements, respectively. Since TTGCγ has two sub-steps, the effective viscosity in variant (ii) will be twice the amount added in the variant (i) case. |G| i contour of the nodal LTA at the shock location i for the variants (i) and (ii) is shown in Fig. 3.15(a)-(b), respectively. A constant Péclet number value of Pe L/R = 0.5 is used for variant (i) and Pe L/R = 0.25 at each sub-step for the variant (ii) case. For the results in the plots, a constant γ value of γ = 0.01 was maintained. Firstly, it can be observed that variant (i) is unstable both in the low and high wavenumber regimes and variant (ii) is unstable only in the low wavenumber regime. In the numerical simulations, as shown in Fig. 3.16, variant (i) produced unstable numerical results. On the contrary, even in the presence of the unstable region at low wavenumber, it was found that variant (ii) was stable even for a stationary shock problem. The numerical stability of variant (ii) becomes clear when one considers the spectral energy cascade of the viscous Burgers' equation shown below [START_REF] Reid | On the transfer of energy in Burgers' model of turbulence[END_REF]:

d U (k) U † (k) dt = d E(k) dt = T (k) -D(k), (3.85)
where U (k) is the complex representation of the Burgers' state u in the wavenumber space, (•) † is the complex conjugate, D(k) = 2νk 2 E(k) represents the rate at which energy is dissipated by wavenumber k, and T (k) is the net rate at which energy is input into the wavenumber k. For problems in 1D, T (k) can be represented using the following form:

T (k) = k * T (k, k * ) = -ikℑ U † (k) U (k * ) U (k -k * ) (3.86)
Note that the term T (k, k * ) is the energy transfer into the wavenumber k due to its interaction with an arbitrary wavenumber k * . A positive value of T (k, k * ) indicates the interactions add energy into k from k * and a negative value indicates the removal of energy from k into k * . For the Burgers' equation, T (k, k * ) is always (i) negative for k * > k, (ii) positive for k * < k and (iii) zero for k * = k [START_REF] Girimaji | Spectrum and energy transfer in steady Burgers turbulence[END_REF]. This means a given scale with wavenumber k always draws energy from larger scales (k * < k) and continuously loses energy to smaller scales (k * > k). Therefore, any numerical instability in scales at the lower wavenumber k (i.e., |G(kh)| > 1) will subsequently be transferred to the scales in the higher wavenumber, which the numerical scheme can dissipate (via the term D(k)). As a result, the energy budget in Eq. (3.85) remains bounded thus providing the stability in variant (ii). For certain critical values of ν, [START_REF] Reid | On the transfer of energy in Burgers' model of turbulence[END_REF] shows that a local or global back-scatter (i.e., reversal of the energy cascade) can occur. By preventing the usual cascade of energy, back-scatter prevents the transfer of energy from lower wavenumbers (k) to higher ones (k * > k); if |G(kh)| > 1 and D(k) is insufficient then this will result in a numerical blowup. In numerical experiments, it was found that for P e > 0.3 the variant (ii) was unstable (see Fig. 3.16(c)) showing such a behavior. The instability at lower kh also makes it impossible to remove over-/undershoots for all shock strengths by simply adding dissipation. Nevertheless, one can manage the instability by bounding E(k) using the dissipation term D(k). [START_REF] Vajjala | Effects of numerical anti-diffusion in closed unsteady flows governed by two-dimensional Navier-Stokes equation[END_REF] applied GSA to the convection-diffusion equation and report a similar route to instability, which they attribute to numerical focusing; a sudden build up of energy in a focused region of the domain. A similar behavior for the Burgers' shock was observed: the build up of errors in lower kh triggers a focused growth of instability near the shock at higher kh. The energy cascade equation is able to clearly explain this numerical focusing or transfer of energy from the unstable low wavenumber region to the higher ones leading to this energy build up.

Using LTA, quantitative and qualitative analyses of numerical schemes on irregular meshes and for non-linear equations were performed. Additionally, LTA suggests that one can optimize numerical schemes either by maximizing power transmission (Sec. 3.2.4.3) or by suppressing spurious waves (Sec. 3.2.4.4). The next chapter is dedicated to exploring the optimization framework resulting from LTA, which can be used to design stable and accurate data-driven numerical schemes.

Chapter 4

Optimizing TTGC-γ L and Machine-Learned TTGC (ML-TTGC)

In the previous chapter, a novel spectral analysis for numerical schemes called the Local Transfer function Analysis (LTA) was introduced. Using LTA, one can formulate a multi-objective optimization framework for designing optimal numerical schemes with stability constraints. This chapter is dedicated to defining and exploring this framework. Firstly, an optimization procedure is used to tune the parameter γ of the original TTGC scheme by introducing element local (optimal) γ values resulting in the new TTGC-γ L scheme (shown in Sec. 3.2.3). However, this procedure is computationally expensive, and its application at every time step is inefficient in practice. A message-passing Graph Neural Network (GNN) model (shown in Sec. 1.1.3) is used as a surrogate model to predict the element local optimal values of γ in the TTGC-γ L . The resulting numerical method is named ML-TTGC. It is used to evolve in time 1D problems governed by the Convection/Burgers' equations. ML-TTGC can diminish dissipation and dispersion errors otherwise produced by the original TTGC scheme. This shows a promising avenue that replaces human-defined rules for numerical methods (often derived for linear systems and on regular grids) with machine-learned discretizations (that will handle non-linear systems and irregular grids).

The optimization problem resulting from LTA

The LTA results derived in Chapter 3 are used to construct a constrained multi-objective optimization framework for numerical schemes. To demonstrate this framework, the considered problem is that of estimating the optimal γ parameter of the TTGC-γ L numerical scheme (defined in Sec. 3.2.3).

The optimization problem is defined as a weighted sum of the two objectives, (i) maximum power output and (ii) minimum reflection as

min γ θ(G L -G † R ) + (1 -θ)(G i+1 -G i ) s.t. |G L | 2 + |G R | 2 ≤ |G phy | 2 (stability condition), (4.1)
where 0 ≤ θ ≤ 1 is the scalarizing weight and θ = 1 implies maximum power while θ = 0 implies minimum reflection. G L/R and G i/i+1 are the local left-/right elemental and nodal impedance in 1D (or an edge in 2D). While LTA enabled to clearly define the design parameters, the correct objective function and constraints for designing numerical schemes, practical enforcement of the impedance matching condition to adjust the dispersion and dissipation characteristics is inconvenient in the spectral space. Since the numerical solution is calculated in the time domain and the objective function is defined in the spectral domain it is convenient to enforce these spectral conditions indirectly in the time domain to optimize the local parameter γ in TTGC-γ L .

Using the results from Sec. 3.2.4.3-3.2.4.4, one can show that the equivalent time domain optimization problem is

min γ θL2(u n+1 ) + (1 -θ) δT V (u n+1 ) + |δT V (u n+1 )| |G L | 2 + |G R | 2 ≤ |G phy | 2 (stability condition), (4.2)
where u n+1 is the numerical solution at time n + 1 obtained from a previous time u n and ūn+1 is a reference or exact solution at time n + 1. The term δT V (u n+1 ) is the difference between the total-variation of the numerical and exact or reference solution, i.e., δT V (u n+1 ) = T V 2(u n+1 )-T V 2(ū n+1 ). Since TVD enforcement is of interest (for minimum reflection impedance matching), i.e., T V 2(u n+1 ) ≤ T V 2(ū n+1 ), only the numerical variations that are higher than the exact or reference solution in the objective are considered. Note that the quantities L2 and T V 2 defined in Eq. (3.75) and Eq. (3.69)

have different mesh scaling of O(h) and 1 O(h) , respectively. Hence, for all the optimization studies the L2 and T V 2 errors are non-dimensionalized using the average element size h avg to bring them to the same order of magnitude, i.e., (h avg ) -1 L2 and h avg T V 2, respectively. The optimization problem should be solved at every time step n to estimate the optimal γ that gives the best numerical solution u n+1 without violating the constraints. Albright and Shashkov (2020) also found that a balance between L1-norm of the error and 1-norm of total variation (i.e., T V measure) is beneficial in their data-driven optimization of artificial viscosity parameters. LTA corroborates theoretically that a similar competing multi-objective measure of the numerical error using L2 and T V is inevitable.

In this work, known exact solutions to the convection problem are used as reference values to obtain the optimal γ. Given the scalar objective in Eq. ( 4.2) and as many γ as the number of finite-elements, the problem is best suited for adjoint-based gradient optimization. By employing the reverse mode of AD (as explained in Sec. 1.2.2), one can obtain exact gradients at a low computational cost, since it is independent of the number of input parameters [START_REF] Griewank | Evaluating Derivatives[END_REF].

These optimal values are then used to build a surrogate model that learns to predict these optimal γ values simply given the solution and discretization. Surrogate models using Machine Learning (ML) have proven to be quite reliable and computationally efficient at representing large dimensional problems [START_REF] Lapointe | Machine learning surrogate models for prediction of point defect vibrational entropy[END_REF]Müller et al., 2019). An important specificity of the data structure here is that unstructured meshes are unsuitable to simple convolutional layers used in many popular neural network architectures, for example, but can instead be represented as a type of graph, as explained in Sec. 1.1.3. Geometric Deep Learning is a dedicated subfield that use graph abstractions to treat relations between entities without prior constraints on their spatial structure, making them most suitable for this task. [START_REF] Battaglia | Relational inductive biases, deep learning, and graph networks[END_REF] recently proposed Graph Neural Network (GNN) models for unstructured mesh simulations to provide physically consistent rollouts of dynamical systems [START_REF] Pfaff | Learning Mesh-Based Simulation with Graph Networks[END_REF][START_REF] Sanchez-Gonzalez | Learning to Simulate Complex Physics with Graph Networks[END_REF]. These GNNs are natural candidates to learn to predict γ directly on the mesh in this study, and are therefore selected employing the representations described later.

A new numerical method named ML-TTGC is proposed, in which the locally tunable TTGC-γ L 's elemental γ values are predicted using the GNN surrogate model. Such a data-driven discretization of the governing equation benefits from improved convergence properties by enforcing stability into the ML architecture using LTA. This is a key differentiation of ML-TTGC from recent works on data-driven discretization [START_REF] Bar-Sinai | Learning data-driven discretizations for partial differential equations[END_REF][START_REF] Kochkov | Machine learning-accelerated computational fluid dynamics[END_REF], where stability cannot be guaranteed or enforced easily. More precisely, the neural network can be constrained to return γ values within the stable regions of TTGC-γ L predicted by the LTA. So for any initial condition, one can enforce stable evolution of the solution. Recent applications of ML models by [START_REF] Stevens | Enhancement of shock-capturing methods via machine learning[END_REF] and [START_REF] Kossaczká | Enhanced fifth order WENO shock-capturing schemes with deep learning[END_REF] to enhance Weighted Essentially Non-Oscillatory (WENO) methods exploited the consistency and stability of the original schemes into the ML model. [START_REF] Melland | Differentiable programming for online training of a neural artificial viscosity function within a staggered grid Lagrangian hydrodynamics scheme[END_REF] replaced the artificial viscosity component of an existing staggered-grid Lagrangian hydrodynamics numerical scheme with a learnable function in the form of an artificial neural network. [START_REF] Magiera | Constraintaware neural networks for Riemann problems[END_REF] introduced a constraint-resolving layer that allows their ML-based surrogate model of a Riemann solver to exactly satisfy the Rankine-Hugoniot condition. An important building block of ML-TTGC is the LTA, which helps identifying the relevant input features for the network. For example, it is known that variations in local CFL number (N c ) and wavenumber spectrum kh of the solution dictate the local optimal value of γ, which in turn changes the impedance. Such a physics-based feature engineering has proven to be useful and effective in the literature. For example, [START_REF] Woo | ETSformer: Exponential Smoothing Transformers for Time-series Forecasting[END_REF] enhanced time-series forecasting performance by extracting dominating seasonal patterns with Fourier transformations.

Numerical implementation and verification of the optimization

This section verifies that the optimization settings defined above leads to a stable evolution of the numerical solution free from spurious waves. The existence of the optimum is critical to the success of the ML surrogate in predicting it. All experiments were conducted in a differentiable solver written in the Julia language [START_REF] Bezanson | Julia: A Fresh Approach to Numerical Computing[END_REF] and the gradients were obtained using the source transformation AD package Zygote (presented in Sec. 1.3.4). Sequential Least-Squares Quadratic Programming (SLSQP) [START_REF] Kraft | Algorithm 733: TOMP-Fortran modules for optimal control calculations[END_REF] gradient-based optimizer in the NLopt [START_REF] Johnson | The NLopt nonlinear-optimization package[END_REF] library was used to solve the constrained multi-objective optimization problem. The optimiza- tion results are compared to the exact solution and to the output of TTGC-γ using the recommended global optimal value of γ = 0.01 by Colin and [START_REF] Rudgyard | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF]. The stability condition |G| < 1 translates into bounded γ values to be set in the optimizer. Non-positive and zero value of γ are unstable, so the lower bound γ min = 0.001. The upper bound of γ is calculated conservatively using GSA (γ max in Fig. 3.6) for the smallest element size that yields a stable result. The optimization procedure is verified for the linear convection of a Wave Packet (WP) with the following form:

u(x, t) = sin [k 0 (x -ct -x 0 )] exp -α(x -ct -x 0 ) 2 (x ∈ [0, 1], t ≥ 0), (4.3)
where k 0 is the central wavenumber of the WP, parameter α controls the spectral spread of the WP and x 0 is the spatial location of the centre of the WP and c is the convection speed.

In the first study (n.1), the TTGC-γ L scheme applied to an equally spaced mesh is optimized. The convection of WP of relatively low wavenumbers (k 0 h < 1.0, inside the GVP region illustrated in Fig. 3.2(e)) on a mesh containing 100 cells (i.e., h = 0.01) is analyzed. Fig. 4.1(a, d, g) show the spatial and spectral content of the considered initial conditions.

The spatio-temporal evolution of the numerical error ε := |u(x, t)u N (x, t)| defined by the absolute difference between the analytical u(x, t) and numerical u N (x, t) solution is plotted for the (i) original TTGC-γ denoted by ε TTGC and (ii) the locally optimized TTGC-γ denoted by ε NLopt as a dedicated subfield of the graphs shown in Fig. 4.1(b,e,h) respectively for indicated kh in the leftmost column. Firstly, one can observe that the optimized scheme has significantly lower error compared to the original scheme thus indicating that an optimal γ exists. Fig. 4.2 plots the order of error convergence with mesh refinement for the two schemes. The optimized scheme follows the same order of convergence of the original scheme but the entire error curve is shifted (reduced) by an order of magnitude (10 -1 ). The thirdorder error convergence of TTGC family of schemes is observed only where the wave packet is inside the GVP region. Outside or near the periphery of the GVP (kh > 1.0) the order of error fails to conform to the third-order behaviour. Qualitatively, the GVP demarcates the smooth or well-resolved region where order of error is valid [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF]. Fig. 4.1(c, f, i) plot the optimal value of γ obtained for the various WP with central wavenumber k 0 h = 0.4, 0.6, 0.8 in the x-t plane. Firstly, one can observe that the optimal γ is highly correlated with the numerical error and hence the numerical solution itself. Moreover, the optimal value of γ is also localized in the time and space and its magnitude increases proportionally with respect to increase in k 0 . These two observations and the error dynamics equation in Eq. (3.31) provide enough evidence that modeling γ based on the spatio-temporal evolution of the numerical solution is feasible.

Multi-objective Pareto optimality

In the second study (n.2), an irregular mesh whose nodes are randomly perturbed is considered. Here, the TTGC-γ scheme generates spurious oscillations even for WP with relatively low wavenumber compared to the Nyquist frequency of the average size (k 0 h = 0.5). Fig. 4.3(a) shows a WP initial condition in low wavenumber regime (k 0 h = 0.5) and its corresponding spectral content. To obtain the irregular mesh node location x * i , perturbations from a uniform random probability distribution (U(0, 1)) toward the positive direction of the x-axis scaled by h 2 are imposed on nodes x i from a uniformly spaced mesh of spacing h, i.e.,

x * i = x i + h 2 U(0, 1) (4.4)
The optimization problem in Eq. (4.2) is solved for different values of parameter θ (the scalarizing weight from Eq. (4.2)) to obtain the Pareto front shown in Fig. 4.3(b). The Pareto front shows that an objective function based purely on the L2-norm of the error is the best choice for evolving waves with low wavenumber content. Note that waves having kh < 1.2 are within the GVP region hence they are well resolved and generate lesser q-waves. Enforcing TVD in this region adds more dissipation and thus cannot provide net improvement to the objective. Fig. 4.3(c-d) compare the numerical solution of the TTGC-γ L (optimized) and the original TTGC-γ at different time instants of the wave convection. For TTGC-γ L , the optimization problem (4.2) is solved for each timestep. The optimized dynamics shows almost negligible spurious oscillations compared to the original TTGC-γ scheme. The optimal values of γ have a high temporal correlation with the numerical solution and exhibit a similar convection pattern. Therefore, it is important to build this correlation information into the surrogate model to improve its predictions.

In the next study (n.3), the convection of a high-frequency WP using a uniform mesh spacing is considered. Here, the original TTGC-γ scheme is . For (n.3) one finds that θ = 0.5 on the Pareto front yields the lowest TVD and L2 error for the wave convection contrary to (n.2) where θ = 1 was found to be better. For kh = 1.2, slightly outside the periphery of the GVP region (where the dissipation as well as the dispersion errors dominate), the optimizer shows that even on regular meshes one can vary γ locally and mitigate dissipation and dispersion errors. For this wavenumber, the original TTGC-γ scheme completely dissipates the wave but the optimized results maintain the wave amplitude as shown by the time evolution in Fig. 4.4(cd). One observes the similar highly correlated convection pattern of the optimal γ values in time for (n.2) as seen in study (n.1).

The final study (n.4) regards the convection of a square wave, which is defined below:

u(x, t) = 1 [x 0 ,x 1 ] (x -ct) (x ∈ [0, 1], t ≥ 0), (4.5)
where x 0 and x 1 are the initial position of the discontinuities of the square wave and 1 is the indicator function defined as

1 A (x) :=    1 if x ∈ A , 0 if x / ∈ A . (4.6)
In gradient-based optimization, discontinuous objective functions are avoided to facilitate smooth descent direction. To ensure a differentiable objective function, a regularized version of the exact solution to the square wave evolution as shown in Eq. (4.7) was used:

u reg (x, t) := {tanh [s(x -ct -x 0 )] -tanh [s(x -ct -x 1 )]} /2, (4.7)
where parameter s is an arbitrary steepening factor, i.e., the higher the value of s steeper is the regularized square wave. In (n.4), a value of s = 150 was used. The Fourier spectrum of the square wave in Fig. 4.5(a) shows that a major proportion of the energy is contained within the GVP region and almost all wavenumbers (broadband) in the Nyquist limit are excited by this waveform. Such a broadband excitation creates spurious waves in the high wavenumber region and even upstream propagating waves, once wavenumbers near the Nyquist limit have negative numerical group velocity (see Fig. L2-norm of the error) yields the best trade-off of L2 error for the lowest TVD error. shows the comparison of the solution evolution of the original TTGC-γ scheme and Pareto optimal one (θ = 1). It can be seen that the latter generates over-/under-shoots near the sharp edges of the square wave of less amplitude than those produced by the TTGC-γ method. For all four studies (n.1-4), one finds the same pattern that the optimal γ is highly correlated temporally and moves with the solution.

In conclusion, by means of the study cases above, the existence of local optimal parameters for the TTGC-γ L scheme could be numerically verified. These optimal values allowed us to address main issues around the original TTGC-γ scheme, namely the generation of spurious oscillations for irregular meshes and the excessive dissipation of waves with relatively high wavenumber content. Despite the encouraging results, this latter case demonstrated that the impedance matching problem cannot always be fully solved (some wave dissipation is still present, see . This is due to two scheme features: stability and resolvable regions. The impedance matching problem is constrained by a stability condition (see Eq. (4.2)). In order to keep the scheme stable, the design space of the γ parameter that could be explored has to be bounded. Along with that, the spectral properties of TTGC-γ indicate resolvable regions (i.e., GVP/DRP zones) that do not significantly change for different values of γ (see Fig. 3.5). This means that high-wavenumber content, out of the resolvable regions for the given scheme, can never be faithfully convected. By means of these optimization studies, not only achievements but also the limits of the proposed approach could be assessed.

Neural network model

The previous section verified that the optimization process on the locally tunable TTGC-γ L scheme leads to improved results on a range of initial conditions for the convection equation. Here, a surrogate model that will provide results of similar quality is sought. For that purpose, an Encode-Process-Decode Graph Neural Network (GNN) architecture inspired by the work of [START_REF] Pfaff | Learning Mesh-Based Simulation with Graph Networks[END_REF] is used. It is composed of a set V of vertices whose relationship is represented using an ensemble of directed edges. These directed edges are represented as a Bipartite graph connecting the sets of so-called sender S and receiver R vertices. S and R participate in a messagepassing scheme that transforms attributes associated with each vertex and Periodic Boundary Figure 4.6: (a) 1D primal/dual mesh (top) and GN based on the dual mesh (bottom); (b) Unstructured primal/dual mesh in 2D (left) and GN based on the dual mesh (right). Arrows denote the GN message-passing, red boxes are the GN vertices (which are equivalent to the dual mesh nodes, indicated by red dots), blue dots are the primal mesh nodes and green dotted lines are the dual mesh edges. edge using neural network layers.

The quantity γ to be predicted is defined element-wise. A natural choice for the set V of vertices would then be the set of finite-elements in the physical mesh, i.e., V := i e i .

For each vertex e i , the following vector of input features (f n i ) is assigned at time step n:

(1D) f n i := [u n i , u n i+1 , γ n-1 i ] ⊤ ,
where u n i and u n i+1 are the solution values at the nodes of the element and γ n-1 i is the cell value of the parameter γ that was predicted at the previous time step. In Sec. 4.2, it was found that the optimal γ are highly correlated temporally (see figs. 4.3(c-d), 4.4(c-d), 4.5(c-d)). Therefore, including γ n-1 i in the input features helps the network to find an optimal γ n i , which should be close to its previous value. By its inclusion, it was observed that convergence of the optimizer to the objective function during the training improved significantly. An edge ϵ ij in this graph corresponds to a pair of a sender vertex e i ∈ S along with a receiver vertex e j ∈ R. Every pair (e i , e j ) of elements sharing a common face is connected by two directed edges, namely ϵ ij and ϵ ji . This is depicted in Fig. 4.6 for 1D and 2D (triangular) meshes. As input edge attribute g ij , the local CFL numbers of the elements e i and e j are considered, i.e., (1D)

g ij := [N c,i , N c,j ] ⊤ .
The above-defined input features are transformed across the Encode-Process-Decode GNN model, as illustrated in Fig. 4.7. This architecture extends Def. 3 of a GNN, which originally only contains message-passing layers. Firstly, an encoding block transports vertices' and edges' features to a higher dimensional latent space (of size L) using independent dense layers. For the applications presented in Sec. 4.4, L = 32. The output of the encoder then feeds a core consisting of message-passing layers (that transforms vertices' attributes based on messages carried by edges). Two message-passing layers are used because the results from the entropy condition (in Appendix D) indicate that second-level neighbours {i -2, i + 2} influence impedance matching problem at a given node i. Finally, a decoding block converts vertex attributes to a feature space compatible with our desired final output, namely the element local γ values. The GNN-predicted γ values feed the TTGC-γ L numerical solver, which outputs an evolved state u n+1 . It is worth noting that the decoder's output is re-scaled using a sigmoid function so that γ n does not violate the stability criteria for the locally tunable TTGC-γ L , which is key to the generalizability of ML-TTGC. During training, the advanced state is compared to the (regularized) exact solution. More precisely, the loss function can be expressed as

L(u 0 ) := 1 N step Nstep k=1 θ L2 u k pred + (1 -θ) δT V u k pred + δT V u k pred , (4.8)
where N step is the number of time steps to predict from a given initial condition u 0 , u k pred is the ML-TTGC output at time k and θ ∈ [0., 1.] is the scalarizing weight from Eq. (4.2). The reader may notice that the GNN model could have been trained using a database of optimal γ generated using NLopt. However, generating such a database for various initial conditions is expensive especially when solving the problem in higher-dimensions.

Figure 4.7: ML-TTGC architecture. The GNN instances trained for the studies presented in Sec. 4.4 use an encoding block, two message-passing layers and a decoding block, which update the vertices and edges attributes using dense layers with a hidden feature space of size L = 32. Each instance has a total of 10, 637 trainable parameters.

Application of ML-TTGC

The previous section introduced ML-TTGC: a framework built upon a GNN model that finely tunes the TTGC-γ L scheme. ML-TTGC will now be used to solve linear and non-linear convection problems. The goal is to assess the proposed architecture's capabilities to provide stable and accurate dynamics for a myriad of initial conditions. More precisely, wiggle-free, less dispersive/dissipative solutions in comparison to those provided by the original TTGC-γ scheme are sought. Tbl. 4.1 below summarizes the set of experiments performed with linear convection and Burgers' equations. 4.1: Sets of initial conditions for ML-TTGC applications.

Problem

Linear convection equation

The evolution of systems following the linear convection equation as defined in Eq. (3.1) is investigated. It is proposed to solve the convection of a WP initial conditions of the form

u(x) = sin [k 0 (x -x 0 )] exp -α(x -x 0 ) 2 (x ∈ [0, 1]), (4.9)
where k 0 is the central frequency, α is the spectral spread and x 0 is the spatial location of the peak of the WP. The GNN γ predictor in ML-TTGC was trained on the forecasting of N step = 16 time steps of a set of 32 samples at wavenumbers k 0 h ∼ U k 0 h (0, 0.6), where U k 0 h is a uniform distribution in the interval (0, 0.6). Here, h denotes the average mesh size of the irregular grids used in training. They are generated according to Eq. ((4.4)). Once only low wavenumber WP are present in the training set, θ = 1 is chosen as impedance matching balance level in the loss function, following findings in Sec. 4.2. Each sample is evolved under a global CFL ∼ U CFL (0.7, 0.9) drawn from the uniform distribution U CFL in the interval (0.7, 0.9), typical to CFD simulations. Fig. 4.8 shows the numerical solution of the ML-TTGC and TTGC-γ scheme for the convection problem of WP at CFL = 0.84, k 0 h = 0.8 (a) and k 0 h = 1.5 (b). Note that these WP are not part of the training data set and were deliberately chosen outside the set to demonstrate robustness. One observes that the amplitude of the spurious oscillations is significantly reduced in ML-TTGC compared to the TTGC-γ scheme where amplitudes of the spurious waves are quite significant. Many flow simulations can exhibit sharp gradients such as these, and such numerical wiggles can lead to unphysical values, notably negative quantities. Therefore, it is critical to mitigate them in the numerical solution. This is most commonly achieved by adding artificial dissipation, which significantly increases the solver runtime.

In this first experiment, the dissipation of high-frequency waves (Fig. 4.8(b)) has proven to be unavoidable, due to the spectral properties of the basis scheme at wavenumbers higher than unity. However, it was found that lowering the CFL values at which the simulations are run can address this to some extent. More precisely, the upper limits for γ under which the scheme is stable are increased as the CFL lowers (following GSA, Fig. 3.6), and hence a larger design space leading to less dissipative solutions can be explored by the neural network. In a second experiment, ML-TTGC was separately trained on 16 samples from the solution sets of WP convection at CFL ∼ U CFL (0.5, 0.7) and k 0 h ∼ U k 0 h (1., 2.). Fig. 4.9 illustrates the high dissipation of TTGC-γ (typical at such range of wavenumbers), which contrasts with ML-TTGC's ability to preserve the wave amplitude for longer time scales. Results shown in Fig. 4.8(b), obtained using an irregular mesh, show higher dissipation of high kh waves; while the results shown in Fig. 4.9(a), obtained using a regular mesh, show less dissipation. The irregularity poses additional challenge as the ML model has to match two competing objective functions, minimum reflection and maximising wave power. Therefore during training the model learns to trade-off dissipation to recover minimum reflection solution (wig-gle free). On a uniform mesh this limitation is non-existent therefore the ML-TTGC scheme learns to be less dissipative.

Inviscid Burgers' equation

For the inviscid Burgers' equation, training database is composed of squarewave initial conditions sampled over irregular meshes:

u(x) = a 1 [x 0 ,x 1 ] (x) (x ∈ [0, 1]), (4.10)
where 1 is the indicator function, as defined in Eq. (4.6), a is the wave amplitude (fixed at 1) and x 0 and x 1 are the initial discontinuities positions that are randomly placed (provided that x 0 < x 1 ). The parameters {x 0 , x 1 } are taken as x 0 = 1 5 (1 + U 0 (0, 1)) and x 1 = 3 25 (5 + U 1 (0, 1)), where U 0 and U 1 are independent random variables following the uniform distribution in the interval (0, 1). Following the method of characteristics, depicted in Fig. D.1(e), their dynamics consist of a progressive rarefaction zone and a moving shock wave that will meet at some point in time. Analytically, the exact solution as a function of time and space can be expressed as

u(x, t) =        (x -x 0 )/t, x ∈ (x 0 , x 0 + at] a, x ∈ [x 0 + at, x 1 + ct] 0, elsewhere (t ∈ (0, t m )), (4.11)
where a is the square wave amplitude, c = a/2 the shock speed and t m is the rarefaction-shock meetup time.

Numerically, the simulation time step for both TTGC-γ and ML-TTGC is defined as ∆t := N global c ∆x min /a, (4.12)

where N global c = 0.84 and ∆x min is the minimum cell size of the considered mesh instance. Results in Fig. 4.10 and Fig. 4.11 compare TTGC-γ and ML-TTGC at the same physical time.

For this particular problem, it was found that the original TTGC-γ produces over-/under-shoots around the shock front and oscillations moving upstream. Using LTA (see Sec. 3.2.5.1), it was shown a mismatch of nodal impedance at the shock with the presence of linear instability was the reason for the over-/under-shoots near the shock. Introducing an artificial dissipation term every sub-step was quite effective in mitigating these numerical issues at the shock. Therefore, it is proposed to add an artificial viscosity term with a locally tunable viscosity parameter into the ML-TTGC scheme.

The viscous terms added after each ML-TTGC sub-step assume the following form: where

D ⋆ e = 1 h 2 i -1 1 1 -1 .
Since the smooth regions of the flow are well evolved by TTGC-γ, it was decided to tune artificial viscosity only near shock, whereas keeping ν = 0 for cells far from it (when at least 10 cells far from the shock front in the upstream direction and 5 cells downstream). A regularized version of the moving shock wave is used as reference, given by u reg, shock (x, t) = a {1 -tanh [s(x -x 1 -ct)]} /2; s = 100.

(4.15)

Due to the relatively limited variability of the training set, it was found that tuning ML-TTGC for the prediction of 20 time steps from 8 initial conditions is enough to give satisfactory results at inference, even at longer time scales, as shown in Fig. 4.10(a). One observes that over-/under-shoots as well as Gibbs' phenomenon are absent by means of the ML-TTGC predicted values for γ and ν. The 'tanh' function was deliberately chosen to regularize the exact solution used in the objective function (Eq. (4.15)) because the exact solution to a moving shock in the viscous Burgers' equation with constant viscosity has a 'tanh' behavior at the shock. It is interesting to notice that the network mimics the viscous Burgers' result by choosing to add constant viscosity (maximum within the set upper bound) in a major part of the shock region. Fig. 4.10(b) illustrates inference outputs on a generalization test. Here, the jump or amplitude of the shock is doubled when compared to the training samples (a = 2). Despite this substantial change to the GNN input features, the resulting dynamics still present significantly fewer undesired oscillations when compared to the original TTGC-γ scheme. Additionally, the GNN model trained only on the square wave was used to predict the dynamics of a sine wave, as shown in Fig. 4.11. The standing shock at the domain boundaries creates a growing numerical instability in the original TTGCγ scheme, leading to numerical blowup (corroborating the LTA findings in Sec. 3.2.5.1). ML-TTGC delivers a stable evolution of the numerical solution, which demonstrates our architectures' excellent extrapolation capability.

A priori knowledge of the shock location allowed to restrict the tuning region for viscosity, which helps the network in finding better optima. In general, one should consider the use of shock sensors to determine strong gradient zones of the flow, where the network will be trained on. Recently, Feng et al. (2020) built shock wave detectors based on machine learning. Alternatively, precise estimates of the local wavenumber could be used to identify shock regions. It was found that the one defined in Appendix C was not accurate enough for the cases studied here, though. Ultimately, one can think of a plug-and-play ML-TTGC architecture composed of some layers trained only around shocks while others specialized in handling smooth regions of the flow. They would be conveniently activated to deliver accurate dynamics in general cases.

The capacities of a GNN architecture used as a surrogate model of the element local parameter γ of the TTGC-γ L scheme were demonstrated. In problems governed by the 1D Convection/Burgers' equations, ML-TTGC outperformed the original TTGC scheme by producing lower levels of dissipation and dispersion errors. In the next chapter, an extension of the ML-TTGC method to higher dimensional problems is proposed. Chapter 5

Extending ML-TTGC to unstructured meshes

Chapter 4 proposed a surrogate model based on graph neural networks that locally tunes a version of the Two-step Taylor Galerkin C (TTGC) scheme used in LES. It was employed in 1D problems governed by the Convection/Burgers' equations. In this chapter, extensions of this machine learning model and training procedure to systems governed by the Convection/Euler equations on higher dimensions are sought. Focus is placed on unstructured irregular grids, typical of industrial CFD. Firstly, the TTGC scheme formulation in higher dimensions is briefly presented. The corresponding ML-TTGC numerical method is then defined and applied in benchmarking problems for numerical schemes on 2D triangular meshes. ML-TTGC's ability to damp spurious oscillations is verified, as observed in the 1D simulations. Finally, an ML-TTGC instance trained only on the convection of vortices following the Euler equations demonstrate capacity to evolve a double shear layer initial condition in a stable fashion, where a numerical blowup otherwise occurs when employing the original TTGC numerical scheme.

TTGC in higher dimensions

By using the divergence theorem, the hyperbolic conservation law (2.1) assumes the following so-called integral form:

∂ ∂t Ω U dV + ∂Ω F • n dS = 0, (5.1)
where Ω is an arbitrary control volume, ∂Ω its boundaries and n is the outward normal to the boundaries. A domain discretization is required to represent U in a finite set of points, and Ω is eventually expressed as the union of polygons (2D) or polyhedrons (3D):

Ω = Ke∈M K e , (5.2)
where M is the mesh object and K e is a cell of this mesh. In what follows, only Triangles (2D) or Tetrahedrons (3D) are considered to be composing a mesh.

In the discretized domain, within a cell K e , Eq. (5.1) reads

∂U e ∂t + R e = 0, (5.3)
where U e is a mean over K e of the state variable and R e (as known as cell residual) is defined as

R e := 1 V e ∂Ke F h • n dS, (5.4)
with F h a numerical approximation of the convective flux F and V e the volume of K e (in 2D, it is the triangle area). Eq. ( 5.3) is solved in different ways following the numerical scheme of choice. On the one hand, the so-called cell-centered schemes store state variables at the center of each cell. On the other hand, cell-vertex schemes store state variables at the mesh nodes. In cell-vertex schemes, one solves a nodal form of Eq. ( 5.3) as follows:

∂U j ∂t + R j = 0, (5.5)
where U j is the value of the state variable at some node j and R j is its associated nodal residual. Fig. 5.1 illustrates the difference between cellcentered and cell-vertex schemes. In particular, the implementation of the TTGC scheme of Colin and Rudgyard (2000) inside the CERFACS solver AVBP is a cell-vertex one. More specifically, it belongs to the class of residual distribution schemes (Deconinck and Ricchiuto, 2017), which compute the nodal residual R j by weighted aggregation of residuals R e from all cells surrounding node j. Such an aggregation is called the scatter operation and is illustrated in Fig. 5.2 below. In order to compute the cell residual R e , it is raised the assumption that convective fluxes vary linearly along faces. In 2D, by considering a triangle cell with nodes {1, 2, 3}, this translates to an approximation of the flux along the triangle face {2 -3} (for instance) as

F h,23 (s) = F 2 + s • F 3 -F 2 ∆s 23 (s ∈ [0, ∆s 23 ]),
where ∆s 23 is the face length. Eq. (5.4) then reads where n ij is the outward normal to the face {i-j} weighted by the face length itself (i.e., n ij := nij ∆s ij , where nij is the outward unit vector normal to the face {i -j}). One can rearrange the flux terms shared among pairs of adjacent faces such that Eq. (5.6) is rewritten as

R e = 1 V e F 2 + F 3 2 • n 23 + F 3 + F 1 2 • n 31 + F 1 + F 2 2 • n 12 , ( 5
R e = 1 2V e [F 1 • (n 31 + n 12 ) + F 2 • (n 23 + n 12 ) + F 3 • (n 31 + n 23 )] . (5.7)
Finally, once the face normals satisfy {i-j} ∈ {{1-2}, {2-3}, {3-1}}

n ij = 0, (5.8)
one is left with the final formula below for cell residuals in 2D triangular meshes:

R e = - 1 2V e k∈Ke, i, j̸ =k F k • n ij .
(5.9)

A relation similar to Eq. (5.9) can be found for 3D tetrahedral meshes. In practice, one works with the more generic formula:

R e = - 1 d V e k∈Ke F k • S k , (5.10)
where S k is the so-called nodal normal defined as

S k := face ∋ k - d n face vert n face ,
with n face vert as the number of vertices in a face (equals 2 for triangles and 3 for tetrahedrons), n face as the outward normal to the face weighted either by its length (in 2D) or surface (in 3D) and d as the number of spatial dimensions (2 in 2D, 3 in 3D).

The nodal residual R j is computed by scattering residuals from all cells which the given node j belongs to: (5.11) where D j is the set of this node's neighboring cells, V j is a control volume around the node defined as

R j := 1 V j e ∈ D j D j,e V e R e ,
V j := e ∈ D j V e n elem vert , ( 5.12)
with n elem vert as the number of nodes in a cell and D j,e is a distribution matrix following the numerical scheme of choice:

• Central Differences (CD) -It is the simplest numerical scheme one can obtain from Eq. (5.11) by choosing

D CD j,e := 1 n elem vert Id,
where Id is the identity matrix. It is unconditionally unstable when used with an Euler explicit first order time stepping scheme.

• Lax-Wendroff (LW) -It is a time-space combined discretization scheme:

forward in time by using first order Euler explicit scheme, and centered in space by using second order central differences (along with an additional diffusive term to ensure stability). Its distribution matrix is expressed as (5.13) where A e is the Jacobian matrix of convective fluxes, averaged across each element.

D LW j,e := 1 n elem vert Id - ∆t 2 d V e A e • S j,e ,
Due to the Galerkin projection, Eq. (5.5) assumes the two-step discrete form below for the TTGC scheme shown in Eq. (3.16)-(3.17)-(3.18) [START_REF] Lamarque | Schémas numériques et conditions limites pour la simulation aux grandes échelles de la combustion diphasique dans les foyers d'hélicoptère[END_REF]:

Ω δ U n ψ j dV = -(0.5 -γ) ∆t L j (U n ) + β ∆t 2 L j (U n ) (5.14) Ω δU n+1 ψ j dV = -∆t L j ( U n ) + γ ∆t 2 L j (U n ), (5.15)
where δ U n := U n -U n and δU n+1 := U n+1 -U n , ψ j is a hat function (as known as P1 element in finite-element terminology) and the operators L j and L j are parts of the nodal residual that can be expressed in terms of the distribution matrices D CD j,e and D LW j,e as follows: (5.17) or, more explicitely,

L j ≡ e∈D j D CD j,e V e R e (5.16) L j ≡ 2 ∆t e∈D j (D CD j,e -D LW j,e ) V e R e ,
L j ≡ 1 n elem vert e∈D j V e R e
(5.18) (5.19) with S j,e as the nodal normal S j with respect to the element e. Finally, the left hand sides of Eq. (5.14)-(5.15) are expressed with the help of a mass matrix M : M δ U n and M δU n+1 , respectively. Let Φ be the resulting nodal residual vector for step 1 of TTGC. The linear system to be solved is

L j ≡ 1 d e∈D j (A e R e ) • S j,e ,
M δ U n = Φ.
(5.20)

The Jacobi iterative method is applied, and is given by [START_REF] Lamarque | Schémas numériques et conditions limites pour la simulation aux grandes échelles de la combustion diphasique dans les foyers d'hélicoptère[END_REF] δ U

(p) = δ U (0) -V -1 (M -V) δ U (p-1)
(5.21) at iteration step p (with V := diag(V j ) and δ U (0) := V -1 Φ).

The next section explains how the numerical method is adapted to accept per-cell values of the parameter γ as done in Sec. 3.2.3 for the 1D case.

ML-TTGC for 2D Convection/Euler

The 2D Convection and Euler equations have convective fluxes respectively given by (see Eq. (2.1))

F Convection (U ) := c X U c Y U , (5.22) F Euler (U ) :=      ρu ρv ρu 2 + p ρuv ρvu ρv 2 + p u(E + p) v(E + p)      , with U :=      ρ ρu ρv E      , (5.23)
where c X and c Y are the components of the convection velocity field along x and y directions, and ρ, u, v, E and p are the scalar fields of density, xvelocity, y-velocity, (total) energy and pressure, respectively. The variables E and p are related by

E = p γ gas -1 + 1 2 ρ(u 2 + v 2 ), (5.24)
where γ gas is the adiabatic index of the considered fluid.

Locally tunable versions of the TTGC scheme exposed in the previous section must be derived (as done in Sec. 3.2.3). Eq. (5.14)-(5.15) are reformulated as include per-cell values of the scheme parameter γ. For the Convection case, these operators can be defined as follows :

Ω δ U n ψ j dV = -∆t L local j (U n ) + β ∆t 2 L j (U n ) (5.25) Ω δU n+1 ψ j dV = -∆t L j ( U n ) + ∆t 2 L local j (U n ), ( 5
L local j ≡ 1 3 e∈D j [(0.5 -γ X ) (V e R e ) X + (0.5 -γ Y ) (V e R e ) Y ] (5.27) L local j ≡ 1 2 e∈D j A Convection e • S j,e V e [γ X (V e R e ) X + γ Y (V e R e ) Y ] , (5.28)
where γ X and γ Y are components of the parameter γ along x and y directions, respectively, and (V e R e ) X := -

1 2 k∈Ke F X k S X k and (V e R e ) Y := - 1 2 k∈Ke F Y k S Y k .
One may notice that (V e R e ) X and (V e R e ) Y result from a splitting of the cell residuals computation that follows Eq. (5.10):

V e R e := - 1 2 k∈Ke F k • S k = - 1 2 k∈Ke F X k S X k + F Y k S Y k ≡ (V e R e ) X + (V e R e ) Y .
(5.29)

For the Euler case,

L local j ≡ 1 3 e∈D j [(0.5 -γ) V e R e ]
(5.30)

L local j ≡ 1 2 e∈D j (γ A Euler e R e ) • S j,e .
(5.31)

Here, a scalar value for the parameter γ is kept for scheme consistency.

Given the locally tunable versions of the TTGC scheme defined above, a Graph Net (GN) model identical to the one described in Sec. 4.3 is coupled to the solver as a surrogate model that provides optimal values of the parameter γ, except that it expects input features of nature and dimensionality specific to the studied problem. The input features of graph vertex (i.e., mesh element) e i at each time step n are (2D Convection)

f n i := U n i1 , U n i2 , U n i3 , γ n-1 i,X , γ n-1 i,Y ⊤ , (5.32) (2D Euler) f n i := [ρ n i1 , ρ n i2 , ρ n i3 , (ρu) n i1 , (ρu) n i2 , (ρu) n i3 , (ρv) n i1 , (ρv) n i2 , (ρv) n i3 , E n i1 , E n i2 , E n i3 , γ n-1 i ] ⊤ .
(5.33)

In other words, γ values from the previous time step and the values of conservative variables at the three nodes composing the mesh cell (a triangle) are provided as information at vertices of the GN. Finally, edges exchanging messages between graph vertices (mesh cells) carry their per-cell defined CFL numbers (N c ), such that the input edge attribute g ij of an edge sending a message from cell i to cell j is

g ij := [N c,i , N c,j ] ⊤ . (5.34)
Training loss function is also similarly defined as a balance between the squared l2-norm of the error over time with respect to the available analytical solution and some Total Variation (TV) error measure specifically set to 2D problems based on the mesh-edge TV. It can be expressed as

L(U 0 ) := 1 N step Nstep k=1 θ L2 U k pred + (1 -θ) δT V U k pred + δT V U k pred , (5.35)
where N step is the number of time steps to predict from a given initial condition U 0 , U k pred is the ML-TTGC output at time k, θ ∈ [0, 1] is some scalarizing weight,

L2(U ) := i U i -U i 2 2
(5.36) with U the reference analytical solution, and δT

V (U n+1 ) := T V 2(U n+1 ) - T V 2(U n+1 ) such that T V 2(U ) := 1 2 i,j ∥U i -U j ∥ 2 2
(5.37) for all pairs {i, j} of adjacent mesh nodes. The reader may notice that the definitions used here for L2 and T V 2 (coming from Eq. (5.36) and Eq. (5.37)) reduce to those from the 1D studies in Chapter 3 (Eq. (3.69) and Eq. (3.75), respectively), apart from the mesh geometry dependency. This is done to avoid related expensive computations during training and is assumed to have little influence on final results once only the values of the solution field U are subject to optimization, not the mesh geometry itself.

In what follows, ML-TTGC training/inference is performed on irregular triangular grids generated from perturbing (interior) regular mesh nodes positions. The node is placed somewhere within a circle around its original position, i.e., ∆x i = R i cos ϕ i and ∆y i = R i sin ϕ i are the displacements imposed to the x and y coordinates of a given node i, respectively. The circle's radius R i is drawn from the uniform distribution U(0, 0.35 h), where h is the (minimal) regular grid spacing, and ϕ i is a random angle drawn from the uniform distribution U(0, 2π). Fig. 5.3 illustrates irregular mesh generation for ML-TTGC evaluation. For the Convection case, the training set was composed of wave packets of the following form : (5.38) where k 0,X and k 0,Y are the components along the x and y directions of some central wavenumber, α is the spectral spread and (x 0 , y 0 ) is the spatial location of the peak of the WP. Training is performed on conditions similar to those of Sec. 4.4.1, with forecasting of N step = 16 time steps of a set of 32 samples at wavenumbers k 0,X h = k 0,Y h drawn from the uniform distribution U(0, 0.8) and (x 0 , y 0 ) = (0.5, 0.5) for all samples. The convection velocity (c X , c Y ) = (1., 1.) is also a constant for the training set. Fig. 5.4 compares inference outputs for one of the training samples from network models trained either on pure l2 loss function (i.e., θ = 1) or on pure TVD one (i.e., θ = 0) against the original scheme and the exact solution. One may notice a significant reduction in the spurious oscillations amplitude as a result of the training procedure, with the pure TVD trained model giving slightly noisier dynamics. The fact that the TTGC scheme is linearly unstable for the junction problem analyzed under the light of LTA in Sec. 3.2.4.1 might explain the reminiscent wiggles on ML-TTGC outputs: wave reflections originated from mesh irregularity will not be fully damped out. For the Euler case, isentropic vortices of varying strengths and radii are chosen to compose the training set1 . It consists of 32 samples with constant vortex core position (x 0 , y 0 ) = (0.5, 0.5), and initialized homogeneous temperature T ∞ = 300 K and pressure p ∞ = 100, 000 P a. The Mach number M ∞ := u ∞ / γ gas R gas T ∞ , the vortex strength β and the vortex radius R are drawn from the uniform distributions U(0.1, 0.6), U(0.1, 0.8) and U(0.05, 0.15), respectively. Values of R gas and γ gas are assumed to remain constant in the considered range of pressure and temperature values. The machine learning model is trained to forecast N step = 16 time steps of the vortex convection, for which the analytical dynamics is given by While effective, these methods require expensive gradient propagation through the adjoint solver for each new configuration. To alleviate this, a (finite) database of cases where the analytical solution is known is generated and a surrogate model is designed by training a Graph Neural Network (GNN) to predict the local optimal parameters for impedance matching. The data-driven numerical method that arises is termed ML-TTGC. A novel aspect of this architecture is that stability is inbuilt into the ML model through the LTA. ML models being approximate surrogates can lead to predictions in the numerically unstable region (for cases outside the training set), which motivated to constrain their output such that ML-TTGC always provides stable dynamics. Applications to the 1D Convection and Burgers' equations demonstrate both the interpolation and extrapolation of the matched impedance values for cases outside the training set.

U 0 (x, y) = sin [k 0,X (x -x 0 ) + k 0,Y (y -y 0 )] exp -α (x -x 0 ) 2 + (y -y 0 ) 2 ((x, y) ∈ [0, 1] × [0, 1]),
U (x, y, t) ≡ U 0 (x -u ∞ t, y) (t > 0). ( 5
Nevertheless, the reader might have noticed that validation sets as defined in Sec. 1.1.4 are absent from ML-TTGC trainings performed in Sec. 4.4 and in Sec. 5.2. This is on purpose, once the capacity to fully solve a restricted set of problems need to be assessed before evaluating generalization features of the coupled ML-PDE system. For instance, it was reported in Sec. 4.4.2 that ML-TTGC without any artificial dissipation term in the numerical scheme side was unable to solve the shock propagation proper to the inviscid Burgers' equation, even if trained on a single sample of a square wave initial condition. Similarly, as shown in Fig. 5.4, tuning of the γ parameter of the TTGC scheme in the convection of 2D wave packets was unable to fully damp all spurious oscillations due to mesh irregularities. Finally, Fig. 5.5 shows that vortex isentropic convection following ML-TTGC still present dispersion errors propagating through the wake of the vortex.

It becomes then clear that the choice of numerical scheme to be tuned by the coupled machine learning model is critical to the success of the proposed approach. Besides, from the experiments reported in Appendix B, it is also clear that an efficient differentiation tool that supports array mutation is required for reasonable performance of the training procedure, both runtime and memory-wise. Despite being able to handle mutating code from a highlevel programming language (Julia), the Enzyme library (Sec. 1.3.5) does not offer the same maturity as older tools like Tapenade (Sec. 1.3.3) and may return unexpected gradient values. Future works will focus on the use of Tapenade to develop data-driven discretizations in the context of highperformance CFD, such that the resulting trained ML models can be easily integrated to mature CFD codes like AVBP for industrial applications.

Recently, large-scale applications of ML models in the field of weather forecasting have produced exciting results. [START_REF] Bi | Accurate medium-range global weather forecasting with 3D neural networks[END_REF] proposed Pangu-Weather: an ML model that can forecast global weekly weather patterns at a significantly faster pace than conventional forecasting methods, while maintaining comparable accuracy. [START_REF] Lam | GraphCast: Learning skillful medium-range global weather forecasting[END_REF] employed GNNs to generate a 10-day forecast in under 60 seconds. Their GNN model was named Graph-Cast. It outperformed the European Centre for Medium-Range Weather Forecasts' operational method on 90 percent of the 2760 variables. Graph-Cast surpassed Pangu-Weather on 99 percent of the 252 targets. Despite the results achieved by these models, they present risks for being exclusively data-driven (Ebert-Uphoff and Hilburn, 2023). One of these risks concerns extreme events, which have an increased likelihood in a changing climate. Besides, these events can lead to highly erratic forecasts, since the performance of fully data-driven ML models usually degrades under conditions not seen during training. Finally, these models do not consider dependencies between forecasting variables, as opposed to the traditional frameworks.

In order to circumvent the aforementioned issues, this work advocates using traditional numerical methods as a basis for the development of ML models in weather forecasting and CFD, in particular, and in physical sciences in general. These hybrid ML models can also produce high-accuracy results with the advantage of being interpretable [START_REF] Murdoch | Definitions, methods, and applications in interpretable machine learning[END_REF]. In regulated industries such as aerospace, compliance with safety standards and regulations is crucial. Interpretable ML models can facilitate compliance by providing explanations and justifications for the predictions provided, making it easier to demonstrate adherence to safety guidelines. The cited works of [START_REF] Bar-Sinai | Learning data-driven discretizations for partial differential equations[END_REF] and [START_REF] Kochkov | Machine learning-accelerated computational fluid dynamics[END_REF] on data-driven spatial discretizations of enforced order of accuracy, [START_REF] Stevens | Enhancement of shock-capturing methods via machine learning[END_REF] and [START_REF] Kossaczká | Enhanced fifth order WENO shock-capturing schemes with deep learning[END_REF] on data-driven perturbations of the parameters of WENO methods, and Bezgin, Schmidt, et al. (2021) on data-driven flux reconstruction for a finite-volume scheme are some examples from the state-of-the-art literature showing the potential of interpretable hybrid ML models in solving PDEs. Still recently, Kossaczká et al. (2023) employed neural networks to approximate the spatial truncation error in finite difference schemes. The distinguishing building block of the work in this thesis is employing spectral analysis on both time and space discretizations to better define stability constraints for the hybrid ML model, and also to interpret the optimization of the model as an impedance-matching problem to be solved at each time step. to build high-level APIs like the NumPy library [START_REF] Harris | Array programming with NumPy[END_REF] on top of C/C++ procedures that execute fast linear algebra computations.

In the realm of machine learning, the same reason led big companies like Google and Meta to launch Python APIs for TensorFlow, JAX and PyTorch. They can enjoy from the fact that most neural network architectures used in image recognition, natural language and time series processing boil down to a relatively small subset of linear algebra kernels whose forward/backward passes can be highly optimized [START_REF] Isaksson | Four Common Types of Neural Network Layers[END_REF]). Nevertheless, one common limitation faced by these popular frameworks is the lack of full support to array mutation, as illustrated by the code snippets below:

TensorFlow :

>>> import tensorflow as tf >>> tf. __ version __ '2.9. File "<stdin>", line 1, in <module> RuntimeError: a view of a leaf Variable that requires grad is being used in an in-place operation. pure functions. " (From JAX -The Sharp Bits -JAX documentation (2022))

An excessive number of reassignments in workflows typical to numerical schemes for unsteady problems (where, for instance, the enforcement of boundary conditions would require array mutation at each time step) eventually leads to high memory pressure, which might be one of the main bottlenecks for the development of differentiable PDE solvers. The literature contains works pointing to memory consumption issues for adjoint codes. In their 2021 paper, [START_REF] Kochkov | Machine learning-accelerated computational fluid dynamics[END_REF] reported using gradient checkpointing2 at each time step of a JAX-based CFD solver3 to avoid prohibitive memory requirements during training. [START_REF] Cardesa | Adjoint computations by algorithmic differentiation of a parallel solver for time-dependent PDEs[END_REF] also used checkpointing to alleviate memory footprint of the adjoint of the CFD solver JAGUAR [START_REF] Cassagne | JAGUAR: a New CFD Code Dedicated to Massively Parallel High-Order LES Computations on Complex Geometry[END_REF]. However, while [START_REF] Cardesa | Adjoint computations by algorithmic differentiation of a parallel solver for time-dependent PDEs[END_REF] could differentiate very long time-stepping sequences (ca. 10 6 time steps) in double precision arithmetics, [START_REF] Kochkov | Machine learning-accelerated computational fluid dynamics[END_REF] were limited to only 32 time steps forecasting during training of their coupled ML-PDE system, still using single precision arithmetics. This is mainly due to the native inability of JAX to handle mutating code. More recently, [START_REF] Bezgin | JAX-Fluids: A fully-differentiable high-order computational fluid dynamics solver for compressible two-phase flows[END_REF] claimed that the largest limitation of their differentiable solver in JAX4 was the available storage space of the GPU, which was an NVIDIA RTX A6000 card with 48 GB of memory, the largest of its series (List of Nvidia graphics processing units 2022).

Sec. B.2-B.3 illustrate the mutation problem with applications of the reverse mode of AD. In each case, a given quantity of interest is computed by two programs: one that employs array mutation whereas the other just performs reassignments. Both programs are written in the Julia programming language [START_REF] Bezanson | Julia: A Fresh Approach to Numerical Computing[END_REF]. Finally, the performances of their adjoints are benchmarked and compared.

B.2 Analysis of a discretization operator

The first problem consists in building a discretization operator of a scalar quantity u over a one dimensional uniform mesh. By means of Taylor expan- where ∆x is the mesh spacing, x i (i ∈ 1, N ) is some of the N mesh nodes and u i := u(x = x i ). Programs that compute the squared l2-norm of the operator du(x)| x=x i := u i-2 -4u i-1 + 3u i (B.1) will be crafted. For the sake of simplicity, the constant 2∆x was ignored.

Here, periodic conditions at the boundaries are considered, so that du(x 1 ) := u N -1 -4u N + 3u 1 , (B.2) du(x 2 ) := u N -4u 1 + 3u 2 .

(B.3) The rolling operator circshift is then applied to the input variable v 0 (i.e., u) to create the arrays u_m1 and u_m2, which will be assigned to the intermediate variables v 1 and v 2 : with the version 1.8.5 of the Julia programming language and the package BenchmarkTools.jl [START_REF] Chen | Robust benchmarking in noisy environments[END_REF]. Fig. B.5(a) shows that the adjoint program du_reassignment_adjoint allocates approximately twice the amount of memory allocated by its primal du_reassignment. Meanwhile, the mutating implementations du_mutation_adjoint and du_mutation exhibit similar memory footprints. Runtime-wise, the mutating implementations present worse performance because their loops are not parallelized.

v 1 ← circshift(v 0 ,
An extension of the problem to higher dimensions is analyzed, motivated by the fact that LES codes handle turbulence, which is inherently a 3D phenomenon. Adjoints of non-mutating and mutating programs that compute Memory-wise, the results are equivalent to those found in the 1D experiments. Runtime-wise, the mutating implementations prove to be faster than the non-mutating ones, as opposed to what was observed in the 1D case.

In the non-mutating 3D implementations, the runtime spent in array allocations becomes dominant when compared to that of the vectorized operations between arrays.

The aftermath of the study conducted in this subsection is that the memory footprint of the adjoint of a non-mutating program when compared to its primal is at least twice that of their mutating counterpart. Along with the elements in the literature already listed [START_REF] Bezgin | JAX-Fluids: A fully-differentiable high-order computational fluid dynamics solver for compressible two-phase flows[END_REF][START_REF] Kochkov | Machine learning-accelerated computational fluid dynamics[END_REF], such a result corroborates with the assumption that building efficient differentiable PDE solvers requires support to array mutation by the AD tool of choice. In the next subsection, a similar analysis is performed on mutating and non-mutating CFD kernels.

B.3 Analysis of a CFD kernel

In order to evaluate the impact of the mutation problem in a CFD configuration, this section considers a function that updates a numerical solution to the 3D Euler equations on tetrahedral meshes using the Two-step Taylor-Galerkin C (TTGC) discretization scheme depicted in Sec. 5.1. No treatment of boundary conditions is performed here. The non-mutating and mutating programs applying a TTGC update will be hereafter referred to as dU_ttgc_reassignment and dU_ttgc_mutation, respectively. They both return the squared l2-norm of the value ∆U by which some input numerical solution U should be incremented over one time step. For the sake of shortness, only parts of the algorithms of the methods dU_ttgc_reassignment and dU_ttgc_mutation are represented in The part of the TTGC scheme highlighted in Fig. B.8 relate to the computation of the nodal residuals Rj ∈ R neq×nnode by scattering of the cell residuals VeRe ∈ R neq×ncell (weighted by the cell volumes), with neq ≡ 5 being the number of equations to be solved (one for density, three for momentum and one for energy), nnode the number of mesh nodes and ncell the number of mesh cells. This is a common operation in the class of residual distribution numerical schemes, which TTGC belongs to [START_REF] Deconinck | Residual Distribution Schemes: Foundations and Analysis[END_REF]. One can notice that the mutating kernel accumulates cells contribution to a zero-initialized array Rj inside a loop over mesh cells (lines 5-17 of dU_ttgc_mutation). By allowing mutation, one avoids to store VeRe for .9 presents performance measurements for the adjoints of both methods, including not only memory usage, but also runtime benchmarks. The computations were performed using an Intel(R) Xeon(R) Gold 6140 CPU. The adjoints of dU_ttgc_reassignment and dU_ttgc_mutation are respectively generated by the Zygote package (which only supports reassignments, as seen in Sec. 1.3.4) and by the Enzyme library (which supports mutation, as seen in Sec. 1.3.5). At the largest problem size under consideration (namely, nnode = 33792), Fig. B.9(a) reveals that the adjoint of dU_ttgc_reassignment allocates around 250 times the amount of memory allocated by the method, whereas the adjoint of dU_ttgc_mutation is only twice as expensive as its primal. Fig. B.9(b) shows that the adjoint of the non-mutating kernel is around 100 times slower than its primal, whereas the adjoint of the mutating kernel is less than 4 times slower than its primal. This means that, in general, one can expect the differentiation of non-mutating PDE kernels to perform much worse than that of mutating implementations, both runtime-and memory-wise.

The quantitative observations from Sec. B.2-B.3 should be sufficient to motivate the development of AD tools that support array mutation for enabling efficient training of ML-PDE coupled systems. Tapenade (Sec. 1.3.3) and other tools developed for Fortran/C languages natively support mutation, but seem out of interest of the ML community for their relatively low-level features when compared to those of the Python APIs of the main ML frameworks (namely, TensorFlow, PyTorch and JAX). Once the field of ML for scientific applications leans toward the development of coupled ML-PDE systems (in a quest for more generalizable predictive capabilities of ML architectures acting on physical models), the mutation problem should be given more attention. Eq. (C.4) closely resembles the shock sensor operator from the widely acclaimed JST dissipation scheme of [START_REF] Jameson | Origins and Further Development of the Jameson-Schmidt-Turkel Scheme[END_REF]. Assuming an equally spaced mesh (i.e., τ = h R /h L ≡ 1) and only right moving waves (kh is positive), the following inequality with the JST sensor (S JST ) holds:

S JST = |∆u i+1/2 -∆u i-1/2 | |∆u i+1/2 | + |∆u i-1/2 | ≥ |∆u i+1/2 | -|∆u i-1/2 | |∆u i+1/2 | + |∆u i-1/2 | . (C.5)
The shock sensor is a conservative upper bound for the local wavenumber in terms of the Nyquist frequency (on equally spaced meshes). In fact, [START_REF] Jameson | Origins and Further Development of the Jameson-Schmidt-Turkel Scheme[END_REF] mentions that this sensor provides robust results near shocks compared to other TVD shock sensors and inspired the CUSP family of schemes.

To the author's knowledge, this is the first work to reveal this mathematical connection of the local wavenumber to the shock sensor. inviscid Burgers' equation, the flux function is f (u) = 1 2 u 2 and the entropy condition can be rewritten as

u 1 + u ≤ 2ξ ≤ u + u 2 .
(D.4)

The problem of shock capturing using discrete finite-elements as illustrated in Fig. 3.8 is considered. It is intended to smooth the shock and capture its location between the two elements e i ∈ [i -1, i] and e i+1 ∈ [i, i + 1], say by matching the impedance to the left and right of node i. The entropy condition in Eq. (D.3) must be satisfied for a valid shock solution between the elements and this simplifies to

u n+1 i-1 ≤ u n+1 i+1 . (D.5)
The non-linear Burgers' equation can be linearized for the time step n to n + 1 and the averaged convection velocity resulting from the linearization [START_REF] Griewank | Evaluating Derivatives[END_REF][START_REF] Ucar | Numerical Solutions of the Modified Burgers' Equation by Finite Difference Methods[END_REF] is assumed constant over the given finite-element. The Fourier transform is applied on both sides of the inequality and simplified to yield the entropy condition in spectral space,

|G i-1 | ≤ |G i+1 |. (D.6)
Note that the impedance G is based on the LTA of linearized convection equation with locally constant (average) convection speed [START_REF] Griewank | Evaluating Derivatives[END_REF][START_REF] Ucar | Numerical Solutions of the Modified Burgers' Equation by Finite Difference Methods[END_REF]. In addition, the nodal impedance G i-1 and G i+1 are functions of the element impedance pairs [G R i-2 , G L i-1 ] and

  (b) Frank Rosenblatt manipulating wired connections of the Mark I Perceptron machine (1960). From dnvdk (2022). . . . . . . . . . . . . . . . . 1.4 Example of linearly separable sets. . . . . . . . . . . . . . . . 1.5 Schematic representation of a 3-layer MLP. . . . . . . . . . . 1.6 Examples from the MNIST dataset. . . . . . . . . . . . . . . 1.7 Schematic comparison between a dense layer (top) and a convolutional layer (bottom). . . . . . . . . . . . . . . . . . . . . 1.8 Schematic representation of operations performed in a graph network layer. . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Performance metrics for training and validation sets at the end of each epoch during the training procedure. . . . . . . . . . 1.10 Absolute error of the ND-evaluated derivative of the function f (x) :

  2.1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Evolution of training and validation loss values across epochs for independent CNN-PAL models trained to predict N step = 16 time steps of the convection of vortices from training datasets composed of (a) N s = 32 samples; (b) N s = 128 samples. Results concern experiment n.2 (see Tbl. 2.2). . . . . . . . . . . . 2.10 Grid convergence error plots for the WENO5 scheme implementation for the Euler equations considering either (a) mean L 2 error or (b) mean L ∞ error of density field over one box turn of the baseline vortex depicted in Tbl. 2.2. The variable h stands for the grid cell size along a single direction. . . . . . 2.11 Outcomes of experiment n.3 performed with two independent CNN-PAL models trained to predict (a) N step = 16 time steps; (b) N step = 32 time steps of the evolution of the double shear layer initial condition depicted in Tbl. 2.2. The training proceure employs the loss function (2.50). (Top) Evolution of training loss values across epochs. (Center) Evolution of L 2 error of density field over time. (Bottom) Absolute error of density field predicted by the trained CNN-PAL models after 40 simulation time steps. . . . . . . . . . . . . . . . . . . . . . 2.12 Outcomes of experiment n.3 performed with a CNN-PAL model trained to predict N step = 32 time steps of the evolution of the double shear layer initial condition depicted in Tbl. 2.2. The training procedure employs the loss function (2.51). (a) Evolution of training loss values across epochs. (b) Evolution of L 2 error of density field over time. (c) Evolution of the sum of L 2 errors of density, components of momentum along the x and y directions, and energy fields over time. . . . . . . . . . .2.13 Evolution of energy values in the simulation of the double shear layer from experiment n.3 (see Tbl. 2.2). The growth in energy values for the trained CNN-PAL models after a certain time demonstrate numerically unstable dynamics. . . . . . . . 3.1 1D linear finite-element shape function (hat function) with nodes and elements enumerated. . . . . . . . . . . . . . . . . . 3.2 Comparison of column-wise (i) |G num |/|G phy |; grey area shows unstable region |G| > 1, (ii) c N /c; grey area shows region where 0.95 ≤ c N /c ≤ 1.05 and dotted line is p = 1.5 and (iii) v gn /v g ; grey area shows GVP region 0.95 ≤ v gn /v g ≤ 1.05 and dotted line is p = 1 for the indicated methods (row-wise (a-e)) tabulated in tbl. 3.1. . . . . . . . . . . . . . . . . . . . . . . . 3.3 Inverse mass-matrix spectral plot ( A = M -1 ) compared with the dissipation D to illustrate its anti-diffusive character and the stabilization effect of β c . . . . . . . . . . . . . . . . . . . . 3.4 Dissipation (solid contour) and GVP region (blue contour) for various β c values in TG-LW scheme; wavy regions are numerically unstable. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 |G num |/|G phy | contours for indicated values of parameter γ in the TTGC-γ scheme; (i) grey shaded area indicates numerically unstable regions (|G num | > 1) and (ii) blue shaded region indicates the group velocity preserving zone (0.95 ≤ v gn /v g ≤ 1.05). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Global CFL number versus maximum stable γ = γ max ; stable region obtained from GSA indicated in grey shaded area. . . . 3.7 Convection of WP on mesh with discontinuity at junction (♦ is the junction) showing upstream reflected waves indicated by ←-and actual wave propagation direction indicated by -→ c; are the mesh nodes. . . . . . . . . . . . . . . . . . . 3.8 Cell and nodal impedance matching problem using local-transfer function; the dash-pot models externally added artificial dissipation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 |G| contours of LTA nodal impedance (a) |G i-1 | (left node), (b) |G i | (junction node) and and (c) |G i+1 | (right); solid shaded region is numerical phase speed within ±5% error of exact phase speed (DRP region), vertical striped region indicates the spectral spread of the WP in terms of local Nyquist frequency p; contours of local instability shown in red (thicker) contour lines (1, 1.1). . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Contour of reflection ratio magnitude (|ρ|) due to nodal impedance mismatch between (a) (i-1)-(i) and (b) (i)-(i+1) for the junction problem for varying CFL (0 < N c < 1) and frequency (0 < p < π/2). . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Space-wavenumber idealization and local neighborhood of nodal and elemental impedance to the continuous initial condition u(x, t n ) leading to the solution at the next step u(x, t n+1 ). . . 3.12 LTA model at a Burgers shock(a) initial condition (solid line) and appearance of solution over-/undershoots (dotted line) created by the reflected waves from impedance mismatch and the individual nodal impedance (b-d) from elemental ones. . . 3.13 Dissipation |G| and phase angle ϕ of the local nodal transfer function G i+1 , G i , and G i-1 ; grey shaded region are numerically unstable (|G| > 1). Note that the phase angle mismatch between G i-1 and G i is much more severe than G i and G i+1 . . 3.14 Time evolution of numerical solution zoomed near the Burgers' shock for the first (left) and second (right) time step from the initial condition. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 |G| contours of local transfer function at shock location i for the TTGC-γ scheme for Burgers' equation (a) for variant (i) with Pe = 0.5, (b) variant (ii) with Pe = 0.25 and (c) variant (ii) with Pe = 0.375. . . . . . . . . . . . . . . . . . . . . . . . 3.16 Time step t = 16∆t for the evolution of a square wave initial condition under Burgers' equation, where ∆t is the simulation time step. (a) for variant (i) with Pe = 0.5, (b) variant (ii) with Pe = 0.25 and (c) variant (ii) with Pe = 0.375. Red lines delimit region near shock front where artificial viscosity is added. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 (a) 1D primal/dual mesh (top) and GN based on the dual mesh (bottom); (b) Unstructured primal/dual mesh in 2D (left) and GN based on the dual mesh (right). Arrows denote the GN message-passing, red boxes are the GN vertices (which are equivalent to the dual mesh nodes, indicated by red dots), blue dots are the primal mesh nodes and green dotted lines are the dual mesh edges. . . . . . . . . . . . . . . . . . . . . . 4.7 ML-TTGC architecture. The GNN instances trained for the studies presented in Sec. 4.4 use an encoding block, two messagepassing layers and a decoding block, which update the vertices and edges attributes using dense layers with a hidden feature space of size L = 32. Each instance has a total of 10, 637 trainable parameters. . . . . . . . . . . . . . . . . . . . . . . . 4.8 (a) Low-frequency WP convection over an irregularly spaced mesh. (b) High-frequency WP convection over irregular mesh. Top: TTGC-γ (Default) and analytical/exact solution. Center: ML-TTGC and analytical/exact solution. Bottom: GNN predicted parameters for locally tunable TTGC-γ L . Initial conditions are not part of the training set. . . . . . . . . . . . 4.9 (a) High-frequency WP convection over regular mesh. Top: TTGC-γ (Default) and analytical/exact solution. Bottom: ML-TTGC and analytical/exact solution. Initial condition is not part of the training set. (b) L2 norm of the error for TTGC-γ (Default) and ML-TTGC, comparing their dissipation levels over time. . . . . . . . . . . . . . . . . . . . . . . . 4.10 Numerical solution obtained at time steps (a) t = 35∆t and (b) t = 40∆t (TTGC-γ and ML-TTGC scheme with predicted parameters γ opt and ν opt ) to the inviscid Burgers' equation for the square-wave initial condition, where ∆t is the simulation time step. Parameter values used are (a, x 0 , x 1 ) = (1., 0.2, 0.65) for (a) on the left and (a, x 0 , x 1 ) = (2., 0.4, 0.75) for (b) on the right; data is not part of the training set. . . . . . . . . 4.11 Numerical solution obtained at time steps (a) t = 40∆t and (b) t = 70∆t (TTGC-γ and ML-TTGC with predicted γ opt and ν opt ) to the inviscid Burgers' equation for the sine wave initial condition, where ∆t is the simulation time step; data is not part of the training set. . . . . . . . . . . . . . . . . . . . 5.1 Different finite volume schemes: (a) cell-centered, (b) cellvertex. Inspired from Crumpton et al. (1993). . . . . . . . . . 5.2 Illustration of cell residuals computation via trapezoidal rule (a) and of their scattering to nodes (b). . . . . . . . . . . . . . 5.3 (a) Unstructured regular triangular grid and (b) its perturbed counterpart. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Surface plots of wave packet whose dynamics are governed by the 2D Convection equation. Time step n. 360. (a) Analytical solution, (b) Output from TTGC scheme, (c) Output from ML-TTGC scheme trained by pure l2 loss function, (d) Output from ML-TTGC scheme trained by pure TVD loss function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 (a, b, c) Surface plots and (d) cut along x = 0.5 of isentropic vortex after completing one box turn (periodic conditions apply). (a) Analytical density field, (b) Output density field from TTGC scheme, (c) Output density field from ML-TTGC scheme trained on pure l2 loss function. . . . . . . . . . . . . . 5.6 Dynamics of fluid's momentum component along the x direction (ρu) for double shear layer problem. (a) TTGC time step n. 600 (just before simulation crash), (b) ML-TTGC time step n. 600, (c) ML-TTGC time step n. 2400, (d) ML-TTGC time step n. 3600. . . . . . . . . . . . . . . . . . . . . . . . . . . . B.1 Illustration of x[0] += 1. operation done either (a) in-place (actual mutation of x) or (b) by re-assignement of x. . . . . . B.2 Illustration of how du is computed either by (a) du_mutation or (b) du_reassignment. . . . . . . . . . . . . . . . . . . . . B.3 Computational graph, forward primal and adjoint evaluation traces for the program du_reassignment. . . . . . . . . . . . B.4 Schematic representation of circshift operation on some array w. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.5 Adjoint programs of (a) du_mutation (Fig. B.2(a)) and (b) du_reassignment (Fig. B.2(b)). Note that the arrays ub and dub respectively corresponding to the reverse derivatives of u and du are initialized to zeros beforehand and given as inputs to each function. . . . . . . . . . . . . . . . . . . . . . . . . . B.6 Benchmarking results for the 1D discretization operator problem regarding the ratios of (a) memory usage; (b) runtime between adjoint and primal codes. . . . . . . . . . . . . . . . B.7 Benchmarking results for the 3D discretization operator problem regarding the ratios of (a) memory usage; (b) runtime between adjoint and primal codes. . . . . . . . . . . . . . . . B.8 Portions of the algorithms of the (a) mutating and (b) nonmutating TTGC scheme implementations (3D Euler equations; tetrahedral meshes). . . . . . . . . . . . . . . . . . . . . . . . B.9 Benchmarking results for the TTGC scheme (3D Euler equations; tetrahedral meshes) regarding the ratios of (a) memory usage; (b) runtime between adjoint and primal codes. . . . . C.1 Comparison of local wavenumber analysis using Eikonal and JST sensor; (a-c) WP function and (d-f) Gaussian pulse function.211 D.1 (a) Thermodynamic cycle of a shock wave, (b) waves overtaking each other, (c) multiple solution at shock, (d) the shock entropy solution and (e) x-t diagram of the shock problem. . .

Figure 1

 1 Figure 1.3: (a) Camera device capturing an image of a letter C drawn on a wall (1962). From Rosenblatt (1962). (b) Frank Rosenblatt manipulating wired connections of the Mark I Perceptron machine (1960). From dnvdk (2022).

Figure 1 . 4 :

 14 Figure 1.4: The sets of blue and red dots are linearly separable, as indicated by the dashed straight lines. Infinitely many separation lines exist and the perceptron learning algorithm cannot distinguish between them.

  Fig. 1.6 shows examples of images from this dataset.

Figure 1 . 6 :

 16 Figure 1.6: Examples from the MNIST dataset. From Lecun et al. (1998).

Figure 1 .

 1 Figure 1.13: (Top) Schematic representation of a dense layer. (Middle and bottom) Computational graph, forward primal and adjoint evaluation traces for a program evaluating the (scalar) output of the dense layer.

Figure 2 . 1 :

 21 Figure 2.1: Representations of the (a) 1D Riemann problem (Def. 4) and (b) 1D periodic Riemann problem (Def. 5).

Figure 2 . 2 :

 22 Figure 2.2: Representations of the analytical solutions for the 1D periodic Riemann problem (see Def. 5) following (a) the Convection equation (2.4) and (b) the Burgers' equation (2.5). In (b), c := (u max + u min ) / 2.

Figure 2

 2 Figure 2.3: Representation of MUSCL numerical fluxes f n i±1/2 using the extrapolations u R, ni±1/2 and u L, n i±1/2 of the state variable at the boundaries of the cell [x i-1/2 , x i+1/2 ], where x i±1/2 := (x i + x i±1 )/2.

Figure 2 Figure 2

 22 Figure 2.4: (a) One box turn of a 1D periodic Riemann problem (see Def. 5) following the Convection equation (2.4) for all schemes depicted in Sec. 2.2. (b) Numerical solution of a 1D periodic Riemann problem following the Burgers' equation (2.5) at time t = 150∆t (where ∆t is the simulation time step) for all schemes depicted in Sec. 2.2.

Fig. 2 .

 2 Fig. 2.6 illustrates the neural network architecture. It consists of a Convolutional Neural Network (CNN) model (see Def.2) followed by a non-trainable layer called the Polynomial Accuracy Layer (PAL). The density (ρ), the components of the momentum along the x and y directions (ρu and ρv, respectively), and the energy (E) fields of a given initial condition are used as input features of the CNN. Subsequently, the CNN output features are split among eight channels. Each channel corresponds to a spatial derivative to be approximated (one per scalar flux, see Eq. (2.32)). The PAL projects the features of each channel on the space of weights that approximate spatial derivatives at a given enforced order of accuracy. Finally, the third-order SSP Runge-Kutta time-stepping method (2.29) is used to evolve the initial condition over time.In order to explain how the PAL works, a 1D field f := f (x) is considered, for the sake of simplicity. The PAL computes a set of weights that is used to approximate the ℓ-th derivative f (ℓ) (x 0 ) of the field f for a given number N of points in the neighborhood of x 0 . The Taylor series of f around the neighbor x n := x 0 + h n (h n > 0) can be expressed as follows:

Figure 2

 2 Figure 2.6: Schematic representation of the CNN-PAL neural network architecture for the data-driven spatial discretization of the Euler equations on 2D cartesian grids. The variables F ⋆,x and F ⋆,y correspond to the components of the flux of the field ⋆ along the x and y directions, respectively.

  and l charact := 1, where M ∞ is the Mach number and c ∞ := γ gas R gas T ∞ is the speed of sound in the initialized homogeneous medium. Tbl. 2.2 illustrates the values of parameters used to set the vortex initial condition. Experiment n.1: Training on a single isentropic vortex Parameter Value Vortex initial position (x 0 , y 0 2: Training on a set of isentropic vortices Parameter Value Vortex initial position (x 0 , y 0 ) (0.5, 0.5) Mach number (M ∞ ) values used in all experiments of Sec. 2.3.2. U(a, b) stands for a continuous uniform distribution on the real interval (a, b). First introduced by Courant, Friedrichs, et al.

Figure 2

 2 Figure 2.7: Outcomes of experiment n.1 performed with two independent CNN-PAL models trained to predict (a) N step = 4 time steps; (b) N step = 16 time steps of the evolution of the vortex depicted in Tbl. 2.2. (Top) Evolution of training loss values across epochs. (Center) Evolution of L 2 error of density field over time. (Bottom) Absolute error of density field predicted by the trained CNN-PAL models after one box turn of the vortex.

Figure 2

 2 Figure 2.9: Evolution of training and validation loss values across epochs for independent CNN-PAL models trained to predict N step = 16 time steps of the convection of vortices from training datasets composed of (a) N s = 32 samples; (b) N s = 128 samples. Results concern experiment n.2 (see Tbl. 2.2).

Figure 2

 2 Figure 2.10: Grid convergence error plots for the WENO5 scheme implementation for the Euler equations considering either (a) mean L 2 error or (b) mean L ∞ error of density field over one box turn of the baseline vortex depicted in Tbl. 2.2. The variable h stands for the grid cell size along a single direction.

Figure 2

 2 Figure 2.11: Outcomes of experiment n.3 performed with two independent CNN-PAL models trained to predict (a) N step = 16 time steps; (b) N step = 32 time steps of the evolution of the double shear layer initial condition depicted in Tbl. 2.2. The training proceure employs the loss function (2.50). (Top) Evolution of training loss values across epochs. (Center) Evolution of L 2 error of density field over time. (Bottom) Absolute error of density field predicted by the trained CNN-PAL models after 40 simulation time steps.

Figure 2

 2 Figure 2.12: Outcomes of experiment n.3 performed with a CNN-PAL model trained to predict N step = 32 time steps of the evolution of the double shear layer initial condition depicted in Tbl. 2.2. The training procedure employs the loss function (2.51). (a) Evolution of training loss values across epochs. (b) Evolution of L 2 error of density field over time. (c) Evolution of the sum of L 2 errors of density, components of momentum along the x and y directions, and energy fields over time.

Figure 2 .

 2 Figure 2.13: Evolution of energy values in the simulation of the double shear layer from experiment n.3 (see Tbl. 2.2). The growth in energy values for the trained CNN-PAL models after a certain time demonstrate numerically unstable dynamics.

Figure 3

 3 Figure 3.1: 1D linear finite-element shape function (hat function) with nodes and elements enumerated.

  3.1. The contour plot of the three ratios for the numerical schemes shown in Tbl. 3.1 are plotted in Fig. 3.2 versus CFL number (N c ) and non-dimensional wavenumber (p = kh). All figures stacked column-wise are the ratios |G num |/|G phy |, c N /c and v gn /v g , while stacked row-wise are the various numerical methods forming a figure matrix.

Figure 3

 3 Figure 3.2: Comparison of column-wise (i) |G num |/|G phy |; grey area shows unstable region |G| > 1, (ii) c N /c; grey area shows region where 0.95 ≤ c N /c ≤ 1.05 and dotted line is p = 1.5 and (iii) v gn /v g ; grey area shows GVP region 0.95 ≤ v gn /v g ≤ 1.05 and dotted line is p = 1 for the indicated methods (row-wise (a-e)) tabulated in tbl. 3.1.

Figure 3 . 3 :

 33 Figure 3.3: Inverse mass-matrix spectral plot ( A = M -1 ) compared with the dissipation D to illustrate its anti-diffusive character and the stabilization effect of β c .

Figure 3

 3 Figure 3.4: Dissipation (solid contour) and GVP region (blue contour) for various β c values in TG-LW scheme; wavy regions are numerically unstable.

Figure 3 . 5 :

 35 Figure 3.5: |G num |/|G phy | contours for indicated values of parameter γ in the TTGC-γ scheme; (i) grey shaded area indicates numerically unstable regions (|G num | > 1) and (ii) blue shaded region indicates the group velocity preserving zone (0.95 ≤ v gn /v g ≤ 1.05).

Figure 3 Figure 3 . 8 :

 338 Figure 3.7: Convection of WP on mesh with discontinuity at junction (♦ is the junction) showing upstream reflected waves indicated by ←-and actual wave propagation direction indicated by -→ c; are the mesh nodes.

  Figure 3.9: |G| contours of LTA nodal impedance (a) |G i-1 | (left node), (b) |G i | (junction node) and and (c) |G i+1 | (right); solid shaded region is numerical phase speed within ±5% error of exact phase speed (DRP region), vertical striped region indicates the spectral spread of the WP in terms of local Nyquist frequency p; contours of local instability shown in red (thicker) contour lines (1, 1.1).

Figure 3

 3 Figure 3.10: Contour of reflection ratio magnitude (|ρ|) due to nodal impedance mismatch between (a) (i -1)-(i) and (b) (i)-(i + 1) for the junction problem for varying CFL (0 < N c < 1) and frequency (0 < p < π/2).

  .73) To minimize the error between the numerical and the exact solution, one should minimize the l2 error objective L2(u n+1 ). A minimum objective is reached when the minimum satisfies the conditions, |G j | ≈ |G phy | and ϕ j = ϕ phy , i.e., the numerical and physical transfer functions match. This naturally means that the wave power is maximized, since |G j | ≈ |G phy | ensures minimum dissipation. Comparing Eq. (3.72) and Eq. (3.73) with the individual impedances of the left and right elements in Eq. (3.60) it is trivial to show that the elemental transfer function has tendency to have opposing elemental phase angles ϕ L ≈ -ϕ R . Therefore, the minimization of L2(u n+1 ) is a good approximation to the maximum power impedance matching G L ≈ G † R . For the minimization of L2(u n+1 ) to perform maximum power matching it is critical to impose the stability constraint on |G L |, |G R |, i.e., |G L | 2 + |G R | 2 ≤ 1 and bound them close to the physical impedance value.

Figure 3 .

 3 Figure 3.12: LTA model at a Burgers shock (a) initial condition (solid line) and appearance of solution over-/undershoots (dotted line) created by the reflected waves from impedance mismatch and the individual nodal impedance (b-d) from elemental ones.

Figure 3 .

 3 Figure 3.13: Dissipation |G| and phase angle ϕ of the local nodal transfer function G i+1 , G i , andG i-1; grey shaded region are numerically unstable (|G| > 1). Note that the phase angle mismatch between G i-1 and G i is much more severe than G i and G i+1 .

Figure 3

 3 Figure 3.14: Time evolution of numerical solution zoomed near the Burgers' shock for the first (left) and second (right) time step from the initial condition.

Figure 3 Figure 3

 33 Figure 3.15: |G| contours of local transfer function at shock location i for the TTGC-γ scheme for Burgers' equation (a) for variant (i) with Pe = 0.5, (b) variant (ii) with Pe = 0.25 and (c) variant (ii) with Pe = 0.375.

Figure 4 . 1 :

 41 Figure 4.1: Study n.1. Spatio-temporal evolution of numerical error in the solution to the WP convection (Eq. (4.3)); WP and its Fast Fourier Transform (FFT) along with the comparison of numerical errors of optimized and unoptimized TTGC and the optimal γ NLopt obtained using NLopt are plotted for central wavenumbers (a-c) k 0 h = 0.4, (d-f) k 0 h = 0.6, and (g-i) k 0 h = 0.8 respectively.

Figure 4

 4 Figure 4.2: (a) Numerical error convergence with mesh refinement (b-e) (finest to coarsest) and the dotted line indicates the third order slope; initial condition and FFT for the mesh refinement levels containing (b) 256, (c) 128, (d) 64 and (e) 32 elements respectively.

Figure 4 Figure 4

 44 Figure 4.3: Study n.2. (a) Left: Periodic initial condition over an irregular mesh of 101 nodes (or 100 finite-elements). Parameters values are (k 0 h, α, x 0 , c) = (0.5, 25, 0.5, 1). Right: Fast-Fourier Transform evaluated over the regular counterpart of the actual mesh; (b) Pareto front for varying θ ∈ [0., 1.] comparing normalized time-averaged values of the two objective functions from the optimization problem (mean over 100 steps); (c, d) Comparison of NLopt and TTGC-γ outputs at steps n. 40 and 120, respectively.

Figure 4

 4 Figure 4.5: Study n.4. (a) Initial condition and its frequency domain counterpart. Parameter a = 1; (b) Pareto front for varying θ ∈ [0., 1.] comparing normalized time-averaged values of the two objective functions from the optimization problem (mean over 100 steps); (c, d) Comparison of exact solution, NLopt and TTGC-γ outputs at steps n. 40 and 120, respectively.

  Fig. 4.8

Figure 4

 4 Figure 4.8: (a) Low-frequency WP convection over an irregularly spaced mesh. (b) High-frequency WP convection over irregular mesh. Top: TTGCγ (Default) and analytical/exact solution. Center: ML-TTGC and analytical/exact solution. Bottom: GNN predicted parameters for locally tunable TTGC-γ L . Initial conditions are not part of the training set.

Figure 4

 4 Figure 4.9: (a) High-frequency WP convection over regular mesh. Top: TTGC-γ (Default) and analytical/exact solution. Bottom: ML-TTGC and analytical/exact solution. Initial condition is not part of the training set. (b) L2 norm of the error for TTGC-γ (Default) and ML-TTGC, comparing their dissipation levels over time.

Figure 4 Figure 4

 44 Figure 4.10: Numerical solution obtained at time steps (a) t = 35∆t and (b) t = 40∆t (TTGC-γ and ML-TTGC scheme with predicted parameters γ opt and ν opt ) to the inviscid Burgers' equation for the square-wave initial condition, where ∆t is the simulation time step. Parameter values used are (a, x 0 , x 1 ) = (1., 0.2, 0.65) for (a) on the left and (a, x 0 , x 1 ) = (2., 0.4, 0.75) for (b) on the right; data is not part of the training set.

Figure 5 . 1 :

 51 Figure 5.1: Different finite volume schemes: (a) cell-centered, (b) cell-vertex. Inspired from Crumpton et al. (1993).

  Figure 5.2: Illustration of cell residuals computation via trapezoidal rule (a) and of their scattering to nodes (b).

Figure 5

 5 Figure 5.3: (a) Unstructured regular triangular grid and (b) its perturbed counterpart.

Figure 5

 5 Figure 5.4: Surface plots of wave packet whose dynamics are governed by the 2D Convection equation. Time step n. 360. (a) Analytical solution, (b) Output from TTGC scheme, (c) Output from ML-TTGC scheme trained by pure l2 loss function, (d) Output from ML-TTGC scheme trained by pure TVD loss function.

Figure 5

 5 Fig.5.5 presents inference results for a vortex from the training set. ML-TTGC was trained on pure l2 loss function only, as reportedly giving smoothier dynamics in the Convection case. Once again, wiggles are significantly damped in the predicted dynamics, approaching it from the analytical one. Despite the limited range of initial conditions explored during training, ML-TTGC also showed generalization capabilities for stabilizing the dynamics of a double shear layer initial condition it was never trained on, as exposed in Fig.5.6.

  diff variables: x:out y:out z:in-zero SUBROUTINEFUNC _ B(n, x, xb, y, yb, z, zb)

  recent call last): File "<stdin>", line 1, in <module> TypeError: 'tensorflow.python.framework.ops.EagerTensor' object does not support item assignment PyTorch : >>> import torch >>> torch. __ version __ '1.11.0' >>> x = torch.rand(10, requires _ grad=True) >>> >>> x[0] += 1. Traceback (most recent call last):

Figure B. 1 :

 1 Figure B.1: Illustration of x[0] += 1. operation done either (a) in-place (actual mutation of x) or (b) by reassignment of x. Represented by red boxes, additional memory in the forward pass is required in (b) when compared to (a).

Fig. B. 2

 2 Fig. B.2 represents the mutating and non-mutating programs. One can notice that the forward pass of du_reassignment (Fig. B.2(b)) requires additional

  Figure B.2: Illustration of how du is computed either by (a) du_mutation or (b) du_reassignment. Represented by red boxes, additional memory in the forward pass is required in (b) when compared to (a). Note that an array du initialized to zeros beforehand is given as input to each function. Furthermore, the instruction sum(abs2, v) computes the squared l2-norm ||v|| 2 2 := i |v i | 2 of a given vector v.

Figure B. 3 :

 3 Figure B.3: Computational graph, forward primal and adjoint evaluation traces for the program du_reassignment.

Figure B. 4 :

 4 Figure B.4: Schematic representation of circshift operation on some array w. The sensitivity of the output array w_shifted is propagated back to the input array w by performing the reverse operation, i.e., circshift(w, shift, w_shifted) := circshift(w_shifted, -shift) (shift ∈ Z).

d

  ∥d x u(x, y, z)∥ 2 2 + ∥d y u(x, y, z)∥ 2 2 + ∥d z u(x, y, z)x u(x, y, z)| x=x i ,y=y j ,z=z k := u i-2,j,k -4u i-1,j,k + 3u i,j,k (B.16) d y u(x, y, z)| x=x i ,y=y j ,z=z k := u i,j-2,k -4u i,j-1,k + 3u i,j,k (B.17) d z u(x, y, z)| x=x i ,y=y j ,z=z k := u i,j,k-2 -4u i,j,k-1 + 3u i,j,k (B.18) (i ∈ 1, N , j ∈ 1, N , k ∈ 1, N , N ∈ N)are benchmarked 8 . Fig.B.7 compares performances of the adjoint programs with respect to their primals for varying problem sizes (N ). The experiments were conducted under conditions similar to those of the 1D case.

  Fig. B.8 9 .

Figure B. 8 :

 8 Figure B.8: Portions of the algorithms of the (a) mutating and (b) nonmutating TTGC scheme implementations (3D Euler equations; tetrahedral meshes).

  Figure C.1: Comparison of local wavenumber analysis using Eikonal and JST sensor; (a-c) WP function and (d-f) Gaussian pulse function.

  2 | -|∆u i-1/2 |τ |∆u i+1/2 | + |∆u i-1/2 |τ . (C.4)

  Figure D.1: (a) Thermodynamic cycle of a shock wave, (b) waves overtaking each other, (c) multiple solution at shock, (d) the shock entropy solution and (e) x-t diagram of the shock problem.

1

  Manuscript's structure. . . . . . . . . . . . . . . . . . . . . .

1.1 Arthur

Samuel and IBM 700 (1956)

. . . . . . . . . . . . . . . 1.2 Schematic representation of the original (photo-)perceptron. . 1.3 (a) Camera device capturing an image of a letter C drawn on a wall

  The symbol u i represents an element of the dataset U which the model learns from, also known as training set. The loss function L compares the model's predictions from the given input u i with a reference value r i . The symbol r i represents the corresponding sample from the dataset R of reference values. Finally, the variable N s stands for the number of samples of the datasets U and R.

	Usually, gradient-based optimization techniques are employed to find W ⋆ .
	More precisely, reverse mode of algorithmic differentiation (see Sec. 1.2.2) 6

end for 1.2 An introduction to Algorithmic Differen- tiation (AD)

  u sample , r sample )

	11:	end for	▷ Estimate loss of the entire batch
	12:
	13:	W ← OPT(W, ∂L/∂W)	▷ Update trainable parameters
	14:
	15:	end for
	16:
	17:	▷ Evaluate the network on the entire validation set every epoch
	18:	L val ← 0
	19:	for v ∈ V do
	20:	L val ← L val + P (W, v)
	21:	end for
	22:
	23:	▷ Evaluate the network on the entire training set every epoch
	24:	L train ← 0
	25:	for u ∈ U do
	26:	L train ← L train + P (W, u)
	27:	end for
	28:
	29:	▷ Check stopping criteria
	30:	if S(L val , L train ) == True then
	31:	save(W)	▷ Save current ML model state
	32:	Break	▷ Interrupt training loop execution
	33:	end if
	34:
	35:

  The ultimate purpose of the works of[START_REF] Bar-Sinai | Learning data-driven discretizations for partial differential equations[END_REF] and[START_REF] Kochkov | Machine learning-accelerated computational fluid dynamics[END_REF] is to replace subgrid-scale models with neural networks that provide suitable discretization of the Navier-Stokes equations.[START_REF] Bar-Sinai | Learning data-driven discretizations for partial differential equations[END_REF] demonstrated the approach on 1D problems governed by the viscous Burgers' equation (2.6), the Korteweg-de Vries equation, which is a mathematical model of waves on shallow water surfaces[START_REF] Korteweg | On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves[END_REF], and the Kuramoto-Sivashinsky equation, which models diffusive-thermal instabilities in laminar flame fronts

Table 2

 2

	Parameter	Value
	Number of convolutional layers	4
	Number of hidden features	16
	Size of convolutional kernel/filter	(3, 3)
	Size of input	(4, 33, 33)
	Size of stencil (N )	5
	Constrained order of accuracy (m)	2

.1: Parameters common to all CNN-PAL models used in experiments.

Hardware-wise, the weights were modelled by variable resistors (also known as potentiometers) to be tuned during the learning process.

In binary classification, two sets of points are linearly separable if and only if there exists a line in the Euclidean plane separating them. This is better illustrated in Fig.1.4.

The S-, A-and R-layers from Fig.1.2 are grouped to compose a single-layer perceptron.

 4 In the modern terminology, it is named activation function.

For example, a CNN employed as a binary classifier to distinguish images of cats from images of dogs will be able to identify a cat's ear independent of its position in the frame.

It was introduced by[START_REF] Linnainmaa | The representation of the cumulative rounding error of an algorithm as a Taylor expansion of the local rounding errors[END_REF] and first applied to neural networks training in 1974 by[START_REF] Werbos | Beyond regression: new tools for prediction and analysis in the behavioral sciences[END_REF] (nevertheless, this work was not cited by[START_REF] Rumelhart | Learning representations by back-propagating errors[END_REF]).

The reverse mode of AD was rediscovered independently by[START_REF] Speelpenning | Compiling fast partial derivatives of functions given by algorithms[END_REF].

The subscript (•) xx = ∂

∂x 2 denote the second-order partial derivative with respect to space.

For the sake of simplicity, the symbol f n will be hereafter employed to denote f (u n ).

The reader is referred to Sec. 1.1.4 for an explanation of the training procedure terminology employed here.

The reader is referred to Sec.

2.3.2 for a description of the vortex and of the double shear layer initial conditions used here.

In JAX, array mutation (an in-place operation) is guaranteed to be performed whenever inside a JIT-compiled function. However, JAX' sharp bits (JAX -The Sharp Bits -JAX

documentation 2022) may turn the coding experience as a whole quite cumbersome.

This technique is not detailed here; the interested reader is referred to Sec. 5.5 of[START_REF] Bartholomew-Biggs | Automatic differentiation of algorithms[END_REF] for a description.

Available at https://github.com/google/jax-cfd

Available at https://github.com/tumaer/jaxfluids

The suffix b stands for bar (• ) here, the symbol that distinguishes a reverse derivative from its primal in the text.

Periodic conditions apply.

The reader can find the adjoint programs in https://github.com/luciano-drozda/ reassignment-vs-mutation/.

Full implementations are made available in https://github.com/luciano-drozda/ reassignment-vs-mutation/.

List of Tables

Conclusions and perspectives

In this work, a new spectral analysis for numerical schemes is introduced based on a local spatio-temporal discretization setting. The Global Spectral Analysis (GSA) (see Sec. 3.2.1) proposed by T. [START_REF] Sengupta | Analysis of central and upwind compact schemes[END_REF] is extended by interpreting the mesh as an electrical circuit whose transfer function is represented by the local numerical amplification factor. Within such an analogy, each mesh element acts as an impedance block providing resistance to the solution propagation. The resulting framework allows to quantify the local spectral properties of a given scheme, applicable to irregular meshes and non-linear systems. It is named Local Transferfunction Analysis (LTA). Taylor series or polynomial order of analysis is only applicable where the order of error makes sense (smooth or well resolved regions) (Brooks and T. J. [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF]. But practical Large-Eddy Simulations (LES) typically employ coarser, under-resolved meshes that target the capturing of the larger flow scales. Therefore, the need to design CFD schemes for bandwidth similar to Roe (2021) is emphasized. LTA being a spectral approach, it uses bandwidth as the basis of measuring accuracy and provides a complete overview of the error behavior even near high gradients and shocks. It provides practical guidance to mitigate errors that might otherwise trigger numerical instability or compromise solution accuracy.

Using the LTA framework, the optimal design of numerical methods can be formulated as a simple impedance matching problem. It was shown that the time-domain counterpart of such a problem corresponds to setting a balance between Total Variation Diminishing (TVD) and minimal norm of error conditions. Following this finding, an optimization problem based on impedance matching to redesign and tune the TTGC scheme of [START_REF] Rudgyard | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF] was set out. Differentiable programming methods to compute gradients required by the optimizer and find local optima for the free parameter of the numerical scheme were employed.

Appendix A

Example of Tapenade-generated adjoint program

Below is a code snippet illustrating the Fortran subroutine FUNC_B as the Tapenade-generated adjoint of the subroutine FUNC which computes

The mutation problem

Chapter 1 introduced ML-related techniques that are employed in the developments proper to data-driven discretizations. Among them, Algorithmic Differentiation (AD) computes the gradients required to train ML models. The main AD frameworks used in ML pipelines, namely, TensorFlow and PyTorch, have found great adoption for their relative ease of use and scaling capabilities. These frameworks fail to deliver performant differentiation for generic control flow proper to coupled ML-PDE systems, though. This is due to their lack of support to array mutation, i.e., modifying the contents of an array after it has been created. This appendix is dedicated to quantifying the impact of this limitation in chains of operations one can find in PDE solvers for CFD. Building differentiable solvers for unstructured CFD in frameworks that lack support to array mutation is prohibitively expensive. To the extent of the author's knowledge, a quantitative study of this problem has not yet been reported in the literature.

B.1 Presentation of the problem

Guido Van Rossum's efforts in late 1980s to create an easily-extensible, highlevel scripting language (artima -The Making of Python 2022) helped turning Python into the most popular programming language as of 2022 (according to the TIOBE Index (2022)). Despite suffering from lack of native support to performant numerical computing, its ease of use has motivated developers Apparent workarounds are provided by certain specific methods, like tf.tensor_scatter_nd_add in TensorFlow, .scatter_add_ in PyTorch, and .at[] in JAX. None of them perform actual mutation of the original array, though 1 . They rather execute re-assignements by creating new arrays of same size. This is illustrated in Fig. B.1. Such a limitation is justified by the fact that reverse mode AD in its simplest form requires all intermediate variables of the forward pass to be stored for the execution of the backward pass. This is what has probably led TensorFlow, PyTorch and, to some extent, JAX to not support mutation from the beginning of their developments:

" Allowing mutation of variables in-place makes program analysis and transformation difficult. JAX requires that programs are

Appendix C Local wavenumber analysis and shock sensor of Jameson

The local wavenumber k is defined as the spatial rate of change of the phase angle of the convecting wave u(x, t) [START_REF] Bracewell | The Fourier Transform and its Applications[END_REF], i.e.,

This quantity is extensively used in the field of magnetism and geology to identify physical sizes of objects like rock formations under the earth by measuring the aberrations in the phase angle of the spatial magnetic or gravitational field [START_REF] Pilkington | The relationship between local wavenumber and analytic signal in magnetic interpretation[END_REF]. The phase angle variation is obtained using the Eikonal E = |∇u|, which represents the surface of constant phase angle. Note that the operator | • | is the Euclidean norm. Thus one can rewrite the local wavenumber using the Eikonal as

The nodal junction i connecting two finite-elements e L [i -1, i] and e R [i, i + 1] is considered. The local gradient at the node junction i can be approximately evaluated using finite-difference formula as

where ∆u i+1/2 = u i+1 -u i and ∆u i-1/2 = u i -u i-1 . If one approximates E at the junction as the average of the values at e L and e R , i.e., E i = determine if the variations are within the resolvable limit of the numerical scheme. When variations are beyond the numerical resolution (unresolved) they can be removed using artificial dissipation. Using the uncertainty principle, one can show that a function u(x) and its Fourier transform Û (k) cannot be supported on arbitrarily small sets as given by the Heisenberg-Pauli-Weyl inequality:

where ∥u∥ 2 is the l2 norm of the function u. Therefore, these spatial estimates are at best approximations to the local wavenumber; when the spatial variations are sharp, the uncertainty in wavenumber grows quite rapidly.

Appendix D Entropy condition and impedance matching

Non-linear convection problems like the Burgers' equation and compressible Euler equations in gas dynamics admit discontinuous solutions with jumps called shock waves. Therefore, the solution that yields the shock wave by enforcing the entropy condition is chosen. There are several ways to formulate this entropy condition, but in this work, the Rankine-Hugoniot (RH) relation is used to enforce the entropy condition. The non-linear conservation law is defined as

Then the shock speed ξ is given by the RH condition as

where u 1 and u 2 are the states before and after the shock (as shown in . With the addition of TVD condition, one can show that the five nodes (i -2, i -1, i, i + 1, and i + 2) influence the impedance matching problem for shock capturing.