# Nonlinear advection-diffusion equation: ADER-DG penalty vs. relaxation schemes

Afaf Bouharguane, Angelo Iollo, Alexis Tardieu

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

## To cite this version:

Afaf Bouharguane, Angelo Iollo, Alexis Tardieu. Nonlinear advection-diffusion equation: ADER-DG penalty vs. relaxation schemes. 2025. hal-05243830

HAL Id: hal-05243830

https://inria.hal.science/hal-05243830v1

Preprint submitted on 7 Sep 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Nonlinear advection-diffusion equation:
ADER-DG penalty vs. relaxation schemes

Afaf Bouharguane<sup>1,2,3</sup>, Angelo Iollo<sup>1,2,3</sup> and Alexis Tardieu<sup>1,2,3*</sup>

<sup>1</sup>Université de Bordeaux, Bordeaux, France.

<sup>2</sup>Centre Inria de l'Université de Bordeaux, Memphis team, Talence, France.

<sup>3</sup>Institut de Mathématiques de Bordeaux, UMR CNRS 5251, Talence, France.

*Corresponding author(s). E-mail(s):
alexis.tardieu@u-bordeaux.fr;

Contributing authors: afaf.bouharguane@u-bordeaux.fr; angelo.iollo@u-bordeaux.fr;

**Abstract**

This paper proposes to solve numerically the two dimensional nonlinear advection-diffusion equation. The space discretization relies on a classical Discontinuous Galerkin (DG) method. This scheme is combined together with an Arbitrary high order DERivatives (ADER) approach to ensure the same high order of accuracy in time compared to the precision in space. More precisely, two different methods are compared regarding the computational cost, the error and the order of convergence: the Symmetric Interior Penalty Galerkin (SIPG) and the Cattaneo relaxation methods. The viscosity of the medium, the mesh and the approximation degree being fixed, we aim at determining whether the penalty or the relaxation scheme is to be preferred. Numerical examples are provided to illustrate and quantify this comparison. We show that both approaches ensure to reach an arbitrary high precision and present an interest from an implementation perspective.

**Keywords:** Nonlinear advection-diffusion, discontinuous Galerkin, ADER approach, penalty method, relaxation method, arbitrary high order schemes.
## 2 ADER-DG schemes for nonlinear advection-diffusion

# 1 Introduction

In mathematical modelling and numerical simulations for the study of partial differential equations (PDEs), many discretization approaches are used, like the Finite Differences (FD), Finite Volumes (FV) or Finite Elements (FE) schemes. For computational fluid dynamics (CFD) applications, one of the main issues is to find a method that can deal with both the advective and diffusive components of the problem under study or, in other words, that can apply to both hyperbolic and elliptic or parabolic operators. In particular, when a high order approximation is required with respect to the space and time dimensions, a great care must be taken when combining the various techniques in order to preserve the effectiveness of each of them.

For the space discretization, if FD schemes rely on local Taylor expansions of the solution around points, FV and FE schemes are known as variational methods, which means that the PDE under study is considered under an integral form, and so is the desired numerical solution. On the one hand, FV schemes [1–5] generally consist of a local integral formulation of the problem, written in each cell of the mesh. The average of the numerical solution in each cell evolves in time under the effect of numerical fluxes exchanged through its interfaces with neighbouring cells or the boundary conditions. On the other hand, FE schemes consist in expanding the numerical solution in a continuous polynomial basis of arbitrary high degree, leading to an arbitrary high accuracy of the method. A hybridization of the FV and FE approaches is the Discontinuous Galerkin (DG) method [6–11], which is based on the expansion of the solution with respect to a discontinuous polynomial basis, locally in each cell. The information is finally recombined between neighbouring cells through well defined numerical fluxes.

For the time discretization, a FD, FV or FE approximation at an arbitrary high order in space is generally combined to the classical Rünge-Kutta (RK) scheme [2, 7] to reach the same order of accuracy in time, which is based on a few iterations of the same scheme at different times. With a completely different approach, the Arbitrary high order DERivatives (ADER) framework [1, 2, 11–21] consists in finding a cell by cell high order space-time approximation in a single step, called a prediction, that will then be corrected to take into account the neighbourhood of the cell. The resulting prediction-correction scheme is then abbreviated ADER-FV or ADER-DG method for instance, if a space FV or DG approximation is chosen respectively.

The ADER combinations provide satisfying results on a wide range of time evolution problems for CFD applications. In particular, the ADER-DG scheme is characterized by a compact stencil and a good capacity to reach a space-time arbitrary high order of accuracy. For instance, the ADER-FV setting has already proven to be efficient on time evolving Chimera meshes for advection-diffusion [1, 18] or incompressible Navier-Stokes equations. With a view to extend this work to nonconforming hierarchical meshes like quadtrees or octrees in two or three dimensions respectively [22–25], the ADER-DG framework thus shows promising perspectives.
In this paper, we propose to compare the Cattaneo relaxation scheme [1, 18] to the Symmetric Interior Penalty Galerkin (SIPG) method [6, 11, 26–29], which is the symmetric variant of the Interior Penalty Discontinuous Galerkin (IPDG) approach. The symmetric version, SIPG, has proven to be always faster and of optimal order of accuracy compared to its incomplete and non-symmetric counterparts, IIPG and NIPG [11]. Both the Cattaneo and SIPG methods are used within the global framework of the ADER-DG space-time prediction-correction scheme and differ in the way the diffusive term of a given PDE is treated. The former consists in rewriting this term under a hyperbolic form, while the latter is derived from a penalty approach to recover the continuity provided by this term. The difference in their approaches results in differences in terms of order of convergence, execution time and asymptotic behaviour. We thus propose to illustrate this comparison on the two dimensional nonlinear advection-diffusion equation and give some theoretical and numerical results.

Section 2 presents all the notations and the space-time tools used in this paper. In section 3, we propose an ADER-DG prediction-correction method with two variants: the SIPG penalty and the Cattaneo relaxation approaches, and we detail the numerical fluxes used at cells interfaces. In section 4, a numerical application is performed and carefully examined. The last section 5 is dedicated to some conclusions and perspectives of the present paper.

# 2 The nonlinear advection-diffusion equation

## 2.1 The governing equation

Let us consider the two dimensional unsteady nonlinear advection-diffusion of $u: \Omega \times [0, T_f] \to \mathbb{R}^2$, where $\Omega \subset \mathbb{R}^2$ is the space domain, and the two components of all the vector quantities are denoted by $\phi := [\phi_x, \phi_y]^T$. The medium is characterized by its kinematic viscosity $\nu: \Omega \times [0, T_f] \to \mathbb{R}_+^*$. In the following, $S: \Omega \times [0, T_f] \to \mathbb{R}^2$ is the source term. We also denote by $F, G$ respectively the physical fluxes in directions $x$ and $y$ that can be expressed as the sum $F := F_a + F_d$ and $G := G_a + G_d$ of their advective and diffusive contributions defined below:

$$F_a(u) := u_x u, \quad F_d(u) := -\nu \partial_x u, \quad G_a(u) := u_y u, \quad G_d(u) := -\nu \partial_y u. \quad (1)$$

The nonlinear advection-diffusion equation then reads:

$$\partial_t u + \nabla \cdot (u \otimes u - \nu \nabla u) = S, \quad (2)$$

where $u \otimes u := u \cdot u^T$ is a square matrix of size 2, and so is $\nabla u$, the gradient of $u$. By using the flux notations, this equation also writes:

$$\partial_t u + \partial_x F(u) + \partial_y G(u) = S. \quad (3)$$
4 ADER-DG schemes for nonlinear advection-diffusion

This time evolution problem is completed with the initial data $u(x, y, 0) = u_0(x, y)$ at time $t = 0$, as well as boundary conditions enforced on the boundary $\partial\Omega$ of the space domain. In this paper, the conditions imposed on $\partial\Omega$ can be of Dirichlet type: $u|_{\partial\Omega} = h_D$, or Neumann type: $H_d(u)|_{\partial\Omega} = h_N$ for $H = F, G$, where $h$ is a given function.

## 2.2 From the physical to reference space-time framework

When using a Galerkin method in order to solve a PDE, whether continuous or discontinuous, it can be helpful to transform all the polynomial expansions, the basis functions or the quadrature rules for numerical integration from the physical coordinates $(x, y, t)$ to reference coordinates $(\xi, \eta, \tau)$. In the present work, we consider a cartesian grid to discretize the space domain $\Omega := [x^-, x^+] \times [y^-, y^+] \subset \mathbb{R}^2$ composed of $N_x N_y$ uniform cells $\Omega_{i,j} := \omega_i^x \times \omega_j^y$ of size $\Delta x \Delta y$. We also define $h := \max(\Delta x, \Delta y)$. The time domain $[0, T_f]$ is made up of $N_t$ intervals $T^n$ of a priori different size $\Delta t^n$. We thus have:

$$ \Omega := \bigcup_{i=1}^{N_x} \bigcup_{j=1}^{N_y} \Omega_{i,j}, \quad [0, T_f] := \bigcup_{n=0}^{N_t-1} T^n, \qquad (4) $$

where we have introduced the notations:

$$ \omega_i^x := [x_{i-1}, x_i], \quad x_i := x^- + i\Delta x, \qquad (5a) $$

$$ \omega_j^y := [y_{j-1}, y_j], \quad y_j := y^- + j\Delta y, \qquad (5b) $$

$$ T^n := [t^n, t^{n+1}], \quad t^n := \Delta t^0 + \dots + \Delta t^{n-1}. \qquad (5c) $$

Our aim is to determine a mapping between any of these three physical variables $(x, y, t)$ to the three reference variables $(\xi, \eta, \tau)$, so that the physical space-time cell $c_{i,j}^n := \Omega_{i,j} \times T^n$ is mapped to the reference space-time cell $\hat{c} := \hat{\omega}^3$, where $\hat{\omega} := [0, 1]$. Consequently, a physical quantity $\phi := [\phi_x, \phi_y]^T$ defined on $c_{i,j}^n$ has a counterpart $\hat{\phi} := [\hat{\phi}_x, \hat{\phi}_y]^T$ defined in $\hat{c}$. This transformation is illustrated on figure 1.

<figure><img src="image_3.png" /><figcaption>Fig. 1 Mapping between the physical and reference frameworks.</figcaption></figure>
A simple choice is to consider the linear bijections $\mathcal{X}_i, \mathcal{Y}_j, \mathcal{T}_n$ below:

$$x = \mathcal{X}_i(\xi) := x_i + \xi \Delta x, \qquad (6a)$$

$$y = \mathcal{Y}_j(\eta) := y_j + \eta \Delta y, \qquad (6b)$$

$$t = \mathcal{T}_n(\tau) := t^n + \tau \Delta t^n. \qquad (6c)$$

In the reference space-time cell $\hat{c}$, the nonlinear advection-diffusion equation (3) then becomes:

$$\partial_{\tau} \hat{u} + \alpha^n \partial_{\xi} \hat{F}(\hat{u}) + \beta^n \partial_{y} \hat{G}(\hat{u}) = \Delta t^n \hat{S}, \qquad (7)$$

by applying the chain rule for the derivative of composite functions, where $\alpha^n := \Delta t^n/\Delta x$ and $\beta^n := \Delta t^n/\Delta y$. The quantities follow the definitions and relations detailed in subsection 2.1 but, for the sake of clarity, we detail the expressions of the advective and diffusive fluxes that become:

$$\hat{F}_a(\hat{u}) := \hat{u}_x \hat{u}, \quad \hat{F}_d(\hat{u}) := -\frac{\nu}{\Delta x} \partial_{\xi} \hat{u}, \quad \hat{G}_a(\hat{u}) := \hat{u}_y \hat{u}, \quad \hat{G}_d(\hat{u}) := -\frac{\nu}{\Delta y} \partial_{\eta} \hat{u}. \qquad (8)$$

The compact forms (3) and (7) of the equation (see [1, 2, 11, 12] for instance) will be used in the rest of this paper, respectively when considering the physical or reference frameworks.

## 2.3 Lagrange basis functions

Contrary to continuous Galerkin methods, the DG approach provides a greater scope to choose the location of the degrees of freedom of the numerical solution. Within the context of a global ADER-DG space-time numerical scheme, both space and time variables are treated equally from a discretization point of vue. The resulting three dimensional space-time structure is obtained as a tensor product of the one dimensional space and time structures.

Let us introduce $\mathbb{P}_D$ the vector space of polynomials of degree $d \le D$, with $D \in \mathbb{N}$, and note $\mathbb{N}(A) := \{1, \dots, A\}$ for all $A \in \mathbb{N}^*$ as well as $S_m := (D+1)^m$ for $m = 1, 2, 3$. In our work, we choose to use the 1D Lagrange basis functions $\vartheta_1, \dots, \vartheta_{S_1} \in \mathbb{P}_D$ passing through the 1D Gauss-Legendre quadrature points $\pi_1, \dots, \pi_{S_1}$. These polynomials are defined as:

$$\vartheta_r(v) := \prod_{k \neq r} \left( \frac{v - \pi_k}{\pi_r - \pi_k} \right), \quad v = \xi, \eta, \tau, \quad r \in \mathbb{N}(S_1). \qquad (9)$$

The 2D and 3D Gauss-Legendre quadrature points are built as the tensor product of the 1D quadrature points previously introduced. Then, the 2D and 3D Lagrange polynomials $\theta_1, \dots, \theta_{S_2}$ and $\Theta_1, \dots, \Theta_{S_3}$ are respectively defined as products of the 1D basis functions:

$$\theta_s(v, v') := \vartheta_k(v)\vartheta_{k'}(v'), \quad v, v' = \xi, \eta, \tau, \quad s \in \mathbb{N}(S_2), \qquad (10)$$
6 ADER-DG schemes for nonlinear advection-diffusion

and:

$$ \Theta_{\ell}(\xi, \eta, \tau) := \vartheta_k(\xi)\vartheta_{k'}(\eta)\vartheta_{k''}(\tau), \quad \ell \in \mathbb{N}(S_3). \quad (11) $$

There are many possibilities for the bijective mappings between local indices $k, k', k''$ and the 2D global index $r$ in (10) or the 3D global index $\ell$ in (11). Here, we propose the following orderings:

$$ \mathbf{m}: \begin{cases} \mathbb{N}(S_1)^2 \to \mathbb{N}(S_2) \\ (k, k') \mapsto r := S_1(k' - 1) + k' \end{cases}, \quad (12) $$

and:

$$ \mathbf{M}: \begin{cases} \mathbb{N}(S_1)^3 \to \mathbb{N}(S_3) \\ (k, k', k'') \mapsto \ell := S_2(k'' - 1) + \mathbf{m}(k, k') \end{cases}. \quad (13) $$

In particular, by applying the function $\phi \mapsto \log(\phi)$ to definition (9), deriving the expression with respect to the variable $v = \xi, \eta, \tau$ and using the usual properties on derivation, one can easily obtain explicitly the first derivative of the $r$-th 1D Lagrange basis function as:

$$ \vartheta'_{r}(v) := \vartheta_{r}(v) \sum_{k \neq r} \left( \frac{1}{v - \pi_{k}} \right), \quad v = \xi, \eta, \tau, \quad r \in \mathbb{N}(S_{1}), \quad (14) $$

for any arbitrary value of $D \in \mathbb{N}$. As the 2D and 3D basis functions have separated variables, their derivative with respect to one variable is immediate, for instance: $\partial_{\xi}\Theta_{\ell}(\xi, \eta, \tau) = \vartheta'_{k}(\xi)\vartheta_{k'}(\eta)\vartheta_{k''}(\tau)$. There are three main advantages of using these basis functions:

* Their implementation simply requires the knowledge of the ($D+1$) quadrature points and weights for a $\mathbb{P}_D$ scheme, and all the Lagrange polynomials and their first derivative can be computed automatically.

* The Lagrange basis is nodal, so it allows to have access to the values of a numerical quantity $\phi$ at the quadrature nodes without additional point-wise evaluations, as they are exactly the coefficients of the local polynomial expansion of $\phi$ in the basis.

* Combining this basis to the Gauss-Legendre quadrature nodes is an efficient choice because the numerical integration can be directly expressed as a function of the same coefficients.

## 2.4 Local Finite Elements definitions

The DG framework and a *fortiori* the ADER-DG method rely on the variational formulations of the equations under consideration. In the present paper, this space-time Galerkin approach calls for introducing the following notations, where we recall that $\hat{\omega} = [0, 1]$ and $\hat{c} = \hat{\omega}^3$:

$$ \langle f, g \rangle := \int_{\hat{c}} f(\xi, \eta, \tau)g(\xi, \eta, \tau)d\xi d\eta d\tau, \quad (15a) $$
$$ \{f, g\}_{\xi} := \int_{\hat{\omega}^2} f(\xi, \eta, \tau)g(\xi, \eta, \tau)d\eta d\tau, \qquad (15b) $$

$$ (f, g)_{\eta} := \int_{\hat{\omega}^2} f(\xi, \eta, \tau)g(\xi, \eta, \tau)d\xi d\tau, \qquad (15c) $$

$$ [f, g]_{\tau} := \int_{\hat{\omega}^2} f(\xi, \eta, \tau)g(\xi, \eta, \tau)d\xi d\eta, \qquad (15d) $$

$$ [f, g] := \int_{\hat{\omega}^2} f(\xi, \eta)g(\xi, \eta)d\xi d\eta. \qquad (15e) $$

The above definitions correspond to classical scalar products in Galerkin methods in general. We complete these integral definitions with the introduction of all the matrices involved in the different steps of the ADER-DG procedure for the nonlinear advection-diffusion. The matrices can be grouped as follows:

* The square matrices belonging to the space $M_{S_3}(\mathbb{R})$:

$$ \underline{M}_{k,\ell} := \langle \Theta_k, \Theta_\ell \rangle, \qquad (16a) $$

$$ \underline{M}^{\tau}_{k,\ell} := [\Theta_k, \Theta_\ell]_{\tau}, \qquad (16b) $$

$$ \underline{K}^{v}_{k,\ell} := \langle \partial_v \Theta_k, \Theta_\ell \rangle, \quad v = \xi, \eta, \tau, \qquad (16c) $$

$$ \underline{G}^{v}_{k,\ell} := \partial_v \Theta_\ell(\pi_s, \pi_{s'}, \pi_{s''}), \quad v = \xi, \eta, \quad k := \mathbf{M}(s, s', s''). \qquad (16d) $$

* The rectangle matrices belonging to the space $M_{S_2,S_3}(\mathbb{R})$:

$$ \underline{P}_{k,\ell} := \langle \theta_k, \Theta_\ell \rangle, \qquad (17a) $$

$$ \underline{N}^{\tau}_{k,\ell} := [\theta_k, \Theta_\ell]_{\tau}, \qquad (17b) $$

$$ \underline{C}^{v}_{k,\ell} := \langle \partial_v \theta_k, \Theta_\ell \rangle, \quad v = \xi, \eta. \qquad (17c) $$

* The square space matrices belonging to the space $M_{S_2}(\mathbb{R})$:

$$ \underline{N}_{k,\ell} := [\theta_k, \theta_\ell], \qquad (18a) $$

$$ (\underline{X}_{\xi}^{p})_{k,\ell} := \{\theta_k(p, \eta), \theta_\ell(\eta, \tau)\}_p, \quad p = 0, 1, \qquad (18b) $$

$$ (\underline{X}_{\eta}^{p})_{k,\ell} := (\theta_k(\xi, q), \theta_\ell(\xi, \tau))_q, \quad q = 0, 1, \qquad (18c) $$

$$ (\underline{Y}_{\xi}^{p})_{k,\ell} := \{\partial_{\xi}\theta_{k}(p, \eta), \theta_{\ell}(\eta, \tau)\}_{p}, \quad p = 0, 1, \qquad (18d) $$

$$ (\underline{Y}_{\eta}^{p})_{k,\ell} := (\partial_{\eta}\theta_{k}(\xi, q), \theta_{\ell}(\xi, \tau))_{q}, \quad q = 0, 1, \qquad (18e) $$

$$ (\underline{V}_{\xi}^{p})_{k,\ell} := \Theta_{\ell}(p, \pi_s, \pi_{s'}), \quad p = 0, 1, \quad k = \mathbf{m}(s, s'), \qquad (18f) $$

$$ (\underline{V}_{\eta}^{q})_{k,\ell} := \Theta_{\ell}(\pi_s, q, \pi_{s'}), \quad q = 0, 1, \quad k = \mathbf{m}(s, s'). \qquad (18g) $$

All the matrices introduced are independent of the characteristics of the problem under study such as the mesh, the initial or boundary conditions or the viscosity of the medium. From an implementation perspective, this allows to build them once at the beginning of the program, which is a gain in the global computational cost.
8 ADER-DG schemes for nonlinear advection-diffusion

# 3 High order ADER-DG method

## 3.1 The ADER-DG prediction-correction scheme

In this paper, the DG setting is used for the space discretization of the nonlinear advection-diffusion equation. It consists in representing all the numerical quantities by piecewise space polynomials of arbitrary high degree $D \in \mathbb{N}$. Their expansions in the nodal Lagrange basis presented in subsection 2.3 is then local to each cell of the cartesian mesh. It is well known that a $\mathbb{P}_D$ approximation supposedly provides a numerical scheme of order of convergence $O(\mathbb{P}_D) = D+1$, but the discussion proposed in [11, 29] will be recalled in the present work.

The time integration, however, has to be performed at the same order of accuracy. We propose here an ADER approach combined to the DG scheme in space, leading to the arbitrary high order space-time ADER-DG method. In its most recent form, the ADER-DG scheme [1, 11–14, 17, 18] is composed of two successive stages: prediction and correction. In the prediction step, we determine a local space-time Galerkin predictor, that is to say in each space-time cell, in order to obtain a first estimation of the solution. The correction step then allows to recombine the information between neighbouring cells through numerical fluxes and to impose the boundary conditions. Both steps are carefully described in the following.

### Prediction stage

In the prediction stage, we evolve the solution locally in each cell as if it had no neighbourhood, leading to consistency errors that will be fixed in the correction stage.

Let us consider the version of the equation written in the reference coordinates (7) and defined on the reference space-time cell $\hat{c} = [0, 1]^3$. Multiplying this equation by a space-time test function $\Theta_k$ for $k \in \mathbb{N}(S_3)$ and integrating over $\hat{c}$ yields the variational formulation below:

$$ \langle \Theta_k, \partial_{\tau} \hat{q} \rangle + \alpha^n \langle \Theta_k, \partial_{\xi} \hat{F}(\hat{q}) \rangle + \beta^n \langle \Theta_k, \partial_{\eta} \hat{G}(\hat{q}) \rangle = \Delta t^n \langle \Theta_k, \hat{S} \rangle, \quad (19) $$

where $\hat{q}(\xi, \eta, \tau)$ is the local Galerkin predictor. Then, we remark that the integration by parts with respect to the reference time variable $\tau$ on the first term of (19) gives:

$$ \langle \Theta_k, \partial_{\tau} \hat{q} \rangle = -\langle \partial_{\tau} \Theta_k, \hat{q} \rangle + [\Theta_k, \hat{q}]_1 - [\Theta_k, \hat{q}]_0. \quad (20) $$

It has to be noted that, on the right hand side of (20), the quantity $\hat{q}(\xi, \eta, 1)$ in the second term is unknown as it corresponds to the value of the numerical solution $q$ at the next time $t^{n+1}$. Yet, we need to enforce the condition $\hat{q}(\xi, \eta, 0) = \hat{u}_0(\xi, \eta)$ corresponding to the numerical solution $u_{i,j}^n$ on cell $\Omega_{i,j}$.
Variational formulation (19) becomes:

$$ \begin{aligned} & [\Theta_k, \hat{q}]_1 - \langle \partial_\tau \Theta_k, \hat{q} \rangle + \alpha^n \langle \Theta_k, \partial_\xi \hat{F}(\hat{q}) \rangle + \beta^n \langle \Theta_k, \partial_\eta \hat{G}(\hat{q}) \rangle \\ & = [\Theta_k, \hat{u}_0]_0 + \Delta t^n \langle \Theta_k, \hat{S} \rangle. \end{aligned} \quad (21) $$

In the 2D space and 3D space-time Lagrange basis respectively, the initial condition $\hat{u}_0$ and the unknowns $\hat{\phi} = \hat{q}, \hat{F}, \hat{G}, \hat{S}$ have the following expansions:

$$ \hat{u}_0(\xi, \eta) := \sum_{\ell=1}^{S_2} (\hat{u}_0)_\ell \theta_\ell(\xi, \eta), \quad \hat{\phi}(\xi, \eta, \tau) := \sum_{\ell=1}^{S_3} \hat{\phi}_\ell \Theta_\ell(\xi, \eta, \tau), \quad (22) $$

and the partial derivative of $\hat{\phi}$ is transferred to the basis functions:

$$ \partial_v \hat{\phi}(\xi, \eta, \tau) := \sum_{\ell=1}^{S_3} \hat{\phi}_\ell \partial_v \Theta_\ell(\xi, \eta, \tau), \quad v = \xi, \eta, \tau. \quad (23) $$

Consequently, and because the Lagrange polynomials form a nodal basis, we can define for $m = 1, 2, 3$ the vector $\mathbf{\Phi} := [\phi_1, \dots, \phi_{S_m}]^T \in \mathbb{R}^{S_m}$ containing the Degrees Of Freedom (DOF) of the $m$D scalar functions $\phi$ with respect to the $m$D Lagrange polynomial basis.

For the sake of clarity, let us continue the computations for each component $z = x, y$ of the problem, that is to say by treating the vector equation (21) of size 2 as two scalar equations. By substituting expansions (22) in (21), and using property (23) as well as matrices presented in subsection 2.4, we show that the prediction stage can be written under the compact matrix form:

$$ \underline{A}\hat{q}_z + \underline{B}^\xi \hat{F}_z + \underline{B}^\eta \hat{G}_z = \hat{c}_z, \quad z = x, y, \quad (24) $$

where we have introduced the intermediate matrices and vectors:

$$ \underline{A} := \underline{M}^1 - \underline{K}^\tau, \quad \underline{B}^\xi := \alpha^n (\underline{K}^\xi)^T, \quad \underline{B}^\eta := \beta^n (\underline{K}^\eta)^T, \quad (25a) $$

$$ \hat{c}_z := (\underline{N}^0)^T \hat{u}_z^0 + \Delta t^n \underline{M} \hat{S}_z, \quad z = x, y. \quad (25b) $$

In the particular case of a linear conservation law of the form (21), that is to say when the physical fluxes $F$ and $G$ are linear with respect to the unknown $q$, it is possible to find matrices $\underline{L}^\xi$ and $\underline{L}^\eta$ such that $\hat{F}_z = \underline{L}^\xi \hat{q}_z$ and $\hat{G}_z = \underline{L}^\eta \hat{q}_z$, for $z = x, y$. Such matrices may depend on the problem under consideration, through the physical parameters involved, or the mesh, for instance. These relations lead to rewrite the previous linear system (24) as $\underline{A}'\hat{q}_z = \hat{c}_z$ that can be inverted directly, where $\underline{A}' := \underline{A} + \underline{B}^\xi \underline{L}^\xi + \underline{B}^\eta \underline{L}^\eta$.

However, in the case of the present article, the PDE studied is nonlinear due to the advective component of the total flux. The two equations of the prediction stage (for $z = x, y$) are thus nonlinear with respect to $\hat{q}_z$ and coupled.
10 ADER-DG schemes for nonlinear advection-diffusion

We propose to use a Picard fixed point algorithm to solve this nonlinear system for $\hat{\mathbf{q}}_x$ and $\hat{\mathbf{q}}_y$. In the following, we introduce the notation $\hat{\mathbf{\Phi}} := [\hat{\mathbf{\Phi}}_x, \hat{\mathbf{\Phi}}_y]$. Then, starting from an initial guess $\hat{\mathbf{q}}^0 := [\hat{\mathbf{q}}_x^0, \hat{\mathbf{q}}_y^0]$, we execute the following instructions until convergence:

1. Compute $\hat{\mathbf{F}}_z^m := \mathbf{F}_z(\hat{\mathbf{q}}^m)$ and $\hat{\mathbf{G}}_z^m := \mathbf{G}_z(\hat{\mathbf{q}}^m)$ for $z = x, y$,

2. Compute $\hat{\mathbf{q}}^{m+1}$ such that $\underline{\mathbf{A}}\hat{\mathbf{q}}_z^{m+1} = \hat{\mathbf{c}}_z - \underline{\mathbf{B}}^\xi \hat{\mathbf{F}}_z^m - \underline{\mathbf{B}}^\eta \hat{\mathbf{G}}_z^m$ for $z = x, y$,

where the initialization can be $\hat{\mathbf{q}}^0 = \mathbf{q}^{n-1}$, that is to say the predictor finally obtained at the end of the Picard loop of the previous iteration. The first initiatization at the physical time $t = 0$ may be done by using the initial condition extruded in the time direction on $[t^0, t^1]$. In practice, a few iterations are enough to achieve convergence in the case of the nonlinear advection-diffusion equation. In our implementation, we chose to perform the Picard instructions until at least one of the two following criteria is satisfied: executing five iterations or reaching a relative error between two consecutive iterations smaller than $10^{-8}$. Yet, more iterations may be needed for more complex nonlinear conservation laws such as the Euler, compressible Navier-Stokes or shallow water problems for example.

Let us close the description of the prediction step of the ADER-DG scheme with additional remarks. Indeed, as the nonlinearity of the total fluxes comes from the advective components, it is also possible to express the linear diffusive fluxes as functions of the predictor $\hat{\mathbf{q}}$, that is to say $(\hat{\mathbf{F}}_d)_z := \underline{\mathbf{L}}^\xi \hat{\mathbf{q}}_z$ and $(\hat{\mathbf{G}}_d)_z := \underline{\mathbf{L}}^\eta \hat{\mathbf{q}}_z$ for $z = x, y$. Here, if we denote by $\underline{\nu}$ the diagonal matrix containing the values of the viscosity function $\underline{\nu}$ at the space-time quadrature nodes, we have $\underline{\mathbf{L}}^\xi := -\underline{\nu}\Delta x^{-1}\underline{\mathbf{G}}^\xi$ and $\underline{\mathbf{L}}^\eta := -\underline{\nu}\Delta y^{-1}\underline{\mathbf{G}}^\eta$, which can be deduced from definitions (8) and expansions (22). As a consequence, the nonlinear system (24) can be rewritten as:

$$ \underline{\mathbf{A}}'\hat{\mathbf{q}}_z + \underline{\mathbf{B}}^\xi(\hat{\mathbf{F}}_a)_z + \underline{\mathbf{B}}^\eta(\hat{\mathbf{G}}_a)_z = \hat{\mathbf{c}}_z, \quad z = x, y, \qquad (26) $$

with the matrix $\underline{\mathbf{A}}'$ already defined, and we can apply the Picard fixed point algorithm previously presented to this new formulation by updating the advective fluxes only. This approach provides the same numerical solution.

Solving the local linear system (24) or (26) of size $S_3$, that is to say of size 8, 27, 64 or 125 for a $\mathbb{P}_D$ approximation for $D \in \{1, \dots, 4\}$, has a reasonable computational cost compared to a global system. However, because of the presence of the coefficients $\alpha^n, \beta^n$ and $\Delta t^n$ that depend on the iteration, a majority of the matrices involved need to be recomputed at each time step of the ADER-DG procedure. Finally, regarding the implementation aspect, the prediction stage can be run in parallel without any difficulty as each cell is independent of its neighbours.

## Correction stage

At the end of the prediction stage, a continuous space-time polynomial $\hat{\mathbf{q}}(\xi, \eta, \tau)$ of arbitrary high degree $D$ has been determined locally in each cell
independently of the others. The correction stage then allows to recombine the information between neighbouring cells through numerical fluxes and to impose the boundary conditions.

We consider, again, the version (7) of the nonlinear advection-diffusion equation. The aim is now to numerically update the unknown $\hat{u}$ with the knowledge of the predictors $\hat{q}$ inside the cells as well as their limits on their faces. The underlying DG scheme consists here in multiplying the PDE by a space test function $\theta_k$ for $k \in \mathbb{N}(S_2)$ and integrating over $\hat{c}$, so that the corresponding DG variational formulation writes:

$$ \langle \theta_k, \partial_{\tau} \hat{u} \rangle + \alpha^n \langle \theta_k, \partial_{\xi} \hat{F}(\hat{q}) \rangle + \beta^n \langle \theta_k, \partial_{\eta} \hat{G}(\hat{q}) \rangle = \Delta t^n \langle \theta_k, \hat{S} \rangle, \quad (27) $$

where $\hat{u}$ has been replaced by $\hat{q}$ in the fluxes terms.

As the test functions do not depend on the reference time variable $\tau$, the integration in time on the first term can be performed exactly. The two fluxes terms, however, are integrated by parts, giving rise to integral terms on the faces where the predictor $\hat{q}$ is not uniquely defined. Indeed, at every face shared by two neighbouring cells, there are two limit values from both sides and there is no reason a priori for these values to be equal. All the quantities $\phi$ related to the predictor in the interface terms, that is to say $\hat{\phi} = \hat{q}, \hat{F}(\hat{q}), \hat{G}(\hat{q})$, are thus replaced by numerical fluxes denoted by $\tilde{\phi}$ respectively. The correction step finally reads:

$$ \begin{aligned} [\theta_k, \hat{u}(\cdot, \cdot, 1)] = [\theta_k, \hat{u}(\cdot, \cdot, 0)] &+ \Delta t^n \langle \theta_k, \hat{S} \rangle \\ &+ \alpha^n \langle \partial_{\xi} \theta_k, \hat{F}(\hat{q}) \rangle - \alpha^n (\{\theta_k, \tilde{F}(\tilde{q})\}_1 - \{\theta_k, \tilde{F}(\tilde{q})\}_0) \\ &+ \beta^n \langle \partial_{\eta} \theta_k, \hat{G}(\hat{q}) \rangle - \beta^n ((\theta_k, \tilde{G}(\tilde{q}))_1 - (\theta_k, \tilde{G}(\tilde{q}))_0), \end{aligned} \quad (28) $$

which, for the sake of brevity, can be shortened as follows:

$$ [\theta_k, \hat{u}(\cdot, \cdot, 1)] = [\theta_k, \hat{u}(\cdot, \cdot, 0)] + \psi_{\text{rhs}}(k, \Delta t^n, \hat{q}, \tilde{q}), \quad (29) $$

with $\psi_{\text{rhs}}$ the quantity regrouping the volume terms and the interface terms of the right hand side (rhs) of equation (28), depending on the k-th test function $\theta_k$ and the time step $\Delta t^n$.

In the DG approximation of the PDE, we expand the numerical solution $\hat{u}$ with respect to the 2D Lagrange space basis. The three space-time volume terms are written in the 3D space-time Lagrange basis, in accordance with the formula (22) for $\hat{\phi} = \hat{S}, \hat{F}, \hat{G}$. In each cell of the cartesian grid, the following expansion holds:

$$ \hat{u}(\xi, \eta, \tau) := \sum_{\ell=1}^{S_2} \hat{u}_{\ell}(\tau) \theta_{\ell}(\xi, \eta). \quad (30) $$

The numerical solution is thus a piecewise polynomial function on the whole domain $\Omega$ at a given iteration $t^n$. The complete definition of the DG scheme
12 ADER-DG schemes for nonlinear advection-diffusion

now depends on the choice of the numerical fluxes $\tilde{\phi}$ at cell interfaces, presented in subsection 3.2. The two next subsections 3.3 and 3.4 propose two approaches respectively: the SIPG penalty and the Cattaneo relaxation methods.

## 3.2 Numerical fluxes at cells interfaces

In a discontinuous finite elements approach like the DG scheme, the numerical solution may be discontinous across an interface separating two space-time cells. In order to give a unique and consistent expression to the solution at a cell boundary, we need to define a conservative numerical flux between the cell and its direct neighbours.

Let us consider the interface $x = x_e$ separating two cells $c_L$ and $c_R$, and the interface $y = y_e$ separating two cells $c_B$ and $c_T$. For any vector function $\phi$, we define its limits from both sides of the given interface by:

$$ \phi_{x_e}^{\pm} := \lim_{\delta \to 0^+} \phi(x_e \pm \delta, y, t), \quad \phi_{y_e}^{\pm} := \lim_{\delta \to 0^+} \phi(x, y_e \pm \delta, t), \quad (31) $$

where the two other variables have been dropped from the expressions. For instance, we have $\phi_{x_e}^{\pm} = \phi_{x_e}^{\pm}(y, t)$, and similarly for $\phi_{y_e}^{\pm}$. The average $\langle \cdot \rangle$ and the jump $\llbracket \cdot \rrbracket$ of $\phi$ can now be expressed at any interface $x = x_e$ or $y = y_e$ common to two neighbouring cells. We have respectively:

$$ \langle \phi \rangle_{x_e} := \frac{1}{2} [\phi_{x_e}^{-} + \phi_{x_e}^{+}], \quad \llbracket \phi \rrbracket_{x_e} := \phi_{x_e}^{-} - \phi_{x_e}^{+}, \quad (32a) $$

$$ \langle \phi \rangle_{y_e} := \frac{1}{2} [\phi_{y_e}^{-} + \phi_{y_e}^{+}], \quad \llbracket \phi \rrbracket_{y_e} := \phi_{y_e}^{-} - \phi_{y_e}^{+}. \quad (32b) $$

In the present work, we choose a Local Lax-Friedrichs (LLF) numerical flux [1, 11], also known as the Rusanov numerical flux. It relies on the averages and jumps of the solution, as well as the stabilization coefficients for the advective components: $\sigma_a^z := |u_z|$ for $z = x, y$, and the unique stabilization coefficient $\sigma_d$ for the diffusive components. The latter will differ depending on whether the penalty or relaxation approach is chosen, and it will be precised in the dedicated subsections 3.3 and 3.4 respectively.

### Interior interfaces

An interior interface is a space-time face which is shared by two adjacent cells. The LLF numerical fluxes $\tilde{H}$ for the fluxes $H = F, G$, and the numerical flux $\tilde{q}$ for the unknown itself $q$, are expressed for an interior interface $x = x_e$ or $y = y_e$ respectively, for each component $c = a, d$, as:

$$ \tilde{F}_c(\tilde{q}_{x_e}) := \langle F_c(q) \rangle_{x_e} + \frac{1}{2} \sigma_c^x(x_e) \llbracket q \rrbracket_{x_e}, \quad \tilde{q}_{x_e} := \frac{1}{2} \llbracket q \rrbracket_{x_e} \vec{n}_{x_e}, \quad (33a) $$

$$ \tilde{G}_c(\tilde{q}_{y_e}) := \langle G_c(q) \rangle_{y_e} + \frac{1}{2} \sigma_c^y(y_e) \llbracket q \rrbracket_{y_e}, \quad \tilde{q}_{y_e} := \frac{1}{2} \llbracket q \rrbracket_{y_e} \vec{n}_{y_e}, \quad (33b) $$
where $\vec{n}_{x_e}$ is equal to +1 for the left cell $c_L$ and -1 for the right cell $c_R$, and the same applies to $\vec{n}_{y_e}$. In order to compact the expressions, we have introduced the notation $\sigma_c^x(x_e) := \sigma_c^x(x_e, y, t)$ for example, and the same remark holds for $\sigma_c^y(y_e)$, for $c = a, d$.

### Dirichlet interfaces

A Dirichlet interface is a space-time face which is located at a part of the boundary $\partial\Omega$ of the computational domain where a Dirichlet boundary condition $q = h_D$ is imposed. In the numerical fluxes, the jump terms have the same form for the advective and diffusive parts, but the average of the physical flux is not the same. We have:

$$ \tilde{F}_a(\tilde{q}_{x_e}) := \frac{1}{2} [F_a(q_{x_e}^{\pm}) + F_a(h_D)] + \frac{1}{2} \sigma_a^x(x_e) [q_{x_e}^{\pm} - h_D] \vec{n}_{x_e}, \quad (34a) $$

$$ \tilde{F}_d(\tilde{q}_{x_e}) := F_d(q_{x_e}^{\pm}) + \frac{1}{2} \sigma_d^x(x_e) [q_{x_e}^{\pm} - h_D] \vec{n}_{x_e}, \quad (34b) $$

$$ \tilde{q}_{x_e} := [q_{x_e}^{\pm} - h_D], \quad (34c) $$

$$ \tilde{G}_a(\tilde{q}_{y_e}) := \frac{1}{2} [G_a(q_{y_e}^{\pm}) + G_a(h_D)] + \frac{1}{2} \sigma_a^y(y_e) [q_{y_e}^{\pm} - h_D] \vec{n}_{y_e}, \quad (34d) $$

$$ \tilde{G}_d(\tilde{q}_{y_e}) := G_d(q_{y_e}^{\pm}) + \frac{1}{2} \sigma_d^y(y_e) [q_{y_e}^{\pm} - h_D] \vec{n}_{y_e}, \quad (34e) $$

$$ \tilde{q}_{y_e} := [q_{y_e}^{\pm} - h_D], \quad (34f) $$

where the same previous remarks and notations hold. For example, in the expression of the first numerical flux $\tilde{F}_a(\tilde{q}_{x_e})$, we have $q_{x_e}^{\pm} = q_{x_e}^{+}$ or $q_{x_e}^{-}$, and $\vec{n}_{x_e} = -1$ or +1, respectively for a Dirichlet interface located on the left or right boundary of the domain $\Omega$.

### Neumann interfaces

A Neumann interface is a space-time face which is located at a part of the boundary $\partial\Omega$ of the computational domain where a Neumann boundary condition $H_d(q) = h_N$ is imposed, for $H = F, G$. The expressions of the numerical fluxes are immediate:

$$ \tilde{F}_d(\tilde{q}_{x_e}) := h_N, \quad \tilde{q}_{x_e} := 0, \quad \tilde{G}_d(\tilde{q}_{y_e}) := h_N, \quad \tilde{q}_{y_e} := 0. \quad (35) $$

At this stage, the DG correction step (28) or equivalently (29) is not completely defined, as we did not precise yet the value of the diffusive stabilization coefficient $\sigma_d$. The two next subsections 3.3 and 3.4 propose two approaches respectively: the SIPG penalty and the Cattaneo relaxation methods, each with its own specific features.

## 3.3 A penalty approach: the SIPG method

The SIPG method is one of the three versions of interior penalty (IP) techniques presented in [11, 29]. This penalty approach consists in modifying the
14 ADER-DG schemes for nonlinear advection-diffusion

DG correction step, for instance in its compact form (29), by adding extra
terms as follows:

$$
\begin{equation}
\begin{aligned}
[\theta_k, \hat{u}(\cdot, \cdot, 1)] = [\theta_k, \hat{u}(\cdot, \cdot, 0)] &+ \psi_{\text{rhs}}(k, \Delta t^n, \hat{q}, \tilde{q}) \\
& + \varepsilon \alpha^n (\{\hat{F}_d(\theta_k), \tilde{q}\}_1 - \{\hat{F}_d(\theta_k), \tilde{q}\}_0) \\
& + \varepsilon \beta^n ((\hat{G}_d(\theta_k), \tilde{q})_1 - (\hat{G}_d(\theta_k), \tilde{q})_0),
\end{aligned}
\tag{36}
\end{equation}
$$

where $\varepsilon \in \{-1, 0, +1\}$ is a discrete parameter leading or not to a symmetry or nonsymmetry with respect to the diffusive interface terms. Indeed, by breaking the total physical fluxes into their advective and diffusive contributions, the right hand side of the equation contains for example the quantity that can be rearranged as:

$$
\pm \alpha^n \left( \{\theta_k, \tilde{F}_d(\tilde{q})\}_p - \varepsilon \{\hat{F}_d(\theta_k), \tilde{q}\}_p \right), \quad p=0,1, \quad (37)
$$

and the same comment holds for the fluxes **G** in the *y* direction, affected by
the coefficient ±*β*ⁿ, and their (non)symmetric couterparts. More precisely:

• the SIPG (symmetric) version corresponds to the case ε = −1,

• the IIPG (incomplete) version corresponds to the case ε = 0,

• the NIPG (nonsymmetric) version corresponds to the case ε = +1.

Of course, the case ε = 0 is, from a formal perspective, strictly equivalent
to the inital problem (29). All three IPDG approaches have in common the
fact that the stabilization coefficients σ<sup>z</sup><sub>d</sub> for z = x, y, in the numerical fluxes
defined in subsection 3.2, have the form:

$$
\sigma_d(x, y, t) := \frac{\gamma^{\text{IP}}(D)}{h} \nu(x, y, t), \quad \gamma^{\text{IP}}(D) > 0, \qquad (38)
$$

with IP = S, I, N for the SIPG, IIPG or NIPG versions respectively. In prac-
tice, we take γ<sup>IP</sup>(D) = 10D<sup>2</sup> for a P<sub>D</sub> approximation, but a broader reflection
can be found in [26]. In the present article, the symmetric version (SIPG)
will be generally preferred as [11] has shown it is faster than the other two.
This SIPG scheme is finally characterized by taking ε = −1 and the diffusive
stabilization coefficient (38).

3.4 A relaxation approach: the Cattaneo method

A relaxation approach consists in considering a modified problem which solu-
tion converges to the one of the original problem in the limit of a vanishing
parameter. In the Cattaneo relaxation scheme, the nonlinear advection-
diffusion equation recalled below:

$$
\partial_t u + \partial_x (\mathbf{F}_a + \mathbf{F}_d) + \partial_y (\mathbf{G}_a + \mathbf{G}_d) = \mathbf{S}, \quad (39)
$$
can be approached by the modified problem:

$$ \left\{ \begin{array}{l} \partial_t \mathbf{u} + \partial_x (\mathbf{F}_a - \nu \mathbf{a}) + \partial_y (\mathbf{G}_a - \nu \mathbf{b}) = \mathbf{S}, \\ \partial_t [\mathbf{a}, \mathbf{b}] = \frac{1}{\delta} (\nabla \mathbf{u} - [\mathbf{a}, \mathbf{b}]), \end{array} \right. \qquad (40) $$

where we recall that $h = \max(\Delta x, \Delta y)$ and we have introduced the small parameter $0 < \delta \ll 1$ and the vector functions $\mathbf{a}, \mathbf{b} : \Omega \times [0, T_f] \to \mathbb{R}^2$. We thus have $[\mathbf{a}, \mathbf{b}] \to \nabla \mathbf{u}$ in the limit of $\delta \to 0$. If we define the following vectors and matrices:

$$ \vec{\mathbf{u}} := \begin{bmatrix} \mathbf{u} \\ \mathbf{a} \\ \mathbf{b} \end{bmatrix}, \quad \vec{\mathbf{S}} := \begin{bmatrix} \mathbf{S} \\ -\mathbf{a}/\delta \\ -\mathbf{b}/\delta \end{bmatrix}, \qquad (41a) $$

$$ \underline{D}^x := \begin{bmatrix} u_x & -\nu & 0 \\ -1/\delta & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad \underline{D}^y := \begin{bmatrix} u_y & 0 & -\nu \\ 0 & 0 & 0 \\ -1/\delta & 0 & 0 \end{bmatrix}, \qquad (41b) $$

then this manipulation leads to the following hyperbolic system:

$$ \partial_t \vec{\mathbf{u}} + \partial_x (\underline{D}^x \vec{\mathbf{u}}) + \partial_y (\underline{D}^y \vec{\mathbf{u}}) = \vec{\mathbf{S}}. \qquad (42) $$

The error due to the DG scheme applied to the first equation of the relaxed system (40) is supposed to be $O(h^{D+1})$ for a $\mathbb{P}_D$ approximation, while the second equation introduces an approximation error of order of magnitude $O(\delta)$. If we want this relaxation error to be negligible compared to the DG approximation error, a natural choice is to take $\delta = O(h^{D+1})$. In the present work, we used the formula proposed in [1]:

$$ \delta = \left( \frac{2^{D+1} - \sqrt{2}}{\sqrt{2} - 1} \right) h^{D+1}. \qquad (43) $$

Our aim now is not to solve this relaxed system (42), which size is three times the one of the initial system (40), even though it may provide the expected results. However, we can identify the nonzero eigenvalues $\lambda^z$ of matrices $\underline{D}^z$ for $z = x, y$. A few computations yield:

$$ \lambda^z := \frac{1}{2} \left( u_z \pm \sqrt{u_z^2 + \frac{4\nu}{\delta}} \right). \qquad (44) $$

An upper bound of their modulus can be determined easily by using the triangular inequality and the concavity of the square root function:

$$ |\lambda^z| \le \frac{1}{2} \left( |u_z| + \sqrt{u_z^2 + \frac{4\nu}{\delta}} \right) \le |u_z| + \sqrt{\frac{\nu}{\delta}}. \qquad (45) $$
16 ADER-DG schemes for nonlinear advection-diffusion

On the one hand, in the particular case of a pure advection equation by taking $\nu \equiv 0$, which is not the object of this article, we recover the advective stabilization coefficients $\sigma_a^z = |u_z|$ in the last bound obtained in (45). On the other hand, for a pure diffusion problem, that is to say by taking $F_a = G_a = 0$, we identify the stabilization coefficients for the diffusive components. Their expressions read:

$$ \sigma_d(x, y, t) := \sqrt{\frac{\nu(x, y, t)}{\delta}} = \sqrt{\frac{\gamma^C(D)}{h^{D+1}}\nu(x, y, t)}, \quad \gamma^C(D) > 0, \quad (46) $$

with $\gamma^C(D) = (\sqrt{2} - 1)/(2^{D+1} - \sqrt{2})$, due to the choice (43). The Cattaneo relaxation scheme can finally be fully described by taking $\varepsilon = 0$ in equation (36), meaning that no additional terms are considered, as well as the stabilization coefficient (46) for the diffusive components. Again, we draw the reader's attention to the fact that the relaxed system is not solved numerically: it is a numerical analysis tool allowing to identify a stabilization coefficient for the diffusion term of the original equation (39).

## 3.5 The CFL condition

As every explicit numerical scheme does, our ADER-DG approach needs to be completed with a well defined CFL condition to ensure the stability of the scheme. For a given characteristic length $h := \max(\Delta x, \Delta y)$ of the mesh, a CFL condition at the iteration $t^n \to t^{n+1}$ always has the form:

$$ \Delta t^n := \kappa \frac{h}{\max|\sigma|^n}, \quad \kappa > 0, \quad (47) $$

where the maximum value $\max|\phi|$ for a given scalar function $\phi$ on the whole space domain $\Omega$ at time $t = t^n$ is defined by:

$$ \max|\phi|^n := \max_{i,j} \max_{\ell=\mathbf{m}(s,s')} \max_{z=x,y} |\phi_{i,j}^z(\mathcal{X}_i(\pi_s), \mathcal{Y}_j(\pi_{s'}), t^n)|. \quad (48) $$

In other words, $\max|\phi|^n$ is the maximum over all the cells $\Omega_{i,j}$ of the quantities $\phi(x, y, t^n)$ for all 2D points $(x, y)$ being the $S_2$ DOFs of the solution on the cell, for the two components $z = x, y$. In the case of the present article, which is dedicated to the nonlinear advection-diffusion with an a priori variable viscosity $\nu$, the advective maximum wave speed $(\sigma_a)_*$ depends on the solution itself, contrary to the diffusive maximum wave speed $(\sigma_d)_*$ obtained by formulas (38) or (46) if the SIPG or Cattaneo approach is chosen. Indeed, using the compact notations, the advective and the SIPG or Cattaneo diffusive maximum wave speeds write respectively:

$$ (\sigma_a)_* := \max|u|, \quad (49) $$

$$ (\sigma_d^S)_* := \frac{\gamma^S(D)}{h} \max|\nu|, \quad \gamma^S(D) := 10D^2, \quad (50) $$
$$ (\sigma_d^C)_* := \sqrt{\frac{\gamma^C(D)}{h^{D+1}}} \max|\nu|, \quad \gamma^C(D) := \frac{\sqrt{2}-1}{2^{D+1}-\sqrt{2}}. \qquad (51) $$

Consequently, the corresponding time steps for advection and diffusion read:

$$ \Delta t_a := \kappa_a(D) \frac{h}{(\sigma_a)_*}, \quad \Delta t_d^M := \kappa_d^M(D) \frac{h}{(\sigma_d^M)_*}, \quad M = S, C, \qquad (52) $$

where $0 < \kappa_a(D) < 1$ and $0 < \kappa_d^M(D) < 1$, for $M = S, C$, are coefficients that can be adjusted numerically. The values we have determined experimentally are shown in table 1. In particular, we notice that the time steps for a DG scheme depend on the approximation degree $D$, which is not the case for finite volumes methods. The reader can refer to [7, 17, 30, 31] for a broader view on this aspect.

<table><caption>**Table 1** Coefficients $\kappa_a(D)$ and $\kappa_d^M(D)$ for $M = S, C$ and $D \in \{0, \dots, 4\}$.</caption><thead><tr><th>D</th><th>κ<sub>a</sub>(D)</th><th>κ<sub>d</sub><sup>S</sup>(D)</th><th>κ<sub>d</sub><sup>C</sup>(D)</th></tr></thead><tbody><tr><td>0</td><td>0.50</td><td>1.0</td><td>1.0</td></tr><tr><td>1</td><td>0.16</td><td>0.20</td><td>0.16</td></tr><tr><td>2</td><td>0.089</td><td>0.095</td><td>0.083</td></tr><tr><td>3</td><td>0.056</td><td>0.056</td><td>0.050</td></tr><tr><td>4</td><td>0.038</td><td>0.037</td><td>0.034</td></tr></tbody></table>

The determination of a suitable CFL condition that will ensure numerically the stability of the scheme for the complete advection-diffusion problem is not immediate, and there is not necessarily a unique approach. In order to be able to pursue the comparison between the SIPG and Cattaneo approaches, we prefer to use the following relation to determine the time step:

$$ \frac{1}{\Delta t^M} := \frac{1}{\Delta t_a} + \frac{1}{\Delta t_d^M}, \qquad (53) $$

which is the classical approach (see [1, 11, 14, 15, 31] for instance), although other options could be considered, such as the one proposed in [32]. It is clear that the advective time step is common to both the SIPG penalty and Cattaneo relaxation methods, while their diffusive time steps differ, leading a *priori* to different computational costs. This aspect will be discussed in section 4 dedicated to numerical results.

## 3.6 The ADER-DG complete algorithm

This subsection is dedicated to writing all the steps of the ADER-DG scheme applied to the nonlinear advection-diffusion equation (2), whether the SIPG penalty or the Cattaneo relaxation scheme is chosen. Indeed, this will emphasize how convenient this ADER-DG approach is from an implementation
18 ADER-DG schemes for nonlinear advection-diffusion

perspective, by using all the matrices (16)-(17)-(18) that can be precomputed once. The complete algorithm is described below:

* Make a loop over the cells $\Omega_{i,j}$ to initialize the unknowns $(\mathbf{u}_z)_{i,j}^0 \in \mathbb{R}^{S_2}$ for $z = x, y$ with the initial condition $\mathbf{u}_0(x, y)$ evaluated at the corresponding DOFs, as well as the predictors $(\mathbf{q}_z)_{i,j}^0 \in \mathbb{R}^{S_3}$ in a similar way:

$$((\mathbf{u}_z)_{i,j}^0)_\ell := (u_z)_0 (\mathcal{X}_i(\pi_k), \mathcal{Y}_j(\pi_{k'})), \quad \ell := \mathbf{m}(k, k'), \quad (54)$$

$$((\mathbf{q}_z)_{i,j}^0)_\ell := (u_z)_0 (\mathcal{X}_i(\pi_k), \mathcal{Y}_j(\pi_{k'})), \quad \ell := \mathbf{M}(k, k', k''). \quad (55)$$

* While the final time $T_f$ is not reached, that is to say $t^n < T_f$ at iteration $n \in \{0, \dots, N_t - 1\}$, follow the next steps:

- Compute the time step $\Delta t^n := t^{n+1} - t^n$ with formula (53) affected by a security coefficient 0.9 to add a safety margin for stability:

$$\Delta t^n := 0.9 \left( \frac{1}{\Delta t_a^n} + \frac{1}{\Delta t_d^n} \right)^{-1}. \quad (56)$$

- Make a loop over the cells $\Omega_{i,j}$ to compute for $z = x, y$ the source terms $(\mathbf{S}_z)_{i,j}^n \in \mathbb{R}^{S_3}$ by evaluating the source term function $\mathbf{S}$ in the space-time cell $c_{i,j}^n$ at the corresponding DOFs:

$$((\mathbf{S}_z)_{i,j}^n)_\ell := S_z (\mathcal{X}_i(\pi_k), \mathcal{Y}_j(\pi_{k'}), \mathcal{T}_n(\pi_{k''})), \quad \ell := \mathbf{M}(k, k', k''). \quad (57)$$

- Make a loop over the cells $\Omega_{i,j}$ to compute the predictors $(\mathbf{q}_z)_{i,j}^n \in \mathbb{R}^3$ in the space-time cells $c_{i,j}^n$ for $z = x, y$ by solving the prediction system (24) or (26) with a Picard iteration algorithm (see subsection 3.1 for details):

$$\underline{A}(\mathbf{q}_z)_{i,j}^n + \underline{B}^\xi(\mathbf{F}_z)_{i,j}^n + \underline{B}^\eta(\mathbf{G}_z)_{i,j}^n = (\mathbf{c}_z)_{i,j}^n. \quad (58)$$

- Make a loop over the faces $z = z_e$ between neighbouring cells or at the boundary $\partial\Omega$ to compute the numerical fluxes $\tilde{\mathbf{H}}_e^n \in \mathbb{R}^{S_2}$ as presented in subsection 3.2, for $z = x, y$ and $\mathbf{H} = \mathbf{F}, \mathbf{G}$ respectively. We give here the example of the numerical fluxes $(\tilde{\mathbf{F}}_z)_{x_i, j}^n$ and $(\tilde{\mathbf{q}}_z)_{x_i, j}^n$ at the interior interface $x = x_i$ shared by cells $c_{i,j}^n$ and $c_{i+1, j}^n$. We note that the space-time fluxes $\mathbf{F}_z^n \in \mathbb{R}^{S_3}$ can be deduced from the predictors $\mathbf{q}_z^n$, as explained in subsection 3.1. We have:

$$(\tilde{\mathbf{F}}_z)^n_{x_i, j} := \frac{1}{2} (\underline{V}_\xi^1(\mathbf{F}_z)^n_{x_i, j} + \underline{V}_\xi^0(\mathbf{F}_z)^n_{i+1, j}) + \frac{1}{2} \sigma^z_{x_i, j} (\underline{V}_\xi^1(\mathbf{q}_z)^n_{x_i, j} - \underline{V}_\xi^0(\mathbf{q}_z)^n_{i+1, j}), \quad (59)$$

$$(\tilde{\mathbf{q}}_z)^n_{x_i, j} := \frac{1}{2} (\underline{V}_\xi^1(\mathbf{q}_z)^n_{x_i, j} + \underline{V}_\xi^0(\mathbf{q}_z)^n_{i+1, j}), \quad (60)$$
where the diagonal matrix of stabilization coefficients is defined by:

$$ (\underline{\underline{\sigma}}_{x_i, j}^z)_{\ell, \ell} := \sigma^z(x_i, \mathcal{Y}_j(\pi_{k'}), \mathcal{T}_n(\pi_{k''})), \quad \ell := \mathbf{m}(k', k''). \quad (61) $$

- Make a loop over the cells $\Omega_{i,j}$ to determine the updated unknowns $(\mathbf{u}_z)_{i,j}^{n+1} \in \mathbb{R}^{S_2}$, for $z = x, y$, by using the DG correction system (36) or (29) if the SIPG or Cattaneo approach is chosen respectively. The two approaches can be cast under the unique form (36) and, after a few computations, this stage writes:

$$ \begin{aligned} \underline{N}(\mathbf{u}_z)_{i,j}^{n+1} = \underline{N}(\mathbf{u}_z)_{i,j}^n &+ \alpha^n \underline{C}^\xi (\mathbf{F}_z)_{i,j}^n + \beta^n \underline{C}^n (\mathbf{G}_z)_{i,j}^n + \Delta t^n \underline{P}(\mathbf{S}_z)_{i,j}^n \\ & - \alpha^n \left( \underline{X}_\xi^1 (\tilde{\mathbf{F}}_z)_{x_i, j}^n - \underline{X}_\xi^0 (\tilde{\mathbf{F}}_z)_{x_{i-1}, j}^n \right) \\ & - \beta^n \left( \underline{X}_\eta^1 (\tilde{\mathbf{G}}_z)_{i, y_j}^n - \underline{X}_\eta^0 (\tilde{\mathbf{G}}_z)_{i, y_{j+1}}^n \right) \\ & - \varepsilon \frac{\alpha^n}{\Delta x} \left( \underline{Y}_\xi^1 \nu_{x_i, j} (\tilde{\mathbf{q}}_z)_{x_i, j}^n - \underline{Y}_\xi^0 \nu_{x_{i-1}, j} (\tilde{\mathbf{q}}_z)_{x_{i-1}, j}^n \right) \\ & - \varepsilon \frac{\beta^n}{\Delta y} \left( \underline{Y}_\eta^1 \nu_{i, y_j} (\tilde{\mathbf{q}}_z)_{i, y_j}^n - \underline{Y}_\eta^0 \nu_{i, y_{j-1}} (\tilde{\mathbf{q}}_z)_{i, y_{j-1}}^n \right). \end{aligned} \quad (62) $$

We draw the reader's attention to the fact that the matrix $\underline{N}$ is diagonal, which is a consequence of our choice of an orthogonal Lagrange polynomial basis. Thus, we can multiply the correction stage (62) by the matrix $\underline{N}^{-1}$ and no linear system is to be solved anymore.

This complete algorithm describing all the steps of the ADER-DG method applied to the nonlinear advection-diffusion equation is then tested against different examples for which an exact solution is known. It is the purpose of the next section 4.

# 4 Numerical results

This section is dedicated to a numerical validation of the ADER-DG space-time method applied to the two dimensional nonlinear advection-diffusion equation. We propose a theoretical and numerical comparison between the penalty and relaxation approaches regarding their orders of convergence and computational cost essentially. Additionally, we investigate what method should be preferred for different meshes, degrees of approximation $D$ or viscosity $\nu$. This complete study is illustrated by numerical validations.

## 4.1 SIPG vs. Cattaneo: theory

In this subsection, we propose a comparative numerical analysis of the SIPG and Cattaneo approaches. We highlight the advantage of using one over the other regarding the execution time, depending on different given parameters.
20 ADER-DG schemes for nonlinear advection-diffusion

From the work detailed in [11], it is clear that, among the three penalty methods, the symmetric version SIPG is always faster than its two incomplete and nonsymmetric counterparts, which guided our choice for the penalty approach.

In subsection 3.6, we remark that the different steps of the ADER-DG algorithm are exactly the same for both the SIPG and Cattaneo methods. The local linear systems to be solved are identical, and the only difference lies in the forms of stabilization coefficients (38) and (46) respectively. This difference is then reflected in the corresponding time steps computed by formula (52). More precisely, the SIPG penalty scheme is supposed to be faster than the Cattaneo relaxation approach if it requires less iterations $N_t$ of the scheme to reach the final time $T_f := \Delta t^0 + \dots + \Delta t^{N_t-1}$ or, equivalently, if it allows a larger time step. If we consider a maximum viscosity $\nu_* := \max(\nu)$, a $\mathbb{P}_D$ approximation with $D \in \mathbb{N}^*$ and a mesh size $h := N^{-1}$, this condition writes:

$$ \Delta t_d^S(\nu_*, D, h) \geq \Delta t_d^C(\nu_*, D, h). \qquad (63) $$

From the formulas (52) for the time steps, where the stabilization coefficients $\sigma_d^M$ for $M = S, C$ can be replaced by their expressions (50)-(51), we have:

$$ \Delta t_d^S(\nu_*, D, h) = \frac{\kappa_d^S(D)}{\gamma^S(D)\nu_*} \cdot h^2, \qquad (64) $$

$$ \Delta t_d^C(\nu_*, D, h) = \frac{\kappa_d^C(D)}{\sqrt{\gamma^C(D)\nu_*}} \cdot h^{(D+3)/2}. \qquad (65) $$

We emphasize that $\Delta t_d^S \propto \nu_*^{-1}$ while $\Delta t_d^C \propto (\sqrt{\nu_*})^{-1}$, where the symbol $\propto$ denotes a proportional relation. Consequently, the following two asymptotic behaviours can be identified for a fixed approximation degree $D$ and mesh size $h$: for arbitrary small (respectively high) values of maximum viscosity $\nu_*$, the SIPG (respectively Cattaneo) method executes faster due to a larger time step.

Thus, using expressions (64)-(65), inequality (63) squared to simplify this study is equivalent to:

$$ \psi_{cfl}^D(\nu_*, N) \leq 1, \qquad (66) $$

where, for a given approximation degree $D$, we have introduced the scalar function $\psi_{cfl}^D$ of two variables:

$$ \psi_{cfl}^D(\nu_*, N) := \alpha(D) \cdot \nu_* \cdot N^{1-D}, \qquad (67) $$

together with the following coefficient:

$$ \alpha(D) := \left( \frac{\kappa_d^C(D)}{\kappa_d^S(D)} \cdot \frac{\gamma^S(D)}{\sqrt{\gamma^C(D)}} \right)^2. \qquad (68) $$

After a few computations, this inequality is equivalent to the following ones, written in terms of number of cells per direction $N := h^{-1}$ or maximum
viscosity $\nu_*$:

$$N \geq \bar{N}(D, \nu_*) := (\alpha(D)\nu_*)^{1/(D-1)}, \quad (69)$$

$$\nu_* \leq \bar{\nu}_*(D, N) := \alpha(D)^{-1} N^{D-1}. \quad (70)$$

Let us comment the above inequalities (69)-(70) for a fixed $\mathbb{P}_D$ approximation, treating the cases $D=1$ and $D \geq 2$ separately:

* For $D=1$, expression (67) simplifies into: $\psi_{\text{cfl}}^1(\nu_*, N) := \alpha(1) \cdot \nu_*$, so the condition (66) can be summed up as: $\nu_* \leq \alpha(1)^{-1} = \bar{\nu}_*(1, N) \approx 2.5 \times 10^{-3}$. It is clear that this condition is satisfied or not independently of the number of cells $N$ per direction $x, y$ composing the mesh. In other words, for a $\mathbb{P}_1$ approximation, if $\nu_* \leq \bar{\nu}_*(1, N)$ then the SIPG is faster regardless of the mesh size used, otherwise this advantage goes to the Cattaneo approach. The next comments only concern the cases where $D \geq 2$.

* On the one hand, for a given maximum viscosity $\nu_*$, the relation (69) suggests that there exists a limit number of cells $\bar{N}(D, \nu_*)$ up to which the Cattaneo method provides a larger time step and thus is faster than the SIPG approach, and then for greater values of $N \geq \bar{N}(D, \nu_*)$ the SIPG method should be favoured.

* On the other hand, for a mesh of characteristic length $h := \max(\Delta x, \Delta y)$, or equivalently a mesh discretized with $N := h^{-1}$ cells per direction $x, y$, the relation (70) expresses that there exists a limit maximum viscosity $\bar{\nu}_*(D, N)$ up to which the SIPG approach is faster than the Cattaneo relaxation scheme, and then for greater values of $\nu_* \geq \bar{\nu}_*(D, N)$ the Cattaneo method becomes faster.

Moreover, both the SIPG and Cattaneo methods are supposed to be of order of convergence $O_{L^2}(\mathbb{P}_D) = D+1$ in the $L^2$ norm, for a $\mathbb{P}_D$ approximation. For the Cattaneo relaxation scheme, this accuracy is clearly obtained thanks to the choice $\delta = O(h^{D+1})$ for the relaxation parameter, as expressed through the choice (43). Still, the SIPG penalty scheme is the only IP version of order $(D+1)$ for any approximation degree $D$. The reader can refer to [11, 29] for more details, where it is precised that the incomplete and nonsymmetric versions, IIPG and NIPG, are of order $(D+1)$ for odd values of $D$ but of order $D$ only for even values of $D$.

In the next subsection 4.2, the SIPG and Cattaneo approaches are tested against two test cases respectively: a pure diffusion problem with constant viscosity and Neumann boundary conditions, and a complete advection-diffusion problem with variable viscosity and Dirichlet boundary conditions.

## 4.2 Application to a manufactured test case

We consider here a manufactured test case to show numerical illustrations of the ADER-DG penalty and relaxation methods. The space domain chosen is $\Omega := [0, 1]^2$ and the final time is $T_f = 0.1$. The exact solution on $\Omega \times [0, T_f]$ is
22 ADER-DG schemes for nonlinear advection-diffusion

given by the following expression:

$$ \mathbf{u}(x, y, t) := \mathbf{u}_0(x, y)e^{-10 \ln(2)t}, \quad \mathbf{u}_0(x, y) := \begin{bmatrix} + \cos(2\pi x) \sin(2\pi y) \\ - \sin(2\pi x) \cos(2\pi y) \end{bmatrix}, \quad (71) $$

and it is illustrated in figure 2.

<figure><img src="image_4.png" /><figcaption>Fig. 2 Analytical solution $\mathbf{u} := [u_x, u_y]^T$ evaluated at the final time $T_f$.</figcaption></figure>

The reader's attention is drawn to the fact that having a manufactured solution induces a source term that compensates all the physical phenomena, namely the advection and the diffusion. In the following, we propose two different applications for different viscosity fields or boundary conditions for example, in order to show the robustness of the arbitrary high order ADER-DG scheme.

### Pure diffusion with constant viscosity and Neumann conditions

The first application will be referred to as EX1 and deals with a pure diffusion problem, as it is the only component for which the SIPG and Cattaneo approaches differ. The constant viscosity is $\nu \equiv 10^{-2}$. The Neumann boundary conditions $H_d(\mathbf{u})|_{\partial\Omega}$ for $H = F, G$ imposed as well as the source term to be enforced $S := \partial_t \mathbf{u} - \nu \Delta \mathbf{u}$, where $\Delta \mathbf{u}$ is the Laplacian of $\mathbf{u}$, are directly deduced from the exact solution (71).

Let us now compare the SIPG penalty and Cattaneo relaxation versions of the ADER-DG method, regarding the execution time $T_{exec}$ (for a simulation in parallel with 6 threads), the $L^2$ relative error $\mathcal{E}_{L^2}$ and the order of convergence $\mathcal{O}_{L^2}$ for the $L^2$ norm. The comparison is made with $\mathbb{P}_D$ approximations for $D \in \{1, \dots, 4\}$ by using four consecutive meshes composed of $N^2$ cells, where $N := 32 \times 2^{s-D+1}$ for $s \in \{0, \dots, 3\}$. This study is grouped in table 2 where, for each approximation degree $D$ and each mesh containing $N^2$ cells, the shortest execution time $T_{exec}$ and the smallest error $\mathcal{E}_{L^2}$ between the numerical and exact solutions are shown in bold. We also report the approximative ratio between the execution times for two successive meshes.
<table><caption>**Table 2** Performances and convergence study for EX1.</caption><thead><tr><th>P<sub>D</sub></th><th>M</th><th>N<sup>2</sup></th><th>T<sub>exec</sub>(s)</th><th>ε<sub>L<sup>2</sup></sub></th><th>O<sub>L<sup>2</sup></sub></th></tr></thead><tbody><tr><td rowspan="12">P<sub>1</sub></td><td rowspan="3">S</td><td>32<sup>2</sup></td><td>1</td><td>4.87 × 10<sup>-4</sup></td><td>—</td></tr><tr><td>64<sup>2</sup></td><td>4</td><td>1.23 × 10<sup>-4</sup></td><td>1.99</td></tr><tr><td>128<sup>2</sup></td><td>55</td><td>13.8</td><td>3.08 × 10<sup>-5</sup></td><td>2.00</td></tr><tr><td rowspan="3">C</td><td>256<sup>2</sup></td><td>863</td><td>15.7</td><td>7.70 × 10<sup>-6</sup></td><td>2.00</td></tr><tr><td>32<sup>2</sup></td><td><b>1</b></td><td>—</td><td><b>2.02 × 10<sup>-4</sup></b></td><td>—</td></tr><tr><td>64<sup>2</sup></td><td><b>3</b></td><td>3.0</td><td><b>4.56 × 10<sup>-5</sup></b></td><td>2.15</td></tr><tr><td rowspan="3">P<sub>2</sub></td><td>128<sup>2</sup></td><td><b>31</b></td><td>10.4</td><td><b>1.11 × 10<sup>-5</sup></b></td><td>2.04</td></tr><tr><td>256<sup>2</sup></td><td><b>451</b></td><td>14.6</td><td><b>2.75 × 10<sup>-6</sup></b></td><td>2.01</td></tr><tr><td>16<sup>2</sup></td><td>1</td><td>—</td><td><b>3.40 × 10<sup>-4</sup></b></td><td>—</td></tr><tr><td rowspan="3">S</td><td>32<sup>2</sup></td><td>9</td><td>9.0</td><td><b>4.17 × 10<sup>-5</sup></b></td><td>3.03</td></tr><tr><td>64<sup>2</sup></td><td>136</td><td>15.2</td><td><b>5.15 × 10<sup>-6</sup></b></td><td>3.02</td></tr><tr><td>128<sup>2</sup></td><td>2 320</td><td>17.1</td><td><b>6.46 × 10<sup>-7</sup></b></td><td>2.99</td></tr><tr><td rowspan="4">P<sub>3</sub></td><td rowspan="3">C</td><td>16<sup>2</sup></td><td><b>1</b></td><td>—</td><td>3.49 × 10<sup>-4</sup></td><td>—</td></tr><tr><td>32<sup>2</sup></td><td><b>6</b></td><td>6.0</td><td>6.28 × 10<sup>-5</sup></td><td>2.47</td></tr><tr><td>64<sup>2</sup></td><td><b>88</b></td><td>14.7</td><td>1.09 × 10<sup>-5</sup></td><td>2.53</td></tr><tr><td>128<sup>2</sup></td><td><b>1 908</b></td><td>21.7</td><td>1.89 × 10<sup>-6</sup></td><td>2.53</td></tr><tr><td rowspan="12">P<sub>4</sub></td><td rowspan="3">S</td><td>8<sup>2</sup></td><td>2</td><td>—</td><td><b>1.65 × 10<sup>-4</sup></b></td><td>—</td></tr><tr><td>16<sup>2</sup></td><td>13</td><td>6.5</td><td><b>9.65 × 10<sup>-6</sup></b></td><td>4.10</td></tr><tr><td>32<sup>2</sup></td><td>216</td><td>16.7</td><td><b>5.87 × 10<sup>-7</sup></b></td><td>4.04</td></tr><tr><td rowspan="3">C</td><td>64<sup>2</sup></td><td><b>3 441</b></td><td>16.0</td><td>3.64 × 10<sup>-8</sup></td><td>4.01</td></tr><tr><td>8<sup>2</sup></td><td><b>1</b></td><td>—</td><td>2.04 × 10<sup>-4</sup></td><td>—</td></tr><tr><td>16<sup>2</sup></td><td><b>6</b></td><td>6.0</td><td>1.04 × 10<sup>-5</sup></td><td>4.29</td></tr><tr><td rowspan="3">P<sub>5</sub></td><td>32<sup>2</sup></td><td><b>146</b></td><td>24.4</td><td>5.98 × 10<sup>-7</sup></td><td>4.12</td></tr><tr><td>64<sup>2</sup></td><td>4 658</td><td>31.9</td><td><b>3.61 × 10<sup>-8</sup></b></td><td>4.05</td></tr><tr><td>4<sup>2</sup></td><td>2</td><td>—</td><td><b>1.70 × 10<sup>-4</sup></b></td><td>—</td></tr><tr><td rowspan="3">S</td><td>8<sup>2</sup></td><td>16</td><td>8.0</td><td><b>5.13 × 10<sup>-6</sup></b></td><td>5.05</td></tr><tr><td>16<sup>2</sup></td><td>221</td><td>13.9</td><td><b>1.54 × 10<sup>-7</sup></b></td><td>5.06</td></tr><tr><td>32<sup>2</sup></td><td><b>4 602</b></td><td>20.9</td><td><b>4.70 × 10<sup>-9</sup></b></td><td>5.03</td></tr><tr><td rowspan="4">P<sub>6</sub></td><td rowspan="3">C</td><td>4<sup>2</sup></td><td><b>2</b></td><td>—</td><td>3.50 × 10<sup>-4</sup></td><td>—</td></tr><tr><td>8<sup>2</sup></td><td><b>5</b></td><td>2.5</td><td>8.57 × 10<sup>-6</sup></td><td>5.35</td></tr><tr><td>16<sup>2</sup></td><td><b>145</b></td><td>29.0</td><td>1.93 × 10<sup>-7</sup></td><td>5.47</td></tr><tr><td>32<sup>2</sup></td><td>6 546</td><td>45.2</td><td>5.41 × 10<sup>-9</sup></td><td>5.16</td></tr></tbody></table>

The first immediate comment is that the numerical solution effectively converges to the manufactured exact solution (71). More precisely, for a $P_D$ approximation, the order of convergence in the $L^2$ norm is indeed $O_{L^2}(P_D) = D + 1$ for both the SIPG and Cattaneo versions, as expected in theory. We can however note that the order of convergence of the Cattaneo method with a $P_2$ scheme is closer to 2.5 rather than 3, but it does not strongly contradict the overall validity of the results obtained. Then in terms of performances, we can comment the execution times and the errors separately:
24 ADER-DG schemes for nonlinear advection-diffusion

• Concerning the errors $\mathcal{E}_L^2$ obtained between the numerical and analytical solutions, SIPG is systematically more accurate than Cattaneo for the viscosity and the range of meshes used, except for $\mathbb{P}_1$ approximations where the opposite is true.

• Regarding the execution times $T_{\text{exec}}$ recorded, the comparison strictly follows the analysis provided in subsection 4.1. Indeed, as $\nu_*=10^{-2} \ge \bar{\nu}_*(1, N) \approx 2.5 \times 10^{-3}$ in our case, the Cattaneo method is always faster independently of the mesh used. Moreover, the limit numbers of cells are $\bar{N}(2, 10^{-2}) = 195$, $\bar{N}(3, 10^{-2}) = 48$ and $\bar{N}(4, 10^{-2}) = 26$ from formula (69), and we recall that SIPG is faster than Cattaneo if and only if $N \ge \bar{N}(D, \nu_*)$. That is why with the meshes used in our convergence study, the SIPG scheme is more advantageous only for the finest mesh for $\mathbb{P}_3$ and $\mathbb{P}_4$ approximations, where $N = 64 > 48$ and $N = 32 > 26$ respectively, but the Cattaneo method is always faster up to the mesh with $N = 128 < 195$ in the $\mathbb{P}_2$ case.

• When a small relative error is required, the higher order schemes should be favoured as they run faster with a lower memory footprint regarding the number of unknowns stored. For instance, if we aim at reaching an error close to $10^{-7}$, it takes about 200 seconds in the $\mathbb{P}_4$ case or 2 000 seconds in the $\mathbb{P}_2$ case.

In conclusion, it is important to carry out a preliminary study before choosing the SIPG or Cattaneo methods. Indeed, if they both provide close errors for a given mesh, their execution times can differ a lot. For instance, with a $\mathbb{P}_1$ approximation and a mesh composed of $N^2 = 256^2$ cells, the execution with SIPG takes about 90% more time than Cattaneo, which has a limited impact because the respective times are around 14 and 7 minutes. However, with a $\mathbb{P}_4$ approximation and a mesh composed of $N^2 = 32^2$ cells, the execution with Cattaneo takes about 40% more time than SIPG, but the corresponding times differ by 35 minutes, which is more substantial. This preliminary study is even more important as the mesh is refined because the ratio between two successive execution times is around 15 for SIPG but keeps increasing for Cattaneo. We close this analysis by recalling that all the behaviours pointed out here directly depend on the maximum viscosity $\nu_*$ through the CFL condition (52).

Advection-diffusion with variable viscosity and Dirichlet conditions

The second application presented here corresponds to a complete advection-diffusion problem, and it will be referred to as EX2. Again, the Cattaneo and SIPG methods are compared in the same way as for EX1. However, for a more challenging test case, we consider a viscosity field which depends on the space variables, and Dirichlet boundary conditions are imposed. The source term $S := \partial_t u + \nabla \cdot (u \otimes u - \nu \nabla u)$ and the boundary conditions are obtained from the exact solution (71).

If we introduce the notation $r(x,y) := \sqrt{(x-0.5)^2 + (y-0.5)^2}$ corresponding to distance of the point $(x,y) \in \Omega$ to the centre point $(0.5, 0.5)$, then
the viscosity function chosen is defined by the smoothed circular step below:

$$ \nu(x, y) := \nu_{\text{out}} + (\nu_{\text{in}} - \nu_{\text{out}}) \exp(-10^3 \cdot r(x, y)^6), \quad (72) $$

with the values $\nu_{\text{in}} := 10^{-2}$ and $\nu_{\text{out}} := 10^{-3}$.

Again, the execution time, the relative error and the order of convergence are carefully examined to compare the Cattaneo and SIPG methods. The same meshes composed of $N := 32 \times 2^{s-D+1}$ cells per direction $x, y$ for $D \in \{1, \dots, 4\}$ and $s \in \{0, \dots, 3\}$ are used on EX2. The numerical results are shown in table 3, where we have written in bold the smallest error and execution time (with 6 threads in parallel) for each $\mathbb{P}_D$ approximation.

The comments for this EX2 are very close to those we made for EX1. In particular, the expected orders of accuracy, that is to say $O_{L^2}(\mathbb{P}_D) = (D+1)$, are globally respected and the numerical solutions converge to the analytical solution (71) for the two methods considered. However, both the SIPG and Cattaneo approaches seem to lead to a super convergence in some cases. For instance, a $\mathbb{P}_1$ scheme with the SIPG version is closer to an order 2.5 or 3 instead of 2, and the Cattaneo method with a $\mathbb{P}_4$ approximation may exceed the seventh order. Nevertheless, the ADER-DG method and its associated orders of convergence are mainly to be studied asymptotically in the limit of mesh refinement. From this perspective, we can see that any $\mathbb{P}_D$ discretization tends to an order $(D+1)$ for the finer meshes. In the particular case $D=4$, the asymptotic region may not be reached yet with the meshes used. Thus, two consecutive meshes may substantially improve the polynomial representation of the stiff viscosity function (72) and lead to the observed super convergence. More precisely, we have simulated this EX2 with a constant viscosity, or without advection, which converged faster to the expected order. All the other comments made for EX1 can be transposed to EX2, which concludes the numerical validation of the ADER-DG scheme presented in its two variants: the SIPG penalty and the Cattaneo relaxation approaches.

# 5 Conclusion and perspectives

In this paper, we carefully presented all the steps of an ADER-DG scheme applied to solve the two dimensional unsteady nonlinear advection-diffusion equation on a cartesian grid. Additionally, we proposed a comparison between a penalty and a relaxation approach, namely the SIPG and Cattaneo methods, in terms of error, execution time and order of accuracy. This study was completed by numerical illustrations through two test cases based on a smooth manufactured solution, and it proved the robustness of our approach against any choice of high approximation degree, viscosity variations or boundary conditions type. More precisely, it appears that for a $\mathbb{P}_1$ approximation, the maximum viscosity $\nu_*$ entirely dictates whether the Cattaneo or SIPG method should be favoured, independently of the mesh. However, for $\mathbb{P}_D$ approximations with $D \ge 2$, the SIPG scheme becomes faster than the Cattaneo method
26 ADER-DG schemes for nonlinear advection-diffusion

<table><caption>**Table 3** Performances and convergence study for EX2.</caption><thead><tr><th>P<sub>D</sub></th><th>M</th><th>N<sup>2</sup></th><th>T<sub>exec</sub></th><th>×</th><th>E<sub>L<sup>2</sup></sub></th><th>O<sub>L<sup>2</sup></sub></th></tr></thead><tbody><tr><td rowspan="12">P<sub>1</sub></td><td rowspan="3">S</td><td>32<sup>2</sup></td><td>2</td><td>—</td><td>3.80 × 10<sup>-3</sup></td><td>—</td></tr><tr><td>64<sup>2</sup></td><td>6</td><td>3.0</td><td><b>5.11 × 10<sup>-4</sup></b></td><td>2.89</td></tr><tr><td>128<sup>2</sup></td><td>89</td><td>14.9</td><td><b>7.53 × 10<sup>-5</sup></b></td><td>2.76</td></tr><tr><td></td><td>256<sup>2</sup></td><td>1 379</td><td>15.5</td><td><b>1.16 × 10<sup>-5</sup></b></td><td>2.70</td></tr><tr><td rowspan="4">C</td><td>32<sup>2</sup></td><td><b>1</b></td><td>—</td><td><b>3.13 × 10<sup>-3</sup></b></td><td>—</td></tr><tr><td>64<sup>2</sup></td><td><b>6</b></td><td>6.0</td><td>5.94 × 10<sup>-4</sup></td><td>2.40</td></tr><tr><td>128<sup>2</sup></td><td><b>52</b></td><td>8.7</td><td>1.40 × 10<sup>-4</sup></td><td>2.09</td></tr><tr><td>256<sup>2</sup></td><td><b>724</b></td><td>14.0</td><td>3.54 × 10<sup>-5</sup></td><td>1.98</td></tr><tr><td rowspan="12">P<sub>2</sub></td><td rowspan="3">S</td><td>16<sup>2</sup></td><td>2</td><td>—</td><td>4.84 × 10<sup>-4</sup></td><td>—</td></tr><tr><td>32<sup>2</sup></td><td>14</td><td>7.0</td><td><b>4.78 × 10<sup>-5</sup></b></td><td>3.34</td></tr><tr><td>64<sup>2</sup></td><td>240</td><td>17.2</td><td><b>5.23 × 10<sup>-6</sup></b></td><td>3.19</td></tr><tr><td></td><td>128<sup>2</sup></td><td>2 936</td><td>12.3</td><td><b>6.26 × 10<sup>-7</sup></b></td><td>3.06</td></tr><tr><td rowspan="4">C</td><td>16<sup>2</sup></td><td><b>1</b></td><td>—</td><td><b>4.47 × 10<sup>-4</sup></b></td><td>—</td></tr><tr><td>32<sup>2</sup></td><td><b>8</b></td><td>8.0</td><td>5.07 × 10<sup>-5</sup></td><td>3.14</td></tr><tr><td>64<sup>2</sup></td><td><b>120</b></td><td>15.0</td><td>6.89 × 10<sup>-6</sup></td><td>2.88</td></tr><tr><td>128<sup>2</sup></td><td><b>2 420</b></td><td>20.2</td><td>1.01 × 10<sup>-6</sup></td><td>2.77</td></tr><tr><td rowspan="12">P<sub>3</sub></td><td rowspan="3">S</td><td>8<sup>2</sup></td><td>3</td><td>—</td><td><b>6.39 × 10<sup>-4</sup></b></td><td>—</td></tr><tr><td>16<sup>2</sup></td><td>19</td><td>6.4</td><td><b>1.83 × 10<sup>-5</sup></b></td><td>5.13</td></tr><tr><td>32<sup>2</sup></td><td>308</td><td>16.3</td><td>6.75 × 10<sup>-7</sup></td><td>4.76</td></tr><tr><td></td><td>64<sup>2</sup></td><td><b>4 614</b></td><td>15.0</td><td><b>3.79 × 10<sup>-8</sup></b></td><td>4.15</td></tr><tr><td rowspan="4">C</td><td>8<sup>2</sup></td><td><b>2</b></td><td>—</td><td>1.21 × 10<sup>-3</sup></td><td>—</td></tr><tr><td>16<sup>2</sup></td><td><b>10</b></td><td>5.0</td><td>2.15 × 10<sup>-5</sup></td><td>5.81</td></tr><tr><td>32<sup>2</sup></td><td><b>221</b></td><td>22.1</td><td><b>5.74 × 10<sup>-7</sup></b></td><td>5.23</td></tr><tr><td>64<sup>2</sup></td><td>5 851</td><td>26.5</td><td>4.05 × 10<sup>-8</sup></td><td>3.83</td></tr><tr><td rowspan="12">P<sub>4</sub></td><td rowspan="3">S</td><td>4<sup>2</sup></td><td>3</td><td>—</td><td><b>2.21 × 10<sup>-3</sup></b></td><td>—</td></tr><tr><td>8<sup>2</sup></td><td>20</td><td>6.7</td><td><b>2.76 × 10<sup>-5</sup></b></td><td>6.32</td></tr><tr><td>16<sup>2</sup></td><td>370</td><td>18.5</td><td>3.80 × 10<sup>-7</sup></td><td>6.18</td></tr><tr><td></td><td>32<sup>2</sup></td><td><b>5 617</b></td><td>15.2</td><td><b>5.67 × 10<sup>-9</sup></b></td><td>6.07</td></tr><tr><td rowspan="4">C</td><td>4<sup>2</sup></td><td><b>2</b></td><td>—</td><td>1.02 × 10<sup>-2</sup></td><td>—</td></tr><tr><td>8<sup>2</sup></td><td><b>10</b></td><td>5.0</td><td>4.71 × 10<sup>-5</sup></td><td>7.76</td></tr><tr><td>16<sup>2</sup></td><td><b>199</b></td><td>19.9</td><td><b>3.67 × 10<sup>-7</sup></b></td><td>7.00</td></tr><tr><td>32<sup>2</sup></td><td>7 601</td><td>38.2</td><td>6.10 × 10<sup>-9</sup></td><td>5.91</td></tr></tbody></table>

for finer mesh sizes or vanishing viscosities, while the Cattaneo method should be preferred for coarser grids or more viscous flows.

The advantages of the ADER-DG framework open up new prospects. On the one hand, the manoeuvrability of the local polynomial expansions of the solution and their easy evaluation at cells interfaces to compute the numerical fluxes suggest that the ADER-DG setting could be extended to quadtree meshes with a reasonable theoretical and implementation effort. On the other hand, the validation of the present work on the nonlinear advection-diffusion equation offers an interesting perspective: solving the incompressible Navier-Stokes equations. Indeed, this system of PDEs is generally solved with the
Chorin-Temam approach [18, 33, 34] where a nonlinear advection-diffusion equation is to be solved, so the present work clearly provides an additional step to reach the aim to solve the incompressible Navier-Stokes equations on quadtree grids.

# Acknowledgements

This work was funded by the University of Bordeaux.

# Declarations

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

# In memoriam

This paper is dedicated to the memory of Prof. Arturo Hidalgo López (*July 03rd 1966 - †August 26th 2024) of the Universidad Politecnica de Madrid, organizer of HONOM 2019 and active participant in many other editions of HONOM. Our thoughts and wishes go to his wife Lourdes and his sister María Jesús, whom he left behind.

# References

[1] Bergmann, M., Carlino, M.G., Iollo, A.: Second Order ADER Scheme for Unsteady Advection-Diffusion on Moving Overset Grids with a Compact Transmission Condition. SIAM Journal on Scientific Computing 44(1), 524–553 (2022). https://doi.org/10.1137/21M1393911

[2] Loubère, R., Dumbser, M., Diot, S.: A New Family of High Order Unstructured MOOD and ADER Finite Volume Schemes for Multidimensional Systems of Hyperbolic Conservation Laws. Communications in Computational Physics 16(3), 718–763 (2014). https://doi.org/10.4208/cicp.181113.140314a

[3] Gaburro, E., Dumbser, M., Castro, M.J.: Direct Arbitrary-Lagrangian-Eulerian finite volume schemes on moving nonconforming unstructured meshes. Computers & Fluids 159, 254–275 (2017). https://doi.org/10.1016/j.compfluid.2017.09.022

[4] Río-Martín, L., Busto, S., Dumbser, M.: A Massively Parallel Hybrid Finite Volume/Finite Element Scheme for Computational Fluid Dynamics. Mathematics 9(18), 2316 (2021). https://doi.org/10.3390/math9182316
28 ADER-DG schemes for nonlinear advection-diffusion

[5] Boscheri, W., Dumbser, M.: Arbitrary-Lagrangian-Eulerian Discontinuous Galerkin schemes with a posteriori subcell finite volume limiting on moving unstructured meshes. Journal of Computational Physics **346**, 449–479 (2017). https://doi.org/10.1016/j.jcp.2017.06.022

[6] Baldassari, C.: Modélisation et simulation numérique pour la migration terrestre par équation d'ondes. Theses, Université de Pau et des Pays de l'Adour (December 2009). https://theses.hal.science/tel-00472810

[7] Cockburn, B., Shu, C.-W.: Runge-kutta discontinuous galerkin methods for convection-dominated problems. Journal of scientific computing **16**, 173–261 (2001)

[8] Cockburn, B., Karniadakis, G.E., Shu, C.-W., Griebel, M., Keyes, D.E., Nieminen, R.M., Roosevelt, D., Schlick, T. (eds.): Discontinuous Galerkin Methods: Theory, Computation and Applications. Lecture Notes in Computational Science and Engineering, vol. 11. Springer, Berlin, Heidelberg (2000). https://doi.org/10.1007/978-3-642-59721-3. https://link.springer.com/10.1007/978-3-642-59721-3

[9] Cockburn, B.: Discontinuous Galerkin methods. ZAMM **83**(11), 731–754 (2003). https://doi.org/10.1002/zamm.200310088

[10] Liu, H., Yan, J.: The Direct Discontinuous Galerkin (DDG) Methods for Diffusion Problems. SIAM Journal on Numerical Analysis **47**(1), 675–698 (2009). https://doi.org/10.1137/080720255

[11] Bergmann, M., Bouharguane, A., Iollo, A., Tardieu, A.: High Order ADER-IPDG Methods for the Unsteady Advection-Diffusion Equation. Communications on Applied Mathematics and Computation **6**(3), 1954–1977 (2024). https://doi.org/10.1007/s42967-023-00355-w

[12] Busto, S., Chiocchetti, S., Dumbser, M., Gaburro, E., Peshkov, I.: High order ader schemes for continuum mechanics. Frontiers in Physics **8** (2020). https://doi.org/10.3389/fphy.2020.00032

[13] Busto, S., Dumbser, M., Escalante, C., Favrie, N., Gavrilyuk, S.: On High Order ADER Discontinuous Galerkin Schemes for First Order Hyperbolic Reformulations of Nonlinear Dispersive Systems. Journal of Scientific Computing **87**(2), 48 (2021). https://doi.org/10.1007/s10915-021-01429-8

[14] Fambri, F., Dumbser, M., Zanotti, O.: Space-time adaptive ADER-DG schemes for dissipative flows: Compressible Navier-Stokes and resistive MHD equations. Computer Physics Communications **220**, 297–318 (2017). https://doi.org/10.1016/j.cpc.2017.08.001
[15] Hidalgo, A., Dumbser, M.: ADER Schemes for Nonlinear Systems of Stiff Advection-Diffusion-Reaction Equations. Journal of Scientific Computing 48(1-3), 173-189 (2011). https://doi.org/10.1007/s10915-010-9426-6

[16] Boscheri, W., Loubère, R.: High Order Accurate Direct Arbitrary-Lagrangian-Eulerian ADER-MOOD Finite Volume Schemes for Non-Conservative Hyperbolic Systems with Stiff Source Terms. Communications in Computational Physics 21(1), 271-312 (2017). https://doi.org/10.4208/cicp.OA-2015-0024

[17] Dumbser, M., Fambri, F., Tavelli, M., Bader, M., Weinzierl, T.: Efficient Implementation of ADER Discontinuous Galerkin Schemes for a Scalable Hyperbolic PDE Engine. Axioms 7(3), 63 (2018). https://doi.org/10.3390/axioms7030063

[18] Bergmann, M., Carlino, M.G., Iollo, A., Telib, H.: ADER scheme for incompressible Navier-Stokes equations on overset grids with a compact transmission condition. Journal of Computational Physics 467, 111414 (2022). https://doi.org/10.1016/j.jcp.2022.111414

[19] Toro, E.F., Montecinos, G.I.: Advection-Diffusion-Reaction Equations: Hyperbolization and High-Order ADER Discretizations. SIAM Journal on Scientific Computing 36(5), 2423-2457 (2014). https://doi.org/10.1137/130937469

[20] Titarev, V.A., Toro, E.F.: Ader: Arbitrary high order godunov approach. Journal of Scientific Computing 17(1/4), 609-618 (2002). https://doi.org/10.1023/A:1015126814947

[21] Busto, S., Toro, E.F., Vázquez-Cendón, M.E.: Design and analysis of ADER-type schemes for model advection-diffusion-reaction equations. Journal of Computational Physics 327, 553-575 (2016). https://doi.org/10.1016/j.jcp.2016.09.043

[22] Cho, H., Park, Y., Kang, M.: Solving incompressible Navier-Stokes equations on irregular domains and quadtrees by monolithic approach. Journal of Computational Physics 463, 111304 (2022). https://doi.org/10.1016/j.jcp.2022.111304. arXiv:2108.12073 [physics]

[23] Bergmann, M., Fondanèche, A., Iollo, A.: An Eulerian finite-volume approach of fluid-structure interaction problems on quadtree meshes. Journal of Computational Physics 471, 111647 (2022). https://doi.org/10.1016/j.jcp.2022.111647

[24] Olshanskii, M.A., Terekhov, K.M., Vassilevski, Y.V.: An octree-based solver for the incompressible Navier-Stokes equations with enhanced stability and low dissipation. Computers & Fluids 84, 231-246 (2013).
30 ADER-DG schemes for nonlinear advection-diffusion

https://doi.org/10.1016/j.compfluid.2013.04.027

[25] Guittet, A., Theillard, M., Gibou, F.: A stable projection method for the incompressible Navier-Stokes equations on arbitrary geometries and adaptive Quad/Octrees. Journal of Computational Physics **292**, 215–238 (2015). https://doi.org/10.1016/j.jcp.2015.03.024

[26] Hartmann, R.: Numerical analysis of higher order discontinuous galerkin finite element methods (2008)

[27] Masri, R., Liu, C., Riviere, B.: A discontinuous Galerkin pressure correction scheme for the incompressible Navier-Stokes equations: Stability and convergence. Mathematics of Computation **91**(336), 1625–1654 (2022). https://doi.org/10.1090/mcom/3731

[28] Girault, V., Rivière, B., Wheeler, M.F.: A splitting method using discontinuous Galerkin for the transient incompressible Navier-Stokes equations. ESAIM: Mathematical Modelling and Numerical Analysis **39**(6), 1115–1147 (2005). https://doi.org/10.1051/m2an:2005048

[29] Rivière, B.: Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation. Frontiers in applied mathematics. SIAM, Society for Industrial and Applied Mathematics, Philadelphia, PA (2008). OCLC: ocn226292048

[30] Fambri, F., Dumbser, M.: Semi-implicit discontinuous Galerkin methods for the incompressible Navier-Stokes equations on adaptive staggered Cartesian grids. Computer Methods in Applied Mechanics and Engineering **324**, 170–203 (2017). https://doi.org/10.1016/j.cma.2017.06.003

[31] Fambri, F.: Discontinuous Galerkin Methods for Compressible and Incompressible Flows on Space-Time Adaptive Meshes: Toward a Novel Family of Efficient Numerical Methods for Fluid Dynamics. Archives of Computational Methods in Engineering **27**(1), 199–283 (2020). https://doi.org/10.1007/s11831-018-09308-6

[32] Gassner, G., Dumbser, M., Hindenlang, F., Munz, C.-D.: Explicit one-step time discretizations for discontinuous Galerkin and finite volume schemes based on local predictors. Journal of Computational Physics **230**(11), 4232–4247 (2011). https://doi.org/10.1016/j.jcp.2010.10.024

[33] Bergmann, M., Iollo, A.: Modeling and simulation of fish-like swimming. Journal of Computational Physics **230**(2), 329–348 (2011). https://doi.org/10.1016/j.jcp.2010.09.017

[34] Bergmann, M., Hovnanian, J., Iollo, A.: An Accurate Cartesian Method for Incompressible Flows with Moving Boundaries. Communications
in Computational Physics 15(5), 1266-1290 (2014). https://doi.org/10.4208/cicp.220313.111013a