# Discontinuous Galerkin methods

Bernardo Cockburn

► To cite this version:

Bernardo Cockburn. Discontinuous Galerkin methods. Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 2003, 83 (11), pp.731-754. 10.1002/zamm.200310088 . hal-01352444

HAL Id: hal-01352444

https://hal.science/hal-01352444

Submitted on 8 Aug 2016

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_1.png" /></figure>
# Discontinuous Galerkin methods

Bernardo Cockburn

University of Minnesota, School of Mathematics, Minneapolis, MN 55455, USA

**Key words** discontinuous Galerkin methods, finite element methods

This paper is a short essay on discontinuous Galerkin methods intended for a very wide audience. We present the discontinuous Galerkin methods and describe and discuss their main features. Since the methods use completely discontinuous approximations, they produce mass matrices that are block-diagonal. This renders the methods highly parallelizable when applied to hyperbolic problems. Another consequence of the use of discontinuous approximations is that these methods can easily handle irregular meshes with hanging nodes and approximations that have polynomials of different degrees in different elements. They are thus ideal for use with adaptive algorithms. Moreover, the methods are locally conservative (a property highly valued by the computational fluid dynamics community) and, in spite of providing discontinuous approximations, stable, and high-order accurate. Even more, when applied to non-linear hyperbolic problems, the discontinuous Galerkin methods are able to capture highly complex solutions presenting discontinuities with high resolution. In this paper, we concentrate on the exposition of the ideas behind the devising of these methods as well as on the mechanisms that allow them to perform so well in such a variety of problems.

## 1 Introduction

The discontinuous Galerkin (DG) methods are locally conservative, stable, and high-order accurate methods which can easily handle complex geometries, irregular meshes with hanging nodes, and approximations that have polynomials of different degrees in different elements. These properties, which render them ideal to be used with hp-adaptive strategies, not only have brought these methods into the main stream of computational fluid dynamics, for example, in gas dynamics [11, 14, 39], compressible [10, 64–66] and incompressible [13, 32, 33] flows, turbomachinery [12], magneto-hydrodynamics [77], granular flows [52, 53], semiconductor device simulation [24, 25], viscoplastic crack growth and chemical transport [21], viscoelasticity [5, 8, 51], and transport of contaminant in porous media [1, 28, 29, 41], but have also prompted their application to a wide variety of problems for which they were not originally intended like, for example, Hamilton-Jacobi equations [59, 60, 62], second-order elliptic problems [4, 6, 20, 22, 31, 38, 67, 70], elasticity [46, 54], and Korteweg-deVries equations [72, 73].

An introduction to DG methods can be found in the short monograph [26]. A history of their development and the state of the art up to 1999 can be found in [34]. Finally, a fairly complete and updated review is given in [40]. This paper is a short essay on DG methods which differs from the above mentioned references in that it is intended for a wider audience and focuses on the exposition of the ideas behind the devising of these methods as well as in the mechanisms that allow them to perform so well in such a variety of problems.

Let us briefly carry out this program in what is perhaps the simplest situation, namely, that of approximating the solution of an ordinary differential equation (ODE). So, consider the following model initial-value problem:

$$ \frac{d}{dt} u(t) = f(t) u(t), \quad t \in (0, T), \quad u(0) = u_0, \qquad (1) $$

and suppose that we want to compute an approximation $u_h$ to $u$ on the interval $(0, T)$ by using a DG method. To do that, we first find a partition of the interval $(0, T)$, $\{t^n\}_{n=0}^N$ and set $I^n = (t^n, t^{n+1})$ for $n = 0, \dots, N-1$. Then we look for a function $u_h$ which, on the interval $I^n$, is the polynomial of degree at most $k^n$ determined by requiring that

$$ - \int_{I^n} u_h(s) \frac{d}{dt} v(s) ds + \hat{u}_h v|_{t^n}^{t^{n+1}} = \int_{I^n} f(s) u(s) v(s) ds, \qquad (2) $$

for all polynomials $v$ of degree at most $k^n$. To complete the definition of the DG method, we still need to define the quantity $\hat{u}_h$. Note that the method establishes a link between the values of $u_h$ in different intervals only through $\hat{u}_h$. Since for the ODE, the information travels “from the past into the future”, it is reasonable to take $\hat{u}_h$ as follows:
$$ \hat{u}_h(t^n) = \begin{cases} u_0, & \text{if } t^n = 0, \\ \lim_{\epsilon \downarrow 0} u_h(t^n - \epsilon) & \text{otherwise.} \end{cases} \tag{3} $$

This completes the definition of the DG method.

In this simple example, we already see the main components of the method, namely,

(i) The use of discontinuous approximations $u_h$,

(ii) The enforcing of the ODE on each interval by means of a Galerkin weak formulation, and

(iii) The introduction and suitable definition of the so-called numerical trace $\hat{u}_h$.

The choice of the numerical trace is perhaps the most delicate and crucial aspect of the definition of the DG methods as it can affect their consistency, stability, and even accuracy. The simple choice we have made is quite natural for this case and gives rise to a very good method; however, other choices can also produce excellent results. Next, we address the question of how to choose the numerical trace $\hat{u}_h$.

Let us begin with the problem of the consistency of the DG method. As it is typical for most finite element methods, the method is said to be consistent if we can replace the approximate solution $u_h$ by the exact solution $u$ in the weak formulation (2). We can immediately see that this is true if and only if $\hat{u} = u$.

Next, let us consider the more subtle issue of the stability of the method. Our strategy is to begin by obtaining a stability property for the ODE (1) which we will then try to enforce for the DG method (2) by a suitable definition of the numerical trace $\hat{u}_h$. If we multiply the ODE by $u$ and integrate over $(0, T)$, we get the equality

$$ \frac{1}{2} u^2(T) - \frac{1}{2} u_0^2 = \int_0^T f(s) u^2(s) ds, $$

from which the $L^\infty$-stability of the solution follows. The result we have to obtain for the DG method is a similar equality.

To do that, it is enough to set $v = u_h$ in the weak formulation (2), integrate by parts and add over $n$. We get

$$ \sum_{n=0}^{N-1} \left( -\frac{1}{2} u_h^2 + \hat{u}_h u_h \right) \Big|_{t^n}^{t^{n+1}} = \frac{1}{2} u_h^2(T^{-}) + \Theta_h(T) - \frac{1}{2} u_0^2 = \int_0^T f(s) u_h^2(s) ds, $$

where

$$ \Theta_h(T) = -\frac{1}{2} u_h^2(T^{-}) + \sum_{n=0}^{N-1} \left( -\frac{1}{2} u_h^2 + \hat{u}_h u_h \right) \Big|_{t^n}^{t^{n+1}} + \frac{1}{2} u_0^2. $$

Note that if $\Theta_h(T)$ were a non-negative quantity, the above equality would imply the stability of the DG method. In other words, the stability of the DG method is guaranteed if we can define the numerical trace $\hat{u}_h$ so that $\Theta_h(T) \ge 0$. Setting

$$ u_h(t) = u_0, \quad t < 0, $$

and using the notation

$$ \{u_h\} = \frac{1}{2}(u_h^{-} + u_h^{+}), \quad \llbracket u_h \rrbracket = u_h^{-} - u_h^{+}, \quad u_h^{\pm}(t) = \lim_{\epsilon \downarrow 0} u_h(t \pm \epsilon), $$

we rewrite $\Theta_h(T)$ as follows:

$$ \begin{aligned} \Theta_h(T) &= -\frac{1}{2}u_h^2(T^{-}) + \left(-\frac{1}{2}u_h^2(T^{-}) + \hat{u}_h(T)u_h(T^{-})\right) + \sum_{n=1}^{N-1}\left(-\frac{1}{2}\llbracket u_h^2 \rrbracket + \hat{u}_h\llbracket u_h \rrbracket\right)(t^n) \\ &\quad -\left(-\frac{1}{2}u_h^2(0^+) + \hat{u}_h(0)u_h(0^+)\right) + \frac{1}{2}u_0^2 \\ &= (\hat{u}_h(T) - u_h(T^{-}))u_h(T^{-}) + \sum_{n=1}^{N-1}\left((\hat{u}_h - \{u_h\})\llbracket u_h \rrbracket\right)(t^n) - (\hat{u}_h(0) - u_0)u_h(0^+) + \frac{1}{2}\llbracket u_h \rrbracket^2(0), \end{aligned} $$

where we used the simple identity

$$ \llbracket u_h^2 \rrbracket = 2 \{u_h\} \llbracket u_h \rrbracket. $$

It is now clear that if we take

$$ \hat{u}_h(t^n) = \begin{cases} u_0, & \text{if } t^n = 0, \\ (\{u_h\} + C^n \llbracket u_h \rrbracket)(t^n), & \text{if } t^n \in (0, T), \\ u_h(T^{-}), & \text{if } t^n = T, \end{cases} $$
where $C^n \ge 0$, we would have, setting $C^0 = 1/2$,

$$ \Theta_h(T) = \sum_{n=0}^{N-1} C^n \llbracket u_h \rrbracket^2(t^n) \ge 0, $$

just as we wanted.

We thus see that it is possible to define the numerical trace $\hat{u}_h$ to enforce the stability of the method. Note that the choice $C^n \equiv 1/2$ corresponds to the numerical trace we chose at the beginning, namely, $\hat{u}_h(t) = u_h(t^-)$. Note that the resulting DG methods are not only stable but consistent for all values of $C^n \ge 0$, as the condition $\hat{u} = u$ is satisfied. Moreover, in our search for stability, we found, in a very natural way, that the numerical trace $\hat{u}_h(t)$ can only depend on both traces of $u_h$ at $t$, that is, on $u_h(t^-)$ and on $u_h(t^+)$.

Finally, let us address the issue of how the choice of the coefficients $C^n$ can affect the accuracy of the method. It can be proven that if we take $C^n = 1/2$, the order of the method at the points $t^n$ is $2k + 1$, [63]. However, if we simply take $C^n = 0$, the order is $2k + 2$, [42]. Of course, in this case we must sacrifice the ability of solving for $u_h$ interval by interval since, if $C^n \ne 1/2$, we must solve for $u_h$ in the whole computational domain $(0, T)$ at once. Thus, if we insist in having a DG method that is consistent, stable and can solve for the approximate solution interval by interval, we must take $C^n \equiv 1/2$!

Next, we want to emphasize three important properties of the DG methods that do carry over to the multi-dimensional case and to all types of problems. The first is that the approximate solution of the DG methods does not have to satisfy any inter-element continuity constraint. As a consequence, the method can be highly parallelizable (when dealing with time-dependent hyperbolic problems) and can easily handle different types of approximations in different elements which renders it ideal for use with hp-adaptivity.

The second is that the DG methods are locally conservative. This is a reflection of the fact that the method enforces the equation element-by-element and of the use of the numerical trace. In our simple setting, this property reads

$$ \hat{u}_h|_{t^n}^{t^{n+1}} = \int_{I_n} f(s) u(s) ds, $$

and is obtained by simply taking $v \equiv 1$ in the weak formulation (2). This a much valued property in computational fluid dynamics.

The third property, which has not been properly discussed in the available literature on DG methods, is the strong relation between the residuals of $u_h$ inside the intervals and its jumps across inter-interval boundaries. To uncover it, let us integrate by parts in (2) to get

$$ \int_{I_n} \frac{d}{dt} u_h(s) v(s) ds + (\hat{u}_h - u_h) v|_{t^n}^{t^{n+1}} = \int_{I_n} f(s) u(s) v(s) ds, $$

or, equivalently,

$$ \int_{I_n} R(s) v(s) ds = (u_h - \hat{u}_h) v|_{t^n}^{t^{n+1}}, $$

where $R$ denotes the residual $(\frac{d}{dt}u_h - f u_h)$. If we now take $v = 1$ and use the definition of the numerical trace $\hat{u}_h$, we obtain

$$ \int_{I_n} R(s) ds = \llbracket u_h \rrbracket(t^n). $$

In other words, the jump of $u_h$ at $t^n$, $\llbracket u_h \rrbracket(t^n)$, is nothing but the integral of the residual over the interval $I^n$.

This means that if the ODE has been very well approximated in the interval $I^n$, the jump $\llbracket u \rrbracket(t^n)$ is very small. On the other hand, if the ODE was not well approximated therein, the jump $\llbracket u_h \rrbracket(t^n)$ is big and the DG method becomes more dissipative as can be seen directly from the stability identity

$$ \frac{1}{2} u_h^2(T^-) + \frac{1}{2} \sum_{n=0}^{N-1} \llbracket u_h \rrbracket^2(t^n) - \frac{1}{2} u_0^2 = \int_0^T f(s) u_h^2(s) ds. $$

Thus, roughly speaking, the jumps $\llbracket u_h \rrbracket$ act as dampers that stabilize the DG method whenever the ride becomes bumpy and the ODE cannot be well approximated. This renders the method very stable without degrading its accuracy. In multi-dimensions, the residual is a more complicated linear functional of the jumps, but the dissipative mechanism just described remains essentially the same. This mechanism is shared with the so-called stabilized finite element methods for convection and for second-order problems.

In what follows, we consider DG methods for a variety of problems. In each of them, we study the properties of consistency, stability, and accuracy of the method, and discuss the relevance of the properties of local conservativity and stabilization by the
jumps. We also comment on some computational issues and establish various links with other numerical schemes. In Sect. 2, devoted to linear problems, we consider (i) the transport equation, (ii) the wave equation, (iii) the Poisson equation, and (iv) the Oseen system. In Sect. 3, we consider the more difficult non-linear hyperbolic problems. We end the paper in Sect. 4, with some concluding remarks.

## 2 Linear problems

### 2.1 The transport equation

In this sub-section, we consider DG methods for the transport equation

$$ 
\begin{aligned}
u_t + \nabla \cdot (a \, u) &= 0 \quad \text{in } \mathbb{R}^d \times (0, T), \\
u(t=0) &= u_0 \quad \text{on } \mathbb{R}^d.
\end{aligned}
$$

We consider only the discretization of this equation in space; the full discretization will be studied when we deal with non-linear hyperbolic problems in Sect. 3.

The objective here is to examine three properties that are especially relevant in this case. The first is that when the DG methods are strongly related to classical finite volume methods like the up-winding and the Lax-Friedrichs methods. The second is that when polynomials of high degree are used, the DG method can achieve high-order accuracy while remaining highly parallelizable. The third is that the artificial viscosity of the method is given by the size of the jumps which, in turn, are associated with the residual inside the elements. As a consequence, as the polynomial degree of its approximate solution increases, the artificial viscosity diminishes even in the presence of discontinuities.

**The DG methods.** To discretize the transport equation in space by using a DG method, we first triangulate the domain $\mathbb{R}^d$; let us denote by $\mathcal{T}_h$ such triangulation. We then seek a discontinuous approximate solution $u_h$ which, in each element $K$ of the triangulation $\mathcal{T}_h$, belongs to the space $V(K)$. There is no restriction in how to choose the space $V(K)$, although a typical choice is the space of polynomials of degree at most $k$, $P^k(K)$. We determine the approximate solution on the element $K$ by weakly enforcing the transport equation as follows:

$$ \int_K (u_h)_t v - \int_K a u_h \cdot \nabla v + \int_{\partial K} \widehat{a u_h} \cdot n \, v \, ds = 0, \quad (4) $$

for all $v \in V(K)$. To complete the definition of the DG method, it only remains to define the numerical trace $\widehat{a u_h}$.

To do that, we proceed as in the ODE considered in the Introduction and begin by obtaining a stability result for the problem under consideration. We thus multiply the transport equation by $u$, and integrate over the space and time to get

$$ \frac{1}{2} \int_{\mathbb{R}^d} u^2(x, T) \, dx + \frac{1}{2} \int_0^T \int_{\mathbb{R}^d} \nabla \cdot a(x) u^2(x, t) \, dx \, dt = \frac{1}{2} \int_{\mathbb{R}^d} u_0^2(x) \, dx. $$

From this equation, a stability result immediately follows if we assume, for example, that $-\nabla \cdot a \le L$.

Now, let us mimic the above procedure for the DG method under consideration. Taking $v = u_h$ in the weak formulation defining the approximate solution $u_h$, and adding on the elements $K$, we get

$$ \frac{1}{2} \int_{\mathbb{R}^d} u_h^2(x, T) \, dx + \frac{1}{2} \int_0^T \int_{\mathbb{R}^d} \nabla \cdot a(x) u_h^2(x, t) \, dx \, dt + \int_0^T \Theta_h(t) \, dt = \frac{1}{2} \int_{\mathbb{R}^d} u_{h,0}^2(x) \, dx, $$

where

$$ \Theta_h(t) = \sum_{K \in \mathcal{T}_h} \left( -\frac{1}{2} \int_K \nabla \cdot (a u_h)(x, t) \, dx + \int_{\partial K} \widehat{a u_h}(x, t) \cdot n \, u_h(x, t) \, ds \right). $$

Next, we investigate if it is possible to define the numerical trace $\widehat{a u_h}$ in such a way as to render $\Theta_h$ non-negative.

At this point, it is convenient to introduce the following notation. Let $x$ be a point on the set $e = \overline{\partial K^+} \cap \overline{\partial K^-}$ and let $n_\pm$ denote the unit outward normal to $\partial K^\pm$ at the point $x$. Let $u_h^\pm(x)$ denote the value $\lim_{\epsilon \downarrow 0} u_h(x - \epsilon n^\pm)$ and set

$$ \{u_h\} = \frac{1}{2} (u_h^+ + u_h^-), \quad \llbracket u_h \rrbracket = u_h^- n^- + u_h^+ n^+. $$

Finally, let $\mathcal{E}_h$ denote the set of sets $e = \overline{\partial K^+} \cap \overline{\partial K^-}$ for all $K^+$ and $K^- \in \mathcal{T}_h$.
We are now ready to rewrite $\Theta_h$ in a suitable way. Indeed, dropping the argument $t$, we have

$$
\begin{align*}
\Theta_h &= \sum_{K \in \mathcal{T}_h} \int_{\partial K} \left( \widehat{a u_h} \cdot n \, u_h - \frac{1}{2} \widehat{a} \, u_h^2 \cdot n \right) \, ds \\
&= \sum_{e \in \mathcal{E}_h} \int_e \left[ \widehat{a u_h} \, u_h - \frac{1}{2} \widehat{a} \, u_h^2 \right] \, ds \\
&= \sum_{e \in \mathcal{E}_h} \int_e \left( \widehat{a u_h} \left[ \llbracket u_h \rrbracket \right] - \frac{1}{2} \left[ \llbracket \widehat{a} \, u_h^2 \rrbracket \right] \right) \, ds \\
&= \sum_{e \in \mathcal{E}_h} \int_e \left( \widehat{a u_h} - \widehat{a} \{ u_h \} \right) \cdot \left[ \llbracket u_h \rrbracket \right] \, ds
\end{align*}
$$

since $\frac{1}{2}[\llbracket u_h^2 \rrbracket] = \{u_h\} [\llbracket u_h \rrbracket]$. Thus, if we take

$$
\widehat{a u_h} = a \{u_h\} + C[\llbracket u_h \rrbracket],
$$

we get

$$
\Theta_h = \sum_{e \in \mathcal{E}_h} \int_e C[\llbracket u_h \rrbracket] \cdot [\llbracket u_h \rrbracket] \, ds \ge 0,
$$

if $C$ is a non-negative definite matrix. This completes the definition of the DG space discretization and shows that it is possible
to define the numerical trace as to render the method consistent and stable.

**Examples of the DG methods.** There are two main examples of DG methods in this case. The first uses the following choice for its artificial viscosity coefficient: $C = \frac{1}{2} |\mathbf{a} \cdot \mathbf{n}| \text{ Id}$. This implies that the numerical trace is

$$
\widehat{a u_h}(\mathbf{x}) = \mathbf{a} \lim_{\epsilon \downarrow 0} u_h(\mathbf{x} - \epsilon \mathbf{a}),
$$

which is nothing but the classical *up-winding numerical flux*.

The second example is when we take $C = \frac{1}{2} |\mathbf{a}| \text{Id}$. For this choice, we have

$$
\widehat{a u_h} = a \{u_h\} + \frac{1}{2} |a| [\llbracket u_h \rrbracket],
$$

which is the so-called *Lax-Friedrichs numerical flux*.

**Some properties.** (i) From the two examples above, we see that the DG methods are strongly related to finite volume methods.
Indeed, the method of lines, that is, the discretization in space, for the up-winding scheme and the local Lax-Friedrichs scheme
*coincide* with the corresponding DG method under consideration when the local space $V(K)$ is taken to be the space of constant
functions.

The DG methods, like finite volume methods, can easily handle complex computational domains. Also like finite volume methods, they have the property of being locally conservative, that is,

$$
\int_K (u_h)_t dx + \int_{\partial K} \widehat{a} u_h \cdot n ds = 0,
$$

provided that the local space $V(K)$ contains the constant functions. Indeed, this property is easily obtained from eq. (4) by
simply taking the test function $v$ to be a constant.

Unlike finite volume methods, however, DG achieve with ease high-order accuracy. Indeed, a theoretical order of convergence of $k+1/2$ can be proven simply by requiring that the local spaces $V(K)$ contain all polynomials of degree at most $k$. Moreover, this is achieved while keeping a high degree of locality since to evolve the degrees of freedom of the approximate solution $u_h$ in an element, only the degrees of freedom of $u_h$ in the immediate neighbors are involved.

(ii) This last property and the fact that the mass matrix is block diagonal, and hence easily invertible, renders the DG methods extremely parallelizable when they are discretized in time by, for example, an explicit Runge-Kutta method.

We give an example taken from [18]. In Table 1 below, we see the solution time and total execution time for the two-
dimensional linear problem

$$
u_t + u_x + u_y = 0, \quad \text{in } (-\pi, \pi)^2 \times (0, T),
$$
<table><caption>**Table 1** Scaled parallel efficiency. Solution times (without I/O) and total execution times measured on the nCUBE/2.</caption><thead><tr><th>Number of processors</th><th>Work (W)</th><th>Solution time (secs.)</th><th>Solution parallel efficiency</th><th>Total time (secs.)</th><th>Total parallel efficiency</th></tr></thead><tbody><tr><td>1</td><td>18,432</td><td>926.92</td><td>-</td><td>927.16</td><td>-</td></tr><tr><td>2</td><td>36,864</td><td>927.06</td><td>99.98 %</td><td>927.31</td><td>99.98 %</td></tr><tr><td>4</td><td>73,728</td><td>927.13</td><td>99.97 %</td><td>927.45</td><td>99.96 %</td></tr><tr><td>8</td><td>147,456</td><td>927.17</td><td>99.97 %</td><td>927.58</td><td>99.95 %</td></tr><tr><td>16</td><td>294,912</td><td>927.38</td><td>99.95 %</td><td>928.13</td><td>99.89 %</td></tr><tr><td>32</td><td>589,824</td><td>927.89</td><td>99.89 %</td><td>929.90</td><td>99.70 %</td></tr><tr><td>64</td><td>1,179,648</td><td>928.63</td><td>99.81 %</td><td>931.28</td><td>99.55 %</td></tr><tr><td>128</td><td>2,359,296</td><td>930.14</td><td>99.65 %</td><td>937.67</td><td>98.88 %</td></tr><tr><td>256</td><td>4,718,592</td><td>933.97</td><td>99.24 %</td><td>950.25</td><td>97.57 %</td></tr></tbody></table>

with initial condition $u(x, y, 0) = \sin(\pi x) \sin(\pi y)$ and periodic boundary conditions. We can see that as we progress from 1 to 256 processors while keeping the work per processor constant, a remarkable parallel efficiency is achieved. Polynomials of degree two and an explicit time-marching scheme were used.

(iii) Next, we show that the dissipation of the DG methods is given by the jumps of their approximate solution. This can be immediately seen when we realize that the DG method has a higher rate of dissipation of the energy, which in this case is nothing but the square of the $L^2$-norm, than the exact solution of the transport equation. The extra rate of dissipation for the DG method is given by

$$ \Theta_h = \sum_{e \in \mathcal{E}_h} \int_e C \llbracket u_h \rrbracket \cdot \llbracket u_h \rrbracket \, ds. $$

In the literature for monotone finite difference schemes for hyperbolic problems, the above term, when $C = \nu$ Id, is introduced by what could be considered to be a term modeling a viscosity effect, with $\nu$ being the viscosity coefficient, artificially inserted to render the scheme stable. That is why it is also called artificial viscosity. We thus see that the artificial viscosity of the DG method solely depends on the jumps of their approximate solution.

Moreover, the jumps and the local residual $(u_h)_t + \nabla \cdot (a u_h)$, which we denote by $R$, are strongly related. To see this, it is enough to carry out a simple integration by parts in the definition of the approximate solution:

$$ \int_K R v = \int_{\partial K} (a u_h \cdot n v - \tilde{a} \tilde{u}_h \cdot n v) \, dx. $$

Note that in the case of the up-winding flux, we get

$$ \int_K R v = \int_{\partial K^{-}} a \cdot \llbracket u_h \rrbracket \, ds, \quad \text{where} \quad \partial K^{-} = \{x \in \partial K : a(x) \cdot n(x) \le 0\}. $$

In other words, the residual of $u_h$ in $K$ is linearly related to the jump of $u_h$ on its inflow boundary $\partial K^{-}$. A similar, but more complicated relation holds for general DG methods.

This indicates that the artificial viscosity of the method depends on the polynomial degree of the approximate solution. If the solution is very smooth, as the polynomial degree increases, we expect to have a smaller residual, hence a smaller jump and hence a smaller artificial viscosity. This is indeed what happens as the method can be proven to be of order $k+1/2$ uniformly in time in the $L^2$-norm when polynomials of degree $k$ are used and the exact solution is smooth. Moreover, the artificial viscosity or, better, the dissipation of the DG methods, does decrease as the polynomial degree increases. This can be seen in Fig. 1 where we compare the exact solution of the transport equation

$$ u_t + u_x = 0, \quad \text{in } (0,1) \times (0,T = 100), $$

with initial condition

$$ u(x, 0) = \begin{cases} 1, & \text{if } x \in (.4, .6), \\ 0, & \text{otherwise,} \end{cases} $$

and periodic boundary conditions, with the DG approximations obtained with polynomial approximations of degree $k = 0, 1, 2$. To march in time, a Runge-Kutta method of order $k+1$ was used.
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Effect of the polynomial degree on the dissipation of the DG method. The exact solution $u$ (solid line) is contrasted against the approximate solution obtained on a mesh of 160 elements with piecewise-constant (dash-point line), piecewise-linear (dotted line) and piecewise-quadratic (dashed line) approximations.</figcaption></figure>

## 2.2 The wave equation

In this subsection, we want to show how to extend the definition of DG methods to linear hyperbolic systems. To do that, we consider DG methods for the wave equation

$$u_{tt} - c^2 \Delta u = 0, \quad \text{in } \mathbb{R}^d \times (0, T),$$

where, for simplicity, we take speed $c$ to be constant, and rewrite it as a first-order hyperbolic problem, namely, as

$$U_t + \nabla \cdot F(U) = 0, \quad \text{in } \mathbb{R}^d \times (0, T),$$

where

$$U = \begin{pmatrix} q_1 \\ q_2 \\ \vdots \\ q_d \\ u \end{pmatrix}, \quad F(U) = -c \begin{pmatrix} u & 0 & \dots & 0 \\ 0 & u & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & u \\ q_1 & q_2 & \dots & q_d \end{pmatrix}.$$

We want to show that the properties of local conservativity, stability, dissipation, and parallelizability of DG methods for the transport equation also hold for the above model hyperbolic problem. We also want to establish a natural link between the DG discretization of the wave equation and the DG discretization of the Laplacian operator.

**The DG methods.** To discretize the wave equation in space by using a DG method, we proceed as follows. After triangulating the domain $\mathbb{R}^d$, for each element $K$ of the triangulation $\mathcal{T}_h$, we choose the local space $\mathcal{U}(K)$ which $U_h|_K$ belongs to. Note, once again, that there is no constraint on how to take the space $\mathcal{U}(K)$; a typical choice, however, is $\mathcal{U}(K) = P^k(K) \times \dots \times P^k(K)$. Then, we determine the approximate solution on the element $K$ by weakly enforcing the conservation law as follows:

$$\int_K (U_{hi})_t V_i - \int_K F_{ij}(U_h) V_{i,j} + \int_{\partial K} \widehat{F_{ij}} n_j V_i dx = 0,$$

for all $V \in \mathcal{U}(K)$. Here, we adopt Einstein summation convention; note that $V_{i,j}$ denotes $\partial_{x_j} V_i$. To complete the definition of the DG method, it only remains to define the numerical trace $\widehat{F_{ij}}$.

To do that, we begin by obtaining the stability result for the problem under consideration. We multiply the hyperbolic system by $U$, integrate over the space and time to get

$$\frac{1}{2} \int_{\mathbb{R}^d} U^2(\boldsymbol{x}, T) dx = \frac{1}{2} \int_{\mathbb{R}^d} U^2(\boldsymbol{x}, 0) dx,$$

where $U^2 = U_i U_i$. Now, let us mimic the above procedure for the DG method under consideration. Taking $V = U_h$ and adding on the elements $K$, we get

$$\frac{1}{2} \int_{\mathbb{R}^d} U_h^2(\boldsymbol{x}, T) dx + \int_0^T \Theta_h(t) dt = \frac{1}{2} \int_{\mathbb{R}^d} U_h^2(\boldsymbol{x}, 0) dx,$$
where

$$ \Theta_h(t) = \sum_{K \in \mathcal{T}_h} \left( - \int_K F_{ij}(U_h)(\mathbf{x}, t) U_{hi,j}(\mathbf{x}, t) dx + \int_{\partial K} \widehat{F}_{ij}(\mathbf{x}, t) \mathbf{n}_j U_{hi}(\mathbf{x}, t) dx \right). $$

Next, we show can define the numerical trace $\widehat{F}_{ij}$ in such a way that the above quantity is non-negative. Dropping the argument $t$, we get

$$ \Theta_h = \left( - \int_K F_{ij} U_{hi,j} dx + \int_{\partial K} \widehat{F}_{ij} \mathbf{n}_j U_{hi} dx \right) = \sum_{K \in \mathcal{T}_h} \int_{\partial K} \left( \widehat{F}_{ij} \mathbf{n}_j U_{hi} - c u \mathbf{q}_j \mathbf{n}_j \right) dx, $$

since $F_{ij} U_{hi,j} = (c u \mathbf{q}_j)_{,j}$. Then, setting $\llbracket V_i \rrbracket_j = V_i^- n_j^- + V_i^+ n_j^+$, we get

$$ \Theta_h = \sum_{e \in \mathcal{E}_h} \int_e \left( \widehat{F}_{ij} \llbracket U_{hi} \rrbracket_j - c \llbracket u \mathbf{q}_j \rrbracket_j \right) dx = \sum_{e \in \mathcal{E}_h} \int_e \left( \widehat{F}_{ij} - \{F_{ij}\} \right) \cdot \llbracket U_{hi} \rrbracket_j dx, $$

since

$$ c \llbracket u \mathbf{q}_j \rrbracket_j = c \{u\} \llbracket \mathbf{q}_j \rrbracket_j + \{\mathbf{q}_j\} \llbracket u \rrbracket_j = c \{u\} \delta_{ij} \llbracket \mathbf{q}_i \rrbracket_j + \{\mathbf{q}_j\} \llbracket u \rrbracket_j = \{F_{ij}\} \llbracket U_{hi} \rrbracket_j. $$

Thus, if we take,

$$ \widehat{F}_{ij} = \{F_{ij}\} + C_{ijk\ell} \llbracket U_{hk} \rrbracket_{\ell}, $$

we obtain

$$ \Theta_h = \sum_{e \in \mathcal{E}_h} \int_e C_{ijk\ell} \llbracket U_{hk} \rrbracket_{\ell} \llbracket U_{hi} \rrbracket_j ds \geq 0, $$

provided $C_{ijk\ell}$ is non-negative definite.

**Examples of DG methods.** We give three examples of DG methods for this case. We only have to determine their numerical trace $\widehat{F}_{ij}$. Two of them are widely known in computational fluids dynamics as the up-winding and Lax-Friedrichs numerical fluxes. The third has been recently discovered in the context of DG methods for second-order elliptic problems. To write them down, we use the following notation:

$$ \llbracket \mathbf{q} \rrbracket = \llbracket \mathbf{q}_j \rrbracket_j. $$

They are as follows:

$$ (i) \quad \widehat{F}_{ij} = \{F_{ij}\} + \frac{|c|}{2} \llbracket \mathbf{q} \rrbracket \delta_{ij} + \frac{|c|}{2} \llbracket u \rrbracket_j \delta_{i(d+1)} \quad \text{(up-winding)}, $$

$$ (ii) \quad \widehat{F}_{ij} = \{F_{ij}\} + \frac{|c|}{2} \llbracket \mathbf{q}_i \rrbracket_j + \frac{|c|}{2} \llbracket u \rrbracket_j \delta_{i(d+1)} \quad \text{(Lax-Friedrichs)}, $$

$$ (iii) \quad \widehat{F}_{ij} = \{F_{ij}\} + (C_{22} \llbracket \mathbf{q} \rrbracket - C_{12} \cdot \llbracket u \rrbracket) \delta_{ij} + ((C_{12})_j \llbracket \mathbf{q} \rrbracket + C_{11} \llbracket u \rrbracket)_j \delta_{i(d+1)}. $$

The last numerical trace is, of course, a generalization of the up-winding numerical trace.

It is not difficult to see that the dissipation term $\Theta_h$ for each of these numerical traces is given by

$$ (i) \quad \Theta_h = \sum_{e \in \mathcal{E}_h} \int_e \left( \frac{|c|}{2} \llbracket \mathbf{q} \rrbracket^2 + \frac{|c|}{2} \llbracket u \rrbracket^2 \right) ds \quad \text{(up-winding)}, $$

$$ (ii) \quad \Theta_h = \sum_{e \in \mathcal{E}_h} \int_e \left( \frac{|c|}{2} \llbracket \mathbf{q}_i \rrbracket \cdot \llbracket \mathbf{q}_i \rrbracket + \frac{|c|}{2} \llbracket u \rrbracket^2 \right) ds \quad \text{(Lax-Friedrichs)}, $$

$$ (iii) \quad \Theta_h = \sum_{e \in \mathcal{E}_h} \int_e (C_{22} \llbracket \mathbf{q} \rrbracket^2 + C_{11} \llbracket u \rrbracket^2) ds. $$

Note that the dissipation of the first two methods is proportional to the speed of propagation $|c|$. Note also that the DG method with Lax-Friedrichs numerical trace is more dissipative than the DG method with the up-winding numerical trace since

$$ \llbracket \mathbf{q}_i \rrbracket \cdot \llbracket \mathbf{q}_i \rrbracket \geq \llbracket \mathbf{q} \rrbracket^2, $$
where the equality holds only if the tangential component of $q$ is continuous. Finally, note that if $C_{22} = C_{11} = |c|/2$, the third method produces the same amount of dissipation than the DG method using the up-winding numerical trace. This indicates that, for that method, the vector-parameter $C_{12}$ does not have any stabilizing effect; it could be used, however, to enhance the accuracy of the method.

**Some properties.** It is not difficult to extend to the case of linear hyperbolic systems, what was discussed for the transport equation. Let us just briefly point out that the DG methods for those systems (i) are related to finite volume methods like the up-winding and the Lax-Friedrichs methods and always are locally conservative, (ii) are high-order accurate (in fact, of order $(k+1/2)$) when polynomials of degree $k$ are used) while remaining highly parallelizable when discretized in time with explicit methods, and (iii) have a dissipative mechanism that solely depends on the jumps which, in turn, are strongly linked with the residuals inside the elements.

A point we are particularly interested in stressing here is that working with the wave equation allows us to see, in a very natural way, how to discretize the Laplace operator $-\Delta$ by using DG methods. Indeed, since the hyperbolic system for the wave equation can be rewritten as

$$q_t - c\nabla u = 0, \quad u_t - c\nabla \cdot q = 0,$$

if we eliminate the term $u_t$ from the second equation and formally replace $q_t$ by $q$, we get

$$q - c\nabla u = 0, \quad -c\nabla \cdot q = 0,$$

which is nothing but a rewriting of the equation $c^2\Delta u = 0$. Thus, we can establish a one-to-one correspondence between the numerical traces used to discretize the hyperbolic system for the wave equation and those used to discretize the Laplace operator. This shows that switching from hyperbolic problems to elliptic ones is, in fact, does not entail a dramatic change as far as DG discretizations are concerned.

## 2.3 Second-order elliptic problems

In this sub-section, we consider DG methods for the model elliptic problem

$$-\Delta u = f \quad \text{in } \Omega, \quad u = 0 \quad \text{on } \partial\Omega,$$

where $\Omega$ is a bounded domain of $\mathbb{R}^d$. Following what was done for the wave equation, to define them, we rewrite our elliptic model problem as

$$q = \nabla u, \quad -\nabla \cdot q = f \quad \text{in } \Omega, \quad u = 0 \quad \text{on } \partial\Omega.$$

If when dealing with hyperbolic problems, we emphasized the DG methods are a generalization of finite volume methods, here we are going to show that DG methods are in fact mixed finite element methods. We also emphasize the fact that also in this context, DG methods are locally conservative methods ideally suited for adaptivity. Finally, we show that the dissipation mechanism of the DG methods, which is associated to the idea of artificial diffusion in the framework of hyperbolic equations, is associated to the idea of *penalization* of the discontinuity jumps in this context.

**The DG methods.** A DG numerical method is obtained as follows. After discretizing the domain $\Omega$, the approximate solution $(q_h, u_h)$ on the element $K$ is taken in the space $Q(K) \times U(K)$ and is determined by requiring that

$$\int_K q_h \cdot v \, dx = - \int_K u_h \nabla \cdot v \, dx + \int_{\partial K} \hat{u}_h v \cdot n \, ds,$$

$$\int_K q_h \cdot \nabla w \, dx - \int_{\partial K} w \hat{q}_h \cdot n \, ds = \int_K f w \, dx,$$

for all $(v, w) \in Q(K) \times U(K)$. Note that now we have two numerical traces, namely, $\hat{u}_h$ and $\hat{q}_h$, that remain to be defined.

To do that, we begin by finding a stability result for the solution of the original equation. To do that, we multiply the first equation by $q$ and integrate over $\Omega$ to get

$$\int_{\Omega} |q|^2 dx - \int_{\Omega} q \cdot \nabla u \, dx = 0.$$

Then, we multiply the second equation by $u$ and integrate over $\Omega$ to obtain

$$- \int_{\Omega} \nabla \cdot q u \, dx = \int_{\Omega} f u \, dx.$$
Adding these two equations, we get

$$\int_{\Omega} |\boldsymbol{q}|^2 dx = \int_{\Omega} f u \, dx.$$

This is the result we sought. Next, we mimic this procedure for the DG method.

We begin by taking $\boldsymbol{v} = \boldsymbol{q}_h$ in the first equation defining the DG method and adding on the elements $K$ to get

$$\int_{\Omega} |\boldsymbol{q}_h|^2 dx - \sum_{K \in \mathcal{T}_h} \left( - \int_K u_h \nabla \cdot \boldsymbol{q}_h \, d\boldsymbol{x} + \int_{\partial K} \hat{u}_h \boldsymbol{q}_h \cdot \boldsymbol{n} \, ds \right) = 0.$$

Next, we take $w = u_h$ in the second equation and add on the elements to obtain

$$\sum_{K \in \mathcal{T}_h} \left( \int_K \boldsymbol{q}_h \cdot \nabla u_h \, d\boldsymbol{x} - \int_{\partial K} u_h \hat{\boldsymbol{q}}_h \cdot \boldsymbol{n} \, ds \right) = \int_{\Omega} f u_h \, d\boldsymbol{x}.$$

Adding the two above equations, we find that

$$\int_{\Omega} |\boldsymbol{q}_h|^2 dx + \Theta_h = \int_{\Omega} f u \, dx,$$

where

$$\Theta_h = - \sum_{K \in \mathcal{T}_h} \left( - \int_K \nabla \cdot (u_h \boldsymbol{q}_h) \, d\boldsymbol{x} + \int_{\partial K} (\hat{u}_h \boldsymbol{q}_h \cdot \boldsymbol{n} + u_h \hat{\boldsymbol{q}}_h \cdot \boldsymbol{n}) \, ds \right).$$

It only remains to show that we can define consistent numerical traces $\hat{u}_h$ and $\hat{\boldsymbol{q}}_h$ that render $\Theta_h$ non-negative. Since,

$$\begin{align*}
\Theta_h &= \sum_{K \in \mathcal{T}_h} \int_{\partial K} (u_h \boldsymbol{q}_h \cdot \boldsymbol{n} - \hat{u}_h \boldsymbol{q}_h \cdot \boldsymbol{n} - u_h \hat{\boldsymbol{q}}_h \cdot \boldsymbol{n}) \, ds \\
&= \sum_{e \in \mathcal{E}_h} \int_e \left[ \left[ u_h \boldsymbol{q}_h - \hat{u}_h \boldsymbol{q}_h - u_h \hat{\boldsymbol{q}}_h \right] \right] ds \\
&= \sum_{e \in \mathcal{E}_{ih}} \int_e \left( \left[ \left[ u_h \boldsymbol{q}_h \right] - \hat{u}_h \left[ \boldsymbol{q}_h \right] - \left[ u_h \right] \cdot \hat{\boldsymbol{q}}_h \right) \right) ds + \int_{\partial \Omega} (u_h \boldsymbol{q}_h - \hat{u}_h \boldsymbol{q}_h \cdot \boldsymbol{n} - u_h \hat{\boldsymbol{q}}_h \cdot \boldsymbol{n}) \, ds \\
&= \sum_{e \in \mathcal{E}_{ih}} \int_e \left( \left( \left\{ u_h \right\} - \hat{u}_h \right) \left[ \boldsymbol{q}_h \right] + \left[ u_h \right] \cdot \left( \left\{ \boldsymbol{q} \right\} - \hat{\boldsymbol{q}}_h \right) \right) ds + \int_{\partial \Omega} (u_h (\boldsymbol{q}_h - \hat{\boldsymbol{q}}_h) \cdot \boldsymbol{n} - \hat{u}_h \boldsymbol{q}_h \cdot \boldsymbol{n}) \, ds,
\end{align*}$$

it is enough to take, inside the domain $\Omega$,

$$\hat{\boldsymbol{q}}_h = \{\boldsymbol{q}_h\} + C_{11} \llbracket u_h \rrbracket + C_{12} \llbracket \boldsymbol{q}_h \rrbracket, \quad \hat{u}_h = \{u_h\} - C_{12} \cdot \llbracket u_h \rrbracket + C_{22} \llbracket \boldsymbol{q}_h \rrbracket,$$

and on its boundary,

$$\hat{\boldsymbol{q}}_h = \boldsymbol{q}_h - C_{11} u_h \boldsymbol{n}, \quad \hat{u}_h = 0,$$

to finally get

$$\Theta_h = \sum_{e \in \mathcal{E}_{ih}} \int_e \left( C_{22} \llbracket \boldsymbol{q}_h \rrbracket^2 + C_{11} \llbracket u_h \rrbracket^2 \right) ds + \int_{\partial \Omega} C_{11} u_h^2 ds \geq 0,$$

provided $C_{11}$ and $C_{22}$ are non-negative. Note how the boundary conditions are imposed weakly through the definition of the numerical traces. This completes the definition of the DG methods.

**Some properties.** (i) Let us show that to guarantee the existence and uniqueness of the approximate solution of the DG methods, the parameter $C_{11}$ has to be greater than zero and the local spaces $\mathcal{U}(K)$ and $\mathcal{Q}(K)$ must satisfy the following compatibility condition:

$$u_h \in \mathcal{U}(K) : \quad \int_K \nabla u_h \boldsymbol{v} \, dx = 0 \quad \forall \boldsymbol{v} \in \mathcal{Q}(K) \quad \text{then} \quad \nabla u_h = \mathbf{0}.$$
Indeed, the approximate solution is well defined if and only if, the only approximate solution to the problem with $f = 0$ is the trivial solution. In that case, our stability identity gives

$$\int_{\Omega} |\mathbf{q}_h|^2 dx + \sum_{e \in \mathcal{E}_{ih}} \int_{e} (C_{22} [[\mathbf{q}_h]]^2 + C_{11} [[u_h]]^2) ds + \int_{\partial \Omega} C_{11} u_h^2 ds = 0,$$

which implies that $\mathbf{q}_h = 0$, $[[u_h]] = 0$ on $\mathcal{E}_{ih}$, and $u_h = 0$ on $\partial\Omega$, provided that $C_{11} > 0$. We can now rewrite the first equation defining the method as follows:

$$\int_K \nabla u_h \mathbf{v} \, dx = 0, \quad \forall \mathbf{v} \in \mathcal{Q}_h,$$

which, by the compatibility condition, implies that $\nabla u_h = 0$. Hence $u_h = 0$, as wanted.

(ii) When all the local spaces contain the polynomials of degree $k$, the orders of convergence of the $L^2$-norms of the errors in $\mathbf{q}$ and $u$ are $k$ and $k + 1$, respectively; see [22] when $C_{11}$ is of order $\mathcal{O}(h^{-1})$. Superconvergence in $\mathbf{q}$ has been proven and numerically observed in Cartesian grids, with a special choice of the numerical fluxes and equal-order elements with $Q^k$-polynomials; see [31].

(iii) Next, we show that DG methods are in fact mixed finite element methods. To see this, let us begin by noting that the DG approximate solution $(\mathbf{q}_h, u_h)$ can be also be characterized as the solution of

$$\begin{aligned} a(\mathbf{q}_h, \mathbf{v}) + b(u_h, \mathbf{v}) &= 0, \\ -b(w, \mathbf{q}_w) + c(u_h, w) &= F(w), \end{aligned}$$

for all $(\mathbf{v}, w) \in \mathcal{Q}_h \times \mathcal{U}_h$ where

$$\mathcal{Q}_h = \{\mathbf{v} : \mathbf{v} \in \mathcal{Q}(K) \ \forall K \in \mathcal{T}_h\}, \quad \mathcal{U}_h = \{w : w \in \mathcal{U}(K) \ \forall K \in \mathcal{T}_h\},$$

and

$$\begin{aligned} a(\mathbf{q}, \mathbf{r}) &:= \int_{\Omega} \mathbf{q} \cdot \mathbf{r} \, dx + \int_{\mathcal{E}_i} C_{22} [[\mathbf{q}]][[\mathbf{r}]] \, ds, \\ b(u, \mathbf{r}) &:= \sum_{K \in \mathcal{T}} \int_K u \nabla \cdot \mathbf{r} \, dx - \int_{\mathcal{E}_i} (\{u\} + C_{12} \cdot [[u]]) [[\mathbf{r}]] \, ds, \\ c(u, \mathbf{v}) &:= \int_{\mathcal{E}_{ih}} C_{11} [[u]] \cdot [[v]] \, ds + \int_{\partial \Omega} C_{11} uv \, ds, \\ F(\mathbf{r}) &:= \int_{\Omega} fv \, dx. \end{aligned}$$

As a consequence, the corresponding matrix equation has the form

$$\begin{pmatrix} \mathbf{A} & -\mathbb{B}^t \\ \mathbb{B} & \mathbf{C} \end{pmatrix} \begin{pmatrix} \mathcal{Q} \\ \mathcal{U} \end{pmatrix} = \begin{pmatrix} 0 \\ \mathbb{F} \end{pmatrix},$$

which is typical of stabilized mixed finite element methods. As it is well known, those methods are not well defined unless the 'stabilizing' form $c(\cdot, \cdot)$, usually associated with residuals, is introduced. For DG methods, the 'stabilizing' form $c(\cdot, \cdot)$ solely depends on the parameter $C_{11}$ and the jumps across elements of the functions in $\mathcal{U}_h$. This is why we could think that this form stabilizes the method by penalizing the jumps, $C_{11}$ being the penalization parameter; thus, what was interpreted to be the artificial dissipation coefficient in hyperbolic problems can now be thought of a penalization parameter.

Finally, let us emphasize that, for DG methods, penalizing the jumps is also a way of introducing stabilization by using residuals. Indeed, just as for the hyperbolic case, the residuals are related to the jumps. To see this, set $\mathbf{R}_1 = \mathbf{q}_h - \nabla u_h$ and $\mathbf{R}_2 = -\nabla \cdot \mathbf{q}_h - f$ and use the weak formulation of the DG method and the definition of its numerical trace to get

$$\begin{aligned} \int_K \mathbf{R}_1 \cdot \mathbf{v} \, dx &= \int_{\partial K} \left( \left( -\frac{1}{2} \mathbf{n} - \mathbf{C}_{12} \right) \cdot [[u_h]] + C_{22} [[\mathbf{q}_h]] \right) \mathbf{v} \cdot \mathbf{n} \, ds, \\ \int_K R_2 w \, dx &= \int_{\partial K} \left( \left( -\frac{1}{2} \mathbf{n} + \mathbf{C}_{12} \right) \cdot [[\mathbf{q}_h]] + C_{11} [[u_h]] \right) \mathbf{n} \, w \, ds, \end{aligned}$$

for all $(\mathbf{v}, w) \in \mathcal{Q}(K) \times \mathcal{U}(K)$.
(iv) Let us now comment on an essential difference between the DG and the classical mixed methods. To solve the system associated to classical mixed methods, namely,

$$ \begin{pmatrix} A & -\mathbb{B}^t \\ \mathbb{B} & C \end{pmatrix} \begin{pmatrix} \mathbb{Q} \\ \mathbb{U} \end{pmatrix} = \begin{pmatrix} 0 \\ \mathbb{F} \end{pmatrix}, $$

we can try to eliminate $Q$ from the equations to obtain

$$ (\mathbb{B}A^{-1}\mathbb{B}^t + \mathbb{C}) U = F. $$

Since the matrix $A$ is not easily invertible, due to continuity constraints on the approximation $q_h$, we can hybridize the method; see [19]. We then obtain a new system of the form

$$ \begin{pmatrix} A & -B^t & -D^t \\ B & C & 0 \\ D & 0 & 0 \end{pmatrix} \begin{pmatrix} \mathbb{Q} \\ \mathbb{U} \\ \mathbb{\Lambda} \end{pmatrix} = \begin{pmatrix} 0 \\ \mathbb{F} \\ 0 \end{pmatrix}, $$

where $\mathbb{\Lambda}$ is the vector of degrees of freedom associated to the so-called Lagrange multiplier $\lambda_h$. As is now well-known, the new vectors of degrees of freedom $Q$ and $U$ actually define the same approximation $(q_h, u_h)$ as the original mixed method. Moreover, since now $A$ is block-diagonal, both $Q$ and $U$ can now be easily eliminated to obtain an equation for the multiplier only, namely,

$$ \mathbb{E}\Lambda = \mathbb{H}, $$

where $\mathbb{E}$ and $\mathbb{H}$ are given by

$$ \begin{aligned} \mathbb{E} &= DA^{-1} (A - B^t (BA^{-1}B^t + C)^{-1} B) A^{-1}D^t, \\ \mathbb{H} &= -DA^{-1}B^t (BA^{-1}B^t + C)^{-1} \mathbb{F}. \end{aligned} \tag{5} $$

In the case of DG methods, the above procedure it not necessary as the matrix $A$ is block diagonal when $C_{22}$ is identically zero. In this case, the DG method remains well-posed and the matrix $A$ can be made equal to the identity, if a suitable basis is used. This allows us to easily eliminate $q_h$ from the equations. Finally, note that, unlike classical and stabilized mixed methods, this can be achieved even if polynomials of different degrees are used in different elements. This renders DG methods ideal for adaptivity.

(v) The methods we have presented are locally conservative. As we saw in the hyperbolic case, this is a reflection of the form of the weak formulation and the fact that the definition of the numerical traces on the face $e$ does not depend on what side of it we are. More general DG methods define the approximate solution by requiring that

$$ \int_K q_h \cdot r \, dx = - \int_K u_h \nabla \cdot r \, dx + \int_{\partial K} \hat{u}_{h,K} r \cdot n_K \, ds, $$

$$ \int_K q_h \cdot \nabla v \, dx = \int_K fv \, dx + \int_{\partial K} v \hat{q}_{h,K} \cdot n_K \, ds, $$

for all $(r, v) \in \mathcal{Q}(K) \times \mathcal{U}(K)$. In this general formulation, the numerical traces $\hat{u}_{h,K}$ and $\hat{q}_{h,K}$ can have definitions that might depend on what side of the element boundaries we are. Hence they are not locally conservative. This is the case for the numerical fluxes in $u$ of the last four schemes in Table 2 taken from [3]. The acronym LDG stands for local DG, and IP for interior penalty.

Finally, let us point out that, in that table, the function $\alpha^r([u_h])$ is a special stabilization term introduced by Bassi and Rebay [12] and later studied by Brezzi et al. [20]; its stabilization properties are equivalent to the one originally presented. In [4], a complete study of these methods is carried out in a single, unified approach.

## 2.4 The Oseen system

In this sub-section, we show how to discretize the Oseen system, namely,

$$ -\Delta u + (a \cdot \nabla) u + \nabla p = f, \quad \nabla \cdot u = 0 \quad \text{in } \Omega, \quad u = 0 \quad \text{on } \partial\Omega, $$

where $\Omega$ is a bounded domain of $\mathbb{R}^d$ and $\nabla \cdot a = 0$.

The objective here is to put to work what we have learned about DG-space discretizations for convection and diffusion operators.
<table><caption>**Table 2** Some DG methods and their numerical fluxes.</caption><thead><tr><th>Method</th><th>q&#770;<sub>e,K</sub></th><th>u&#770;<sub>h,K</sub></th></tr></thead><tbody><tr><td>Bassi-Rebay [10]</td><td>{q<sub>h</sub>}</td><td>{u<sub>h</sub>}</td></tr><tr><td>LDG [38]</td><td>{q<sub>h</sub>} + C<sub>11</sub> [[u<sub>h</sub>]] - C<sub>12</sub>[[q<sub>h</sub>]]</td><td>{u<sub>h</sub>} + C<sub>12</sub> · [[u<sub>h</sub>]]</td></tr><tr><td>DG [22]</td><td>{q<sub>h</sub>} + C<sub>11</sub> [[u<sub>h</sub>]] - C<sub>12</sub>[[q<sub>h</sub>]]</td><td>{u<sub>h</sub>} + C<sub>12</sub> · [[u<sub>h</sub>]] + C<sub>22</sub>[[q<sub>h</sub>]]</td></tr><tr><td>Brezzi et al. [20]</td><td>{q<sub>h</sub>} - α<sup>r</sup>([[u<sub>h</sub>]])</td><td>{u<sub>h</sub>}</td></tr><tr><td>IP [45]</td><td>{∇u<sub>h</sub>} + C<sub>11</sub> [[u<sub>h</sub>]]</td><td>{u<sub>h</sub>}</td></tr><tr><td>Bassi-Rebay [12]</td><td>{∇u<sub>h</sub>} - α<sup>r</sup>([[u<sub>h</sub>]])</td><td>{u<sub>h</sub>}</td></tr><tr><td>Baumann-Oden [15]</td><td>{∇u<sub>h</sub>}</td><td>{u<sub>h</sub>} - n<sub>K</sub> · [[u<sub>h</sub>]]</td></tr><tr><td>NIPG [70]</td><td>{∇u<sub>h</sub>} + C<sub>11</sub>[[u<sub>h</sub>]]</td><td>{u<sub>h</sub>} - n<sub>K</sub> · [[u<sub>h</sub>]]</td></tr><tr><td>Babuška-Zlámal [7]</td><td>C<sub>11</sub>[[u<sub>h</sub>]]</td><td>u<sub>h</sub>|<sub>K</sub></td></tr><tr><td>Brezzi et al. [20]</td><td>-α<sup>r</sup>([[u<sub>h</sub>]])</td><td>u<sub>h</sub>|<sub>K</sub></td></tr></tbody></table>

**The DG methods.** To define a DG method for the Oseen system, we begin by rewriting it as a first-order system,

$$ \sigma_i = \nabla u_i, \quad -\nabla \cdot \sigma_i + \nabla \cdot (a u_i) + \partial_i p = f_i, \quad 1 \le i \le d, \quad \nabla \cdot u = 0 \quad \text{in } \Omega, \quad u = 0 \quad \text{on } \partial\Omega. $$

where $u_i$ denotes the $i$-th component of the velocity $u$. To discretize the above equations, we take the approximate solution $(\sigma_h, u_h, p_h)$ on the element $K$ to be in the space $S(K)^d \times U(K)^d \times P(K)$ and determine it by imposing that, for $1 \le i \le d$, and for all $(\tau, v, w) \in S(K) \times U(K) \times P(K)$,

$$ \int_K \sigma_{i,h} \cdot \tau \, dx = - \int_K u_{i,h} \nabla \cdot \tau \, dx + \int_{\partial K} \hat{u}_{\sigma_{i,h}} \tau \cdot n \, ds, $$

$$ \int_K (\sigma_{i,h} \cdot \nabla v - u_{i,h} a \cdot \nabla v - p_h \partial_i v) \, dx + \int_{\partial K} (-\hat{\sigma}_{i,h} \cdot n v + \widehat{a u_{i,h}} \cdot n v + \hat{p}_h v n_i) \, ds = \int_K f_i v \, dx, $$

$$ - \int_K u_h \cdot \nabla q \, dx + \int_{\partial K} \hat{u}_{p,h} \cdot n \, q \, ds = 0. $$

It remains to find numerical traces that ensure that the method is well defined and stable. As usual, we begin by finding the stability equality for the continuous case. To do that, we first multiply the equation defining $\sigma_i$ by $\sigma_i$ sum over $i$ and integrate over $\Omega$ to get

$$ \int_{\Omega} |\sigma|^2 dx = \int_{\Omega} \sigma_i \cdot \nabla u_i dx, $$

where $|\sigma|^2 = \sigma_i \cdot \sigma_i$; we are using the Einstein summation convention. Now, we multiply the second equation of the Oseen system by $u_i$, sum over $i$ and integrate over $\Omega$ to get

$$ \int_{\Omega} \left( -\nabla \cdot \sigma_i u_i - \frac{1}{2} \nabla \cdot a |u|^2 + u \cdot \nabla p \right) dx = \int_{\Omega} f \cdot u \, dx, $$

where $|u|^2 = u_i u_i$. Adding the last two equations, and integrating by parts, we get

$$ \int_{\Omega} |\sigma|^2 dx - \int_{\Omega} p \nabla \cdot u \, dx = \int_{\Omega} f \cdot u \, dx, $$

and using the incompressibility of the velocity $u$,

$$ \int_{\Omega} |\sigma|^2 dx = \int_{\Omega} f \cdot u \, dx. $$

Next, we mimic the above procedure for the DG method. First, we take $\tau = \sigma_{i,h}$, add over $i$ and then over $K$, to get

$$ \int_{\Omega} |\sigma_h|^2 dx = - \sum_{K \in \mathcal{T}_h} \int_K u_{i,h} \nabla \cdot \sigma_{i,h} dx + \sum_{K \in \mathcal{T}_h} \int_{\partial K} \hat{u}_{\sigma_{i,h}} \sigma_{i,h} \cdot n \, ds. $$
Now, take $v = u_{i,h}$, add over $i$ and then over $K$ to obtain

$$ \sum_{K \in \mathcal{T}_h} \int_K \left( \boldsymbol{\sigma}_{i h} \cdot \nabla u_{i h} - \frac{1}{2} \boldsymbol{a} \cdot \nabla |\boldsymbol{u}_h|^2 - p_h \nabla \cdot \boldsymbol{u} \right) d x \\ + \sum_{K \in \mathcal{T}_h} \int_{\partial K} \left( -\hat{\boldsymbol{\sigma}}_{i h} \cdot \boldsymbol{n} u_{i h} + \widehat{\boldsymbol{a} \cdot \boldsymbol{u}_{i h}} \cdot \boldsymbol{n} u_{i h} + \hat{p}_h u_h \cdot \boldsymbol{n} \right) d s = \int_{\Omega} \boldsymbol{f} \cdot \boldsymbol{u}_h d x. $$

Adding the last two equations, we get

$$ \int_{\Omega} |\boldsymbol{\sigma}_h|^2 d x + \sum_{K \in \mathcal{T}_h} \int_K \left( \nabla \cdot (\boldsymbol{\sigma}_{i h} u_{i h}) - \frac{1}{2} \boldsymbol{a} \cdot \nabla |\boldsymbol{u}_h|^2 - p_h \nabla \cdot \boldsymbol{u} \right) d x \\ + \sum_{K \in \mathcal{T}_h} \int_{\partial K} \left( -\hat{u}_{\sigma, i h} \boldsymbol{\sigma}_{i h} \cdot \boldsymbol{n} - \hat{\boldsymbol{\sigma}}_{i h} \cdot \boldsymbol{n} u_{i h} + \widehat{\boldsymbol{a} \cdot \boldsymbol{u}_{i h}} \cdot \boldsymbol{n} u_{i h} + \hat{p}_h u_h \cdot \boldsymbol{n} \right) d s = \int_{\Omega} \boldsymbol{f} \cdot \boldsymbol{u}_h d x. $$

Finally, we take $q = p_h$, add over $K$ and add the result to the above equation to get

$$ \int_{\Omega} |\boldsymbol{\sigma}|^2 d x + \Theta_h = \int_{\Omega} \boldsymbol{f} \cdot \boldsymbol{u} d x, $$

where

$$ \Theta_h = \Theta_h^1 + \Theta_h^2, $$

and

$$ \Theta_h^1 = \sum_{K \in \mathcal{T}_h} \int_K \nabla \cdot \left( \boldsymbol{\sigma}_{i h} u_{i h} - \frac{1}{2} \boldsymbol{a} |\boldsymbol{u}_h|^2 - p_h \boldsymbol{u}_h \right) d x \\ = \sum_{e \in \mathcal{E}_h} \int_e \left[ \boldsymbol{\sigma}_{i h} u_{i h} - \frac{1}{2} \boldsymbol{a} |\boldsymbol{u}_h|^2 - p_h \boldsymbol{u}_h \right] d x, $$

and

$$ \Theta_h^2 = \sum_{K \in \mathcal{T}_h} \int_{\partial K} \left( -\hat{u}_{\sigma, i h} \boldsymbol{\sigma}_{i h} \cdot \boldsymbol{n} - \hat{\boldsymbol{\sigma}}_{i h} \cdot \boldsymbol{n} u_{i h} + \widehat{\boldsymbol{a} \cdot \boldsymbol{u}_{i h}} \cdot \boldsymbol{n} u_{i h} + \hat{p}_h u_h \cdot \boldsymbol{n} + \hat{\boldsymbol{u}}_{p,h} \cdot \boldsymbol{n} p_h \right) d s \\ = \sum_{e \in \mathcal{E}_h} \int_e \left[ -\hat{u}_{\sigma, i h} \boldsymbol{\sigma}_{i h} - \hat{\boldsymbol{\sigma}}_{i h} u_{i h} + \widehat{\boldsymbol{a} \cdot \boldsymbol{u}_{i h}} u_{i h} + \hat{p}_h u_h + \hat{\boldsymbol{u}}_{p,h} p_h \right] d s. $$

Adding and rearranging terms, we get

$$ \Theta_h = \sum_{e \in \mathcal{E}_{ih}} \int_e \left( \left( \{u_{ih}\} - \hat{u}_{\sigma,ih} \right) \llbracket \boldsymbol{\sigma}_{ih} \rrbracket + \left( \{\boldsymbol{\sigma}_{ih}\} - \hat{\boldsymbol{\sigma}}_{ih} \right) \llbracket u_{ih} \rrbracket - \left( \{\boldsymbol{a} u_{ih}\} - \widehat{\boldsymbol{a} u_{ih}} \right) \llbracket u_{ih} \rrbracket \right. \\ \left. - \left( \{p_h\} - \hat{p}_h \right) \llbracket \boldsymbol{u}_h \rrbracket - \left( \{\boldsymbol{u}_h\} - \hat{\boldsymbol{u}}_{p,h} \right) \llbracket p_h \rrbracket \right) d s \\ + \int_{\partial \Omega} \left( -\hat{u}_{\sigma,ih} \boldsymbol{\sigma}_{ih} + (\boldsymbol{\sigma}_{ih} - \hat{\boldsymbol{\sigma}}_{ih}) u_{ih} - \left( \frac{1}{2} \boldsymbol{a} u_{ih} - \widehat{\boldsymbol{a} u_{ih}} \right) u_{ih} - (p_h - \hat{p}_h) \boldsymbol{u}_h + \hat{\boldsymbol{u}}_{p,h} p_h \right) \cdot \boldsymbol{n} d s. $$

We now take the numerical traces as to render $\Theta_h$ non-negative. On the interior of the domain, we take

$$ \hat{\boldsymbol{\sigma}}_{i h} = \{\boldsymbol{\sigma}_{i h}\} - C_{11} \llbracket u_{i h} \rrbracket + C_{12} \llbracket \boldsymbol{\sigma}_{i h} \rrbracket, \\ \widehat{\boldsymbol{a} \cdot \boldsymbol{u}_{i h}}(\boldsymbol{x}) = \boldsymbol{a} \lim_{\epsilon \downarrow 0} u_{i,h}(\boldsymbol{x} - \epsilon \boldsymbol{a}), \\ \hat{\boldsymbol{u}}_{\sigma, i h} = \{\boldsymbol{u}_{i h}\} - C_{12} \llbracket \boldsymbol{u}_{i h} \rrbracket, $$

and, on the boundary,

$$ \hat{\boldsymbol{\sigma}}_{i h} = \boldsymbol{\sigma}_{i h} - C_{11} u_{i h} \boldsymbol{n}, \quad \hat{\boldsymbol{u}}_{\sigma, h} = \mathbf{0}. $$
The numerical traces associated with the incompressibility constraint, $\hat{u}_{p,h}$ and $\hat{p}_h$, are defined by using an analogous recipe. In the interior of $\Omega$, we take

$$ \hat{u}_{p,h} = \{u_h\} + D_{11} \llbracket p_h \rrbracket - D_{12} \llbracket u_h \rrbracket, $$

$$ \hat{p}_h = \{p_h\} + D_{12} \cdot \llbracket p_h \rrbracket, $$

and on the boundary, we take

$$ \hat{u}_{p,h} = 0, \quad \hat{p}_h = p_h^+. $$

With the above choice, we get

$$ \Theta_h = \sum_{e \in \mathcal{E}_{ih}} \int_e ((C_{11} + |\mathbf{a} \cdot \mathbf{n}|) |\llbracket \mathbf{u}_h \rrbracket|^2 + D_{11} \llbracket p_h \rrbracket^2) ds + \frac{1}{2} \int_{\partial \Omega} |\mathbf{a} \cdot \mathbf{n}| |\mathbf{u}_h|^2 ds \geq 0, $$

whenever the parameters $C_{11}$ and $D_{11}$ are non-negative. This completes the definition of the DG method for the Oseen system.

**Some properties.** (i) Let us show that the above DG method is well defined when the stabilization parameters $C_{11}$ and $D_{11}$ are positive and when the two following compatibility conditions on the local spaces are satisfied:

$$ u_h \in \mathcal{U}(K) : \quad \int_K \tau \cdot \nabla u_h \, dx = 0 \quad \forall \tau \in S^d(K) \quad \text{then} \quad \nabla u_h = 0, $$

and

$$ p_h \in \mathcal{P}(K) : \quad \int_K v \cdot \nabla p_h \, dx = 0 \quad \forall v \in \mathcal{U}^d(K) \quad \text{then} \quad \nabla p_h = 0. $$

Thus, we must show that when $f = 0$, the only solution is $u_h = \sigma_{ih} = 0$ and $p_h$ is a constant. From the stability result, we get that

$$ \int_{\Omega} |\boldsymbol{\sigma}|^2 dx + \sum_{e \in \mathcal{E}_{ih}} \int_e ((C_{11} + |\boldsymbol{a} \cdot \boldsymbol{n}|) |\llbracket \boldsymbol{u}_h \rrbracket|^2 + D_{11} \llbracket p_h \rrbracket^2) ds + \frac{1}{2} \int_{\partial \Omega} |\boldsymbol{a} \cdot \boldsymbol{n}| |\boldsymbol{u}_h|^2 ds = 0, $$

which implies that $\sigma_{ih} = 0$, $\llbracket u_{ih} \rrbracket = 0$, $u_{ih} = 0$ on $\partial\Omega$, and $\llbracket p_h \rrbracket = 0$.

With this information, the first equation defining the DG method becomes

$$ \int_K \tau \cdot \nabla u_{ih} \, dx = 0 \quad \forall \tau \in S^d, $$

which, by the first compatibility condition implies that $\nabla u_{ih} = 0$ on each element $K$. Since $u_h$ is continuous and zero in the border, we must have that $u_h = u$.

Finally, the second equation defining the DG method becomes

$$ \int_K v \cdot \nabla p_h \, dx = 0 \quad \forall v \in \mathcal{U}^d(K), $$

which, by the second compatibility condition, implies that $\nabla p_h = 0$ on each element $K$. Since $p_h$ is continuous, this shows that $p_h$ is a constant.

(ii) In obtaining error estimates, one of the main issues is the so-called inf-sup condition. It turns out that, when $D_{11} > 0$, it is possible to circumvent that condition and obtain optimal error estimates. Indeed, if we take $S(K) = U(K) = P(K)$ to be the space of polynomials of degree $k$, the L$^2$-norm of the error in $p$ and $\sigma_i$ converge with order $k$, and the L$^2$-norm of the velocity with order $k+1$. Moreover, if polynomials of degree $k-1$ are used to approximate the pressure $p$ and the stress tensor $\sigma_i$, the above mentioned orders of convergence remain the same. This method is, however, not more efficient; see [33].

The fact that the inf-sup condition can be circumvented is typical of stabilized numerical methods. In the case of DG methods, this happens because it is possible to obtain a much weaker, generalized inf-sup condition. Moreover, in some cases, this inf-sup condition can be proven even if we take $D_{11} = 0$. In such a case, the conditions on the velocity and pressure spaces are more restrictive, of course; see [71, 76].

(iii) Let us end this section by pointing out that the method is locally conservative and that all the local residuals can be easily computed solely in terms of the jumps of the approximate solution. Moreover, the method is robust with respect to the Reynolds number. To support this claim, we quote a result from [32] about the Kovasznay flow with different Reynolds numbers. We use quadrilateral meshes generated by consecutive refinements. In each refinement step, each grid cell is divided into four similar cells by connecting the edge midpoints. Therefore, grid level $L$ corresponds to a mesh-size $h_L = 2^{1-L}$. All the unknowns are approximated with tensor product polynomials of degree $k \ge 1$. The discontinuity and pressure stabilization functions $\sigma$ and $\delta$ are chosen of the order $O(1/h)$ and $O(h)$, respectively. Fig. 2 shows robustness of the discretization with respect to the Reynolds number.
<figure><img src="image_3.png" /><figcaption>Fig. 2 Scaled L²-errors in $u$ with $k = 2$ for different Reynolds numbers.</figcaption></figure>

## 3 Non-linear hyperbolic conservation laws

In all the previous section, we have considered DG-space discretizations for linear problems ranging from ODEs to the Oseen system. In computational fluid dynamics, the next natural step would be to consider the incompressible case. However, this topic is currently being vigorously addressed by several researchers; see the pioneering work [9, 61].

In this section, we switch to the numerical solution of non-linear hyperbolic equations. One of the main applications is to devise a locally conservative, stable and high-order accurate method to discretize the Euler equations of gas dynamics, as this is the bottleneck to deal with the compressible Navier-Stokes equations.

### 3.1 The RKDG methods: Introduction

The DG methods we consider are called the Runge Kutta discontinuous Galerkin (RKDG) methods, [30, 35–37, 39]. To describe these methods, we use the simple model problem of the non-linear hyperbolic scalar conservation law

$$u_t + \nabla f(u) = 0.$$

The RKDG methods are obtained in three steps:

**Step 1: The DG space discretization.** First, the conservation law is discretized in space by using a DG method. A discontinuous approximate solution $u_h$ is sought such that when restricted to the element $K$, it belongs to the finite dimensional space $\mathcal{U}(K)$. It is defined by imposing that, for all $v_h \in \mathcal{U}(K)$,

$$\int_K (u_h)_t v_h \, dx - \int_K f(u_h) \cdot \nabla v_h \, dx + \int_{\partial K} \hat{f}(u_h) \cdot n_K v_h \, ds = 0.$$

Here, as we have seen in the previous section, the proper definition of the numerical trace $\hat{f}(u_h)$ (called, in this kind of problems, the approximate Riemann solver or also the numerical flux) is essential for the stability and convergence of the method.

**Step 2: The RK time discretization.** Then, we discretize the resulting system of ordinary differential equations, $\frac{d}{dt}u_h = L(u_h)$, by using special explicit Runge-Kutta (RK) methods:

1. Set $u_h^{(0)} = u_h^n$.

2. For $i = 1, \dots, K$ compute the intermediate functions:

$$u_h^{(i)} = \sum_{l=0}^{i-1} \alpha_{il} w_h^{il}, \quad w_h^{il} = u_h^{(l)} + \frac{\beta_{il}}{\alpha_{il}} \Delta t^n L_h(u_h^{(l)}).$$

3. Set $u_h^{n+1} = u_h^K$.
The distinctive feature of these RK methods is that their stability follows from the stability of the mapping $u_h^{(l)} \mapsto w_h^{il}$ defining the intermediate steps. Unfortunately, this mapping is not necessarily stable and, as a consequence, the method requires another component to enforce stability.

**Step 3: The generalized slope limiter.** This component is nothing but the so-called generalized slope limiter $\Lambda\Pi_h$. This non-linear projection operator, is devised in such a way that if $u_h^{(l)} = \Lambda\Pi_h v_h$ for some function $v_h$, then the mapping $u_h^{(l)} \mapsto w_h^{il}$ is stable.

Thus, we incorporate the generalized slope limiter in the above time-marching algorithm as follows:

1. Set $u_h^{(0)} = u_h^n$.

2. For $i = 1, \dots, K$ compute the intermediate functions:

$$u_h^{(i)} = \Lambda\Pi_h \left( \sum_{l=0}^{i-1} \alpha_{il} w_h^{il} \right), \quad w_h^{il} = u_h^{(l)} + \frac{\beta_{il}}{\alpha_{il}} \Delta t^n L_h (u_h^{(l)})$$

3. Set $u_h^{n+1} = u_h^K$.

This is the general form of the RKDG methods

In what follows, we elaborate on each of the above points and pay special attention to how to enforce the stability of the method. This is by far more delicate than what was done in the previous section because, as we can see, the stability of the method relies on the stability of the forward Euler step $u_h^{(l)} \mapsto w_h^{il}$. Unfortunately, this operator is always unstable in the $L^2$-norm. Nevertheless, we are going to show that it is possible to find a weaker stability property for the forward Euler and then enforce it with the generalized slope limiter to obtain a stable method. The remarkable fact is that this can be achieved while maintaining the high-order accuracy of the method.

## 3.2 The DG-space discretization

With respect to the DG-space discretization, we do not have to say too many things that are different from the DG-space discretization of the transport equation. However, let us very briefly emphasize, that the method is locally conservative, and has order $k+1/2$ when polynomials of degree $k$ are used. It is highly parallelizable since its mass matrix is block diagonal and since, as it is typical of DG methods, the numerical trace $\hat{f}$ solely depends on the traces of $u_h$ on both sides of the inter-element boundary. The classical examples are the following:

(i) The Godunov flux:

$$\hat{f}^G(a, b) = \begin{cases} \min_{a \le u \le b} f(u), & \text{if } a \le b, \\ \max_{b \le u \le a} f(u), & \text{otherwise.} \end{cases}$$

(ii) The Engquist-Osher flux:

$$\hat{f}^{EO}(a, b) = \int_{0}^{b} \min(f'(s), 0) ds + \int_{0}^{a} \max(f'(s), 0) ds + f(0).$$

(iii) The Lax-Friedrichs flux:

$$\hat{f}^{LF}(a, b) = \frac{1}{2} [f(a) + f(b) - C(b - a)], \quad C = \max_{\inf u^0(x) \le s \le \sup u^0(x)} |f'(s)|.$$

When piecewise constant approximate solutions are taken and when the forward Euler method is used to march in time, a monotone scheme is obtained. Monotone schemes are not only very stable but converge to the physically relevant solution, the so-called entropy solution. Unfortunately, monotone schemes are at most first order accurate. By using high-degree polynomials in a DG-space discretization, the accuracy of the scheme is raised. However, a RK time-discretization that has the same accuracy in time and renders the scheme stable has to be devised.

## 3.3 The RK time discretization

The RK method we consider is required to satisfy the following conditions:

(i) If $\beta_{il} \ne 0$ then $\alpha_{il} \ne 0$.

(ii) $\alpha_{il} \ge 0$.

(iii) $\sum_{l=0}^{i-1} \alpha_{il} = 1$.
Note that, by the first property, we can express the RK method in terms of the functions $w_h^{il}$. Note also that if we assume that, for some semi-norm $|\cdot|$, we have that $|w_h^{il}| \le |u_h^{(l)}|$, then

$$
\begin{align*}
|u_h^{(i)}| &= \left| \sum_{l=0}^{i-1} \alpha_{il} w_h^{il} \right|, \\
&\le \sum_{l=0}^{i-1} \alpha_{il} |w_h^{il}|, && \text{by the positivity property (ii)}, \\
&\le \sum_{l=0}^{i-1} \alpha_{il} |u_h^{(l)}|, && \text{by the stability assumption,} \\
&\le \max_{0 \le l \le i-1} |u_h^{(l)}|, && \text{by the consistency property (iii)},
\end{align*}
$$

which readily implies that $|u_h^n| \le |u_h^0|$, $\forall n \ge 0$, by a simple induction argument. In other words, the stability of the Euler
forward step $u_h^{(l)} \mapsto w_h^{il}$ implies the stability of the RK method!

We must also keep in mind that, the RK method being explicit, we must ensure that the round-off errors are not amplified.
For DG-space discretizations using polynomials of degree $k$ and a $(k+1)$-stage RK method of order $k+1$ (which give rise
to a $(k+1)$-st order accurate method), a von Neumann stability analysis for the one-dimensional linear case $f(u) = cu$ gives
us the stability

$$
|c| \frac{\Delta t}{\Delta x} \le \frac{1}{2k+1}.
$$

This is a condition to be respected, since otherwise the round-off errors will be amplified.

3.4 The stability of the step $u_h \mapsto w_h = u_h + \delta L_h(u_h)$

It is not difficult to prove that the forward Euler step is *not* stable in the $L^2$-norm, except in the case in which polynomials of degree 0 are used. If polynomials of degree $k > 0$ are used, the Euler step can be rendered stable if $\Delta t/\Delta x$ is proportional to $(\Delta x)^p(k)$, where $p(k) > 0$; for example, $p(1) = 1/2$, see [23]. This is clearly an unacceptable situation for hyperbolic problems.

This negative result prompted the search of weaker norms, or semi-norms, for which the forward Euler step would be stable.
To do that, recall that, when using polynomials of degree zero, the DG methods are nothing but monotone schemes, which are
stable methods for the $L^\infty$ norm in several space dimensions and in the total-variation semi-norm (in one space dimension).
Thus, the idea is to see if these stability properties remain invariant when high-order degree polynomials are used.

To address this issue as clearly as possible, let us restrict ourselves to the one dimensional case. In this case, the forward
Euler step for the method reads

$$
\int_{I_j} \frac{(w_h - u_h)}{\delta t} v_h \, dx - \int_{I_j} f(u_h) (v_h)_x \, dx + \hat{f}(u_h) v_h \Big|_{x_{j-1/2}}^{x_{j+1/2}} = 0.
$$

Taking $v_h \equiv 1$, we obtain,

$$
(\bar{w}_j - \bar{u}_j) / \delta + \left( \hat{f} \left( u_{j+1/2}^{-}, u_{j+1/2}^{+} \right) - \hat{f} \left( u_{j-1/2}^{-}, u_{j-1/2}^{+} \right) \right) / \Delta_j = 0, \quad (6)
$$

where $\bar{u}_j$ denotes the mean of $u_h$ on the interval $I_j$, $u_{j+1/2}^{-}$ denotes the limit from the left and $u_{j+1/2}^{+}$ the limit from the right.
When the approximate solution is piecewise-constant, we obtain a monotone scheme for small enough values of $|\delta|$ and, as a
consequence, we do have that the scheme is total variation diminishing (TVD), that is, that

$$
| \bar{w}_h |_{TV(0,1)} \leq | \bar{u}_h |_{TV(0,1)},
$$

where

$$
| \bar{u}_h |_{TV(0,1)} \equiv \sum_{1 \le j \le N} | \bar{u}_{j+1} - \bar{u}_j |
$$

is the total variation of the local means.
<figure><img src="image_4.png" /><figcaption>Fig. 3 Example of slope limiters: The MUSCL limiter (left) and the less restrictive $\Lambda\Pi_h^1$ limiter (right). Displayed are the local means of $u_h$ (thick line), the linear function $u_h$ in the element of the middle before limiting (dotted line), and the resulting function after limiting (solid line).</figcaption></figure>

For approximate solutions that are not necessarily piecewise-constant, the above result still holds provided that the following sign conditions are satisfied:

$$ \operatorname{sign}(u_{j+1/2}^{+} - u_{j-1/2}^{+}) = \operatorname{sign}(\bar{u}_{j+1} - \bar{u}_{j}), $$

$$ \operatorname{sign}(u_{j+1/2}^{-} - u_{j-1/2}^{-}) = \operatorname{sign}(\bar{u}_{j} - \bar{u}_{j-1}). $$

Since these conditions are not necessarily satisfied it is necessary to enforce them by means of what will be called a generalized slope limiter, $\Lambda\Pi_h$.

We see that, unlike what happened in the previous section, the introduction of the numerical traces is not enough to guarantee the stability of the DG method. For non-linear hyperbolic problems, the use of a generalized slope limiter is indispensable, as was originally shown for the so-called high-resolution methods. See also [27] for a motivation of the introduction of this operator.

## 3.5 The generalized slope limiter

For piecewise-linear functions

$$ v_h|_{I_j} = \bar{v}_j + (x - x_j) v_{x,j}, $$

we define $u_h = \Lambda\Pi_h^1(v_h)$, [68], as follows:

$$ u_h|_{I_j} = \bar{v}_j + (x - x_j) m \left( v_{x,j}, \frac{\bar{v}_{j+1} - \bar{v}_j}{\Delta_j/2}, \frac{\bar{v}_j - \bar{v}_{j-1}}{\Delta_j/2} \right), $$

where the minmod function $m$ is defined by

$$ m(a_1, a_2, a_3) = \begin{cases} s \min_{1 \le n \le 3} |a_n| & \text{if } s = \operatorname{sign}(a_1) = \operatorname{sign}(a_2) = \operatorname{sign}(a_3), \\ 0 & \text{otherwise.} \end{cases} $$

Note that this projection is non-linear and can be rewritten as follows:

$$ u_{j+1/2}^{-} = \bar{v}_{j} + m \left( v_{j+1/2}^{-} - \bar{v}_{j}, \bar{v}_{j} - \bar{v}_{j-1}, \bar{v}_{j+1}^{-} - \bar{v}_{j} \right) \quad (7) $$

$$ u_{j-1/2}^{+} = \bar{v}_{j} - m \left( \bar{v}_{j} - v_{j-1/2}^{+}, \bar{v}_{j} - \bar{v}_{j-1}, \bar{v}_{j+1}^{-} - \bar{v}_{j} \right). \quad (8) $$

Next, we define a generalized slope limiter $\Lambda\Pi_h$ for general discontinuous functions. Let us denote by $v_h^1$ the L²-projection of $v_h$ into the space of piecewise-linear functions. We then define $u_h = \Lambda\Pi_h(v_h)$ on the interval $I_j$, as follows:

(i) Compute $u_{j+1/2}^{-}$ and $u_{j-1/2}^{+}$ by using (7) and (8).

(ii) If $u_{j+1/2}^{-} = v_{j+1/2}^{-}$ and $u_{j-1/2}^{+} = v_{j-1/2}^{+}$ set $u_h|_{I_j} = v_h|_{I_j}$.

(iii) If not, take $u_h|_{I_j}$ equal to $\Lambda\Pi_h^1(v_h^1)$.
Since the above generalized slope limiter enforces the sign conditions, the forward Euler step is stable, that is, we have the following result.

**Proposition 1 (The TVBM property).** Suppose that for $j = 1, \dots, N$

$$|\delta| \left( \frac{|\hat{f}(a, \cdot)|_{Lip}}{\Delta_{j+1}} + \frac{|\hat{f}(\cdot, b)|_{Lip}}{\Delta_j} \right) \le 1/2.$$

Then, if $u_h = \Lambda\Pi_h v_h$, for some $v_h$, then

$$|\bar{w}_h|_{TV(0,1)} \le |\bar{u}_h|_{TV(0,1)}.$$

## 3.6 The stability of the RKDG method

For this method, we have the following stability result.

**Theorem 2 (TVBM-stability of the RKDG method).** Let each time step $\Delta t^n$ satisfy the following CFL condition:

$$\max_{il} \left| \frac{\beta_{il}}{\alpha_{il}} \right| \Delta t^n \left( \frac{|\hat{f}(a, \cdot)|_{Lip}}{\Delta_{j+1}} + \frac{|\hat{f}(\cdot, b)|_{Lip}}{\Delta_j} \right) \le 1/2. \qquad (9)$$

Then we have

$$|\bar{u}_h^n|_{TV(0,1)} \le |u_0|_{TV(0,1)} \quad \forall n \in \mathbb{N}.$$

## 3.7 Extensions

We must point out that with the above generalized slope limiter, there is loss of accuracy when the exact solution displays critical points. It is possible, however, to modify the above limiter to completely overcome this difficulty. It is also possible to extend the above results to the multi-dimensional case if we use the $L^\infty$-norm of the local means instead of their total variation. Finally, let us point that, although there are no stability proofs for DG methods applied to non-linear hyperbolic problems, the methods work extremely well.

To show this, we show some computational results from [39]. We display the contours of the density for the so-called double Mach reflection problem. We see that both the strong shocks as well as the contacts are well approximated and that no spurious oscillations appear when we go from linear to quadratic approximations. Moreover, the contact discontinuities seem to be better approximated by using quadratic polynomials. In fact, as argued in [39], even though the use of higher degree polynomials entails a more restrictive CFL condition, the enhanced quality of the approximation more than off-sets the increase of cost per mesh point.

# 4 Concluding remarks

In this paper, we have displayed the main ingredients (discontinuous approximations, element-by-element Galerkin weak formulations, and the numerical traces) and properties (local conservativity, high-order accuracy, and dissipation or stabilization through the jumps) of the DG methods as applied to a wide variety of problem ranging from ODEs to non-linear hyperbolic systems. We have also shown some of the particularities of these methods when applied to different problems: High accuracy and parallelizability for hyperbolic problems, easy elimination of the auxiliary unknowns for elliptic problems, and robustness with respect to the Reynolds numbers for the Oseen problem.

In [40], some of the challenges that the development of these methods present are described. Let us end this paper by pointing out that a very important issue that has not been touched in this paper is the ease with which these methods can handle hp-adaptive strategies [16–18, 43, 44, 47–50, 55–58, 74, 75] and can be coupled with already existing methods [2, 29, 69]. This is already having an important impact on the development of hp-adaptive algorithms and will also facilitate the handling of multi-physics models.

**Acknowledgements** The author would like to thank Prof. Dr. Ronald H. W. Hoppe and Prof. Dr. H.-J. Bungartz for the kind invitation to deliver a plenary lecture at the Annual Meeting of GAMM at Augsburg in March of 2002. The content of paper is based on that lecture. The author would also like to thank Ryuta Suzuki for discussions leading to a better presentation of the material in this paper.
<figure><img src="image_5.png" /><figcaption>Rectangles P1, $\Delta x = \Delta y = 1/480$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Rectangles P2, $\Delta x = \Delta y = 1/480$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4** Euler equations of gas dynamics: Double Mach reflection problem. Isolines of the density around the double Mach stems. Quadratic polynomials on squares $\Delta x = \Delta y = 1/240$ (top); linear polynomials on squares $\Delta x = \Delta y = 1/480$ (middle); and quadratic polynomials on squares $\Delta x = \Delta y = 1/480$ (bottom).</figcaption></figure>
References

[1] V. Aizinger, C. Dawson, B. Cockburn, and P. Castillo, Local discontinuous Galerkin method for contaminant transport, Adv. Water Resour. **24**, 73–87 (2000).

[2] P. Alotto, A. Bertoni, I. Perugia, and D. Schötzau, Discontinuous finite element methods for the simulation of rotating electrical machines, COMPEL - Int. J. Comput. Math. Electr. Electron. Eng. (UK) **20**, 448–462 (2001).

[3] D. Arnold, F. Brezzi, B. Cockburn, and D. Marini, Discontinuous Galerkin methods for elliptic problems, in: Discontinuous Galerkin Methods. Theory, Computation and Applications, edited by B. Cockburn, G. Karniadakis, and C. W. Shu, Lecture Notes in Computational Science and Engineering, Vol. 11 (Springer Verlag, Berlin, 2000), pp. 89–101.

[4] D. Arnold, F. Brezzi, B. Cockburn, and D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. **39**, 1749–1779 (2001).

[5] F. Baaijens, Application of low-order discontinuous Galerkin methods to the analysis of viscoelastic flows, J. Non-Newton. Fluid Mech. **52**, 37–57 (1994).

[6] I. Babuška, C. Baumann, and J. Oden, A discontinuous hp-finite element method for diffusion problems: 1-D analysis, Comput. Appl. Math. **37**, 103–122 (1999).

[7] I. Babuška and M. Zlámal, Nonconforming elements in the finite element method with penalty, SIAM J. Numer. Anal. **10**, 863–875 (1973).

[8] A. Bahhar, J. Baranger, and D. Sandri, Galerkin discontinuous approximation of the transport equation and viscoelastic fluid flow on quadrilaterals, Numer. Methods Partial Differ. Equ. **14**, 97–114 (1998).

[9] G. Baker, W. Jureidini, and O. Karakashian, Piecewise solenoidal vector fields and the Stokes problem, SIAM J. Numer. Anal. **27**, 1466–1485 (1990).

[10] F. Bassi and S. Rebay, A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations, J. Comput. Phys. **131**, 267–279 (1997).

[11] F. Bassi and S. Rebay, High-order accurate discontinuous finite element solution of the 2D Euler equations, J. Comput. Phys. **138**, 251–285 (1997).

[12] F. Bassi, S. Rebay, G. Mariotti, S. Pedinotti, and M. Savini, A high-order accurate discontinuous finite element method for inviscid and viscous turbomachinery flows, in: 2nd European Conference on Turbomachinery Fluid Dynamics and Thermodynamics, edited by R. Decuypere and G. Dibelius. Technologisch Instituut, Antwerpen, Belgium, March 5–7 1997, pp. 99–108.

[13] C. Baumann and J. Oden, A discontinuous hp-finite element method for the Navier-Stokes equations, in: Proceedings of the 10th International Conference on Finite Element in Fluids, 1998.

[14] C. Baumann and J. Oden, A discontinuous hp-finite element method for the solution of the Euler equation of gas dynamics, in: Proceedings of the 10th International Conference on Finite Element in Fluids, 1998.

[15] C. Baumann and J. Oden, A discontinuous hp-finite element method for convection-diffusion problems, Comput. Methods Appl. Mech. Eng. **175**, 311–341 (1999).

[16] K. Bey, An hp-adaptive discontinuous Galerkin method for hyperbolic conservation laws, PhD thesis, The University of Texas at Austin (1994).

[17] K. Bey and J. Oden, hp-version discontinuous Galerkin methods for hyperbolic conservation laws, Comput. Methods Appl. Mech. Eng. **133**, 259–286 (1996).

[18] R. Biswas, K. Devine, and J. Flaherty, Parallel, adaptive finite element methods for conservation laws, Appl. Numer. Math. **14**, 255–283 (1994).

[19] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods (Springer Verlag, New York, 1991).

[20] F. Brezzi, G. Manzini, D. Marini, P. Pietra, and A. Russo, Discontinuous Galerkin approximations for elliptic problems, Numer. Methods Partial Differ. Equ. **16**, 365–378 (2000).

[21] F. Carranza, B. Fang, and R. Haber, An adaptive discontinuous Galerkin model for coupled viscoplastic crack growth and chemical transport, in: Discontinuous Galerkin Methods. Theory, Computation and Applications, edited by B. Cockburn, G. Karniadakis, and C. W. Shu, Lecture Notes in Computational Science and Engineering, Vol. 11 (Springer Verlag, Berlin, 2000), pp. 277–283.

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

[23] G. Chavent and B. Cockburn, The local projection $P^0 P^1$-discontinuous-Galerkin finite element method for scalar conservation laws, RAIRO Modél. Math. Anal. Numér. **23**, 565–592 (1989).

[24] Z. Chen, B. Cockburn, C. Gardner, and J. Jerome, Quantum hydrodynamic simulation of hysteresis in the resonant tunneling diode, J. Comput. Phys. **117**, 274–280 (1995).

[25] Z. Chen, B. Cockburn, J. Jerome, and C. W. Shu, Mixed-RKDG finite element methods for the 2-D hydrodynamic model for semiconductor device simulation, VLSI Des. **3**, 145–158 (1995).

[26] B. Cockburn, Discontinuous Galerkin methods for convection-dominated problems, in: High-Order Methods for Computational Physics, edited by T. Barth and H. Deconink, Lecture Notes in Computational Science and Engineering, Vol. 9 (Springer Verlag, Berlin, 1999), pp. 69–224.

[27] B. Cockburn, Devising discontinuous Galerkin methods for non-linear hyperbolic conservation laws, J. Comput. Appl. Math. **128**, 187–204 (2001).

[28] B. Cockburn and C. Dawson, Some extensions of the local discontinuous Galerkin method for convection-diffusion equations in multidimensions, in: The Proceedings of the Conference on the Mathematics of Finite Elements and Applications: MAFELAP X, edited by J. Whiteman (Elsevier, Amsterdam, 2000), pp. 225–238.

[29] B. Cockburn and C. Dawson, Approximation of the velocity by coupling discontinuous Galerkin and mixed finite element methods for flow problems, Special issue: Locally Conservative Numerical Methods for Flow in Porous Media Comput. Geosci. **6**, 502–522 (2002).

[30] B. Cockburn, S. Hou, and C. W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: The multidimensional case, Math. Comput. **54**, 545–581 (1990).
[31] B. Cockburn, G. Kanschat, I. Perugia, and D. Schötzau, Superconvergence of the local discontinuous Galerkin method for elliptic problems on Cartesian grids, SIAM J. Numer. Anal. **39**, 264–285 (2001).

[32] B. Cockburn, G. Kanschat, and D. Schötzau, Local discontinuous Galerkin methods for the Oseen equations, Math. Comput., to appear.

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

[34] B. Cockburn, G. Karniadakis, and C. W. Shu, The development of discontinuous Galerkin methods, in: Discontinuous Galerkin Methods. Theory, Computation and Applications, edited by B. Cockburn, G. Karniadakis, and C. W. Shu, Lecture Notes in Computational Science and Engineering, Vol. 11 (Springer Verlag, Berlin, 2000), pp. 3–50.

[35] B. Cockburn, S. Lin, and C. W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One dimensional systems, J. Comput. Phys. **84**, 90–113 (1989).

[36] B. Cockburn and C. W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws II: General framework, Math. Comput. **52**, 411–435 (1989).

[37] B. Cockburn and C. W. Shu, The Runge-Kutta local projection $P^1$-discontinuous Galerkin method for scalar conservation laws, RAIRO Modél. Math. Anal. Numér. **25**, 337–361 (1991).

[38] B. Cockburn and C. W. Shu, The local discontinuous Galerkin method for time-dependent convection-diffusion systems, SIAM J. Numer. Anal. **35**, 2440–2463 (1998).

[39] B. Cockburn and C. W. Shu, The Runge-Kutta discontinuous Galerkin finite element method for conservation laws V: Multidimensional systems, J. Comput. Phys. **141**, 199–224 (1998).

[40] B. Cockburn and C. W. Shu, Runge-Kutta discontinuous Galerkin methods for convection-dominated problems, J. Sci. Comput. **16**, 173–261 (2001).

[41] C. Dawson, V. Aizinger, and B. Cockburn, Local discontinuous Galerkin methods for problems in contaminant transport, in: Discontinuous Galerkin Methods. Theory, Computation and Applications, edited by B. Cockburn, G. Karniadakis, and C. W. Shu, Lecture Notes in Computational Science and Engineering, Vol. 11 (Springer Verlag, Berlin, 2000), pp. 309–314.

[42] M. Delfour, W. Hager, and F. Trochu, Discontinuous Galerkin methods for ordinary differential equations, Math. Comput. **36**, 455–473 (1981).

[43] K. Devine and J. Flaherty, Parallel adaptive hp-refinement techniques for conservation laws, Appl. Numer. Math. **20**, 367–386 (1996).

[44] K. Devine, J. Flaherty, R. Loy, and S. Wheat, Parallel partitioning strategies for the adaptive solution of conservation laws, in: Modeling, Mesh Generation, and Adaptive Numerical Methods for Partial Differential Equations, edited by I. Babuška, W. Henshaw, J. Hopcroft, J. Oliger, and T. Tezduyar (Springer Verlag, New York, 1995), pp. 215–242.

[45] J. Douglas, jr. and T. Dupont, Interior Penalty Procedures for Elliptic and Parabolic Galerkin Methods, Lecture Notes in Physics, Vol. 58 (Springer-Verlag, Berlin, 1976).

[46] G. Engel, K. Garikipati, T. Hughes, M. Larson, L. Mazzei, and R. Taylor, Continuous/discontinuous finite element approximations of fourth-order elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity, Comput. Methods Appl. Mech. Eng. **191**, 3669–3750 (2002).

[47] J. Flaherty, R. Loy, M. Shephard, and J. Teresco, Software for the parallel adaptive solution of conservation laws by a discontinuous Galerkin method, in: Discontinuous Galerkin Methods. Theory, Computation and Applications, edited by B. Cockburn, G. Karniadakis, and C. W. Shu, Lecture Notes in Computational Science and Engineering, Vol. 11 (Springer Verlag, Berlin, 2000), pp. 113–123.

[48] J. Flaherty, R. Loy, C. Özturan, M. Shephard, B. Szymanski, J. Teresco, and L. Ziantz, Parallel structures and dynamic load balancing for adaptive finite element computation, Appl. Numer. Math. **26**, 241–265 (1998).

[49] J. Flaherty, R. Loy, M. Shephard, M. Simone, B. Szymanski, J. Teresco, and L. Ziantz, Distributed octree data structures and local refinement method for the parallel solution of three-dimensional conservation laws, in: Grid Generation and Adaptive Algorithms, edited by M. Bern, J. Flaherty, and M. Luskin, The IMA Volumes in Mathematics and its Applications, Vol. 113. Institute for Mathematics and its Applications (Springer, Minneapolis, 1999), pp. 113–134.

[50] J. Flaherty, R. Loy, M. Shephard, B. Szymanski, J. Teresco, and L. Ziantz, Adaptive local refinement with octree load-balancing for the parallel solution of three-dimensional conservation laws, J. Parallel Distrib. Comput. **47**, 139–152 (1997).

[51] M. Fortin and A. Fortin, New approach for the finite element method simulation of viscoelastic flows, J. Non-Newton. Fluid Mech. **32**, 295–310 (1989).

[52] P. Gremaud and J. Matthews, On the computation of steady hopper flows I. Stress determination for coulomb materials, J. Comput. Phys. **166**, 63–83 (2001).

[53] P. A. Gremaud and J. Matthews, Simulation of gravity flow of granular materials in silos, in: Discontinuous Galerkin Methods. Theory, Computation and Applications, edited by B. Cockburn, G. Karniadakis, and C. W. Shu, Lecture Notes in Computational Science and Engineering, Vol. 11 (Springer Verlag, Berlin, 2000), pp. 125–134.

[54] P. Hansbo and M. Larson, Discontinuous finite element methods for incompressible and nearly incompressible elasticity by use of Nitsche's method, Comput. Methods Appl. Mech. Eng. **191**, 1895–1908 (2002).

[55] P. Houston, C. Schwab, and E. Süli, Discontinuous hp finite element methods for advection-diffusion problems, Tech. Rep. NA 00-15, Oxford University Computing Laboratory (Oxford, UK, 2000)

[56] P. Houston, C. Schwab, and E. Süli, Stabilized hp-finite element methods for hyperbolic problems, SIAM J. Numer. Anal. **37**, 1618–1643 (2000).

[57] P. Houston, C. Schwab, and E. Süli, hp-adaptive discontinuous Galerkin finite element methods for hyperbolic problems, SIAM J. Math. Anal. **23**, 1226–1252 (2001).

[58] P. Houston and E. Süli, Stabilized hp-finite element approximation of partial differential equations with nonnegative characteristic form, Computing **60**, 99–119 (2001).

[59] C. Hu, O. Lepsky, and C. W. Shu, The effect of the least square procedure for discontinuous Galerkin methods for Hamilton-Jacobi equations, in: Discontinuous Galerkin Methods. Theory, Computation and Applications, edited by B. Cockburn, G. Karniadakis, and C. W. Shu, Lecture Notes in Computational Science and Engineering, Vol. 11 (Springer Verlag, Berlin, 2000), pp. 343–348.
[60] C. Hu and C. W. Shu, A discontinuous Galerkin finite element method for Hamilton-Jacobi equations, SIAM J. Sci. Comput. **21**, 666-690 (1999).

[61] O. Karakashian and W. Jureidini, A nonconforming finite element method for the stationary Navier-Stokes equations, SIAM J. Numer. Anal. **35**, 93-120 (1998).

[62] O. Lepsky, C. Hu, and C. W. Shu, Analysis of the discontinuous Galerkin method for Hamilton-Jacobi equations, Appl. Numer. Math. **33**, 423-434 (2000).

[63] P. Lesaint and P. Raviart, On a finite element method for solving the neutron transport equation, in: Mathematical Aspects of finite Elements in Partial Differential Equations, edited by C. deBoor (Academic Press, New York, 1974), pp. 89-145.

[64] I. Lomtev and G. Karniadakis, Simulations of viscous supersonic flows on unstructured hp-meshes, in: Proceedings of the 35th Aerospace Sciences Meeting, Reno, Nevada (1997), AIAA-97-0754.

[65] I. Lomtev and G. Karniadakis, A discontinuous Galerkin method for the Navier-Stokes equations, Int. J. Numer. Methods Fluids **29**, 587-603 (1999).

[66] I. Lomtev, C. Quillen, and G. Karniadakis, Spectral/hp methods for viscous compressible flows on unstructured 2D meshes, J. Comput. Phys. **144**, 325-357 (1998).

[67] J. Oden, I. Babuška, and C. Baumann, A discontinuous hp-finite element method for diffusion problems, J. Comput. Phys. **146**, 491-519 (1998).

[68] S. Osher, Convergence of generalized MUSCL schemes, SIAM J. Numer. Anal. **22**, 947-961 (1984).

[69] I. Perugia and D. Schötzau, On the coupling of local discontinuous Galerkin and conforming finite element methods, J. Sci. Comput. **16**, 411-433 (2002).

[70] B. Rivière, M. Wheeler, and V. Girault, Improved energy estimates for interior penalty, constrained and discontinuous Galerkin methods for elliptic problems. Part I, Comput. Geosci. 1999 (3), 337-360 (1999).

[71] D. Schötzau, C. Schwab, and A. Toselli, hp-DGFEM for incompressible flows, SIAM J. Numer. Anal. **40**, 2171-2194 (2003).

[72] C. W. Shu and J. Yan, A local discontinuous Galerkin method for KdV type equations, SIAM J. Numer. Anal. **40**, 769-791 (2002) (electronic).

[73] C. W. Shu and J. Yan, Local discontinuous Galerkin for partial differential equations with higher order derivatives, J. Sci. Comput. **17**, 27-47 (2002).

[74] E. Süli, P. Houston, and C. Schwab, hp-finite element methods for hyperbolic problems, in: The Proceedings of the Conference on the Mathematics of Finite Elements and Applications: MAFELAP X, edited by J. Whiteman (Elsevier, Amsterdam, 2000), pp. 143-162.

[75] E. Süli, C. Schwab, and P. Houston, hp-DGFEM for partial differential equations with non-negative characteristic form, in: Discontinuous Galerkin Methods. Theory, Computation and Applications, edited by B. Cockburn, G. Karniadakis, and C. W. Shu, Lecture Notes in Computational Science and Engineering, Vol. 11 (Springer Verlag, Berlin, 2000), pp. 221-230.

[76] A. Toselli, hp-discontinuous Galerkin approximations for the Stokes problem, Math. Models Methods Appl. Sci. **12**, 1565-1616 (2002).

[77] T. Warburton and G. Karniadakis, A discontinuous Galerkin method for the viscous MHD equations, J. Comput. Phys. **152**, 1-34 (1999).