*2.3. The Discontinuous Galerkin Discretization*

In this section, we briefly introduce the semi-discrete form of the Discontinuous Galerkin finite element method (DG) for compressible inviscid flows. The compressible Euler equations were derived from the Navier–Stokes equations by neglecting diffusive terms. They still provide a model for the conservation of mass, momentum, and energy in the fluid and can be described in vectorial notation by:

$$
\partial\_t \mathbf{u} + \nabla \cdot \mathbf{F}(\mathbf{u}) = 0,\tag{12}
$$

equipped with suitable initial and boundary conditions. Here, **u** is a vector of the conservative variables, and the flux function **F**(**u**)=(**f**(**u**), **g**(**u**))*<sup>T</sup>* for two spatial dimensions is given by:

$$\mathbf{u} = \begin{bmatrix} \rho \\ \rho u \\ \rho v \\ \rho E \end{bmatrix}, \mathbf{f}(\mathbf{u}) = \begin{bmatrix} \rho u \\ \rho u^2 + p \\ \rho uv \\ (\rho E + p)u \end{bmatrix}, \mathbf{g}(\mathbf{u}) = \begin{bmatrix} \rho v \\ \rho uv \\ \rho v^2 + p \\ (\rho E + p)v \end{bmatrix}.$$

where *ρ*, **v** = (*u*, *v*)*T*, *E*, *p* denotes the density, velocity vector, specific total energy, and pressure, respectively. The system is closed by the equation of state assuming the fluid obeys the ideal gas law with pressure defined as *p* = (*γ* − 1)*ρ <sup>e</sup>* <sup>−</sup> <sup>1</sup> <sup>2</sup> (*u*<sup>2</sup> + *<sup>v</sup>*2) ; where *γ* = *cp cv* is the ratio of specific heat capacities and *e* is the total internal energy per unit mass.

The discontinuous Galerkin formulation of the above equation was obtained by multiplying it with a test function *ψ* and integrating it over the domain Ω. Thereafter, integration by parts was used to obtain the following weak formulation:

$$
\int\_{\Omega} \psi \frac{\partial \mathbf{u}}{\partial t} d\Omega + \oint\_{\partial \Omega} \psi \mathbf{F}(\mathbf{u}) \cdot \mathbf{n} ds - \int\_{\Omega} \nabla \psi \cdot \mathbf{F}(\mathbf{u}) d\Omega = 0, \forall \psi,\tag{13}
$$

where *ds* denotes the surface integral. A discrete analogue of the above equation was obtained by considering a tessellation of the domain Ω into *n* closed, non-overlapping elements given by *<sup>T</sup>* = {Ω*i*|*<sup>i</sup>* = 1, 2, ... , *<sup>n</sup>*}, such that <sup>Ω</sup> = ∪*<sup>n</sup> <sup>i</sup>*=1Ω*<sup>i</sup>* and Ω*<sup>i</sup>* ∩ Ω*<sup>j</sup>* = ∅∀*i* = *j*. We define a finite element space consisting of discontinuous polynomial functions of degree *m* ≥ 0 given by:

$$P^{m} = \{ f \in [L^{2}(\Omega)]^{m} : f\_{|\Omega\_{k}} \in \mathbb{P}^{m}(\Omega\_{k}) \forall \Omega\_{k} \in \Omega \} \tag{14}$$

where P*m*(Ω*k*) is the space of polynomials with largest degree *<sup>m</sup>* on element <sup>Ω</sup>*k*. With the above definition, we can write the approximate solution **u***h*(**x**, *t*) within each element using a polynomial of degree *m*:

$$\mathbf{u}\_h(\mathbf{x}, t) = \sum\_{i=1}^m \vartheta\_i \phi\_i, \quad \psi\_h(\mathbf{x}) = \sum\_{i=1}^m \vartheta\_i \phi\_i. \tag{15}$$

where the expansion coefficients *u*ˆ*<sup>i</sup>* and *v*ˆ*<sup>i</sup>* denote the degrees of freedom of the numerical solution and the test function, respectively. Notice that there is no global continuity requirement for **u***<sup>h</sup>* and *ψ<sup>h</sup>* in the previous definition. Splitting the integrals in Equation (13) into a sum of integrals over elements Ω*i*, we obtain the space-discrete variational formulation:

$$\sum\_{i=1}^{n} \frac{\partial}{\partial t} \int\_{\Omega\_{i}} \psi\_{\hbar} \mathbf{u}\_{\hbar} d\Omega + \oint\_{\partial \Omega\_{i}} \psi\_{\hbar} \mathbf{F}(\mathbf{u}\_{\hbar}) \cdot \mathbf{n} ds - \int\_{\Omega\_{i}} \nabla \psi\_{\hbar} \cdot \mathbf{F}(\mathbf{u}\_{\hbar}) d\Omega = 0, \forall \psi\_{\hbar} \in P^{n}. \tag{16}$$

Due to the element local support of the numerical representation, the flux term is not uniquely defined at the element interfaces. The flux function is, therefore, replaced by a numerical flux function **F**∗(**u**− *<sup>h</sup>* , **<sup>u</sup>**<sup>+</sup> *<sup>h</sup>* , **n**), where **u**<sup>−</sup> *<sup>h</sup>* and **<sup>u</sup>**<sup>+</sup> *<sup>h</sup>* are the interior and exterior traces at the element face in the direction **n** normal to the interface. A choice of appropriate numerical flux can then be selected from several numerical flux schemes. For our simulations, we used the Lax–Friedrichs scheme for numerical flux.

For simplicity, we can re-write the equation above in matrix vector notation and obtain:

$$\frac{\partial}{\partial t}\hat{\mathbf{u}} = M^{-1}\left(\mathbb{S}\cdot\mathbb{F}(\hat{\mathbf{u}}) - M^{\mathrm{F}}\cdot\mathbb{F}(\hat{\mathbf{u}})\right) =: rhs(\hat{\mathbf{u}}).\tag{17}$$

where *M*, *S* denote the mass and the stiffness matrices and *M<sup>F</sup>* are the so-called face mass lumping matrices. The above obtained ordinary differential Equation (17) can be solved in time using any standard time-stepping method, e.g., a Runge–Kutta method.

In our implementation, we exploited the fact that we only used cubical elements. This choice of simple elements allowed for a tensor-product notation in the multi-dimensional basis functions. The symmetry of the elements enabled efficient dimension-by-dimension algorithms in the computation.
