AN ADAPTIVE DISCONTINUOUS GALERKIN TECHNIQUE WITH AN
ORTHOGONAL BASIS APPLIED TO COMPRESSIBLE FLOW PROBLEMS
JEAN-FRANÇOIS REMACLE
, JOSEPH E. FLAHERTY
, AND MARK S. SHEPHARD
Abstract. We present a high-order formulation for solving hyperbolic conservation laws using the Discon-
tinuous Galerkin Method (DGM). We introduce an orthogonal basis for the spatial discretization and use explicit
Runge-Kutta time discretization. Some results of higher-order adaptive refinement calculations are presented for in-
viscid Rayleigh Taylor flow instability and shock reflexion problems. The adaptive procedure uses an error indicator
that concentrates the computational effort near discontinuities.
Key words. Discontinuous Galerkin, adaptive meshing , Orthogonal Basis.
1. Introduction. The Discontinuous Galerkin Method (DGM) was initially introduced
by Reed and Hill in 1973 [16] as a technique to solve neutron transport problems. Lesaint [13]
presented the first numerical analysis of the method for a linear advection equation. However,
the technique lay dormant for several years and has only recently become popular as a method
for solving fluid dynamics or electromagnetic problems [4]. The DGM is somewhere between
a finite element and a finite volume method and has many good features of both.
Finite element methods (FEMs), for example, involve a double discretization. First, the
physical domain
is discretized into a collection of
elements
(1.1)
called a mesh. Then, the continuous function space
containing the solution of the
problem is approximated on each element
of the mesh, defining a finite-dimensional space
. The DGM is a finite element method in the sense that both geometrical and functional
discretizations define the finite-dimensional approximation space
.
The accuracy of a finite element discretization depends both on geometrical and func-
tional discretizations. Adaptivity seeks an optimal combination of these two ingredients:
p-refinement is the expression used for functional enrichment and h-refinement for mesh en-
richment.
Classical continuous FEMs typically use conforming meshes where elements share only
complete boundary segments. Thus, spatial discretizations like those shown in Figure 1.1
would, normally, not be allowed. Since the approximation space
is also constrained to
be a subspace of a continuous function space, e.g.,
, the basis (shape functions) for
are
typically associated with element vertices, edges, faces, or interiors. These simplify the impo-
sition of continuity requirements but limit choices. The DGM allows more general mesh con-
figurations and discontinuous bases (see Figure 1.1) which simplify both h- and p-refinement.
For example, non-conforming meshes and arbitrary bases for functional approximation [20]
may be used. In particular, we use a
-orthogonal basis that yields a diagonal mass matrix.
The DGM can also be regarded as an extension of Finite Volume Methods (FVMs) to
arbitrary orders of accuracy without the need to construct complex stencils for high-order
reconstruction. Indeed, the DGM stencil remains invariant for all polynomial degrees. This
greatly simplifies parallel implementation for methods of all polynomial orders. Finally, the
This work was supported ASCI Flash Center at the University of Chicago, under contract B341495, by the U.S.
Army Research Office through grant DAAG55-98-1-0200, and by the National Science Foundation through grant
DMS-0074174
Scientific Computation Research Center, Rensselaer Polytechnic Institute, Troy, New York, USA
1
!
"
$#
%
%'&
)(
%*&
+
%*&
,
FIG. 1.1. Non-conforming mesh with discontinuous approximations
-. /10
on elements
243
.
DGM has aspects in common with finite difference schemes in that it may use fluxes associ-
ated with Riemann problems [3].
Herein, we concentrate on DGM formulations for hyperbolic conservation laws (§2).
For the spatial discretization, we choose an orthogonal basis that diagonalizes the mass matrix
and, thus, simplifies its evaluation (§2.1). Time discretization is performed by an explicit total
variation diminishing Runge-Kutta scheme [3]. To improve the performance of the explicit
integration, we use a new local time stepping procedure similar to one used by Flaherty et
al. [7] and which will be explained in a forthcoming paper [18].
We present procedures to perform adaptive computations where the discretization space
changes in time. Because of the flexibility of the DGM, we are able to change both
mesh and elementary polynomial orders often, e.g., several thousand times with very little
computational overhead.
Transient computation of unstable flows provide an application where adaptivity in time
is crucial. The instability of an interface separating miscible fluids of different densities
subject to gravity is known as a Rayleigh-Taylor Instability (RTI). Bubbles (spikes) of lighter
(heavier) fluid penetrate into the heavier (lighter) fluid, leaving behind a region where the
two fluids are mixed. This mixing region quickly becomes irregular and may provide an
understanding of turbulence since the flow there has chaotic features [11, 22].
Young et al. [22] solved an incompressible RTI problem governed by the Boussinesq
equations using spectral methods. We likewise believe that the complex structure of the mix-
ing zone could be efficiently represented by high-order polynomials. Fryxell [9] used a piece-
wise parabolic method [21] with adaptive h-refinement to solve compressible RTI problems in
two and three dimensions. Without explicit interface tracking [10], h-adaptivity will certainly
be necessary to accurately represent the complex evolution of bubbles and spikes [11].
We present solutions of a standard two-dimenstional RTI problem using h- and p-refinement.
Increasing the polynomial degree
5
improves the quality of the solution. However, p-refinement
alone is not effective for capturing the fine scale structures near discontinuities. Using an er-
ror indicator based on solution jumps, we present results for the same problem using adaptive
h-refinement and compare computations with those using adaptive p-refinement. Finally, an
2
adaptive hp-refinement computation is performed which is shown to be the best of these RTI
calculations.
2. Discontinuous Finite Element Formulation for Conservation Laws. Consider an
open set
7698
#
whose boundary
:;
is Lipschitz continuous with a normal
<
=
that is defined
everywhere. We seek to determine
>?A@)B)DCE8
#GF
8IHJ
LKM7
as the solution of a
system of conservation laws
:EN>PORQSUT
<
V
U>*WYXEZ
(2.1)
Here
Q[S\T]]_^a`cbA@dZeZdZe@f^g`hbi
is the vector valued divergence operator and
<
V
U>*j]
<
k
U>*l@dZeZeZd@
<
k
K
_>[f
is the flux vector with the
m
th component
<
k[n
_>[oC[f
_f
K
HpU^g`hbq@r
. Function space
U^a`cbs@4
consists of square integrable vector valued functions whose divergence is also
square integrable i.e.,
U^a`cbs@4jutv<
wx&
<
wPy
#
,
^g`hbG<
wy
{z|Z
With the aim of constructing a Galerkin form of (2.1), let
)}h@d} 4~
and
L}c@e})~
respectively
denote the standard
and
U:;
scalar products. Multiply equation (2.1) by a test
function
y
_
, integrate over
and use the divergence theorem to obtain the following
variational formulation
_:
N
>?@))~9
<
V
U>*l@4XgQGL~xO7
<
V
_>['}<
=
@)f~\X@fL~W@
y
lZ
(2.2)
Finite element methods (FEMs) involve a double discretization. First, the physical domain
is discretized into a collection of
elements like in (1.1) The continuous function space
containing the solution of (2.2) is approximated on each element
of the mesh to
define a finite-dimensional space
. With discontinuous finite elements,
is a “broken”
function space that consists in the direct sum of elementary approximations
>
(we use here
a polynomial basis
*
of order
):
M$>
&
>
y
_
K
@4>
y
K
7
@
y
d
Z
(2.3)
Because all approximation are disconnected, we can solve the conservation laws on each
element to obtain
U:ENL>
@)
9
<
V
U>
l@4XgQ|
O7
V
@f
$
]UXE@f
@
y
lZ
(2.4)
Now, a discontinuous basis implies that the normal trace
V
<
V
U>*j};<
=
is not defined
on
:
. In this situation, a numerical flux
V{
U>
@4>
L
is usually used on each portion
:
)
of
:
shared by element
and neighboring element
. Here,
>
and
>
)
are the restrictions of
solution
>
, respectively, to element
and element
. This numerical flux must be continuous,
so
<
V
y
U^a`cbs@4)K
, and be consistent, so
Vi
_>@4>[?
<
V
U>*}E<
=
. With such a numerical flux,
equation (2.4) becomes
U:
N
>[@))D
<
V
_>[el@fXgQofjO
V
U>*$@4>[Ll@)f$Lo]UX@))i@
y
@
(2.5)
3
¡
¢¤£)¥
¦§g¨
£
§¦
¥
¨
FIG. 2.1. Reference triangular element
where
=
is the number of faces of element
. Only the normal traces have to be defined on
:
and several operators are possible [12, 21]. It is usual to define the trace as the solution
of a Riemann problem across
:
©
. Herein, when we consider problems with strong shocks
[6,21], an exact Riemann solver is used to compute the numerical fluxes and a slope limiter [2]
is used to produce monotonic solutions when polynomial degrees
«ª¬
. For Rayleigh-Taylor
instabilities, Roe’s flux linearization [19] is used with a physical limiter that we describe in
§4.
2.1. Spatial discretization. Even if DGM solutions do not depend on the choice of basis
(because they all span
*$
), some of them are more convenient and computationally effi-
cient than others. We construct an orthogonal basis of
j$
with respect to the
scalar
product. As a result, an explicit time integration scheme will neither necessitate “lumping”
nor inversion of the mass matrix. Another advantage, which is, perhaps, more important, is
that the orthogonal basis makes p-refinement trivial.
In two dimensions, consider a right-triangular reference element as shown in Figure 2.1.
Without a need to maintain inter-element continuity, consider a basis of
5
NU®
degree monomials
in
\¯a@f°g
, i.e.,
±
²³
n
@)m'²´ @eZeZdZe@4µ
´@f¯a@)°@)¯
@)¯°@)°
@)¯
#
@dZeZdZe@)°
@
(2.6)
with
µ
_OM´eU{O¶ f·©¶
. This basis is said to be hierarchical in the sense that, if
¸
n
is
the space of
m
NU®
-degree polynomials in
\¯a@)°g
, then
±
¹º»½¼U
'¾
ºd¿?»]¼_
¾
¿»¹}e}e}»]¼_
lÀ
'¾
¿'Z
(2.7)
Any field
Á
is approximated in
*
as
Á*U¯a@)°gÃÂ*Ä
n
Á
n
³
n
\¯a@f°aZ
(2.8)
Let us define the scalar product on the reference element as
³
n
@4³fÅdLÆ{MÇ
º
Ç
ÀÈ
º
³
n
³fÅ?É ¯©É °
(2.9)
and the induced norm
Ê
Á
Ê
Æ
½_Á@rÁ;LÆZ
(2.10)
4
ËiÌÍiÎ'ÏÏ
1.414213562373095E+00
Ë
ÌÍiÎÐfÏ
-2.00000000000000E+00
ËiÌÍiÎÐÐ
6.000000000000000E+00
ËiÌÍiÎÑ
Ï
-3.464101615137754E+00
Ë
ÌÍiÎÑ
Ð
3.464101615137750E+00
ËiÌÍiÎ
ÑÑ
6.928203230275512E+00
ËiÌÍiÎÒ4Ï
2.449489742783153E+00
Ë
ÌÍiÎ
ÒLÐ
-1.959591794226528E+01
Ë
ÌÍiÎÒ
Ñ
0.000000000000000E+00
ËiÌÍiÎÒÒ
2.449489742783160E+01
ËiÌÍiÎÓfÏ
4.242640687119131E+00
Ë
ÌÍiÎ
ÓÐ
-2.545584412271482E+01
Ë
ÌÍiÎÓ
Ñ
-8.485281374238392E+00
ËiÌÍiÎÓ1Ò
2.121320343559552E+01
ËiÌÍiÎÓÓ
4.242640687119219E+01
Ë
ÌÍiÎÔ
Ï
5.477225575051629E+00
ËiÌÍiÎ
Ô
Ð
-1.095445115010309E+01
Ë
ÌÍiÎ
ÔÑ
-3.286335345030997E+01
ËiÌÍiÎ
Ô
Ò
5.477225575051381E+00
Ë
ÌÍiÎ
Ô
Ó
3.286335345031001E+01
Ë
ÌÍiÎÔÔ
3.286335345030994E+01
TAB L E 2.1
Second-order basis coefficients on the reference triangle of Figure 2.1.
We seek an alternate basis
ÕÖpd×
n
@fmvÖ´@dZeZeZe@fµ
of
[
which is orthonormal, i.e.,
×
n
@L×
Å
ÙØ
n
Å
. For this purpose, we apply Gram-Schmidt orthogonalization to basis
±
and construct
×
n
n
Åf[ÚqÛ[ÜqÝ
n
Å
³fÅ
(2.11)
with
Ú
Û[ÜoÝ
Þd×E@eZdZeZd@L×
Âß
a triangular matrix representing the change of basis from
±
to
Õ
. The
à
NU®
column of
Ú
Û*ÜqÝ
is the coordinates of
×Å
in basis
±
. Any shape function
³
may be expressed as
³
I¯ágâ
lã
° äâ
lã
with exponents
å?æg
and
ç_æ
depending on
æ
. Scalar
products
_³
n
@r³fÅ)
are calculated as
³
n
@r³fÅd)MÇ
º
Ç
ÀÈ
º
¯
K
°
ɯ©É °è
´
=
OI´
Ç
º
¯
K
L´W¯
©é
É ¯
´
=
O¹´
©é
ê
ºë
ê
Ç
º
¯
K
é
ê
É ¯q
´
=
OI´
©é
ê
º
ë
ê
ì
O¤íO¹´
(2.12)
with
ì
½åD\mL'O9åDhà
and
=
½ç\mL'Oçhà
. This simple result (2.12) avoids the need for
numerical integration in the Gram-Schmidt process so that any order shape functions can be
computed without a loss of precision. In Table 2.1, we give the transformation
Ú
Û*ÜqÝ
for a
complete second-order basis (q=2).
Integration of shape functions are usually not done in the parametric coordinates
\¯a@f°a
of the element but in the actual coordinates
!
@
"
.
orthogonality of
×
n
’s will only
be preserved if the mapping from the actual to the parametric coordinates is linear, i.e., the
Jacobian of the mapping is constant. Curved elements, which are essential for higher-order
5