Discontinuous Galerkin Method

§1 Introduction

Discontinuous Galerkin method is a finite element method using completely discontinuous piecewise polynomial space for the numerical solution as well as the test functions. It was first developed in 1973 by Reed and Hill [17] to solve the neutron transport equation, which is a linear hyperbolic PDE. LeSaint and Raviart [14] published its first mathematical analysis in 1974, proving that the method converges with a rate (Δx)^{k} in the L^{2}-norm for general triangulations using polynomials of degree k and with a rate (Δx)^{k+1} for tensor products of polynomials of degree k in each variable on Cartesian meshes. In 1986, Johnson and Pitkaränta [13] proved a rate of convergence of (Δx)^{k+1/2} for general triangulations and in 1991, Peterson [16] numerically confirmed this rate to be optimal.

The discontinuous Galerkin method for solving nonlinear hyperbolic problems started with a series of papers by Cockburn, Shu and their collaborators starting in the late 80s, [5, 6, 7, 8], which incorporated important ideas of high resolution finite difference and finite volume methods such as monotone fluxes or approximate Riemann solvers, high order TVD Runge-Kutta time discretizations, and nonlinear limiters for discontinuities in the solutions, into the framework of discontinuous Galerkin methods.

The original discontinuous Galerkin method is designed solely for solving hyperbolic problems. However recently the method has been successfully generalized to solve other types of PDEs such as the parabolic and elliptic problems [9, 4].

The applications of the discontinuous Galerkin method in the literature include areas as diverse as meteorology, weather-forecasting, oceanography, gas dynamics, turbomachinery, turbulent flows, granular flows, oil recovery simulation, modeling of shallow waters, transport of contaminant in porous media, viscoelastic flows, semiconductor device simulation, magnetohydrodynamics, and electro-magnetism, among many others.

In 1998 the first International Symposium on Discontinuous Galerkin methods was held in Newport, Rhode Island, USA. The symposium proceedings was published by Springer and contained papers in algorithm development, analysis, application, and parallel implemenation issues for the method, as well as an overview paper written by the organizers, Cockburn, Karniadakis and Shu [10]. This volume can be used as a genearl reference if more details about the method is desired.

§2 A simple example which carries most of the ideas

We shall first describe the essential points of the discontinuous Galerkin method using the one dimensional linear equation

u_{t} − u_{x} = 0  (1)

as an example. We shall solve (1) for x ∈ [0,1] with periodic boundary conditions and with an initial condition u(x, 0) = sin(2πx).

Let’s denote I_{j} =[x_{j− \frac{1}{2}} , x_{j+ \frac{1}{2}} ], for j = 1, ..., N, as a mesh for [0,1], where x_{\frac{1}{2}} = 0 and x_{N+ \frac{1}{2}} = 1. We denote the center of each cell by x_{j} = \frac{1}{2} (x_{j− \frac{1}{2}} + x_{j+ \frac{1}{2}} ) and the size of each cell by Δx_{j} = x_{j+ \frac{1}{2}} − x_{j− \frac{1}{2}} . The cells do not need to be uniform but for simplicity we will perform calculations only with uniform meshes and will denote the uniform mesh size by Δx.

If we multiply (1) by an arbitrary test function v(x), integrate over the interval I_{j}, and integrate by parts, we get

∫_{I_{j}} u_{t}vdx + ∫_{I_{j}} uv_{x}dx − u(x_{j+ \frac{1}{2}} , t)v(x_{j+ \frac{1}{2}} ) + u(x_{j− \frac{1}{2}} , t)v(x_{j− \frac{1}{2}} ) = 0.  (2)

This is the starting point for designing the discontinuous Galerkin method. We replace both the solution u and the test function v by piecewise polynomials of degree at most k. That is, u, v ∈ V_{Δx} where

V_{Δx} = {v : v is a polynomial of degree at most k for x ∈ I_{j}, j = 1, ..., N} .  (3)

With this choice, there is an ambiguity in (2) in the last two terms involving the boundary values at x_{j± \frac{1}{2}} , as both the solution u and the test function v are discontinuous exactly at these boundary points. This is a nuisance but also an opportunity: one could cleverly design these terms so that the resulting numerical method is stable and accurate. To motivate the ideas, let’s look at the simplest case k = 0. That is, the solution as well as the test functions are piecewise constants. If we denote by u_{j} the value of u (which is constant in each cell) in the cell I_{j}, (2) would become the familiar first order upwind finite volume scheme

\frac{d}{dt}u_{j} − \frac{1}{Δx_{j}} (u_{j+1} − u_{j}) = 0

if we perform the following in (2):

1. Replace the boundary terms u(x_{j± \frac{1}{2}} , t) by single valued numerical fluxes û_{j± \frac{1}{2}} = û(u_{j± \frac{1}{2}}^{−} , u_{j± \frac{1}{2}}^{+} ). These fluxes in general depend both on the left limit (e.g. u_{j+ \frac{1}{2}}^{−} = lim_{x→x_{j+ \frac{1}{2}}^{−}} u(x, t)) and on the right limit (e.g. u_{j+ \frac{1}{2}}^{+} = lim_{x→x_{j+ \frac{1}{2}}^{+}} u(x, t)). For the equation (1), the flux û_{j+ \frac{1}{2}} is taken as u_{j+ \frac{1}{2}}^{+} according to upwinding.

2. Replace the test function v at the boundaries by the values taken from inside the cell I_{j}, namely v_{j+ \frac{1}{2}}^{−} and v_{j− \frac{1}{2}}^{+} .

The scheme now becomes: find u ∈ V_{Δx} such that, for all test functions v ∈ V_{Δx},

∫_{I_{j}} u_{t}vdx + ∫_{I_{j}} uv_{x}dx − û_{j+ \frac{1}{2}} v_{j+ \frac{1}{2}}^{−} + û_{j− \frac{1}{2}} v_{j− \frac{1}{2}}^{+} = 0  (4)

where the numerical flux û_{j+ \frac{1}{2}} = u_{j+ \frac{1}{2}}^{+} .

After picking a local basis and inverting a local (k + 1) × (k + 1) mass matrix (by hand), the scheme (4) can be written as

\frac{d}{dt}u_{j} + \frac{1}{Δx_{j}} (Au_{j} + Bu_{j+1}) = 0  (5)

���� 1: L^{2} and L^{∞} errors and orders of accuracy for the discontinuous Galerkin method (4) applied to the linear equation (1) with an initial condition u(x, 0) = sin(2πx), t = 1. Third order Runge-Kutta in time.�

k = 1  k = 2
Δx L^{2} error order L^{∞} error order L^{2} error order L^{∞} error order
2π/20 4.63E-03  — 1.34E-02  — 1.14E-04  — 5.10E-04  —
2π/40 1.09E-03  2.08 3.75E-03  1.84 1.42E-05  3.00 6.44E-05  2.98
2π/80 2.69E-04  2.02 9.84E-04  1.93 1.77E-06  3.00 8.08E-06  3.00
2π/160 6.69E-05  2.01 2.52E-04  1.97 2.21E-07  3.00 1.01E-06  3.00

where u_{j} is a small vector of length k + 1 containing the coefficients of the solution u in the local basis inside cell I_{j}, and A and B are (k + 1) × (k + 1) constant matrices which can be computed once and for all and stored at the beginning of the code.

Scheme (5) can then be easily discretized in time by the nonlinearly stable high order Runge-Kutta methods in [20]. We use the third order version in [20] for the accuracy tests.

We remark that the method (5) is extremely simple to code and easy to parallelize.

For illustration purpose we show in Table 1 the L^{2} and L^{∞} errors and numerically observed orders of accuracy for the two cases k = 1 and 2 (piecewise linear and piecewise quadratic cases) for t = 1 (after one time period). We can clearly see that an order of k + 1 is achieved.

§3 DG method for scalar 1D hyperbolic problems

In this section, we give more details of the DG method for the non-linear 1D hyperbolic conservation laws.

u_{t} + f(u)_{x} = 0, in (0, 1) × (0, T),  u(x, 0) = u_{0}(x), ∀ x ∈ (0, 1),  (6)

with periodic boundary conditions.

As is well known, the main difficulty of a numerical solution to (6) is the appearance of shocks even if the initial condition u_{0}(x) is smooth. A good scheme for (6) would hopefully have the following properties:

• It is locally conservative.

• It is high order accurate in smooth regions of the solution.

• It has sharp and monotone (non-oscillatory) shock transitions.

• The numerical solution should be self-similar, that is, it should remain invariant when both space x and time t are scaled by the same constant. Notice that this self-similarity is an important property held by the exact solution of (6) and should be maintained by the numerical solution whenever possible. We will discuss in detail RKDG method in terms of these requirements.

§3.1 The DG space discretization

We seek an approximation u_{h} whose restriction to each element I_{j} = (x_{j−1/2}, x_{j+1/2}) is, for each value of the time variable, an element of the local space U(I_{j}); typically, U(I_{j}) is the space of polynomials of degree at most k ≥ 0. A reasonable way to define the initial data u_{h}(·, 0) on the element I_{j} is to take the L^{2}-projection of u_{0} on the local space U(I_{j}), that is, for all v_{h} ∈ U(I_{j}),

∫_{I_{j}} u_{h}(x, 0) v_{h}(x) dx = ∫_{I_{j}} u_{0}(x) v_{h}(x) dx.  (7)

To determine the approximate solution for t > 0, we enforce the non-linear conservation law element-by-element by means of a Galerkin method. Thus, on each interval I_{j} = (x_{j−1/2}, x_{j+1/2}), we require that, for all v_{h} ∈ U(I_{j}),

∫_{I_{j}} (u_{h}(x, t))_{t} v_{h}(x) dx − ∫_{I_{j}} f(u_{h}(x, t)) (v_{h}(x))_{x} dx + f̂(u_{h}(·, t)) v_{h}|_{x_{j−1/2}}^{x_{j+1/2}} = 0,  (8)

where f̂(u_{h}) is the numerical flux. Note that u_{h} is a well defined function since there are as many equations per element as unknowns. The integral

∫_{I_{j}} f(u_{h}(x, t)) (v_{h}(x))_{x} dx,

could either be computed exactly or approximately by using suitable numerical quadratures or other methods.

To define the numerical flux we can borrow from the successful development in the socalled high resolution finite volume schemes for conservation laws. The numerical flux depends on the two values of the approximate solution u_{h} at the discontinuities, that is,

f̂(u_{h})(x_{j+1/2}) = f̂(u_{h}(x_{j+1/2}^{−}), u_{h}(x_{j+1/2}^{+})).

To guarantee that, when the approximate solution u_{h} is piecewise-constant, the DG space discretization gives rise to a monotone finite volume scheme, which is stable and convergent to the exact entropy solution, we only need to require that the flux is monotone, namely a ↦ f̂(a, ·) is non-decreasing and b ↦ f̂(·, b) is non-increasing. A few examples of numerical fluxes satisfying the above properties are the following:

(i) The Godunov flux:

f̂^{G}(a, b) = \begin{cases}min_{a≤u≤b} f(u),  if a ≤ b, \\ max_{b≤u≤a} f(u),  otherwise. \end{cases}

(ii) The Engquist-Osher flux:

f̂^{EO}(a, b) = ∫_{0}^{b} min(f′(s), 0) ds + ∫_{0}^{a} max(f′(s), 0) ds + f(0);

(iii) The Lax-Friedrichs flux:

f̂^{LF} (a, b) = \frac{1}{2} [f(a) + f(b) − C (b − a)],  C = max_{inf u^{0}(x)≤s≤sup u^{0}(x)} |f′(s)|.

This completes the definition of the DG space approximation.

Several comments about this DG space discretization are in order:

• In the linear case f(u) = c u, all the above numerical fluxes coincide with the so-called upwind numerical flux, namely, f̂(a, b) = c a if c ≥ 0, and f̂(a, b) = c b if c < 0. The first DG method was introduced in 1973 by Reed and Hill [17] for a linear equation modeling transport of neutrons; they used the upwind numerical flux, even though the concept of monotone schemes and numerical fluxes had not been introduced yet.

• We mentioned briefly in the introduction section the theoretical analysis carried out for the DG method in linear case. The non-linear case is much more difficult to study. So far there is no convergence analysis of the DG space discretization for non-linear scalar conservation laws with non-smooth solutions except for the piecewise-constant case. Notice that the DG space discretization is linear (i.e. the scheme is linear for linear PDEs), hence it is expected to be oscillatory for problems with shocks except for the piecewise-constant case (the well known “Godunov Theorem”. However, comparing with other high order linear schemes such as finite difference and finite volume schemes, the method of lines version of the DG space discretization (and also certain implicit time discretization of it, such as backward Euler or Crank-Nicholson) satisfies a remarkably stronger provable stability property: a local cell entropy inequality for the square entropy as proven by Jiang and Shu in 1994 [12]. This result is valid for any order of accuracy and any triangulation in any spatial dimensions. It trivially implies a L^{2} norm non-increasing with time for the numerical solution and enforces any limit solution to be the correct entropy solution when f(u) is convex. The only comparable result for finite difference or finite volume schemes is for one space dimension, second order accurate only, nonlinear schemes. However, unlike for the linear case, numerical evidence suggests that the DG method with explicit time stepping might be unstable in the general case with strong shocks.

• For non-linear problems, the best choice of numerical flux is the Godunov flux f̂^{G} since it is well-known that this is the numerical flux that produces the smallest amount of artificial viscosity. The local Lax-Friedrichs flux

f̂^{LLF} (a, b) = \frac{1}{2} [f(a) + f(b) − C(b − a)],  C = max_{min(a,b)≤s≤max(a,b)} |f′(s)|,

produces more artificial viscosity than the Godunov flux, but their performances are remarkably similar. Of course, if f is too complicated, we can always use the simple Lax-Friedrichs flux. Numerical experience suggests that as the degree k of the approximate solution increases, the choice of the numerical flux does not have a significant impact on the quality of the approximations.

• In the general case, local mass matrix can be easily inverted, by means of a symbolic manipulator, for example, since its order is equal to the dimension of the local spaces. We thus can always obtain a system of ordinary differential equations for the degrees of freedom that we can write as follows:

\frac{d}{dt} u_{h} = L_{h}(u_{h}),  in (0, T),  u_{h}(t = 0) = P_{h}u_{0},

where P_{h} denotes the L^{2}-projection and the function L_{h}(u_{h}) is, of course, the approximation to −f(u)_{x} provided by the DG-space discretization.

§3.2 The RK time discretization

We discretize in time our system of ordinary differential equations by using the following RK method:

1. Set u_{h}^{(0)} = u_{h}^{n};

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

u_{h}^{(i)} = ∑_{l=0}^{i−1} α_{il} w_{h}^{il}, w_{h}^{il} = u_{h}^{(l)} + \frac{β_{il}}{α_{il}} Δt^{n}L_{h}(u_{h}^{(l)});

3. Set u_{h}^{n+1} = u_{h}^{K}. When the coefficients α_{il} and β_{il} are all non-negative, the method is a simple convex combination of forward Euler steps and its stability follows from that of the forward Euler. These TVD Runge-Kutta methods were first developed in [20]. Some of the RK methods in this class are displayed in the Table 2.

���� 2: TVD-RK time discretization parameters.�

order  α_{il}  β_{il}  max{β_{il}/α_{il}}
2  1  1  1
\frac{1}{2} \frac{1}{2}  0 \frac{1}{2}
1  1
3  \frac{3}{4} \frac{1}{4}  0 \frac{1}{4}  1
\frac{1}{3} 0 \frac{2}{3}  0 0 \frac{2}{3}

It is essential to carry out a von Neumann stability analysis of the method for the linear case f(u) = c u in order to know for what values of the number CFL_{L^{2}}, the condition

|c|\frac{Δt}{Δx} ≤ CFL_{L^{2}},

ensures its L^{2}-stability. This condition has to be respected even for non-linear functions f since only under this condition the round-off errors are not amplified.

For example, for DG discretizations using polynomials of degree k and a k + 1 stage RK method of order k + 1 (which give rise to an (k + 1)-th order accurate method), we can take in practice

CFL_{L^{2}} = \frac{1}{2k + 1}.

Indeed, this can be trivially proven for k = 0 and was proven for k = 1 in [?]. Moreover, for k ≥ 2, the number \frac{1}{2k+1} is less than 5% smaller than numerically-obtained estimates of CFL_{L^{2}}. In Table 3, we display these CFL_{L^{2}} numbers for a wide variety of time and space discretizations. The symbol ‘⋆’ indicates that the method is unstable when the ratio Δt/Δx is held constant.

���� 3: The CFL_{L^{2}} numbers for polynomials of degree k and RK methods of order ν.�

\begin{matrix}k & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\ ν = 1 & 1.000 & ⋆ & ⋆ & ⋆ & ⋆ & ⋆ & ⋆ & ⋆ & ⋆ \\ ν = 2 & 1.000 & 0.333 & ⋆ & ⋆ & ⋆ & ⋆ & ⋆ & ⋆ & ⋆ \\ ν = 3 & 1.256 & 0.409 & 0.209 & 0.130 & 0.089 & 0.066 & 0.051 & 0.040 & 0.033 \\ ν = 4 & 1.392 & 0.464 & 0.235 & 0.145 & 0.100 & 0.073 & 0.056 & 0.045 & 0.037 \\ ν = 5 & 1.608 & 0.534 & 0.271 & 0.167 & 0.115 & 0.085 & 0.065 & 0.052 & 0.042 \\ ν = 6 & 1.776 & 0.592 & 0.300 & 0.185 & 0.127 & 0.093 & 0.072 & 0.057 & 0.047 \\ ν = 7 & 1.977 & 0.659 & 0.333 & 0.206 & 0.142 & 0.104 & 0.080 & 0.064 & 0.052 \\ ν = 8 & 2.156 & 0.718 & 0.364 & 0.225 & 0.154 & 0.114 & 0.087 & 0.070 & 0.057 \\ ν = 9 & 2.350 & 0.783 & 0.396 & 0.245 & 0.168 & 0.124 & 0.095 & 0.076 & 0.062 \\ ν = 10 & 2.534 & 0.844 & 0.428 & 0.264 & 0.182 & 0.134 & 0.103 & 0.082 & 0.067 \\ ν = 11 & 2.725 & 0.908 & 0.460 & 0.284 & 0.195 & 0.144 & 0.111 & 0.088 & 0.072 \\ ν = 12 & 2.911 & 0.970 & 0.491 & 0.303 & 0.209 & 0.153 & 0.118 & 0.094 & 0.077 \end{matrix}

§3.3 Limiters

Even though the RKDG method can be applied to problems with smooth solutions as well as most linear problems even with discontinuous solutions, for nonlinear problems with strong discontinuities a nonlinear limiter is usually needed to enforce stability and to control spurious oscillations. We now describe the limiter operator ΛΠ_{h}; set u_{h} = ΛΠ_{h}v_{h}. We begin by noting that for piecewise linear approximate solutions, that is,

v_{h}|_{I_{j}} = \overline{v}_{j} + (x − x_{j}) v_{x,j},

van Leer [21, 22] introduced the following slope limiter in the construction of his MUSCL schemes:

u_{h}|_{I_{j}} = \overline{v}_{j} + (x − x_{j}) m (v_{x,j}, \frac{\overline{v}_{j+1} − \overline{v}_{j}}{Δ_{j}} , \frac{\overline{v}_{j} − \overline{v}_{j−1}}{Δ_{j}}),

where the minmod function m is defined as follows:

m (a_{1}, a_{2}, a_{3}) = \begin{cases}s min_{1≤n≤3} | a_{n} | if s = sign(a_{1}) = sign(a_{2}) = sign(a_{3}), \\ 0 otherwise. \end{cases}

We use a less restrictive slope limiter, denoted by ΛΠ_{h}^{1}, due to Osher [15], which is defined as follows:

u_{h}|_{I_{j}} = \overline{v}_{j} + (x − x_{j}) m (v_{x,j}, \frac{\overline{v}_{j+1} − \overline{v}_{j}}{Δ_{j}/2} , \frac{\overline{v}_{j} − \overline{v}_{j−1}}{Δ_{j}/2} ),

which can be rewritten as follows:

u_{j+1/2}^{−} = \overline{v}_{j} + m ( v_{j+1/2}^{−} − \overline{v}_{j}, \overline{v}_{j} − \overline{v}_{j−1}, \overline{v}_{j+1} − \overline{v}_{j})  (9)
u_{j−1/2}^{+} = \overline{v}_{j} − m ( \overline{v}_{j} − v_{j−1/2}^{+}, \overline{v}_{j} − \overline{v}_{j−1}, \overline{v}_{j+1} − \overline{v}_{j}).  (10)

For general functions v_{h}, we can define a generalized slope limiter ΛΠ_{h} in a very simple way. To do that, let us denote by v_{h}^{1} the L^{2}-projection of v_{h} into the space of piecewise-linear functions. We then define u_{h} = ΛΠ_{h}(v_{h}) on the interval I_{j}, as follows:

(i) Compute u_{j+1/2}^{−} and u_{j−1/2}^{+} by using (9) and (10), (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 ΛΠ_{h}^{1}(v_{h}^{1}).

Let us discuss some important points about this limiter:

• The above recipe is remarkably simple as it can be applied to any type of approximate solution v_{h} by using the minmod function m at most three times per element.

• This generalized slope limiter can be efficiently parallelized since to compute it on the element I_{j}, the only information needed which is not associated with this element are the means on the two neighboring elements.

• The original high-order accuracy of the RKDG method is maintained away from smooth extrema (places where the first derivative of the solution is zero) but is lost near them. Thus, methods using a TVD limiter is only first order accurate in the maximum norm and second order in L^{1}. To correct this deficiency, we can use the idea in Shu [18] and modify the definition of the generalized slope limiters by simply replacing the minmod function m by the corrected minmod function \overline{m}_{j} defined as follows:

\overline{m}_{j} (a_{1}, a_{2}, a_{3}) = \begin{cases}a_{1} & if |a_{1}| ≤ MΔ_{j}^{2}, \\ m (a_{1}, a_{2}, a_{3}) & otherwise, \end{cases}

where M is, of course, an upper bound of the absolute value of the second-order derivative of the solution at local extrema. In the case of the non-linear conservation laws under consideration, it is easy to prove that, if the initial data is piecewise C^{2}, we can take

M = C sup{ | (u_{0})_{xx}(y) |, y : (u_{0})_{x}(y) = 0},

where for a uniform mesh we could take C = 2/3; see [5]. Fortunately, in practice, the numerical results are not very sensitive to the choice of this constant which can be taken fairly big without degrading the quality of the results. It is remarkable that this small change in the limiter actually makes all the difference: the method now becomes fully high order in smooth regions in the same way as the version without limiters.

• We have used the total variation of the local means to devise our generalized slope limiter but stability in the L^{∞}-norm of the means is also enforced by this limiter; see [5] for details.

§3.4 Extension to bounded domains

In all this section, we have only considered periodic boundary conditions. To extend our results to the bounded domain case, see the formulation of the corresponding initial-boundary value problem by Bardos, LeRoux and Nédélec [2], we proceed in three steps:

(i) To extend the DG space discretization, we simply have to replace the numerical fluxes at the boundaries, namely,

f̂(u_{h}(0^{−}, t), u_{h}(0^{+}, t))  and  f̂(u_{h}(1^{−}, t), u_{h}(1^{+}, t)),

by

f̂(a(t), u_{h}(0^{+}, t))  and  f̂(u_{h}(1^{−}, t), b(t)),

respectively, where a(t) and b(t) are the boundary data.

(ii) The TDV-RK time discretization has now to take into consideration the boundary data; see the work by Shu [19] for details.

(iii) To extend the generalized slope limiter, we simply have to define the quantities

\overline{u}_{0} = 2 a(t) − \overline{u}_{1}  and  \overline{u}_{N+1} = 2 b(t) − \overline{u}_{N},

and proceed as usual.

§4 DG method for multi-dimensional hyperbolic systems

In this section, we consider the extension of the RKDG method to multi-dimensional systems:

u_{t} + ∇ · f(u) = 0,  in Ω × (0, T),  (11)
u(x, 0) = u_{0}(x),  ∀ x ∈ Ω.  (12)

For simplicity, we assume that Ω is the d-dimensional unit cube. The RKDG method for multidimensional systems has the same structure it has for one-dimensional scalar conservation laws; we only need to describe the DG-space discretization and the generalized slope limiter ΛΠ_{h}.

§4.1 The Discontinuous Galerkin space discretization

To discretize in space our multi-dimensional system (11), we simply proceed component by component; thus, it is enough to show how to do this in the case in which u is a scalar.

For this case, we seek an approximate solution u_{h} whose restriction to the element K of the triangulation T_{h} of Ω is, for each value of the time variable, in the local space U(K). Just as done in the one-dimensional case, we take u_{h}(0) ≡ u_{h}(·, 0) on the element K to be the L^{2}-projection of the data on U(K), that is, for all v_{h} ∈ U(K),

∫_{K} u_{h}(0) v_{h} dx = ∫_{K} u_{0} v_{h} dx.  (13)

We now determine the approximate solution for t > 0 on each element K of our triangulation by imposing that, for all v_{h} ∈ U(K),

∫_{K}(u_{h})_{t} v_{h} dx − ∫_{K} f(u_{h}) · ∇v_{h} dx + ∫_{∂K} { \atop f · n_{K}}(u_{h}) v_{h} ds = 0,  (14)

where n_{K} is the outward unit normal to the boundary of K.

To complete the definition of the DG space discretization, it only remains to define the numerical flux { \atop f · n_{K}}. This is in effect just a one dimensional flux we have discussed in the previous section, in the normal direction of the edge. However, to explain it clearly, we need to introduce some notation. For two adjacent elements K^{+} and K^{−} of the triangulation T_{h} and a point x of their common boundary at which the vectors n_{K^{±}} are well defined, we set

u_{h}^{±}(x) = lim_{ϵ↓0} u_{h}(x − ϵ n_{K^{±}}),

and call these values the traces of u_{h} from the interior of K^{±}. Now, just like for the onedimensional case, we take the numerical flux at x to be solely a function of the traces u_{h}^{±}(x), that is,

{ \atop f · n_{K^{−}}}(u_{h})(x) = { \atop f · n_{K^{−}}}(u_{h}^{−}(x), u_{h}^{+}(x)),

and require that it be consistent with the non-linearity f · n_{K^{−}}, which in this case amounts to ask that { \atop f · n_{K^{−}}}(a, a) = f(a) · n_{K^{−}}. Another criterion to pick our numerical fluxes is that when a piecewise-constant approximation is taken, the DG space discretization should give rise to a monotone finite volume scheme. This is ensured if we ask that our numerical flux be conservative, that is, that

{ \atop f · n_{K^{−}}}(u_{h}^{−}(x), u_{h}^{+}(x)) + { \atop f · n_{K^{+}}}(u_{h}^{+}(x), u_{h}^{−}(x)) = 0,

and that the mapping a ↦ { \atop f · n_{K^{−}}}(a, ·) be non-decreasing. The main examples of numerical fluxes satisfying all the above requirements are the following:

(i) The Godunov flux:

f̂ · n^{G}(a, b) = \begin{cases}min_{a≤u≤b} f · n (u),  if a ≤ b, \\ max_{b≤u≤a} f · n (u),  otherwise. \end{cases}

(ii) The Engquist-Osher flux:

f̂ · n^{EO}(a, b) = ∫_{0}^{b} min(f′ · n(s), 0) ds + ∫_{0}^{a} max(f′ · n(s), 0) ds + f · n(0);

(iii) The Lax-Friedrichs flux:

f̂ · n^{LF} (a, b) = \frac{1}{2} (f(a) + f(b)) · n − \frac{C}{2} (b − a),

where C = max_{inf u_{0}(x)≤s≤sup u_{0}(x)} |n · f′(s)|.

In other words, to define the multi-dimensional DG discretization, we can use simple onedimensional numerical fluxes.

Before discussing the DG discretization under consideration, we introduce some notations. Thus, we define the mean values · and jumps · by

u_{h} := \frac{1}{2}(u_{h}^{+} + u_{h}^{−}), u_{h} := u_{h}^{+}n_{K^{−}} + u_{h}^{−}n_{K^{+}},

we realize that we have the identity f̂ · n^{LF} = f̂^{LF} · n where

f̂^{LF} (u_{h}^{−}, u_{h}^{+}) = f(u_{h}) − \frac{C}{2} u_{h}.

The Godunov and the Engquist-Osher numerical fluxes do not satisfy a similar identity. Next, we discuss a few important points concerning this discretization:

• Just like in the one dimensional case, the mass matrix is block-diagonal; the block associated with the element K is a square matrix of order equal to the dimension of the local space U(K) and hence can be easily inverted. Moreover, for a variety of elements and spaces U(K), a basis can be found which is orthonormal in L^{2}. This is the case, for example, of rectangles and tensor product polynomials, in which case the orthonormal basis is a properly scaled tensor product of Legendre polynomials. Another remarkable example is that of simplexes and polynomials of a given total degree, case for which the so-called Dubiner basis is the orthonormal basis; see Dubiner [11].

Thus, after performing the DG space discretization, and just like for the one-dimensional case, the resulting equations can be rewritten in ODE form as \frac{d}{dt}u_{h} = L_{h}(u_{h}) where L_{h}(u_{h}) denotes the approximation to −∇ · f(u) provided by the DG method.

• In practice, the integrals appearing in the weak formulation (14) need to be approximated by quadrature rules. It was proven in [7] that

∥L_{h}(u) + ∇ · f(u)∥_{L^{∞}(K)} ≤ C h^{k+1}|f(u)|_{W^{k+2,∞}(K)},

if the quadrature rules over each of the faces of the border of the element K are exact for polynomials of degree 2k + 1, and if the one over the element is exact for polynomials of degree 2k. In fact, these requirements are also necessary, as we have verified numerically; moreover, the method is more sensitive to the quality of the quadrature rules used on the boundary of the elements than to that used in their interior.

Finally, let us point out that a quadrature-free version of the method was devised by Atkins and Shu [1] which results in a very efficient method for linear problems and certain nonlinear problems such as Euler equations of gas dynamics where the nonlinearity in the flux is mainly low order polynomials and perhaps one or two divisions of the components of the independent variable u.

• When dealing with multi-dimensional hyperbolic systems, the local Lax-Friedrichs numerical flux is a particularly convenient choice of numerical flux because it can be easily applied to any non-linear hyperbolic system, it is simple to compute, and because it gives good results.

This numerical flux is defined as follows. First, note that for multi-dimensional systems, u is a vector-valued function and f(u) is a matrix whose rows will be denoted by f_{j}(u); as a consequence, f̂^{LLF} is also a matrix whose j-th. row is given by

f̂_{j}^{LLF} (u_{h}^{−}, u_{h}^{+}) = f_{j}(u_{h}) − \frac{C}{2} (u_{h})_{j},

where C = C(K^{±}) is the larger one of the largest eigenvalue (in absolute value) of \frac{∂}{∂u}f(u_{h}^{±})·n_{K^{±}}. In practice, one could also determine C = C(K^{±}) to be the larger one of the largest eigenvalue (in absolute value) of \frac{∂}{∂u}f(\overline{u}_{K^{±}}) · n_{K^{±}} where \overline{u}_{K^{±}} are the means of the approximate solution u_{h} in the elements K^{±}.

• The DG space discretization can be applied to any high-order hyperbolic equation by simply rewriting it as a first order system of equations. For example, the wave equation

u_{tt} − c^{2}Δu = 0,

which is a second-order hyperbolic equation, can be rewritten as follows:

U_{t} + ∇ · F (U) = 0,

where,

U = (\begin{matrix}q_{1} \\ q_{2} \\ … \\ q_{d} \\ u \end{matrix}) , F (U) = − (\begin{matrix}u & 0 & … & 0 \\ 0 & u & … & 0 \\ … & … & … & … \\ 0 & 0 & … & u \\ c^{2} q_{1} & c^{2} q_{2} & … & c^{2} q_{d} \end{matrix}) .

The DG space discretization can now be easily applied to this system.

• Let us finally point out that since the wave equation can also be rewritten as

q_{t} − ∇u = 0,
u_{t} − c^{2}∇ · q = 0,

the DG space discretization of the hyperbolic system for U, can be rewritten in terms of (u, q) as follows: Find (u_{h}, q_{h}) such that its restriction to the element K belongs to the local space U(K) × U^{d}(K) and is such that, for all (v_{h}, r_{h}) ∈ U(K) × U^{d}(K),

∫_{K}(q_{h})_{t} · r_{h} dx + ∫_{K} u_{h} ∇ · r_{h} dx − ∫_{∂K} û_{h} r_{h} · n_{K} ds = 0,
∫_{K}(u_{h})_{t} v_{h} dx + ∫_{K} c^{2}q_{h} · ∇v_{h} dx − ∫_{∂K} ĉ^{2} q_{h} · n_{K} v_{h} ds = 0,

where the numerical fluxes û_{h} and ĉ^{2}q_{h} can be easily written in terms of the numerical flux F̂(U_{h}). For example, the Lax-Friedrichs flux for F̂(U_{h}) corresponds to

−û_{h} = −u_{h} − \frac{| c |}{2}q_{h}, −ĉ^{2} q_{h} = −c^{2}q_{h} − \frac{| c |}{2}u_{h},

where

u_{h} = \frac{1}{2}(u_{h}^{+} + u_{h}^{−}),  u_{h} = u_{h}^{+}n_{K^{−}} + u_{h}^{−}n_{K^{+}},
q_{h} = \frac{1}{2}(q_{h}^{+} + q_{h}^{−}),  q_{h} = q_{h}^{+} · n_{K^{−}} + q_{h}^{−} · n_{K^{+}}.

§4.2 The generalized slope limiter ΛΠ_{h}

The role of the generalized slope limiter ΛΠ_{h} in multi space dimensions is to enforce the maximum norm stability. In [7] Cockburn, Hou and Shu constructed a generalized slope limiter that enforces a local maximum principle without degrading the accuracy of the numerical scheme; this property holds for approximate solutions of arbitrary shapes and quite general meshes.

A practically simpler and effective limiter is advocated in [8]. There is however no rigorous proof that the use of this limiter does enforce the stability of the method, although numerically it seems to perform well. Let v_{h} be the function to which we are going to apply the limiter and let u_{h} be the result; let also v_{h}^{1} be its L^{2}-projection into the space of piecewise linear functions. Inspired by the construction of the one-dimensional limiter described in the previous section, we first construct a slope limiter for piecewise linear functions, ΛΠ_{h,M}^{1}. Then, we construct a limiter for general function as follows:

(i) Detect the spurious oscillations in v_{h}|_{K}, (ii) If there is no spurious oscillation, set u_{h}|_{K} = v_{h}|_{K}, (iii) If not, take u_{h}|_{K} equal to ΛΠ_{h,M}^{1}v_{h}^{1}. It remains now to decide how to ‘detect the spurious oscillations’. To do that, we assume that spurious oscillations are present in v_{h}|_{K} on the element K only if they are present in v_{h}^{1}|_{K} and by this we mean that v_{h}^{1}|_{K} ≠ ΛΠ_{h,M}v_{h}^{1}|_{K}. Thus, our generalized slope limiter is defined on the element K as follows:

(i) Compute r_{h}|_{K} = ΛΠ_{h,M}^{1}v_{h}^{1}|_{K}, (ii) If r_{h}|_{K} = v_{h}^{1}|_{K}, set u_{h}|_{K} = v_{h}|_{K}, (iii) If not, set u_{h}|_{K} = r_{h}|_{K}. It only remains to define the slope limiter ΛΠ_{h,M}^{1}. To construct it for triangular elements, we proceed as follows. Consider the triangles in Fig. 1, where m_{1} is the mid-point of the edge on the boundary of K_{0} and b_{i} denotes the barycenter of the triangle K_{i} for i = 0, 1, 2, 3. Since

m_{1} − b_{0} = α_{1} (b_{1} − b_{0}) + α_{2} (b_{2} − b_{0}),

for some nonnegative coefficients α_{1}, α_{2} which depend only on m_{1} and the geometry, we can write, for any linear function v_{h},

v_{h}(m_{1}) − v_{h}(b_{0}) = α_{1} (v_{h}(b_{1}) − v_{h}(b_{0})) + α_{2} (v_{h}(b_{2}) − v_{h}(b_{0})),

and since

\overline{v}_{K_{i}} = \frac{1}{|K_{i}|} ∫_{K_{i}} v_{h} = v_{h}(b_{i}), i = 0, 1, 2, 3,

we have that

ṽ_{h}(m_{1}, K_{0}) ≡ v_{h}(m_{1}) − \overline{v}_{K_{0}} = α_{1} (\overline{v}_{K_{1}} − \overline{v}_{K_{0}}) + α_{2} (\overline{v}_{K_{2}} − \overline{v}_{K_{0}}) ≡ Δ\overline{v}(m_{1}, K_{0}).

                           NO FILE: fig_19_chap4.eps

���� 1: Illustration of limiting.�

Now, we are ready to describe the slope limiting. Let us consider a piecewise linear function v_{h}, and let m_{i}, i = 1, 2, 3, be the three mid-points of the edges of the triangle K_{0}. We then can write, for (x, y) ∈ K_{0},

v_{h}(x, y) = ∑_{i=1}^{3} v_{h}(m_{i}) φ_{i}(x, y) = \overline{v}_{K_{0}} + ∑_{i=1}^{3} ṽ_{h}(m_{i}, K_{0}) φ_{i}(x, y).

To compute ΛΠ_{h}^{1}v_{h}, we first compute the quantities

Δ_{i} = \overline{m}( ṽ_{h}(m_{i}, K_{0}), ν Δ\overline{v}(m_{i}, K_{0}) ),

where \overline{m} is the TVB modified minmod function (without its third argument) and ν is an auxiliary parameter which we took equal to 2 in the one-dimensional case. Then, if ∑_{i=1}^{3} Δ_{i} = 0, we simply set

ΛΠ_{h,M}^{1}v_{h}(x, y) = \overline{v}_{K_{0}} + ∑_{i=1}^{3} Δ_{i} φ_{i}(x, y).

Note that if v_{h} is a linear function, then ṽ_{h}(m_{i}, K_{0}) = Δ\overline{v}(m_{i}, K_{0}) and Δ_{i} = ṽ_{h}(m_{i}, K_{0}) provided ν > 1; in this case we have ΛΠ_{h}^{1}(v_{h}) = v_{h}. This ensures that there is no degradation of accuracy after the application of the slope limiter away of critical points; when there are critical points, the suitable choice of the parameter M, hidden in the definition of the TVB modified minmod function, ensures the same effect.

If ∑_{i=1}^{3} Δ_{i} ≠ 0, we compute

pos = ∑_{i=1}^{3} max(0, Δ_{i}), neg = ∑_{i=1}^{3} max(0, −Δ_{i}),

and set

θ^{+} = min (1, \frac{neg}{pos} ) , θ^{−} = min (1, \frac{pos}{neg} ) .

Then, we define

ΛΠ_{h}v_{h}(x, y) = \overline{v}_{K_{0}} + ∑_{i=1}^{3} Δ̂_{i} φ_{i}(x, y),

where

Δ̂_{i} = θ^{+} max(0, Δ_{i}) − θ^{−} max(0, −Δ_{i}).

For systems, limiting in the local characteristic variables gives remarkably superior results than doing it component-by-component. Thus, to limit the vector ṽ_{h}(m_{i}, K_{0}) in the element K_{0}, we proceed as follows:

• Find the matrix R and its inverse R^{−1}, which diagonalizes the Jacobian

J = \frac{∂}{∂u}f(\overline{v}_{K_{0}}) · \frac{m_{i} − b_{0}}{|m_{i} − b_{0}|},

that is, R^{−1} J R = Λ, where Λ is a diagonal matrix containing the eigenvalues of J. Notice that the columns of R are the right eigenvectors of J and the rows of R^{−1} are the left eigenvectors.

• Transform ṽ_{h}(m_{i}, K_{0}) and Δ\overline{v}(m_{i}, K_{0}) to the characteristic fields. This is achieved by left multiplying these vectors by R^{−1}.

• Apply the scalar limiter to each of the components of the transformed vectors.

• The result is transformed back to the original space by multiplying R on the left.

���������

[1] H.L. Atkins and C.-W. Shu, Quadrature-free implementation of discontinuous Galerkin methods for hyperbolic equations, AIAA Journal 36 (1998), 775–782.

[2] C. Bardos, A.Y. LeRoux, and J.C. Nédélec, First order quasilinear equations with boundary conditions, Comm. in P.D.E. 4 (1979), 1017–1034.

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

[4] 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., to appear.

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

[6] B. Cockburn, S.-Y. 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 (1989), pp. 90–113.

[7] 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. Comp., 54 (1990), pp. 545–581.

[8] B. Cockburn and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws V: multidimensional systems, J. Comput. Phys., 141 (1998), pp. 199-224.

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

[10] B. Cockburn, G. Karniadakis and C.-W. Shu, The development of discontinuous Galerkin methods, in Discontinuous Galerkin Methods: Theory, Computation and Applications, B. Cockburn, G. Karniadakis and C.-W. Shu, editors, Lecture Notes in Computational Science and Engineering, volume 11, Springer, 2000, Part I: Overview, pp. 3-50.

[11] M. Dubiner, Spectral methods on triangles and other domains, J. Sci. Comp. 6 (1991), 345–390.

[12] G. Jiang and C.-W. Shu, On cell entropy inequality for discontinuous Galerkin methods, Math. Comp., 62 (1994), pp. 531–538.

[13] C. Johnson and J. Pitkäranta, An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation, Math. Comp. 46 (1986), 1–26.

[14] P. LeSaint and P.A. Raviart, On a finite element method for solving the neutron transport equation, Mathematical aspects of finite elements in partial differential equations (C. de Boor, ed.), Academic Press, 1974, pp. 89–145.

[15] S. Osher, Convergence of generalized muscl schemes, SIAM J. Numer. Anal. 22 (1985), 947–961.

[16] T. Peterson, A note on the convergence of the discontinuous Galerkin method for a scalar hyperbolic equation, SIAM J. Numer. Anal. 28 (1991), 133–140.

[17] W. H. Reed and T. R. Hill, Triangular mesh methods for the neutron transport equation, Tech. Report LA-UR-73-479, Los Alamos Scientific Laboratory, 1973.

[18] C.-W. Shu, TVB uniformly high order schemes for conservation laws, Math. Comp. 49 (1987), 105–121.

[19] C.-W. Shu, TVB boundary treatment for numerical solutions of conservation laws, Math. Comp. 49 (1987), 123–134.

[20] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes, J. Comput. Phys., 77 (1988), pp. 439–471.

[21] B. van Leer, Towards the ultimate conservation difference scheme, II, J. Comput. Phys. 14 (1974), 361–376.

[22] , Towards the ultimate conservation difference scheme, V, J. Comput. Phys. 32 (1979), 1–136.