# High order ADER-IPDG methods for the unsteady advection-diffusion equation

Michel Bergmann, Afaf Bouharguane, Angelo Iollo, Alexis Tardieu

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

## To cite this version:

Michel Bergmann, Afaf Bouharguane, Angelo Iollo, Alexis Tardieu. High order ADER-IPDG methods for the unsteady advection-diffusion equation. 2024. hal-04032238v2

HAL Id: hal-04032238

https://inria.hal.science/hal-04032238v2

Preprint submitted on 15 Jan 2024

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
High order ADER-IPDG methods for the
unsteady advection-diffusion equation

Michel Bergmann<sup>2,3</sup>, 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: michel.bergmann@u-bordeaux.fr; afaf.bouharguane@u-bordeaux.fr; angelo.iollo@u-bordeaux.fr;

**Abstract**

We present a high-order Galerkin method in both space and time for the one-dimensional unsteady linear advection-diffusion equation. Three Interior Penalty Discontinuous Galerkin (IPDG) schemes are detailed for the space discretization, while the time integration is performed at the same order of accuracy thanks to an Arbitrary high order DERivatives (ADER) method. The orders of convergence of the three ADER-IPDG methods are carefully examined through numerical illustrations, showing that the approach is consistent, accurate and efficient. The numerical results indicate that the symmetric version of IPDG is typically more accurate and more efficient compared to the other approaches.

**Keywords:** Advection-diffusion, Galerkin, ADER approach, IPDG, high-order schemes, empirical convergence rates.
High order ADER-IPDG methods for the
unsteady advection-diffusion equation

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

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

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

**Abstract**

We present a high-order Galerkin method in both space and time for the one-dimensional unsteady linear advection-diffusion equation. Three Interior Penalty Discontinuous Galerkin (IPDG) schemes are detailed for the space discretization, while the time integration is performed at the same order of accuracy thanks to an Arbitrary high order DERivatives (ADER) method. The orders of convergence of the three ADER-IPDG methods are carefully examined through numerical illustrations, showing that the approach is consistent, accurate and efficient. The numerical results indicate that the symmetric version of IPDG is typically more accurate and more efficient compared to the other approaches.

**Keywords:** Advection-diffusion, Galerkin, ADER approach, IPDG, high-order schemes, empirical convergence rates.

<footnote n="1"><sup>1</sup>Université de Bordeaux, Bordeaux, France.</footnote>

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

<footnote n="3"><sup>3</sup>Institut de Mathématiques de Bordeaux, UMR CNRS 5251, Talence, France.</footnote>
## 2 *High order ADER-IPDG methods for advection-diffusion*

# 1 Introduction

The numerical integration of the unsteady advection-diffusion equation has been widely studied by using a large variety of schemes. Finite Differences (FD), Finite Volumes (FV) or Finite Elements (FE) methods are the most common space discretization techniques, each with advantages and drawbacks especially for reaching a high order accuracy with compact stencils. The time integration is generally achieved thanks to Runge-Kutta (RK) techniques. A compromise between the upwind FV schemes and the FE methods is the Residual Distribution Scheme which can be coupled to a RK integration in time [1, 2]. Halfway between the high accuracy of continuous polynomial basis offered by classical FE methods and the local formulation of a FV scheme over a space-time cell, the Discontinuous Galerkin (DG) schemes [3–5] represent the numerical solution by a high degree polynomial in each cell.

Thanks to its local space-time high-order formulation, the Arbitrary high-order DERivatives-FV (ADER-FV) scheme [6, 7] proved accurate and efficient on a large variety of problems. In addition, the Multidimensional Optimal Order Detection (MOOD) paradigm has been developed in the last decade. It allows to detect the optimal order to ensure the preservation of the physical properties of the solution, giving birth to high order ADER-MOOD schemes for hyperbolic systems of conservation laws [8, 9]. On the other hand, the ADER-DG combination is particularly suitable to reach a desirable accuracy at the interface between overset grids or at level jumps with hierarchical meshes like quadtrees or octrees. In this sense, the work presented in [10] and [11] has shown promising results respectively for the advection-diffusion equation and the incompressible Navier-Stokes equations using a second or third order ADER-FV scheme on Chimera meshes. In our ongoing research, the plan is to set an additional step to extend that work to improve accuracy at quadtrees or octrees level jumps (see for example [12]) in hierarchical meshes. However, to attain a higher order precision at level jumps with a compact stencil for improved parallelism, we switch from the ADER-FV to the ADER-DG framework.

The Interior Penalty Discontinuous Galerkin (IPDG) scheme is a combination of the DG framework and the Interior Penalty (IP) methods, the latter being detailed in [13, 14] and applied to a wide range of equations [13, 15]. As in the DG framework, the solution is approximated within each cell using a set of basis functions. Since the solution is possibly discontinuous across the edges of the cells, the IPDG method introduces a penalty term at the interfaces between the cells to ensure the continuity of the solution and its gradient.

The novelty of the approach presented in this paper is to combine the IPDG method to handle the diffusion term and the ADER method to perform the space-time integration. This can result in a method that is accurate, efficient, and can deal with a wide range of advection-diffusion problems. More precisely, this article proposes a comparison between the classical ADER-DG scheme used in [9, 16–22], which is equivalent to the IIPG method, and its symmetric and anti-symmetric variants, named respectively SIPG and NIPG methods.
## 2 *High order ADER-IPDG methods for advection-diffusion*

# 1 Introduction

The numerical integration of the unsteady advection-diffusion equation has been widely studied by using a large variety of schemes. Finite Differences (FD), Finite Volumes (FV) or Finite Elements (FE) methods are the most common space discretization techniques, each with advantages and drawbacks especially for reaching a high order accuracy with compact stencils. The time integration is generally achieved thanks to Runge-Kutta (RK) techniques. A compromise between the upwind FV schemes and the FE methods is the Residual Distribution Scheme which can be coupled to a RK integration in time [1, 2]. Halfway between the high accuracy of continuous polynomial basis offered by classical FE methods and the local formulation of a FV scheme over a space-time cell, the Discontinuous Galerkin (DG) schemes [3–5] represent the numerical solution by a high degree polynomial in each cell.

Thanks to its local space-time high-order formulation, the Arbitrary high-order DERivatives-FV (ADER-FV) scheme [6, 7] proved accurate and efficient on a large variety of problems. In addition, the Multidimensional Optimal Order Detection (MOOD) paradigm has been developed in the last decade. It allows to detect the optimal order to ensure the preservation of the physical properties of the solution, giving birth to high order ADER-MOOD schemes for hyperbolic systems of conservation laws [8, 9]. On the other hand, the ADER-DG combination is particularly suitable to reach a desirable accuracy at the interface between overset grids or at level jumps with hierarchical meshes like quadtrees or octrees. In this sense, the work presented in [10] and [11] has shown promising results respectively for the advection-diffusion equation and the incompressible Navier-Stokes equations using a second or third order ADER-FV scheme on Chimera meshes. In our ongoing research, the plan is to set an additional step to extend that work to improve accuracy at quadtrees or octrees level jumps (see for example [12]) in hierarchical meshes. However, to attain a higher order precision at level jumps with a compact stencil for improved parallelism, we switch from the ADER-FV to the ADER-DG framework.

The Interior Penalty Discontinuous Galerkin (IPDG) scheme is a combination of the DG framework and the Interior Penalty (IP) methods, the latter being detailed in [13, 14] and applied to a wide range of equations [13, 15]. As in the DG framework, the solution is approximated within each cell using a set of basis functions. Since the solution is possibly discontinuous across the edges of the cells, the IPDG method introduces a penalty term at the interfaces between the cells to ensure the continuity of the solution and its gradient.

The novelty of the approach presented in this paper is to combine the IPDG method to handle the diffusion term and the ADER method to perform the space-time integration. This can result in a method that is accurate, efficient, and can deal with a wide range of advection-diffusion problems. More precisely, this article proposes a comparison between the classical ADER-DG scheme used in [9, 16–22], which is equivalent to the IIPG method, and its symmetric and anti-symmetric variants, named respectively SIPG and NIPG methods.
Section 2 introduces the important notations as well as the space-time structure used in this paper. In section 3, we present three ADER-IPDG prediction-correction schemes carefully describing the numerical fluxes chosen at interfaces between cells. In section 4, we validate our implementations by numerical illustrations. Section 5 is devoted to some conclusions and perspectives of our work.

# 2 Preliminary definitions and notations

## 2.1 The advection-diffusion equation

Let us consider a function $u : \Omega \times [0, T_f] \to \mathbb{R}$, where $\Omega := [x^{-}, x^{+}] \subset \mathbb{R}$, as well as the constants $a \in \mathbb{R}$ and $\lambda \in \mathbb{R}^{+}$ standing for the advection velocity and the diffusion coefficient respectively. On the same space-time domain as the unknown $u$, we denote by $S$ the source term and $F := F_a + F_d$ the flux function that is expressed as the sum of its advective and diffusive components, respectively defined as

$$F_a(u) := au, \quad F_d(u) := -\lambda \partial_x u. \qquad (1)$$

Thus, the advection-diffusion equation completed with an initial condition as well as Dirichlet boundary conditions reads

$$\left\{ \begin{array}{ll} \partial_t u + \partial_x F(u) = S & \text{on } \Omega \times [0, T_f], \\ u(x, 0) = u_0(x) & \text{on } \Omega, \\ u(x^{\pm}, t) = g_{\pm}(t) & \text{on } \partial\Omega \times [0, T_f]. \end{array} \right. \qquad (2)$$

The space and time domains $\Omega$ and $[0, T_f]$ are discretized respectively into $N_x$ uniform intervals $\omega_i$ of size $\Delta x$ and $N_t$ uniform intervals $T^n$ of size $\Delta t$ such that

$$\Omega := \bigcup_{i=1}^{N_x} \omega_i, \qquad \omega_i := [x_{i-1}, x_i], \qquad x_i := x^{-} + i\Delta x, \qquad (3a)$$

$$[0, T_f] := \bigcup_{n=0}^{N_t-1} T^n, \qquad T^n := [t^n, t^{n+1}], \qquad t^n := n\Delta t. \qquad (3b)$$

On a given space-time cell $c_i^n := \omega_i \times T^n$, we consider the usual change of variables

$$x(\xi) := x_i + \xi \Delta x, \qquad \xi \in \hat{\omega} := [0, 1], \qquad (4a)$$

$$t(\tau) := t^n + \tau \Delta t, \qquad \tau \in \hat{T} := [0, 1], \qquad (4b)$$

which allows to rewrite equation (2) with respect to the reference variables $\xi \in \hat{\omega}$ and $\tau \in \hat{T}$ instead of the physical variables $x \in \omega_i$ and $t \in T^n$. By
Section 2 introduces the important notations as well as the space-time structure used in this paper. In section 3, we present three ADER-IPDG prediction-correction schemes carefully describing the numerical fluxes chosen at interfaces between cells. In section 4, we validate our implementations by numerical illustrations. Section 5 is devoted to some conclusions and perspectives of our work.

# 2 Preliminary definitions and notations

## 2.1 The advection-diffusion equation

Let us consider a function $u: \Omega \times [0, T_f] \to \mathbb{R}$, where $\Omega := [x^{-}, x^{+}] \subset \mathbb{R}$, as well as the constants $a \in \mathbb{R}$ and $\lambda \in \mathbb{R}^{+}$ standing for the advection velocity and the diffusion coefficient respectively. On the same space-time domain as the unknown $u$, we denote by $S$ the source term and $F := F_a + F_d$ the flux function that is expressed as the sum of its advective and diffusive components, respectively defined as

$$F_a(u) := au, \quad F_d(u) := -\lambda \partial_x u. \qquad (1)$$

Thus, the advection-diffusion equation completed with an initial condition as well as Dirichlet boundary conditions reads

$$\left\{ \begin{array}{ll} \partial_t u + \partial_x F(u) = S & \text{on } \Omega \times [0, T_f], \\ u(x, 0) = u_0(x) & \text{on } \Omega, \\ u(x^{\pm}, t) = g_{\pm}(t) & \text{on } \partial\Omega \times [0, T_f]. \end{array} \right. \qquad (2)$$

The space and time domains $\Omega$ and $[0, T_f]$ are discretized respectively into $N_x$ uniform intervals $\omega_i$ of size $\Delta x$ and $N_t$ uniform intervals $T^n$ of size $\Delta t$ such that

$$\Omega := \bigcup_{i=1}^{N_x} \omega_i, \qquad \omega_i := [x_{i-1}, x_i], \qquad x_i := x^{-} + i\Delta x, \qquad (3a)$$

$$[0, T_f] := \bigcup_{n=0}^{N_t-1} T^n, \qquad T^n := [t^n, t^{n+1}], \qquad t^n := n\Delta t. \qquad (3b)$$

On a given space-time cell $c_i^n := \omega_i \times T^n$, we consider the usual change of variables

$$x(\xi) := x_i + \xi \Delta x, \qquad \xi \in \hat{\omega} := [0, 1], \qquad (4a)$$

$$t(\tau) := t^n + \tau \Delta t, \qquad \tau \in \hat{T} := [0, 1], \qquad (4b)$$

which allows to rewrite equation (2) with respect to the reference variables $\xi \in \hat{\omega}$ and $\tau \in \hat{T}$ instead of the physical variables $x \in \omega_i$ and $t \in T^n$. By
4 High order ADER-IPDG methods for advection-diffusion

using the chain rule for the derivative of composite functions, and adding a hat on quantities that are defined in the reference coordinates, we can define the following advective and diffusive fluxes

$$
\hat{F}_a(\hat{u}) := a\hat{u}, \quad \hat{F}_d(\hat{u}) := -\frac{\lambda}{\Delta x} \partial_{\xi} \hat{u}, \qquad (5)
$$

and we still have $\hat{F} := \hat{F}_a + \hat{F}_d$. As a consequence, the advection-diffusion defined on $c_i^n$ can be expressed in the reference variables on $\hat{c} := \hat{\omega} \times \hat{T}$ as follows

$$
\partial_{\tau} \hat{u} + \left( \frac{\Delta t}{\Delta x} \right) \partial_{\xi} \hat{F}(\hat{u}) = \Delta t \hat{S}. \qquad (6)
$$

The main advantage of working on this formulation at this early stage (see
[10]) is that the space and time discretizations can be entirely developed in
the reference space-time cell $\hat{c}$, which avoids making change of variables in
integrals as it is generally done for classical finite elements methods.

2.2 Reference space-time structure

Generally speaking, Galerkin methods are based on a high degree piecewise
polynomial representation of the numerical solution on the space domain Ω,
while the high order in time is reached by using RK schemes for instance. In the
particular case of the DG framework, the solution may be discontinuous across
interfaces which implies that the constraints on the localization of the degrees
of freedom disappear. As a consequence, a remarkable interest in combining the
DG discretization in space with the ADER approach in time is the possibility
to build a space-time structure as a tensor product of the polynomial basis
with respect to the space and time variables. In this paper, we use the 1D
Lagrange polynomials θ₁, ..., θK ∈ PK−₁ (the vector space of polynomials of
degree d ≤ K − 1) passing through the 1D Gauss-Legendre quadrature points
denoted by p₁, ..., pK. They write

$$
\theta_s(p) := \prod_{k \neq s} \left( \frac{p - p_k}{p_s - p_k} \right), \quad p = \xi \in \hat{\omega}, \tau \in \hat{T}. \qquad (7)
$$

The 2D space-time Lagrange basis functions Θ₁, ..., Θ_L (with L := K²) pass-
ing through the 2D quadrature points π₁, ..., π_L are then defined as products
of the 1D polynomials introduced in (7)

$$
\Theta_{\ell}(\xi, \tau) := \theta_{s}(\xi)\theta_{s'}(\tau), \quad \ell := \mathbf{m}(s, s'), \qquad (8)
$$

where the mapping **m** between local 1D indices (s, s') and the 2D global index
ℓ is defined as

$$
\mathbf{m} : \begin{cases} \{1, \dots, K\}^2 & \rightarrow \{1, \dots, L\} \\ (s, s') & \mapsto \ell := K(s' - 1) + s. \end{cases} \tag{9}
$$
4 High order ADER-IPDG methods for advection-diffusion

using the chain rule for the derivative of composite functions, and adding a hat on quantities that are defined in the reference coordinates, we can define the following advective and diffusive fluxes

$$
\hat{F}_a(\hat{u}) := a\hat{u}, \quad \hat{F}_d(\hat{u}) := -\frac{\lambda}{\Delta x} \partial_{\xi} \hat{u}, \qquad (5)
$$

and we still have $\hat{F} := \hat{F}_a + \hat{F}_d$. As a consequence, the advection-diffusion defined on $c_i^n$ can be expressed in the reference variables on $\hat{c} := \hat{\omega} \times \hat{T}$ as follows

$$
\partial_{\tau} \hat{u} + \left( \frac{\Delta t}{\Delta x} \right) \partial_{\xi} \hat{F}(\hat{u}) = \Delta t \hat{S}. \qquad (6)
$$

The main advantage of working on this formulation at this early stage (see
[10]) is that the space and time discretizations can be entirely developed in
the reference space-time cell $\hat{c}$, which avoids making change of variables in
integrals as it is generally done for classical finite elements methods.

2.2 Reference space-time structure

Generally speaking, Galerkin methods are based on a high degree piecewise
polynomial representation of the numerical solution on the space domain Ω,
while the high order in time is reached by using RK schemes for instance. In the
particular case of the DG framework, the solution may be discontinuous across
interfaces which implies that the constraints on the localization of the degrees
of freedom disappear. As a consequence, a remarkable interest in combining the
DG discretization in space with the ADER approach in time is the possibility
to build a space-time structure as a tensor product of the polynomial basis
with respect to the space and time variables. In this paper, we use the 1D
Lagrange polynomials θ₁, ..., θK ∈ PK−₁ (the vector space of polynomials of
degree d ≤ K − 1) passing through the 1D Gauss-Legendre quadrature points
denoted by p₁, ..., pK. They write

$$
\theta_s(p) := \prod_{k \neq s} \left( \frac{p - p_k}{p_s - p_k} \right), \quad p = \xi \in \hat{\omega}, \tau \in \hat{T}. \qquad (7)
$$

The 2D space-time Lagrange basis functions Θ₁, ..., Θ_L (with L := K²) pass-
ing through the 2D quadrature points π₁, ..., π_L are then defined as products
of the 1D polynomials introduced in (7)

$$
\Theta_{\ell}(\xi, \tau) := \theta_{s}(\xi)\theta_{s'}(\tau), \quad \ell := \mathbf{m}(s, s'), \qquad (8)
$$

where the mapping **m** between local 1D indices (s, s') and the 2D global index
ℓ is defined as

$$
\mathbf{m} : \begin{cases} \{1, \dots, K\}^2 & \rightarrow \{1, \dots, L\} \\ (s, s') & \mapsto \ell := K(s' - 1) + s. \end{cases} \tag{9}
$$
For instance, the $l$-th 2D quadrature point is defined by $\pi_l := (p_s, p_{s'}) \in \hat{c}$ where we have taken $l := m(s, s')$. The case $K = 3$ is illustrated in figure 1.

<figure><img src="image_1.png" /><figcaption>Fig. 1 2D space-time quadrature points $\pi_l$ obtained as tensor products of the space and time 1D quadrature points $p_k$.</figcaption></figure>

This structure is particularly convenient as its extension from 1D to 2D or 3D problems in space is straightforward for structured grids. Indeed, the 2D or 3D basis functions in space can be defined and implemented as the products of the 1D basis functions with respect to the two or three space variables, the time variable being treated exactly in the same way as for the 1D case. This remark on structured grids stands not only for regular meshes but also for quadtree or octree grids, which are the target of this work.

## 2.3 Local scalar products and matrix definitions

Let us introduce the following notations corresponding to various scalar products defined by integrals

$$ \langle f, g \rangle := \int_{\hat{\omega}} \int_{\hat{T}} f(\xi, \tau)g(\xi, \tau)d\xi d\tau, \quad (10a) $$

$$ [f, g]_{\tau} := \int_{\hat{\omega}} f(\xi, \tau)g(\xi, \tau)d\xi, \quad (10b) $$

$$ \{f, g\}_{\xi} := \int_{\hat{T}} f(\xi, \tau)g(\xi, \tau)d\tau, \quad (10c) $$

$$ [f, g] := \int_{\hat{\omega}} f(p)g(p)dp, \quad p = \xi, \tau. \quad (10d) $$

This allows us to define some classical FE matrices below.
For instance, the $\ell$-th 2D quadrature point is defined by $\pi_\ell := (p_s, p_{s'}) \in \hat{c}$ where we have taken $\ell := \mathbf{m}(s, s')$. The case $K = 3$ is illustrated in figure 1.

<figure><img src="image_3.png" /><figcaption>Fig. 1 2D space-time quadrature points $\pi_\ell$ obtained as tensor products of the space and time 1D quadrature points $p_k$.</figcaption></figure>

This structure is particularly convenient as its extension from 1D to 2D or 3D problems in space is straightforward for structured grids. Indeed, the 2D or 3D basis functions in space can be defined and implemented as the products of the 1D basis functions with respect to the two or three space variables, the time variable being treated exactly in the same way as for the 1D case. This remark on structured grids stands not only for regular meshes but also for quadtree or octree grids, which are the target of this work.

## 2.3 Local scalar products and matrix definitions

Let us introduce the following notations corresponding to various scalar products defined by integrals

$$ \langle f, g \rangle := \int_{\hat{\omega}} \int_{\hat{T}} f(\xi, \tau)g(\xi, \tau)d\xi d\tau, \quad (10a) $$

$$ [f, g]_{\tau} := \int_{\hat{\omega}} f(\xi, \tau)g(\xi, \tau)d\xi, \quad (10b) $$

$$ \{f, g\}_{\xi} := \int_{\hat{T}} f(\xi, \tau)g(\xi, \tau)d\tau, \quad (10c) $$

$$ [f, g] := \int_{\hat{\omega}} f(p)g(p)dp, \quad p = \xi, \tau. \quad (10d) $$

This allows us to define some classical FE matrices below.
6 High order ADER-IPDG methods for advection-diffusion

• The space-time matrices belonging to the space M_L(ℝ)

$$
\underline{M}_{k, \ell} := \langle \Theta_k, \Theta_\ell \rangle, \qquad \underline{M}_{k, \ell}^{\tau} := [\Theta_k, \Theta_\ell]_{\tau}, \qquad (11a)
$$

$$
\underline{K}_{k, \ell}^{\tau} := \langle \partial_{\tau} \Theta_k, \Theta_{\ell} \rangle, \qquad (11b)
$$

• The mixed matrices belonging to the space M_{K,L}(ℝ)

$$
\underline{N}_{k, \ell}^{\tau} := [\theta_k, \Theta_{\ell}]_{\tau}, \qquad \underline{P}_{k, \ell} := \langle \theta_k, \Theta_{\ell} \rangle, \qquad (12a)
$$

$$
\underline{C}_{k, \ell} := \langle \partial_{\xi} \theta_k, \Theta_{\ell} \rangle, \qquad \underline{E}_{k, \ell} := \langle \partial_{\xi} \theta_k, \partial_{\xi} \Theta_{\ell} \rangle, \qquad (12b)
$$

$$
\underline{Q}_{k, \ell}^{p,q} := [\theta_k(p), \Theta_{\ell}(q, \cdot)], \qquad \underline{R}_{k, \ell}^{p,q} := [\theta_k(p), \partial_{\xi} \Theta_{\ell}(q, \cdot)], \qquad (12c)
$$

$$
\underline{W}_{k, \ell}^{p,q} := [\partial_{\xi} \theta_k(p), \Theta_{\ell}(q, \cdot)]. \qquad (12d)
$$

• The space matrix belonging to the space M_K(ℝ)

$$
\underline{N}_{k,l} := [\theta_k, \theta_l]. \qquad (13a)
$$

For the sake of clarity, for instance, the notation $\langle f, g \rangle$ denotes the space-time integral of the quantity $f(\xi, \tau)g(\xi, \tau)$ when both $f$ and $g$ depend on two variables, but it also holds for the quantity $f(\xi)g(\xi, \tau)$ if the function $f$ only depends on $\xi$ (we make the same comment on $g$).

The use of all these quantities will simplify the notations and help us to compactly write numerical schemes where all the contributions between neighbouring cells are easy to identify.

3 The ADER-IPDG scheme

3.1 The ADER-DG prediction-correction scheme

The ADER approach to build a high-order scheme in time is combined to a DG method in space where the solution is composed of piecewise polynomials per cell at each iteration in time. The ADER technique has been used for the resolution of a wide range of problems, especially in the framework of compressible flows [16–18, 21–23], even with a diffusion term [19, 24, 25].

Initially, the ADER discretization was performed by using a Cauchy-
Kowalewski procedure [26, 27], whether it was combined with a FV or a DG
scheme in space. In its current form, the ADER-DG prediction-correction
scheme is composed of two stages that we describe below in the particular case
of a linear unsteady advection-diffusion problem.

Predictor step

For any function $\hat{\phi} = \hat{q}, \hat{F}, \hat{S}$ defined on the reference space-time cell $\hat{c}$, let us introduce its polynomial expansions with respect to the space-time Lagrange basis functions (8). By linearity, the expansion of any derivative is
6 High order ADER-IPDG methods for advection-diffusion

• The space-time matrices belonging to the space M_L(ℝ)

$$
\underline{M}_{k, \ell} := \langle \Theta_k, \Theta_\ell \rangle, \qquad \underline{M}_{k, \ell}^{\tau} := [\Theta_k, \Theta_\ell]_{\tau}, \qquad (11a)
$$

$$
\underline{K}_{k, \ell}^{\tau} := \langle \partial_{\tau} \Theta_k, \Theta_{\ell} \rangle, \qquad (11b)
$$

• The mixed matrices belonging to the space M_{K,L}(ℝ)

$$
\underline{N}_{k, \ell}^{\tau} := [\theta_k, \Theta_{\ell}]_{\tau}, \qquad \underline{P}_{k, \ell} := \langle \theta_k, \Theta_{\ell} \rangle, \qquad (12a)
$$

$$
\underline{C}_{k, \ell} := \langle \partial_{\xi} \theta_k, \Theta_{\ell} \rangle, \qquad \underline{E}_{k, \ell} := \langle \partial_{\xi} \theta_k, \partial_{\xi} \Theta_{\ell} \rangle, \qquad (12b)
$$

$$
Q_{k, \ell}^{p,q} := [\theta_k(p), \Theta_{\ell}(q, \cdot)], \qquad \underline{R}_{k, \ell}^{p,q} := [\theta_k(p), \partial_{\xi} \Theta_{\ell}(q, \cdot)], \qquad (12c)
$$

$$
\underline{W}_{k, \ell}^{p,q} := [\partial_{\xi} \theta_k(p), \Theta_{\ell}(q, \cdot)]. \qquad (12d)
$$

• The space matrix belonging to the space M_K(ℝ)

$$
\underline{N}_{k,l} := [\theta_k, \theta_l]. \qquad (13a)
$$

For the sake of clarity, for instance, the notation $\langle f, g \rangle$ denotes the space-time integral of the quantity $f(\xi, \tau)g(\xi, \tau)$ when both $f$ and $g$ depend on two variables, but it also holds for the quantity $f(\xi)g(\xi, \tau)$ if the function $f$ only depends on $\xi$ (we make the same comment on $g$).

The use of all these quantities will simplify the notations and help us
to compactly write numerical schemes where all the contributions between
neighbouring cells are easy to identify.

3 The ADER-IPDG scheme

3.1 The ADER-DG prediction-correction scheme

The ADER approach to build a high-order scheme in time is combined to a
DG method in space where the solution is composed of piecewise polynomials
per cell at each iteration in time. The ADER technique has been used for
the resolution of a wide range of problems, especially in the framework of
compressible flows [16–18, 21–23], even with a diffusion term [19, 24, 25].

Initially, the ADER discretization was performed by using a Cauchy-
Kowalewski procedure [26, 27], whether it was combined with a FV or a DG
scheme in space. In its current form, the ADER-DG prediction-correction
scheme is composed of two stages that we describe below in the particular case
of a linear unsteady advection-diffusion problem.

Predictor step

For any function $\hat{\phi} = \hat{q}, \hat{F}, \hat{S}$ defined on the reference space-time cell $\hat{c}$,
let us introduce its polynomial expansions with respect to the space-time
Lagrange basis functions (8). By linearity, the expansion of any derivative is
straightforward

$$ \hat{\phi}(\xi, \tau) := \sum_{\ell=1}^{L} \hat{\phi}_{\ell} \Theta_{\ell}(\xi, \tau), \quad \partial_{p} \hat{\phi}(\xi, \tau) := \sum_{\ell=1}^{L} \hat{\phi}_{\ell} \partial_{p} \Theta_{\ell}(\xi, \tau), \quad p = \xi, \tau. \quad (14) $$

In the prediction stage, we multiply equation (6) by a test function $\Theta_k(\xi, \tau)$ and integrate by parts in time to determine the 2D local space-time Galerkin predictor $\hat{q}(\xi, \tau)$ starting from a 1D polynomial $\hat{u}^0(\xi, 0)$ at fictious time $\tau = 0$ (see for instance [22]). Using expansions (14), this procedure writes

$$ \langle \Theta_k, \partial_{\tau} \hat{q} \rangle + \left( \frac{\Delta t}{\Delta x} \right) \langle \Theta_k, \partial_{\xi} \hat{F}(\hat{q}) \rangle = \Delta t \langle \Theta_k, \hat{S} \rangle, \quad (15) $$

$$ \begin{aligned} & \Leftrightarrow \sum_{\ell=1}^{L} \left( \left[ \Theta_k, \Theta_\ell \right]_1 - \langle \partial_\tau \Theta_k, \Theta_\ell \rangle \right) \hat{q}_\ell + \left( \frac{\Delta t}{\Delta x} \right) \sum_{\ell=1}^{L} \langle \Theta_k, \partial_\xi \Theta_\ell \rangle \hat{F}_\ell \\ & = \sum_{\ell=1}^{K} \left[ \Theta_k, \theta_\ell \right]_0 \hat{u}_\ell^0 + \Delta t \sum_{\ell=1}^{L} \langle \Theta_k, \Theta_\ell \rangle \hat{S}_\ell, \end{aligned} \quad (16) $$

$$ \Leftrightarrow (\underline{M}^1 - \underline{K}^\tau) \hat{q} + \left(\frac{\Delta t}{\Delta x}\right) (\underline{K}^\xi)^T \hat{\mathbf{F}} = (\underline{N}^0)^T \hat{\mathbf{u}}^0 + \Delta t \underline{M} \underline{S}, \quad (17) $$

where we reused the notations (11) and $\mathbf{\Phi} := [\phi_1, \dots, \phi_*]^T \in \mathbb{R}^*$ for $\star = K, L$ is the vector containing the Degrees Of Freedom (DOF) of the 1D or 2D function $\phi$ with respect to the corresponding Lagrange basis. In particular, $\hat{\mathbf{u}}^0 \in \mathbb{R}^K$ contains the DOF of the numerical solution on a physical cell $\omega_i$ at time $t^n$.

Introducing the matrices and the vector

$$ \underline{A} := \underline{M}^1 - \underline{K}^\tau, \quad \underline{B} := \left(\frac{\Delta t}{\Delta x}\right) (\underline{K}^\xi)^T, \quad \hat{\mathbf{c}} := (\underline{N}^0)^T \hat{\mathbf{u}}^0 + \Delta t \underline{M} \hat{\mathbf{S}}, \quad (18) $$

the prediction step consists in solving the equation

$$ \underline{A}\hat{\mathbf{q}} + \underline{B}\hat{\mathbf{F}} = \hat{\mathbf{c}}. \quad (19) $$

When the physical flux is nonlinear, this system can be solved for $\hat{\mathbf{q}}$ with a Picard fixed-point algorithm for example. As the advective and diffusive fluxes given in (5) are linear with respect to the solution or its derivative, one can express the DOF of the total flux $\hat{F}(\hat{q})$ as a linear combination of the DOF of the predictor $\hat{q}$

$$ \hat{F} = a\hat{q} - \frac{\lambda}{\Delta x}\partial_{\xi}\hat{q} \Leftrightarrow \hat{\mathbf{F}} = \left(a\underline{I} - \frac{\lambda}{\Delta x}\underline{G}\right) \hat{\mathbf{q}} := \underline{H}\hat{\mathbf{q}}, \quad (20) $$

where the gradient matrix is defined by $\underline{G}_{k,\ell} := \partial_{\xi} \Theta_{\ell}(\pi_k)$ and $\underline{I} \in \mathcal{M}_L(\mathbb{R})$ stands for the identity matrix. Indeed, we exploit the fact that the physical
straightforward

$$ \hat{\phi}(\xi, \tau) := \sum_{\ell=1}^{L} \hat{\phi}_{\ell} \Theta_{\ell}(\xi, \tau), \quad \partial_{p} \hat{\phi}(\xi, \tau) := \sum_{\ell=1}^{L} \hat{\phi}_{\ell} \partial_{p} \Theta_{\ell}(\xi, \tau), \quad p = \xi, \tau. \quad (14) $$

In the prediction stage, we multiply equation (6) by a test function $\Theta_k(\xi, \tau)$ and integrate by parts in time to determine the 2D local space-time Galerkin predictor $\hat{q}(\xi, \tau)$ starting from a 1D polynomial $\hat{u}^0(\xi, 0)$ at fictious time $\tau = 0$ (see for instance [22]). Using expansions (14), this procedure writes

$$ \langle \Theta_k, \partial_{\tau} \hat{q} \rangle + \left( \frac{\Delta t}{\Delta x} \right) \langle \Theta_k, \partial_{\xi} \hat{F}(\hat{q}) \rangle = \Delta t \langle \Theta_k, \hat{S} \rangle, \quad (15) $$

$$ \begin{aligned} & \Leftrightarrow \sum_{\ell=1}^{L} \left( [\Theta_k, \Theta_\ell]_1 - \langle \partial_\tau \Theta_k, \Theta_\ell \rangle \right) \hat{q}_\ell + \left( \frac{\Delta t}{\Delta x} \right) \sum_{\ell=1}^{L} \langle \Theta_k, \partial_\xi \Theta_\ell \rangle \hat{F}_\ell \\ & = \sum_{\ell=1}^{K} [\Theta_k, \theta_\ell]_0 \hat{u}_\ell^0 + \Delta t \sum_{\ell=1}^{L} \langle \Theta_k, \Theta_\ell \rangle \hat{S}_\ell, \end{aligned} \quad (16) $$

$$ \Leftrightarrow (\underline{M}^1 - \underline{K}^\tau) \hat{q} + \left(\frac{\Delta t}{\Delta x}\right) (\underline{K}^\xi)^T \hat{\mathbf{F}} = (\underline{N}^0)^T \hat{\mathbf{u}}^0 + \Delta t \underline{M} \hat{\mathbf{S}}, \quad (17) $$

where we reused the notations (11) and $\mathbf{\Phi} := [\phi_1, \dots, \phi_*]^T \in \mathbb{R}^*$ for $\star = K, L$ is the vector containing the Degrees Of Freedom (DOF) of the 1D or 2D function $\phi$ with respect to the corresponding Lagrange basis. In particular, $\hat{\mathbf{u}}^0 \in \mathbb{R}^K$ contains the DOF of the numerical solution on a physical cell $\omega_i$ at time $t^n$.

Introducing the matrices and the vector

$$ \underline{A} := \underline{M}^1 - \underline{K}^\tau, \quad \underline{B} := \left(\frac{\Delta t}{\Delta x}\right) (\underline{K}^\xi)^T, \quad \hat{\mathbf{c}} := (\underline{N}^0)^T \hat{\mathbf{u}}^0 + \Delta t \underline{M} \hat{\mathbf{S}}, \quad (18) $$

the prediction step consists in solving the equation

$$ \underline{A}\hat{\mathbf{q}} + \underline{B}\hat{\mathbf{F}} = \hat{\mathbf{c}}. \quad (19) $$

When the physical flux is nonlinear, this system can be solved for $\hat{\mathbf{q}}$ with a Picard fixed-point algorithm for example. As the advective and diffusive fluxes given in (5) are linear with respect to the solution or its derivative, one can express the DOF of the total flux $\hat{F}(\hat{q})$ as a linear combination of the DOF of the predictor $\hat{q}$

$$ \hat{F} = a\hat{q} - \frac{\lambda}{\Delta x}\partial_{\xi}\hat{q} \Leftrightarrow \hat{\mathbf{F}} = \left(a\underline{I} - \frac{\lambda}{\Delta x}\underline{G}\right) \hat{\mathbf{q}} := \underline{H}\hat{\mathbf{q}}, \quad (20) $$

where the gradient matrix is defined by $\underline{G}_{k,\ell} := \partial_{\xi} \Theta_{\ell}(\pi_k)$ and $\underline{I} \in \mathcal{M}_L(\mathbb{R})$ stands for the identity matrix. Indeed, we exploit the fact that the physical
8 *High order ADER-IPDG methods for advection-diffusion*

flux $\hat{F}$ is linear, so we use the following expansion for its evaluation at the k-th space-time quadrature point (i.e. the k-th DOF of the vector $\hat{\mathbf{F}}$)

$$ \hat{F}(\hat{q}(\pi_k)) = \sum_{\ell=1}^{L} \left( a\Theta_{\ell}(\pi_k) - \frac{\lambda}{\Delta x} \partial_{\xi}\Theta_{\ell}(\pi_k) \right) \hat{q}_{\ell}, \quad (21) $$

as presented in eq. (14). In particular, we observe that the matrix $\underline{B}$ represents the divergence operator (through matrix $\underline{K}^{\xi}$, which is computed with quadratures), while the matrix $\underline{G}$ represents the gradient operator (through a pointwise evaluation). This implies that the second order derivatives of the basis functions never need to be computed.

Finally, the space-time Galerkin predictor is completely determined by solving the linear system of size $L$

$$ (\underline{A} + \underline{B}\underline{H}) \hat{\mathbf{q}} = \hat{\mathbf{c}}. \quad (22) $$

The resolution of this small local linear system is not expensive, especially if one uses a uniform mesh and constant values for $a$ and $\lambda$. Indeed, the matrix $(\underline{A} + \underline{B}\underline{H})$ does not depend on the space-time cell $c_i^n$ on which the prediction step is solved, so it can be pre-computed once to lower the total computational cost. We also remark that this prediction stage is embarrassingly parallel because no communications are needed between the discretization cells.

### Corrector step

Within the DG framework, the solution is represented by piecewise polynomials of arbitrary degree ($K-1$) in space, i.e., each cell contains $K$ local DOF relative to $\hat{u}$, defined on the reference space-time cell $\hat{c}$. Let us introduce its polynomial expansion with respect to the space Lagrange basis functions (7)

$$ \hat{u}(\xi, \tau) := \sum_{\ell=1}^{K} \hat{u}_{\ell}(\tau)\theta_{\ell}(\xi), \quad (23) $$

the DOF $\hat{u}_{\ell}(\tau)$ being updated at each iteration.

In the correction stage, we multiply equation (6) by a space test function $\hat{\psi} \equiv \theta_k(\xi)$ and integrate by parts in both time and space to couple the DOF between neighbouring cells, and we replace the solution $\hat{u}$ by the predictor $\hat{q}$ in the fluxes. In the following, we denote by $\tilde{\phi}$ the numerical flux of the quantity $\hat{\phi}$ standing for the physical fluxes $\hat{F}(\hat{q})$ or the solution itself $\hat{q}$. This step then reads

$$ \begin{aligned} & \langle \hat{\psi}, \partial_{\tau} \hat{u} \rangle + \left( \frac{\Delta t}{\Delta x} \right) \langle \hat{\psi}, \partial_{\xi} \hat{F}(\hat{q}) \rangle = \Delta t \langle \hat{\psi}, \hat{S} \rangle, \\ \Leftrightarrow & [\hat{\psi}, \hat{u}(\cdot, 1)] = [\hat{\psi}, \hat{u}(\cdot, 0)] + \left( \frac{\Delta t}{\Delta x} \right) \langle \partial_{\xi} \hat{\psi}, \hat{F}(\hat{q}) \rangle + \Delta t \langle \hat{\psi}, \hat{S} \rangle \end{aligned} $$
8 High order ADER-IPDG methods for advection-diffusion

flux $\hat{F}$ is linear, so we use the following expansion for its evaluation at the k-th space-time quadrature point (i.e. the k-th DOF of the vector $\hat{\mathbf{F}}$)

$$ \hat{F}(\hat{q}(\pi_k)) = \sum_{\ell=1}^{L} \left( a \Theta_{\ell}(\pi_k) - \frac{\lambda}{\Delta x} \partial_{\xi} \Theta_{\ell}(\pi_k) \right) \hat{q}_{\ell}, \qquad (21) $$

as presented in eq. (14). In particular, we observe that the matrix $\underline{B}$ represents the divergence operator (through matrix $\underline{K}^{\xi}$, which is computed with quadratures), while the matrix $\underline{G}$ represents the gradient operator (through a pointwise evaluation). This implies that the second order derivatives of the basis functions never need to be computed.

Finally, the space-time Galerkin predictor is completely determined by solving the linear system of size $L$

$$ (\underline{A} + \underline{B}\underline{H}) \hat{\mathbf{q}} = \hat{\mathbf{c}}. \qquad (22) $$

The resolution of this small local linear system is not expensive, especially if one uses a uniform mesh and constant values for $a$ and $\lambda$. Indeed, the matrix $(\underline{A} + \underline{B}\underline{H})$ does not depend on the space-time cell $c_i^n$ on which the prediction step is solved, so it can be pre-computed once to lower the total computational cost. We also remark that this prediction stage is embarrassingly parallel because no communications are needed between the discretization cells.

### Corrector step

Within the DG framework, the solution is represented by piecewise polynomials of arbitrary degree ($K-1$) in space, i.e., each cell contains $K$ local DOF relative to $\hat{u}$, defined on the reference space-time cell $\hat{c}$. Let us introduce its polynomial expansion with respect to the space Lagrange basis functions (7)

$$ \hat{u}(\xi, \tau) := \sum_{\ell=1}^{K} \hat{u}_{\ell}(\tau) \theta_{\ell}(\xi), \qquad (23) $$

the DOF $\hat{u}_{\ell}(\tau)$ being updated at each iteration.

In the correction stage, we multiply equation (6) by a space test function $\hat{\psi} \equiv \theta_k(\xi)$ and integrate by parts in both time and space to couple the DOF between neighbouring cells, and we replace the solution $\hat{u}$ by the predictor $\hat{q}$ in the fluxes. In the following, we denote by $\tilde{\phi}$ the numerical flux of the quantity $\hat{\phi}$ standing for the physical fluxes $\hat{F}(\hat{q})$ or the solution itself $\hat{q}$. This step then reads

$$ \begin{aligned} & \langle \hat{\psi}, \partial_{\tau} \hat{u} \rangle + \left( \frac{\Delta t}{\Delta x} \right) \langle \hat{\psi}, \partial_{\xi} \hat{F}(\hat{q}) \rangle = \Delta t \langle \hat{\psi}, \hat{S} \rangle, \\ & \Leftrightarrow [\hat{\psi}, \hat{u}(\cdot, 1)] = [\hat{\psi}, \hat{u}(\cdot, 0)] + \left( \frac{\Delta t}{\Delta x} \right) \langle \partial_{\xi} \hat{\psi}, \hat{F}(\hat{q}) \rangle + \Delta t \langle \hat{\psi}, \hat{S} \rangle \end{aligned} $$
$$ - \left( \frac{\Delta t}{\Delta x} \right) \left( \{\hat{\psi}, \tilde{F}(\tilde{q})\}_1 - \{\hat{\psi}, \tilde{F}(\tilde{q})\}_0 \right), \quad (24) $$

where we keep in mind that $\hat{F} = \hat{F}_a + \hat{F}_d$. In the following, the terms $\langle \partial_{\xi} \hat{\psi}, \hat{F}(\hat{q}) \rangle$ and $\{\hat{\psi}, \tilde{F}(\tilde{q})\}_{\xi}$ will be refered to as volume and interface terms respectively. The three different possibilities to solve this correction stage in this paper are developed in the next section dedicated to IPDG methods applied to time-dependent problems.

In figure 2, we summarize the ADER procedure. The prediction stage consists in determining a local continuous space-time Galerkin predictor from the knowledge of the initial data, which is a polynomial of high degree in the DG approach. The correction stage recombines the information between neighbouring cells through the computation of numerical fluxes.

<figure><img src="image_2.png" /><figcaption>Fig. 2 Top: representation of the solution by discontinuous piecewise polynomials. Bottom: illustration of the ADER-DG prediction-correction scheme.</figcaption></figure>

The choice of Lagrange polynomials passing through the Gauss-Legendre quadrature points as nodal basis functions allows to reduce the global computational cost. Indeed, all the integrals (10) are computed by directly using the values of the unknowns at the Gauss-Legendre nodes at no additional cost.

## 3.2 A family of IPDG methods

In the previous variational formulation (24), let us break the total flux into its advective and diffusive components and replace them by their expressions in the volume terms. Since the solution may be discontinuous across cell interfaces, there is not a unique definition for the solution and its derivative at a
$$ - \left( \frac{\Delta t}{\Delta x} \right) \left( \{\hat{\psi}, \tilde{F}(\tilde{q})\}_1 - \{\hat{\psi}, \tilde{F}(\tilde{q})\}_0 \right), \quad (24) $$

where we keep in mind that $\tilde{F} = \tilde{F}_a + \tilde{F}_d$. In the following, the terms $\langle \partial_{\xi} \hat{\psi}, \tilde{F}(\tilde{q}) \rangle$ and $\{\hat{\psi}, \tilde{F}(\tilde{q})\}_{\xi}$ will be refered to as volume and interface terms respectively. The three different possibilities to solve this correction stage in this paper are developed in the next section dedicated to IPDG methods applied to time-dependent problems.

In figure 2, we summarize the ADER procedure. The prediction stage consists in determining a local continuous space-time Galerkin predictor from the knowledge of the initial data, which is a polynomial of high degree in the DG approach. The correction stage recombines the information between neighbouring cells through the computation of numerical fluxes.

<figure><img src="image_4.png" /><figcaption>Fig. 2 Top: representation of the solution by discontinuous piecewise polynomials. Bottom: illustration of the ADER-DG prediction-correction scheme.</figcaption></figure>

The choice of Lagrange polynomials passing through the Gauss-Legendre quadrature points as nodal basis functions allows to reduce the global computational cost. Indeed, all the integrals (10) are computed by directly using the values of the unknowns at the Gauss-Legendre nodes at no additional cost.

## 3.2 A family of IPDG methods

In the previous variational formulation (24), let us break the total flux into its advective and diffusive components and replace them by their expressions in the volume terms. Since the solution may be discontinuous across cell interfaces, there is not a unique definition for the solution and its derivative at a
10 *High order ADER-IPDG methods for advection-diffusion*

given interface between two neighbouring cells. Thus, we add the subscript *i* to the quantities that refer to the cell $\omega_i$ and are uniquely defined. Now, for all interfaces, we introduce the discrete parameter $\varepsilon \in \{-1, 0, +1\}$ and consider the modified variational formulation below

$$ 
\begin{aligned}
[\hat{\psi}_i, \hat{u}_i(\cdot, 1)] = & [\hat{\psi}_i, \hat{u}_i(\cdot, 0)] + \Delta t \langle \hat{\psi}_i, \hat{S}_i \rangle \\
& + \left(\frac{\Delta t}{\Delta x}\right) \langle \partial_{\xi} \hat{\psi}_i, a\hat{q}_i \rangle - \left(\frac{\Delta t}{\Delta x^2}\right) \langle \partial_{\xi} \hat{\psi}_i, \lambda \partial_{\xi} \hat{q}_i \rangle \\
& - \left(\frac{\Delta t}{\Delta x}\right) \left( \left\{ \hat{\psi}_i, \tilde{F}_a(\tilde{q}) \right\}_1 - \left\{ \hat{\psi}_i, \tilde{F}_a(\tilde{q}) \right\}_0 \right) \\
& - \left(\frac{\Delta t}{\Delta x}\right) \left( \left\{ \hat{\psi}_i, \tilde{F}_d(\tilde{q}) \right\}_1 - \left\{ \hat{\psi}_i, \tilde{F}_d(\tilde{q}) \right\}_0 \right) \\
& + \varepsilon \left(\frac{\Delta t}{\Delta x}\right) \left( \left\{ \tilde{q}, \hat{F}_d(\hat{\psi}_i) \right\}_1 - \left\{ \tilde{q}, \hat{F}_d(\hat{\psi}_i) \right\}_0 \right).
\end{aligned} 
\quad (25) $$

Clearly, the case $\varepsilon = 0$ corresponds to the previous equation (24); the case $\varepsilon = +1$ creates an anti-symmetry with respect to $\tilde{q}$ and $\hat{\psi}$ on the diffusive interface terms; finally, the case $\varepsilon = -1$ provides a total symmetry on these diffusive interface terms. The interface terms involving the advective flux are already symmetric by construction. This manipulation leads to the definition of the three different IPDG methods studied in this paper, as shown in table 1.

<table><caption>**Table 1** Three IPDG methods depending on the value of the parameter $\varepsilon$.</caption><thead><tr><th>Parameter &epsilon;</th><th>IPDG method</th><th>Abbreviation</th></tr></thead><tbody><tr><td>-1</td><td>Symmetric Interior Penalty Galerkin</td><td>SIPG</td></tr><tr><td>0</td><td>Incomplete Interior Penalty Galerkin</td><td>IIPG</td></tr><tr><td>+1</td><td>Nonsymmetric Interior Penalty Galerkin</td><td>NIPG</td></tr></tbody></table>

Many details on these three IPDG methods can be found in [13] for theoretical or implementation aspects, with Euler or RK schemes for the time integration. A presentation of some classical IP techniques can be found in [5, 14]. We draw the reader's attention to the fact that the ADER-IIPG scheme actually corresponds to the classical ADER-DG scheme used in [17, 18, 22] for example.

Independently of the IPDG method used, the quantities $\hat{q}$ and $\hat{F}(\hat{q})$ need to be replaced by a consistent and conservative numerical flux which is detailed in the next section.
10 *High order ADER-IPDG methods for advection-diffusion*

given interface between two neighbouring cells. Thus, we add the subscript *i* to the quantities that refer to the cell $\omega_i$ and are uniquely defined. Now, for all interfaces, we introduce the discrete parameter $\varepsilon \in \{-1, 0, +1\}$ and consider the modified variational formulation below

$$ 
\begin{aligned}
[\hat{\psi}_i, \hat{u}_i(\cdot, 1)] = [\hat{\psi}_i, \hat{u}_i(\cdot, 0)] &+ \Delta t \langle \hat{\psi}_i, \hat{S}_i \rangle \\
& + \left(\frac{\Delta t}{\Delta x}\right) \langle \partial_{\xi} \hat{\psi}_i, a\hat{q}_i \rangle - \left(\frac{\Delta t}{\Delta x^2}\right) \langle \partial_{\xi} \hat{\psi}_i, \lambda\partial_{\xi} \hat{q}_i \rangle \\
& - \left(\frac{\Delta t}{\Delta x}\right) \left(\left\{\hat{\psi}_i, \tilde{F}_a(\tilde{q})\right\}_1 - \left\{\hat{\psi}_i, \tilde{F}_a(\tilde{q})\right\}_0\right) \\
& - \left(\frac{\Delta t}{\Delta x}\right) \left(\left\{\hat{\psi}_i, \tilde{F}_d(\tilde{q})\right\}_1 - \left\{\hat{\psi}_i, \tilde{F}_d(\tilde{q})\right\}_0\right) \\
& + \varepsilon \left(\frac{\Delta t}{\Delta x}\right) \left(\left\{\tilde{q}, \hat{F}_d(\hat{\psi}_i)\right\}_1 - \left\{\tilde{q}, \hat{F}_d(\hat{\psi}_i)\right\}_0\right).
\end{aligned} 
\quad (25) $$

Clearly, the case $\varepsilon = 0$ corresponds to the previous equation (24); the case $\varepsilon = +1$ creates an anti-symmetry with respect to $\tilde{q}$ and $\hat{\psi}$ on the diffusive interface terms; finally, the case $\varepsilon = -1$ provides a total symmetry on these diffusive interface terms. The interface terms involving the advective flux are already symmetric by construction. This manipulation leads to the definition of the three different IPDG methods studied in this paper, as shown in table 1.

<table><caption>**Table 1** Three IPDG methods depending on the value of the parameter $\varepsilon$.</caption><thead><tr><td>Parameter ε</td><td>IPDG method</td><td>Abbreviation</td></tr></thead><tbody><tr><td>-1</td><td>Symmetric Interior Penalty Galerkin</td><td>SIPG</td></tr><tr><td>0</td><td>Incomplete Interior Penalty Galerkin</td><td>IIPG</td></tr><tr><td>+1</td><td>Nonsymmetric Interior Penalty Galerkin</td><td>NIPG</td></tr></tbody></table>

Many details on these three IPDG methods can be found in [13] for theoretical or implementation aspects, with Euler or RK schemes for the time integration. A presentation of some classical IP techniques can be found in [5, 14]. We draw the reader's attention to the fact that the ADER-IIPG scheme actually corresponds to the classical ADER-DG scheme used in [17, 18, 22] for example.

Independently of the IPDG method used, the quantities $\hat{q}$ and $\hat{F}(\hat{q})$ need to be replaced by a consistent and conservative numerical flux which is detailed in the next section.
## 3.3 Cell average, jump and numerical flux

At the interface $x_i$ separating the two cells $\omega_i$ and $\omega_{i+1}$, we define the limits of any function $\phi$ from the left and the right sides at a given time $t$ as

$$ \phi_i^-(t) := \lim_{\delta \to 0^+} \phi(x_i - \delta, t), \quad \phi_i^+(t) := \lim_{\delta \to 0^+} \phi(x_i + \delta, t). \qquad (26) $$

Thus the average and the jump of the quantity $\phi$ across an interior interface $x_i$ for $i = 1, \dots, N_x - 1$ read respectively

$$ \langle \phi(t) \rangle_i := \frac{1}{2} [\phi_i^-(t) + \phi_i^+(t)], \qquad (27) $$

$$ \llbracket \phi(t) \rrbracket_i := \phi_i^-(t) - \phi_i^+(t). \qquad (28) $$

The numerical fluxes are generally designed on the same basis which is a linear combination of the average of the physical flux and the jump of the solution, such as the Lax-Friedrichs or the Rusanov numerical fluxes. In this work, we use a Rusanov-like numerical flux both for the advective and the diffusive components at an interface $e := \omega_i \cap \omega_{i \pm 1}$.

* For interior interfaces

$$ \tilde{F}_a(\tilde{q}_e(\tau)) := \langle \hat{F}_a(\hat{q}(\tau)) \rangle_i + \sigma_a \llbracket \hat{q}(\tau) \rrbracket_i, \qquad (29) $$

$$ \tilde{F}_d(\tilde{q}_e(\tau)) := \langle \hat{F}_d(\hat{q}(\tau)) \rangle_i + \sigma_d \llbracket \hat{q}(\tau) \rrbracket_i, \qquad (30) $$

$$ \tilde{q}_e(\tau) := \frac{1}{2} \llbracket \hat{q}(\tau) \rrbracket_i \vec{n}. \qquad (31) $$

* For Dirichlet boundary interfaces

$$ \tilde{F}_a(\tilde{q}_e(\tau)) := \frac{1}{2} [\hat{F}(\hat{q}^{\mp}(\tau)) + \hat{F}(g_{\pm}(\tau))] + \sigma_a (\hat{q}^{\mp}(\tau) - g_{\pm}(\tau)) \vec{n}, \qquad (32) $$

$$ \tilde{F}_d(\tilde{q}_e(\tau)) := \hat{F}_d(\hat{q}^{\mp}(\tau)) + \sigma_d (\hat{q}^{\mp}(\tau) - g_{\pm}(\tau)) \vec{n}, \qquad (33) $$

$$ \tilde{q}_e(\tau) := \hat{q}^{\mp}(\tau) - g_{\pm}(\tau). \qquad (34) $$

When we consider the cell $\omega_i$, the outward unit vector $\vec{n}$ in 1D comes down to $-1$ (respectively $+1$) on the left (respectively right) interface of the cell. Moreover, recalling that $(K - 1)$ is the degree of the polynomials used in the ADER-IPDG prediction-correction scheme, the coefficients $\sigma_a, \sigma_d > 0$ are defined as

$$ \sigma_a := \frac{|a|}{2}, \quad \sigma_d := \mu \frac{\lambda}{2\Delta x}, \quad \mu := \gamma(K - 1)^2, \qquad (35) $$

where $\gamma > 0$ is a penalization parameter for the diffusive flux (see [14] for discussions on the choice for $\mu$). Finally, replacing the numerical fluxes $\tilde{q}$ and $\tilde{F}$ in variational formulation (25) by their expressions (29) and (32) allows to rewrite explicitly the complete scheme according to the different unknowns as well as the boundary conditions.
## 3.3 Cell average, jump and numerical flux

At the interface $x_i$ separating the two cells $\omega_i$ and $\omega_{i+1}$, we define the limits of any function $\phi$ from the left and the right sides at a given time $t$ as

$$ \phi_i^-(t) := \lim_{\delta \to 0^+} \phi(x_i - \delta, t), \quad \phi_i^+(t) := \lim_{\delta \to 0^+} \phi(x_i + \delta, t). \quad (26) $$

Thus the average and the jump of the quantity $\phi$ across an interior interface $x_i$ for $i = 1, \dots, N_x - 1$ read respectively

$$ \langle \phi(t) \rangle_i := \frac{1}{2} [\phi_i^-(t) + \phi_i^+(t)], \quad (27) $$

$$ \llbracket \phi(t) \rrbracket_i := \phi_i^-(t) - \phi_i^+(t). \quad (28) $$

The numerical fluxes are generally designed on the same basis which is a linear combination of the average of the physical flux and the jump of the solution, such as the Lax-Friedrichs or the Rusanov numerical fluxes. In this work, we use a Rusanov-like numerical flux both for the advective and the diffusive components at an interface $e := \omega_i \cap \omega_{i \pm 1}$.

* For interior interfaces

$$ \tilde{F}_a(\tilde{q}_e(\tau)) := \langle \hat{F}_a(\hat{q}(\tau)) \rangle_i + \sigma_a \llbracket \hat{q}(\tau) \rrbracket_i, \quad (29) $$

$$ \tilde{F}_d(\tilde{q}_e(\tau)) := \langle \hat{F}_d(\hat{q}(\tau)) \rangle_i + \sigma_d \llbracket \hat{q}(\tau) \rrbracket_i, \quad (30) $$

$$ \tilde{q}_e(\tau) := \frac{1}{2} \llbracket \hat{q}(\tau) \rrbracket_i \vec{n}. \quad (31) $$

* For Dirichlet boundary interfaces

$$ \tilde{F}_a(\tilde{q}_e(\tau)) := \frac{1}{2} [\hat{F}(\hat{q}^{\mp}(\tau)) + \hat{F}(g_{\pm}(\tau))] + \sigma_a (\hat{q}^{\mp}(\tau) - g_{\pm}(\tau)) \vec{n}, \quad (32) $$

$$ \tilde{F}_d(\tilde{q}_e(\tau)) := \hat{F}_d(\hat{q}^{\mp}(\tau)) + \sigma_d (\hat{q}^{\mp}(\tau) - g_{\pm}(\tau)) \vec{n}, \quad (33) $$

$$ \tilde{q}_e(\tau) := \hat{q}^{\mp}(\tau) - g_{\pm}(\tau). \quad (34) $$

When we consider the cell $\omega_i$, the outward unit vector $\vec{n}$ in 1D comes down to $-1$ (respectively $+1$) on the left (respectively right) interface of the cell. Moreover, recalling that $(K - 1)$ is the degree of the polynomials used in the ADER-IPDG prediction-correction scheme, the coefficients $\sigma_a, \sigma_d > 0$ are defined as

$$ \sigma_a := \frac{|a|}{2}, \quad \sigma_d := \mu \frac{\lambda}{2\Delta x}, \quad \mu := \gamma(K - 1)^2, \quad (35) $$

where $\gamma > 0$ is a penalization parameter for the diffusive flux (see [14] for discussions on the choice for $\mu$). Finally, replacing the numerical fluxes $\tilde{q}$ and $\tilde{F}$ in variational formulation (25) by their expressions (29) and (32) allows to rewrite explicitly the complete scheme according to the different unknowns as well as the boundary conditions.
12 High order ADER-IPDG methods for advection-diffusion

This choice of numerical fluxes for the advective and diffusive components comes with the definition of a suitable CFL condition. A synthesis of the approaches proposed in [19–21, 24] leads us to compute the time step using the following CFL restriction

$$ \Delta t := \alpha(k) \frac{\Delta x_{\star}^2}{2\lambda + |a|\Delta x_{\star}}, \quad \Delta x_{\star} := \frac{\Delta x}{2k+1}. \qquad (36) $$

Contrary to FV schemes, the time step for DG schemes strongly depends on the degree $k$ of the approximation polynomials. It is possible to apply an additional stabilization through the parameter $\alpha(k)$ which may depend on $k$, as explored in [5] and reported in [26].

## 3.4 Construction of the local linear systems

In the particular case of the advection-diffusion equation, the physical flux is linear with respect to the unknown and its derivative. Consequently, the numerical scheme can be cast into the form of a linear combination of the unknown DOF on each cell of the mesh. Starting from the variational formulation (25), the case of interior or boundary cells needs to be treated separately (see appendix A for the computations). To simplify the final expression of the numerical schemes, we introduce the notations

$$ \underline{Z}_{\star} := \left( \frac{1}{\Delta x} \right) \underline{Z}_{\star}^{a} - \left( \frac{1}{\Delta x^2} \right) \underline{Z}_{\star}^{d}, \qquad \star = L, C, R, \qquad (37) $$

$$ \mathbf{b}_{\star} := \left( \frac{1}{\Delta x} \right) \mathbf{b}_{\star}^{a} - \left( \frac{1}{\Delta x^2} \right) \mathbf{b}_{\star}^{d}, \qquad \star = L, R, \qquad (38) $$

the former quantity standing for matrices and the latter for vectors. The letters $L, C, R$ refer respectively to the left, center and right positions.

Interior cells. The ADER-IPDG scheme for interior cells reads

$$ \underline{N} \mathbf{u}_{i}^{n+1} = \underline{N} \mathbf{u}_{i}^{n} + \Delta t (\underline{Z}_{L} \mathbf{q}_{i-1}^{n} + \underline{Z}_{C} \mathbf{q}_{i}^{n} + \underline{Z}_{R} \mathbf{q}_{i+1}^{n}) + \Delta t \underline{P} \mathbf{S}_{i}^{n}, \qquad (39) $$

where the advective and diffusive local space-time matrices are defined as

$$ \underline{Z}_{L}^{a} := a^{+} \underline{Q}^{0,1}, \qquad (40) $$

$$ \underline{Z}_{C}^{a} := a \underline{C} - a^{+} \underline{Q}^{1,1} - a^{-} \underline{Q}^{0,0}, \qquad (41) $$

$$ \underline{Z}_{R}^{a} := a^{-} \underline{Q}^{1,0}, \qquad (42) $$

$$ \underline{Z}_{L}^{d} := \frac{\lambda}{2} [\underline{R}^{0,1} + \varepsilon \underline{W}^{0,1} - \mu \underline{Q}^{0,1}], \qquad (43) $$

$$ \underline{Z}_{C}^{d} := \frac{\lambda}{2} [2\underline{E} + \underline{R}^{0,0} - \underline{R}^{1,1} - \varepsilon (\underline{W}^{0,0} - \underline{W}^{1,1}) + \mu (\underline{Q}^{0,0} + \underline{Q}^{1,1})], \qquad (44) $$

$$ \underline{Z}_{R}^{d} := \frac{\lambda}{2} [-\underline{R}^{1,0} - \varepsilon \underline{W}^{1,0} - \mu \underline{Q}^{1,0}]. \qquad (45) $$
12 High order ADER-IPDG methods for advection-diffusion

This choice of numerical fluxes for the advective and diffusive components comes with the definition of a suitable CFL condition. A synthesis of the approaches proposed in [19–21, 24] leads us to compute the time step using the following CFL restriction

$$ \Delta t := \alpha(k) \frac{\Delta x_{\star}^2}{2\lambda + |a|\Delta x_{\star}}, \quad \Delta x_{\star} := \frac{\Delta x}{2k+1}. \qquad (36) $$

Contrary to FV schemes, the time step for DG schemes strongly depends on the degree $k$ of the approximation polynomials. It is possible to apply an additional stabilization through the parameter $\alpha(k)$ which may depend on $k$, as explored in [5] and reported in [26].

## 3.4 Construction of the local linear systems

In the particular case of the advection-diffusion equation, the physical flux is linear with respect to the unknown and its derivative. Consequently, the numerical scheme can be cast into the form of a linear combination of the unknown DOF on each cell of the mesh. Starting from the variational formulation (25), the case of interior or boundary cells needs to be treated separately (see appendix A for the computations). To simplify the final expression of the numerical schemes, we introduce the notations

$$ \underline{Z}_{\star} := \left( \frac{1}{\Delta x} \right) \underline{Z}_{\star}^{a} - \left( \frac{1}{\Delta x^2} \right) \underline{Z}_{\star}^{d}, \qquad \star = L, C, R, \qquad (37) $$

$$ \mathbf{b}_{\star} := \left( \frac{1}{\Delta x} \right) \mathbf{b}_{\star}^{a} - \left( \frac{1}{\Delta x^2} \right) \mathbf{b}_{\star}^{d}, \qquad \star = L, R, \qquad (38) $$

the former quantity standing for matrices and the latter for vectors. The letters $L, C, R$ refer respectively to the left, center and right positions.

Interior cells. The ADER-IPDG scheme for interior cells reads

$$ \underline{N} \mathbf{u}_{i}^{n+1} = \underline{N} \mathbf{u}_{i}^{n} + \Delta t (\underline{Z}_{L} \mathbf{q}_{i-1}^{n} + \underline{Z}_{C} \mathbf{q}_{i}^{n} + \underline{Z}_{R} \mathbf{q}_{i+1}^{n}) + \Delta t \underline{P} \mathbf{S}_{i}^{n}, \qquad (39) $$

where the advective and diffusive local space-time matrices are defined as

$$ \underline{Z}_{L}^{a} := a^{+} \underline{Q}^{0,1}, \qquad (40) $$

$$ \underline{Z}_{C}^{a} := a \underline{C} - a^{+} \underline{Q}^{1,1} - a^{-} \underline{Q}^{0,0}, \qquad (41) $$

$$ \underline{Z}_{R}^{a} := a^{-} \underline{Q}^{1,0}, \qquad (42) $$

$$ \underline{Z}_{L}^{d} := \frac{\lambda}{2} [\underline{R}^{0,1} + \varepsilon \underline{W}^{0,1} - \mu \underline{Q}^{0,1}], \qquad (43) $$

$$ \underline{Z}_{C}^{d} := \frac{\lambda}{2} [2\underline{E} + \underline{R}^{0,0} - \underline{R}^{1,1} - \varepsilon (\underline{W}^{0,0} - \underline{W}^{1,1}) + \mu (\underline{Q}^{0,0} + \underline{Q}^{1,1})], \qquad (44) $$

$$ \underline{Z}_{R}^{d} := \frac{\lambda}{2} [-\underline{R}^{1,0} - \varepsilon \underline{W}^{1,0} - \mu \underline{Q}^{1,0}]. \qquad (45) $$
*Left boundary cell.* The ADER-IPDG scheme for the left cell reads

$$ \underline{N} \mathbf{u}_1^{n+1} = \underline{N} \mathbf{u}_1^n + \Delta t (\underline{Z}_C \mathbf{q}_1^n + \underline{Z}_R \mathbf{q}_2^n) + \Delta t \underline{P} \mathbf{S}_1^n + \Delta t \mathbf{b}_L^n, \quad (46) $$

where the local space-time matrices are exactly the same as for interior cells, except for $\underline{Z}_C^d$ in which a coefficient 2 has to be added in front of the matrices $\underline{R}^{0,0}$ and $\underline{W}^{0,0}$. The boundary conditions for the advective and diffusive components are

$$ (\mathbf{b}_L^a)^n := a^+ \mathbf{d}^0, \quad (\mathbf{b}_L^d)^n := \frac{\lambda}{2} (-\mu \mathbf{d}^0 + 2\varepsilon \mathbf{e}^0). \quad (47) $$

*Right boundary cell.* The ADER-IPDG scheme for the right cell reads

$$ \underline{N} \mathbf{u}_{N_x}^{n+1} = \underline{N} \mathbf{u}_{N_x}^n + \Delta t (\underline{Z}_L \mathbf{q}_{N_x-1}^n + \underline{Z}_C \mathbf{q}_{N_x}^n) + \Delta t \underline{P} \mathbf{S}_{N_x}^n + \Delta t \mathbf{b}_R^n, \quad (48) $$

where, again, all the local matrices are the same as for interior cells, except for $\underline{Z}_C^d$ in which a coefficient 2 needs to be added in front of the matrices $\underline{R}^{1,1}$ and $\underline{W}^{1,1}$. The advective and diffusive components for the boundary conditions are

$$ (\mathbf{b}_R^a)^n := a^- \mathbf{d}^1, \quad (\mathbf{b}_R^d)^n := -\frac{\lambda}{2} (\mu \mathbf{d}^1 + 2\varepsilon \mathbf{e}^1). \quad (49) $$

The remark on the computational cost that we made on the prediction step also holds for this correction stage, regardless the IPDG method chosen. Indeed, as we consider that $a$ and $\lambda$ are constants and we use a uniform mesh, all the matrices involved are constant and can be pre-computed only once. Moreover, this linear system to be solved for $\mathbf{u}_i^n$ is of small size $K$ and the matrix $\underline{N}$ is diagonal, which is a consequence of our choice of an orthogonal Lagrange polynomial basis. On the other hand, and contrary to FV methods that need a large stencil to reconstruct the information with a high precision, this DG scheme is compact (it is enough to know the DOF of the cell $\omega_i$ and of its two immediate neighbours) and thus the parallelization is simple and efficient.

# 4 Numerical validation of the high-order ADER-IPDG scheme

In this section, we present some numerical results on two different test cases in one dimension. The first one is a classical test: a gaussian that remains far from the boundaries. The second one is a bit more challenging: a smooth step initially located out of the domain and which crosses the boundary. Finally, for the third example, we propose a two dimensional extension.

## 4.1 First example in 1D

In this section, we present a manufactured test case. We take an advection velocity $a = 1$ and a diffusion coefficient $\lambda = 1$ on a computational domain
Left boundary cell. The ADER-IPDG scheme for the left cell reads

$$ \underline{N} \mathbf{u}_1^{n+1} = \underline{N} \mathbf{u}_1^n + \Delta t (\underline{Z}_C \mathbf{q}_1^n + \underline{Z}_R \mathbf{q}_2^n) + \Delta t \underline{P} \mathbf{S}_1^n + \Delta t \mathbf{b}_L^n, \quad (46) $$

where the local space-time matrices are exactly the same as for interior cells, except for $\underline{Z}_C^d$ in which a coefficient 2 has to be added in front of the matrices $\underline{R}^{0,0}$ and $\underline{W}^{0,0}$. The boundary conditions for the advective and diffusive components are

$$ (\mathbf{b}_L^a)^n := a^+ \mathbf{d}^0, \quad (\mathbf{b}_L^d)^n := \frac{\lambda}{2} (-\mu \mathbf{d}^0 + 2\varepsilon \mathbf{e}^0). \quad (47) $$

Right boundary cell. The ADER-IPDG scheme for the right cell reads

$$ \underline{N} \mathbf{u}_{N_x}^{n+1} = \underline{N} \mathbf{u}_{N_x}^n + \Delta t (\underline{Z}_L \mathbf{q}_{N_x-1}^n + \underline{Z}_C \mathbf{q}_{N_x}^n) + \Delta t \underline{P} \mathbf{S}_{N_x}^n + \Delta t \mathbf{b}_R^n, \quad (48) $$

where, again, all the local matrices are the same as for interior cells, except for $\underline{Z}_C^d$ in which a coefficient 2 needs to be added in front of the matrices $\underline{R}^{1,1}$ and $\underline{W}^{1,1}$. The advective and diffusive components for the boundary conditions are

$$ (\mathbf{b}_R^a)^n := a^- \mathbf{d}^1, \quad (\mathbf{b}_R^d)^n := -\frac{\lambda}{2} (\mu \mathbf{d}^1 + 2\varepsilon \mathbf{e}^1). \quad (49) $$

The remark on the computational cost that we made on the prediction step also holds for this correction stage, regardless the IPDG method chosen. Indeed, as we consider that $a$ and $\lambda$ are constants and we use a uniform mesh, all the matrices involved are constant and can be pre-computed only once. Moreover, this linear system to be solved for $\mathbf{u}_i^n$ is of small size $K$ and the matrix $\underline{N}$ is diagonal, which is a consequence of our choice of an orthogonal Lagrange polynomial basis. On the other hand, and contrary to FV methods that need a large stencil to reconstruct the information with a high precision, this DG scheme is compact (it is enough to know the DOF of the cell $\omega_i$ and of its two immediate neighbours) and thus the parallelization is simple and efficient.

# 4 Numerical validation of the high-order ADER-IPDG scheme

In this section, we present some numerical results on two different test cases in one dimension. The first one is a classical test: a gaussian that remains far from the boundaries. The second one is a bit more challenging: a smooth step initially located out of the domain and which crosses the boundary. Finally, for the third example, we propose a two dimensional extension.

## 4.1 First example in 1D

In this section, we present a manufactured test case. We take an advection velocity $a = 1$ and a diffusion coefficient $\lambda = 1$ on a computational domain
14 High order ADER-IPDG methods for advection-diffusion

$\Omega := [0, 3]$. The Dirichlet boundary conditions and the initial condition are imposed by using the values of the exact solution at the boundaries of $\Omega$ or at time $t = 0$ respectively, and the final time is $T_f = 1$. In all our simulations, we take $\gamma = 10$ in the diffusive numerical flux. We enforce the source term $S := \partial_t u + \partial_x F(u)$ corresponding to the exact solution below

$$u(x, t) := u_0(x - at)e^{-\lambda t}, \quad u_0(x) := e^{-30(x-1)^2}. \qquad (50)$$

For instance, the results obtained using an ADER-SIPG method on a space mesh composed of $N_x = 30$ cells with space-time $\mathbb{P}_4$ polynomials are presented in figure 3, where the dots represent the values of the numerical solution at each Gauss-Legendre quadrature point on every cell of the mesh.

<figure><img src="image_3.png" /><figcaption>Fig. 3 Test case 1: numerical and exact solutions obtained by the ADER-SIPG scheme with $\mathbb{P}_4$ polynomials on $N_x = 30$ cells.</figcaption></figure>

To validate numerically the theoretical orders of convergence of the global ADER-IPDG schemes, we need to examine the convergence rates of the three methods presented.

## Convergence study

On the same test case based on manufactured solution (50), we now perform a study of the mesh convergence of the SIPG, IIPG and NIPG methods. It is well known that the order of convergence of the DG scheme for pure advection is always $(k + 1)$ for $\mathbb{P}_k$ polynomials. However, we know from [13] that, in the $L^2$ norm, the orders of convergence of the different IPDG methods for the pure diffusion problem are supposed to be

* For IIPG and NIPG: order $(k + 1)$ for odd $k$ and order $k$ for even $k$.

* For SIPG: order $(k + 1)$ no matter if $k$ is odd or even.
14 High order ADER-IPDG methods for advection-diffusion

$\Omega := [0, 3]$. The Dirichlet boundary conditions and the initial condition are imposed by using the values of the exact solution at the boundaries of $\Omega$ or at time $t = 0$ respectively, and the final time is $T_f = 1$. In all our simulations, we take $\gamma = 10$ in the diffusive numerical flux. We enforce the source term $S := \partial_t u + \partial_x F(u)$ corresponding to the exact solution below

$$u(x, t) := u_0(x - at)e^{-\lambda t}, \quad u_0(x) := e^{-30(x-1)^2}. \qquad (50)$$

For instance, the results obtained using an ADER-SIPG method on a space mesh composed of $N_x = 30$ cells with space-time $\mathbb{P}_4$ polynomials are presented in figure 3, where the dots represent the values of the numerical solution at each Gauss-Legendre quadrature point on every cell of the mesh.

<figure><img src="image_5.png" /><figcaption>Fig. 3 Test case 1: numerical and exact solutions obtained by the ADER-SIPG scheme with $\mathbb{P}_4$ polynomials on $N_x = 30$ cells.</figcaption></figure>

To validate numerically the theoretical orders of convergence of the global ADER-IPDG schemes, we need to examine the convergence rates of the three methods presented.

## Convergence study

On the same test case based on manufactured solution (50), we now perform a study of the mesh convergence of the SIPG, IIPG and NIPG methods. It is well known that the order of convergence of the DG scheme for pure advection is always $(k + 1)$ for $\mathbb{P}_k$ polynomials. However, we know from [13] that, in the $L^2$ norm, the orders of convergence of the different IPDG methods for the pure diffusion problem are supposed to be

* For IIPG and NIPG: order $(k + 1)$ for odd $k$ and order $k$ for even $k$.

* For SIPG: order $(k + 1)$ no matter if $k$ is odd or even.
As a consequence, for the complete advection-diffusion equation, the resulting order of convergence of a given scheme is the minimum of the orders of con-vergence for the pure advection or pure diffusion problems, which comes down to the one for pure diffusion as it is the most restrictive. For more details on this theoretical result, the reader may refer to [14].

We thus compare the three IPDG methods by considering a series of meshes composed of $N_x := 16 \times 2^m$ cells for $m \in \{0, \dots, 3\}$, with $\mathbb{P}_k$ polynomials for all $k \in \{1, \dots, 4\}$. The resulting mesh sizes range from $3 \times 16^{-1}$ to $3 \times 128^{-1}$. The numerical orders of convergence are presented on figure 4.

<figure><img src="image_4.png" /><figcaption>Fig. 4 Test case 1: orders of convergence of the three ADER-IPDG methods.</figcaption></figure>

From the theory, we expect the order of convergence for an advection-diffusion problem to always be $(k+1)$ for $\mathbb{P}_k$ polynomials when using the SIPG scheme. On the contrary, this order may depend on the parity of $k$ with the IIPG or NIPG methods. In our results, the numerical orders of convergence are 2 and 4 respectively for $\mathbb{P}_1$ and $\mathbb{P}_3$ polynomials for all the methods. For a $\mathbb{P}_2$ representation of the solution, the SIPG method is of order 3 as expected, while the two others are slightly better than order 2. On the other hand, for $\mathbb{P}_4$ polynomials, if the fifth order is reached by the SIPG approach, its advantage over the NIPG and especially the IIPG methods is less evident in this case. The convergence study is thus globally compatible with the theoretical results.
As a consequence, for the complete advection-diffusion equation, the resulting order of convergence of a given scheme is the minimum of the orders of convergence for the pure advection or pure diffusion problems, which comes down to the one for pure diffusion as it is the most restrictive. For more details on this theoretical result, the reader may refer to [14].

We thus compare the three IPDG methods by considering a series of meshes composed of $N_x := 16 \times 2^m$ cells for $m \in \{0, \dots, 3\}$, with $\mathbb{P}_k$ polynomials for all $k \in \{1, \dots, 4\}$. The resulting mesh sizes range from $3 \times 16^{-1}$ to $3 \times 128^{-1}$. The numerical orders of convergence are presented on figure 4.

<figure><img src="image_6.png" /><figcaption>Fig. 4 Test case 1: orders of convergence of the three ADER-IPDG methods.</figcaption></figure>

From the theory, we expect the order of convergence for an advection-diffusion problem to always be $(k+1)$ for $\mathbb{P}_k$ polynomials when using the SIPG scheme. On the contrary, this order may depend on the parity of $k$ with the IIPG or NIPG methods. In our results, the numerical orders of convergence are 2 and 4 respectively for $\mathbb{P}_1$ and $\mathbb{P}_3$ polynomials for all the methods. For a $\mathbb{P}_2$ representation of the solution, the SIPG method is of order 3 as expected, while the two others are slightly better than order 2. On the other hand, for $\mathbb{P}_4$ polynomials, if the fifth order is reached by the SIPG approach, its advantage over the NIPG and especially the IIPG methods is less evident in this case. The convergence study is thus globally compatible with the theoretical results.
16 *High order ADER-IPDG methods for advection-diffusion*

Now, the performances of the different ADER-IPDG methods need to be further analyzed, especially in terms of execution time to reach a certain precision and with a convergence study.

*Performances*

For $\mathbb{P}_k$ polynomials for $k \in \{1, \dots, 4\}$, we propose to determine the number of cells $N_x$ required so that the $\ell^2$ relative error is less than $5 \times 10^{-2}\%$ or $5 \times 10^{-3}\%$ respectively. This study allows to compare the different $\mathbb{P}_k$ approximations regarding the total number $N_{\text{DOF}} := (k+1)N_x$ of DOF on the whole mesh as well as the execution time (on 6 threads) in order to reach a given precision. The results are presented in tables 2 and 3 respectively.

<table><caption>**Table 2** Test case 1: total number of DOF and execution time of the three IPDG methods to reach a precision of $5 \times 10^{-2}\%$.</caption><thead><tr><th>Poly.</th><th>IP method</th><th>N<sub>x</sub></th><th>N<sub>DOF</sub></th><th>Time (s)</th></tr></thead><tbody><tr><td rowspan="3">P<sub>1</sub></td><td>SIPG</td><td>306</td><td>612</td><td>144</td></tr><tr><td>IIPG</td><td>247</td><td>494</td><td>73</td></tr><tr><td>NIPG</td><td>199</td><td>398</td><td>48</td></tr><tr><td rowspan="3">idrule<br>P<sub>2</sub></td><td>SIPG</td><td>54</td><td>162</td><td>11</td></tr><tr><td>IIPG</td><td>70</td><td>210</td><td>24</td></tr><tr><td>NIPG</td><td>89</td><td>267</td><td>64</td></tr><tr><td rowspan="3">P<sub>3</sub></td><td>SIPG</td><td>26</td><td>104</td><td>10</td></tr><tr><td>IIPG</td><td>26</td><td>104</td><td>11</td></tr><tr><td>NIPG</td><td>27</td><td>108</td><td>13</td></tr><tr><td rowspan="3">P<sub>4</sub></td><td>SIPG</td><td>16</td><td>80</td><td>11</td></tr><tr><td>IIPG</td><td>16</td><td>80</td><td>13</td></tr><tr><td>NIPG</td><td>16</td><td>80</td><td>17</td></tr></tbody></table>

<table><caption>**Table 3** Test case 1: total number of DOF and execution time of the three IPDG methods to reach a precision of $5 \times 10^{-3}\%$.</caption><thead><tr><th>Poly.</th><th>IP method</th><th>N<sub>x</sub></th><th>N<sub>DOF</sub></th><th>Time (s)</th></tr></thead><tbody><tr><td rowspan="3">P<sub>2</sub></td><td>SIPG</td><td>116</td><td>348</td><td>109</td></tr><tr><td>IIPG</td><td>202</td><td>606</td><td>594</td></tr><tr><td>NIPG</td><td>275</td><td>825</td><td>1,623</td></tr><tr><td rowspan="3">P<sub>3</sub></td><td>SIPG</td><td>47</td><td>188</td><td>56</td></tr><tr><td>IIPG</td><td>46</td><td>184</td><td>56</td></tr><tr><td>NIPG</td><td>47</td><td>188</td><td>71</td></tr><tr><td rowspan="3">P<sub>4</sub></td><td>SIPG</td><td>25</td><td>125</td><td>51</td></tr><tr><td>IIPG</td><td>26</td><td>130</td><td>63</td></tr><tr><td>NIPG</td><td>26</td><td>130</td><td>68</td></tr></tbody></table>
16 *High order ADER-IPDG methods for advection-diffusion*

Now, the performances of the different ADER-IPDG methods need to be further analyzed, especially in terms of execution time to reach a certain precision and with a convergence study.

*Performances*

For $\mathbb{P}_k$ polynomials for $k \in \{1, \dots, 4\}$, we propose to determine the number of cells $N_x$ required so that the $\ell^2$ relative error is less than $5 \times 10^{-2}\%$ or $5 \times 10^{-3}\%$ respectively. This study allows to compare the different $\mathbb{P}_k$ approximations regarding the total number $N_{\text{DOF}} := (k+1)N_x$ of DOF on the whole mesh as well as the execution time (on 6 threads) in order to reach a given precision. The results are presented in tables 2 and 3 respectively.

<table><caption>**Table 2** Test case 1: total number of DOF and execution time of the three IPDG methods to reach a precision of $5 \times 10^{-2}\%$.</caption><thead><tr><th>Poly.</th><th>IP method</th><th>N<sub>x</sub></th><th>N<sub>DOF</sub></th><th>Time (s)</th></tr></thead><tbody><tr><td rowspan="3">P<sub>1</sub></td><td>SIPG</td><td>306</td><td>612</td><td>144</td></tr><tr><td>IIPG</td><td>247</td><td>494</td><td>73</td></tr><tr><td>NIPG</td><td>199</td><td>398</td><td>48</td></tr><tr><td rowspan="3">idrule<br>P<sub>2</sub></td><td>SIPG</td><td>54</td><td>162</td><td>11</td></tr><tr><td>IIPG</td><td>70</td><td>210</td><td>24</td></tr><tr><td>NIPG</td><td>89</td><td>267</td><td>64</td></tr><tr><td rowspan="3">P<sub>3</sub></td><td>SIPG</td><td>26</td><td>104</td><td>10</td></tr><tr><td>IIPG</td><td>26</td><td>104</td><td>11</td></tr><tr><td>NIPG</td><td>27</td><td>108</td><td>13</td></tr><tr><td rowspan="3">P<sub>4</sub></td><td>SIPG</td><td>16</td><td>80</td><td>11</td></tr><tr><td>IIPG</td><td>16</td><td>80</td><td>13</td></tr><tr><td>NIPG</td><td>16</td><td>80</td><td>17</td></tr></tbody></table>

<table><caption>**Table 3** Test case 1: total number of DOF and execution time of the three IPDG methods to reach a precision of $5 \times 10^{-3}\%$.</caption><thead><tr><th>Poly.</th><th>IP method</th><th>N<sub>x</sub></th><th>N<sub>DOF</sub></th><th>Time (s)</th></tr></thead><tbody><tr><td rowspan="3">P<sub>2</sub></td><td>SIPG</td><td>116</td><td>348</td><td>109</td></tr><tr><td>IIPG</td><td>202</td><td>606</td><td>594</td></tr><tr><td>NIPG</td><td>275</td><td>825</td><td>1,623</td></tr><tr><td rowspan="3">P<sub>3</sub></td><td>SIPG</td><td>47</td><td>188</td><td>56</td></tr><tr><td>IIPG</td><td>46</td><td>184</td><td>56</td></tr><tr><td>NIPG</td><td>47</td><td>188</td><td>71</td></tr><tr><td rowspan="3">P<sub>4</sub></td><td>SIPG</td><td>25</td><td>125</td><td>51</td></tr><tr><td>IIPG</td><td>26</td><td>130</td><td>63</td></tr><tr><td>NIPG</td><td>26</td><td>130</td><td>68</td></tr></tbody></table>
A fair remark is that when a smaller relative error is required, the schemes of higher orders seem preferable for their speed and low memory footprint regarding the total number of DOF, which is shown in table 3 for $P_3$ and especially $P_4$ polynomials. More precisely, the execution times show that when requiring a relative error of $5 \times 10^{-2}\%$, $P_1$ schemes are not competitive while higher degree approximations remain comparable. If we ask for an error divided by a factor 10, the same comment can be made on $P_2$ approximations of the solution. In addition, the SIPG variant appears to be generally more competitive compared to the other IP methods regarding the execution time needed to reach a certain precision. It also always needs less DOF, which is an advantage in terms of memory footprint.

## 4.2 Second example in 1D

Here, we propose another manufactured test case that reproduces the same physical behaviours, with an advective velocity $a = 1$ and a diffusion coefficient $\lambda = 1$. The computational domain is $\Omega := [0, 1]$ and the final time of the simulation is $T_f = 1$. We take $\gamma = 10$ in the diffusive numerical flux. The Dirichlet boundary conditions, the initial condition and the source term $S$ are imposed by using the values of the exact solution below

$$u(x,t) := u_0(x - at)e^{-\lambda t}, \quad u_0(x) := 1 - \tanh(10(x + 0.4)). \quad (51)$$

The results of the ADER-SIPG scheme with a space mesh composed of $N_x = 20$ and $P_4$ polynomials are presentd in figure 5, where the dots represent the values of the numerical solution at the Gauss-Legendre quadrature points in every cell.

<figure><img src="image_5.png" /><figcaption>Fig. 5 Test case 2: numerical and exact solutions obtained by the ADER-SIPG scheme with $P_4$ polynomials on $N_x = 20$ cells.</figcaption></figure>
A fair remark is that when a smaller relative error is required, the schemes of higher orders seem preferable for their speed and low memory footprint regarding the total number of DOF, which is shown in table 3 for $\mathbb{P}_3$ and especially $\mathbb{P}_4$ polynomials. More precisely, the execution times show that when requiring a relative error of $5 \times 10^{-2}\%$, $\mathbb{P}_1$ schemes are not competitive while higher degree approximations remain comparable. If we ask for an error divided by a factor 10, the same comment can be made on $\mathbb{P}_2$ approximations of the solution. In addition, the SIPG variant appears to be generally more competitive compared to the other IP methods regarding the execution time needed to reach a certain precision. It also always needs less DOF, which is an advantage in terms of memory footprint.

## 4.2 Second example in 1D

Here, we propose another manufactured test case that reproduces the same physical behaviours, with an advective velocity $a = 1$ and a diffusion coefficient $\lambda = 1$. The computational domain is $\Omega := [0, 1]$ and the final time of the simulation is $T_f = 1$. We take $\gamma = 10$ in the diffusive numerical flux. The Dirichlet boundary conditions, the initial condition and the source term $S$ are imposed by using the values of the exact solution below

$$u(x,t) := u_0(x - at)e^{-\lambda t}, \quad u_0(x) := 1 - \tanh(10(x + 0.4)). \quad (51)$$

The results of the ADER-SIPG scheme with a space mesh composed of $N_x = 20$ and $\mathbb{P}_4$ polynomials are presentd in figure 5, where the dots represent the values of the numerical solution at the Gauss-Legendre quadrature points in every cell.

<figure><img src="image_7.png" /><figcaption>Fig. 5 Test case 2: numerical and exact solutions obtained by the ADER-SIPG scheme with $\mathbb{P}_4$ polynomials on $N_x = 20$ cells.</figcaption></figure>
18 *High order ADER-IPDG methods for advection-diffusion*

A similar convergence study is performed on this test case. The goal is to
determine if the boundary conditions (especially on the left boundary of the
domain) deteriorate the global order of accuracy of the different ADER-IPDG
schemes or if the method remains consistent to the expected precision.

Convergence study

On the same test case based on manufactured solution (51), we carry out a
convergence study of the SIPG, IIPG and NIPG methods, keeping in mind
that, in the L² norm, the orders of convergence of the different IPDG methods
for the complete advection-diffusion problem should be (k + 1) for all of them
except for IIPG and NIPG with even k for which the order is supposed to be k.

In order to compare the three IPDG methods, we consider a series of meshes composed of $N_x := 8 \times 2^m$ cells for $m \in \{0, \dots, 3\}$, with $\mathbb{P}_k$ polynomials for all $k \in \{1, \dots, 4\}$. The resulting mesh sizes now range from $2 \times 16^{-1}$ to $2 \times 128^{-1}$, which is comparable to the sizes of the cells for the first test case as the domain $\Omega$ is three times smaller here. Again, the numerical orders of convergence are shown in Figure 6.

<figure><img src="image_6.png" /><figcaption>**Fig. 6** Test case 2: orders of convergence of the three ADER-IPDG methods.</figcaption></figure>

The conclusion here is the same as for the first test case: the numerical
orders of convergence are coherent with those expected. More precisely, the
18 *High order ADER-IPDG methods for advection-diffusion*

A similar convergence study is performed on this test case. The goal is to
determine if the boundary conditions (especially on the left boundary of the
domain) deteriorate the global order of accuracy of the different ADER-IPDG
schemes or if the method remains consistent to the expected precision.

Convergence study

On the same test case based on manufactured solution (51), we carry out a
convergence study of the SIPG, IIPG and NIPG methods, keeping in mind
that, in the L² norm, the orders of convergence of the different IPDG methods
for the complete advection-diffusion problem should be (k + 1) for all of them
except for IIPG and NIPG with even k for which the order is supposed to be k.

In order to compare the three IPDG methods, we consider a series of meshes composed of $N_x := 8 \times 2^m$ cells for $m \in \{0, \dots, 3\}$, with $\mathbb{P}_k$ polynomials for all $k \in \{1, \dots, 4\}$. The resulting mesh sizes now range from $2 \times 16^{-1}$ to $2 \times 128^{-1}$, which is comparable to the sizes of the cells for the first test case as the domain $\Omega$ is three times smaller here. Again, the numerical orders of convergence are shown in Figure 6.

<figure><img src="image_8.png" /><figcaption>**Fig. 6** Test case 2: orders of convergence of the three ADER-IPDG methods.</figcaption></figure>

The conclusion here is the same as for the first test case: the numerical
orders of convergence are coherent with those expected. More precisely, the
SIPG method is always of optimal order ($k+1$) when the mesh is refined. The IIPG and NIPG methods depend on the parity of $k$, which is well documented (see [13, 24, 28]). Taking finer meshes may allow to better show the advantage of the SIPG method compared to the other two in the case of $\mathbb{P}_4$ polynomials, in relation to their respective orders of convergence.

## 4.3 Third example: 2D case

This section is dedicated to the two dimensional (2D) version of the first manufactured test case. As we use a cartesian mesh, writing both the predictor and the corrector step of the ADER-IPDG scheme is a straightforward extension of the 1D study presented in this paper.

We take the advection velocities $a = b = 1$ with respect to the $x$ and $y$ directions, and a diffusion coefficient $\lambda = 1$ on a computational domain $\Omega := [0, 3]^2$. The final time is $T_f = 1$ and we keep the value $\gamma = 10$ in the diffusive numerical flux. The Dirichlet boundary conditions, the initial condition and the source term imposed are those corresponding to the following exact solution

$$u(x, y, t) := u_0(x - at, y - bt)e^{-\lambda t}, \quad u_0(x, y) := e^{-30[(x-1)^2+(y-1)^2]}. \quad (52)$$

A convergence study is carried out on this 2D example in order to show that the 2D extension remains coherent with respect to the theoretical results announced in 1D.

### Convergence study

On test case (52), we perform a study of the mesh convergence of the same three ADER-IPDG methods. We recall that the order of convergence of the SIPG scheme is expected to be still ($k+1$) for $\mathbb{P}_k$ polynomials, while it should be $k$ or ($k+1$) if $k$ is even or odd respectively.

We thus compare the three IPDG variants by considering a series of meshes composed of $N_x = N_y := 32 \times 2^{m-k+1}$ cells for $m \in \{0, \dots, 3\}$, with $\mathbb{P}_k$ polynomials for all $k \in \{1, \dots, 4\}$. This choice of shifting the range of meshes used for higher degree polynomials allows to reduce the computational cost of the study. In this sense, we also fix the duration of the simulation to $T_f = 0.01$. The numerical orders of convergence are shown on figure 7, where we have defined $h = \Delta x = \Delta y$.
SIPG method is always of optimal order ($k+1$) when the mesh is refined. The IIPG and NIPG methods depend on the parity of $k$, which is well documented (see [13, 24, 28]). Taking finer meshes may allow to better show the advantage of the SIPG method compared to the other two in the case of $\mathbb{P}_4$ polynomials, in relation to their respective orders of convergence.

## 4.3 Third example: 2D case

This section is dedicated to the two dimensional (2D) version of the first manufactured test case. As we use a cartesian mesh, writing both the predictor and the corrector step of the ADER-IPDG scheme is a straightforward extension of the 1D study presented in this paper.

We take the advection velocities $a = b = 1$ with respect to the $x$ and $y$ directions, and a diffusion coefficient $\lambda = 1$ on a computational domain $\Omega := [0, 3]^2$. The final time is $T_f = 1$ and we keep the value $\gamma = 10$ in the diffusive numerical flux. The Dirichlet boundary conditions, the initial condition and the source term imposed are those corresponding to the following exact solution

$$u(x, y, t) := u_0(x - at, y - bt)e^{-\lambda t}, \quad u_0(x, y) := e^{-30[(x-1)^2+(y-1)^2]}. \quad (52)$$

A convergence study is carried out on this 2D example in order to show that the 2D extension remains coherent with respect to the theoretical results announced in 1D.

### Convergence study

On test case (52), we perform a study of the mesh convergence of the same three ADER-IPDG methods. We recall that the order of convergence of the SIPG scheme is expected to be still ($k+1$) for $\mathbb{P}_k$ polynomials, while it should be $k$ or ($k+1$) if $k$ is even or odd respectively.

We thus compare the three IPDG variants by considering a series of meshes composed of $N_x = N_y := 32 \times 2^{m-k+1}$ cells for $m \in \{0, \dots, 3\}$, with $\mathbb{P}_k$ polynomials for all $k \in \{1, \dots, 4\}$. This choice of shifting the range of meshes used for higher degree polynomials allows to reduce the computational cost of the study. In this sense, we also fix the duration of the simulation to $T_f = 0.01$. The numerical orders of convergence are shown on figure 7, where we have defined $h = \Delta x = \Delta y$.
20 High order ADER-IPDG methods for advection-diffusion

<figure><img src="image_7.png" /><figcaption>**Fig. 7** Test case 3: orders of convergence of the three ADER-IPDG methods.</figcaption></figure>

For the $P_1$ and $P_3$ approximations, the three ADER-IPDG methods provide the expected order of accuracy, that is to say 2 and 4 respectively. The same observation can be made concerning the $P_2$ schemes: the numerical convergence rates are coherent with the theoretical expectations. In the case of $P_4$ polynomials, no significant difference between the SIPG method and the other two can be noticed yet. We expect that by using finer meshes, one can recover the results obtained in 1D. Indeed, in this 2D convergence study, the finer mesh is around $h = 10^{-1}$: at this level of refinement, the 1D curves for $P_4$ approximations did not reach the asymptotic behaviour yet.

As we did for the 1D version of this test case, we now present a study of the performance of the three ADER-IPDG schemes.

## Performances

For $P_k$ polynomials for $k \in \{1, \dots, 4\}$, we want to determine the total number $N^2$ of cells (i.e. $N_x = N_y = N$ cells in each space direction) to reach an $\ell^2$ relative error of 0.1% or 0.05% respectively. We do not ask for the same errors as in the 1D case because of the increased computational cost in 2D. Thus, we propose to compare the different $P_k$ approximations regarding the total number $N_{\text{DOF}} := (k+1)^2 N^2$ of DOF and the execution time (again, on 6 threads) to reach a prescribed precision. The results are shown in tables 4 and 5 respectively.
20 High order ADER-IPDG methods for advection-diffusion

<figure><img src="image_9.png" /><figcaption>Fig. 7 Test case 3: orders of convergence of the three ADER-IPDG methods.</figcaption></figure>

For the $P_1$ and $P_3$ approximations, the three ADER-IPDG methods provide the expected order of accuracy, that is to say 2 and 4 respectively. The same observation can be made concerning the $P_2$ schemes: the numerical convergence rates are coherent with the theoretical expectations. In the case of $P_4$ polynomials, no significant difference between the SIPG method and the other two can be noticed yet. We expect that by using finer meshes, one can recover the results obtained in 1D. Indeed, in this 2D convergence study, the finer mesh is around $h = 10^{-1}$: at this level of refinement, the 1D curves for $P_4$ approximations did not reach the asymptotic behaviour yet.

As we did for the 1D version of this test case, we now present a study of the performance of the three ADER-IPDG schemes.

## Performances

For $P_k$ polynomials for $k \in \{1, \dots, 4\}$, we want to determine the total number $N^2$ of cells (i.e. $N_x = N_y = N$ cells in each space direction) to reach an $\ell^2$ relative error of 0.1% or 0.05% respectively. We do not ask for the same errors as in the 1D case because of the increased computational cost in 2D. Thus, we propose to compare the different $P_k$ approximations regarding the total number $N_{\text{DOF}} := (k+1)^2 N^2$ of DOF and the execution time (again, on 6 threads) to reach a prescribed precision. The results are shown in tables 4 and 5 respectively.
<table><caption>**Table 4** Test case 3: total number of DOF and execution time of the three IPDG methods to reach a precision of 0.1%.</caption><thead><tr><th>Poly.</th><th>IP method</th><th>N<sup>2</sup></th><th>N<sub>DOF</sub></th><th>Time (s)</th></tr></thead><tbody><tr><td rowspan="3">P<sub>1</sub></td><td>SIPG</td><td>207<sup>2</sup></td><td>171,396</td><td>539</td></tr><tr><td>IIPG</td><td>165<sup>2</sup></td><td>108,900</td><td>238</td></tr><tr><td>NIPG</td><td>131<sup>2</sup></td><td>68,644</td><td>100</td></tr><tr><td rowspan="3">P<sub>2</sub></td><td>SIPG</td><td>48<sup>2</sup></td><td>20,736</td><td>68</td></tr><tr><td>IIPG</td><td>54<sup>2</sup></td><td>26,244</td><td>107</td></tr><tr><td>NIPG</td><td>65<sup>2</sup></td><td>38,025</td><td>273</td></tr><tr><td rowspan="3">P<sub>3</sub></td><td>SIPG</td><td>23<sup>2</sup></td><td>8,464</td><td>76</td></tr><tr><td>IIPG</td><td>24<sup>2</sup></td><td>9,216</td><td>96</td></tr><tr><td>NIPG</td><td>24<sup>2</sup></td><td>9,216</td><td>88</td></tr><tr><td rowspan="3">P<sub>4</sub></td><td>SIPG</td><td>16<sup>2</sup></td><td>6,400</td><td>271</td></tr><tr><td>IIPG</td><td>16<sup>2</sup></td><td>6,400</td><td>297</td></tr><tr><td>NIPG</td><td>16<sup>2</sup></td><td>6,400</td><td>326</td></tr></tbody></table>

<table><caption>**Table 5** Test case 3: total number of DOF and execution time of the three IPDG methods to reach a precision of 0.05%.</caption><thead><tr><th>Poly.</th><th>IP method</th><th>N<sup>2</sup></th><th>N<sub>DOF</sub></th><th>Time (s)</th></tr></thead><tbody><tr><td rowspan="3">P<sub>2</sub></td><td>SIPG</td><td>60<sup>2</sup></td><td>32,400</td><td>158</td></tr><tr><td>IIPG</td><td>72<sup>2</sup></td><td>46,656</td><td>372</td></tr><tr><td>NIPG</td><td>89<sup>2</sup></td><td>71,289</td><td>1,026</td></tr><tr><td rowspan="3">P<sub>3</sub></td><td>SIPG</td><td>28<sup>2</sup></td><td>12,544</td><td>149</td></tr><tr><td>IIPG</td><td>28<sup>2</sup></td><td>12,544</td><td>170</td></tr><tr><td>NIPG</td><td>29<sup>2</sup></td><td>13,456</td><td>247</td></tr><tr><td rowspan="3">P<sub>4</sub></td><td>SIPG</td><td>16<sup>2</sup></td><td>6,400</td><td>271</td></tr><tr><td>IIPG</td><td>16<sup>2</sup></td><td>6,400</td><td>297</td></tr><tr><td>NIPG</td><td>16<sup>2</sup></td><td>6,400</td><td>326</td></tr></tbody></table>

From the two tables, the numerical results show that the SIPG method is almost always faster than its two variants to reach a given precision. In this 2D framework, it also has to be noticed that higher order schemes provide a reasonable advantage regarding the memory footprint and the execution time. From table 4, the $P_4$ schemes may not seem competitive, but table 5 shows that they even allow to reach an error divided by a factor 2 (0.05% instead of 0.1%) without refining the mesh. It is due to the high order of the method that provides a strongly decreasing error when the mesh is refined. Thus, the 2D performance study finally joins the analysis made for the 1D case.
<table><caption>**Table 4** Test case 3: total number of DOF and execution time of the three IPDG methods to reach a precision of 0.1%.</caption><thead><tr><th>Poly.</th><th>IP method</th><th>N<sup>2</sup></th><th>N<sub>DOF</sub></th><th>Time (s)</th></tr></thead><tbody><tr><td rowspan="3">P<sub>1</sub></td><td>SIPG</td><td>207<sup>2</sup></td><td>171,396</td><td>539</td></tr><tr><td>IIPG</td><td>165<sup>2</sup></td><td>108,900</td><td>238</td></tr><tr><td>NIPG</td><td>131<sup>2</sup></td><td>68,644</td><td>100</td></tr><tr><td rowspan="3">P<sub>2</sub></td><td>SIPG</td><td>48<sup>2</sup></td><td>20,736</td><td>68</td></tr><tr><td>IIPG</td><td>54<sup>2</sup></td><td>26,244</td><td>107</td></tr><tr><td>NIPG</td><td>65<sup>2</sup></td><td>38,025</td><td>273</td></tr><tr><td rowspan="3">P<sub>3</sub></td><td>SIPG</td><td>23<sup>2</sup></td><td>8,464</td><td>76</td></tr><tr><td>IIPG</td><td>24<sup>2</sup></td><td>9,216</td><td>96</td></tr><tr><td>NIPG</td><td>24<sup>2</sup></td><td>9,216</td><td>88</td></tr><tr><td rowspan="3">P<sub>4</sub></td><td>SIPG</td><td>16<sup>2</sup></td><td>6,400</td><td>271</td></tr><tr><td>IIPG</td><td>16<sup>2</sup></td><td>6,400</td><td>297</td></tr><tr><td>NIPG</td><td>16<sup>2</sup></td><td>6,400</td><td>326</td></tr></tbody></table>

<table><caption>**Table 5** Test case 3: total number of DOF and execution time of the three IPDG methods to reach a precision of 0.05%.</caption><thead><tr><th>Poly.</th><th>IP method</th><th>N<sup>2</sup></th><th>N<sub>DOF</sub></th><th>Time (s)</th></tr></thead><tbody><tr><td rowspan="3">P<sub>2</sub></td><td>SIPG</td><td>60<sup>2</sup></td><td>32,400</td><td>158</td></tr><tr><td>IIPG</td><td>72<sup>2</sup></td><td>46,656</td><td>372</td></tr><tr><td>NIPG</td><td>89<sup>2</sup></td><td>71,289</td><td>1,026</td></tr><tr><td rowspan="3">P<sub>3</sub></td><td>SIPG</td><td>28<sup>2</sup></td><td>12,544</td><td>149</td></tr><tr><td>IIPG</td><td>28<sup>2</sup></td><td>12,544</td><td>170</td></tr><tr><td>NIPG</td><td>29<sup>2</sup></td><td>13,456</td><td>247</td></tr><tr><td rowspan="3">P<sub>4</sub></td><td>SIPG</td><td>16<sup>2</sup></td><td>6,400</td><td>271</td></tr><tr><td>IIPG</td><td>16<sup>2</sup></td><td>6,400</td><td>297</td></tr><tr><td>NIPG</td><td>16<sup>2</sup></td><td>6,400</td><td>326</td></tr></tbody></table>

From the two tables, the numerical results show that the SIPG method is almost always faster than its two variants to reach a given precision. In this 2D framework, it also has to be noticed that higher order schemes provide a reasonable advantage regarding the memory footprint and the execution time. From table 4, the $P_4$ schemes may not seem competitive, but table 5 shows that they even allow to reach an error divided by a factor 2 (0.05% instead of 0.1%) without refining the mesh. It is due to the high order of the method that provides a strongly decreasing error when the mesh is refined. Thus, the 2D performance study finally joins the analysis made for the 1D case.
# 5 Conclusion and perspectives

In this paper, we proposed an original association of the IPDG framework with an ADER technique for the time integration to solve the unsteady linear advection-diffusion equation. This ADER-IPDG combination allowed to reach an arbitrary high order of convergence on smooth test cases and showed that higher order schemes are competitive to meet certain requirements regarding both the desired precision and the execution time with a compact stencil. Our numerical studies suggest that the symmetric version, SIPG, is preferable in most of the cases: it is always of optimal order of convergence, and its execution is always faster to reach a desired precision.

In [10] and [11], an ADER-FV scheme was developed with $\mathbb{P}_2$ polynomials for the advection-diffusion equation. As a consequence, the contribution from the diffusive part of the total flux was evaluated by using a relaxation equation to approximate the gradient of the solution. Although this approach granted the global consistency and accuracy of the ADER-FV scheme, the relaxation parameter deteriorated the CFL condition in the limit of mesh refinement. Our different high order ADER-IPDG methods overcome this difficulty at no additional cost.

The results obtained in this paper constitute a first step towards the solution of more challenging problems like the incompressible Navier-Stokes equations. In this sense, the compactness and precision of the IPDG schemes presented combined with the adaptivity of the ADER procedure to the deforming cells of Chimera meshes show promising perspectives.

# Appendix A Additional Definitions

Let us recall that the test functions considered in (25) are the space basis functions $\hat{\psi} \equiv \theta_k(\xi)$. In order to present all the necessary computations in a legible way, we first define the quantities

$$ (D_i^\tau)_k := [\theta_k, \hat{u}_i(\cdot, \tau)], \qquad (J_i)_k := \langle \theta_k, \hat{S}_i \rangle, \qquad (A1a) $$

$$ (G_i^a)_k := \langle \partial_\xi \theta_k, a \hat{q}_i \rangle, \qquad (G_i^d)_k := \langle \partial_\xi \theta_k, \lambda \partial_\xi \hat{q}_i \rangle, \qquad (A1b) $$

$$ (H_i^\xi)_k^a := \{\theta_k, \tilde{F}_a(\tilde{q})\}_\xi, \qquad (H_i^\xi)_k^d := \{\theta_k, \tilde{F}_d(\tilde{q})\}_\xi, \qquad (A1c) $$

$$ (\bar{H}_i^\xi)_k^d := \{\tilde{q}, \hat{F}_d(\theta_k)\}_\xi. \qquad (A1d) $$

The ADER-IPDG variational formulation (25) written for all test functions $\theta_k$ for $k \in \{1, \dots, K\}$ now reads

$$ \mathbf{D}_i^1 = \mathbf{D}_i^0 + \left(\frac{\Delta t}{\Delta x}\right) \mathbf{G}_i^a - \left(\frac{\Delta t}{\Delta x^2}\right) \mathbf{G}_i^d + \Delta t \mathbf{J}_i - \left(\frac{\Delta t}{\Delta x}\right) ((\mathbf{H}_i^1)^a - (\mathbf{H}_i^0)^a) $$
# 5 Conclusion and perspectives

In this paper, we proposed an original association of the IPDG framework with an ADER technique for the time integration to solve the unsteady linear advection-diffusion equation. This ADER-IPDG combination allowed to reach an arbitrary high order of convergence on smooth test cases and showed that higher order schemes are competitive to meet certain requirements regarding both the desired precision and the execution time with a compact stencil. Our numerical studies suggest that the symmetric version, SIPG, is preferable in most of the cases: it is always of optimal order of convergence, and its execution is always faster to reach a desired precision.

In [10] and [11], an ADER-FV scheme was developed with $P_2$ polynomials for the advection-diffusion equation. As a consequence, the contribution from the diffusive part of the total flux was evaluated by using a relaxation equation to approximate the gradient of the solution. Although this approach granted the global consistency and accuracy of the ADER-FV scheme, the relaxation parameter deteriorated the CFL condition in the limit of mesh refinement. Our different high order ADER-IPDG methods overcome this difficulty at no additional cost.

The results obtained in this paper constitute a first step towards the solution of more challenging problems like the incompressible Navier-Stokes equations. In this sense, the compactness and precision of the IPDG schemes presented combined with the adaptivity of the ADER procedure to the deforming cells of Chimera meshes show promising perspectives.

# Appendix A Additional Definitions

Let us recall that the test functions considered in (25) are the space basis functions $\hat{\psi} \equiv \theta_k(\xi)$. In order to present all the necessary computations in a legible way, we first define the quantities

$$ (D_i^\tau)_k := [\theta_k, \hat{u}_i(\cdot, \tau)], \qquad (J_i)_k := \langle \theta_k, \hat{S}_i \rangle, \qquad (A1a) $$

$$ (G_i^a)_k := \langle \partial_\xi \theta_k, a \hat{q}_i \rangle, \qquad (G_i^d)_k := \langle \partial_\xi \theta_k, \lambda \partial_\xi \hat{q}_i \rangle, \qquad (A1b) $$

$$ (H_i^\xi)_k^a := \{\theta_k, \tilde{F}_a(\tilde{q})\}_\xi, \qquad (H_i^\xi)_k^d := \{\theta_k, \tilde{F}_d(\tilde{q})\}_\xi, \qquad (A1c) $$

$$ (\bar{H}_i^\xi)_k^d := \{\tilde{q}, \hat{F}_d(\theta_k)\}_\xi. \qquad (A1d) $$

The ADER-IPDG variational formulation (25) written for all test functions $\theta_k$ for $k \in \{1, \dots, K\}$ now reads

$$ \mathbf{D}_i^1 = \mathbf{D}_i^0 + \left(\frac{\Delta t}{\Delta x}\right) \mathbf{G}_i^a - \left(\frac{\Delta t}{\Delta x^2}\right) \mathbf{G}_i^d + \Delta t \mathbf{J}_i - \left(\frac{\Delta t}{\Delta x}\right) ((\mathbf{H}_i^1)^a - (\mathbf{H}_i^0)^a) $$
$$ - \left( \frac{\Delta t}{\Delta x} \right) \left( (\mathbf{H}_i^1)^d - (\mathbf{H}_i^0)^d \right) + \varepsilon \left( \frac{\Delta t}{\Delta x} \right) \left( (\bar{\mathbf{H}}_i^1)^d - (\bar{\mathbf{H}}_i^0)^d \right). \quad (A2) $$

The goal now is to express the different terms $\mathbf{D}_i$, $\mathbf{G}_i$, $\mathbf{H}_i$, $\bar{\mathbf{H}}_i$ as functions of the DOF of the Galerkin predictors on cells $\omega_j$ for $j \in \{i-1, i, i+1\}$ denoted by $\mathbf{q}_j$ in order to build the local linear systems to be solved. Let us reuse the polynomial expansions (14)–(23), the numerical fluxes (29)–(32) as well as the matrices (11)–(12)–(13). In the following, the quantities $\mathbf{H}_i$ and $\bar{\mathbf{H}}_i$ will be refered to as original terms and extra terms respectively.

The term of DOF.

$$ (D_i^T)_k = \sum_{\ell=1}^{K} [\theta_k, \theta_\ell] \hat{u}_\ell^i(\tau) \quad \Leftrightarrow \quad \mathbf{D}_i^T = \underline{N} \mathbf{u}_i^{n+\tau}. \quad (A3) $$

The advective volume term.

$$ (G_i^a)_k = a \sum_{\ell=1}^{L} [\partial_\xi \theta_k, \Theta_\ell] \hat{q}_\ell^i \quad \Leftrightarrow \quad \mathbf{G}_i^a = a \underline{C} \mathbf{q}_i^n. \quad (A4) $$

The diffusive volume term.

$$ (G_i^d)_k = \lambda \sum_{\ell=1}^{L} [\partial_\xi \theta_k, \partial_\xi \Theta_\ell] \hat{u}_\ell^i \quad \Leftrightarrow \quad \mathbf{G}_i^d = \lambda \underline{E} \mathbf{q}_i^n. \quad (A5) $$

The source term.

$$ (J_i)_k = \sum_{\ell=1}^{L} [\theta_k, \Theta_\ell] \hat{S}_\ell^i \quad \Leftrightarrow \quad \mathbf{J}_i = \underline{P} \mathbf{S}_i^n. \quad (A6) $$

We now introduce the notations $\mathbf{d}_k^\xi := [\theta_k(\xi), g_\pm(\cdot)]$ and $\mathbf{e}_k^\xi := [\partial_\xi \theta_k(\xi), g_\pm(\cdot)]$, where $g_\pm$ is the Dirichlet boundary condition. For any real number $z$, we define its positive and negative parts respectively by $z^+ := \frac{1}{2}(|z| + z) \ge 0$ and $z^- := \frac{1}{2}(|z| - z) \ge 0$. Then, we can compute the interface terms on a given cell.

The left advective original term.

• For interior interfaces

$$ \begin{aligned} (H_i^0)_k^a &= \left[ \theta_k(0), \frac{1}{2} \left( \hat{F}_a(\hat{q}_{i-1}(1, \cdot)) + \hat{F}_a(\hat{q}_i(0, \cdot)) \right) + \sigma_a(\hat{q}_{i-1}(1, \cdot) - \hat{q}_i(0, \cdot)) \right], \\ \Leftrightarrow (H_i^0)_k^a &= (a^+ \underline{Q}^{0,1}) \mathbf{q}_{i-1}^n + (-a^- \underline{Q}^{0,0}) \mathbf{q}_i^n. \end{aligned} \quad (A7) $$
$$ - \left( \frac{\Delta t}{\Delta x} \right) \left( (\mathbf{H}_i^1)^d - (\mathbf{H}_i^0)^d \right) + \varepsilon \left( \frac{\Delta t}{\Delta x} \right) \left( (\bar{\mathbf{H}}_i^1)^d - (\bar{\mathbf{H}}_i^0)^d \right). \qquad (A2) $$

The goal now is to express the different terms $\mathbf{D}_i$, $\mathbf{G}_i$, $\mathbf{H}_i$, $\bar{\mathbf{H}}_i$ as functions of the DOF of the Galerkin predictors on cells $\omega_j$ for $j \in \{i-1, i, i+1\}$ denoted by $\mathbf{q}_j$ in order to build the local linear systems to be solved. Let us reuse the polynomial expansions (14)–(23), the numerical fluxes (29)–(32) as well as the matrices (11)–(12)–(13). In the following, the quantities $\mathbf{H}_i$ and $\bar{\mathbf{H}}_i$ will be refered to as original terms and extra terms respectively.

The term of DOF.

$$ (D_i^T)_k = \sum_{\ell=1}^{K} [\theta_k, \theta_\ell] \hat{u}_\ell^i(\tau) \quad \Leftrightarrow \quad \mathbf{D}_i^T = \underline{N} \mathbf{u}_i^{n+\tau}. \qquad (A3) $$

The advective volume term.

$$ (G_i^a)_k = a \sum_{\ell=1}^{L} [\partial_\xi \theta_k, \Theta_\ell] \hat{q}_\ell^i \quad \Leftrightarrow \quad \mathbf{G}_i^a = a \underline{C} \mathbf{q}_i^n. \qquad (A4) $$

The diffusive volume term.

$$ (G_i^d)_k = \lambda \sum_{\ell=1}^{L} [\partial_\xi \theta_k, \partial_\xi \Theta_\ell] \hat{u}_\ell^i \quad \Leftrightarrow \quad \mathbf{G}_i^d = \lambda \underline{E} \mathbf{q}_i^n. \qquad (A5) $$

The source term.

$$ (J_i)_k = \sum_{\ell=1}^{L} [\theta_k, \Theta_\ell] \hat{S}_\ell^i \quad \Leftrightarrow \quad \mathbf{J}_i = \underline{P} \mathbf{S}_i^n. \qquad (A6) $$

We now introduce the notations $\mathbf{d}_k^\xi := [\theta_k(\xi), g_\pm(\cdot)]$ and $\mathbf{e}_k^\xi := [\partial_\xi \theta_k(\xi), g_\pm(\cdot)]$, where $g_\pm$ is the Dirichlet boundary condition. For any real number $z$, we define its positive and negative parts respectively by $z^+ := \frac{1}{2}(|z| + z) \ge 0$ and $z^- := \frac{1}{2}(|z| - z) \ge 0$. Then, we can compute the interface terms on a given cell.

The left advective original term.

• For interior interfaces

$$ \begin{aligned} (H_i^0)_k^a &= \left[ \theta_k(0), \frac{1}{2} \left( \hat{F}_a(\hat{q}_{i-1}(1, \cdot)) + \hat{F}_a(\hat{q}_i(0, \cdot)) \right) + \sigma_a(\hat{q}_{i-1}(1, \cdot) - \hat{q}_i(0, \cdot)) \right], \\ \Leftrightarrow (\mathbf{H}_i^0)^a &= (a^+ \underline{Q}^{0,1}) \mathbf{q}_{i-1}^n + (-a^- \underline{Q}^{0,0}) \mathbf{q}_i^n. \end{aligned} \qquad (A7) $$
24 High order ADER-IPDG methods for advection-diffusion

• For a left Dirichlet boundary interface

$$
\begin{align}
(H_0^0)_k^a &= \left[ \theta_k(0), \frac{1}{2} \left( \hat{F}_a(g_-(\cdot)) + \hat{F}_a(\hat{q}_1(0, \cdot)) \right) + \sigma_a (g_-(\cdot) - \hat{q}_1(0, \cdot)) \right], \nonumber \\
\Leftrightarrow (\mathbf{H}_0^0)^a &= a^+ \mathbf{d}^0 + (-a^- \underline{Q}^{0,0}) \mathbf{q}_1^n. \tag{A8}
\end{align}
$$

The right advective original term.

• For interior interfaces

$$
\begin{align}
(H_i^1)_k^a &= \left[ \theta_k(1), \frac{1}{2} \left( \hat{F}_a(\hat{q}_i(1, \cdot)) + \hat{F}_a(\hat{q}_{i+1}(0, \cdot)) \right) + \sigma_a (\hat{q}_i(1, \cdot) - \hat{q}_{i+1}(0, \cdot)) \right], \nonumber \\
\Leftrightarrow (\mathbf{H}_i^1)^a &= (a^+ \underline{Q}^{1,1}) \mathbf{q}_i^n + (-a^- \underline{Q}^{1,0}) \mathbf{q}_{i+1}^n. \tag{A9}
\end{align}
$$

• For a right Dirichlet boundary interface

$$
\begin{align}
(H_{N_x}^1)_k^a &= \left[ \theta_k(1), \frac{1}{2} \left( \hat{F}_a(\hat{q}_{N_x}(1, \cdot)) + \hat{F}_a(g_+(\cdot)) \right) + \sigma_a (\hat{q}_{N_x}(1, \cdot) - g_+(\cdot)) \right], \nonumber \\
\Leftrightarrow (\mathbf{H}_{N_x}^1)^a &= (a^+ \underline{Q}^{1,1}) \mathbf{q}_{N_x}^n - a^- \mathbf{d}^1. \tag{A10}
\end{align}
$$

The left diffusive original term.

• For interior interfaces

$$
\begin{align}
(H_i^0)_k^d &= \left[ \theta_k(0), \frac{1}{2} \left( \hat{F}_d(\hat{q}_{i-1}(1, \cdot)) + \hat{F}_d(\hat{q}_i(0, \cdot)) \right) + \sigma_d (\hat{q}_{i-1}(1, \cdot) - \hat{q}_i(0, \cdot)) \right], \nonumber \\
\Leftrightarrow (\mathbf{H}_i^0)^d &= \frac{\lambda}{2\Delta x} \left[ (-\underline{R}^{0,1} + \mu \underline{Q}^{0,1}) \mathbf{q}_{i-1}^n + (-\underline{R}^{0,0} - \mu \underline{Q}^{0,0}) \mathbf{q}_i^n \right]. \tag{A11}
\end{align}
$$

• For a left Dirichlet boundary interface

$$
\begin{align}
(H_1^0)_k^d &= \left[ \theta_k(0), \hat{F}_d(\hat{q}_1(0, \cdot)) + \sigma_d(g_-(\cdot) - \hat{q}_1(0, \cdot)) \right], \nonumber \\
\Leftrightarrow (\mathbf{H}_1^0)^d &= \frac{\lambda}{2\Delta x} \left[ \mu \mathbf{d}^0 + (-2\underline{R}^{0,0} - \mu \underline{Q}^{0,0}) \mathbf{q}_1^n \right]. \tag{A12}
\end{align}
$$

The right diffusive original term.

• For interior interfaces

$$
\begin{align}
(H_i^1)_k^d &= \left[ \theta_k(1), \frac{1}{2} \left( \hat{F}_d(\hat{q}_i(1, \cdot)) + \hat{F}_d(\hat{q}_{i+1}(0, \cdot)) \right) + \sigma_d (\hat{q}_i(1, \cdot) - \hat{q}_{i+1}(0, \cdot)) \right], \nonumber \\
\Leftrightarrow (\mathbf{H}_i^1)^d &= \frac{\lambda}{2\Delta x} \left[ (-\underline{R}^{1,1} + \mu \underline{Q}^{1,1}) \mathbf{q}_i^n + (-\underline{R}^{1,0} - \mu \underline{Q}^{1,0}) \mathbf{q}_{i+1}^n \right]. \tag{A13}
\end{align}
$$

• For a right Dirichlet boundary interface

$$
(H_{N_x}^1)_k^d = \left[ \theta_k(1), \hat{F}_d(\hat{q}_{N_x}(1, \cdot)) + \sigma_d(\hat{q}_{N_x}(1, \cdot) - g_+(\cdot)) \right],
$$
24 High order ADER-IPDG methods for advection-diffusion

• For a left Dirichlet boundary interface

$$
\begin{align}
(H_0^0)_k^a &= \left[ \theta_k(0), \frac{1}{2} \left( \hat{F}_a(g_-(\cdot)) + \hat{F}_a(\hat{q}_1(0, \cdot)) \right) + \sigma_a (g_-(\cdot) - \hat{q}_1(0, \cdot)) \right], \nonumber \\
\Leftrightarrow (\mathbf{H}_0^0)^a &= a^+ \mathbf{d}^0 + (-a^- \underline{Q}^{0,0}) \mathbf{q}_1^n. \tag{A8}
\end{align}
$$

The right advective original term.

• For interior interfaces

$$
\begin{align}
(H_i^1)_k^a &= \left[ \theta_k(1), \frac{1}{2} \left( \hat{F}_a(\hat{q}_i(1, \cdot)) + \hat{F}_a(\hat{q}_{i+1}(0, \cdot)) \right) + \sigma_a (\hat{q}_i(1, \cdot) - \hat{q}_{i+1}(0, \cdot)) \right], \nonumber \\
\Leftrightarrow (\mathbf{H}_i^1)^a &= (a^+ \underline{Q}^{1,1}) \mathbf{q}_i^n + (-a^- \underline{Q}^{1,0}) \mathbf{q}_{i+1}^n. \tag{A9}
\end{align}
$$

• For a right Dirichlet boundary interface

$$
\begin{align}
(H_{N_x}^1)_k^a &= \left[ \theta_k(1), \frac{1}{2} \left( \hat{F}_a(\hat{q}_{N_x}(1, \cdot)) + \hat{F}_a(g_+(\cdot)) \right) + \sigma_a (\hat{q}_{N_x}(1, \cdot) - g_+(\cdot)) \right], \nonumber \\
\Leftrightarrow (\mathbf{H}_{N_x}^1)^a &= (a^+ \underline{Q}^{1,1}) \mathbf{q}_{N_x}^n - a^- \mathbf{d}^1. \tag{A10}
\end{align}
$$

The left diffusive original term.

• For interior interfaces

$$
\begin{align}
(H_i^0)_k^d &= \left[ \theta_k(0), \frac{1}{2} \left( \hat{F}_d(\hat{q}_{i-1}(1, \cdot)) + \hat{F}_d(\hat{q}_i(0, \cdot)) \right) + \sigma_d (\hat{q}_{i-1}(1, \cdot) - \hat{q}_i(0, \cdot)) \right], \nonumber \\
\Leftrightarrow (\mathbf{H}_i^0)^d &= \frac{\lambda}{2\Delta x} \left[ (-\underline{R}^{0,1} + \mu \underline{Q}^{0,1}) \mathbf{q}_{i-1}^n + (-\underline{R}^{0,0} - \mu \underline{Q}^{0,0}) \mathbf{q}_i^n \right]. \tag{A11}
\end{align}
$$

• For a left Dirichlet boundary interface

$$
\begin{align}
(H_1^0)_k^d &= \left[ \theta_k(0), \hat{F}_d(\hat{q}_1(0, \cdot)) + \sigma_d(g_-(\cdot) - \hat{q}_1(0, \cdot)) \right], \nonumber \\
\Leftrightarrow (\mathbf{H}_1^0)^d &= \frac{\lambda}{2\Delta x} \left[ \mu \mathbf{d}^0 + (-2\underline{R}^{0,0} - \mu \underline{Q}^{0,0}) \mathbf{q}_1^n \right]. \tag{A12}
\end{align}
$$

The right diffusive original term.

• For interior interfaces

$$
\begin{align}
(H_i^1)_k^d &= \left[ \theta_k(1), \frac{1}{2} \left( \hat{F}_d(\hat{q}_i(1, \cdot)) + \hat{F}_d(\hat{q}_{i+1}(0, \cdot)) \right) + \sigma_d (\hat{q}_i(1, \cdot) - \hat{q}_{i+1}(0, \cdot)) \right], \nonumber \\
\Leftrightarrow (\mathbf{H}_i^1)^d &= \frac{\lambda}{2\Delta x} \left[ (-\underline{R}^{1,1} + \mu \underline{Q}^{1,1}) \mathbf{q}_i^n + (-\underline{R}^{1,0} - \mu \underline{Q}^{1,0}) \mathbf{q}_{i+1}^n \right]. \tag{A13}
\end{align}
$$

• For a right Dirichlet boundary interface

$$
(H_{N_x}^1)_k^d = \left[ \theta_k(1), \hat{F}_d(\hat{q}_{N_x}(1, \cdot)) + \sigma_d(\hat{q}_{N_x}(1, \cdot) - g_+(\cdot)) \right],
$$
$$ \Leftrightarrow (\mathbf{H}_{N_x}^1)^d = \frac{\lambda}{2\Delta x} [(-2\underline{R}^{1,1} + \mu\underline{Q}^{1,1}) \mathbf{q}_{N_x}^n - \mu\mathbf{d}^1]. \quad (A14) $$

The left diffusive extra term.

* For interior interfaces

$$ \begin{aligned} (\bar{H}_i^0)_k^d &= \left[-\frac{\lambda}{\Delta x} \partial_\xi \theta_k(0), \frac{1}{2}(\hat{q}_i(0, \cdot) - \hat{q}_{i-1}(1, \cdot))\right], \\ \Leftrightarrow (\bar{H}_i^0)^d &= \frac{\lambda}{2\Delta x} \left[(\underline{W}^{0,1}) \mathbf{q}_{i-1}^n + (-\underline{W}^{0,0}) \mathbf{q}_i^n\right]. \end{aligned} \quad (A15) $$

* For a left Dirichlet boundary interface

$$ \begin{aligned} (\bar{H}_1^0)_k^d &= \left[-\frac{\lambda}{\Delta x} \partial_\xi \theta_k(0), (\hat{q}_1(0, \cdot) - g_-(\cdot))\right], \\ \Leftrightarrow (\bar{H}_1^0)^d &= \frac{\lambda}{2\Delta x} \left[2\mathbf{e}^0 + (-2\underline{W}^{0,0}) \mathbf{q}_1^n\right]. \end{aligned} \quad (A16) $$

The right diffusive extra term.

* For interior interfaces

$$ \begin{aligned} (\bar{H}_i^1)_k^d &= \left[-\frac{\lambda}{\Delta x} \partial_\xi \theta_k(1), \frac{1}{2}(\hat{q}_i(1, \cdot) - \hat{q}_{i+1}(0, \cdot))\right], \\ \Leftrightarrow (\bar{H}_i^1)^d &= \frac{\lambda}{2\Delta x} \left[(-\underline{W}^{1,1}) \mathbf{q}_i^n + (\underline{W}^{1,0}) \mathbf{q}_{i+1}^n\right]. \end{aligned} \quad (A17) $$

* For a right Dirichlet boundary interface

$$ \begin{aligned} (\bar{H}_{N_x}^1)_k^d &= \left[-\frac{\lambda}{\Delta x} \partial_\xi \theta_k(1), (\hat{q}_{N_x}(1, \cdot) - g^+(\cdot))\right], \\ \Leftrightarrow (\bar{H}_{N_x}^1)^d &= \frac{\lambda}{2\Delta x} \left[(-2\underline{W}^{1,1}) \mathbf{q}_{N_x}^n + 2\mathbf{e}^1\right]. \end{aligned} \quad (A18) $$

Substituting these different quantities that have been expressed in terms of the unknown DOF in each cell provides a complete and fully explicit ADER-IPDG scheme.

## Declarations

On behalf of all authors, the corresponding author states that there is no conflict of interest.
$$ \Leftrightarrow (\mathbf{H}_{N_x}^1)^d = \frac{\lambda}{2\Delta x} [(-2\underline{R}^{1,1} + \mu\underline{Q}^{1,1}) \mathbf{q}_{N_x}^n - \mu\mathbf{d}^1]. \quad (A14) $$

The left diffusive extra term.

* For interior interfaces

$$ \begin{aligned} (\bar{H}_i^0)_k^d &= \left[-\frac{\lambda}{\Delta x} \partial_\xi \theta_k(0), \frac{1}{2} (\hat{q}_i(0, \cdot) - \hat{q}_{i-1}(1, \cdot))\right], \\ \Leftrightarrow (\bar{H}_i^0)^d &= \frac{\lambda}{2\Delta x} \left[(\underline{W}^{0,1}) \mathbf{q}_{i-1}^n + (-\underline{W}^{0,0}) \mathbf{q}_i^n\right]. \end{aligned} \quad (A15) $$

* For a left Dirichlet boundary interface

$$ \begin{aligned} (\bar{H}_1^0)_k^d &= \left[-\frac{\lambda}{\Delta x} \partial_\xi \theta_k(0), (\hat{q}_1(0, \cdot) - g_-(\cdot))\right], \\ \Leftrightarrow (\bar{H}_1^0)^d &= \frac{\lambda}{2\Delta x} \left[2\mathbf{e}^0 + (-2\underline{W}^{0,0}) \mathbf{q}_1^n\right]. \end{aligned} \quad (A16) $$

The right diffusive extra term.

* For interior interfaces

$$ \begin{aligned} (\bar{H}_i^1)_k^d &= \left[-\frac{\lambda}{\Delta x} \partial_\xi \theta_k(1), \frac{1}{2} (\hat{q}_i(1, \cdot) - \hat{q}_{i+1}(0, \cdot))\right], \\ \Leftrightarrow (\bar{H}_i^1)^d &= \frac{\lambda}{2\Delta x} \left[(-\underline{W}^{1,1}) \mathbf{q}_i^n + (\underline{W}^{1,0}) \mathbf{q}_{i+1}^n\right]. \end{aligned} \quad (A17) $$

* For a right Dirichlet boundary interface

$$ \begin{aligned} (\bar{H}_{N_x}^1)_k^d &= \left[-\frac{\lambda}{\Delta x} \partial_\xi \theta_k(1), (\hat{q}_{N_x}(1, \cdot) - g^+(\cdot))\right], \\ \Leftrightarrow (\bar{H}_{N_x}^1)^d &= \frac{\lambda}{2\Delta x} \left[(-2\underline{W}^{1,1}) \mathbf{q}_{N_x}^n + 2\mathbf{e}^1\right]. \end{aligned} \quad (A18) $$

Substituting these different quantities that have been expressed in terms of the unknown DOF in each cell provides a complete and fully explicit ADER-IPDG scheme.

## Declarations

On behalf of all authors, the corresponding author states that there is no conflict of interest.
References

[1] Abgrall, R.: Residual distribution schemes: Current status and future trends. Computers & Fluids 35(7), 641–669 (2006). https://doi.org/10.1016/j.compfluid.2005.01.007

[2] Ricchiuto, M., Abgrall, R.: Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case. Journal of Computational Physics 229(16), 5653–5691 (2010). https://doi.org/10.1016/j.jcp.2010.04.002

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

[4] 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

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

[6] 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

[7] 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

[8] 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

[9] 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

[10] Bergmann, M., Carlino, M.G., Iollo, A.: Second Order ADER Scheme for
References

[1] Abgrall, R.: Residual distribution schemes: Current status and future trends. Computers & Fluids 35(7), 641–669 (2006). https://doi.org/10.1016/j.compfluid.2005.01.007

[2] Ricchiuto, M., Abgrall, R.: Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case. Journal of Computational Physics 229(16), 5653–5691 (2010). https://doi.org/10.1016/j.jcp.2010.04.002

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

[4] 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

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

[6] 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

[7] 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

[8] 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

[9] 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

[10] 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

[11] 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

[12] 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

[13] 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

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

[15] 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

[16] 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

[17] 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

[18] Dumbser, M., Munz, C.-D.: Building Blocks for Arbitrary High Order Discontinuous Galerkin Schemes. Journal of Scientific Computing **27**(1-3), 215–230 (2006). https://doi.org/10.1007/s10915-005-9025-0

[19] 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
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

[11] 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

[12] 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

[13] 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

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

[15] 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

[16] 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

[17] 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

[18] Dumbser, M., Munz, C.-D.: Building Blocks for Arbitrary High Order Discontinuous Galerkin Schemes. Journal of Scientific Computing **27**(1-3), 215–230 (2006). https://doi.org/10.1007/s10915-005-9025-0

[19] 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
[20] 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

[21] 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

[22] Dumbser, M., Balsara, D.S., Toro, E.F., Munz, C.-D.: A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes. Journal of Computational Physics 227(18), 8209-8253 (2008). https://doi.org/10.1016/j.jcp.2008.05.025

[23] 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

[24] 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

[25] 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

[26] 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

[27] Zahran, Y.H.: Central ADER schemes for hyperbolic conservation laws. Journal of Mathematical Analysis and Applications 346(1), 120-140 (2008). https://doi.org/10.1016/j.jmaa.2008.05.032

[28] 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
[20] 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

[21] 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

[22] Dumbser, M., Balsara, D.S., Toro, E.F., Munz, C.-D.: A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes. Journal of Computational Physics **227**(18), 8209–8253 (2008). https://doi.org/10.1016/j.jcp.2008.05.025

[23] 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

[24] 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

[25] 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

[26] 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

[27] Zahran, Y.H.: Central ADER schemes for hyperbolic conservation laws. Journal of Mathematical Analysis and Applications **346**(1), 120–140 (2008). https://doi.org/10.1016/j.jmaa.2008.05.032

[28] 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