A Comparison of Hybridized and
Standard DG Methods for
Target-Based hp-Adaptive
Simulation of Compressible Flow
Michael Woopen, Aravind Balan, Georg May, and Jochen Schütz

Aachen Institute
for Advanced Study in
Computational Engineering Science
Financial support from the
Deutsche Forschungsgemeinschaft (German Research Foundation)
through grant GSC 111 is gratefully acknowledged.

A COMPARISON OF HYBRIDIZED AND STANDARD DG
METHODS FOR TARGET-BASED hp-ADAPTIVE SIMULATION
OF COMPRESSIBLE FLOW

MICHAEL WOOPEN, ARAVIND BALAN, GEORG MAY
Aachen Institute for Advanced Study in Computational Engineering Science,
RWTH Aachen University, Aachen, Germany

JOCHEN SCHÜTZ
Institut für Geometrie und Praktische Mathematik, RWTH Aachen University,
Aachen, Germany

Abstract. We present a comparison between hybridized and non-hybridized
discontinuous Galerkin methods in the context of target-based hp-adaptation
for compressible flow problems. The aim is to provide a critical assessment of
the computational efficiency of hybridized DG methods.
Hybridization of finite element discretizations has the main advantage, that
the resulting set of algebraic equations has globally coupled degrees of freedom only on the skeleton of the computational mesh. Consequently, solving for
these degrees of freedom involves the solution of a potentially much smaller system. This not only reduces storage requirements, but also allows for a faster
solution with iterative solvers. Using a discrete-adjoint approach, sensitivities with respect to output functionals are computed to drive the adaptation.
From the error distribution given by the adjoint-based error estimator, h- or
p-refinement is chosen based on the smoothness of the solution which can be
quantified by properly-chosen smoothness indicators.
Numerical results are shown for subsonic, transonic, and supersonic flow
around the NACA0012 airfoil. hp-adaptation proves to be superior to pure
h-adaptation if discontinuous or singular flow features are involved. In all
cases, a higher polynomial degree turns out to be beneficial. We show that for
polynomial degree of approximation p = 2 and higher, and for a broad range
of test cases, HDG performs better than DG in terms of runtime and memory
requirements.

1. Introduction
During the last years, discontinuous Galerkin (DG) methods (see, e.g., [1, 2, 3])
have become increasingly popular. This is indisputable due to their advantages
— high-order accuracy on unstructured meshes, a variational setting, and local
conservation, just to name a few.
However, the use of discontinuous function spaces is at the same time the reason
for a major disadvantage: unlike in continuous Galerkin (CG) methods, degrees
1

2

WOOPEN, BALAN, MAY, SCHÜTZ

of freedom are not shared between elements. As a consequence, the number of
unknowns is substantially higher compared to a CG discretization. Especially for
implicit time discretization this imposes large memory requirements, and potentially leads to increased time-to-solution.
In order to avoid these disadvantages, a technique called hybridization may be
utilized (see [4, 5, 6, 7, 8, 9, 10]), resulting in hybridized discontinuous Galerkin
(HDG) methods. Here, the globally coupled unknowns have support on the mesh
skeleton, i.e. the element interfaces, only. This reduces the size of the global system
and coincidentally improves the sparsity pattern.
However, aiming at industry applications, e.g. turbulent flow around a complete airplane or within an aircraft engine, hybridization alone does most likely
not provide a sufficiently successful overall algorithm. In these applications one
is usually interested in certain quantities only, for example lift or drag coefficients
in aerospace, instead of the solution quality per se. Thus, it might be beneficial
to distribute the degrees of freedom within the computational domain in such a
way that the solution to the discretized problem is close to optimal with respect
to the accuracy of these quantities. To achieve this goal, target-based error control
methods have been developed (see [11, 12, 13, 14, 15]). One such method is based
on the adjoint solution of the original governing equations with respect to the target functional. In this method, an additional linear system of partial differential
equations is solved which then gives an estimate on the spatial error distribution
contributing to the error in the target functional. This estimate can be used as
a criterion for local adaptation. Within the context of low order schemes, mesh
refinement is used for adaptation [12, 13]. Using DG (or HDG), however, offers the
additional possibility of varying the polynomial degree within each element. For
smooth solutions, this is more efficient compared to mesh refinement, as it yields
exponential convergence. In the context of wave problems, Giorgiani et al. [16]
showed the benefit of using p-adaptation within an HDG-framework. Combining
both mesh- and order-refinement results in so-called hp-adaptation.
In [17], we presented a discretization method for nonlinear convection-diffusion
equations. The method is based on a discontinuous Galerkin discretization for
convection terms, and a mixed method using H(div) spaces for the diffusive terms.
Furthermore, hybridization is used to reduce the number of globally coupled degrees
of freedom. Adjoint consistency was shown in [18]. In [19, 20], we extended our
computational framework to include HDG schemes, as well as adjoint-based hand hp-adaptation. In the current paper, we compare our HDG method with a
standard DG method in the context of hp-adaptation for stationary compressible
flow, mainly with the aim to assess the efficiency of both methods.
This paper is structured as follows. We briefly cover the governing equations,
namely the compressible Euler and Navier-Stokes equations, in Sec. 2. After that
we introduce our discretization and describe the concept of hybridization in Sec. 3.
In Sec. 4 we establish the adjoint formulation and show how hybridization can be
applied to the dual problem. Then we show its efficiency and robustness with
examples from compressible flow, including the subsonic, transonic, and supersonic
regime, in Sec. 5. Finally, we offer conclusions and outlook on future work in Sec. 6.

A COMPARISON OF HDG AND DG FOR TARGET-BASED hp-ADAPTATION

3

2. Governing Equations
We consider systems of partial differential equations
(1)

∇ · (fc (w) − fv (w, ∇w)) = s (w, ∇w)

with convective and diffusive fluxes
(2)

fc : Rm → Rm×d

and fv : Rm × Rm×d → Rm×d ,

respectively, and a state-dependent source term
(3)

s : Rm × Rm×d → Rm

on domain Ω ⊂ Rd . Potentially, some of these quantities could be zero. We
denote the spatial dimension by d and the number of conservative variables by m.
Boundary conditions can be applied either to the conservative variables w ∈ Rm
and their gradient ∇w ∈ Rm×d or directly to the fluxes fc and fv .
2.1. Two-Dimensional Euler Equations. The Euler equations are comprised
of the inviscid compressible continuity, momentum and energy equations. They are
given in conservative form as
(4)

∇ · fc (w) = 0

with the vector of conserved variables
(5)

w = (ρ, ρv, E)T

where ρ is the density, v is the velocity vector v := (vx , vy )T , and E the total
energy. The convective flux is given by
(6)

T

fc = (ρv, p Id + v ⊗ v, v(E + p)) .

Pressure is related to the conservative flow variables w by the equation of state


1
(7)
p = (γ − 1) E − ρv · v
2
where γ = cp /cv is the ratio of specific heats, generally taken as 1.4 for air.
Along wall boundaries we apply the slip boundary condition
(8)

vn (w) := v · n = 0.

We also define a boundary function which satisfies vn (w∂Ω (w)) = 0 as


1
0T
0
(9)
w∂Ω (w) = 0 Id − n ⊗ n 0 w.
0
0T
1
Prescribing boundary conditions at the far-field can be realized with the aid of
characteristic upwinding [21]. Here, the normal convective flux Jacobian is decomposed as
(10)

fc0 (w) · n = Q(w, n) · Λ(w, n) · Q−1 (w, n)

with Λ(w, n) being a diagonal matrix containing the eigenvalues of fc0 (w) · n. The
corresponding right eigenvectors can be found in the columns of Q(w, n). The interior and far-field states in characteristic variables are then given by wc = Q(w, n)w

4

WOOPEN, BALAN, MAY, SCHÜTZ

and wc,∞ = Q(w, n)w∞ , respectively. Finally, depending on the sign of (Λ(w, n))i,i ,
we can construct a boundary state
(
(Q(w, n)wc )i ,
(Λ(w, n))i,i ≥ 0
(11)
(w∂Ω (w))i =
(Q(w, n)wc,∞ )i , (Λ(w, n))i,i < 0.
2.2. Two-Dimensional Navier-Stokes Equations. The Navier-Stokes equations
in conservative form are given by
∇ · (fc (w) − fv (w, ∇w)) = 0.

(12)

The convective part fc of the Navier-Stokes equations coincides with the Euler
equations. The viscous flux is given by
T

(13)

fv = (0, τ , τ v + k∇T ) .

The temperature is defined via the ideal gas law


µγ
E
1
1
p
(14)
T =
− v·v =
k · Pr ρ
2
(γ − 1)cv ρ
µc

where Pr = kp is the Prandtl number, which for air at moderate conditions can be
taken as a constant with a value of Pr ≈ 0.72. k denotes the thermal conductivity
coefficient. For a Newtonian fluid, the stress tensor is defined as


2
T
(15)
τ = µ ∇v + (∇v) − (∇ · v) Id .
3
The variation of the molecular viscosity µ as a function of temperature is determined by Sutherland’s law as
(16)

µ=

C1 T 3/2
T + C2

√
with C1 = 1.458 × 10−6 mskg
and C2 = 110.4 K.
K
Along wall boundaries, we apply the no-slip boundary condition, i.e.

(17)

v=0

with corresponding boundary function
(18)

T

w∂Ω (w) = (ρ, 0, E) .

Furthermore, one has to give boundary conditions for the temperature. In the
present work we use the adiabatic wall condition, i.e.
(19)

∇T · n = 0.

Combining both no-slip and adiabatic wall boundary conditions, yields a condition
for the viscous flux, namely
T
(20)
fv (w∂Ω , q∂Ω ) = 0 τ 0 .

A COMPARISON OF HDG AND DG FOR TARGET-BASED hp-ADAPTATION

5

3. Discretization
3.1. Notation. We tesselate the domain
S Ω into a collection of non-overlapping
elements, denoted by Th , such that K∈Th K = Ω. For the element edges we
consider two different kinds of sets, ∂Th and Γh , which are element-oriented and
edge-oriented, respectively.
(21)
(22)

∂Th := { ∂K\∂Ω : K ∈ Th },
Γh := { e : e = K ∩ K 0 for K, K 0 ∈ Th ; measd−1 (e) 6= 0 }.

The first is the collection of all element boundaries, which means that every edge
appears twice. The latter, however, includes every edge just once. The reason for
this distinction will become clear later. Please note that neither of these sets shall
include edges lying on the domain boundary; the set of boundary edges is denoted
by Γbh .
We denote by Πp (D) the set of polynomials of degree at most p on some domain D. We will need discontinuous function spaces for the domain and the mesh
skeleton:
(23)
(24)
(25)

Vh = {v ∈ L2 (Ω) : v|K ∈ ΠpK (K),
2

Wh = {w ∈ L (Ω) : w|K ∈ Π

pK

(K),

Mh = {µ ∈ L2 (Γh ) : µ|e ∈ Πpe (e),

K ∈ Th }m×d
K ∈ Th }m
e ∈ Γh }m .

Thus, v ∈ Vh , w ∈ Wh and µ ∈ Mh are piecewise polynomials of degree p which
can be discontinuous across edges (for v, w) or vertices (for µ), respectively.
Usually, the polynomial degree between elements and interfaces does not vary.
In the case of varying polynomial degrees pK − and pK + , we follow Cockburn
et al. [6] and choose the polynomial degree for the interface e = K − ∩ K + as
pe = max {pK − , pK + }. This approach is commonly applied for p-adaptive HDG
discretizations (see Chen and Cockburn [22], Cockburn and Gopalakrishnan [23]
and Giorgiani et al. [16]).
We will distinguish between element-oriented inner products (defined with Th ),
i.e.
X Z
(v, w)Th :=
(26)
vw dx,
K∈Th

(27)

(v, w)Th :=

K∈Th

(28)

hv, wi∂Th :=

K

X Z

v · w dx,

K

X Z
K∈Th

vw dσ,

∂K

and edge-oriented inner products (defined with Γh ), i.e.
XZ
(29)
hv, wiΓh :=
vw dσ.
e∈Γh

e

3.2. Weak Formulation. We can rewrite general convection-diffusion equations
as a first-order system by introducing an additional unknown representing the gradient of the solution
q = ∇w
(30)
∇ · (fc (w) − fv (w, q)) = s (w, q) .

6

WOOPEN, BALAN, MAY, SCHÜTZ

By multiplying the strong, mixed form (30) with appropriate test functions
(τ h , ϕh ) ∈ Vh × Wh and integrating by parts, we obtain a standard DG discretization in mixed formulation of the problem, i.e.:
Find (qh , wh ) ∈ Vh × Wh s.t. ∀(τ h , ϕh ) ∈ Vh × Wh
(31)
(32)
(33)
(34)

0 = NhDG (qh , wh ; τ h , ϕh )
:= (τ h , qh )Th + (∇ · τ h , wh )Th − hτ h · n, wi
b ∂Th
− (∇ϕh , fc (wh ) − fv (wh , qh ))Th
D
E
+ ϕh , fbc − fbv
− (ϕh , s(wh , qh ))Th
∂Th

(35)

+

DG
Nh,∂Ω

DG
(qh , wh ; τ h , ϕh ) + Nh,sc
(qh , wh ; ϕh ) .

Here the numerical trace w
b and the numerical fluxes fbc , fbv have to be chosen
appropriately to define a stable and consistent method. Furthermore, boundary
conditions and shock-capturing are incorporated using appropriate operators, here
DG
DG
denoted by Nh,∂Ω
(qh , wh ; τ h , ϕh ) and Nh,sc
(qh , wh ; ϕh ).
In contrast to a DG discretization, where the numerical trace w
b is defined explicitly in terms of wh and qh , it is treated as an additional unknown in an HDG
method. This additional unknown is called λh and has support on the skeleton of
the mesh only. In order to close the system the continuity of the numerical fluxes
across edges is required in a weak sense, resulting in an additional equation.
The weak formulation of the hybrid system, comprised of equations for the gradient qh , the solution wh itself and its trace on the mesh skeleton λh , is then given
by:
Find (qh , wh , λh ) ∈ Xh := Vh × Wh × Mh s.t. ∀(τ h , ϕh , µh ) ∈ Xh
(36)
(37)
(38)
(39)

0 = Nh (qh , wh , λh ; τ h , ϕh , µh )
:= (τ h , qh )Th + (∇ · τ h , wh )Th − hτ h · n, λh i∂Th
− (∇ϕh , fc (wh ) − fv (wh , qh ))Th
D
E
+ ϕh , fbc − fbv
− (ϕh , s(wh , qh ))Th
∂Th

(40)
(41)

+ Nh,∂Ω (qh , wh ; τ h , ϕh ) + Nh,sc (qh , wh ; ϕh )
D
r
zE
+ µh , fbc − fbv
.
Γh

Please note the use of ∂Th in the weak formulation of the mixed form (see
Eq. (37)–(40)) and Γh in Eq. (41) defining λh . This perfectly resembles the character of these equations, being element- and edge-oriented, respectively. The terms
tested against τ h and ϕh are called local solvers, meaning they do not depend
on the solution within neighboring elements but only on the trace of the solution
which is approximated by λh . The coupling between elements is then introduced
by weakly enforcing the normal continuity of the numerical fluxes across interfaces
(see Eq. (41)).
We choose numerical fluxes similar to the Lax-Friedrich flux and to the LDG
flux for the convective and diffusive flux, respectively, i.e.
(42)

fbc (λh , wh ) = fc (λh ) · n − αc (λh − wh )

(43)

fbv (λh , wh , qh ) = fv (λh , qh ) · n + αv (λh − wh )

A COMPARISON OF HDG AND DG FOR TARGET-BASED hp-ADAPTATION

7

which can be combined into
(44)

fbc − fbv = (fc (λh ) − fv (λh , qh )) · n − α (λh − wh )

where α = αc + αv . The stabilization introduced can be given by a tensor; in
our work, however, we restrict ourselves to a constant scalar α which seems to be
sufficient for a wide range of test cases.
3.2.1. Boundary Conditions. In order to retrieve an adjoint-consistent scheme, special care has to be taken when discretizing the boundary conditions (see Schütz and
May [18]). The boundary conditions have to be incorporated by using the boundary
states w∂Ω (wh ) and gradients q∂Ω (wh , qh ), i.e.
(45)

Nh,∂Ω (qh , wh ; τ h , ϕh )
:= hτ h · n, w∂Ω iΓb + hϕh , (fc (w∂Ω ) − fv (w∂Ω , q∂Ω )) · niΓb .
h

h

We would like to emphasize that λh does not occur in this boundary term.
3.2.2. Shock-Capturing. In non-smooth parts of the solution, for example shocks in
compressible flows, a stabilization term has to be introduced. We adopt the shockcapturing approach by Hartmann and Houston [11] where an artificial viscosity
term, given by ∇ · ( (w, ∇w) ∇w), is used. The viscosity  is given by the L1 norm of the strong residual ∇ · fc (w) in every element. In order to accelerate the
convergence of this term to zero with mesh refinement, it is premultiplied with an
K
. The latter resembles the actual resolution within an
effective mesh size h̃K := hpK
element. Furthermore, a user-defined factor 0 is introduced which can be reliably
tuned for a rather large range of test cases. Finally, the artificial viscosity is given
by
Z
0 h̃2−β
K
(46)
|K :=
d(w) dx
|K|
K
where the strong residual is given by
(47)

d(w) :=

m
X

|(∇ · fc (w))i | .

i=1

In the discretization of this shock-capturing term the interface integral is neglected so that only the volume contribution is considered, i.e.
(48)

Nh,sc (wh ; ϕh ) := (∇ϕh ,  (wh , ∇wh ) ∇wh )Th .

This obviates the need for introducing qh in purely convective problems (e.g. the
compressible Euler equations). In the viscous case, where the gradient is explicitly
given, ∇wh can be replaced by qh yielding
(49)

Nh,sc (qh , wh ; ϕh ) := (∇ϕh ,  (wh , qh ) qh )Th .

Please note, that this term enters only the local part of the discretization.
Furthermore, it is noteworthy that this shock-capturing term is not only consistent and conservative but also yields an asymptotically adjoint-consistent scheme.

8

WOOPEN, BALAN, MAY, SCHÜTZ

3.3. Relaxation. In order to solve the nonlinear system of equations that defines
the HDG method, the Newton-Raphson
method is applied. Beginning with an

initial guess x0h := q0h , wh0 , λ0h , one iteratively solves the linear system
Nh0 [xnh ] (δ xnh ; yh ) = −Nh (xnh ; yh )

(50)

∀yh ∈ Xh

and updates the solution as

xn+1
= xnh + δ xnh
h

(51)

until the residual Nh (xnh ; yh ) has reached a certain threshold. Please note, that
we have grouped qh , wh and λh , and the test functions into xh := (qh , wh , λh )
and yh := (τ h , ϕh , µh ), respectively. Nh0 denotes the Fréchet derivative of Nh with
respect to xh .
This routine can, however, lead to stability problems if the starting value x0h is
too far away from the solution xh . Therefore, an artificial time is introduced and a
backward Euler method is applied, which yields a slight modification of the linear
system given in Eq. (50), namely ∀yh ∈ Xh


(52)

ϕh ,

1
δwhn
∆tn



n
n
+ Nh0 [xn
h ] (δ xh ; yh ) = −Nh (xh ; yh ) .
Th

Please note that by choosing ∆tn → ∞, a pure Newton-Raphson method is obtained. Usually the time step is kept finite for a few initial steps to ensure stability.
As soon as the residual is lower than a certain threshold, i.e. the current approximation xnh is thought to be sufficiently close to the solution xh , we let the time
step go towards infinity.
As we are interested in steady-state problems, we can use local time-stepping in
order to accelerate the computation. For each element, we apply a time step based
on a global CFL number, the element volume, and an approximation of the flux
Jacobian’s spectral radius, i.e.
(53)

∆tnK = CFLn

|K|
.
λc + 4λv

Here, λc and λv represent approximations to the maximum eigenvalue of the convective and diffusive flux, respectively (see Mavriplis and Jameson [24]).
The next question is, how to choose CFLn in each iteration such that for xnh →
xh , CFLn → ∞. Therefore, we choose a sequence (CFLn )n to consist of a ramping
phase followed by a very fast increase, i.e.
(54)

   
 3 
2


− 2 nn
c0 3 nn0
0


 n−1 
CFLn :=
kNh
k2

n−1

1 + c1 max 0, log
CFL
kN n k
h

2

n ≤ n0
otherwise

Strictly speaking, the parameters c0 , c1 and n0 depend on the actual problem;
experiments, however, showed that they can be reliably tuned for a wide range of
problems.
3.4. Hybridization. Using an appropriate polynomial expansion for δqh , δwh and
δλh , the linearized global system Eq. (50) is given in matrix form as


 

A B R
δQ
F
 C D S   δW  =  G 
(55)
L M N
δΛ
H

A COMPARISON OF HDG AND DG FOR TARGET-BASED hp-ADAPTATION

9

T

where the vector [δQ, δW, δΛ] contains the expansion coefficients of δxh with respect to the chosen basis.
In order to carry on with the derivation of the hybridized method, we want to
formulate that system in terms of δΛ only. Therefore we split it into


 
 

A B
δQ
F
R
(56)
=
−
δΛ
C D
δW
G
S
and
(57)



L

M





δQ
δW


+ N δΛ = H.

Substituting Eq. (56) into Eq. (57) yields the hybridized system
(58)

K δΛ = E.

where the system matrix and right-hand side vector are given by

 
!

 A B −1 R
(59)
K= N− L M
C D
S

−1 


 A B
F
E=H− L M
(60)
.
C D
G
The workflow is as follows: First, the hybridized system is assembled and then
being solved for δΛ. Subsequently, δQ and δW can be reconstructed inside the
elements via Eq. (56). It is very important to note that it is not necessary to solve
the large system given by Eq. (56). In fact, the matrix in Eq. (56) can be reordered
to be block diagonal. Each of these blocks is associated to one element. Thus, both
the assembly of the hybridized matrix in Eq. (58) and the reconstruction of δQ and
δW can be done in an element-wise fashion. In order to save computational time,
the solutions to Eq. (56) can be saved after the assembly of the hybridized system
and reused during the reconstruction of δQ and δW .
The hybridized matrix is an nf × nf block matrix, where nf = |Γh | is the
number of interior edges. In each block row there is one block on the diagonal
and 2d off-diagonal blocks in the case of simplex elements. These blocks represent
the edges of the
 neighboring elements of one edge. Each block is dense and has
O m2 · p2(d−1) entries. Please recall that p is the polynomial degree of the ansatz
functions, d is the spatial dimension of the domain Ω and m is the number of
partial differential equations (m = 4 for the 2-dimensional Euler or Navier-Stokes
equations). This structure is very similar to that
 of a normal DG discretization,
whereas the blocks in the latter have O m2 · p2d entries and thus are considerably
bigger for higher polynomial order p. The size of the system matrix does not only
play a big role in terms of memory consumption but also for the iterative solver.
Here, a major portion of the overall workload goes into matrix-vector products
which can obviously be performed faster, if the problem dimensions are smaller. In
our code we use an ILU(n)-preconditioned GMRES which is available through the
PETSc library [25, 26].
4. Adaptation Procedure
In the context of adjoint-based (also referred to as target- or output-based) error
estimation, one is interested in quantifying the error of a specific target functional

10

WOOPEN, BALAN, MAY, SCHÜTZ

Jh : Xh → R, i.e.
eh := Jh (x) − Jh (xh ) ,

(61)

where xh is the approximation to x in Xh . This target functional can for example
represent lift or drag coefficients in aerospace applications. In general, the target
functional is an integrated value, where integration can be both on a volume or
along the boundary. For the derivation of the adjoint-based error estimate we
expand the target functional in a Taylor series as follows

(62)
Jh (x) − Jh (xh ) = Jh0 [xh ] δ xh + O kδ xh k2
where we defined δ xh := x − xh .
We proceed in a similar manner with the error in the residual, i.e.

Nh (x; yh ) − Nh (xh ; yh ) = Nh0 [xh ] (δ xh ; yh ) + O kδ xh k2 .

(63)

As our discretization is consistent the first term Nh (x; yh ) vanishes.
Substituting Eq. (63) into Eq. (62) and neglecting the quadratic terms yields
eh ≈ η := −Nh (xh ; zh )

(64)
where

zh

is defined by the so-called adjoint equation

e h.
Nh0 [xh ] (yh ; zh ) = Jh0 [xh ] (yh )
∀yh ∈ X


eh ∈ X
e h represents the link between variaThe adjoint solution zh = qeh , w
eh , λ
tions in the residual and in the target functional.
The global error estimate η can then be restricted to a single element to yield a
local indicator to drive an adaptation procedure, i.e.
(65)

(66)

eh , w
eh , 0)
ηK := Nh (qh , wh , λh ; q

K

.

Here, we want to emphasize that, in contrast to the global error estimate, we ignore
eh . We found that by taking the
the contribution of the hybrid adjoint variable λ
whole adjoint into account, jumps across element interfaces were overly penalized.
This deserves a more in-depth analysis.
Please note, that the functionals Nh and Jh and their jacobians have to be
e h ⊃ Xh . Otherwise the
evaluated in a somewhat richer space than Xh , namely X
weighted residual Nh (xh ; zh ) would be identical zero as
(67)

Nh (xh ; yh ) = 0

∀yh ∈ Xh .

This can be achieved by either mesh refinement or a higher polynomial degree of
the ansatz functions. In our setting, especially when using a hierarchical basis, the
latter is advantageous with respect to implementational effort and efficiency.
4.1. Hybridization. In matrix form, the
follows

T 
A B R
 C D S  
(68)

L M N

adjoint system (see Eq. (65)) reads as
 
e
Q

f 
W
=
e
Λ


Fe
e 
G

e
H

e = 0 as λh is not defined on the boundary
Please note, that in our formulation H
and thus the target functional depends only on wh and qh .

A COMPARISON OF HDG AND DG FOR TARGET-BASED hp-ADAPTATION

11

As the overall structure of the adjoint equation is similar to the primal system
(see Eq. (58)), one can also apply static condensation to the adjoint system which
then yields its hybridized form:
e=E
e
KT Λ

(69)

where the right-hand side is given by
(70)

e=−
E



RT

ST





A
C

B
D

−T "

Fe
e
G

#
.

It is interesting to note that the hybridized adjoint system matrix is also the
transpose of the hybridized primal system matrix (for a higher polynomial order).
This is very beneficial for the implementation as the routines for the assembly of
this matrix are already available.
The adjoint solution within each element can then be computed with the aid of
the adjoint local system, given by
#
# "

T "
e

T
Q
A B
Fe
e
− L M
Λ
(71)
=
e
f
C D
G
W
where the matrix is also the transpose of the primal local matrix (see Eq. (56)).
4.2. Marking Elements for Refinement. After having obtained a localized error estimate, we have to choose a set of elements to be refined. This can be done
in many different ways. We adopt a strategy proposed by Dörfler [27]. The aim
of this marking strategy is to find the smallest set M ⊆ Th such that the error
contributed by this set represents a certain fraction of the total error, i.e.
(72)

ηM ≥ (1 − θ) ηTh .

The user-defined parameter θ is of course problem dependent. It can, however, be
tuned for a large rangePof test cases. Please note, that we define the error of any
2
2
subset of Th as ηM
:= K∈M ηK
.
4.3. Choosing between h- and p-Adaptation. The final step in the adaptation
procedure is the decision between mesh refinement and order enrichment. There
exist several ways to make this decision. Ceze and Fidkowski [28] solve local adjoint
problems for both options and then decide which one is more efficient with respect
to degrees of freedom or the non-zero entries in the system matrix. We, however,
adopt the strategy by Wang and Mavriplis [29]. They used a smoothness sensor
originally devised by Persson and Peraire [30] for an artificial viscosity approach.
This sensor exploits the fact that the decay of the expansion coefficients of the
solution is closely linked to its smoothness. For smoother solutions, decay is faster.
This can be exploited to check the regularity of the solution. On each element, the
smoothness sensor is defined as
(w − w? , w − w? )K
(73)
SK :=
(w, w)K
where w? is the element-wise projection of w to the next smaller polynomial space,
given by
(74)

(ϕh , w? )K = (ϕh , w)K

∀ϕh ∈ ΠpK −1 (K).

12

WOOPEN, BALAN, MAY, SCHÜTZ

(a) Coarse mesh

(b) Adapted mesh

Figure 1. Smoothness sensor for a transonic test case

Hence, w − w? represents the higher order components of the solution (see Fig. 1).
As we use a hierarchical basis [31], this projection is very cheap. By introducing a
threshold S , a decision between mesh-refinement and p-enrichment can be made,
i.e.
(
(75)

SK

< S
≥ S

p-enrichment
mesh-refinement

First promising results for this approach are given in [20].

A COMPARISON OF HDG AND DG FOR TARGET-BASED hp-ADAPTATION

13

5. Results
In the following we compare our in-house HDG and DG solvers in terms of degrees of freedom and runtime. The DG discretization is based on the Lax-Friedrich
and the BR2 [32] fluxes for convective and viscous terms, respectively. Boundary
conditions and target functionals are evaluated in an adjoint-consistent manner
[33]. Both solvers share the same computational framework, so we believe that our
comparison is meaningful.
We apply both solvers to compressible flow problems, including inviscid subsonic,
transonic and supersonic, as well as subsonic laminar flow. In all cases we show
results for pure mesh-adaptation (p = 1 . . . 4) and hp-adaptation (p = 2 . . . 5).
Please note, that we choose the same relaxation parameters for all test cases,
namely, c0 = 106 , c1 = 103 , n0 = 4 (see Eq. (54)). If artificial viscosity is necessary,
we use 0 = 0.2 and β = 0. The parameters for the adaptation procedure are chosen
as θ = 0.05 and S = 10−6 . This set of parameters seems to be robust for both
DG and HDG for a broad range of test cases. In order to approximate the error
in drag, reference solutions on hp-adapted meshes with more than 2 · 106 degrees
of freedom are used. (Please note, that we refer to ndof w whenever we speak of
degrees of freedom in the following as this is a good measure for the resolution.)
Before we turn our attention to the adaptive computations, we compare runtimes for both methods on a fixed mesh for several polynomial orders. This way,
we can learn which improvement can be expected. In Fig. 2 timings for the assembly
procedure and the iterative solver are given. We show Euler and a Navier-Stokes
computation on a mesh with 2396 elements and 3544 interior faces. We used polynomial orders from p = 0 to p = 6. We want to emphasize that the necessary
Newton and GMRES iterations for both HDG and DG were comparable. As there
are more faces than elements, DG is faster than HDG for p = 0 and p = 1. However, already for p = 2 HDG outruns DG. At p = 6 there is a ratio of 2.5 for the
Euler test case and 2.1 for the Navier-Stokes test case. For the Euler test case it
is interesting to note, that the iterative solver dominates the computational time
for DG. For HDG it is the other way around. This has two reasons: firstly, the
assembly is more involved due to the local solves; and secondly, the global system
is considerably smaller for higher polynomial degree. In the case of a Navier-Stokes
computation, the DG assembly takes over the dominating part as the lifting operators are very expensive to compute. The time for the HDG assembly is also
increased which is among other things due to the introduction of the gradient. For
both HDG and DG, the time spent in the iterative solver is comparable for both
Euler and Navier-Stokes.
5.1. Subsonic Inviscid Flow over the NACA 0012 Airfoil. In the first test
case, we consider subsonic inviscid flow over the NACA 0012 airfoil which is defined
by

√
y = ±0.6 0.2969 x − 0.1260x − 0.3516x2 + 0.2843x3 − 0.1015x4

with x ∈ [0, 1]. Using this definition, the airfoil would have a finite trailing edge
thickness of .252 %. In order to obtain a sharp trailing edge we modify the x4
coefficient, i.e.

√
y = ±0.6 0.2969 x − 0.1260x − 0.3516x2 + 0.2843x3 − 0.1036x4 .

The flow is characterized by a free stream Mach number of Ma∞ = 0.5 and an
angle of attack of α = 2◦ . In Fig. 3 the baseline mesh for the Euler test cases

14

WOOPEN, BALAN, MAY, SCHÜTZ

Assembly

40

Assembly

200

Linear Solve

Linear Solve

t

t

150
100

20

50
0

0

1

2

3
p

4

5

6

(a) Euler

0

0

1

2

3
p

4

5

6

(b) Navier-Stokes

Figure 2. Runtime comparison of the hybridized and nonhybridized DG method for a fixed mesh and varying polynomial
degree (Solid: HDG, Shaded: DG).
(subsonic and transonic) can be seen. It consists of 719 elements and its far field is
over a 1000 chords away.
Admissible target functionals defined on the boundary for the Euler equations
are given by the weighted pressure along wall boundaries, i.e.
Z
(76)
J(w) =
ψ · (pn) dσ
∂Ω

where n is the outward pointing normal. By using ψ =
1
C∞

T

1
C∞

T

(cos α, sin α)

or

ψ =
(− sin α, cos α) along wall boundaries and 0 otherwise, the functional
represents the pressure drag coefficient cD or the pressure lift coefficient cL , respectively. C∞ is a normalized reference value defined by C∞ = 12 γMa2∞ p∞ l. Here, l is
the chord length of the airfoil.
In Fig. 4a, a purely h-adapted mesh can be seen. The most refined regions are
the leading and trailing edge. The former is of importance as the flow experiences
high gradients towards the stagnation point. Refinement of the latter is necessary
due to the sharp trailing edge and the slip-wall boundary conditions. As soon as
the error in these two regions is sufficiently low, other elements close to the airfoil
get refined as well. For the hp-adapted mesh (see Fig. 4b) the leading and trailing
edge are refined as well. All other regions, however, undergo mostly p-enrichment.
In terms of degrees of freedom, both HDG and DG show similar results. For
all computations it takes some adaptations until the critical regions, leading and
trailing edge, are resolved. From this point on, one can see the benefit of a higher
order discretization: the error drops significantly faster with respect to degrees of
freedom and computational time (see Fig. 5 and 6).
In Tbl. 1, we give the runtime ratios for a fixed error level (we always choose
the minimum level attained). For p = 1, HDG and DG are comparable in both
runtime and nonzero entries. From p = 2 on, HDG is faster than DG and needs
less memories for the global system. Here, it is important to note that the adjoint is approximated with p + 1, so that already for lower orders HDG is faster.

A COMPARISON OF HDG AND DG FOR TARGET-BASED hp-ADAPTATION

15

Figure 3. Baseline mesh with 719 elements for inviscid computations
p
1
2
3
4
tDG /tHDG
1.293 3.799 3.551 3.092
nnz,DG /nnz,HDG 1.123 2.150 3.248 4.281
Table 1. Runtime and nonzero ratios for a fixed
(Ma∞ = 0.5, α = 2◦ )

hp
2.209
4.601
error level

Furthermore, we list the ratio of nonzero entries in the system matrix in order to
compare the memory requirements of both methods. As expected, the ratio grows
with increasing polynomial degree in favor of HDG.
5.2. Transonic Inviscid Flow over the NACA 0012 Airfoil. Next, we turn
our attention to transonic flow which develops more complex features (e.g. compression shocks) compared to the subsonic regime. The flow is characterized by a
free stream Mach number of Ma∞ = 0.8 and an angle of attack of α = 1.25◦ .
In Fig. 7a a purely h-adapted mesh can be seen. The adjoint sensor detects all
regions of relevance for the drag: the upper shock, the leading and trailing edge,
and the lower weak shock. Further refinement is added upstream of the shock,
where the adjoint has steep gradients and thus needs higher resolution. In the case
of hp-adaptation, the mesh-refinement is stronger confined to the shock region and
the trailing edge. The other features undergo p-enrichment.
As expected, both methods show a similar accuracy for a given number of degrees
of freedom. The computations with p = 2 . . . 4 outperform p = 1 but are comparable
to each other. hp-adaptation shows very good results which is due to the accurate
prediction of the solution smoothness (see Fig. 8 and 9).
For this test case, HDG is faster than DG from p = 1 on (see Tbl. 2). The hpadaptive run is more than 2.5 times faster. The ratio of necessary nonzero entries
attains its highest value for p = 4. In the case of hp-adaptation, this ratio is not as
high, as in the shock region a lot of elements with p = 2 exist.

16

WOOPEN, BALAN, MAY, SCHÜTZ

(a) Pure h-adapation (p = 2)

(b) hp-adapation (p = 2 . . . 5)

Figure 4. Adapted meshes for the subsonic Euler test case
(Ma∞ = 0.5, α = 2◦ )
p
1
2
3
4
tDG /tHDG
2.636 1.607 2.154 2.362
nnz,DG /nnz,HDG 1.237 1.807 3.098 4.437
Table 2. Runtime and nonzero ratios for a fixed
(Ma∞ = 0.8, α = 1.25◦ )

hp
2.628
2.411
error level

A COMPARISON OF HDG AND DG FOR TARGET-BASED hp-ADAPTATION

10−3

10−5

∆cD

∆cD

10−3

17

p=1
p=2
p=3
p=4
hp

10−7

10−5

10−7

103

104
ndof w

105

p=1
p=2
p=3
p=4
hp

103

104
ndof w

(a) HDG

105

(b) DG

Figure 5. Drag convergence with respect to degrees of freedom
(Ma∞ = 0.5, α = 2◦ )

10−5

10−7
10−2

10−3

∆cD

∆cD

10−3

p=1
p=2
p=3
p=4
hp

10−5

10−7

10−1

100
t/t0

(a) HDG

101

102

10−2

p=1
p=2
p=3
p=4
hp
10−1

100
t/t0

101

102

(b) DG

Figure 6. Drag convergence with respect to time (Ma∞ = 0.5, α = 2◦ )

5.3. Supersonic Inviscid Flow over the NACA 0012 Airfoil. Now, we consider supersonic, inviscid flow with a free stream Mach number of Ma = 1.4. We
choose the angle of attack to be α = 0◦ . In contrast to the transonic test case, a
so-called bow shock develops in front of the airfoil.
A purely h-adapted mesh is given in Fig. 10a. Several refinement regions can
be identified: the bow shock, leading and trailing edge, and a feature in the adjoint solution emanating from the trailing edge in the upstream direction. During
hp-adaptation, the latter feature and the leading edge are being treated with penrichment.
Both HDG and DG show a similar behavior in terms of error reduction with
respect to degrees of freedom (see Fig. 11). A comparison between different polynomial orders is rather inconclusive which might be due to the lack of smoothness
of this test case. Furthermore, we note that the convergence behavior is not monotonic (see for example p = 3 in Fig. 11). The error reduction with respect to

18

WOOPEN, BALAN, MAY, SCHÜTZ

(a) Pure h-adaptation (p = 2)

(b) hp-adaptation (p = 2 . . . 5)

Figure 7. Adapted meshes for the transonic Euler test case
(Ma∞ = 0.8, α = 1.25◦ )
computational time is similar (see Fig. 12). In both measures, the hp-adaptive run
gives the overall best solution.
Considering the time to solution, HDG is approximately 50% faster than DG
(see Tbl. 3). The hp-adaptive run is more than 3.5 times faster. Similar to the
transonic test case, the ratio of necessary nonzero entries attains its highest value
for p = 4. In the case of hp-adaptation, this ratio is not as high, as in the shock
region a lot of elements with p = 2 exist.
5.4. Subsonic Laminar Flow over the NACA 0012 Airfoil. Finally, we consider viscous flow in the subsonic regime. The free stream Mach number is Ma∞ =

10−1

10−1

10−2

10−2

10−3

10−4

10−5
103

∆cD

∆cD

A COMPARISON OF HDG AND DG FOR TARGET-BASED hp-ADAPTATION

p=1
p=2
p=3
p=4
hp

10−3

10−4

104
ndof w

10−5
103

105

19

p=1
p=2
p=3
p=4
hp
104
ndof w

(a) HDG

105

(b) DG

10−1

10−1

10−2

10−2

10−3

10−4

10−5
10−1

∆cD

∆cD

Figure 8. Drag convergence with respect to degrees of freedom
(Ma∞ = 0.8, α = 1.25◦ )

p=1
p=2
p=3
p=4
hp

10−3

10−4

101

100
t/t0

(a) HDG

102

10−5
10−1

p=1
p=2
p=3
p=4
hp
101

100

102

t/t0

(b) DG

Figure 9. Drag convergence with respect to time (Ma∞ = 0.8,
α = 1.25◦ )
p
1
2
3
4
tDG /tHDG
2.817 1.335 1.751 1.441
nnz,DG /nnz,HDG 1.212 2.132 3.429 4.030
Table 3. Runtime and nonzero ratios for a fixed
(Ma∞ = 1.4, α = 0◦ )

hp
3.554
3.806
error level

0.5, the angle of attack α = 1◦ and the Reynolds number Re = 5000. Due to the
latter, a thin boundary layer develops around the airfoil.
The baseline mesh for the Navier-Stokes test case is more refined around the
airfoil such that the boundary layer is correctly captured (see Fig. 13). It consists
of 1781 elements and its far field is over a 1000 chords away.

20

WOOPEN, BALAN, MAY, SCHÜTZ

(a) Pure h-adapation (p = 2)

(b) hp-adapation (p = 2 . . . 5)

Figure 10. Adapted meshes for the supersonic Euler test case
(Ma∞ = 1.4, α = 0◦ )
Admissible target functionals defined on the boundary for the Navier-Stokes
equations are given by the weighted boundary flux (including pressure and shear
forces) along wall boundaries, i.e.
Z
(77)
J (w, ∇w) =
ψ · (pn − τ n) dσ
∂Ω

where n is the outward pointing normal. Here, ψ is nonzero only on wall boundaries.
T
T
By using ψ = C1∞ (cos α, sin α) or ψ = C1∞ (− sin α, cos α) along wall boundaries
and otherwise 0, the functional represents the viscous drag coefficient cD or the
viscous lift coefficient cL , respectively.

10−2

10−2

10−3

10−3

10−4

10−5

10−6
103

∆cD

∆cD

A COMPARISON OF HDG AND DG FOR TARGET-BASED hp-ADAPTATION

p=1
p=2
p=3
p=4
hp

10−4

10−5

104
ndof w

10−6
103

105

(a) HDG

21

p=1
p=2
p=3
p=4
hp
104
ndof w

105

(b) DG

Figure 11. Drag convergence with respect to degrees of freedom
(Ma∞ = 1.4, α = 0◦ )

p=1
p=2
p=3
p=4
hp

∆cD

10−3

10−4

10−5

10−6
10−1

10−2

p=1
p=2
p=3
p=4
hp

10−3
∆cD

10−2

10−4

10−5

100

101
t/t0

(a) HDG

102

10−6
10−1

100

101

102

t/t0

(b) DG

Figure 12. Drag convergence with respect to time (Ma∞ = 1.4, α = 0◦ )
Both the h-adapted mesh (see Fig. 14a) and the hp-adapted mesh (see Fig. 14b)
undergo refinement within the boundary layer and the wake region. The mesh
refinement for the hp-adaptive run is however more restricted to the leading edge
region where the boundary layer develops. Further downstream, p-enrichment is
used as soon as the necessary mesh-resolution is reached.
In terms of accuracy versus degrees of freedom, HDG and DG perform comparably well. The higher the polynomial degree the more accurate and efficient the
computations are for both HDG and DG (see Fig. 15 and 16). The difference between hp, p = 3 and p = 4 is not as big, though. This might lead to the conclusion
that isotropic mesh refinement is not longer efficiently applicable in cases involving
strong gradients. Anisotropic refinement might be able to remedy this problem so
that the full potential of hp-adaptation can be attained.
Concerning the timings, we can see a similar trend as in the previous test cases
(see Tbl. 4). For p = 1 . . . 4, HDG is more than twice as fast. The hp-adaptive

22

WOOPEN, BALAN, MAY, SCHÜTZ

Figure 13. Baseline mesh with 1781 elements for viscous computations
p
1
2
3
4
tDG /tHDG
2.244 2.572 2.397 2.288
nnz,DG /nnz,HDG 1.196 2.235 3.472 4.531
Table 4. Runtime and nonzero ratios for a fixed
(Ma∞ = 0.5, α = 1◦ , Re = 5000)

hp
3.122
3.463
error level

HDG computation is even three times as fast compared to the DG run. For p ≥ 2
the savings in nonzero entries for HDG become significant.
6. Conclusion and Outlook
We presented an adjoint-based hp-adaptation methodology and compared it for
hybridized and non-hybridized discontinuous Galerkin methods. hp-adaptation
proved to be superior to pure h-adaptation if discontinuous or singular flow features
were involved. In all cases, a higher polynomial degree turned out to be beneficial.
We showed that one can expect HDG to be better than DG in terms of runtime and
memory requirements from p = 2 on for a broad range of test cases. The results
are promising enough to pursue further comparisons with more complex test cases.
We plan to extend our computational framework to three dimensional problems
[34]. Then, adaptivity will play an even more crucial role, as the problem size
increases drastically compared to the two dimensional case.
Furthermore, compressible flows are often dominated by anisotropic features,
such as shocks or very thin boundary layers. Thus, taking this anisotropy into
account during adaptation is crucial. We plan to incorporate anisotropic adaptation
within future projects.
Acknowledgments
Financial support from the Deutsche Forschungsgemeinschaft (German Research
Association) through grant GSC 111 is gratefully acknowledged.

A COMPARISON OF HDG AND DG FOR TARGET-BASED hp-ADAPTATION

23

(a) Pure h-adaptation (p = 2)

(b) hp-adaptation (p = 2 . . . 5)

Figure 14. Adapted meshes for the Navier-Stokes test case
(Ma∞ = 0.5, α = 1◦ , Re = 5000)
References
[1] F. Bassi, S. Rebay, A High-Order Accurate Discontinuous Finite-Element
Method for the Numerical Solution of the Compressible Navier-Stokes Equations, Journal of Computational Physics 131 (1997) 267–279.
[2] B. Cockburn, C.-W. Shu, TVB Runge-Kutta Local Projection Discontinuous
Galerkin Finite Element Method for Conservation Laws II: General Framework, Mathematics of Computation 52 (1988) 411–435.
[3] D. Arnold, F. Brezzi, B. Cockburn, L. Marini, Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems, SIAM Journal on Numerical
Analysis 39 (2002) 1749–1779.

WOOPEN, BALAN, MAY, SCHÜTZ

10−2

10−2

10−3

10−3

10−4

∆cD

∆cD

24

p=1
p=2
p=3
p=4
hp

10−5

10−6

10−4

10−5

104

10−6

105

ndof w

p=1
p=2
p=3
p=4
hp
104

(a) HDG

105

ndof w

(b) DG

10−2

10−2

10−3

10−3

10−4

10−5

10−6
10−1

∆cD

∆cD

Figure 15. Drag convergence with respect to degrees of freedom
(Ma∞ = 0.5, α = 1◦ , Re = 5000)

p=1
p=2
p=3
p=4
hp

10−4

10−5

100

101
t/t0

(a) HDG

102

103

10−6
10−1

p=1
p=2
p=3
p=4
hp
100

101
t/t0

102

103

(b) DG

Figure 16. Drag convergence with respect to time (Ma∞ = 0.5,
α = 1◦ , Re = 5000)

[4] B. F. de Veubeke, Displacement and Equilibrium Models in the Finite Element
Method, chap. 9, Wiley New York, 145–197, 1965.
[5] B. Cockburn, J. Gopalakrishnan, A Characterization of Hybridized Mixed
Methods for Second Order Elliptic Problems, SIAM Journal on Numerical
Analysis 42 (1) (2004) 283–301.
[6] B. Cockburn, J. Gopalakrishnan, R. Lazarov, Unified Hybridization of Discontinuous Galerkin, Mixed, and Continuous Galerkin Methods for Second
Order Elliptic Problems, SIAM Journal on Numerical Analysis 47 (2) (2009)
1319–1365.
[7] N. Nguyen, J. Peraire, B. Cockburn, An Implicit High-Order Hybridizable
Discontinuous Galerkin Method for Linear Convection-Diffusion Equations,
Journal of Computational Physics 228 (9) (2009) 3232–3254.

A COMPARISON OF HDG AND DG FOR TARGET-BASED hp-ADAPTATION

25

[8] N. Nguyen, J. Peraire, B. Cockburn, An Implicit High-Order Hybridizable
Discontinuous Galerkin Method for Nonlinear Convection-Diffusion Equations,
Journal of Computational Physics 228 (2009) 8841–8855.
[9] J. Peraire, N. Nguyen, B. Cockburn, A Hybridizable Discontinuous Galerkin
Method for the Compressible Euler and Navier-Stokes Equations, AIAA Paper
2010-0363, 2010.
[10] A. Huerta, A. Angeloski, X. Roca, J. Peraire, Efficiency of High-Order Elements for Continuous and Discontinuous Galerkin Methods, Int. J. Numer.
Methods Eng. 96 (9) (2013) 529–560.
[11] R. Hartmann, P. Houston, Adaptive Discontinuous Galerkin Finite Element
Methods for the Compressible Euler Equations, Journal of Computational
Physics 183 (2) (2002) 508–532.
[12] D. Venditti, D. Darmofal, Grid Adaptation for Functional Outputs: Application to Two-Dimensional Inviscid Flows, Journal of Computational Physics
176 (1) (2002) 40–69.
[13] D. Venditti, D. Darmofal, Anisotropic Grid Adaptation for Functional Outputs: Application to Two-Dimensional Viscous Flows, Journal of Computational Physics 187 (1) (2003) 22–46.
[14] R. Hartmann, Adaptive Discontinuous Galerkin Methods with ShockCapturing for the Compressible Navier–Stokes Equations, International Journal for Numerical Methods in Fluids 51 (9-10) (2006) 1131–1156.
[15] K. Fidkowski, D. Darmofal, Review of Output-Based Error Estimation and
Mesh Adaptation in Computational Fluid Dynamics, AIAA Journal 49 (4)
(2011) 673–694.
[16] G. Giorgiani, S. Fernández-Méndez, A. Huerta, Hybridizable Discontinuous
Galerkin p-Adaptivity for Wave Propagation Problems, International Journal
for Numerical Methods in Fluids .
[17] J. Schütz, G. May, A Hybrid Mixed Method for the Compressible Navier–
Stokes Equations, Journal of Computational Physics 240 (2013) 58–75.
[18] J. Schütz, G. May, An Adjoint Consistency Analysis for a Class of Hybrid
Mixed Methods, IMA Journal of Numerical Analysis .
[19] M. Woopen, G. May, J. Schütz, Adjoint-Based Error Estimation and Mesh
Adaptation for Hybridized Discontinuous Galerkin Methods, arXiv preprint
arXiv:1309.3649 .
[20] A. Balan, M. Woopen, G. May, Adjoint-Based hp-Adaptation for a Class of
High-Order Hybridized Finite Element Schemes for Compressible Flows, AIAA
Paper 2013-2938, 2013.
[21] E. Godlewski, P. Raviart, Numerical Approximation of Hyperbolic Systems of
Conservation Laws, no. Nr. 118 in Applied Mathematical Sciences, Springer,
ISBN 9780387945293, 1996.
[22] Y. Chen, B. Cockburn, Analysis of Variable-Degree HDG Methods for
Convection–Diffusion Equations. Part I: General Nonconforming Meshes, IMA
Journal of Numerical Analysis 32 (4) (2012) 1267–1293.
[23] B. Cockburn, J. Gopalakrishnan, Error Analysis of Variable Degree Mixed
Methods for Elliptic Problems via Hybridization, Mathematics of Computation
74 (252) (2005) 1653–1677.
[24] D. Mavriplis, A. Jameson, Multigrid Solution of the Navier-Stokes Equations
on Triangular Meshes, AIAA Journal 28 (8) (1990) 1415–1425.

26

WOOPEN, BALAN, MAY, SCHÜTZ

[25] S. Balay, J. Brown, K. Buschelman, W. Gropp, D. Kaushik,
M. Knepley, L. McInnes, B. Smith, H. Zhang, PETSc Web Page,
http://www.mcs.anl.gov/petsc, 2013.
[26] S. Balay, J. Brown, K. Buschelman, V. Eijkhout, W. Gropp, D. Kaushik,
M. Knepley, L. McInnes, B. Smith, H. Zhang, PETSc Users Manual, Tech.
Rep. ANL-95/11 - Revision 3.4, Argonne National Laboratory, 2013.
[27] W. Dörfler, A Convergent Adaptive Algorithm for Poisson’s Equation, SIAM
Journal on Numerical Analysis 33 (3) (1996) 1106–1124.
[28] M. Ceze, K. Fidkowski, Output-Driven Anisotropic Mesh Adaptation for Viscous Flows using Discrete Choice Optimization, AIAA Paper 2010-0170, 2010.
[29] L. Wang, D. J. Mavriplis, Adjoint-Based h-p Adaptive Discontinuous Galerkin
Methods for the 2D Compressible Euler Equations, Journal of Computational
Physics 228 (2009) 7643–7661.
[30] P. Persson, J. Peraire, Sub-cell Shock Capturing for Discontinuous Galerkin
Methods, AIAA Paper 2006-0112, 2006.
[31] M. Dubiner, Spectral Methods on Triangles and Other Domains, J. Sci. Comput. 6 (4) (1991) 345–390.
[32] F. Bassi, A. Crivellini, S. Rebay, M. Savini, Discontinuous Galerkin Solution of
the Reynolds-Averaged Navier-Stokes and k-ω Turbulence Model Equations,
Computers & Fluids 34 (4) (2005) 507–540.
[33] R. Hartmann, Adjoint Consistency Analysis of Discontinuous Galerkin Discretizations, SIAM Journal on Numerical Analysis 45 (6) (2007) 2671–2696.
[34] M. Woopen, A. Balan, G. May, A Hybridized Discontinuous Galerkin Method
for Three-Dimensional Compressible Flow Problems, AIAA Paper 2014-0938,
2014.

