High order exactly divergence-free Hybrid Discontinuous Galerkin Methods
for unsteady incompressible flows
Christoph Lehrenfelda , Joachim Schöberlb
a Institute

for Computational and Applied Mathematics, WWU Münster, Münster, Germany
for Analysis and Scientific Computing, TU Wien, Wien, Austria

arXiv:1508.04245v2 [math.NA] 26 Apr 2016

b Institute

Abstract
In this paper we present an efficient discretization method for the solution of the unsteady incompressible
Navier-Stokes equations based on a high order (Hybrid) Discontinuous Galerkin formulation. The crucial
component for the efficiency of the discretization method is the disctinction between stiff linear parts and
less stiff non-linear parts with respect to their temporal and spatial treatment.
Exploiting the flexibility of operator-splitting time integration schemes we combine two spatial discretizations which are tailored for two simpler sub-problems: a corresponding hyperbolic transport problem and
an unsteady Stokes problem.
For the hyperbolic transport problem a spatial discretization with an Upwind Discontinuous Galerkin
method and an explicit treatment in the time integration scheme is rather natural and allows for an efficient
implementation. The treatment of the Stokes part involves the solution of linear systems. In this case a discretization with Hybrid Discontinuous Galerkin methods is better suited. We consider such a discretization
for the Stokes part with two important features: H(div)-conforming finite elements to garantuee exactly
divergence-free velocity solutions and a projection operator which reduces the number of globally coupled
unknowns. We present the method, discuss implementational aspects and demonstrate the performance on
two and three dimensional benchmark problems.
Keywords: Navier-Stokes equations, Hybrid Discontinuous Galerkin Methods, H(div)-conforming Finite
Elements, exactly divergence-free, Operator-Splitting, reduced stabilization

1. Introduction
1.1. Problem statement
We consider the numerical solution of the unsteady incompressible Navier Stokes equations in a velocitypressure formulation:
 ∂
in Ω
∂t u + div(−ν∇u + u⊗u + pI) = f
(1)
div u = 0 in Ω
with boundary conditions u = uD on ΓD ⊂ ∂Ω and (ν∇u − pI) · n = 0 on Γout = ∂Ω \ ΓD . Here, ν = const
is the kinematic viscosity, u the velocity, p the pressure, and f is an external body force. We consider
a discretization to (1) with high order (Hybrid) Discontinuous Galerkin (DG) finite element methods for
complex geometries with underlying meshes consisting of possibly curved tetrahedra, hexahedra, prisms and
pyramids.
Email addresses: christoph.lehrenfeld@uni-muenster.de (Christoph Lehrenfeld), joachim.schoeberl@tu-wien.ac.at
(Joachim Schöberl)

Preprint submitted to Comp. Meth. Appl. Mech. Eng.

April 27, 2016

1.2. Literature
Since its introduction in the paper of Reed and Hill [1], DG methods have been developed and used
for hyperbolic problems [2, 3, 4, 5, 6] and later extended to second-order elliptic (and parabolic) problems
[7, 8, 9, 10]. DG methods are specifically popular for flow problems [11, 12]. For incompressible flows
different finite element methods have been discussed in the literature [13, 14, 15]. Among these methods
several DG formulations have been considered, e.g. [16, 17, 18, 19, 20].
DG methods provide a flexibility which can be utilized for different purposes. In our case, the motivation
to consider this type of discretizations is twofold. First, relevant flow problems that can be modeled with
the incompressible Navier-Stokes equations are often convection dominated. Using the Upwind mechanism
DG methods offer a natural way to devise stable discretizations of (dominating) convection. Secondly,
abandoning H 1 -conforming finite element spaces, DG methods allow to consider (only) H(div)-conforming
finite elements (see [21]) for (Navier-)Stokes problems. This is attractive as it facilitates the design of a
discretization with exactly solenoidal solution which in turn implies energy-stability of the discretization for
the Navier-Stokes problem, cf. [19, 20, 22].
Compared to Continuous Galerkin (CG) methods the number of degrees of freedom of Discontinuous
Galerkin methods increases significantly. This drawdack is often outwayed by the advantages of the method.
However, when it comes to solving linear systems Discontinuous Galerkin methods suffer most from drastically increased globally coupled degrees of freedom. An approach to compensate for this is the concept of
Hybridization where additional unknowns on element interfaces are introduced. This increases the number
of degrees of freedom, but introduces two advantages. First, the global couplings are reduced and secondly,
the structure of the couplings allows to apply static condensation for the element unknowns. The concept has originally been introduced in the context of mixed finite element methods, cf. [21]. In the last
decade Hybridization has also been applied to Discontinuous Galerkin method for a variety of problems,
cf. [23, 24, 25, 26, 27]. We also mention, that recently similar concepts are also known in the literature
under the name “Hybridized Weak Galerkin” methods [28] where a slightly different framework is used to
derive the methods. With the emphasis to treat general polyhedral meshes “Hybrid High-Order” (HHO)
methods [29, 30, 31] have recently been introduced where a combination of Hybrid (or hybridized) methods
and higher order spaces is considered.
Concerning HDG discretizations for incompressible flows, we also mention the papers [32, 26, 33]. Further, different approaches to implement exact incompressible finite element solutions to the Stokes problem
using Hybridization have been investigated in [34, 35, 36]. In [37] the Stokes-Brinkman problem, the Stokes
problem with an additional zero order term, has been discretized using a hybridized H(div)-conforming
finite element formulation.
An interesting and often praised aspect of Hybrid mixed finite element methods is the fact that a postprocessing step can be used to reconstruct interior unknowns of an increased higher order accuracy for
elliptic problems [21]. The same is also possible for HDG methods in mixed formulation which is another
advantage over conventional DG methods, see e.g. [25]. For elliptic problems an accuracy of order k + 2
in the volume can be achieved if polynomials of order k are used on the element interfaces. One main
contribution in this paper is the introduction of a projection operator which achieves the same effect. In
[22] we already discussed this operator which has recently also been addressed under the name “reduced
stabilization” in [38, 39, 40, 41]. In these papers the construction of the corresponding operator is only
possible in two dimensions using special integration rules. We explain how the operator can be implemented
in a more general setting.
The efficiency and practicability of HDG methods is rarely addressed in the literature. Interesting
exception are [42, 43] where the computational cost of the method is compared to CG and DG methods.
The nature of diffusive (viscous) and convective terms appearing in convection-diffusion-type problems
have a substantially different character. Discretizations of problems involving only one of the two mechanisms
would typically lead to different spatial and temporal discretizations. It is therefore often desirable to consider
operator-splitting time discretization schemes which allow for a separate treatment of the different operators
as in [44, 45, 46, 47, 48, 49]. In these papers the spatial discretization for the stiff (diffusion/Stokes) and the
non-stiff (convection) operators is typically the same. We consider a different treatment of the operators
with respect to their temporal and spatial discretization.
2

1.3. The concept: Efficiency through operator splitting
A crucial ingredient in the considered discretization is the fact that the convection and the Stokes
problem are separated by means of an operator-splitting method. The operator-splitting method is chosen
such that only operator evaluations of the convection operator are required so that the time integration
scheme is explicit in terms of the convection operator. This is often affordable as the time step restrictions
following from the convection operator are the least restrictive ones in the Navier-Stokes problem. Moreover,
the convection operator is non-linear, s.t. the set up of linear systems of equations and corresponding
preconditioners or solvers would have to take place every time step which renders implicit approaches for
the convection very expensive. The Stokes operator is dealt with implicitly, i.e. it appears in linear systems
in every time step. Due to the differential-algebraic structure of the Navier-Stokes equations this is necessary
w.r.t. the pressure and the incompressibility constraint. As completely explicit handling of the viscosity
terms would introduce severe time step restrictions it is also advisable to treat viscous forces implicit. Note
that the Stokes operator is time-independent such that the setup of linear systems and preconditioners or
solvers can be done once and re-used in every time step.
The spatial discretizations for the different operators are designed differently as the treatment of the
convection term can be optimized for operator evaluations while the discretization of the Stokes operator
has to provide an efficient handling of implicit solution steps, i.e. the solution of corresponding linear
systems. This different treatment reflects in the use of two different finite element spaces which are used:
An H(div)-conforming Hybrid DG space for the velocity-pressure-pair for solving Stokes problems and a
DG space for handling convection.
1.4. Main contributions and structure of this paper
In this paper we introduce a new discretization method for the incompressible Navier-Stokes equations.
The discretization is based on a decomposition of the problem into the (unsteady) Stokes problem and
a hyperbolic transport problem. This decomposition reflects in the use of appropriate operator splitting
methods in the time discretization and the use of different finite element spaces for the different spatial
operators.
While we use a rather standard DG formulation for the spatial discretization of the hyperbolic transport
problem, we consider a new HDG formulation for the spatial discretization of the Stokes-type problem. We
use H(div)-conforming functions of element-wise polynomials of degree k for the velocity and discontinuous
pressure functions of degree k − 1. Due to the introduction of additional unknowns of degree k − 1 on
the facets for the approximation of the tangential trace of the velocity static condensation allows to reduce
the unknowns and their couplings significantly. The resulting globally coupled unknowns correspond to the
approximation of the normal trace of the velocity on the facets by order k functions, the tangential trace of
the velocity by order k − 1 functions and the approximation of the pressure field by piecewise constants.
The main contributions of this paper are:
1. Introduction of a new H(div)-conforming high order accurate HDG method with a projected jumps
formulation (also known as reduced stabilization or reduced-order HDG) for the solution of Stokes-type
problems.
2. Presentation of a combined spatial discretization for the Navier-Stokes equations based on a standard
Upwind DG formulation for the hyperbolic transport problem with the new HDG method for StokesBrinkman problems.
3. Discussion of operator-splitting time integration schemes which restrict solutions of linear systems to
Stokes-Brinkman problems.
In section 2 we discuss the discretization of spatial operators. The discussion is divided into two parts, the
discretization of the Stokes-Brinkman problem, the problem which involves all relevant spatial operators
except for the convection and the discretization of the convection operator. For the former part we consider
an HDG discretization, for the latter part a standard DG discretization. In section 3 operator-splitting
time integration schemes are discussed which are tailored for such a situation. Finally, in section 4 we give
numerical examples which demonstrate the accuracy and the performance of the method.
3

1.5. Preliminaries
Before describing the methods, we introduce some basic notation and assumptions: Ω is an open bounded
domain in Rd with a Lipschitz boundary Γ. It is decomposed into a shape regular partition Th of Ω consisting
of elements T which are (curved) simplices, quadrilaterals, hexahedrals, prisms or pyramids. For ease of
presentation we assume that all elements are not curved and consider only homogeneous Dirichlet boundary
conditions. The element interfaces and element boundaries coinciding
boundary are called
S with the domain
S
facets. The set of those facets F is denoted by Fh and there holds T ∈Th ∂T = F ∈Fh F .
In the sequel we distinguish functions with support only on facets indicated by a subscript F and those
with support also on the volume elements which is indicated by a subscript T . Compositions of both types
are used for the HDG discretization of the velocity which is denoted by u = (uT , uF ).
For vector-valued functions the superscripts t denotes the application of the tangential projection: v t =
v − (v·n) · n ∈ Rd . The index k which describes the polynomial degree of the finite element approximation
at many places through out the paper is an arbitrary but fixed positive integer number.
We identify finite
functions u ∈ Xh with their representation in terms of coefficient vectors
PNelement
X
ui ϕi for a corresponding basis {ϕi } of Xh and NX = dim(Xh ). A (generic) bilinear
u ∈ RNX , s.t. u = i=1
Y
form Gh : Xh × Yh → R is identified with the matrix G ∈ RNY ×NX , s.t. Gi,j = Gh (ϕX
j , ϕi ).
2. DG/HDG spatial discretization
In this section we introduce the spatial discretization for the Stokes operator, the convection operator
and transfer operations between both. First, we introduce the H(div)-conforming Hybrid DG discretization
of the Stokes-Brinkman problem, i.e. the stationary Navier-Stokes problem without convection and an
additional zero order term in section in section 2.1. This discretization is improved significantly. Further
on, a modification of the discretization using the idea of projected jumps (also known as reduced stabilization
or reduced-order HDG) is presented in section 2.2 including an explanation of how the projection operator
can be realized. For the convection part of the Navier-Stokes problem we consider a DG discretization using
standard approaches. This is discussed in section 2.3. As the discretization spaces for the Stokes part and
the convection part are different, we present transfer operations between the spaces in section 2.4 which
allow us to finally formulate the semi-discrete problem in section 2.5.
2.1. H(div)-conforming HDG formulation of the Stokes-Brinkman problem
In this part we consider the discretization of the Stokes part of the Navier-Stokes problem. For simplicity
we restrict the discussion to homogeneous Dirichlet boundary conditions. We present the method and
elaborate on important properties. For an error analysis of the method we refer to [22]. The reaction
term τ −1 corresponding to an inertia term stemming from an implicit time integration scheme is further
incorporated. The resulting problem is known as the (stationary) Stokes-Brinkman problem:
 −1
τ u + div(−ν∇u + pI) = f
in Ω
(2)
div u = 0
in Ω
We first introduce the finite element spaces, followed by the definition of the bilinear forms corresponding
to the involved operators.
2.1.1. H(div)-conforming Finite Elements for Stokes.
Following [19] a DG formulation for the incompressible Navier Stokes equations which is locally conservative and energy-stable at the same time has to provide discrete solutions which are exactly divergence-free.
This can be achieved with a suitable pair of finite element spaces. We consider the use of H(div)-conforming
Finite Element spaces for the velocity u. H(div)-conformity requires that every discrete function uh is in
H(div, Ω) = {u ∈ [L2 (Ω)]d : div u ∈ L2 (Ω)}
4

We use piecewise polynomials, so that on each element the functions are in H(div, T ). For global conformity,
continuity of the normal component is necessary, resulting in
Y
Wh := {uT ∈
[P k (T )]d , [[uT ·n]]F = 0 ∀ F ∈ Fh } ⊂ H(div, Ω), NW := dim(Wh ),
(3)
T ∈Th

with [[·]]F the usual jump operator and P k the space of polynomials up to degree k. We refer to [21] for
details on the construction of H(div)-conforming finite elements such as the Brezzi-Douglas-Marini finite
element.
The appropriate Finite Element space for the pressure is the space of piecewise polynomials which are
discontinuous and of one degree less:
Y
Qh :=
P k−1 (T ), NQ := dim(Qh ).
(4)
T ∈Th

This velocity-pressure pair Wh /Qh fulfills
div(Wh ) = Qh .
The crucial point of this choice of the velocity-pressure pair is the property:
If a velocity uT ∈ Wh is weakly incompressible, it is also strongly incompressible:
Z
div(uT )q dx = 0 ∀qh ∈ Qh ⇔ div(uT ) = 0 in Ω.

(5)

(6)

Ω

The benefit of (6) is twofold: First, it allows to show energy-stability for the incompressible Navier-Stokes
equations, cf. section 2.5. Secondly, error estimates for the velocity error can be derived which are independent of the pressure field, we comment on this in remark 5 below.
As solutions of the incompressible Navier Stokes equations are [H 1 (Ω)]d × L2 (Ω)-regular and tangential
continuity is not imposed as an essential condition on the finite element space the discrete formulation has
to incorporate the tangential continuity weakly. We do this with a corresponding (Hybrid) DG formulation
for the tangential components across element interfaces
Remark 1 (Reduced spaces). Due to (6) solutions of the discretized (Navier)-Stokes problem will be
exactly divergence-free velocity fields. This a priori knowledge can be exploited. The basis for the space Wh
can be constructed in such a way such that we can discard certain higher order basis functions (with non-zero
divergence) that will have no contribution. A corresponding basis is introduced in [50] in the context of a
set of higher order basis functions which fulfill an exact sequence property. The resulting space Whred has
div(Whred ) = Qred
h := {p|T = const : T ∈ Th } such that also most of the degrees of freedom of the pressure
space can be discarded. The reduction of basis functions for Wh and Qh is explained in more detail in [22,
Chapter 2.2].
2.1.2. The HDG space for the velocity.
To (weakly) enforce continuity we apply a discontinuous Galerkin (DG) formulation such as the Interior
Penalty method [51]. However, to avoid the full coupling of degrees of freedom of neighboring elements,
we introduce additional unknowns on the skeleton, the facet unknowns, which represent an approximation
of the tangential trace of the solution. The DG formulation is then replaced with a corresponding HDG
formulation, s.t. degrees of freedom of neighboring elements couple only through the facet unknowns. We
note that the resulting space is only “hybrid” in the sense of the tangential unknowns. Further, we do not
consider a hybridization with respect to the pressure, cf. also remark 7.
As normal continuity is already implemented in Wh we only need a DG enforcement of continuity in the
tangential direction and hence only introduce the facet unknown for the tangential direction of the trace:
Y
Fh := {uF ∈
[P k (F )]d , uF · n = 0 }, NF := dim(Fh ).
(7)
F ∈Fh

5

H 1 -conforming (CG)

DG

H(div)-conforming DG

HDG

H(div)-conforming HDG

Figure 1: Tangential and normal continuity for different finite element spaces for the velocity.

Functions in Fh have normal component zero. For the discretization of the velocity field we use the composite
space
Uh := Wh × Fh , NU := dim(Uh ).
(8)
Remark 2 (The role of the facet space). We note that the space Fh is only introduced to allow for a
more efficient handling of linear systems. In section 4.3 we consider a numerical example where a vectorvalued Poisson problem is considered for the HDG space Uh in comparison to other (H(div)-conforming) DG
spaces to illustrate the impact of Hybridization. In case that only explicit applications of discrete operators
are used, the introduction of the space Fh entails no advantages. In this paper, facet variables appear only in
the discretization of the viscous forces. Although the facet variable has no contribution in the discretization
of inertia and pressure forces or the incompressibilty constraint we define the corresponding operations for
u ∈ Uh instead of uT ∈ Wh to simplify the presentation.
2.1.3. Viscous forces.
For ease of presentation we consider a hybridized version of the Interior Penalty method for the discretization of the viscous term. In remark 3 we comment on alternatives. In the hybridized version of the
Interior Penalty method, the usual jump across element interfaces of the Interior Penalty method is replaced
with jumps between element interior and facet unknown (in tangential direction) [[ut ]] = utT − uF from both
sides of a facet. The bilinear form corresponding to the HDG discretization of viscous forces is
Z
X Z
∂uT t
Ah (u, v) :=
ν∇uT : ∇vT dx −
ν
[[v ]] ds
∂n
T
∂T
T ∈T
(9)
Zh
Z
∂vT t
α
−
ν
[[u ]] ds +
ν [[ut ]][[v t ]] ds, u = (uT ,uF ), v = (vT ,vF ) ∈ Uh
∂n
∂T
∂T h
In this bilinear form the four terms have different functions. While the first two terms ensure consistency
(in the sense of Galerkin orthogonality) the third and fourth term are tailored to ensure symmetry (adjoint
consistency) and stability, respectively. Due to continuity of the solution ([[·t ]] = 0) the latter terms also
preserve consistency. For a more detailed introduction we refer to [22, section 2.3.1]. With respect to the
discrete norm
)
(
2
X
∂u
1
T
|||u|||21 :=
k∇uT k2T + k[[ut ]]k2∂T + h
h
∂n ∂T
T ∈Th

which is an appropriately modified version of the discrete norm typically used in the analysis of Standard
Interior Penalty methods, the bilinearform Ah is (for a sufficiently large α) consistent, bounded and coercive.
The coupling through facet unknowns enforces the same kind of continuity as in the Interior Penalty
DG formulation while preserving the following structure in the sparsity pattern: element unknowns are only
coupled with unknowns associated with the same element or unknowns associated with aligned facets.
6

k
quadrilateral
triangle

4

8

16

32

0.305
0.167

0.313
0.190

0.315
0.201

0.315
0.205

Table 1: LBB constant in dependence of k for a single element.

Remark 3 (Interior Penalty and alternatives). A drawback of the (Hybrid) Interior Penalty method
is the fact that the stabilization parameter α depends on the shape regularity of the mesh. Often α is chosen
on the safe side, but as was pointed out in [52] the condition number of arising linear systems increases
with α. In the same paper a hybridized variant of the Bassi-Rebay stabilization method (cf. [4, 5, 53]) for
a scalar Poisson equation has been proposed, see also remark 12. Such a variant is used in the numerical
examples, but as it does not have any further consequences for the remainder of this work, we stick to the
well-known (Hybrid) Interior Penalty method for ease of presentation.
2.1.4. Mass bilinear form.
The HDG mass matrix is defined as
Z
uT vT dx,
MU
(u,
v)
:=
h

u = (uT , uF ), v = (vT , vF ) ∈ Uh .

Ω

(10)

Note that we defined the mass matrix for u ∈ Uh although uF has no contribution, cf. remark 2.
2.1.5. Pressure force and incompressibility constraint.
For the pressure force and the incompressibility constraint we define the bilinearform
X Z
Dh (u, p) :=
−
p div uT dx for u = (uT , uF ) ∈ Uh , p ∈ Qh
T ∈Th

for which the LBB-condition
sup
u∈Uh

(11)

T

Dh (u, p)
≥ cLBB kpkL2 ,
|||u|||1

∀ p ∈ Qh

(12)

holds true for a cLBB independent of the mesh size h, cf. [22, Proposition 2.3.5].
Remark 4 (Robustness in polynomial degree k). In numerical experiments we observed that the LBB
constant cLBB in (12) is also robust in k. To this
R end we computed the LBB constant of one element with
Dirichlet boundary conditions and the condition Ω p dx = 0 on the pressure. The results are shown in Table
1 for varying the polynomial degree k. From the boundedness of the LBB constant on one element and the
h-robustness in (12), the robustness in h and k on the global spaces follows. This will is in the forthcoming
master’s thesis of Philip Lederer.
2.1.6. Discretization of the Brinkman-Stokes problem
With the introduced discretizations of the spatial operators the discrete Brinkman-Stokes problem can
be written as:
Find u ∈ Uh and p ∈ Qh , s.t.
 −1
τ Mh (u, v) + Ah (u, v) + Dh (v, p) = hf, vi ∀ v ∈ Uh ,
(13)
Dh (u, q) =
0 ∀ q ∈ Qh .
Due to coercivity of Ah (respectively τ −1 Mh +Ah ), the LBB-condition of Dh and consistency and continuity
of all bilinear forms Brezzi’s famous theorem (see [21]) can be applied to obtain optimal order a priori error
estimates. For a discussion of the coupling structure of this discretization we refer to Remark 7 below.

7

Remark 5 (Pressure-independence of velocity error). Classical error estimates for mixed problems
result in error estimates which are formulated in the compound norm of the velocity and pressure space.
This has the disadvantage that the discretization error in the velocity depends on the approximation error
in the pressure. As was pointed out in [54], ideally this should not be the case. Due to the fact that discrete
solutions to (13) are exactly divergence-free an error estimate for the velocity field which does not involve
the pressure can be derived, cf. [22, Lemma 2.3.13].
2.2. Projected jumps: An enhancement of the HDG Stokes discretization
The proposed HDG formulation for viscous forces in section 2.1.3 can also be derived as a hybrid mixed
method with a modified flux. This is done in detail in [25] (see also [22, section 1.2.2] or [33]). In that
setting the unknowns for the primal variable (uT ) are approximated with the same polynomial degree k as
the facet unknowns (uF ). Afterwards, in a postprocessing step approximations for the primal unknown, of
one degree higher, k + 1, are reconstructed in an element-by-element fashion. This approach ends up with
a higher order approximation than the previously introduced HDG method considering the use of the same
polynomial degree k for the facet unknowns. This sub-optimality can be overcome by means of a projection
operator which leads to the projected jumps formulation. We first introduce the method and explain how
the method can be implemented afterwards.
2.2.1. Projected jumps: The method.
The idea of the projected jumps formulation is to reduce the polynomial degree of the facet unknowns in
Fh to k − 1,
Y
[P k−1 (F )]d , uF · n = 0, ∀ F ∈ Fh },
(14)
Fh → F h := {uF ∈
F ∈Fh

while keeping the polynomial degree k in Wh . In order to do this in a consistent fashion we have to modify
the bilinearform Ah . First, we introduce the L2 projection Π for a fixed facet F ∈ Fh :
Z
Z
Π : [P k (F )]d → [P k−1 (F )]d ,
(Π f ) v dx =
f v dx ∀ v ∈ [P k−1 (F )]d
(15)
F

F

k−1
T
Due to the fact that ∂u
(F ) on affine linearly mapped elements, the test functions vT and vF in
∂n ∈ P
the first integral on the element boundary in (9) can be replaced by their L2 (F ) projections ΠvT and ΠvF .
If we additionally use Π[[ut ]] instead of [[ut ]] to symmetrize and stabilize the formulation in (9) we end up
with a modified version of the Hybrid Interior Penalty method:
Z
X Z
∂uT
Arh (u, v) :=
ν∇uT : ∇vT dx −
ν
Π[[v t ]] ds
∂n
T ∈T
(16)
Zh T
Z ∂T
∂vT
α
t
−
ν
Π[[u ]] ds +
ν Π[[ut ]]Π[[v t ]] ds, u, v ∈ Wh × F h
∂n
∂T
∂T h

Note that the first two boundary integrals are just reformulated while only the last integral is really changed.
This modification preserves the important properties of Ah : It is still consistent and bounded. Further, a
deeper look into the coercivity proof (see [22]) reveals that coercivity can easily be shown in the modified
(weaker) norm
(
)
2
X
1
∂uT
2
2
t 2
|||u|||1,∗ :=
k∇uT kT + k[[Π u ]]k∂T + h
.
(17)
h
∂n ∂T
T ∈Th

Notice that, as we do not modify Wh , the normal component of uT ∈ Wh is still a polynomial of degree k.
The idea of applying such a reduced stabilization is also discussed and analyzed in [39, 40]. However, only
for the two-dimensional case the realization of the projection Π is discussed (by means of Gauss quadrature,
cf. [39, section 3.4]). In the next section a simple way to implement this projection operator is presented.

8

Remark 6 (Interplay with operator-splitting). If a Hybrid DG formulation for a problem involving
diffusion (viscosity) and convection is applied, such a reduction of the facet unknowns is only appropriate if
diffusion is dominating. One premise of the operator-splitting in this paper is that convection is not involved
in implicit solution steps. Hence, the facet variable will never be involved in the discretization of convection,
s.t. even if the physical problem is convection dominated, the projected jumps modification can be applied
without loss of accuracy.
Remark 7 (Comparison to a fully hybridized formulation). We note that the formulation in (13)
(and the corresponding reduced formulation with Arh (·, ·)) is not a “hybridized” formulation in the usual
sense, see e.g. [25]. Only for the tangential component of the velocity a new unknown is introduced, such
that we only have a “tangentially hybridized” velocity discretization. The pressure is not hybridized. We
comment on the resulting coupling of this HDG method and compare it to a HDG formulation where only
Hybrid unknowns on the facet exist, for the velocity and the pressure. The velocity unknowns in Uh = Wh ×Fh
(or Wh × F h in the reduced case) can be reduced to unknowns on the skeleton by static condensation. The
remaining facet unknowns are the usual unknowns for the normal component of the H(div)-conforming
finite element space and the (tangential) unknowns in Fh (or F h ). Due to the pressure unknowns, static
condensation of the formulation in (13) does not yield a global system which only involves facet unknowns,
but also element unknowns for the pressure. However, except for the constant pressure all pressure unknowns
can be eliminated by static condensation. We summarize: The globally reduced system only involves velocity
unknowns on the facets and a constant pressure on each element. Alternatively one could formulate a
related HDG formulation which only involves facet unknowns for the velocity and the pressure. In such a
formulation the pressure plays the role of the lagrange multiplier for the normal continuity and only for
the (weak) tangential continuity velocity unknowns have to be added on the facet. To obtain an accurate
(and exactly divergencefree) order k approximation of the velocity and an order k − 1 approximation of the
pressure, the facet unknowns would be an order k approximation of the pressure and an order k (or order k−1
if projected jumps or corresponding postprocessing techniques are applied) approximation of the tangential
velocities. Overall the number of unknowns would be smaller compared to the formulation proposed before
by one (constant) pressure unknown per element. However, the price for this is the fact, that the Lagrange
multiplier space (the pressure unknowns on the facets) is much larger and hence the structure of the arising
linear system is different. It very much depends on the linear solver strategies which of the two Hybrid DG
formulations gives the better perfomance.
2.2.2. Projected jumps: Realization.
The following way to implement projected jumps needed for Arh in (16) relies on an L2 -orthogonal basis
for the facet functions in Fh . To obtain an L2 -orthogonal basis, we take a local coordinate system on
each facet (spanned by d − 1 aligned edges) and an L2 -orthogonal Dubiner basis (see [55]) for each vector
component. We consider the three dimensional case, and denote the vectors of the local coordinate system
by e1 and e2 . Translation to the two dimensional case is then obvious.
On each facet F ∈ Fh , the space of (vector-valued) polynomials up to degree k, span(e1 , e2 ) · P k (F ), can
be split into the orthogonal subspaces
Vk−1 := span{e1 , e2 } · P k−1 (F )

and

⊥
Vk−1
:= span{e1 , e2 } · [P k (F ) ∩ P k−1 (F )

⊥

].

(18)

⊥
On each facet F the facet unknown uF can thus be written as uF = ūF +λ with unique ūF ∈ Vk−1 , λ ∈ Vk−1
.
0
We now replace the highest order function λ with two copies λT and λT , each of these functions is associated
with one of the two neighbouring elements T and T 0 , the functions λT are only defined element-local. λT
can thus be eliminated after the computation of the element matrix and finally, only uT ∈ Wh and ūF ∈ F h
appear in the global system.
λT is only responsible for implicitly realising the projection operator. To see this we now consider one
facet F ∈ Fh and one of the neighboring elements T ∈ Th . We use the decomposition corresponding to (18)
for trial and test functions
uF = ūF + λT and vF = v̄F + µT ,
(19)

9

⊥
with ūF , v̄F ∈ Vk−1 and λT , µT ∈ Vk−1
where λT , µT are only supported on T . In (13) we choose the test
⊥
function v such that vT = v̄F = 0, µT ∈ Vk−1
on T and µT 0 = 0 on every other element T 0 . This yields
Z
Z
Z
∂uT
α
α
⊥
ν
µT ds +
.
(20)
ν (utT − ūF − λT )µT ds =
ν (utT − λT )µT ds = 0, ∀ µT ∈ Vk−1
∂n
F
F h
F h

Hence, there holds λT = (I − Π)utT and we have
[[ut ]] = utT − ūF − λT = Π(utT ) − ūF = Π(utT − ūF ) = Π[[ut ]]

(21)

We conclude that it is sufficient to consider the HDG bilinear form Ah as before, where the local element
matrices are computed according to Fh . For each element matrix the degrees of freedom corresponding to
⊥
Vk−1
are then eliminated forming a corresponding Schur complement. This yiels a final stiffness matrix
which is only set up with respect to space unknowns of Wh × F h .
2.3. DG formulation for the convection
For the discretization of the convection part we consider a standard DG finite element space:
Vh := {u : u ∈ [P k (T )]d ∀ T ∈ Th },

NV := dim(Vh ).

(22)

We define the mass bilinear form
MVh

Z
(w, z) :=

w z dx,
Ω

w, z ∈ Vh .

(23)

Using an L2 -orthogonal basis on each element renders the associated mass matrix MV : RNV → RNV
diagonal. A stable spatial discretization of (1) with respect to the convection part is achieved using a
Standard Upwind DG trilinearform. We assume that the convection velocity uT is exactly divergence-free.
Z
X Z
Ch (uT ; w, z) :=
− w ⊗ uT : ∇z dx +
uT ·n ŵ z ds, uT ∈ Wh , w, z ∈ Vh .
(24)
T

T

∂T

where ŵ denotes the upwind value ŵ = limε&0 w(x − εuT (x)). Standard Upwind DG formulations are stable
in the sense that with ∂Ωin := {x ∈ ∂Ω, uT (x) · n < 0} there holds
Z
1
|uT · n| w2 ds + Ch (uT ; w, w) ≥ 0,
∀ w ∈ Vh , u ∈ Uh , div(uT ) = 0.
(25)
2 ∂Ωin
2.4. Transfer operations and embeddings
The discrete convection and Stokes operators have been defined on different spaces. In order to combine
both we introduce transfer operations between the (finite element) spaces to make both discretizations
compatible. We restrict ourselves to two types of transfer operations which are based on embeddings:
I: With Wh ⊂ Vh we have a canonical embedding of Uh in Vh with the embedding operator I : (uT , uF ) ∈
Uh → uT ∈ Vh . Note that the corresponding operation in terms of coefficient vectors, denoted by
I : RNU → RNV , is not an identity.
IT : The embedding operator I implies the canonical embedding I 0 : Vh0 → Uh0 , which maps functionals
I 0 : f ∈ Vh0 → [(uT , uF ) ∈ Uh → f (uT )] ∈ Uh0 . The corresponding matrix representation is IT : RNV →
RNU .
To realize the operator I we consider the equivalent L2 problem for u ∈ Uh .
Z
Z
(Iu) v dx =
uT v dx ∀ v ∈ Vh .
Ω

Ω

10

(26)

In terms of coefficient vectors this reads as
MV (Iu) = MU,V u, ∀ u ∈ RNU

=⇒

with the mixed mass matrix
Z
U,V
V
Mi,j =
ϕW
j ϕi dx, i = 1, . . . , NV , j = 1, . . . , NW
Ω

I = (MV )−1 MU,V

(27)

and MU,V
i,j = 0, j > NW .

Note that MV is diagonal (for affine linear transformations) such that (MV )−1 can be evaluated very
efficiently. The overall cost of the transfer operator I is essentially that of one sparse matrix multiplication.
Let CV : RNU × RNV → RNV denote the discrete convection operator corresponding to the trilinearform
Ch in (25).
With the operator I we can formulate applications of the convection and the mass operations CV (u), MV :
NV
R
→ RNV with respect to functions in the HDG space Uh and denote the corresponding operators by
U
C (u), MU : RNU → RNU ,
CU (u) := IT CV (u)I,

MU := IT MV I.

(28)

Remark 8 (Restriction on time integration scheme). Note that the restriction to these two transfer
operations implies that we do not allow to apply any part of the Stokes operator to a function in Vh and that
no functional on Uh can be used in solution steps involving the convection operator. This is a restriction on
the time integration scheme.
Remark 9 (Curved elements). If curved elements are considered we no longer have Wh ⊂ Vh due to the
Piola transform usually applied to construct H(div)-conforming finite elements. In this case I : Uh → Vh
is not an embedding, but the L2 projection. Nevertheless, MV is block diagonal with blocks which are not
diagonal matrices only on curved elements. Hence, applications of (MV )−1 are still cheap.
2.5. The semidiscrete formulation
With the definitions of the bilinear forms in the previous sections we arrive at
discrete DAE problem: Find u(t) ∈ Uh and p(t) ∈ Qh , such that

∂
hf, vi
∀ v ∈ Uh ,
 MU
h ( ∂t u, v) + Ah (u, v) + Ch (uT ; u, v)+ Dh (v, p) =
Dh (u, q) =
0
∀ q ∈ Qh ,

Mh (u, v) = Mh (u0 , v) ∀ v ∈ Uh ,

the following spatially
t ∈ [0, T ],
t ∈ [0, T ],
t = 0.

(29)

Here, we implicitly used the embedding I to define Ch (uT ; ·, ·) on Uh . Due to (25) we have with
(

∂
1 d
d
u, u)L2 =
kuk2L2 = kukL2 kukL2 = hf, ui − Ah (u, u) − Ch (uT ; u, u) − Dh (u, p) ≤ kf kL2 kukL2 (30)
| {z } |
{z
} | {z }
∂t
2 dt
dt
≥0

≥0

=0

the stability of the kinetic energy:

d
kukL2 ≤ kf (t)kL2
dt
In the next section we discuss operator splitting time integration methods to solve (29) efficiently.

(31)

3. Operator-splitting time integration
In this section we are faced with the problem of solving the semi-discrete Navier-Stokes problem with a
proper time integration scheme. For ease of presentation we neglect external forces (f = 0) in the following

11

and consider the problem in terms of discrete operators corresponding to the bilinear (trilinear) forms
introduced in the previous section: Find u(t) ∈ RNU and p(t) ∈ RNQ , such that

U
in [0, T ],
 MU ∂u
∂t + Au + C (u) u +∆t D p = f
(32)
DT u = 0 in [0, T ],

u(t = 0) = u0 .
In the time integration scheme we want to explicitly exploit the properties of the spatial discretization,
i.e. the convection operator C should only be involved explicitly in terms of operator evaluation. Due
to the DAE-structure we require time integration schemes which are stiffly accurate. Hence, the solution
of a Stokes-Brinkman problem should conclude every time step so that the incompressibility constraint is
ensured. For operator splittings of this kind different approaches exist. We briefly discuss three approaches.
In section 3.1 we discuss additive decomposition methods, like the famous class of IMplicit EXplicit (IMEX)
schemes. Product decomposition methods, sometimes also called exponential factor splittings, are discussed
in section 3.2 in the framework of operator-integration-factor splittings introduced in [48]. We discuss the
advantages and disadvantages of both approaches and motivate the consideration of a different approach.
An operator-splitting modification of the famous fractional step method, cf. [56], which eliminates the
most important disadvantages of the additive and multiplicative decomposition methods is then introduced
and discussed in section 3.3. We want to stress that the considered operator splitting approaches are of
convection-diffusion type and should not be confused with projection methods like the Chorin splitting [57].
3.1. Additive decomposition using IMEX schemes
Additive decomposition methods distinguish spatial operators that are treated implicitly and those that
are treated explicitly. The decomposition is additive in the sense that every solution (sub-)step involves
both operators. This is used by IMplicit EXplicit (IMEX) schemes (see [44, 45, 46]). The simplest of these
schemes is the semi-implicit Euler method for which one time step of size ∆t reads as:

(MU + ∆tA)un+1 + ∆tDpn+1 = MU un − ∆tCU (un )un
(33)
DT un+1 = 0
Convection only appears on the r.h.s. so that only operator evaluations occur for the convection operator,
while the remainder appears fully implicit. This scheme is obviously only first order accurate and only conditionally stable. For higher order methods of this decomposition type (e.g. using Multistep or partitioned
Runge-Kutta schemes) we refer to [44, 45, 46].
The major drawback of this type of operator splitting methods is the fact that the time step size of the
explicit and the implicit part of the decomposition have to coincide. Thereby the stability restriction caused
by the explicit treatment of the convection part dictates not only the number of explicit evaluations but also
- which is typically more expensive - the number of solution steps for the implicit part. The decomposition
methods considered in the sections 3.2 and 3.3 overcome this issue.
3.2. Product decomposition with operator-integration-factor splitting
The disadvantage of IMEX methods can be avoided with product decomposition methods, where sequences of separated problems are solved successively. The separated problems then only involve the Stokes
or the convection operator at the same time. The major benefit of this is the fact, that the numerical solution of the sub-problems (e.g. the size of the time steps) can be chosen completely different. The derivation
of the so called operator-integration-factor splitting approach has been derived in [48, Section 2.1]. We use
this idea to obtain an operator decoupling for the DAE which allows to formally rewrite (32) as
 ∂ t→t∗
∗
u) + Qt→t M−1 (Au + Dp) = 0 in [0, T ],
∂t (Q
(34)
DT u = 0 in [0, T ],
∗

for some arbitrary t∗ ∈ R with the integration factor Qt→t , the propagation operator, specified later in
section 3.2.2.
12

Be aware that M−1 is to be understood only formally for now. We are able to apply any suitable (i.e.
implicit and stiffly accurate) time integration method for (34). We do this for a first order method here to
explain the procedure and refer to the literature [48] for higher order variants.
3.2.1. Implicit Euler.
We employ the implicit Euler method on (34) and arrive at
 1
n
∗
n+1
∗
tn+1 →t∗ n+1
u
− Qt →t un ) + Qt →t M−1 (Aun+1 + Dpn+1 ) = 0,
∆t (Q
DT un+1 = 0.
After setting t∗ = tn+1 and multiplication with M this simplifies to

n
n+1
(M + ∆tA) un+1 + ∆tDpn+1 ) = MV,U Qt →t Iun ,
T n+1
D u
= 0,

(35)

(36)

which is the solution of a Stokes-Brinkman problem as in (33), but with a different right hand side.
3.2.2. The propagation operator.
1
2
1
2
The propagation operator Qt →t is defined as Qt →t w = v(t2 ) where v(s) solves
∂v
+ (MV )−1 CV (s)v(s) = 0,
∂s

∀ s ∈ (t1 , t2 ],

v(t1 ) = w.

(37)

Here CV (s) := CV (ú(s)) with ú(s) an extrapolation of divergence-free solutions of previous time steps.
Note, that the extrapolation of convection velocities renders the problem (37) linear hyperbolic and ensures
stability in the sense of (25). After replacing Q by a numerical time integrator Q∆t for (37), the time
integration method is completely specified. The order of accuracy of the extrapolation ú and the time
integrator Q∆t should coincide with the order of the time integration method applied on (34).
3.2.3. Properties of product decompositions.
At this point the advantage of product decomposition methods over IMEX schemes is evident: As the
time integration of (37) is independent of the one for (34) the (not necessarily) explicit time integrator can
deal with stability restrictions (typically using multiple time steps) without influencing the time step for
(34). This separation of the two problems also introduces the biggest disadvantage of product decomposition
methods: an additional consistency error. Even if the DAE (32) has a stable stationary solution, product
decomposition methods may not reach it. This is not the case for monolithic or additive decomposition
approaches. Nevertheless, this splitting error is controlled by the time discretization error.
Remark 10 ((Marchuk-)Yanenko splitting). If the implicit Euler discretization as in (36) is combined
with an Euler method for (37), the famous Yanenko splitting method is recovered.
3.3. A modified fractional-step-θ-scheme
In this section we introduce an approach to circumvent severe CFL-restrictions and splitting errors at
the same time. We no longer ask for sub-problems that only involve the Stokes or the convection operator
(as in section 3.2), but ask for sub-problems which only involve one of both implicitly. The method is based
on [49] where the well-known fractional-step-θ-scheme, cf. [56, Chapter II, section 10], is modified. The
resulting method is an additive decomposition method without the time step restrictions of IMEX schemes.
We start with formally writing down an operator-splitting version of the fractional-step-θ-scheme. The
scheme is divided into three steps, the first and the last step treat the Stokes part implicitly and the

13

convection part explicitly as in (33) while the second step treats convection implicitly and viscosity forces
explicitly.

(M + θ∆tA)un1 + θ∆t D pn1 = Mun0 − θ∆tCun0
Step 1(tn0 → tn1 ) :
(38a)
DT un1 = 0
= Mun1 − θ∗ ∆t(Aun1 + Dpn1 )
Step 2 (tn1 → tn2 ) : { (M + θ∗ ∆tC)un2

(M + θ∆tA)un3 + θ∆t D pn3 = Mun2 − θ∆tCun2
Step 3 (tn2 → tn3 ) :
DT un3 = 0

(38b)
(38c)

The time stages are labeled by the superscripts n0 , n1 , n2 , n3 , respectively, where n0 denotes initial data
∗
and n3 the final time stage. The time steps size
√ for the first and the last step is θ∆t and θ ∆t in the
∗
middle step with θ = 1 − 2θ. Here θ = 1 − 1/ 2. We note that specifications for M and C with respect
to the considered spaces and the convection velocity for C are still missing at this points. The scheme is
second order accurate. In contrast to the unsplit fractional-step-θ-scheme the stability analysis of this time
integration scheme is an open problem. In our experience, however, the stability restrictions are much less
restrictive than those of a comparable IMEX schemes.
3.3.1. Sub-steps in different spaces.
Initially, we stated that we want to avoid solving linear systems involving convection, for efficiency
reasons. This seems to be contradictory to what is formulated in (38b). Nevertheless, as only convection is
involved implicitly an efficient numerical solution is still possible. We apply a simple iterative scheme which
only involves explicit operator evaluations of the convection to do so. We explain this in section 3.3.2. To
do this efficiently, we want Step 2 to be formulated in the space Vh while Step 1 and 3 are to be formulated
in the space (Uh , Qh ). This poses problems the solution of which we discuss in this section.
In Step 1 and Step 3 the adjustments are obvious: Step 1 only depends on initial data in Uh . The
initial data for Step 3 is in Vh but appears only in terms of functionals (Mv and Cv) such that the transfer
operations are clear. Step 2 is more involved, cf. remark 8. Functionals in Uh0 (such as Au) are in general
not functionals in Vh0 . We use the first equation in (38a) to formally define a different representation of the
functionals required in Step 2:
θ∆t gn1 = −θ∆t(Aun1 + Dpn1 ) = M(un1 − un0 ) + θ∆tCun0

(39)

Now we replace M and C with operations suitable for a setting in Vh :
RNV 3 gn1 :=

1
MU,V (un1 − un0 ) + θ∆tCV (un0 )Iun0
θ∆t

We arrive at the modified fractional-step-θ-scheme:

(MU+ θ∆tA)un1 + θ∆tD pn1
Step 1 :
DT un1

Step 2 : (MV + θ∗ ∆tCV (tn2 )vn2

(MU+ θ∆tA)un3 + θ∆tD pn3
Step 3 :
DT un3

(40)

= MU un0 − θ∆t IT CV (un0 ) I un0
=0

(41a)

= MU,Vun1 − θ∗ ∆tgn1

(41b)

V,U

=M
=0

n2

T

V

n2

v − θ∆t I C (t )v

n2

(41c)

where we replaced the generic operators M, C with suitable ones. We recall the definition of the extrapolated
convection operator CV (tn2 ) := CV (ú(tn2 )) the convection velocity of which is (linearly) extrapolated from
exactly divergence-free velocities.
3.3.2. Iterative solution of the implicit convection problem.
In (41b) we need to solve a problem of the form
v + τ ∗ (MV )−1 CV v = g∗
14

(42)

for given τ ∗ , g∗ and constant convection CV . We do this by means of a pseudo time-stepping method, i.e.
we formulate (42) as the stationary solution to
∂
v(s) + v(s) + τ ∗ CV v(s) = g∗ ,
∂s

s ∈ [0, ∞).

(43)

Note that a stationary solution exists as Id + τ ∗ CV only has eigenvalues λ, with Re(λ) > 1. This stationary
solution is approximated with a few explicit Euler time step with an artificial time step size ∆s.
vi+1 = vi + ∆s(g∗ − vi + τ ∗ (MV )−1 CV vi ).

(44)

This procedure can be interpreted as a Richardson iteration applied to (42). With a time step size which is
tailored for stability (as in the numerical solution to (37)) we iterate (44) until the initial residual is reduced
by a prescribed factor. To solve the convection step, Step 2, we only require operator evaluations for the
convection. The time step size ∆s in this iteration is decoupled from the time step size ∆t used of the
overall scheme.
4. Numerical examples
In this section we consider different test problems which essentially purpose three different goals:
1. The validation of convergence properties of the HDG method with and without the projected jumps.
2. The investigation and quantification of the dependency of the sparsity pattern of arising linear systems
on the choice of the discrete velocity spaces.
3. The evaluation of the performance of the space and time discretization on benchmark problem.
The examples in section 4.1 and section 4.2 aim at goal 1 and are two-dimensional test cases with known
exact solutions for a stationary Stokes and Navier-Stokes problem, respectively. These cases allow for a
thorough investigation of the convergence history of the proposed spatial discretizations in the usual norms.
The test case in section 4.3 aims at goal 2 and is concerned with the impact of the choice of discretization
spaces on the complexity of linear systems. For this purpose, a three-dimensional vector-valued Poisson
problem is considered using different H(div)-conforming DG and HDG spaces. The impact of hybridization
and the improvement due to the projected jumps modification is compared to other DG methods. Finally,
we approach goal 3 by considering two challenging transient benchmark problems in sections 4.4 and 4.5.
The benchmark problems have been defined within the DFG Priority Research Programm ’Flow Simulation
on High Performance Computers’ and are formulated in [58]. In section 4.4 a two-dimensional benchmark
problem is considered to demonstrate the accuracy of our method for a demanding test case and to compare
the discussed time integration methods. Finally, in section 4.5 we discuss a three-dimensional, and hence
computationally demanding, benchmark problem from [58]. We compare accuracy and run-time performance
to the data of the studies in [58, 59, 60].
The methods discussed in this paper have been implemented in the add-on package ngsflow [61] for the
high order finite element library NGSolve [62]. The computations in this section have also been carried out
with this software. Throughout this section we only consider direct solvers and comment on linear solvers
below, in remark 12.
For the computations of the benchmark problem in sections 4.4 and 4.5 we used the reduced H(div)conforming space Whred , cf. remark 1 and the projected jumps modification, cf. section 2.2.
4.1. Stokes flow around obstacle
We consider the Stokes problem in the domain Ω = [−2, 2]2 \ Ω− with the circular obstacle Ω− :=
{kxk ≤ 1} and viscosity ν = 1. On the whole boundary except for the outflow boundary {x = 2} we
prescribe Dirichlet boundary conditions, on {x = 2} we prescribe Neumann-kind boundary conditions such
1
that the solution to the problem is given by u = (∂y Ψ, −∂x Ψ) with the potential Ψ = y(1 − x2 +y
2 ). Note
−
that Ψ|∂Ω− is constant so that u · n = ∇Ψ × n = 0. On ∂Ω we have u · n = 0, but u 6= 0, i.e. we have a
slip on the obstacle. Further we have that p = 0 in Ω.
15

k∇(u − uh )kL2 (Ω)

ku − uh kL2 (Ω)

kp − ph kL2 (Ω)

100

10−1

10−2

10−2

10−4

10−4

10−6
10

10−3

k=1
k=2
k=3
k=4
k=5

−8

10−10
1

2

10−6
O(h6 )

3
refinement level

10−8
4

5

10−5

k=1
k=2
k=3
k=4
k=5
1

2

k=1
k=2
k=3
k=4
k=5

10−7
O(h5 )

10
3

refinement level

4

5

−9

1

2

O(h5 )

3

4

5

refinement level

Figure 2: Convergence of divergence-conforming HDG method with (solid) and without (dotted) projected jumps in different
norms for the example in section 4.1. The dotted lines in gray correspond to the HDG method without projected jumps, the
solid lines in color correspond to the HDG method with projected jumps. Except for the error in the pressure the results are
hardly distinguishable.

On an initially coarse unstructured grid with only 72 triangles we obtain the results displayed in Figure
2 after succesive uniform mesh refinements. We compare the exact solution u, p and the computed solution
uh , ph and investigate the convergence of the error in the norms ku − uh kL2 (Ω) , k∇(u − uh )kL2 (Ω) and
kp − ph kL2 (Ω) . We consider the discretization with and without the projected jumps modification. In all
norms we observe optimal order convergence, i.e. ku − uh kL2 (Ω) = O(hk+1 ), k∇(u − uh )kL2 (Ω) = O(hk ) and
kp − ph kL2 (Ω) = O(hk ) where k is the order of the velocity field inside the elements. Note that after static
condensation the remaining degrees of freedoms are the unknowns corresponding to order k polynomials
for the tangential and order k polynomials for the normal component on the facets and one constant per
element for the pressure. For the projected jumps formulation the tangential unknowns are reduced by one
order. The difference between the error of both formulations — with and without projected jumps — is only
marginal in the velocity. In the pressure field we observe a difference between the methods, however both
converge with optimal order and the difference decreases for increasing polynomial degree k. For the impact
of the projected jumps modification concerning the computational effort, we refer to section 4.3 where this
aspect is discussed in more detail.
4.2. Kovasznay flow
As a test case for the stationary Navier-Stokes equations, we consider the famous example by [63]. On
the boundary of the domain Ω = [− 12 , 32 ] × [0, 2] we again prescribe inhomogeneous Dirichlet data, so that
the exact solution to the Navier-Stokes equations (with ν = 1) is


1
1 − eλx cos(2πy)
u(x, y) =
, p(x, y) = − e2λx + p̄ with p̄ ∈ R
λ λx
sin(2πy)
2
2π e
R
2
Here, λ = ν −1 +√−8π
and p̄ so that Ω p dx = 0.
ν −2 +64π 2
To approximate the stationary solution we initially solve the Stokes problem corresponding to the boundary data and use the simple IMEX scheme of first order, cf. (33), to progress to t = 1000 with 50 time steps
of size ∆t = 20. This choice of time discretization parameters gives stable solutions for all considered spatial
discretizations but at the same time provides time discretization errors which are neglegible compared to
the spatial errors. On an initally coarse unstructured grid with only 18 triangles we obtain the results
displayed in Figure 3 after succesive uniform mesh refinements. We observe that the errors converge optimal
in the considered norms for the pressure and the velocity. Further, we observe that the difference between
the results obtained with and without the projected jumps formulation, e.g. with and without a reduction
16

k∇(u − uh )kL2 (Ω)

ku − uh kL2 (Ω)

kp − ph kL2 (Ω)

100

101

101

10−2

10−1

10−1

10−4
10

10−3

k=1
k=2
k=3
k=4
k=5
k=6

−6

10−8
10−10
1

2

k=1
k=2
k=3
k=4
k=5
k=6

10−5
O(h7 )

10−7
3

refinement level

4

5

1

2

k=1
k=2
k=3
k=4
k=5
k=6

10−3
10−5

O(h6 )

10−7
3

4

5

refinement level

1

2

O(h6 )

3

4

5

refinement level

Figure 3: Convergence of divergence-conforming HDG method with (solid) and without (dotted) projected jumps in different
norms for the example in section 4.2. The results are almost identical.

of the degrees of freedoms at the facets, is again only marginal. In fact one only observes a (very small)
difference in the case k = 1.
4.3. Linear systems - A comparison between DG and HDG methods
We consider the comparably simple vector-valued Poisson problem
− ∆u = f in Ω,

u = 0 on ∂Ω,

(45)

discretized with four different methods. The three-dimensional domain Ω is the same as the one in section 4.5
with an unstructured tetrahedral mesh consisting of 3487 elements. The problem (45) leads to discretizations
with symmetric positive definite matrices A. Only the lower triangular part of the sparse matrix has to be
stored and a sparse Cholesky factorization algorithm can be used. In this comparison, we used the sparse
direct solver PARDISO, cf. [64, 65]. We are only concerned with the sparsity pattern of different methods
before and after static condensation and do not compare accuracy or conditioning of the discretizations.
Moreover, for the sake of simplicity we do not apply the reduction of the space Wh as discussed in remark 1.
The following quantities of interest for the linear systems arising from discretizations of (45) are displayed in
#dof[K] : number of unknowns (in thousands)
#cdof[K] : number of unknowns after static condensation (in thousands)
Table 2:
#nzeA[K] : number of non-zero entries in the system matrix A (lower triangle) (in thousands)
#nzeL[K] : number of non-zero entries in the Cholesky factor L (in thousands)
Four different methods are considered on Wh or Uh with varying polynomial degree between k = 1 and
k = 6:
1. HDG: The HDG method proposed in section 2 without projected jumps.
2. PHDG: The HDG method with projected jumps, cf. section 2.2.
3. Std.DG.: A standard DG method using the space Wh where the basis is constructed such that all
degrees of freedom from one element couple with all degrees of freedom from adjacent elements.
4. N.DG: A nodal DG method using the space Wh where basis functions are assumed to be constructed
such that degrees of freedom associated to one element couple only with degrees of freedom from
adjacent elements which have support on the shared facet. At the same time we assume that basis
functions are constructed such that the number of basis functions with support on a facets is minimized,
cf. [11]. In terms of the sparsity pattern this nodal DG method represents the best case for a DG
method without hybridization.
17

Std.DG

N.DG

HDG PHDG

Std.DG

k=1
#dof[K]
#cdof[K]
#nzeA[K]
#nzeL[K]

23
23
732
10 113

23
23
676
10 208

69
69
2 037
17 768

146
146
21 686
261 977

146
146
16 051
260 416

#dof[K]
#cdof[K]
#nzeA[K]
#nzeL[K]

453
453
184 035
2 099 690

453
411
101 473
1 741 072

298
229
22 443
194 524

38
38
637
5 569

67
67
5 073
64 463

PHDG

158
137
8 103
70 138

112
91
3 628
31 415

500
343
50 412
435 129

423
267
30 588
264 183

67
67
4 177
69 831
k=4

237
168
12 124
104 496

271
271
69 525
814 168

271
261
46 912
731 253

681
389
64 816
557 798

702
702
424 764
4 752 072

702
597
200 813
3 519 061

k=5
773
480
98 696
847 913

HDG

k=2

k=3
#dof[K]
#cdof[K]
#nzeA[K]
#nzeL[K]

N.DG

k=6
1 128
640
175 321
1 502 558

1 022
533
121 944
1 045 444

Table 2: Comparison of different DG methods for the vector-valued reaction Poisson problem (45).

In Table 2 we observe that for small polynomial degree k the amount of additional unknowns required for
the HDG formulation is quite large as are the nonzero entries in the system matrix. Nevertheless except
for k = 1, the number of nonzero entries in the Cholesky factor are comparable (k = 2) to the Standard
DG methods or less (k > 2). For high order, i.e. k ≥ 4 the HDG method performs significantly better as it
has less nonzero entries in A and L. The HDG space with the projected jumps modification improves the
situation dramatically. It essentially compensates the overhead of the HDG method for small k. But even
for k = 6 the effect is still significant. We note, that the difference between the first three methods increases
for increasing polynomial degree k whereas the difference between the last two method, PHDG and HDG,
decreases. In all cases the HDG method with projected jumps outperforms all alternatives.
4.4. A two-dimensional benchmark problem
In this section we consider the benchmark problem denoted as “2D-2Z” in [58] where a laminar flow
around a circle-shaped obstacle is considered. The Reynolds number is moderately high (Re = 100) and
results in a periodic vortex street behind the obstacle. We briefly introduce the problem (for more details
we refer to [58]) and the numerical setup to investigate spatial and temporal discretization errors. Finally,
we discuss the obtained results.
4.4.1. Geometrical setup and boundary conditions.
The domain is a rectangular channel without an almost vertically centered circular obstacle, cf. Figure
4,
Ω := [0, 2.2]×[0, 0.41] \ {kx − (0.2, 0.2)k2 ≤ 0.05}.
(46)

The boundary is decomposed into Γin := {x = 0}, the inflow boundary, Γout := {x = 2.2}, the outflow
boundary and ΓW := ∂Ω\(Γin ∪Γout ), the wall boundary. On Γout we prescribe natural boundary conditions
(−ν∇u + pI) · n = 0, on ΓW homogeneous Dirichlet boundary conditions for the velocity and on Γin the
inflow Dirichlet boundary conditions
u(0, y, t) = uD = (3/2 · ū) · 4 · y(dy − y)/d2y · (1, 0, 0).

Here, ū = 1 and the viscosity is fixed to ν = 10−3 which results in a Reynolds number Re = 100.
4.4.2. Drag and Lift.
The quantities of interest in this example are the (maximal and minimal) drag and lift forces cD , cL that
act on the disc. These are defined as


Z 
Z 
1
∂u
1
∂u
cD := 2
ν
− pn · ex ds,
cL := 2
ν
− pn · ey ds.
ū r Γ◦
∂n
ū r Γ◦
∂n
18

Here ex , ey denote the unit vectors in x and y direction, r = 0.05 is the radius of the obstacle, ū is the
average inflow velocity (ū = 1) and Γ◦ denotes the surface of the obstacle.
4.4.3. Numerical setup.
We use an unstructured triangular grid with an additional layer of quadrilaterals around the disk which
is anisotropically refined towards the disk once. In Figure 4 the geometry, the mesh and a typical solution
is depicted.
2

0
Figure 4: Sketch of the mesh and the solution (color coding corresponding to velocity magnitude kuk2 ) to the problem considered
in section 4.4 at a fixed time t (left) and zoom-in on the boundary layer mesh (right).

In order to be able to neglect time discretization errors, when investigating the spatial accuracy, we
consider the use of a well-known stiffly accurate second order Runge-Kutta-IMEX scheme, taken from [45,
section 2.6], with an extremely small time step size 3.125 · 10−5 which means that roughly 10000 time steps
are used to resolve one full period. The duration of one full cycle is roughly 1/3s. We also fix the mesh and
consider a pure p-refinement, i.e. variations of the polynomial degree k.
For the investigations of the temporal discretization error, we consider a fixed polynomial degree k = 6,
with the same mesh. For the considered example the stability restriction of the second order IMEX scheme is
severe. A time step of below 10−3 has to be considered. The intention of the discussion of operator-splitting
methods in section 3 has been to present strategies to circumvent or relax these severe conditions. We
compare the performance of a product decomposition method, a second order operator-integration-factor
splitting version of the BDF2 method, and the modified fractional-step-θ-scheme discussed in section 3.3.
Note that these methods allows to consider much larger time steps than the IMEX scheme. As references
we give the values from the literature [66] and from the IMEX scheme with ∆t = 10−3 (stability limit) and
the reference solution with ∆t = 3.125 · 10−5 .
4.4.4. Numerical results: spatial discretization.
In Table 3 the quantities of interest are shown for varying polynomial degree k. As a reference we also
show the result obtained by FEATFLOW [67] with a discretization using quadrilateral meshes and continuous
second order finite elements for the velocity with a discontinuous piecewise linear pressure (Q2 /P1disc ). These
results have been made accessible on [66]. We observe a rapid convergence for the p-refinement, i.e. the
#dof
k
k
k
k
k
k

=1
=2
=3
=4
=5
=6

ref. [66]

max cD

min cD

max cL

min cL

211
148
558
441
797
626

2.52594
3.22841
3.23184
3.22714
3.22759
3.22757

2.47871
3.16260
3.16842
3.16401
3.16432
3.16430

0.65728
1.00571
0.98822
0.98431
0.98578
0.98580

-0.81672
-1.03894
-1.02427
-1.01906
-1.02053
-1.02053

167 232
667 264

3.22662
3.22711

3.16351
3.16426

0.98620
0.98658

-1.02093
-1.02129

2
4
6
9
12
16

Table 3: Accuracy of the spatial discretization: results for different polynomial degrees.

increase of the polynomial degree. Compared to the results from the literature [66] the same order of
accuracy is achieved with a lot less degrees of freedoms.
19

4.4.5. Numerical results: temporal discretization.
The results for the temporal discretization are shown in Table 4. First of all, we observe that the second
order IMEX scheme is already very accurate at its stability limit. But the method does not allow to choose
larger time steps. This is in contrast to the alternative methods discussed here. For these methods we
can consider much large time steps and observe a second order convergence. In this example the product
decomposition method is more accurate than the modified fractional step method by one (time) level. We
note that one major concern with this method is however, that splitting errors appear also if a stationary to
the flow problem exists. This is not the case for additive decomposition methods or the proposed modified
fractional step method.

2nd order IMEX

operator-integrationfactor splitting
(BDF2)
modified fractionalstep-θ-scheme

1/∆t

max cD

min cD

max cL

min cL

<1 000
1 000
32 000

3.22754
3.22757

unstable
3.16437 0.98486
3.16431 0.98580

-1.01957
-1.02053

125
250
500
1 000

3.38456
3.25564
3.23239
3.22819

3.28279
3.18643
3.16836
3.16394

1.16052
1.01725
0.98864
0.98320

-1.19939
-1.05276
-1.02378
-1.02013

125
250
500
1 000

3.38272
3.28713
3.25036
3.23656

3.29673
3.21483
3.18127
3.17094

1.07667
1.00937
0.99172
0.98721

-1.12227
-1.04840
-1.02807
-1.02227

Table 4: Accuracy of operator-splitting time integration methods

4.5. A three-dimensional benchmark problem
Finally, we consider a three dimensional benchmark problem, the problem “3D-3Z” in [58]. In contrast
to the problem in section 4.4 the inflow velocity is varied over time and the observed time interval is fixed.
The maximal Reynolds number in this configuration is also Re = 100 as in the previous section. The focus
in this section is on the study of performance in the sense of computational effort over accuracy.
4.5.1. Geometrical setup and boundary conditions.
The geometrical setup is a generalization of the problem in the previous section. A cylindrical-shaped
obstacle is places in a cuboid-shaped channel slightly above the vertical center:
Ω := [0, 2.5]×[0, 0.41] ×[0, 0.41] \ {k(x1 , x2 ) − (0.5, 0.2)k2 ≤ 0.05}.

(47)

Boundary conditions are chosen as in the previous example except for a change to unsteady inflow boundary
conditions:
u(0, y, z, t) = uD (t) = (9/4 · ū(t)) · 16 · y(dy −y)/d2y · z(dz −z)/d2z · (1, 0, 0),
with ū(t) the average inflow velocity which is time-dependent, ū(t) = sin(πt/8). The considered time interval
is [0, 8s], the viscosity is set to ν = 10−3 s.t. the Reynolds number varies between Re = 0 and Re = 100.
4.5.2. Drag and Lift.
Again, the quantities of interest in this example are the (maximal and minimal) drag and lift forces cD ,
cL that act on the disc. These are defined as


Z 
Z 
∂u
1
∂u
1
ν
− pn · ex ds
cL := 2
ν
− pn · ey ds
cD := 2
ūmax rh Γ◦
∂n
ūmax rh Γ◦
∂n
with ūmax = maxt∈[0,8] ū(t) = 1, r = 0.05, h = 0.41 and Γ◦ the surface of the obstacle.
20

Figure 5: Used mesh (left) and solution at t = 0.4 (right) for the benchmark problem in section 4.5.

4.5.3. Numerical setup.
We use an unstructured tetrahedral mesh consisting of 5922 elements, cf. Figure 5. For the time
discretization we use the modified fractional-step-θ-scheme discussed in section 3.3. To compensate for the
increase in the spatial accuracy for increasing polynomial degree k, we adapt the number of time steps
accordingly. The computations were carried out on a shared-memory computer with 24 cores. We comment
on details of the computing times in remark 11.

k
k
k
k

=2
=3
=4
=5

#ndof [K]

max cD

max cL

min cL

169
343
595
939

3.43046
3.29331
3.29853
3.29798

0.00262
0.00277
0.00278
0.00278

-0.016289
-0.011099
-0.010762
-0.011054

0.0080
0.0080
0.0040
0.0040

0.0028
0.0028

-0.010992
-0.010999

0.01
0.005

ref. [59]

11 432
89 760

3.2963
3.2978

ref. [60]

7 036

3.2968

ref. [58]

[3.2,3.3]

∆t [s]

comp. time [s]
492
964
3 087
6 670

×
×
×
×

24
24
24
24

35 550 × 24
214 473 × 48

-0.011
[0.002,0.004]

Table 5: Numerical results for the benchmark problem “3D3Z” in [58]. Results obtained with different polynomial degrees and
reference values.

4.5.4. Numerical results.
In Table 5 the results obtained are compared with the literature in terms of accuracy and computing
time. We observe that we can achieve the same level of accuracy as the results in [59] (and [60]) with a
computing time which is dramatically smaller. We note that the computing time per degree of freedom is
actually worse than in [59]. Nevertheless, the same accuracy is achieved with much less degrees of freedoms
using the high order method (k > 2), s.t. our computations exceed the performance results in the literature.
In the study [59] one of the conclusions is that their third order method (Q2 /P1disc ) is much more efficient
compared to lower order methods. We extend this conclusion in the sense that the use of even higher order
discretizations, i.e. k > 2, increases efficiency even further. Moreover, high order discretizations can be
implemented efficiently. One important component for the efficient handling of the Navier-Stokes equations
with our high order discretization is the time integration using operator-splitting.
Remark 11 (Computation times). We remark on the computation for the case k = 5. In that computation approximately 65% of the computing time has been spend on the solution of linear systems for Stokes-type
problems, 30% on convection operator evaluations and 5% on the setup and remaining operations. To solve
the implicit convection problem (Step 2 in (41)) an average of 20 iterations has been applied.
Remark 12 (Linear systems). In the test cases in this section we only applied direct solvers which is
possible due to the (comparably) small size of the arising linear systems. For problems with increasing
21

complexity efficient linear solvers are mandatory. The development of suitable preconditioners of the Stokes
problem is based on efficient preconditioning of the bilinear form Ah . For the scalar problem we could
show poly-logarithmic bounds in k for the condition number of standard p-version domain decomposition
preconditioners, cf. [52]. We plan to investigate suitable generalizations of this preconditioner for the Stokes
problem in the future.
5. Conclusion
We presented and discussed a combined DG/HDG discretization tailored for efficiency. We summarize
the core components. We split the Navier-Stokes problem into linear Stokes-type problems and hyperbolic
transport problems by means of operator-splitting time integration. For the Stokes-type problem we use an
H(div)-conforming Hybrid DG formulation with a new modification: the projected jumps formulation. The
Hybrid DG formulation facilitates the efficient solution of linear systems compared to other DG methods.
The projected jumps formulation improves its efficiency even more. For the hyperbolic transport problem we
apply a standard DG formulation. In numerical test cases we demonstrated the performance of the method.
Acknowledgements
The authors greatly appreciate the contribution of Philip Lederer related to the calculation of LBBconstants, cf. remark 4.
References
[1] Reed WH, Hill T. Triangular mesh methods for the neutron transport equation. Los Alamos Report LA-UR-73-479 1973;
.
[2] Lesaint P, Raviart PA. On a finite element method for solving the neutron transport equation. Mathematical aspects of
finite elements in partial differential equations 1974; 33:89–123.
[3] Johnson C, Pitkäranta J. An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation. Mathematics
of computation 1986; 46(173):1–26.
[4] Bassi F, Rebay S. High-order accurate discontinuous finite element solution of the 2D Euler equations. Journal of computational physics 1997; 138(2):251–285.
[5] Bassi F, Rebay S, Mariotti G, Pedinotti S, Savini M. A high-order accurate discontinuous finite element method for
inviscid and viscous turbomachinery flows. Proceedings of 2nd European Conference on Turbomachinery, Fluid Dynamics
and Thermodynamics, Technologisch Instituut, Antwerpen, Belgium, 1997; 99–108.
[6] Cockburn B, Shu CW. The Runge–Kutta discontinuous Galerkin method for conservation laws V: multidimensional
systems. Journal of Computational Physics 1998; 141(2):199–224.
[7] Arnold DN, Brezzi F, Cockburn B, Marini LD. Unified analysis of discontinuous Galerkin methods for elliptic problems.
SIAM journal on numerical analysis 2002; 39(5):1749–1779.
[8] Houston P, Schwab C, Süli E. Discontinuous hp-finite element methods for advection-diffusion-reaction problems. SIAM
Journal on Numerical Analysis 2002; 39(6):2133–2163.
[9] Rivière B. Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation. Society
for Industrial and Applied Mathematics, 2008.
[10] Di Pietro DA, Ern A. Mathematical Aspects of Discontinuous Galerkin Methods, vol. 69. Springer Science & Business
Media, 2012.
[11] Hesthaven JS, Warburton T. Nodal discontinuous Galerkin methods: algorithms, analysis, and applications. Springer
Science & Business Media, 2007.
[12] Karniadakis G, Sherwin S. Spectral/hp element methods for computational fluid dynamics. Oxford University Press, 2013.
[13] Girault V, Raviart PA. Finite element methods for Navier-Stokes equations: theory and algorithms, vol. 5. Springer
Science & Business Media, 2012.
[14] Elman HC, Silvester DJ, Wathen AJ. Finite elements and fast iterative solvers: with applications in incompressible fluid
dynamics. Oxford University Press, 2014.
[15] Donea J, Huerta A. Finite element methods for flow problems. John Wiley & Sons, 2003.
[16] Toselli A. hp discontinuous Galerkin approximations for the Stokes problem. Mathematical Models and Methods in Applied
Sciences 2002; 12(11):1565–1597.
[17] Schötzau D, Schwab C, Toselli A. Mixed hp-DGFEM for incompressible flows. SIAM Journal on Numerical Analysis 2002;
40(6):2171–2194.
[18] Girault V, Rivière B, Wheeler M. A discontinuous Galerkin method with nonoverlapping domain decomposition for the
Stokes and Navier-Stokes problems. Mathematics of Computation 2005; 74(249):53–84.

22

[19] Cockburn B, Kanschat G, Schötzau D. A locally conservative LDG method for the incompressible Navier-Stokes equations.
Mathematics of Computation 2005; 74(251):1067–1095.
[20] Cockburn B, Kanschat G, Schötzau D. A note on discontinuous Galerkin divergence-free solutions of the Navier–Stokes
equations. Journal of Scientific Computing 2007; 31(1-2):61–73.
[21] Brezzi F, Fortin M. Mixed and hybrid finite element methods, vol. 15. Springer Science & Business Media, 2012.
[22] Lehrenfeld C. Hybrid discontinuous Galerkin methods for solving incompressible flow problems. Rheinisch-Westfalischen
Technischen Hochschule Aachen 2010; .
[23] Egger H, Schöberl J. A hybrid mixed discontinuous Galerkin method for convection–diffusion problems. J. Numer. Anal
2008; .
[24] Nguyen NC, Peraire J, Cockburn B. An implicit high-order hybridizable discontinuous Galerkin method for linear
convection–diffusion equations. Journal of Computational Physics 2009; 228(9):3232–3254.
[25] Cockburn B, Gopalakrishnan J, Lazarov R. Unified hybridization of discontinuous Galerkin, mixed, and continuous
Galerkin methods for second order elliptic problems. SIAM Journal on Numerical Analysis 2009; 47(2):1319–1365.
[26] Cockburn B, Gopalakrishnan J, Nguyen N, Peraire J, Sayas FJ. Analysis of HDG methods for Stokes flow. Mathematics
of Computation 2011; 80(274):723–760.
[27] Cesmelioglu A, Cockburn B, Nguyen NC, Peraire J. Analysis of HDG methods for Oseen equations. Journal of Scientific
Computing 2013; 55(2):392–431.
[28] Zhai Q, Zhang R, Wang X. A hybridized weak Galerkin finite element scheme for the Stokes equations. Science China
Mathematics 2015; :1–18.
[29] Pietro DAD, Ern A, Lemaire S. An arbitrary-order and compact-stencil discretization of diffusion on general meshes based
on local reconstruction operators. Comp. Meth. Appl. Math. 2014; 14(4):461–472.
[30] Di Pietro D, Ern A, Lemaire S. A review of hybrid high-order methods: formulations, computational aspects, comparison
with other methods. Technical Report HAL-01163569 September 2015.
[31] Cockburn B, Pietro DD, Ern A. Bridging the hybrid high-order and hybridizable discontinuous galerkin methods. Technical
Report HAL-01115318 July 2015.
[32] Egger H, Waluga C. hp analysis of a hybrid DG method for Stokes flow. IMA Journal of Numerical Analysis 2012; :drs018.
[33] Nguyen N, Peraire J, Cockburn B. A hybridizable discontinuous Galerkin method for Stokes flow. Computer Methods in
Applied Mechanics and Engineering 2010; 199(9):582–597.
[34] Carrero J, Cockburn B, Schötzau D. Hybridized globally divergence-free LDG methods. part i: The Stokes problem.
Mathematics of computation 2006; 75(254):533–563.
[35] Cockburn B, Gopalakrishnan J. Incompressible finite elements via hybridization. part i: The Stokes system in two space
dimensions. SIAM Journal on Numerical Analysis 2005; 43(4):1627–1650.
[36] Cockburn B, Gopalakrishnan J. Incompressible finite elements via hybridization. part ii: The Stokes system in three space
dimensions. SIAM Journal on Numerical Analysis 2005; 43(4):1651–1672.
[37] Könnö J, Stenberg R. Numerical computations with H(div)-finite elements for the Brinkman problem. Computational
Geosciences 2012; 16(1):139–158.
[38] Qiu W, Shi K. An HDG method for linear elasticity with strong symmetric stresses. arXiv preprint arXiv:1312.1407 2013;
.
[39] Oikawa I. A hybridized discontinuous Galerkin method with reduced stabilization. Journal of Scientific Computing 2014;
:1–14.
[40] Oikawa I. A reduced HDG method for the Stokes equations. arXiv preprint arXiv:1502.01833 2015; .
[41] Qiu W, Shi K. A superconvergent HDG method for the incompressible Navier-Stokes equations on general polyhedral
meshes. CoRR 2015; abs/1506.07543.
[42] Kirby RM, Sherwin SJ, Cockburn B. To CG or to HDG: a comparative study. Journal of Scientific Computing 2012;
51(1):183–212.
[43] Huerta A, Angeloski A, Roca X, Peraire J. Efficiency of high-order elements for continuous and discontinuous Galerkin
methods. International Journal for Numerical Methods in Engineering 2013; 96(9):529–560.
[44] Ascher UM, Ruuth SJ, Wetton BT. Implicit-explicit methods for time-dependent partial differential equations. SIAM
Journal on Numerical Analysis 1995; 32(3):797–823.
[45] Ascher UM, Ruuth SJ, Spiteri RJ. Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations.
Applied Numerical Mathematics 1997; 25(2):151–167.
[46] Kanevsky A, Carpenter MH, Gottlieb D, Hesthaven JS. Application of implicit–explicit high order Runge–Kutta methods
to discontinuous-Galerkin schemes. Journal of Computational Physics 2007; 225(2):1753–1781.
[47] Strang G. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis 1968; 5(3):506–
517.
[48] Maday Y, Patera AT, Rønquist EM. An operator-integration-factor splitting method for time-dependent problems: application to incompressible fluid flow. Journal of Scientific Computing 1990; 5(4):263–292.
[49] Chrispell J, Ervin V, Jenkins E. A fractional step θ-method for convection–diffusion problems. Journal of Mathematical
Analysis and Applications 2007; 333(1):204–218.
[50] Zaglmayr S. High order finite element methods for electromagnetic field computation. PhD Thesis, JKU Linz 2006.
[51] Arnold DN. An interior penalty finite element method with discontinuous elements. SIAM Journal on Numerical Analysis
1982; 19(4):742–760.
[52] Schöberl J, Lehrenfeld C. Domain decomposition preconditioning for high order hybrid discontinuous Galerkin methods
on tetrahedral meshes. Advanced finite element methods and applications. Springer, 2013; 27–56.
[53] Brezzi F, Manzini G, Marini D, Pietra P, Russo A. Discontinuous finite elements for diffusion problems. Atti Convegno in

23

onore di F. Brioschi (Milano 1997), Istituto Lombardo, Accademia di Scienze e Lettere 1999; :197–217.
[54] Linke A. On the role of the Helmholtz decomposition in mixed methods for incompressible flows and a new variational
crime. Computer Methods in Applied Mechanics and Engineering 2014; 268:782–800.
[55] Dubiner M. Spectral methods on triangles and other domains. Journal of Scientific Computing 1991; 6(4):345–390.
[56] Glowinski R. Finite element methods for incompressible viscous flow. Handbook of numerical analysis 2003; 9:3–1176.
[57] Chorin AJ. The numerical solution of the Navier-Stokes equations for an incompressible fluid. Bulletin of the American
Mathematical Society 1967; 73(6):928–931.
[58] Schäfer M, Turek S, Durst F, Krause E, Rannacher R. Benchmark computations of laminar flow around a cylinder. Flow
simulation with high-performance computers II 1996; :547–566.
[59] Bayraktar E, Mierka O, Turek S. Benchmark computations of 3D laminar flow around a cylinder with CFX, OpenFOAM
and FeatFlow. International Journal of Computational Science and Engineering 2012; 7(3):253–266.
[60] John V. On the efficiency of linearization schemes and coupled multigrid methods in the simulation of a 3d flow around a
cylinder. International Journal for Numerical Methods in Fluids 2006; 50(7):845–862.
[61] ngsflow: Flow solver package for NGSolve. www.sourceforge.net/projects/ngsflow/.
[62] Schöberl J. C++11 implementation of finite elements in NGSolve. Technical Report ASC-2014-30, Institute for Analysis and Scientific Computing September 2014. URL http://www.asc.tuwien.ac.at/~schoeberl/wiki/publications/
ngs-cpp11.pdf.
[63] Kovasznay L. Laminar flow behind a two-dimensional grid. Mathematical Proceedings of the Cambridge Philosophical
Society, vol. 44, Cambridge Univ Press, 1948; 58–62.
[64] PARDISO Solver project. www.pardiso-project.org.
[65] Schenk O, Gärtner K. Solving unsymmetric sparse systems of linear equations with PARDISO. J. of Future Generation
Computer Systems 2004; 20:475–487.
[66] FEATFLOW Finite element software for the incompressible Navier-Stokes equations. www.featflow.de.
[67] Turek S, Becker C. FEATFLOW: Finite element software for the incompressible Navier-Stokes equations, User Manual
Release 1.1. Preprint, Heidelberg 1998; 4.

24

