# A high-order discontinuous Galerkin Method using a mixture of Gauss-Legendre and Gauss-Lobatto quadratures for improved efficiency

Stéphanie Chaillat, Régis Cottereau, Ruben Sevilla

<figure><img src="image_1.png" /></figure>

## To cite this version:

Stéphanie Chaillat, Régis Cottereau, Ruben Sevilla. A high-order discontinuous Galerkin Method using a mixture of Gauss-Legendre and Gauss-Lobatto quadratures for improved efficiency. 2022. [hal-03695573](https://hal.science/hal-03695573)

HAL Id: hal-03695573

https://hal.science/hal-03695573

Preprint submitted on 15 Jun 2022

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

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

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A high-order discontinuous Galerkin Method using a mixture of Gauss-Legendre and Gauss-Lobatto quadratures for improved efficiency

S. Chaillat<fnref n="1" />, R. Cottereau<fnref n="2" />, and R. Sevilla<fnref n="3" />

In discontinuous Galerkin spectral element methods (DGSEM), the two most common approaches to numerically integrate the terms of the weak form are either using Gauss-Legendre or Gauss-Lobatto quadratures. The former yields more accurate results but at a higher computational cost, so that a priori it is not clear whether one approach is more efficient than the other. In this paper, it is shown (theoretically for a particular case and numerically for the general case) that using Gauss-Lobatto quadrature for the convection matrix actually introduces a negligible error. In contrast, using Gauss-Lobatto quadratures for the evaluation of the jump term in the element faces introduces a sizeable error. This leads to the proposal of a new DG approach, where the convection matrix is evaluated using Gauss-Lobatto quadratures, whereas the face mass matrices are integrated using Gauss-Legendre quadratures. For elements with constant Jacobian and constant coefficients, a formal proof shows that no numerical integration error is actually introduced in the evaluation of the residual, even though both the mass and the convection matrices are not computed exactly with Gauss-Lobatto quadratures. For elements with non-constant Jacobian and/or non-constant coefficients, the impact of numerical integration error on the overall error is evaluated through a series of numerical tests, showing that this is also negligible. In addition, the computational cost associated to the matrix-vector products required to evaluate the residual is evaluated precisely for the different cases considered. The proposed approach is particularly attractive in the most general case, since the use of Gauss-Lobatto quadratures significantly speeds-up the evaluation of the residual.

Keywords discontinuous Galerkin method, spectral element method, high-order methods, wave propagation

## 1 Introduction

High-order discontinuous Galerkin (DG) methods (Bey et al. 1996; Hu et al. 1999; Warburton et al. 1999; Rasetarinera et al. 2000; Sherwin 2000; Houston et al. 2002; Ainsworth 2004a; Cockburn et al. 2005) have proven very efficient to approximate the solution of wave propagation over large complex media. They come in different flavours, depending in particular on the choice of quadrature for integration, Gauss-Legendre (GLegQ) or Gauss-Lobatto (GLobQ), on whether the interpolation and integrations points are the same or not, on whether the weak formulation is integrated by parts once or twice, and on the choice of flux that weakly enforces the continuity between elements.

When considering collocation approaches, using the nodes of either the GLegQ or the GLobQ, the mass matrix becomes diagonal and very efficient strategies arise in combination with time marching algorithms and fast matrix-vector products in the evaluation of the residual (Deville et al. 2002) (see also Section 3.2). At first, using the GLegQ seems like a better option in terms of accuracy because that quadrature can integrate exactly higher order polynomials. However, the fluxes then need to be interpolated on the faces (Canuto et al. 2007), which induces an additional cost at each time step to compute the residual (evaluated at 15% of the total simulation cost in Kopriva et al. (2010)). These competing effects make it difficult to know which method is better in terms of accuracy and cost, although most authors (Castel et al. 2009; Kopriva et al. 2010; Gassner et al. 2011) seem to (slightly) prefer GLegQ over GLobQ, in particular with arguments of stability (Castel et al. 2009; Gassner et al. 2011).

<footnote n="1">¹ Laboratoire POEMS, CNRS-ENSTA-INRIA, IP Paris, France.</footnote>

<footnote n="2">² Aix-Marseille University, CNRS, Centrale Marseille, LMA UMR 7031, Marseille, France.</footnote>

<footnote n="3">³ Zienkiewicz Centre for Computational Engineering, Faculty of Science and Engineering, Swansea University, Bay Campus, Fabian Way, Crymlyn Burrows, SA1 8EN, Swansea, Wales, UK.</footnote>
Another flavour of the high-order DG method is the so-called nodal version (Giraldo et al. 2002; Hesthaven et al. 2008b; Giraldo et al. 2008), where the interpolation points are the nodes of the GLobQ, but integration follows the GLegQ. In that version, a good accuracy is expected because integrals are accurately evaluated, but the cost of computing the residual at each time step is potentially high because the mass matrix is not diagonal. Here again, the combination of these advantages and defects makes it very difficult to draw a definite conclusion on the interest of this version in terms of accuracy and cost. In this paper, a comparison and an analysis is provided of the nodal and the GLobQ versions discussed above, where the GLobQ is used for both integration and interpolation. These two versions will be referred as the DG method with GLegQ and GLobQ respectively because the only difference lies in the quadrature formula used for numerical integration. The interpolation polynomials are the same in both approaches and defined using the nodes of the GLobQ.

The first objective of this paper consists of evaluating precisely the complexity of both DG methods in terms of operation count required to perform the matrix-vector products required when employing an explicit time marching scheme. Additionally, accuracy of both methods will be evaluated by using a series of numerical simulations of increasing complexity. The second objective of the paper is to propose a novel version of the high-order DG method, bridging the gap between the DG methods with GLegQ and GLobQ, where the convection part of the residual is integrated with the GLobQ, while the face jump term is evaluated with the GLegQ. It will be shown that this simple change allows to improve greatly the accuracy of the method, while retaining its numerical efficiency. In particular, we show analytically that in the case of a constant Jacobian, the use of GLobQ to integrate the convection part of the residual induces absolutely no error with respect to the GLegQ. In the non-constant Jacobian case, a series of numerical examples illustrate that this property is approximately conserved.

The remainder of the paper is organised as follows. Section 2 introduces the classical DG formulation for a linear conservation law, its computational complexity in the most general case for tensor-product elements and a brief discussion on the different choices to perform the numerical integration. Note that the discussion throughout the paper is limited to tensorial elements (quadrilaterals and hexaedra). Section 3 presents alternative formulations of the residual for the general case and for the particular case when constant coefficients and meshes with affine elements (i.e., constant Jacobian of the isoparametric mapping) are considered. It also provides the formulation of the residual for the DG method with GLobQ. In all cases the cost required to compute the matrix-vector products required during the explicit time marching process is detailed. In Section 4, the new DG method with mixed quadratures is proposed. First, a formal proof is provided showing that, under certain hypotheses, the contribution to the residual given by the convection term is identical if GLegQ or GLobQ are used. This motivates the introduction of the new method that bridges the gap between the DG methods with GLegQ and GLobQ. Section 5 shows a series of numerical examples to evaluate the accuracy of the DG approaches considered in this work. The examples involve problems with constant and non-constant velocity fields both in two and three dimensions. In addition, the use of Cartesian and non-Cartesian grids is also studied, to evaluate the influence of having elements with non-constant Jacobian of the isoparametric mapping. The discussion also considers the computational cost of the DG approaches considered. Finally, Section 6 presents the conclusions of the work.

## 2 DG formulation for conservation laws

This section introduces the DG formulation of a scalar conservation law and the spatial and temporal discretisations. The particular choice considered here of tensor product elements is recalled and the cost of evaluating the residual of the semi-discrete system is presented. Finally the different choices for the numerical integration of the terms appearing in the weak form are briefly discussed.
## 2.1 Model problem

Let us consider the general linear conservation law in a open bounded domain $\Omega \subset \mathbb{R}^{nsd}$ with boundary $\partial\Omega$, where $n_{sd}$ denotes the number of spatial dimensions,

$$ \frac{\partial u}{\partial t} + \nabla \cdot f(\mathbf{x}, u) = 0, \quad \text{on } \Omega \times (0, T]. \qquad (1) $$

where $u(\mathbf{x}, t)$ is the conserved variable, $f(\mathbf{x}, u)$ is the hyperbolic flux vector and $T$ is the final time.

It is often convenient to write the hyperbolic flux vector as

$$ f(\mathbf{x}, u) = a(\mathbf{x})u, \qquad (2) $$

where the velocity field vector $a = (a_1, \dots, a_{n_{sd}}) \in [\mathcal{L}^{\infty}(\Omega)]^{n_{sd}}$ is assumed to be divergence free.

The boundary of the domain is assumed to be partitioned as $\partial\Omega = \bar{\Gamma}_{in} \cup \bar{\Gamma}_{out}$ where $\Gamma_{in}$ and $\Gamma_{out}$ are the inflow and outflow parts of the boundary respectively, defined as

$$ \Gamma_{in} = \{\mathbf{x} \in \partial\Omega \mid a_n(\mathbf{x}) < 0\}, \quad \Gamma_{out} = \{\mathbf{x} \in \partial\Omega \mid a_n(\mathbf{x}) > 0\}, \qquad (3) $$

where $a_n(\mathbf{x}) = a(\mathbf{x}) \cdot n(\mathbf{x})$ is the normal component of the velocity vector and $n(\mathbf{x})$ is the outward unit normal vector at $\mathbf{x} \in \partial\Omega$.

As usual for hyperbolic problems, boundary conditions can only be imposed on the inflow part of the boundary (LeVeque 2002). Here a Dirichlet boundary condition is considered, namely

$$ u(\mathbf{x}, t) = u_D(\mathbf{x}, t) \quad \text{on } \Gamma_{in} \times (0, T]. \qquad (4) $$

Finally, the strong form of the problem must be completed with an initial condition

$$ u(\mathbf{x}, t_0) = u_0(\mathbf{x}) \quad \text{in } \Omega. \qquad (5) $$

## 2.2 DG weak formulation

The domain is partitioned in $n_{e1}$ disjoint elements $\Omega_e$ such that

$$ \bar{\Omega} = \bigcup_{e=1}^{n_{e1}} \bar{\Omega}_e. \qquad (6) $$

The weak formulation for a generic element $\Omega_e$ is obtained after multiplying Eq. (1) by a test function, $v$ and integrating in $\Omega_e$

$$ \int_{\Omega_e} v \frac{\partial u}{\partial t} d\Omega - \int_{\Omega_e} \nabla v \cdot f(\mathbf{x}, u) d\Omega + \int_{\partial \Omega_e} v \hat{f}_n(\mathbf{x}, u, u^{\text{out}}) d\Gamma = 0, \qquad (7) $$

where the integration by parts has been already performed on the divergence term and the boundary term features the numerical normal flux, $\hat{f}_n$.

A natural choice to define the numerical flux, for the linear hyperbolic equation considered here, is to employ a flux splitting technique (Hesthaven et al. 2008a). First, the physical normal flux $f_n(\mathbf{x}, u) = f(\mathbf{x}, u) \cdot n(\mathbf{x})$ is split as

$$ f_n(\mathbf{x}, u) = f_n^-(\mathbf{x}, u) + f_n^+(\mathbf{x}, u), \qquad (8) $$

where

$$ f_n^-(\mathbf{x}, u) = a_n^-(\mathbf{x})u, \quad f_n^+(\mathbf{x}, u) = a_n^+(\mathbf{x})u, \qquad (9) $$

are the inflow and outflow normal fluxes, respectively, and

$$ a_n^+(\mathbf{x}) := \frac{1}{2}(a_n(\mathbf{x}) + |a_n(\mathbf{x})|), \quad a_n^-(\mathbf{x}) := \frac{1}{2}(a_n(\mathbf{x}) - |a_n(\mathbf{x})|). \qquad (10) $$
The numerical normal flux, evaluated in terms of the trace of the solution on element $\Omega_e$ and the trace of the solution $u^{\text{out}}$ on the neighbouring element, is then defined as

$$ \hat{f}_n(\mathbf{x}, u, u^{\text{out}}) = f_n^+(\mathbf{x}, u) + f_n^-(\mathbf{x}, u^{\text{out}}) = a_n^+(\mathbf{x})u + a_n^-(\mathbf{x})u^{\text{out}}. \quad (11) $$

The DG weak formulation of Eq. (1) with flux splitting can be written as

$$ \int_{\Omega_e} \nu \frac{\partial u}{\partial t} d\Omega + \int_{\Omega_e} \nu (a(\mathbf{x}) \cdot \nabla u) d\Omega + \int_{\partial \Omega_e} \nu a_n^-(\mathbf{x}) \llbracket u \rrbracket d\Gamma = 0. \quad (12) $$

where a second integration by parts has been performed after introducing the numerical normal flux and the operator $\llbracket u \rrbracket := u^{\text{out}} - u$ denotes the jump of the trace of the solution on the element boundary.

## 2.3 Spatial and temporal discretisations

The semi-discrete formulation is obtained by approximating the solution $u$ as

$$ u(\xi, t) \approx u^h(\xi, t) = \sum_{J=1}^{n_{\text{en}}} u_J(t) N_J(\xi) \quad (13) $$

on a reference element, with local coordinates $\xi = (\xi_1, \dots, \xi_{n_{\text{sd}}})$. Here $u_J(t)$ denotes the value of the solution at node $x_J$ at time $t$, $N_J$ is the shape function associated with the $J$-th node and $n_{\text{en}}$ is the total number of nodes on the element. The shape functions considered here are standard Lagrange polynomials of order $p$, in each direction, defined on the nodes of the GLobQ.

The isoparametric mapping is employed to link the reference element $\hat{\Omega}$ and a physical element $\Omega_e$, namely

$$ \varphi : \hat{\Omega} \subset \mathbb{R}^{n_{\text{sd}}} \rightarrow \Omega_e \subset \mathbb{R}^{n_{\text{sd}}} \\ \xi \mapsto \varphi(\xi) := \sum_{J=1}^{n_{\text{en}}} x_J N_J(\xi), \quad (14) $$

where $\{x_J\}_{J=1,\dots,n_{\text{en}}}$ are the nodal coordinates of the physical element $\Omega_e$ (Zienkiewicz et al. 2000).

Similarly, an isoparametric mapping is used to link the reference face $\hat{\Gamma}$ and a physical face $\Gamma_e^f \subset \partial \Omega_e$

$$ \psi : \hat{\Gamma} \subset \mathbb{R}^{n_{\text{sd}}-1} \rightarrow \Gamma_e^f \subset \mathbb{R}^{n_{\text{sd}}} \\ \eta \mapsto \psi(\eta) := \sum_{J=1}^{n_{\text{fn}}} x_J^f N_J^f(\eta), \quad (15) $$

where $\{x_J^f\}_{J=1,\dots,n_{\text{fn}}}$ are the nodal coordinates of the subset of nodes that belong to the physical face $\Gamma_e^f$ and $N_J^f$ is the shape function associated to $J$-th face node, defined in local coordinates $\eta = (\eta_1, \dots, \eta_{n_{\text{sd}}-1})$.

Introducing the approximation of the solution in the weak formulation of Eq. (12) and selecting the space of the weighting functions to be the same as the space of the interpolation functions, leads to the following system of ordinary differential equations (ODEs)

$$ M \frac{dU}{dt} + CU + \sum_{f=1}^{n_{\text{fa}}} M^f \llbracket U^f \rrbracket = 0, \quad (16) $$

where $n_{\text{fa}}^e$ is the number of faces of the element $\Omega_e$. It is worth noting that $U^f$ is obtained by restricting $U$ to the degrees of freedom associated to the nodes of the faces of $\Omega_e$.

The elemental mass and convection matrices are given, respectively, by

$$ M_{IJ} = \int_{\hat{\Omega}} N_I N_J |J| d\Omega, \quad C_{IJ} = \sum_{\ell=1}^{n_{\text{sd}}} \int_{\hat{\Omega}} N_I G_{\ell} \frac{\partial N_J}{\partial \xi_{\ell}} d\Omega, \quad (17) $$
for $I, J = 1, \dots, n_{en}$, where

$$G_{\ell}(\xi) = |J(\xi)| \sum_{k=1}^{n_{sd}} a_k(\varphi(\xi)) J_{kl}^{-1}(\xi) = \sum_{k=1}^{n_{sd}} a_k(\varphi(\xi)) H_{kl}(\xi), \quad (18)$$

for $\ell = 1, \dots, n_{sd}$, $J = \partial\varphi/\partial\xi$ denotes the Jacobian of the isoparametric mapping and $H = \text{adj}(J)$ is the adjoint of $J$.

Similarly, the face mass matrix is given by

$$M_{IJ}^{f} = \begin{cases} \int_{\hat{\Gamma}} a_{n}^{-} N_{I}^{f} N_{J}^{f} \|J^{f}\| d\Gamma, & \text{if } \mathbf{x}_{I} \in \Gamma_{e}^{f} \\ 0 & \text{otherwise,} \end{cases} \quad (19)$$

for $I = 1, \dots, n_{en}$ and $J = 1, \dots, n_{fn}$, where $J^f$ is the Jacobian of the restriction of the isoparametric mapping to an element face. In two dimensions $J^f = \partial\psi/\partial\eta_1$ and in three dimensions $J^f = \partial\psi/\partial\eta_1 \times \partial\psi/\partial\eta_2$.

**Remark 1** The size of the mass and convection matrices is $n_{en} \times n_{en}$, whereas the face mass matrix $M^f$ is defined as a rectangular matrix of size $n_{en} \times n_{fn}$ only to ensure consistency in Eq. (16). Therefore, when analysing the cost of the different approaches presented in this work (in Sections 2.5, 3.1.3 and 3.2.3), non necessary operations are not accounted for.

When the system of ODEs given by Eq. (16) is advanced in time using an explicit time marching algorithm, the scheme requires the solution of a linear system of equations element by element at each time step, avoiding the assembly and solution of a large sparse linear system in each time step. Here a classical fourth order explicit Runge-Kutta scheme is considered, namely

$$U^{n+1} = U^n + \frac{\Delta t}{6} (R^{(1)} + 2R^{(2)} + 2R^{(3)} + R^{(4)}). \quad (20)$$

The four stages are given by

$$\begin{align*} R^{(1)} &= R(U^n, t^n) \\ R^{(2)} &= R(U^n + R^{(1)}\Delta t/2, t^n + \Delta t/2) \\ R^{(3)} &= R(U^n + R^{(2)}\Delta t/2, t^n + \Delta t/2) \\ R^{(4)} &= R(U^n + R^{(3)}\Delta t, t^n + \Delta t). \end{align*} \quad (21)$$

where $R$ denotes the residual of the system of ODEs

$$R := - \left\{ C U + \sum_{f=1}^{n_{fa}^{e}} M^{f} [U^{f}] \right\}. \quad (22)$$

where $C := M^{-1}C$ and $M^f := M^{-1}M^f$.

This time marching algorithm is known to be conditionally stable, so the time step must be selected below the stability limit (Hairer et al. 2006). In all the numerical examples the time step is chosen small enough so that the error is dominated by the spatial discretisation. However, as the stability can influence the overall complexity by changing the number of time steps necessary to reach a final time, some comments and numerical tests will be provided in Section 5.

**Remark 2** Although the discussion presented in this work is based on this particular fourth order explicit Runge-Kutta scheme, the conclusions of the paper apply to any explicit time integrator.
## 2.4 Tensorial basis and tensor products

This work focuses on tensor-product elements, namely quadrilateral and hexahedral elements in two and three dimensions, respectively. The reference element is taken as $\hat{\Omega} = [-1, 1]^{nsd}$ and the reference face is $\hat{\Gamma} = [-1, 1]^{nsd-1}$. Assuming for simplicity that the interpolation polynomials have the same order $p$ in each direction, the number of element nodes is $n_{en} = (p+1)^{nsd}$ and the number of face nodes are $n_{fn} = (p+1)^{nsd-1}$, for each face. Throughout the document, lowercase (both for matrices and indices) denotes one dimensional quantities whereas uppercase is used for multi-dimensional quantities.

In this scenario, the shape functions {$N_j$}$_{j=1,...,n_{en}}$ are defined as a product of one-dimensional shape functions, denoted by {$n_j$}$_{j=1,...,p+1}$. Using the lexicographic order, multi-dimensional shape functions are defined as

$$N_J(\xi) = \prod_{k=1}^{n_{sd}} n_{jk}(\xi_k), \quad \text{with} \quad J = 1 + \sum_{k=1}^{n_{sd}} (j_k - 1)(p+1)^{k-1}. \qquad (23)$$

## 2.5 Operations count to evaluate the residual

Once the elemental matrix $C$, of size $n_{en} \times n_{en}$ is computed for an element, the evaluation of the term of the residual in Eq. (22) associated to an element requires $n_{op}^e = (2n_{en} - 1)n_{en}$ operations to perform the matrix-vector product. The leading term in this evaluation is thus of the order of $O(2n_{en}^2)$.

Similarly, once the matrix $M^f$ is computed for a face, the evaluation of the term of the residual in Eq. (22) associated to an internal face reduces to the evaluation of a matrix-vector product. Since each $M^f$ is a matrix of size $n_{en} \times n_{fn}$, the evaluation requires a total of $n_{op}^f = (2n_{fn} - 1)n_{en}$ operations. The leading term in this evaluation is of the order of $O(2n_{fn}n_{en})$.

The overall cost of evaluating the residual, once the matrices $M^f$ and $C$ are computed, is $n_{op} = n_{e1}(2n_{en} - 1)n_{en} + n_{fa}(2n_{fn} - 1)n_{en}$, where $n_{e1}$ is the number of mesh elements and $n_{fa}$ is the total number of mesh edges or faces in two or three dimensions respectively.

For a regular mesh of tensor-product elements, i.e. quadrilateral or hexahedral elements, in two and three dimensions respectively, the number of mesh edges/faces can be written in terms of the total number of mesh elements, namely $n_{fa} = n_{sd}n_{e1}$, by assuming that the mesh is large enough so that the number of exterior faces can be neglected. This leads to the following expression for the total number of operations, only in terms of the number of elements, $n_{e1}$, and the degree of approximation, $p$,

$$n_{op} = n_{e1}(p+1)^{nsd} [2(p + n_{sd} + 1)(p+1)^{nsd-1} - n_{sd} - 1], \qquad (24)$$

with leading term $O(2n_{e1}p^{2nsd})$. As the number of elements can be factorised, in the remainder of the paper the total number of operations per element will be considered. The number of operations required to compute the matrix-vector products in Eq. (22) for a tensor-product element are given in Table 1 for future reference.

<table><caption>**Table 1:** Number of operations required (per element) to compute the matrix-vector products in Eq. (22) using the GLegQ for numerical integration.</caption><thead><tr><th></th><th>Operations</th><th>Leading term</th></tr></thead><tbody><tr><td>$n_{op}^e$</td><td>$[2(p+1)^{nsd} - 1](p+1)^{nsd}$</td><td>$O(2p^{2nsd})$</td></tr><tr><td>$n_{op}^f$</td><td>$[2(p+1)^{nsd-1} - 1](p+1)^{nsd}$</td><td>$O(2p^{2nsd-1})$</td></tr></tbody></table>

**Remark 3** This approach requires to store the matrices $C$ and $M^f$. The cost to assemble these matrices and the memory requirements are thus both of the order of $O(n_{e1}p^{2nsd})$. Depending on the available resources, it can be interesting to save these matrices but that does not change the main complexity of the algorithm.
## 2.6 Choices of quadrature for integration

The computation of the elemental mass and convection matrices M and C in Eq. (17) is performed using a numerical quadrature defined on the reference element $\hat{\Omega}$, with $n_{ip}$ points $\{\xi_G\}_{G=1,...,n_{ip}}$ and weights $\{\omega_G\}_{G=1,...,n_{ip}}$, namely

$$ 
\begin{aligned}
M_{IJ} &\approx \sum_{G=1}^{n_{ip}} N_I(\xi_G) N_J(\xi_G) |J(\xi_G)| \omega_G, \\
C_{IJ} &\approx \sum_{G=1}^{n_{ip}} N_I(\xi_G) \omega_G \sum_{\ell=1}^{n_{sd}} G_{\ell}(\xi_G) \frac{\partial N_J(\xi_G)}{\partial \xi_{\ell}}.
\end{aligned} 
\quad (25) $$

Similarly the face mass matrix $M^f$ in Eq. (19) is computed by using a numerical quadrature defined on the reference face $\hat{\Gamma}$, with $n_{ip}^f$ points $\{\eta_G\}_{G=1,...,n_{ip}^f}$ and weights $\{\partial_G\}_{G=1,...,n_{ip}^f}$, namely

$$ M_{IJ}^f \approx \begin{cases} \sum_{G=1}^{n_{ip}^f} a_n^-(\psi(\eta_G)) N_I(\eta_G) N_J(\eta_G) \|J^f(\eta_G)\| \partial_G & \text{if } \mathbf{x}_I \in \Gamma_e^f \\ 0 & \text{otherwise,} \end{cases} \quad (26) $$

In this work two numerical quadratures, namely Gauss-Legendre (GLegQ) and Gauss-Lobatto (GLobQ) quadratures, are considered and compared. For integration over the interval $[-1, 1]$, GLegQ are based on the roots of the Legendre polynomials, of order $p+1$, and they integrate exactly polynomials of order $2p+1$. In contrast, GLobQ are based on the extremes of the segment and the roots of the derivative of the Legendre polynomials, of order $p-1$, (for a total of $p+1$ nodes) and integrates exactly polynomials of order $2p-1$.

In general, i.e. even when the Jacobian is not constant, it is possible to select the number of integration points to integrate exactly the terms of the elemental mass matrix. For an approximation with polynomials of degree $p$, the entries of the mass matrix require the integration of a polynomial of degree $(n_{sd} + 2)p - 1$, which can be computed exactly with a GLegQ with $2p-1$ points. The convection matrix can be computed exactly only if the convection velocity is a polynomial. For instance, if the convection velocity is a polynomial of degree $q$, the entries of the mass matrix require the integration of a polynomial of degree $(n_{sd} + q + 1)p - 1$. Finally, the face mass matrix cannot be computed exactly in general due to the non-polynomial nature of the norm of $J^f$ (Sevilla et al. 2011).

## 3 Alternative formulations of the residual

This section presents an efficient strategy to evaluate the residual of the ODE system of Eq. (22) for tensor-product elements. First, the particular case involving a constant velocity field and meshes with elements with constant Jacobian is considered. The more general case, where the velocity is not constant and/or the Jacobian of the isoparametric mapping is not constant is later discussed.

The main property used in the following is that the tensor product $C = a \otimes b$ of two matrices of size $n_a \times m_a$ and $n_b \times m_b$ is a matrix of size $(n_a n_b) \times (m_a m_b)$ given by

$$ C_{IJ} = a_{i_1 j_1} b_{i_2 j_2}, \quad \text{with} \quad I = i_2 + (i_1 - 1)n_b, \quad J = j_2 + (j_1 - 1)m_b. \quad (27) $$

In addition, the derivations presented in this section use the following properties of matrix multiplication, assuming compatible sizes,

$$ (a \otimes b)(f \otimes g) = af \otimes bg, \quad (28) $$

and inversion, assuming compatible sizes and that a and b are invertible,

$$ (a \otimes b)^{-1}(f \otimes g) = a^{-1}f \otimes b^{-1}g. \quad (29) $$
These properties also apply to tensor products in higher dimensions (see Deville et al. 2002, for more details).

Finally, let us introduce the one dimensional mass and convection elemental matrices

$$ \hat{\mathbf{m}}_{ij} = \int_{-1}^{1} n_i(\xi) n_j(\xi) d\xi, \quad \hat{\mathbf{c}}_{ij} = \int_{-1}^{1} n_i(\xi) \frac{\partial n_j}{\partial \xi}(\xi) d\xi. \qquad (30) $$

### 3.1 Case of constant Jacobian and velocity field

When the isoparametric mapping of Eq. (14) is affine, its Jacobian becomes constant for each element. If, additionally, the vector field $\mathbf{a}(\mathbf{x})$ is constant on the element, this can be used to accelerate the code substantially by adopting a quadrature-free implementation of the DG method (Atkins et al. 1998a; Sevilla et al. 2014). In this situation, the elemental matrices are computed as

$$ \mathbf{M} = |\mathbf{J}|\hat{\mathbf{M}}, \quad \mathbf{C} = \sum_{l=1}^{n_{sd}} G_l \hat{\mathbf{C}}^l, \quad \mathbf{M}^f = a_n^\bar{\|} \mathbf{J}^f \| \hat{\mathbf{M}}^f, \qquad (31) $$

where the functions $G_l$, given in Eq. (18), become constant in each element. The elemental matrices

$$ \hat{\mathbf{M}}_{IJ} = \int_{\hat{\Omega}} N_I N_J d\Omega, \quad \hat{\mathbf{C}}_{IJ}^{\ell} = \int_{\hat{\Omega}} N_I \frac{\partial N_J}{\partial \xi^{\ell}} d\Omega, \qquad (32) $$

for $I, J = 1, \dots, n_{en}$, are computed only once, in the reference element and stored. Similarly, the face matrix

$$ \hat{\mathbf{M}}_{IJ}^{f} = \begin{cases} \int_{\hat{\Gamma}} N_{I}^{f} N_{J}^{f} d\Gamma, & \text{if } \mathbf{x}_{I} \in \Gamma_{e}^{f} \\ 0 & \text{otherwise,} \end{cases} \qquad (33) $$

for $I = 1, \dots, n_{en}$ and $J = 1, \dots, n_{fn}$ is computed only once, in the reference face and stored.

This strategy leads to a very efficient implementation of the DG method where, at each time step, the evaluation of the residual of the system of ODEs, given by Eq. (22), reduces to a matrix scaling and matrix-vector product operations. It is also worth noting that in this situation, it is possible to compute the integrals of all the matrices in Eq. (32) exactly by using a numerical quadrature of order $2p$, for instance using a GLegQ.

#### 3.1.1 Residual evaluation for quadrilateral elements

Using the tensor product nature of the shape functions described in Section 2.4, the mass and convection matrices of the reference element given in Eq. (32) can be written as

$$ \hat{\mathbf{M}} = \hat{\mathbf{m}} \otimes \hat{\mathbf{m}}, \quad \hat{\mathbf{C}}^{1} = \hat{\mathbf{c}} \otimes \hat{\mathbf{m}}, \quad \hat{\mathbf{C}}^{2} = \hat{\mathbf{m}} \otimes \hat{\mathbf{c}}. \qquad (34) $$

Using these tensorised expressions of the elemental mass and convection matrices, the contribution of the convection term to the residual in Eq. (22) can be evaluated as

$$ \mathbf{C} = G_1 (\hat{\mathbf{m}}^{-1} \hat{\mathbf{c}}) \otimes \mathbf{i} + G_2 \mathbf{i} \otimes (\hat{\mathbf{m}}^{-1} \hat{\mathbf{c}}), \qquad (35) $$

where $\mathbf{i}$ is the identity matrix of size $(p+1) \times (p+1)$.

The mass matrices corresponding to the four edges of a quadrilateral can be written as

$$ \hat{\mathbf{M}}^{1} = \mathbf{e}_{1} \otimes \hat{\mathbf{m}}, \quad \hat{\mathbf{M}}^{2} = \mathbf{e}_{2} \otimes \hat{\mathbf{m}}, \quad \hat{\mathbf{M}}^{3} = \hat{\mathbf{m}} \otimes \mathbf{e}_{1}, \quad \hat{\mathbf{M}}^{4} = \hat{\mathbf{m}} \otimes \mathbf{e}_{2}, \qquad (36) $$

where the vectors $\mathbf{e}_1 = [1, 0, \dots, 0]^T$ and $\mathbf{e}_2 = [0, \dots, 0, 1]^T$ are of dimension $p+1$.

Figure 1 shows the reference quadrilateral element with the assumed local numbering of the nodes for the particular case of $p = 2$. The numbering of the edges assumed is such that the first edge contains vertices 1 and 3, the second edge contains vertices 7 and 9, the third face contains vertices 1 and 7 and the fourth face contains vertices 3 and 9.
<figure><img src="image_3.png" /><figcaption>Figure 1: Reference quadrilateral element, showing the local numbering of the nodes for $p = 2$.</figcaption></figure>

The contribution to the residual corresponding to the four edges of a quadrilateral can then be written in tensorised form as

$$ \begin{aligned} \mathcal{M}^1 &= a_n^{-} \frac{\|J^1\|}{|J|} (\hat{\mathbf{m}}^{-1} \mathbf{e}_1) \otimes \mathbf{i}, & \mathcal{M}^2 &= a_n^{-} \frac{\|J^2\|}{|J|} (\hat{\mathbf{m}}^{-1} \mathbf{e}_2) \otimes \mathbf{i}, \\ \mathcal{M}^3 &= a_n^{-} \frac{\|J^3\|}{|J|} \mathbf{i} \otimes (\hat{\mathbf{m}}^{-1} \mathbf{e}_1), & \mathcal{M}^4 &= a_n^{-} \frac{\|J^4\|}{|J|} \mathbf{i} \otimes (\hat{\mathbf{m}}^{-1} \mathbf{e}_2), \end{aligned} \quad (37) $$

### 3.1.2 Residual evaluation for hexahedral elements

The tensorial decompositions described above generalize to hexahedral elements in 3D. The mass and convection matrices of the reference element given in Eq. (32) can be written as

$$ \hat{\mathbf{M}} = \hat{\mathbf{m}} \otimes \hat{\mathbf{m}} \otimes \hat{\mathbf{m}}, \quad \hat{\mathbf{C}}^1 = \hat{\mathbf{c}} \otimes \hat{\mathbf{m}} \otimes \hat{\mathbf{m}}, \quad \hat{\mathbf{C}}^2 = \hat{\mathbf{m}} \otimes \hat{\mathbf{c}} \otimes \hat{\mathbf{m}}, \quad \hat{\mathbf{C}}^3 = \hat{\mathbf{m}} \otimes \hat{\mathbf{m}} \otimes \hat{\mathbf{c}}, \quad (38) $$

and the contribution of the convection term to the residual in Eq. (22) can be evaluated as

$$ \mathbf{C} = G_1 (\hat{\mathbf{m}}^{-1} \hat{\mathbf{c}}) \otimes \mathbf{i} \otimes \mathbf{i} + G_2 \mathbf{i} \otimes (\hat{\mathbf{m}}^{-1} \hat{\mathbf{c}}) \otimes \mathbf{i} + G_3 \mathbf{i} \otimes \mathbf{i} \otimes (\hat{\mathbf{m}}^{-1} \hat{\mathbf{c}}). \quad (39) $$

The mass matrices corresponding to the six faces of a hexahedral element can be written as

$$ \begin{aligned} \hat{\mathbf{M}}^1 &= \mathbf{e}_1 \otimes \hat{\mathbf{m}} \otimes \hat{\mathbf{m}}, & \hat{\mathbf{M}}^2 &= \mathbf{e}_2 \otimes \hat{\mathbf{m}} \otimes \hat{\mathbf{m}}, & \hat{\mathbf{M}}^3 &= \hat{\mathbf{m}} \otimes \mathbf{e}_1 \otimes \hat{\mathbf{m}} \\ \hat{\mathbf{M}}^4 &= \hat{\mathbf{m}} \otimes \mathbf{e}_2 \otimes \hat{\mathbf{m}}, & \hat{\mathbf{M}}^5 &= \hat{\mathbf{m}} \otimes \hat{\mathbf{m}} \otimes \mathbf{e}_1, & \hat{\mathbf{M}}^6 &= \hat{\mathbf{m}} \otimes \hat{\mathbf{m}} \otimes \mathbf{e}_2, \end{aligned} \quad (40) $$

leading to the following tensorised expressions for the face contributions to the residual:

$$ \begin{aligned} \mathcal{M}^1 &= a_n^{-} \frac{\|J^1\|}{|J|} (\hat{\mathbf{m}}^{-1} \mathbf{e}_1) \otimes \mathbf{i} \otimes \mathbf{i}, & \mathcal{M}^2 &= a_n^{-} \frac{\|J^2\|}{|J|} (\hat{\mathbf{m}}^{-1} \mathbf{e}_2) \otimes \mathbf{i} \otimes \mathbf{i}, \\ \mathcal{M}^3 &= a_n^{-} \frac{\|J^3\|}{|J|} \mathbf{i} \otimes (\hat{\mathbf{m}}^{-1} \mathbf{e}_1) \otimes \mathbf{i}, & \mathcal{M}^4 &= a_n^{-} \frac{\|J^4\|}{|J|} \mathbf{i} \otimes (\hat{\mathbf{m}}^{-1} \mathbf{e}_2) \otimes \mathbf{i}, \quad (41) \\ \mathcal{M}^5 &= a_n^{-} \frac{\|J^5\|}{|J|} \mathbf{i} \otimes \mathbf{i} \otimes (\hat{\mathbf{m}}^{-1} \mathbf{e}_1), & \mathcal{M}^6 &= a_n^{-} \frac{\|J^6\|}{|J|} \mathbf{i} \otimes \mathbf{i} \otimes (\hat{\mathbf{m}}^{-1} \mathbf{e}_2). \end{aligned} $$

Figure 2 shows the reference hexahedral element with the assumed local numbering of the nodes for the particular case of $p = 2$. The numbering of the faces assumed is such that the first face contains vertices 1, 3, 9 and 7, the second face contains vertices 19, 21, 27 and 25, the third face contains vertices 7, 9, 27 and 25, the fourth face contains vertices 1, 3, 21 and 19, the fifth face contains vertices 1, 7, 25 and 19 and the sixth face contains vertices 3, 9, 27 and 21.

### 3.1.3 Operations count to evaluate the residual

Using the tensorisation proposed above, the cost to compute the residual term is re-evaluated. We recall that the cost to compute the action on a vector of the tensorised matrix $\mathbf{B} = \mathbf{i} \otimes \mathbf{a}$, where $\mathbf{a}$ is a matrix of size $n \times n$ and $\mathbf{i}$ is the identity matrix of size $m \times m$, is $mn(2n - 1)$. We also recall that $\mathbf{i} \otimes \mathbf{i}$ is just the identity matrix of size $m^2 \times m^2$, and that the action of $\mathbf{i} \otimes \mathbf{a} \otimes \mathbf{i}$ has the same complexity as the action of $\mathbf{i} \otimes \mathbf{i} \otimes \mathbf{a}$ on that same vector.

It follows that, once the elemental matrices $G_\ell (\hat{\mathbf{m}}^{-1} \hat{\mathbf{c}})$, of size $(p+1) \times (p+1)$, are computed for an element, the evaluation of the term of the residual in Eq. (22) associated to such an
<figure><img src="image_4.png" /><figcaption>Figure 2: Reference hexahedral element, showing the local numbering of the nodes for $p = 2$.</figcaption></figure>

elemental matrix requires $(p+1)^{nsd} [2(p+1)-1]$ operations to perform the matrix-vector product. There are $n_{sd}$ such matrix-vector products to perform leading to $n_{op}^e = n_{sd}(p+1)^{nsd} [2(p+1)-1]$. The leading term in this evaluation is thus of the order of $O(2n_{sd}p^{nsd+1})$.

Considering vectors now, the general cost to compute the action on a vector of the tensorised matrix $B = i \otimes w$, where $w$ is a vector of size $n$ and $i$ is still the identity matrix of size $m \times m$, is $mn$. The cost to evaluate the term of the residual in Eq. (22) associated to an internal face now reduces to the application of $n_{sd}$ tensorised matrices per element (as explained in Section 2.5). The vector $(\hat{m}^{-1}e_i)$ of size $(p+1)$ is tensorised with an identity matrix of size $(p+1)^{nsd-1} \times (p+1)^{nsd-1}$, so the cost to apply a matrix $M^i$ is thus $n_{op}^f = (p+1)^{nsd}$.

These numbers of operations are summarized in Table 2.

<table><caption>Table 2: Number of operations required (per element) to compute the matrix-vector products in Eq. (22) using the GLegQ for numerical integration, assuming a constant Jacobian and a constant velocity field.</caption><thead><tr><th></th><th>Operations</th><th>Leading term</th></tr></thead><tbody><tr><td>$n_{op}^e$</td><td>$n_{sd}(p+1)^{nsd}[2(p+1)-1]$</td><td>$O(2n_{sd}p^{nsd+1})$</td></tr><tr><td>$n_{op}^f$</td><td>$n_{sd}(p+1)^{nsd}$</td><td>$O(n_{sd}p^{nsd})$</td></tr></tbody></table>

**Remark 4** The use of a tensorised version of the operators drastically reduces the memory requirements. It is thus necessary in this context to save only the the $n_{sd}$ elemental matrices $G_\ell (\hat{m}^{-1}\hat{c})$ of size $(p+1) \times (p+1)$ and the $n_{sd}$ vectors $(\hat{m}^{-1}e_i)$ of size $(p+1)$. The total memory requirements are thus of the order of $O(n_{e1}p^2)$ advocating in favor of an approach where these matrices are precomputed.

## 3.2 General case using Gauss-Lobatto quadratures

This section considers the general case, where the Jacobian and/or the velocity field are not constant. Despite the tensorised expressions obtained in the previous section are no longer valid in the general case, it is still possible to obtain efficient formulas for estimating the residual, but only when using the GLobQ for numerical integration. In order to simplify the comparisons in the next section, the mass, convection, and face mass matrices will be denoted by $M$, $C$ and $M^f$, respectively, when they have been integrated using GLobQ.

Using the delta Kronecker property of the shape functions, that is $N_I(\xi_J) = \delta_{IJ}$, along with the GLobQ numerical quadrature, simplified expressions for the elemental and face matrices in Eqs. (25) and (26) can be obtained. Indeed, the elemental mass matrix can be written as

$$M_{IJ} \approx \sum_{G=1}^{n_{ip}} N_I(\xi_G) N_J(\xi_G) \mathcal{J}(\xi_G) |\omega_G = \delta_{IJ} \mathcal{J}(\xi_I) |\omega_I, \quad (42)$$
with the trivial expression for the inverse

$$ \mathbb{M}_{IJ}^{-1} = \frac{\delta_{IJ}}{|\mathbf{J}(\xi_I)|\omega_I}. \qquad (43) $$

Similarly, the elemental convection matrix can be written as

$$ \mathbb{C}_{IJ} \approx \sum_{G=1}^{n_{1p}} N_I(\xi_G) \omega_G \sum_{\ell=1}^{n_{sd}} G_\ell(\xi_G) \frac{\partial N_J(\xi_G)}{\partial \xi_\ell} = \sum_{\ell=1}^{n_{sd}} G_\ell(\xi_I) \frac{\partial N_J(\xi_I)}{\partial \xi_\ell} \omega_I. \qquad (44) $$

Finally, the face mass matrix can be written as

$$ \mathbb{M}_{IJ}^{f} \approx \begin{cases} \sum_{G=1}^{n_{1p}} a_{n}^{-}(\xi_{G}) N_{I}(\xi_{G}) N_{J}(\xi_{G}) \|\mathbf{J}^{f}(\xi_{G})\|_{\omega_{G}} = \delta_{IJ} a_{n}^{-}(\xi_{I}) \|\mathbf{J}^{f}(\xi_{I})\|_{\omega_{I}} & \text{if } \mathbf{x}_{I} \in \Gamma_{e}^{f} \\ 0 & \text{otherwise.} \end{cases} \qquad (45) $$

### 3.2.1 Residual evaluation for quadrilateral elements

Observing in the 2D case that

$$ \frac{\partial N_J}{\partial \xi_1}(\xi_1^*) = n'_{j_1}(\xi_{i_1}) n_{j_2}(\xi_{i_2}), \quad \frac{\partial N_J}{\partial \xi_2}(\xi_1^*) = n_{j_1}(\xi_{i_1}) n'_{j_2}(\xi_{i_2}), \quad \omega_I = \omega_{i_1} \omega_{i_2}, \qquad (46) $$

the matrices $n_{\ell j} := \omega_\ell n_j(\xi_\ell)$ and $d_{\ell j} := \omega_\ell n'_j(\xi_\ell)$, of size $(p+1) \times (p+1)$ are introduced. These matrices are used to define the 2D derivative matrices, namely

$$ \mathbf{D}_1 := \mathbf{d} \otimes \mathbf{n}, \quad \mathbf{D}_2 := \mathbf{n} \otimes \mathbf{d}. \qquad (47) $$

Introducing the diagonal matrix $\hat{\mathbf{G}}_\ell$ (of size $n_{en} \times n_{en}$), whose elements are $G_\ell(\xi_I)$, the convection matrix can then be computed as

$$ \mathbb{C} = \hat{\mathbf{G}}_1(\mathbf{d} \otimes \mathbf{n}) + \hat{\mathbf{G}}_2(\mathbf{n} \otimes \mathbf{d}), \qquad (48) $$

and

$$ \mathbb{M}^{-1}\mathbb{C} = \mathbb{M}^{-1}\hat{\mathbf{G}}_1(\mathbf{d} \otimes \mathbf{n}) + \mathbb{M}^{-1}\hat{\mathbf{G}}_2(\mathbf{n} \otimes \mathbf{d}), \qquad (49) $$

where $\mathbb{M}^{-1}$ is the diagonal matrix defined in Eq. (43). It is worth noting that each term in the sum above is the product of a diagonal matrix ($\mathbb{M}^{-1}\hat{\mathbf{G}}_i$) by a matrix that has a tensorised form.

For the evaluation of the part of the residual related to edges, the following diagonal formula is obtained

$$ \mathbb{M}^{-1} \mathbb{M}_{IJ}^{f} \approx \begin{cases} \delta_{IJ} a_{n}^{-}(\xi_{I}) \frac{\|\mathbf{J}^{f}(\xi_{I})\|_{\omega_{I}}}{|\mathbf{J}(\xi_{I})|\omega_{I}} & \text{if } \mathbf{x}_{I} \in \Gamma_{e}^{f} \\ 0 & \text{otherwise.} \end{cases} \qquad (50) $$

### 3.2.2 Residual evaluation for hexahedral elements

The decompositions described above generalize to 3D tensor-product elements in the following way for the convection matrix,

$$ \mathbb{C} = \hat{\mathbf{G}}_1(\mathbf{d} \otimes \mathbf{n} \otimes \mathbf{n}) + \hat{\mathbf{G}}_2(\mathbf{n} \otimes \mathbf{d} \otimes \mathbf{n}) + \hat{\mathbf{G}}_3(\mathbf{n} \otimes \mathbf{n} \otimes \mathbf{d}), \qquad (51) $$

and the contribution of the convection matrices to the residual of the ODE system,

$$ \mathbb{M}^{-1}\mathbb{C} = \mathbb{M}^{-1}\hat{\mathbf{G}}_1(\mathbf{d} \otimes \mathbf{n} \otimes \mathbf{n}) + \mathbb{M}^{-1}\hat{\mathbf{G}}_2(\mathbf{n} \otimes \mathbf{d} \otimes \mathbf{n}) + \mathbb{M}^{-1}\hat{\mathbf{G}}_3(\mathbf{n} \otimes \mathbf{n} \otimes \mathbf{d}). \qquad (52) $$

Finally, the contribution of the face mass matrices to the ODE residual is the same as in Eq. (50).
### 3.2.3 Operations count to evaluate the residual

To evaluate the cost of computing the action of the tensorised matrix $a \otimes b$ on a vector, where $a$ and $b$ are matrices of size $n \times n$, the decomposition $a \otimes b = (a \otimes i)(i \otimes b)$ is introduced, where $i$ is the identity matrix of size $n \times n$. As previously shown, the cost of evaluating $(i \otimes b)u$, which is also the cost of $(a \otimes i)u$, is equal to $n^2(2n - 1)$. The total cost to compute the action of the tensorised matrix $a \otimes b$ on a vector is thus $2n^2(2n - 1)$. Similarly, in 3D the decomposition $a \otimes b \otimes c = (a \otimes i \otimes i)(i \otimes b \otimes i)(i \otimes i \otimes c)$ is used to obtain that the cost to compute the action of the tensorised matrix $a \otimes b \otimes c$ on a vector is $3n^3(2n - 1)$.

Using Eq. (49), the action of the matrix $M^{-1}C$ on a vector can be therefore decomposed into: (i) the action of $n_{sd}$ tensorised matrices on that vector, with cost $n_{sd}(p + 1)^{n_{sd}}(2p + 1)$ per matrix; and (ii) the action of $n_{sd}$ diagonal matrices of size $(p + 1)^{n_{sd}}$, with a cost of $n_{sd}(p + 1)^{n_{sd}}$. The total cost to evaluate the residual part associated with the convection is therefore of order $n_{sd}^2(p + 1)^{n_{sd}}(2p + 1) + n_{sd}(p + 1)^{n_{sd}}$, with a leading order similar to the constant case.

The cost to evaluate the term of the residual in Eq. (22) associated to an internal face now reduces to the application of $n_{sd}$ diagonal matrices of size $(p + 1)^{n_{sd}-1} \times (p + 1)^{n_{sd}-1}$, with a cost of $n_{sd}(p + 1)^{n_{sd}-1}$.

These numbers of operations are summarized in Table 3.

<table><caption>**Table 3:** Number of operations required (per element) to compute the matrix-vector products in Eq. (22) using the GLobQ for numerical integration.</caption><thead><tr><th></th><th>Operations</th><th>Leading term</th></tr></thead><tbody><tr><td>$n_{op}^e$</td><td>$n_{sd}^2(p + 1)^{n_{sd}}(2p + 1) + n_{sd}(p + 1)^{n_{sd}}$</td><td>$O(n_{sd}^2p^{n_{sd}+1})$</td></tr><tr><td>$n_{op}^f$</td><td>$n_{sd}(p + 1)^{n_{sd}-1}$</td><td>$O(n_{sd}p^{n_{sd}-1})$</td></tr></tbody></table>

**Remark 5** Similarly to the case with a constant Jacobian and velocity the use of a tensorized version of the operators drastically reduces the memory requirements. Due to the use of Gauss-Lobatto quadrature rules, most of the matrices are diagonal ones. It is thus necessary in this context to save only the matrices $n$ and $d$ of size $(p + 1) \times (p + 1)$, the $n_{sd}$ diagonal matrices $\hat{G}_l$ of size $(p + 1)^{n_{sd}} \times (p + 1)^{n_{sd}}$, the diagonal matrices $M^{-1}$ together with the diagonal matrices $M^f$ of size $(p + 1)^{n_{sd}-1} \times (p + 1)^{n_{sd}-1}$. The total memory requirements are thus of the order of $O(p^{n_{sd}})$ advocating in favor of an approach where these matrices are precomputed.

## 4 A DG method with mixed quadratures

In the previous section, it was concluded that, in general, the most expensive part of the computation of the residual of Eq. (22) corresponds to the convection contribution. In this section, a formal proof is provided showing that under certain hypothesis (explained further down), the contribution to the residual given by the convection term is identical if GLegQ or GLobQ are used. This result is in principle not expected because the GLobQ does not provide an exact integration for the entries of the mass matrix. Furthermore, numerical integration with GLobQ does not provide exact integration for the convection matrix under the hypothesis considered. So even if both the mass and convection matrices are not integrated exactly, the result provided here shows that the contribution to the residual can be computed exactly.

This is the main motivation for proposing a new DG method using mixed quadratures (MixQ). The proposed approach uses a GLobQ for computing the contribution to the residual due to convection whereas a GLegQ is used to compute the face contribution. Numerical examples in Section 5 will be used to numerically demonstrate that, in the most general case, it is indeed the inaccurate integration of the face mass contribution (with the GLobQ) that influences accuracy most drastically, rather than the inaccurate integration of the convection contribution. Therefore, the proposed approach is considered attractive as it will improve the performance of the code by reducing the cost of the most expensive part of the computation, i.e., the convection part of the residual, whilst maintaining the overall accuracy of the standard DG method with GLegQ.
## 4.1 A surprising equality for residual evaluation

This section formally proves that, somewhat surprisingly, the evaluation of the convection part of the residual is exact (given an hypothesis that will be explained below) with GLobQ even when both the convection and mass matrices contain errors due to the numerical integration.

Let us first consider the mass and convection matrices obtained by integration with the GLobQ. Using Eqs. (43), (25) and (18), the expression for the elemental mass and convection matrices can be written as

$$ \mathbb{M}_{IJ}^{-1} = \frac{\delta_{IJ}}{|\mathbb{J}(\xi_I)|\omega_I}, \quad \mathbb{C}_{IJ} = \omega_I |\mathbb{J}(\xi_I)| \sum_{k, \ell=1}^{n_{sd}} a_k(\varphi(\xi_I)) J_{kl}^{-1}(\xi_I) \frac{\partial N_J(\xi_I)}{\partial \xi_\ell}, \quad (53) $$

so that the convection part of the residual is given by

$$ (\mathbb{M}^{-1}\mathbb{C})_{IJ} = \sum_{k, \ell=1}^{n_{sd}} a_k(\varphi(\xi_I)) J_{kl}^{-1}(\xi_I) \frac{\partial N_J(\xi_I)}{\partial \xi_\ell}. \quad (54) $$

We first note that any polynomial $P(\xi)$ over $\hat{\Omega}$, of order $p$ in each direction, can be exactly represented on the basis of functions $N_I(\xi)$ in the following manner:

$$ P(\xi) = \sum_{I=0}^{n_{en}} N_I(\xi) \int_{\hat{\Omega}} N_I(\xi) P(\xi) d\Omega. \quad (55) $$

Secondly, if the order of that polynomial is less than or equal to $p-1$ in each direction, numerical integration with the GLobQ is exact (because $N_I(\xi)P(\xi)$ is then a polynomial of order less than $2p-1$), so that we further have

$$ P(\xi) = \sum_{I=0}^{n_{en}} N_I(\xi) P(\xi_I). \quad (56) $$

Therefore, if, for any $1 \le I, J \le n_{en}$, $\sum_{k, \ell=1}^{n_{sd}} a_k(\varphi(\xi)) J_{kl}^{-1}(\xi) \frac{\partial N_J(\xi)}{\partial \xi_\ell}$ is a polynomial of order less than $p-1$ in each direction,

$$ \sum_{k, \ell=1}^{n_{sd}} a_k(\varphi(\xi)) J_{kl}^{-1}(\xi) \frac{\partial N_J(\xi)}{\partial \xi_\ell} = \sum_{I=0}^{n_{en}} N_I(\xi) (\mathbb{M}^{-1}\mathbb{C})_{IJ}. \quad (57) $$

If we now consider the convection and mass matrices integrated exactly (for instance using GLegQ with enough integration points), as in Eq. (17) and (18), we obtain, using also Eq. (57),

$$ \sum_{J=1}^{n_{en}} \mathbf{M}_{IJ} (\mathbb{M}^{-1}\mathbb{C})_{JK} = \int_{\hat{\Omega}} N_I(\xi) |\mathbb{J}(\xi)| \left( \sum_{J=1}^{n_{en}} N_J(\xi) (\mathbb{M}^{-1}\mathbb{C})_{JK} \right) d\Omega = C_{IK}, \quad (58) $$

or, equivalently,

$$ \mathbb{M}^{-1}\mathbb{C} = \mathbb{M}^{-1}\mathbb{C}. \quad (59) $$

Therefore, provided that for any $1 \le I, J \le n_{en}$, $\sum_{k, \ell=1}^{n_{sd}} a_k(\varphi(\xi)) J_{kl}^{-1}(\xi) \frac{\partial N_J(\xi)}{\partial \xi_\ell}$ is a polynomial of order less than $p-1$ in each direction, the product $\mathbb{M}^{-1}\mathbb{C}$ is evaluated exactly, even though both the matrices $\mathbb{M}$ and $\mathbb{C}$ may be wrongly integrated. Since this product is the only quantity required during the evaluation of the convective part of the residual, this means that using GLobQ essentially brings no additional error with respect to GLegQ. Of course, this says nothing of the face mass part or the evaluation of the residual.

Note that the constraint on the polynomial order is quite strong, but is at least verified for a constant velocity field and a constant Jacobian.
## 4.2 The DG method with mixed quadratures

Combined with the numerical tests of next section (in both cases when the polynomial hypothesis is verified or not), the previous section essentially states that integration error bears little weight in the evaluation of the convective part of the residual. Otherwise said, when aiming for a certain overall accuracy, it is less important to integrate exactly the convective part of the residual than the face mass part. This leads us to propose a novel DG method, with the following characteristics:

* Interpolation points are the nodes of the GLobQ;

* Residual is evaluated using Eq. (22) where

- The convective part is $CU = M^{-1}CU$, where the elements of $M^{-1}C$ are integrated using GLobQ, as in Eq. (54);

- The face mass part is $\sum_{f=1}^{n_{fa}^e} M^f[\llap{\text{U}}f] = \sum_{f=1}^{n_{fa}^e} M^{-1}M^f[\llap{\text{U}}f]$ where the elements of $M^{-1}M^f$ are integrated exactly (for instance using GLegQ in Eq. (25) and (26) with enough integration points).

The overall operations count for the evaluation of the residual is therefore $[2(p + 1)^{n_{sd} - 1} + n_{sd}^2(2p + 1) + n_{sd} - 1](p + 1)^{n_{sd}}$, with leading order $O(2p^{2n_{sd} - 1})$ or $O(2n_{sd}^2p^{n_{sd} + 1})$ depending on the dimension $n_{sd}$ and $p$. It is interesting to note here that except for $n_{sd} = 3$ and high order $p$, the contributions of the two terms in the residual are on the same order. In the worst case, the proposed method is therefore slightly more expensive than the full DGSEM, but next section will show that it is also much more accurate. With respect to the DG method with exact integration, the proposed method is more efficient, and the examples of next section will show that it is as accurate.

Table 4 summarises the number of operations of the DG methods with GLegQ, GLobQ and with the proposed MixQ.

<table><caption>**Table 4:** Orders of magnitude of the numbers of operations required (per element) to compute the matrix-vector products to evaluate the ODE residual in Eq. (16) for the different DG approaches considered. The results, from Tables 1, 2 and 3, are particularised for two and three dimensions and displayed separately for the particular case of constant velocity and constant Jacobian alongside the general case.</caption><thead><tr><th colspan="4">Constant Jacobian and velocity</th></tr><tr><th></th><th colspan="2">n<sub>sd</sub> = 2</th><th colspan="2">n<sub>sd</sub> = 3</th></tr><tr><th></th><th>n<sub>op</sub><sup>e</sup></th><th>n<sub>op</sub><sup>f</sup></th><th>n<sub>op</sub><sup>e</sup></th><th>n<sub>op</sub><sup>f</sup></th></tr></thead><tbody><tr><th>GLegQ</th><td>O(4p<sup>3</sup>)</td><td>O(2p<sup>2</sup>)</td><td>O(6p<sup>4</sup>)</td><td>O(3p<sup>3</sup>)</td></tr><tr><th>GLobQ</th><td>O(4p<sup>3</sup>)</td><td>O(2p)</td><td>O(9p<sup>4</sup>)</td><td>O(3p<sup>2</sup>)</td></tr><tr><th>MixQ</th><td>O(4p<sup>3</sup>)</td><td>O(2p<sup>2</sup>)</td><td>O(9p<sup>4</sup>)</td><td>O(3p<sup>3</sup>)</td></tr><tr><th colspan="4">Non-constant Jacobian and/or velocity</th></tr><tr><th></th><th colspan="2">n<sub>sd</sub> = 2</th><th colspan="2">n<sub>sd</sub> = 3</th></tr><tr><th></th><th>n<sub>op</sub><sup>e</sup></th><th>n<sub>op</sub><sup>f</sup></th><th>n<sub>op</sub><sup>e</sup></th><th>n<sub>op</sub><sup>f</sup></th></tr><tr><th>GLegQ</th><td>O(2p<sup>4</sup>)</td><td>O(2p<sup>3</sup>)</td><td>O(2p<sup>6</sup>)</td><td>O(2p<sup>5</sup>)</td></tr><tr><th>GLobQ</th><td>O(4p<sup>3</sup>)</td><td>O(2p)</td><td>O(9p<sup>4</sup>)</td><td>O(3p<sup>2</sup>)</td></tr><tr><th>MixQ</th><td>O(4p<sup>3</sup>)</td><td>O(2p<sup>3</sup>)</td><td>O(9p<sup>4</sup>)</td><td>O(2p<sup>5</sup>)</td></tr></tbody></table>

# 5 Numerical examples

This section presents a set of numerical tests to evaluate the accuracy of the DG approaches considered in this work. The examples involve problems with constant and non-constant velocity fields both in two and three dimensions. In addition, the use of Cartesian and non-Cartesian grids is also studied, to evaluate the influence of having elements with non-constant Jacobian of the isoparametric mapping. The discussion also considers the computational cost of the DG approaches considered, as well as the impact of stability.
<figure><img src="image_5.png" /><figcaption>Figure 3: Cartesian meshes of quadrilateral elements.</figcaption></figure>

<table><caption>Table 5: Relative error in the $L^2(\Omega)$ norm ($\epsilon$), and rate of convergence ($r$) for the DG methods with GLegQ and GLobQ for the 2D problem with constant velocity using Cartesian meshes.</caption><thead><tr><th>p</th><th>h</th><th>ε<sub>GLegQ</sub> = ε<sub>MixQ</sub></th><th>r<sub>GLegQ</sub> = r<sub>MixQ</sub></th><th>ε<sub>GLobQ</sub></th><th>r<sub>GLobQ</sub></th><th>ε<sub>GLobQ</sub>/ε<sub>GLegQ</sub></th></tr></thead><tbody><tr><td>1</td><td>0.5000</td><td>4.5 × 10<sup>-1</sup></td><td>–</td><td>9.7 × 10<sup>-1</sup></td><td>–</td><td>2.1</td></tr><tr><td>1</td><td>0.2500</td><td>1.4 × 10<sup>-1</sup></td><td>1.7</td><td>6.0 × 10<sup>-1</sup></td><td>0.7</td><td>4.2</td></tr><tr><td>1</td><td>0.1250</td><td>3.3 × 10<sup>-2</sup></td><td>2.1</td><td>2.3 × 10<sup>-1</sup></td><td>1.4</td><td>6.9</td></tr><tr><td>1</td><td>0.0625</td><td>7.7 × 10<sup>-3</sup></td><td>2.1</td><td>6.6 × 10<sup>-2</sup></td><td>1.8</td><td>8.6</td></tr><tr><td>2</td><td>0.5000</td><td>8.8 × 10<sup>-2</sup></td><td>–</td><td>3.1 × 10<sup>-1</sup></td><td>–</td><td>3.5</td></tr><tr><td>2</td><td>0.2500</td><td>1.0 × 10<sup>-2</sup></td><td>3.1</td><td>4.9 × 10<sup>-2</sup></td><td>2.7</td><td>5.0</td></tr><tr><td>2</td><td>0.1250</td><td>1.2 × 10<sup>-3</sup></td><td>3.0</td><td>6.2 × 10<sup>-3</sup></td><td>3.0</td><td>5.1</td></tr><tr><td>2</td><td>0.0625</td><td>1.5 × 10<sup>-4</sup></td><td>3.0</td><td>7.6 × 10<sup>-4</sup></td><td>3.0</td><td>5.1</td></tr><tr><td>3</td><td>0.5000</td><td>1.3 × 10<sup>-2</sup></td><td>–</td><td>4.9 × 10<sup>-2</sup></td><td>–</td><td>4.0</td></tr><tr><td>3</td><td>0.2500</td><td>7.9 × 10<sup>-4</sup></td><td>4.0</td><td>3.4 × 10<sup>-3</sup></td><td>3.9</td><td>4.3</td></tr><tr><td>3</td><td>0.1250</td><td>4.8 × 10<sup>-5</sup></td><td>4.0</td><td>2.2 × 10<sup>-4</sup></td><td>4.0</td><td>4.5</td></tr><tr><td>3</td><td>0.0625</td><td>3.0 × 10<sup>-6</sup></td><td>4.0</td><td>1.4 × 10<sup>-5</sup></td><td>4.0</td><td>4.5</td></tr><tr><td>4</td><td>0.5000</td><td>1.7 × 10<sup>-3</sup></td><td>–</td><td>6.5 × 10<sup>-3</sup></td><td>–</td><td>3.9</td></tr><tr><td>4</td><td>0.2500</td><td>5.1 × 10<sup>-5</sup></td><td>5.0</td><td>2.1 × 10<sup>-4</sup></td><td>4.9</td><td>4.1</td></tr><tr><td>4</td><td>0.1250</td><td>1.6 × 10<sup>-6</sup></td><td>5.0</td><td>6.8 × 10<sup>-6</sup></td><td>5.0</td><td>4.2</td></tr><tr><td>4</td><td>0.0625</td><td>5.0 × 10<sup>-8</sup></td><td>5.0</td><td>2.1 × 10<sup>-7</sup></td><td>5.0</td><td>4.2</td></tr></tbody></table>

## 5.1 Two dimensional examples

The first two dimensional example considers the propagation of a sinusoidal wave with a constant velocity field in a Cartesian mesh. The domain is $\Omega = [0, 1]^2$ and the velocity field is given by $\mathbf{a} = [\cos(\pi/6), \sin(\pi/6)]$. The exact solution of the problem is given by

$$u(\mathbf{x}, t) = \sin(2\pi(\mathbf{a} \cdot \mathbf{x} - t)). \qquad (60)$$

The final time for the simulations is taken as $T = 1$ and the time step is selected small enough so that the error induced by the time marching scheme is negligible when compared to the error induced by the spatial discretisation.

A set of Cartesian meshes is first utilised to test the optimal convergence properties of the DG methods with Gauss-Legendre (GLegQ) and Gauss-Lobatto (GLobQ) quadratures. As detailed in Section 4.1, the proposed approach with mixed quadratures (MixQ) leads to identical results when compared to the DG method with GLegQ.

The mesh corresponding to the $i$-th level of refinement is denoted by $M^i$ and contains $2^i \times 2^i$ quadrilateral elements with characteristic element size $h = 2^{-i}$, as shown in figure 3.

Table 5 shows the relative error, measured in the $L^2(\Omega)$ norm and denoted by $\epsilon$, for the DG methods with GLegQ and GLobQ for a polynomial approximation of order ranging from $p = 1$ up to $p = 4$. The table also shows the rates of convergence, denoted by $r$ and the ratio between the error of both approaches.

The results show that both approaches provide an optimal rate of convergence, approximately $p + 1$, in the $L^2(\Omega)$ norm. The results also reveal the extra accuracy provided by the DG method
<figure><img src="image_6.png" /><figcaption>**Figure 4:** Influence, in the case of constant velocity and constant jacobian, of time step on the $L^2(\Omega)$-error at $T = 1$ for different polynomial orders, meshes and method.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5:** Non-Cartesian meshes of quadrilateral elements.</figcaption></figure>

with GLegQ, when compared to the results of the DG method with GLobQ on the same mesh. It is worth emphasising that the computation of the convection term in both cases is identical because the velocity is constant and the Jacobian of the isoparametric mapping is constant for the selected meshes. Therefore, the extra accuracy is due to the more accurate integration provided by the GLegQ for the contribution to the residual of the face term. For the lowest order case, the use of GLegQ provides an error nine times lower than that produced by the DG method with GLobQ. For higher orders, the use of GLegQ consistently provides an accuracy between four to five times higher than the DG method with GLobQ.

In terms of the computational cost, the DG method with GLegQ requires $O(4p^3)$ and $O(2p^2)$ operations for computing the matrix-vector products corresponding to the convection and face contributions to the residual, respectively. The DG method with GLobQ requires $O(4p^3)$ and $O(2p)$ operations for computing the matrix-vector products corresponding to the convection and face contributions to the residual, respectively. Therefore, it is apparent that the approach with GLegQ is beneficial as the dominant term, corresponding to the convection term, has the same order of operations.

For completeness, we propose an additional series of numerical tests to try and evaluate the influence of stability on this conclusion. In Fig. 4 is represented the influence of the time step on the error (at $T = 1$) for all three methods, different polynomial orders, and the refined meshes ($M^3$ and $M^4$). The instability appears as an explosion of the error above a given critical time step. As already discussed, that time step is the same for GLegQ and MixQ. We observe that the stability time is consistently twice as large for GLobQ than GLegQ/MixQ. This means that the conclusion of the previous paragraph is somehow lessened as larger time steps can be considered (and hence twice less computational effort) with GLobQ than the other two method. However, MixQ remain the overall best method, with 2 to 4 times less error for the same computational cost as GLobQ, and the same error for a smaller computational cost as GLegQ.

To evaluate the influence of having elements with non-constant Jacobian, the same problem is now solved using the non-Cartesian meshes shown in Figure 5. The meshes have been obtained by randomly perturbing the position of the interior nodes of the Cartesian meshes of Figure 3.

Table 6 shows the relative error, measured in the $L^2(\Omega)$ norm, for the DG methods with GLegQ, GLobQ and MixQ. The table also shows the rates of convergence and the ratio between
the error of the DG methods with GLobQ and MixQ with respect to the GLegQ.

<table><caption>**Table 6:** Relative error in the $L^2(\Omega)$ norm ($\epsilon$), and rate of convergence ($r$) for the DG methods with GLegQ, GLobQ and MixQ for the 2D problem with constant velocity using non-Cartesian meshes.</caption><thead><tr><th>p</th><th>h</th><th>ε<sub>GLegQ</sub></th><th>r<sub>GLegQ</sub></th><th>ε<sub>GLobQ</sub></th><th>r<sub>GLobQ</sub></th><th>ε<sub>MixQ</sub></th><th>r<sub>MixQ</sub></th><th>ε<sub>GLobQ</sub> / ε<sub>GLegQ</sub></th><th>ε<sub>MixQ</sub> / ε<sub>GLegQ</sub></th></tr></thead><tbody><tr><td>1</td><td>0.5737</td><td>4.7 × 10<sup>-1</sup></td><td>—</td><td>1.0</td><td>—</td><td>4.8 × 10<sup>-1</sup></td><td>—</td><td>2.1</td><td>1.0</td></tr><tr><td>1</td><td>0.2991</td><td>1.5 × 10<sup>-1</sup></td><td>1.8</td><td>6.7 × 10<sup>-1</sup></td><td>0.6</td><td>1.5 × 10<sup>-1</sup></td><td>1.8</td><td>4.4</td><td>1.0</td></tr><tr><td>1</td><td>0.1699</td><td>3.5 × 10<sup>-2</sup></td><td>2.6</td><td>2.8 × 10<sup>-1</sup></td><td>1.6</td><td>3.8 × 10<sup>-2</sup></td><td>2.4</td><td>7.9</td><td>1.1</td></tr><tr><td>1</td><td>0.0871</td><td>8.4 × 10<sup>-3</sup></td><td>2.1</td><td>7.8 × 10<sup>-2</sup></td><td>1.9</td><td>1.2 × 10<sup>-2</sup></td><td>1.8</td><td>9.3</td><td>1.4</td></tr><tr><td>2</td><td>0.5737</td><td>9.6 × 10<sup>-2</sup></td><td>—</td><td>2.6 × 10<sup>-1</sup></td><td>—</td><td>9.8 × 10<sup>-2</sup></td><td>—</td><td>2.7</td><td>1.0</td></tr><tr><td>2</td><td>0.2991</td><td>1.2 × 10<sup>-2</sup></td><td>3.1</td><td>3.9 × 10<sup>-2</sup></td><td>2.9</td><td>1.3 × 10<sup>-2</sup></td><td>3.1</td><td>3.1</td><td>1.0</td></tr><tr><td>2</td><td>0.1699</td><td>1.5 × 10<sup>-3</sup></td><td>3.7</td><td>5.3 × 10<sup>-3</sup></td><td>3.5</td><td>1.6 × 10<sup>-3</sup></td><td>3.7</td><td>3.5</td><td>1.0</td></tr><tr><td>2</td><td>0.0871</td><td>1.9 × 10<sup>-4</sup></td><td>3.1</td><td>6.2 × 10<sup>-4</sup></td><td>3.2</td><td>2.4 × 10<sup>-4</sup></td><td>2.8</td><td>3.3</td><td>1.3</td></tr><tr><td>3</td><td>0.5737</td><td>1.4 × 10<sup>-2</sup></td><td>—</td><td>4.2 × 10<sup>-2</sup></td><td>—</td><td>1.4 × 10<sup>-2</sup></td><td>—</td><td>3.0</td><td>1.0</td></tr><tr><td>3</td><td>0.2991</td><td>8.5 × 10<sup>-4</sup></td><td>4.3</td><td>2.4 × 10<sup>-3</sup></td><td>4.4</td><td>8.5 × 10<sup>-4</sup></td><td>4.3</td><td>2.9</td><td>1.0</td></tr><tr><td>3</td><td>0.1699</td><td>5.9 × 10<sup>-5</sup></td><td>4.7</td><td>1.7 × 10<sup>-4</sup></td><td>4.7</td><td>6.0 × 10<sup>-5</sup></td><td>4.7</td><td>2.9</td><td>1.0</td></tr><tr><td>3</td><td>0.0871</td><td>4.0 × 10<sup>-6</sup></td><td>4.0</td><td>1.2 × 10<sup>-5</sup></td><td>4.0</td><td>4.2 × 10<sup>-6</sup></td><td>4.0</td><td>2.9</td><td>1.1</td></tr><tr><td>4</td><td>0.5737</td><td>2.1 × 10<sup>-3</sup></td><td>—</td><td>5.3 × 10<sup>-3</sup></td><td>—</td><td>2.1 × 10<sup>-3</sup></td><td>—</td><td>2.6</td><td>1.0</td></tr><tr><td>4</td><td>0.2991</td><td>8.2 × 10<sup>-5</sup></td><td>4.9</td><td>2.1 × 10<sup>-4</sup></td><td>4.9</td><td>8.3 × 10<sup>-5</sup></td><td>5.0</td><td>2.6</td><td>1.0</td></tr><tr><td>4</td><td>0.1699</td><td>2.6 × 10<sup>-6</sup></td><td>6.1</td><td>7.2 × 10<sup>-6</sup></td><td>6.0</td><td>2.6 × 10<sup>-6</sup></td><td>6.1</td><td>2.7</td><td>1.0</td></tr><tr><td>4</td><td>0.0871</td><td>7.9 × 10<sup>-8</sup></td><td>5.2</td><td>2.2 × 10<sup>-7</sup></td><td>5.2</td><td>8.0 × 10<sup>-8</sup></td><td>5.2</td><td>2.7</td><td>1.0</td></tr></tbody></table>

The results show the optimal rate of convergence for the three DG approaches. In terms of the accuracy, the DG method with GLegQ and the proposed DG method with MixQ provide almost identical results. This shows that the extra error present in the DG method with MixQ that is induced by the numerical integration of the part of the residual corresponding to the convection is almost negligible. In contrast when GLobQ are used for both the convection and the face integrals, the error increases by a factor between 2.7 and 3.3 for $p > 1$.

When comparing the computational costs of the three approaches, the most expensive is the DG method with GLegQ as the number of operations associated to the matrix-vector products of the convection and face contributions to the residual are $O(2p^4)$ and $O(2p^3)$, respectively. For the DG method with GLobQ the number of operations associated to the matrix-vector products of the convection and face contributions to the residual are $O(4p^3)$ and $O(2p)$, respectively. With the proposed DG method with MixQ $O(4p^3)$ and $O(2p^3)$ operations are required to compute the matrix-vector products of the convection and face contributions to the residual, respectively. In all cases the cost is dominated by the convection term, which requires $O(2p^4)$, $O(4p^3)$, and $O(4p^3)$ operations for the DG method with GLegQ, GLobQ and MixQ, respectively. In fact it is worth noting that the cost of the approach with MixQ is only slightly higher than the cost of the DG method with GLobQ, but as a result the accuracy is approximately three times higher. It is also worth noting that the proposed DG method with MixQ is more efficient than the DG method with GLegQ. In terms of accuracy both methods provide almost identical results but in terms of the cost, the proposed method has a dominant cost of $O(p^3)$ whereas the DG method with GLegQ has a dominant cost of $O(p^4)$.

As in the previous case, the conclusions are lessened when considering stability (see Fig. 6), but the MixQ remains 1.5 to 3 times more efficient than the GLobQ, and as precise as the GLeqQ for the same computational cost.

The next example considers a deformational flow, similar to the tests employed to assess the accuracy of advection schemes for climate and weather applications (Kent 2019). The velocity field is given by

$$a(\mathbf{x}, t) = \cos(\pi t/T) \begin{bmatrix} 0.1 \sin^2(\pi x_1) \sin(2\pi x_2) \\ 0.1 \sin(2\pi x_1) \sin^2(\pi x_2) \end{bmatrix}^T \quad (61)$$
<figure><img src="image_8.png" /><figcaption>**Figure 6:** Influence, in the case of constant velocity and non-constant jacobian, of time step on the $L^2(\Omega)$-error at $T = 1$ for different polynomial orders, meshes and method.</figcaption></figure>

and the final time is $T = 1$.

The initial condition is given by $u_0(x) = 0.5 + 0.5 \sin(2\pi x_1) \sin(2\pi x_2)$ and periodic boundary conditions are considered on the boundary of the domain $\partial\Omega$. It is worth noting that the problem has no known analytical solution but the cosine term in Eq. (61) ensures that the flow reverses at time $t = T/2$ and reaches the same state as the initial condition at time $t = T$. This idea was originally proposed in (Leveque 1996) to test the accuracy of numerical schemes in problems with no known analytical solution.

Table 7 shows the relative error, measured in the $L^2(\Omega)$ norm, for the three DG approaches considered in the previous section using Cartesian meshes.

The results show that all methods achieve the optimal rate of convergence, but it is worth noting that the DG method with GLobQ takes substantially longer to reach the asymptotic range of convergence compared to the DG methods with GLegQ and MixQ. Once the asymptotic rate of convergence is reached, the DG methods with GLegQ and MixQ provide between 2 and 2.5 times more accurate results than the DG method with GLobQ for $p > 1$. The results also show that in some cases the proposed DG method with MixQ is marginally more accurate than the DG method with GLegQ. This is attributed to the non-exact integration of the convection term, which involves a non-polynomial velocity field.

As the cost associated to the matrix-vector products when the velocity field or the Jacobian are not constant, the conclusions are very similar to the ones obtained in the previous example. The dominant cost, associated to the convection term, is $O(2p^4)$, $O(4p^3)$, and $O(4p^3)$ operations for the DG method with GLegQ, GLobQ and MixQ, respectively. Therefore, the proposed DG method with MixQ has the same dominant cost as the GLobQ with an accuracy more than two times higher. The proposed DG method with MixQ is also more efficient than the DG method with GLegQ because the accuracy is almost identical but the dominant cost is lower.

The last two dimensional example involves the solution of the deformational flow with non-constant velocity by using non-Cartesian meshes. Table 8 shows the relative error, measured in the $L^2(\Omega)$ norm, for the three DG approaches considered.

The results are both, qualitatively and quantitatively, very similar to the ones obtained in the previous example with a non-constant velocity field using Cartesian meshes. In all cases the proposed DG method with MixQ produces almost identical results when compared to the DG method with GLegQ, despite the extra accuracy of the numerical quadratures used for the convection term in the DG method with GLegQ. This denotes that even if the velocity and Jacobian are both not constant, the error induced by a numerical integration of the convection term by a Gauss-Lobatto quadrature compared to the error induced by a Gauss-Legendre quadrature is below the discretisation error.

This time however, the conclusions are not so clear when considering stability (see Fig. 7), as the higher precision of the MixQ with respect to GLobQ is somewhat compensated by the higher stability of the GlobQ, so that both methods seem equivalent in terms of efficiency in the case of non-constant velocity.

As a conclusion for this series of 2D numerical tests, the MixQ method is an attractive alternative to the DG methods with GLegQ and GLobQ. At least for constant velocity, the
<table><caption>**Table 7:** Relative error in the $L^2(\Omega)$ norm ($\epsilon$), and rate of convergence ($r$) for the DG methods with GLegQ, GLobQ and MixQ for the 2D problem with non-constant velocity using Cartesian meshes.</caption><thead><tr><th>p</th><th>h</th><th>ε<sub>GLegQ</sub></th><th>r<sub>GLegQ</sub></th><th>ε<sub>GLobQ</sub></th><th>r<sub>GLobQ</sub></th><th>ε<sub>MixQ</sub></th><th>r<sub>MixQ</sub></th><th>ε<sub>GLobQ</sub><br>ε<sub>GLegQ</sub></th><th>ε<sub>MixQ</sub><br>ε<sub>GLegQ</sub></th></tr></thead><tbody><tr><td>1</td><td>0.5000</td><td>3.5 × 10<sup>-1</sup></td><td>–</td><td>3.5 × 10<sup>-1</sup></td><td>–</td><td>3.5 × 10<sup>-1</sup></td><td>–</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.2500</td><td>1.5 × 10<sup>-1</sup></td><td>1.2</td><td>1.5 × 10<sup>-1</sup></td><td>1.2</td><td>1.5 × 10<sup>-1</sup></td><td>1.2</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.1250</td><td>4.4 × 10<sup>-2</sup></td><td>1.8</td><td>4.4 × 10<sup>-2</sup></td><td>1.8</td><td>4.4 × 10<sup>-2</sup></td><td>1.8</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.0625</td><td>1.2 × 10<sup>-2</sup></td><td>1.9</td><td>1.2 × 10<sup>-2</sup></td><td>1.9</td><td>1.2 × 10<sup>-2</sup></td><td>1.9</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.0313</td><td>3.0 × 10<sup>-3</sup></td><td>2.0</td><td>3.2 × 10<sup>-3</sup></td><td>1.9</td><td>3.0 × 10<sup>-3</sup></td><td>2.0</td><td>1.1</td><td>1.0</td></tr><tr><td>1</td><td>0.0156</td><td>7.7 × 10<sup>-4</sup></td><td>2.0</td><td>9.2 × 10<sup>-4</sup></td><td>1.8</td><td>7.7 × 10<sup>-4</sup></td><td>2.0</td><td>1.2</td><td>1.0</td></tr><tr><td>1</td><td>0.0078</td><td>1.9 × 10<sup>-4</sup></td><td>2.0</td><td>2.6 × 10<sup>-4</sup></td><td>1.8</td><td>1.9 × 10<sup>-4</sup></td><td>2.0</td><td>1.3</td><td>1.0</td></tr><tr><td>2</td><td>0.5000</td><td>4.0 × 10<sup>-2</sup></td><td>–</td><td>4.0 × 10<sup>-2</sup></td><td>–</td><td>4.0 × 10<sup>-2</sup></td><td>–</td><td>1.0</td><td>1.0</td></tr><tr><td>2</td><td>0.2500</td><td>1.2 × 10<sup>-2</sup></td><td>1.7</td><td>1.2 × 10<sup>-2</sup></td><td>1.7</td><td>1.2 × 10<sup>-2</sup></td><td>1.7</td><td>1.0</td><td>1.0</td></tr><tr><td>2</td><td>0.1250</td><td>1.6 × 10<sup>-3</sup></td><td>2.9</td><td>1.7 × 10<sup>-3</sup></td><td>2.9</td><td>1.6 × 10<sup>-3</sup></td><td>2.9</td><td>1.0</td><td>1.0</td></tr><tr><td>2</td><td>0.0625</td><td>2.1 × 10<sup>-4</sup></td><td>2.9</td><td>2.3 × 10<sup>-4</sup></td><td>2.8</td><td>2.1 × 10<sup>-4</sup></td><td>2.9</td><td>1.1</td><td>1.0</td></tr><tr><td>2</td><td>0.0313</td><td>2.2 × 10<sup>-5</sup></td><td>3.2</td><td>4.6 × 10<sup>-5</sup></td><td>2.3</td><td>2.2 × 10<sup>-5</sup></td><td>3.2</td><td>2.0</td><td>1.0</td></tr><tr><td>2</td><td>0.0156</td><td>2.8 × 10<sup>-6</sup></td><td>3.0</td><td>5.6 × 10<sup>-6</sup></td><td>3.0</td><td>2.7 × 10<sup>-6</sup></td><td>3.0</td><td>2.0</td><td>1.0</td></tr><tr><td>2</td><td>0.0078</td><td>3.5 × 10<sup>-7</sup></td><td>3.0</td><td>7.3 × 10<sup>-7</sup></td><td>2.9</td><td>3.4 × 10<sup>-7</sup></td><td>3.0</td><td>2.1</td><td>1.0</td></tr><tr><td>3</td><td>0.5000</td><td>2.2 × 10<sup>-2</sup></td><td>–</td><td>2.2 × 10<sup>-2</sup></td><td>–</td><td>2.2 × 10<sup>-2</sup></td><td>–</td><td>1.0</td><td>1.0</td></tr><tr><td>3</td><td>0.2500</td><td>1.0 × 10<sup>-3</sup></td><td>4.4</td><td>1.1 × 10<sup>-3</sup></td><td>4.4</td><td>1.1 × 10<sup>-3</sup></td><td>4.4</td><td>1.0</td><td>1.0</td></tr><tr><td>3</td><td>0.1250</td><td>8.5 × 10<sup>-5</sup></td><td>3.6</td><td>9.1 × 10<sup>-5</sup></td><td>3.5</td><td>8.2 × 10<sup>-5</sup></td><td>3.7</td><td>1.1</td><td>1.0</td></tr><tr><td>3</td><td>0.0625</td><td>6.4 × 10<sup>-6</sup></td><td>3.7</td><td>8.7 × 10<sup>-6</sup></td><td>3.4</td><td>6.3 × 10<sup>-6</sup></td><td>3.7</td><td>1.4</td><td>1.0</td></tr><tr><td>3</td><td>0.0313</td><td>4.0 × 10<sup>-7</sup></td><td>4.0</td><td>8.0 × 10<sup>-7</sup></td><td>3.5</td><td>4.0 × 10<sup>-7</sup></td><td>4.0</td><td>2.0</td><td>1.0</td></tr><tr><td>3</td><td>0.0156</td><td>2.3 × 10<sup>-8</sup></td><td>4.2</td><td>5.2 × 10<sup>-8</sup></td><td>3.9</td><td>2.3 × 10<sup>-8</sup></td><td>4.1</td><td>2.3</td><td>1.0</td></tr><tr><td>3</td><td>0.0078</td><td>1.2 × 10<sup>-9</sup></td><td>4.2</td><td>2.9 × 10<sup>-9</sup></td><td>4.2</td><td>1.2 × 10<sup>-9</sup></td><td>4.2</td><td>2.3</td><td>1.0</td></tr><tr><td>4</td><td>0.5000</td><td>8.8 × 10<sup>-4</sup></td><td>–</td><td>7.5 × 10<sup>-4</sup></td><td>–</td><td>7.6 × 10<sup>-4</sup></td><td>–</td><td>0.8</td><td>0.9</td></tr><tr><td>4</td><td>0.2500</td><td>9.6 × 10<sup>-5</sup></td><td>3.2</td><td>9.5 × 10<sup>-5</sup></td><td>3.0</td><td>8.9 × 10<sup>-5</sup></td><td>3.1</td><td>1.0</td><td>0.9</td></tr><tr><td>4</td><td>0.1250</td><td>3.6 × 10<sup>-6</sup></td><td>4.7</td><td>3.9 × 10<sup>-6</sup></td><td>4.6</td><td>3.2 × 10<sup>-6</sup></td><td>4.8</td><td>1.1</td><td>0.9</td></tr><tr><td>4</td><td>0.0625</td><td>1.2 × 10<sup>-7</sup></td><td>4.9</td><td>1.9 × 10<sup>-7</sup></td><td>4.3</td><td>1.1 × 10<sup>-7</sup></td><td>4.9</td><td>1.6</td><td>0.9</td></tr><tr><td>4</td><td>0.0313</td><td>3.8 × 10<sup>-9</sup></td><td>5.0</td><td>8.0 × 10<sup>-9</sup></td><td>4.6</td><td>3.3 × 10<sup>-9</sup></td><td>5.0</td><td>2.1</td><td>0.9</td></tr><tr><td>4</td><td>0.0156</td><td>1.2 × 10<sup>-10</sup></td><td>5.0</td><td>2.6 × 10<sup>-10</sup></td><td>5.0</td><td>1.0 × 10<sup>-10</sup></td><td>5.0</td><td>2.2</td><td>0.9</td></tr><tr><td>4</td><td>0.0078</td><td>3.5 × 10<sup>-12</sup></td><td>5.1</td><td>8.4 × 10<sup>-12</sup></td><td>4.9</td><td>3.3 × 10<sup>-12</sup></td><td>5.0</td><td>2.4</td><td>0.9</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 7:** Influence, in the case of non-constant velocity and non-constant jacobian, of time step on the $L^2(\Omega)$-error at $T = 1$ for different polynomial orders, meshes and method.</figcaption></figure>

proposed approach provides the same accuracy as the DG method with GLegQ, with a lower cost induced by the matrix-vector products. In addition, the leading term of the cost for the DG methods with GLegQ and MixQ is identical but the proposed approach provides results twice as accurate. For non-constant velocity, the efficiency of GLobQ and MixQ seems to be equivalent, because the improved stability of the GLobQ method compensates for the improved precision of
<table><caption>**Table 8:** Relative error in the $L^2(\Omega)$ norm ($\epsilon$), and rate of convergence ($r$) for the DG methods with GLegQ, GLobQ and MixQ for the 2D problem with non-constant velocity using non-Cartesian meshes.</caption><thead><tr><th>p</th><th>h</th><th>ε<sub>GLegQ</sub></th><th>r<sub>GLegQ</sub></th><th>ε<sub>GLobQ</sub></th><th>r<sub>GLobQ</sub></th><th>ε<sub>MixQ</sub></th><th>r<sub>MixQ</sub></th><th>ε<sub>GLobQ</sub><br>ε<sub>GLegQ</sub></th><th>ε<sub>MixQ</sub><br>ε<sub>GLegQ</sub></th></tr></thead><tbody><tr><td>1</td><td>0.5737</td><td>3.6 × 10<sup>-1</sup></td><td>–</td><td>3.6 × 10<sup>-1</sup></td><td>–</td><td>3.6 × 10<sup>-1</sup></td><td>–</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.2991</td><td>1.7 × 10<sup>-1</sup></td><td>1.2</td><td>1.7 × 10<sup>-1</sup></td><td>1.2</td><td>1.7 × 10<sup>-1</sup></td><td>1.2</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.1699</td><td>4.9 × 10<sup>-2</sup></td><td>2.2</td><td>4.9 × 10<sup>-2</sup></td><td>2.2</td><td>4.9 × 10<sup>-2</sup></td><td>2.2</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.0871</td><td>1.3 × 10<sup>-2</sup></td><td>2.0</td><td>1.3 × 10<sup>-2</sup></td><td>2.0</td><td>1.3 × 10<sup>-2</sup></td><td>2.0</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.0440</td><td>3.3 × 10<sup>-3</sup></td><td>2.0</td><td>3.5 × 10<sup>-3</sup></td><td>1.9</td><td>3.3 × 10<sup>-3</sup></td><td>2.0</td><td>1.1</td><td>1.0</td></tr><tr><td>1</td><td>0.0224</td><td>8.6 × 10<sup>-4</sup></td><td>2.0</td><td>1.0 × 10<sup>-3</sup></td><td>1.8</td><td>8.6 × 10<sup>-4</sup></td><td>2.0</td><td>1.2</td><td>1.0</td></tr><tr><td>1</td><td>0.0111</td><td>2.2 × 10<sup>-4</sup></td><td>2.0</td><td>3.2 × 10<sup>-4</sup></td><td>1.7</td><td>2.2 × 10<sup>-4</sup></td><td>1.9</td><td>1.5</td><td>1.0</td></tr><tr><td>2</td><td>0.5737</td><td>4.7 × 10<sup>-2</sup></td><td>–</td><td>4.7 × 10<sup>-2</sup></td><td>–</td><td>4.7 × 10<sup>-2</sup></td><td>–</td><td>1.0</td><td>1.0</td></tr><tr><td>2</td><td>0.2991</td><td>1.4 × 10<sup>-2</sup></td><td>1.9</td><td>1.4 × 10<sup>-2</sup></td><td>1.9</td><td>1.4 × 10<sup>-2</sup></td><td>1.9</td><td>1.0</td><td>1.0</td></tr><tr><td>2</td><td>0.1699</td><td>2.0 × 10<sup>-3</sup></td><td>3.4</td><td>2.1 × 10<sup>-3</sup></td><td>3.3</td><td>2.0 × 10<sup>-3</sup></td><td>3.4</td><td>1.0</td><td>1.0</td></tr><tr><td>2</td><td>0.0871</td><td>2.7 × 10<sup>-4</sup></td><td>3.0</td><td>2.9 × 10<sup>-4</sup></td><td>2.9</td><td>2.7 × 10<sup>-4</sup></td><td>3.0</td><td>1.1</td><td>1.0</td></tr><tr><td>2</td><td>0.0440</td><td>3.5 × 10<sup>-5</sup></td><td>3.0</td><td>6.0 × 10<sup>-5</sup></td><td>2.3</td><td>3.5 × 10<sup>-5</sup></td><td>3.0</td><td>1.7</td><td>1.0</td></tr><tr><td>2</td><td>0.0224</td><td>4.1 × 10<sup>-6</sup></td><td>3.2</td><td>8.2 × 10<sup>-6</sup></td><td>2.9</td><td>4.1 × 10<sup>-6</sup></td><td>3.2</td><td>2.0</td><td>1.0</td></tr><tr><td>2</td><td>0.0111</td><td>5.2 × 10<sup>-7</sup></td><td>2.9</td><td>1.1 × 10<sup>-6</sup></td><td>2.9</td><td>5.3 × 10<sup>-7</sup></td><td>2.9</td><td>2.0</td><td>1.0</td></tr><tr><td>3</td><td>0.5737</td><td>2.4 × 10<sup>-2</sup></td><td>–</td><td>2.4 × 10<sup>-2</sup></td><td>–</td><td>2.4 × 10<sup>-2</sup></td><td>–</td><td>1.0</td><td>1.0</td></tr><tr><td>3</td><td>0.2991</td><td>1.5 × 10<sup>-3</sup></td><td>4.3</td><td>1.5 × 10<sup>-3</sup></td><td>4.3</td><td>1.4 × 10<sup>-3</sup></td><td>4.3</td><td>1.0</td><td>1.0</td></tr><tr><td>3</td><td>0.1699</td><td>1.2 × 10<sup>-4</sup></td><td>4.5</td><td>1.2 × 10<sup>-4</sup></td><td>4.4</td><td>1.1 × 10<sup>-4</sup></td><td>4.5</td><td>1.1</td><td>1.0</td></tr><tr><td>3</td><td>0.0871</td><td>9.0 × 10<sup>-6</sup></td><td>3.8</td><td>1.2 × 10<sup>-5</sup></td><td>3.5</td><td>8.9 × 10<sup>-6</sup></td><td>3.8</td><td>1.3</td><td>1.0</td></tr><tr><td>3</td><td>0.0440</td><td>5.7 × 10<sup>-7</sup></td><td>4.0</td><td>1.0 × 10<sup>-6</sup></td><td>3.5</td><td>5.7 × 10<sup>-7</sup></td><td>4.0</td><td>1.8</td><td>1.0</td></tr><tr><td>3</td><td>0.0224</td><td>3.5 × 10<sup>-8</sup></td><td>4.1</td><td>7.5 × 10<sup>-8</sup></td><td>3.9</td><td>3.5 × 10<sup>-8</sup></td><td>4.1</td><td>2.1</td><td>1.0</td></tr><tr><td>3</td><td>0.0111</td><td>2.1 × 10<sup>-9</sup></td><td>4.0</td><td>4.4 × 10<sup>-9</sup></td><td>4.0</td><td>2.1 × 10<sup>-9</sup></td><td>4.0</td><td>2.1</td><td>1.0</td></tr><tr><td>4</td><td>0.5737</td><td>1.7 × 10<sup>-3</sup></td><td>–</td><td>1.5 × 10<sup>-3</sup></td><td>–</td><td>1.5 × 10<sup>-3</sup></td><td>–</td><td>0.9</td><td>0.9</td></tr><tr><td>4</td><td>0.2991</td><td>1.3 × 10<sup>-4</sup></td><td>4.0</td><td>1.2 × 10<sup>-4</sup></td><td>3.8</td><td>1.2 × 10<sup>-4</sup></td><td>3.9</td><td>1.0</td><td>0.9</td></tr><tr><td>4</td><td>0.1699</td><td>5.8 × 10<sup>-6</sup></td><td>5.4</td><td>6.5 × 10<sup>-6</sup></td><td>5.2</td><td>5.4 × 10<sup>-6</sup></td><td>5.4</td><td>1.1</td><td>0.9</td></tr><tr><td>4</td><td>0.0871</td><td>2.1 × 10<sup>-7</sup></td><td>5.0</td><td>2.9 × 10<sup>-7</sup></td><td>4.7</td><td>1.9 × 10<sup>-7</sup></td><td>5.0</td><td>1.4</td><td>0.9</td></tr><tr><td>4</td><td>0.0440</td><td>7.4 × 10<sup>-9</sup></td><td>4.9</td><td>1.4 × 10<sup>-8</sup></td><td>4.4</td><td>6.9 × 10<sup>-9</sup></td><td>4.8</td><td>1.9</td><td>0.9</td></tr><tr><td>4</td><td>0.0224</td><td>2.4 × 10<sup>-10</sup></td><td>5.1</td><td>4.5 × 10<sup>-10</sup></td><td>5.1</td><td>2.2 × 10<sup>-10</sup></td><td>5.1</td><td>1.9</td><td>0.9</td></tr><tr><td>4</td><td>0.0111</td><td>7.9 × 10<sup>-12</sup></td><td>4.9</td><td>1.5 × 10<sup>-11</sup></td><td>4.9</td><td>7.7 × 10<sup>-12</sup></td><td>4.8</td><td>1.9</td><td>1.0</td></tr></tbody></table>

the MixQ.

## 5.2 Three dimensional examples

The first three dimensional example considers the propagation of a sinusoidal wave at constant velocity in a Cartesian mesh. The domain is $\Omega = [0, 1]^3$ and the velocity field is given by $a = [\cos(\pi/6) \sin(\pi/4), \sin(\pi/6) \sin(\pi/4), \cos(\pi/4)]$. As in the previous examples, the final time is $T = 1$ and the time step is selected small enough so that the error induced by the time marching scheme is negligible when compared to the error induced by the spatial discretisation.

A set of Cartesian meshes is utilised in three dimensions. The mesh corresponding to the $i$-th level of refinement is denoted by $M^i$ and contains $2^i \times 2^i \times 2^i$ hexahedral elements with characteristic element size $h = 2^{-i}$, as shown in Figure 8.

Table 9 shows the relative error, measured in the $L^2(\Omega)$ norm and denoted by $\epsilon$, for the DG methods with GLegQ and GLobQ for a polynomial approximation of order ranging from $p = 1$ up to $p = 4$. The table also shows the rates of convergence, denoted by $r$ and the ratio between the error of both DG approaches.

The optimal rate of convergence is shown for both DG approaches, that is a rate of approximately $p+1$, in the $L^2(\Omega)$ norm. As in the two dimensional example, the extra accuracy provided by the DG method with GLegQ is clearly observed. For the lowest order case, the use of GLegQ
<figure><img src="image_10.png" /><figcaption>Figure 8: Cartesian meshes of hexahedral elements.</figcaption></figure>

<table><caption>Table 9: Relative error in the $L^2(\Omega)$ norm ($\epsilon$), and rate of convergence ($r$) for the DG methods with GLegQ and GLobQ for the 3D problem with constant velocity using Cartesian meshes.</caption><thead><tr><th>p</th><th>h</th><th>ε<sub>GLegQ</sub> = ε<sub>MixQ</sub></th><th>r<sub>GLegQ</sub> = r<sub>MixQ</sub></th><th>ε<sub>GLobQ</sub></th><th>r<sub>GLobQ</sub></th><th>ε<sub>GLobQ</sub>/ε<sub>GLegQ</sub></th></tr></thead><tbody><tr><td>1</td><td>0.5000</td><td>3.1 × 10<sup>-1</sup></td><td>-</td><td>8.6 × 10<sup>-1</sup></td><td>-</td><td>2.8</td></tr><tr><td>1</td><td>0.2500</td><td>9.8 × 10<sup>-2</sup></td><td>1.6</td><td>4.7 × 10<sup>-1</sup></td><td>0.9</td><td>4.8</td></tr><tr><td>1</td><td>0.1250</td><td>2.5 × 10<sup>-2</sup></td><td>2.0</td><td>1.6 × 10<sup>-1</sup></td><td>1.6</td><td>6.4</td></tr><tr><td>1</td><td>0.0625</td><td>6.1 × 10<sup>-3</sup></td><td>2.0</td><td>4.3 × 10<sup>-2</sup></td><td>1.9</td><td>7.1</td></tr><tr><td>2</td><td>0.5000</td><td>5.3 × 10<sup>-2</sup></td><td>-</td><td>2.1 × 10<sup>-1</sup></td><td>-</td><td>3.9</td></tr><tr><td>2</td><td>0.2500</td><td>6.2 × 10<sup>-3</sup></td><td>3.1</td><td>3.0 × 10<sup>-2</sup></td><td>2.8</td><td>4.9</td></tr><tr><td>2</td><td>0.1250</td><td>7.6 × 10<sup>-4</sup></td><td>3.0</td><td>3.8 × 10<sup>-3</sup></td><td>3.0</td><td>5.0</td></tr><tr><td>2</td><td>0.0625</td><td>9.5 × 10<sup>-5</sup></td><td>3.0</td><td>4.8 × 10<sup>-4</sup></td><td>3.0</td><td>5.1</td></tr><tr><td>3</td><td>0.5000</td><td>6.3 × 10<sup>-3</sup></td><td>-</td><td>2.6 × 10<sup>-2</sup></td><td>-</td><td>4.1</td></tr><tr><td>3</td><td>0.2500</td><td>3.9 × 10<sup>-4</sup></td><td>4.0</td><td>1.7 × 10<sup>-3</sup></td><td>3.9</td><td>4.4</td></tr><tr><td>3</td><td>0.1250</td><td>2.5 × 10<sup>-5</sup></td><td>4.0</td><td>1.1 × 10<sup>-4</sup></td><td>4.0</td><td>4.5</td></tr><tr><td>3</td><td>0.0625</td><td>1.5 × 10<sup>-6</sup></td><td>4.0</td><td>6.9 × 10<sup>-6</sup></td><td>4.0</td><td>4.5</td></tr><tr><td>4</td><td>0.5000</td><td>6.7 × 10<sup>-4</sup></td><td>-</td><td>2.6 × 10<sup>-3</sup></td><td>-</td><td>3.9</td></tr><tr><td>4</td><td>0.2500</td><td>2.1 × 10<sup>-5</sup></td><td>5.0</td><td>8.6 × 10<sup>-5</sup></td><td>4.9</td><td>4.2</td></tr><tr><td>4</td><td>0.1250</td><td>6.5 × 10<sup>-7</sup></td><td>5.0</td><td>2.7 × 10<sup>-6</sup></td><td>5.0</td><td>4.2</td></tr><tr><td>4</td><td>0.0625</td><td>2.0 × 10<sup>-8</sup></td><td>5.0</td><td>8.5 × 10<sup>-8</sup></td><td>5.0</td><td>4.2</td></tr></tbody></table>

provides an error almost nine times lower than that produced by the DG method with GLobQ. For higher orders, the use of GLegQ consistently provides an accuracy between four to five times higher than the DG method with GLobQ. It is worth noticing that for $p > 1$ the improvement due to the use of GLegQ is the same as for the two dimensional example with constant velocity and using Cartesian meshes.

In terms of the computational cost, the DG method with GLegQ requires $O(6p^4)$ and $O(3p^3)$ operations for computing the matrix-vector products corresponding to the convection and face contributions to the residual, respectively. The DG method with GLobQ requires $O(9p^4)$ and $O(3p^2)$ operations for computing the matrix-vector products corresponding to the convection and face contributions to the residual, respectively. Therefore, it is clear that the approach with GLegQ is beneficial as the dominant term, corresponding to the convection term, has a lower order of operations and provides more accurate results. As in two dimensions, note that stability tests (not shown here because the results are very similar in three dimensions) indicate that part of the efficiency of the DG method with GLegQ is lost because a time step twice smaller needs to be used.

As done in two dimensions, the next three dimensional example considers a problem with constant velocity by using non-Cartesian meshes. The meshes considered are shown in Figure 9 and are, again, generated by randomly perturbing the internal nodes of the Cartesian meshes employed in the previous three dimensional example.
<figure><img src="image_11.png" /><figcaption>Figure 9: Non-Cartesian meshes of hexahedral elements.</figcaption></figure>

Table 10 shows the relative error, measured in the $L^2(\Omega)$ norm, for the DG methods with GLegQ, GLobQ and MixQ. The table also shows the rates of convergence and the ratio between the error of the DG methods with GLobQ and MixQ with respect to the GLegQ.

<table><caption>Table 10: Relative error in the $L^2(\Omega)$ norm ($\epsilon$), and rate of convergence ($r$) for the DG methods with GLegQ, GLobQ and MixQ for the 3D problem with constant velocity using non-Cartesian meshes.</caption><thead><tr><th>p</th><th>h</th><th>ε<sub>GLegQ</sub></th><th>r<sub>GLegQ</sub></th><th>ε<sub>GLobQ</sub></th><th>r<sub>GLobQ</sub></th><th>ε<sub>MixQ</sub></th><th>r<sub>MixQ</sub></th><th>ε<sub>GLobQ</sub>/ε<sub>GLegQ</sub></th><th>ε<sub>MixQ</sub>/ε<sub>GLegQ</sub></th></tr></thead><tbody><tr><td>1</td><td>0.5750</td><td>3.1 × 10<sup>-1</sup></td><td>–</td><td>8.6 × 10<sup>-1</sup></td><td>–</td><td>3.1 × 10<sup>-1</sup></td><td>–</td><td>2.8</td><td>1.0</td></tr><tr><td>1</td><td>0.3482</td><td>1.0 × 10<sup>-1</sup></td><td>2.2</td><td>4.9 × 10<sup>-1</sup></td><td>1.1</td><td>1.1 × 10<sup>-1</sup></td><td>2.0</td><td>4.8</td><td>1.1</td></tr><tr><td>1</td><td>0.1758</td><td>2.7 × 10<sup>-2</sup></td><td>1.9</td><td>1.9 × 10<sup>-1</sup></td><td>1.4</td><td>3.6 × 10<sup>-2</sup></td><td>1.6</td><td>6.8</td><td>1.3</td></tr><tr><td>1</td><td>0.0897</td><td>6.8 × 10<sup>-3</sup></td><td>2.1</td><td>5.5 × 10<sup>-2</sup></td><td>1.8</td><td>1.3 × 10<sup>-2</sup></td><td>1.5</td><td>8.1</td><td>1.9</td></tr><tr><td>2</td><td>0.5750</td><td>6.5 × 10<sup>-2</sup></td><td>–</td><td>2.4 × 10<sup>-1</sup></td><td>–</td><td>6.8 × 10<sup>-2</sup></td><td>–</td><td>3.6</td><td>1.0</td></tr><tr><td>2</td><td>0.3482</td><td>7.4 × 10<sup>-3</sup></td><td>4.3</td><td>3.4 × 10<sup>-2</sup></td><td>3.9</td><td>7.8 × 10<sup>-3</sup></td><td>4.3</td><td>4.6</td><td>1.1</td></tr><tr><td>2</td><td>0.1758</td><td>9.9 × 10<sup>-4</sup></td><td>2.9</td><td>4.8 × 10<sup>-3</sup></td><td>2.9</td><td>1.1 × 10<sup>-3</sup></td><td>2.8</td><td>4.9</td><td>1.1</td></tr><tr><td>2</td><td>0.0897</td><td>1.2 × 10<sup>-4</sup></td><td>3.1</td><td>6.0 × 10<sup>-4</sup></td><td>3.1</td><td>1.5 × 10<sup>-4</sup></td><td>2.9</td><td>4.9</td><td>1.2</td></tr><tr><td>3</td><td>0.5750</td><td>6.5 × 10<sup>-3</sup></td><td>–</td><td>2.6 × 10<sup>-2</sup></td><td>–</td><td>6.5 × 10<sup>-3</sup></td><td>–</td><td>4.0</td><td>1.0</td></tr><tr><td>3</td><td>0.3482</td><td>4.6 × 10<sup>-4</sup></td><td>5.3</td><td>2.0 × 10<sup>-3</sup></td><td>5.1</td><td>4.8 × 10<sup>-4</sup></td><td>5.2</td><td>4.3</td><td>1.0</td></tr><tr><td>3</td><td>0.1758</td><td>3.3 × 10<sup>-5</sup></td><td>3.9</td><td>1.5 × 10<sup>-4</sup></td><td>3.8</td><td>3.6 × 10<sup>-5</sup></td><td>3.8</td><td>4.5</td><td>1.1</td></tr><tr><td>3</td><td>0.0897</td><td>2.2 × 10<sup>-6</sup></td><td>4.0</td><td>1.0 × 10<sup>-5</sup></td><td>4.0</td><td>2.6 × 10<sup>-6</sup></td><td>3.9</td><td>4.6</td><td>1.2</td></tr><tr><td>4</td><td>0.5750</td><td>7.9 × 10<sup>-4</sup></td><td>–</td><td>3.2 × 10<sup>-3</sup></td><td>–</td><td>8.0 × 10<sup>-4</sup></td><td>–</td><td>4.0</td><td>1.0</td></tr><tr><td>4</td><td>0.3482</td><td>2.8 × 10<sup>-5</sup></td><td>6.6</td><td>1.2 × 10<sup>-4</sup></td><td>6.6</td><td>2.9 × 10<sup>-5</sup></td><td>6.6</td><td>4.2</td><td>1.0</td></tr><tr><td>4</td><td>0.1758</td><td>9.9 × 10<sup>-7</sup></td><td>4.9</td><td>4.2 × 10<sup>-6</sup></td><td>4.9</td><td>1.0 × 10<sup>-6</sup></td><td>4.9</td><td>4.2</td><td>1.0</td></tr><tr><td>4</td><td>0.0897</td><td>3.4 × 10<sup>-8</sup></td><td>5.0</td><td>1.4 × 10<sup>-7</sup></td><td>5.0</td><td>3.6 × 10<sup>-8</sup></td><td>5.0</td><td>4.2</td><td>1.1</td></tr></tbody></table>

The results reveal that the proposed DG method with MixQ offers again almost identical results when compared to the DG method with GLegQ. This indicates that the numerical integration error introduced by the MixQ approach is below the discretisation error even when three dimensional elements with non-constant Jacobian are considered. In addition, the results in Table 10 indicate that the DG methods with GLegQ and MixQ offer between four and five times more accurate results when compared to the approach with GLobQ. This indicates a higher gain than what was observed in the two dimensional example.

In terms of the computational cost, the DG method with GLegQ requires $O(2p^6)$ and $O(2p^5)$ operations for computing the matrix-vector products corresponding to the convection and face contributions to the residual, respectively. The DG method with GLobQ requires $O(9p^4)$ and $O(3p^2)$ operations for computing the matrix-vector products corresponding to the convection and face contributions to the residual, respectively. With the proposed DG method with MixQ $O(9p^4)$ and $O(2p^5)$ operations are required to compute the matrix-vector products of the convection and face contributions to the residual, respectively. It is worth noting that contrary to the observations made in two dimensions, for the DG method with MixQ the dominant term is now associated to the face term. However, when compared to the DG method with GLegQ, the proposed DG
method still has a leading cost lower and provides the same accuracy. As in two dimensions, larger time steps need to be used so that the method is overall two to three times more efficient.

The last example considers a deformational flow in three dimensions. The non-constant velocity field is given by

$$ a(\mathbf{x}, t) = \cos(\pi t/T) \begin{bmatrix} 0.1 \sin^2(\pi x_1) \sin(2\pi x_2) \sin(\pi x_3) \cos(\pi x_3) \\ 0.2 \sin(2\pi x_1) \sin^2(\pi x_2) \sin(\pi x_3) \cos(\pi x_3) \\ 0.05\pi \sin(2\pi x_1) \sin(2\pi x_2) \sin^2(\pi x_3) \end{bmatrix}^T \quad (62) $$

and the final time is $T = 1$.

The initial condition is given by $u_0(\mathbf{x}) = 0.5 + 0.5 \sin(2\pi x_1) \sin(2\pi x_2) \sin(2\pi x_3)$ and periodic boundary conditions are considered on the boundary of the domain $\partial\Omega$. Analogously to the two dimensional deformation flow considered in Section 5.1, the cosine term in Eq. (62) ensures that the flow reverses at time $t = T/2$ and reaches the same state as the initial condition at time $t = T$.

Table 11 shows the relative error, measured in the $L^2(\Omega)$ norm, for the three DG approaches considered using Cartesian meshes.

<table><caption>**Table 11:** Relative error in the $L^2(\Omega)$ norm ($\epsilon$), and rate of convergence ($r$) for the DG methods with GLegQ, GLobQ and MixQ for the 3D problem with non-constant velocity using Cartesian meshes.</caption><thead><tr><th>p</th><th>h</th><th>ε<sub>GLegQ</sub></th><th>r<sub>GLegQ</sub></th><th>ε<sub>GLobQ</sub></th><th>r<sub>GLobQ</sub></th><th>ε<sub>MixQ</sub></th><th>r<sub>MixQ</sub></th><th>ε<sub>GLobQ</sub>/ε<sub>GLegQ</sub></th><th>ε<sub>MixQ</sub>/ε<sub>GLegQ</sub></th></tr></thead><tbody><tr><td>1</td><td>0.5000</td><td>2.3 × 10<sup>-1</sup></td><td>–</td><td>2.3 × 10<sup>-1</sup></td><td>–</td><td>2.3 × 10<sup>-1</sup></td><td>–</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.2500</td><td>1.5 × 10<sup>-1</sup></td><td>0.6</td><td>1.5 × 10<sup>-1</sup></td><td>0.6</td><td>1.5 × 10<sup>-1</sup></td><td>0.6</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.1250</td><td>4.8 × 10<sup>-2</sup></td><td>1.7</td><td>4.8 × 10<sup>-2</sup></td><td>1.7</td><td>4.8 × 10<sup>-2</sup></td><td>1.7</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.0625</td><td>1.3 × 10<sup>-2</sup></td><td>1.9</td><td>1.3 × 10<sup>-2</sup></td><td>1.9</td><td>1.3 × 10<sup>-2</sup></td><td>1.9</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.0313</td><td>3.3 × 10<sup>-3</sup></td><td>2.0</td><td>3.4 × 10<sup>-3</sup></td><td>1.9</td><td>3.3 × 10<sup>-3</sup></td><td>2.0</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.0156</td><td>8.4 × 10<sup>-4</sup></td><td>2.0</td><td>9.3 × 10<sup>-4</sup></td><td>1.9</td><td>8.4 × 10<sup>-4</sup></td><td>2.0</td><td>1.1</td><td>1.0</td></tr><tr><td>1</td><td>0.0078</td><td>2.1 × 10<sup>-4</sup></td><td>2.0</td><td>2.5 × 10<sup>-4</sup></td><td>1.9</td><td>2.1 × 10<sup>-4</sup></td><td>2.0</td><td>1.2</td><td>1.0</td></tr><tr><td>2</td><td>0.5000</td><td>3.8 × 10<sup>-2</sup></td><td>–</td><td>3.8 × 10<sup>-2</sup></td><td>–</td><td>3.8 × 10<sup>-2</sup></td><td>–</td><td>1.0</td><td>1.0</td></tr><tr><td>2</td><td>0.2500</td><td>1.2 × 10<sup>-2</sup></td><td>1.7</td><td>1.2 × 10<sup>-2</sup></td><td>1.7</td><td>1.2 × 10<sup>-2</sup></td><td>1.7</td><td>1.0</td><td>1.0</td></tr><tr><td>2</td><td>0.1250</td><td>1.5 × 10<sup>-3</sup></td><td>3.0</td><td>1.5 × 10<sup>-3</sup></td><td>3.0</td><td>1.5 × 10<sup>-3</sup></td><td>3.0</td><td>1.0</td><td>1.0</td></tr><tr><td>2</td><td>0.0625</td><td>1.9 × 10<sup>-4</sup></td><td>3.0</td><td>2.2 × 10<sup>-4</sup></td><td>2.8</td><td>1.8 × 10<sup>-4</sup></td><td>3.0</td><td>1.2</td><td>1.0</td></tr><tr><td>2</td><td>0.0313</td><td>2.2 × 10<sup>-5</sup></td><td>3.1</td><td>3.6 × 10<sup>-5</sup></td><td>2.6</td><td>2.2 × 10<sup>-5</sup></td><td>3.1</td><td>1.6</td><td>1.0</td></tr><tr><td>2</td><td>0.0156</td><td>2.6 × 10<sup>-6</sup></td><td>3.1</td><td>5.1 × 10<sup>-6</sup></td><td>2.8</td><td>2.6 × 10<sup>-6</sup></td><td>3.1</td><td>1.9</td><td>1.0</td></tr><tr><td>2</td><td>0.0078</td><td>3.2 × 10<sup>-7</sup></td><td>3.0</td><td>6.7 × 10<sup>-7</sup></td><td>2.9</td><td>3.1 × 10<sup>-7</sup></td><td>3.0</td><td>2.1</td><td>1.0</td></tr><tr><td>3</td><td>0.5000</td><td>2.2 × 10<sup>-2</sup></td><td>–</td><td>2.2 × 10<sup>-2</sup></td><td>–</td><td>2.2 × 10<sup>-2</sup></td><td>–</td><td>1.0</td><td>1.0</td></tr><tr><td>3</td><td>0.2500</td><td>9.6 × 10<sup>-4</sup></td><td>4.5</td><td>9.7 × 10<sup>-4</sup></td><td>4.5</td><td>9.6 × 10<sup>-4</sup></td><td>4.5</td><td>1.0</td><td>1.0</td></tr><tr><td>3</td><td>0.1250</td><td>7.6 × 10<sup>-5</sup></td><td>3.7</td><td>8.2 × 10<sup>-5</sup></td><td>3.6</td><td>7.3 × 10<sup>-5</sup></td><td>3.7</td><td>1.1</td><td>1.0</td></tr><tr><td>3</td><td>0.0625</td><td>5.4 × 10<sup>-6</sup></td><td>3.8</td><td>7.6 × 10<sup>-6</sup></td><td>3.4</td><td>5.4 × 10<sup>-6</sup></td><td>3.8</td><td>1.4</td><td>1.0</td></tr><tr><td>3</td><td>0.0313</td><td>3.6 × 10<sup>-7</sup></td><td>3.9</td><td>6.6 × 10<sup>-7</sup></td><td>3.5</td><td>3.6 × 10<sup>-7</sup></td><td>3.9</td><td>1.8</td><td>1.0</td></tr><tr><td>3</td><td>0.0156</td><td>2.2 × 10<sup>-8</sup></td><td>4.0</td><td>4.6 × 10<sup>-8</sup></td><td>3.8</td><td>2.2 × 10<sup>-8</sup></td><td>4.0</td><td>2.1</td><td>1.0</td></tr><tr><td>3</td><td>0.0078</td><td>1.2 × 10<sup>-9</sup></td><td>4.1</td><td>2.8 × 10<sup>-9</sup></td><td>4.1</td><td>1.2 × 10<sup>-9</sup></td><td>4.1</td><td>2.2</td><td>1.0</td></tr><tr><td>4</td><td>0.5000</td><td>8.2 × 10<sup>-4</sup></td><td>–</td><td>7.1 × 10<sup>-4</sup></td><td>–</td><td>7.2 × 10<sup>-4</sup></td><td>–</td><td>0.9</td><td>0.9</td></tr><tr><td>4</td><td>0.2500</td><td>8.9 × 10<sup>-5</sup></td><td>3.2</td><td>9.0 × 10<sup>-5</sup></td><td>3.0</td><td>8.3 × 10<sup>-5</sup></td><td>3.1</td><td>1.0</td><td>0.9</td></tr><tr><td>4</td><td>0.1250</td><td>3.3 × 10<sup>-6</sup></td><td>4.7</td><td>3.7 × 10<sup>-6</sup></td><td>4.6</td><td>2.9 × 10<sup>-6</sup></td><td>4.9</td><td>1.1</td><td>0.9</td></tr><tr><td>4</td><td>0.0625</td><td>1.2 × 10<sup>-7</sup></td><td>4.8</td><td>1.7 × 10<sup>-7</sup></td><td>4.4</td><td>1.0 × 10<sup>-7</sup></td><td>4.8</td><td>1.4</td><td>0.8</td></tr><tr><td>4</td><td>0.0313</td><td>3.8 × 10<sup>-9</sup></td><td>5.0</td><td>6.9 × 10<sup>-9</sup></td><td>4.7</td><td>3.2 × 10<sup>-9</sup></td><td>5.0</td><td>1.8</td><td>0.8</td></tr><tr><td>4</td><td>0.0156</td><td>1.2 × 10<sup>-10</sup></td><td>5.0</td><td>2.5 × 10<sup>-10</sup></td><td>4.8</td><td>1.0 × 10<sup>-10</sup></td><td>5.0</td><td>2.1</td><td>0.9</td></tr></tbody></table>

As observed in previous examples, all methods show the optimal rate of convergence but the DG method with GLobQ requires very fine meshes to exhibit the asymptotic convergence. In this example the DG methods with GLegQ and MixQ offer an error two times lower than the
DG method with GLobQ. The results of the proposed approach, with MixQ, are in some cases
marginally more accurate than the DG method with GLegQ, but this is attributed by the inexact
integration of the complex velocity field considered.

It is worth noting that in three dimensions the dominant cost associated to the proposed
approach is of lower order when compared to the DG method with GLegQ but the accuracy is
almost identical. When comparing the proposed DG method with MixQ and GLobQ, the cost of
the MixQ is only slightly higher but the extra accuracy is clearly observed. However, as in two
dimensions, for non-constant velocity, it seems that the overall efficiency is the same for DG
methods with MixQ and GLegQ because the weaker stability compensates for the lower cost for
computing matrix-vector products.

6 Concluding remarks

A novel DG method for the solution of linear conservation laws has been presented. The approach uses a Gauss-Lobatto distribution of nodes to define the polynomial approximation and a mixture of Gauss-Lobatto and Gauss-Legendre quadratures for the numerical integration. The convection part of the residual of the system of ordinary differential equations is evaluated using Gauss-Lobatto quadratures, whereas the term corresponding to the face integrals is evaluated using Gauss-Legendre quadratures.

A formal proof is provided to show that in the case of constant Jacobian and constant coefficients (i.e., velocity field), the proposed method provides exactly the same results as a standard DG method with Gauss-Legendre quadratures, despite both the mass and convection matrices are underintegrated. For the general case, the numerical results show that the accuracy is also preserved, meaning that the numerical integration error is below the spatial discretisation error. This brings the advantages of both commonly used DG methods with Gauss-Lobatto and Gauss-Legendre quadratures. The proposed approach is as accurate as the standard DG method with Gauss-Legendre quadratures but at a lower cost, which is only slightly more expensive than the DG method with Gauss-Lobatto quadratures in three dimensions and as expensive in two dimensions. This conclusion remains generally correct when adding stability to the discussion, but, in some cases, the improved stability of the DG method with Gauss-Lobatto quadratures compensates for the improved precision of the method with mixed quadratures.

References

Ainsworth, M. (2004a). "Dispersive and dissipative behavior of high order discontinuous Galerkin finite element methods". *Journal of Computational Physics* 198.1, pp. 106–130. doi: 10.1016/j.jcp.2004.01.004

Ainsworth, M. (2004b). "Dispersive and dissipative behaviour of high order discontinuous Galerkin finite element methods". *J. Comput. Phys.* 198.1, pp. 106–130

Ainsworth, M. and H. Abdul Wajid (2009). "Dispersive and dissipative behavior of the Spectral Element Method". *SIAM Journal of Numerical Analysis* 47.5, pp. 3910–3937. doi: 10.1137/080724976

Antonietti, P. F., I. Mazzieri, A. Quarteroni, and F. Rapetti (2012). “Non-conforming high order approximations of the elastodynamics equation”. *Computer Methods in Applied Mechanics and Engineering* 209-212, pp. 212–238. doi: 10.1016/j.cma.2011.11.004

Arnold, D. N. (1982). "An interior penalty finite element method with discontinuous elements". *SIAM J. Numer. Anal.* 19.4, pp. 742–760

Arnold, D. N., F. Brezzi, B. Cockburn, and L. D. Marini (2002). “Unified analysis of discontinuous Galerkin methods for elliptic problems”. *SIAM J. Numer. Anal.* 39.5, pp. 1749–1779

Atkins, H. L. and C. W. Shu (1998a). “Quadrature-free implementation of discontinuous Galerkin method for hyperbolic equations”. *AIAA Journal* 36.5, pp. 775–782

Atkins, H. L. and C. W. Shu (1998b). “Quadrature-free implementation of discontinuous Galerkin method for hyperbolic equations”. *AIAA Journal* 36.5, pp. 775–782

Balanis, C. A. (1989). *Advanced Engineering Electromagnetics*. New York: John Wiley and Sons
Beck, A. D., T. Bolemann, D. Flad, H. Frank, G. J. Gassner, F. Hindenlang, and C.-D. Munz (2014). "High-order Discontinuous Galerkin Spectral Element Methods for transitional and turbulent flow simulations". *International Journal for Numerical Methods in Fluids* 76.8, pp. 522–548. doi: 10.1002/fld.3943

Bey, K. S. and J. T. Oden (1996). "hp-version discontinuous Galerkin methods for hyperbolic conservation laws". *Computer Methods in Applied Mechanics and Engineering* 133.3-4, pp. 259–286. doi: 10.1016/0045-7825(95)00944-2

Brezzi, F., G. Manzini, D. Marini, P. Pietra, and A. Russo (2000). "Discontinuous Galerkin approximations for elliptic problems". *Numer. Methods Partial Differential Equations* 16.4, pp. 365–378

Bui-Thanh, T. and O. Ghattas (2012). "Analysis of an hp-nonconforming Discontinuous Galerkin Spectral Element Method for wave propagation". *SIAM Journal of Numerical Analysis* 50.3, pp. 1801–1826. doi: 10.1137/110828010

Canuto, C., M. Y. Hussaini, A. Quarteroni, and T. A. Zang (2007). *Spectral methods: evolution to complex geometries and applications to fluid dynamics*. Springer

Castel, N., G. Cohen, and M. Duruflé (2009). "Application of Discontinuous Galerkin spectral method on hexahedral elements for aeroacoustic". *Journal of Computational Acoustics* 17.2, pp. 175–196. doi: 10.1142/S0218396X09003914

Castillo, P. (2002). "Performance of discontinuous Galerkin methods for elliptic PDEs". *SIAM J. Sci. Comp.* 24.2, pp. 524–547

Castillo, P. (2003). "A superconvergence result for Discontinuous Galerkin methods applied to elliptic problems". *Comput. Methods Appl. Mech. Engrg.* 192.41-42, pp. 4675–4685

Castillo, P. (2005). "An a posteriori error estimate for the local Discontinuous Galerkin method". *J. Sci. Comp.* 22.1, pp. 187–204

Castillo, P., B. Cockburn, I. Perugia, and D. Schötzau (2000). "An a priori error analysis of the local discontinuous Galerkin method for elliptic problems". *SIAM J. Numer. Anal.* 38.5, pp. 1676–1706

Castillo, P., B. Cockburn, D. Schotzau, and C. Schwab (2002). "Optimal a priori error estimates for the hp-version of the local Discontinuous Galerkin method for convection-diffusion problems". *Math. Comp.* 71.238, pp. 455–478

Chavent, G. and G. Salzano (1982). "A finite element method for the 1D water flooding problem with gravity". *J. Comput. Phys.* 45.3, pp. 307–344

Chavent, G. and B. Cockburn (1989). "The local projection $P^0 P^1$-Discontinuous Galerkin finite element method for scalar conservation laws". *M$^2$AN* 23.4, pp. 565–592

Chen, Z., B. Cockburn, J. W. Jerome, and C. W. Shu (1995). "Mixed-RKDG finite element methods for the 2-D hydrodynamics model for semiconductor device simulation". *VLSI Design* 3.2, pp. 145–158

Cockburn, B. (1999). "Discontinuous Galerkin methods for convection-dominated problems". *High-Order methods for computational physics*. Ed. by T. J. Barth and H. Deconinck. Vol. 9. Lecture Notes in Computational Science and Engineering. Springer, pp. 69–224. doi: 10.1007/978-3-662-03882-6_2

Cockburn, B. (2003). "Discontinuous Galerkin methods. Plenary lecture presented at the 80th Annual GAMM Conference, Augsburg, 25-28 March 2002". *Z. Angew. Math. Mech.* 83.11, pp. 731–754

Cockburn, B. and C. W. Shu (1989a). "TVB Runge-Kutta Local Projection discontinuous Galerkin Finite Element Method For Conservation-Laws II. General Framework". *Math. Comp.* 52.186, pp. 411–435

Cockburn, B., S. Y. Lin, and C. W. Shu (1989b). "TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite-Element Method For Conservation-Laws III. One-Dimensional Systems". *J. Comput. Phys.* 84.1, pp. 90–113

Cockburn, B., S. C. Hou, and C. W. Shu (1990). "The Runge-Kutta Local Projection Discontinuous Galerkin Finite-Element Method For Conservation-Laws IV. The Multidimensional Case". *Math. Comp.* 54.190, pp. 545–581

Cockburn, B. and C. W. Shu (1991a). "The Runge-Kutta local projection $P^1$-Discontinuous Galerkin method for scalar conservation laws". *M$^2$AN* 25, pp. 337–361
Cockburn, B. and C. W. Shu (1991b). "The Runge-Kutta Local Projection Rho-1-Discontinuous Galerkin Finite-Element Method For Scalar Conservation-Laws". *RAIRO-Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique* 25.3, pp. 337-361

Cockburn, B. and C. W. Shu (1998a). "The local Discontinuous Galerkin method for time-dependent convection-diffusion systems". *SIAM J. Numer. Anal.* 35.6, pp. 2440-2463

Cockburn, B. and C. W. Shu (1998b). "The Runge-Kutta Discontinuous Galerkin method for conservation laws V. Multidimensional systems". *J. Comput. Phys.* 141.2, pp. 199-224

Cockburn, B., G. Kanschat, D. Schötzau, and C. Schwab (2002). "Local discontinuous Galerkin methods for the Stokes system". *SIAM J. Numer. Anal.* 40, pp. 319-343

Cockburn, B., G. Karniadakis, and C. Shu (2005). *Discontinuous Galerkin methods. Theory, computation and applications*. Springer

Cockburn, B. (1998). "An introduction to the discontinuous Galerkin method for convection-dominated problems". *Advanced numerical approximation of nonlinear hyperbolic equations*, pp. 151-268

Cockburn, B., G. E. Karniadakis, and C.-W. Shu (2000). "The development of discontinuous Galerkin methods". *Discontinuous Galerkin Methods*. Springer, pp. 3-50

Delcourte, S. and N. Glinsky (2015). "Analysis of a high-order space and time Discontinuous Galerkin Method for elastodynamic equations. Applications to 3D wave propagation". *ESAIM: Mathematical Modelling and Numerical Analysis* 49.4, pp. 1085-1126. doi: 10.1051/m2an/2015001

Deville, M. O., P. F. Fischer, P. F. Fischer, E. Mund, et al. (2002). *High-order methods for incompressible fluid flow*. Vol. 9. Cambridge university press

Donea, J. and A. Huerta (2005). *Finite Element Methods for Flow Problems*. Wiley

Flad, D., A. D. Beck, and P. Guthke (2020). "A large eddy simulation method for DGSEM using non-linearly optimized relaxation filters". *Journal of Computational Physics* 408.109303, pp. 1-12. doi: 10.1016/j.jcp.2020.109303

Friedrich, L., A. R. Winters, D. C. Del Rey Fernández, G. J. Gassner, M. Parsani, and M. H. Carpenter (2018). "An entropy stable h/p non-conforming Discontinuous Galerkin Method with the summation-by-parts property". *Journal of Scientific Computing* 77, pp. 689-725. doi: 10.1007/s10915-018-0733-7

Gassner, G. J. and D. A. Kopriva (2011). "A comparison of the dispersion and dissipation errors of Gauss and Gauss-Lobatto Discontinuous Galerkin Spectral Element Methods". *SIAM Journal of Scientific Computing* 33.5, pp. 2560-2579. doi: 10.1137/100807211

Gassner, G. J. and A. D. Beck (2013). "On the accuracy of high-order discretizations for underresolved turbulence simulations". *Theoretical and Computational Fluid Dynamics* 27, pp. 221-237. doi: 10.1007/s00162-011-0253-7

Giraldo, F. X., J. S. Hesthaven, and T. Warburton (2002). "Nodal high-order Discontinuous Galerkin methods for the spherical shallow water equations". *Journal of Computational Physics* 181.2, pp. 499-525. doi: 10.1006/jcph.2002.7139

Giraldo, F. X. and M. Restelli (2008). "A study of spectral element and discontinuous Galerkin methods for the Navier-Stokes equations in nonhydrostatic mesoscale atmospheric modeling: equation sets and test cases". *Journal of Computational Physics* 227.8, pp. 3849-3877. doi: 10.1016/j.jcp.2007.12.009

Givoli, D. (1992). *Numerical Methods for Problems in Infinite Domains*. Amsterdam: Elsevier

Hairer, E., C. Lubich, and M. Roche (2006). *The numerical solution of differential-algebraic systems by Runge-Kutta methods*. Vol. 1409. Springer

Hesthaven, J. S. and T. Warburton (2002). "Nodal high-order methods on unstructured grids. I. Time-domain solution of Maxwell's equations". *Journal of Computational Physics* 181.1, pp. 186-221

Hesthaven, J. S. and T. Warburton (2008a). "Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications". Vol. 54. Texts in Applied Mathematics. New York: Springer

Hesthaven, J. S. and T. Warburton (2008b). *Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications*. Vol. 54. Texts in Applied Mathematics. New York: Springer
Hesthaven, J. and T Warburton (2004). "Discontinuous Galerkin methods for the time-domain Maxwell's equations". *ACES Newsletter* 19, pp. 10–29

Houston, P., C. Schwab, and E. Süli (2002). "Discontinuous hp-Finite Element Methods for advection-diffusion-reaction problems". *SIAM Journal of Numerical Analysis* 39.6, pp. 2133–2163. DOI: 10.1137/S0036142900374111

Hu, F. Q., M. Y. Hussaini, and P. Rasetarinera (1999). "An analysis of the Discontinuous Galerkin Method for wave propagation problems". *Journal of Computational Physics* 151.2, pp. 921–946. DOI: 10.1006/jcph.1999.6227

Jameson, A., W. Schmidt, and E. Turkel (1981). "Numerical solution of the Euler equations using Runge-Kutta time-stepping schemes". *Proceedings of the 14th AIAA Fluid and Plasma Dynamic Conference*. AIAA. Palo Alto, California

Johnson, C. and J. Pitkaranta (1986). "An Analysis Of The discontinuous Galerkin Method For A Scalar Hyperbolic Equation". *Math. Comp.* 46.173, pp. 1–26

Kabakian, A. V., V. Shankar, and W. F. Hall (2004). "Unstructured grid-based discontinuous Galerkin method for broadband electromagnetic simulations". *Journal of Scientific Computing* 20.3, pp. 405–431

Kent, J. (2019). "Capturing the cross-terms in multidimensional advection schemes". *International Journal for Numerical Methods in Fluids* 91.2, pp. 49–62

Kopriva, D. A. and G. Gassner (2010). "On the quadrature and weak form choices in collocation type discontinuous Galerkin spectral element methods". *SIAM Journal of Scientific Computing* 44, pp. 136–155. DOI: 10.1007/s10915-010-9372-3

Laskowski, W., A. M. Rueda-Ramirez, G. Rubio, E. Valero, and E. Ferrer (2020). "Advantages of static condensation in implicit compressible Navier Stokes DGSEM solvers". *Computers & Fluids* 209.104646, pp. 1–17. DOI: 10.1016/j.compfluid.2020.104646

Lesaint, P. and P. A. Raviart (1974). "On a finite element method for solving the neutron transport equation". *Mathematical aspects of finite elements in PDE*. Ed. by C. de Boor. Academic Press, pp. 89–123

LeVeque, R. J. (1992). *Numerical methods for conservation laws*. Second. Lectures in Mathematics ETH Zürich. Basel: Birkhäuser Verlag, pp. x+214

LeVeque, R. J. (1996). "High-resolution conservative algorithms for advection in incompressible flow". *SIAM Journal on Numerical Analysis* 33.2, pp. 627–665

LeVeque, R. J. (2002). *Finite volume methods for hyperbolic problems*. Cambridge university press

Li, B. Q. (2006). *Discontinuous Finite Elements in Fluid Dynamics and Heat Transfer (Computational Fluid and Solid Mechanics)*. Springer

Lilienthal, M. (2015). "Error controlled hp-adaptive Finite Element Methods for the time-dependent Maxwell equations". PhD thesis. Germany: Technische Universität Darmstadt

Malm, J., P. Schlatter, P. F. Fischer, and D. S. Henningson (2013). "Stabilization of the spectral element method in convection dominated flows by recovery of skew-symmetry". *Journal of Scientific Computing* 57, pp. 254–277

Noventa, G., F. Massa, F. Bassi, A. Colombo, N. Franchina, and A. Ghidoni (2016). "A high-order Discontinuous Galerkin solver for unsteady incompressible turbulent flows". *Computers & Fluids* 139, pp. 248–260

Osher, S. (1984). "Riemann Solvers, The Entropy Condition, And Difference Approximations". *SIAM J. Numer. Anal.* 21.2, pp. 217–235

Peraire, J. and P.-O. Persson (2008). "The Compact Discontinuous Galerkin (CDG) Method for Elliptic Problems". *SIAM J. Sci. Comput.* 30.4, pp. 1806–1824

Persson, P.-O. and J. Peraire (2006a). "An efficient low memory implicit DG algorithm for time dependent problems". *Proceedings of the 44th AIAA Aerospace Sciences Meeting and Exhibit*. AIAA

Persson, P.-O. and J. Peraire (2006b). "Sub-cell shock capturing for discontinuous Galerkin Methods". *Proceedings of the 44th AIAA Aerospace Sciences Meeting and Exhibit*. AIAA

Persson, P.-O. and J. Peraire (2009). "Curved Mesh Generation and Mesh Refinement using Lagrangian Solid Mechanics". *Proceedings of the 47th AIAA Aerospace Sciences Meeting and Exhibit*. AIAA
Peterson, T. (1991). "A Note on the Convergence of the Discontinuous Galerkin Method for a Scalar Hyperbolic Equation". *SIAM J. Numer. Anal.* 28.1, pp. 133-140

Pomponet Oliveira, S. and S. Angelozi Leite (2018). "Error analysis of the spectral element method with Gauss-Lobatto-Legendre points for the acoustic wave equation in heterogeneous media". *Appl. Numer. Math.* 129, pp. 39-57. DOI: 10.1016/j.apnum.2018.02.007

Rasetarinera, P., M. Y. Hussaini, and F. Q. Hu (2000). "Some remarks on the accuracy of a Discontinuous Galerkin Method". *Discontinuous Galerkin methods. Theory, computation and applications*. Ed. by B. Cockburn, G. Karniadakis, and C. W. Shu. Springer-Verlag, pp. 407-412

Reed, W. H. and T. R. Hill (1973). *Triangular mesh methods for the neutron transport equation*. Tech. rep. LA-UR-73-479. Los Alamos: Los Alamos Scientific Laboratory

Renac, F. (2019). "Entropy stable DGSEM for nonlinear hyperbolic systems in nonconservative form with application to two-phase flows". *Journal of Computational Physics* 382, pp. 1-26. DOI: 10.1016/j.jcp.2018.12.035

Rider, W. J. and R. B. Lowrie (2002). "The use of classical Lax-Friedrichs Riemann solvers with discontinuous Galerkin methods". *Int. J. Numer. Meth. Engng. Fluids* 40.3-4, pp. 479-486

Sármány, D., M. A. Botchev, and J. J. W. van der Vegt (2007). "Dispersion and Dissipation Error in High-Order Runge-Kutta Discontinuous Galerkin Discretisations of the Maxwell Equations". *J. Sci. Comp.* 33.1, pp. 47-74

Sevilla, R., S. Fernández-Méndez, and A. Huerta (2011). "Comparison of high-order curved finite elements". *International Journal for Numerical Methods in Engineering* 87.8, pp. 719-734

Sevilla, R., O. Hassan, and K. Morgan (2014). "The Use of Hybrid Meshes to Improve the Efficiency of a Discontinuous Galerkin Method for the Solution of Maxwell's Equations". *Computers & Structures* 137, pp. 2-13

Sherwin, S. (2000). "Dispersion analysis of the continuous and Discontinuous Galerkin formulations". *Discontinuous Galerkin methods. Theory, computation and applications*. Ed. by B. Cockburn, G. Karniadakis, and C. W. Shu. Springer-Verlag, pp. 425-431

Sherwin, S. J., R. M. Kirby, J. Peiro, L. R. Taylor, and O. C. Zienkiewicz (2006). "On 2D elliptic discontinuous Galerkin methods". *Int. J. Numer. Meth. Engng.* 65.5, pp. 752-784

Shu, C. W. (2003). "High-order finite difference and finite volume WENO schemes and Discontinuous Galerkin methods for CFD". *Int. J. Comput. Fluid Dyn.* 17.2, pp. 107-118

Terrana, S., J.-P. Vilotte, and L. Guillot (2018). "A spectral hybridizable Discontinuous Galerkin Method for elastic-acoustic wave propagation". *Geophysical Journal International* 213.1, pp. 574-602. DOI: 10.1093/gji/ggx557

Ven, H. V. der and J. J. W. V. der Vegt (2002). "Space-time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flows II. Efficient flux quadrature". *Comput. Methods Appl. Mech. Engrg.* 191.41-42, pp. 4747-4780

Warburton, T. C., I. Lomtev, Y. Du, S. J. Sherwin, and G. E. Karniadakis (1999). "Galerkin and discontinuous Galerkin spectral/hp methods". *Computer Methods in Applied Mechanics and Engineering* 175.3-4, pp. 343-359. DOI: 10.1016/S0045-7825(98)00360-0

Wilcox, L. C., G. Stadler, C. Burstedde, and O. Ghattas (2010). "A high-order Discontinuous Galerkin Method for wave propagation through coupled elastic-acoustic media". *Journal of Computational Physics* 229.24, pp. 9373-9396. DOI: 10.1016/j.jcp.2010.09.008

Winters, A. R. and D. A. Kopriva (2014). "ALE-DGSEM approximation of wave reflection and transmission from a moving medium". *Journal of Computational Physics* 263, pp. 233-267. DOI: 10.1016/j.jcp.2014.01.022

Zienkiewicz, O. C. and R. L. Taylor (2000). *The Finite Element Method*. Fifth. Vol. 1. The basis. Butterwort-Heinemann

Zienkiewicz, O. C., R. L. Taylor, S. J. Sherwin, and J. Peiro (2003). "On discontinuous Galerkin methods". *Int. J. Numer. Meth. Engng.* 58.8, pp. 1119-1148

Zienkiewicz, O. C., R. L. Taylor, O. C. Zienkiewicz, and R. L. Taylor (1977). *The finite element method*. Vol. 3. McGraw-hill London