HDG for incompressible flow - RWTH Aachen …...Chapter 2translates the discretization discussed in...

Hybrid Discontinuous Galerkin methodsfor solving incompressible flow problems

Diplomarbeitzur Erlangung des akademischen Grades

Diplomingenieurin der Studienrichtung

Computational Engineering Sciencean der

Rheinisch-Westfälischen Technischen Hochschule Aachen

Autor:Christoph Lehrenfeld

Betreuer:Prof. Dr. Joachim Schöberl

Prof. Dr. Arnold Reusken

June 6, 2010

1

Abstract

This thesis deals with a higher order discretization of incompressible flow problems using the HybridDiscontinuous Galerkin Method. It aims to introduce an appropriate and computationally efficientHybrid Discontinous Galerkin formulation for the most important model problems of incompressiblefluid flow, namely the convection diffusion equation and the incompressible Navier-Stokes equations.

The main contribution is the derivation, discussion and analysis of the Hybrid (exactly divergence-free)Discontinuous Galerkin method which results from the combination of the following three modern andpromissing concepts in the field of (Discontinuous) Finite Element Methods:

The first of those concepts goes back to [CGL08] where the hybridization of Discontinuous GalerkinMethods is proposed resulting in the Hybrid Discontinuous Galerkin Method which by constructioncompensates for the computational disadvantages of Discontinuous Galerkin Method and thereby al-lows for a more efficient use its advantages. Linear systems of equations stemming from a HybridDiscontinuous Galerkin formulation can be reduced to interface unknowns similar to the static con-densation approach. This reduces the computational effort of Discontinuous Galerkin Method.

The second concept is motivated by [CKS05] where it was shown that Discontinuous Galerkin Methodsfor the Navier-Stokes equations should provide exactly divergence-free solutions in order to be stableand locally conservative. This leads to the idea of using divergence-conforming Finite Elements, whichare normal-continuous across element interfaces, instead of using completely discontinous ones.

In [Zagl06] and [SZ05] a strategy for the construction of conforming higher order Finite Elements forH1, H(curl), H(div) and L2 based on the exact sequence property of those spaces is proposed anddiscussed. This approach leads - beside other advantages - to a natural separation of the H(div)-conforming Finite Elements into basis functions which are relevant for the discrete solution of theincompressible Navier-Stokes equations and those which can be neglected (a priori). This leads to areduction of the number of unknowns.

2

Acknowledgements

At this point I want to, first of all, thank Prof. Dr. Joachim Schöberl for supervising,inspiring and supporting this diploma thesis as well as many other projects during thelast two years. At the same time, I want to thank Prof. Dr. Arnold Reusken fortaking the official supervision of my diploma thesis on him.

I also want to thank Armin Westerkamp for many hours of proof-reading this workwhich led to several enhancements and corrections. Many errors in english grammarand especially punctuation wouldn’t have been found without him.

The groups of the course “Advanced Finite Element Methods” in the summer semester2009, which I was allowed to supervise, helped me to get to know another perspectiveon several aspects of the topics discussed in this diploma thesis. I want to thank allparticipants for a fruitful cooperation.

Furthermore, I want to acknowledge the scientific environment and hostage of theInstitute MathCCES, chaired by Prof. Dr. Joachim Schöberl and Prof. Dr. MartinFrank.

3

Eidesstattliche Erklärung

Ich versichere hiermit, die vorgelegte Arbeit in dem gemeldeten Zeitraum ohne fremdeHilfe verfaßt und mich keiner anderen als der angegebenen Hilfsmittel und Quellenbedient zu haben.Diese Diplomarbeit wurde bisher weder im In- noch im Ausland als Prüfungsarbeitvorgelegt.

Aachen, den June 6, 2010

(Vorname, Nachname)

4 CONTENTS

Contents

Introduction 6

1 Hybrid DG for steady convection diffusion problems 9

1.1 The scalar convection diffusion problem . . . . . . . . . . . . . . . . . . . . . . . . . 9

1.1.1 State of the Art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

1.2 Hybrid DG for second order elliptic problems . . . . . . . . . . . . . . . . . . . . . . 11

1.2.1 Introducing the method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.2.2 Characterizations of the method . . . . . . . . . . . . . . . . . . . . . . . . . 15

1.2.3 A priori error analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

1.2.4 Conservation property of the Hybrid DG method . . . . . . . . . . . . . . . . 32

1.2.5 Numerical Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

1.3 Hybrid DG for (linear) hyperbolic problems . . . . . . . . . . . . . . . . . . . . . . . 36




1.4 Adding diffusion and convection together . . . . . . . . . . . . . . . . . . . . . . . . 48

1.4.1 Joining the limits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48



1.5 Computational aspects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

1.5.1 Elimination of inner degrees of freedom . . . . . . . . . . . . . . . . . . . . . 55

1.5.2 Sparsity Pattern of HDG compared to DG methods . . . . . . . . . . . . . . . 56

2 Divergence-free Hybrid DG for Navier-Stokes equations 58

2.1 A semi-conforming Finite Element space for Navier-Stokes . . . . . . . . . . . . . . . 59

2.2 H(div)-conforming Finite Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

CONTENTS 5

2.2.1 Lowest Order H(div)-conforming Finite Elements . . . . . . . . . . . . . . . . 61

2.2.2 H(div)-conforming Transformation (Piola Transformation) . . . . . . . . . . . 63

2.2.3 The Exact sequence and the DeRham Complex . . . . . . . . . . . . . . . . . 64

2.2.4 Construction of Higher Order H(div)-conforming Finite Elements . . . . . . . 66

2.2.5 H(div)-conforming Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . 70

2.3 Divergence-free Hybrid DG for Stokes Equations . . . . . . . . . . . . . . . . . . . . 71



2.4 Divergence-free Hybrid DG for Oseen Equations . . . . . . . . . . . . . . . . . . . . 83



2.5 Divergence-free Hybrid DG for Navier-Stokes Equations . . . . . . . . . . . . . . . . 86



3 Time Integration for incompressible flow problems 89

3.1 Semi-discrete form of convection diffusion type problems . . . . . . . . . . . . . . . . 89

3.2 Approaches for Time Integration: explicit vs. implicit . . . . . . . . . . . . . . . . . . 90

3.2.1 Explicit methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

3.2.2 Implicit methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

3.2.3 Semi-Implicit methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

3.2.4 Stability vs. Efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

3.3 Higher Order methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

3.4 IMEX Time Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

3.4.1 IMEX Runge-Kutta methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

3.4.2 Multistep IMEX methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

A Notation, elementary ineq. and orthogonal polynomials 98

A.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

A.2 Elementary Inequalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

A.3 Orthogonal polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

B Selected Proofs 104

Bibliography 108

6 CONTENTS

Introduction

Motivation

The numerical solution of partial differential equations (PDE) related to (incompressible) flow problemsis an important issue in many engineering applicitions. Facing the task of solving partial differentialequations numerically, one has to choose between quite many different available methods. In Com-putational Fluid Dynamics the Finite Volume Method as well as the Standard (Continuous Galerkin)Finite Element Methods (CG FEM) are often used as they - in contrast to the Finite Difference Method- also work on unstructured meshes.The concept of Standard (Continuous Galerkin) Finite Element Methods (CG FEM) is based upon thevariational formulation of the PDE. An approximation to the PDE is achieved by replacing the infinitedimensional space in which the variational formulation is posed by a finite dimensional subspace whichnormally uses (element-) piecewise polynomials which are continuous across element interfaces. Thisallows for increasing the order of approximation quite easily also on unstructured meshes as only thepolynomial degree has to be increased (given that the solution is sufficiently smooth). However, theCG FEM suffers from stability issues for convection-dominated flows and the lack of a conservationproperty. Finite Volume Methods (FVM) on the other hand, which are frequently used in engineeringpractice and are build upon the balance of in- and outfluxes on a control volumes, have neither ofthese problems. But here the problem is that providing higher order approximations with Finite VolumeMethods is hardly possible.A modern approach which combines the features and overcomes the draw-backs of both methods is theuse of Discontinuous Galerkin Finite Element Methods (DG FEM). It is also based upon a variationalformulation but this time the finite dimensional space is not taken to be a subspace of the infiniteone. Instead, the finite one uses discontinuous functions, which are not allowed in the continuoussetting. As the DG FEM approximation with piecewise constant functions coincides with the use of aFinite Volume Method the Finite Volume method is, in a sense, a subclass and thereby its conservationproperty can be inherited to the DG FEM. Although all stated problems of FVM and CG FEM can besolved with DG FEM, DG FEM has one (big) disadvantage:It introduces considerably more unknowns as CG FEM or FVM. Especially linear systems arising fromdiscretizations with DG FEM are considerably larger and less sparse as a lot of couplings betweenunknowns of different elements are introduced.One approach to tackle this problem is the concept of hybridization which has already been usedfor discretization such as mixed methods (see [BF91]). The general concept which also allows thehybridization of DG FEM was proposed in [CGL08]. This approach leads to the Hybrid DiscontinuousGalerkin Finite Element Method (HDG FEM) which is the major subject of this diploma thesis. As wewill see, HDG FEM is actually a subset of DG FEM, but allows for a more efficient solution of linearsystems stemming from the discretization.

CONTENTS 7

For the incompressible Navier-Stokes equations it was shown in [CKS05] that DG FEM discretiza-tions have to provide exactly divergence-free solutions to ensure a locally conservative, energy-stableand optimally convergent method. This motivates the use of divergence-conforming DiscontinuousGalerkin Finite Elements and, to combine this with the aspects mentioned before, to transform it toa HDG FEM formulation. Furthermore the construction of higher order divergence-conforming FiniteElements possesses another potential. The construction of those Finite Elements can be based on theexact sequence property of the Sobolev spaces H1, H(curl), H(div) and L2 which allows for furtherimprovement of the HDG FEM formulation for the incompressible Navier-Stokes equations. For thistype of construction, a natural separation of the Finite Element space takes place. One subset of thisFinite Element space can then be removed without consequences on the approximation of (strongly)divergence-free velocity fields. This reduces the total number of unknowns compared to other dis-cretizations.

The above adressed issues are only problems in the spatial discretization and adding time integrationwith the method of lines seems to be clear. Here you normally use either fully implicit or explicitmethods, like Runge-Kutta schemes. But the structure of convection diffusion equations - and theNavier-Stokes equations belong to this class of equations - can be considered for by semi-implicit timeintegration methods. Those combine the features of implicit and explicit methods resulting in a methodwhich has to solve linear systems of equations which are far more easy to handle than those arising infully implicit schemes but whose time step restriction - basically - only depends on the convection termof the PDE which is a mild condition compared to time step restrictions for fully explicit methods.As this approach is not as famous as the fully explicit and fully implicit counterparts, we will at leastpresent higher order time integration schemes of this kind, which have been proposed in [ARW95] and[ARS97], the so called IMEX (implicit-explicit) schemes.

As two model problems on which we’ll apply the above indicated strategies we consider the scalarconvection-diffusion equation and the incompressible Navier-Stokes equations.

Organization of the work

The thesis is organized as follows:

Chapter 1 derives an HDG FEM formulation for the steady scalar convection diffusion equation byconsidering an HDG FEM formulation for each limit case, i.e. the pure convective one (linear transportequation) and the pure diffusive one (Poisson’s equation), first and joining them afterwards. Com-parisons and relations to other (e.g. Discontinuous or Mixed) Finite Element Methods are discussedand a rather complete a priori error analysis is carried out. In addition, numerical examples and thediscussion of properties and computational aspects substantiate the potential of Hybrid DiscontinuousGalerkin Methods.

Chapter 2 translates the discretization discussed in chapter 1 to the steady incompressible Navier-Stokes equations. Before the new divergence-conforming HDG formulation is introduced an extrasection about the constructing of higher order H(div)-conforming Finite Elements is given. This sec-tion includes the discussion of the construction proposed and discussed in [SZ05, Zagl06] which allowsfor the reduction of the number of unknowns. Starting from the Stokes problem the new discretization

8 CONTENTS

of the steady incompressible Navier-Stokes equations is analyzed and discussed.

Chapter 3, which is about time integration (especially of semi-implicit type) for Navier-Stokes andconvection diffusion equations, brings the diploma thesis to the conclusion. Therein the IMEX timeintegration method is derived and two different approaches to achieve higher order time integration ofsemi-implicit type are suggested.

The Appendix includes a summary of the notation used through out the thesis, some elementaryinequality and the definition of some orthogonal polynomials. Furthermore two (more involved) proofswhich were moved from chapter 2 can be found at the end of the appendix.

Requirements on the reader

The reader is presumed to have basic knowledge of the Finite Element Method including basic knowl-edge of functional analytic concepts. No knowledge of nonconforming Finite Elements like DG FEMis required.

Implementations

All implementations used for the numerical examples in chapter 1 were carried out with the opensource software package NETGEN/NGSolve. NETGEN is an automatic mesh generator which is availableat

http://sourceforge.net/projects/netgen-mesher

The Finite Element Library NGSolve which comes as an Add-On to NETGEN is also open source andavailable at

http://sourceforge.net/projects/ngsolve

It already includes the proposed HDG discretization for the scalar problem.Implementations of the Stokes and Navier-Stokes discretization discussed in chapter 2 and also thetime integration schemes of chapter 3 can be found in the NGSolve-Add-On package ngsflow whichis developped by Joachim Schöberl and the author. It is available under

http://sourceforge.net/projects/ngsflow

and is able to solve the steady and unsteady incompressible Navier-Stokes problem on unstructured(hybrid) 2D and 3D meshes with arbitrary polynomial degree.

9

Chapter 1

The Hybrid Discontinuous Galerkin methodfor steady scalar convection diffusionproblems

1.1 The scalar convection diffusion problem

Convection diffusion equations model the transport of a scalar quantity - denoted by u - which maybe a concentration of a species or the temperature. It is a conservation law which accounts for twosuperimposed mechanisms: The first one, the convection, is the transport of the scalar quantity alonga given stream b and is a directed phenomena (information moves downwind, information speed islimited). The corresponding convective flux is simply bu. Convection is superimposed by diffusionwhich describes the exchange of momentum, concentration or temperature of neighbouring particleson a molecular level and so is not directed (information goes in all directions and propagates infinitelyfast). The corresponding diffusive flux is −ε∇u, where ε is a scalar for isotropic diffusion, which is thestandard case. The right hand side of the convection diffusion equation describes sources and sinks forthe scalar quantity. So if you balance the change of the scalar in time and the total in- and outgoingfluxes with sources and sinks on an arbitrary control volume you end up with the convection diffusionequation:

∂u

∂t+ div(−ε∇u+ bu) = f

This type of equation arises in many flow problems and other physical applications like semiconductors.The steady Navier-Stokes equations which will be discussed later in chapter 2 are a vector-valued versionof the convection diffusion problem, where the momentum ρu is the quantity which is transported withthe velocity u and superimposed by viscous effects which are of diffusive character.Our objective in this part is to solve the following scalar convection diffusion equation (with appropriateboundary and initial conditions) numerically. The problem reads

∂u∂t

+ div(−ε∇u+ bu) = f in Ω× [t0, T ], div(b) = 0u = uD on ΓD × [t0, T ]∂u∂n

= g on ΓN × [t0, T ]∂u∂n

+ βu = h on ΓR × [t0, T ]u = u0 on Ω× t0

(1.1.1)

10 CHAPTER 1. HYBRID DG FOR STEADY CONVECTION DIFFUSION PROBLEMS

where Ω is a bounded open domain in Rd, d = 1, 2, 3 with boundary ∂Ω = ΓD ∪ ΓN ∪ ΓR. Theboundary parts ΓD, ΓN and ΓR denote the part of the boundary where Dirichlet, Neumann and Robinboundary conditions are prescribed respectively. [t0, T ] ⊂ R describes the time interval, where t0 isthe start time where the initial condition u0 is prescribed and T is some later time.The spatial domain is decomposed into a shape regular partition Th of Ω consisting of simplices T .The element interfaces and element boundaries coinciding with the domain boundary are called facets.Furthermore the set of those facets E is denoted by Fh and there holds ⋃T∈Th ∂T = ⋃

E∈Fh E.Note that we always assume that the convective velocity b is divergence-free which is reasonable if weassume that it comes from an incompressible flow field.In this chapter we will only consider steady problems, i.e. ∂u

∂t= 0 for now. The structure of this

chapter is as follows:First we will look at one limit of the convection diffusion equation, the poisson’s equation, where theconvective velocity is set to zero (‖b‖ = 0) and only the purely diffusive phenomenon is described.The Hybrid DG formulation is introduced, discussed and a complete a priori error analysis is carriedout. The other limit of the convection diffusion equation (ε = 0), which is the pure convective (linearand hyperbolic) one, will be the topic of the next section where the discretization of the convectionwith hybrid DG is presented and analyzed. Both cases will be joined afterwards to handle the fullconvection diffusion equation and the a priori error analysis for this case will be reconstructed fromjoining the analysis of both limits. As the major advantages and key properties of the Hybrid DGmethod are the same in all three parts, at least from a computational point of view, computationalaspects will be discussed in more detail after the a priori error analysis. We conclude each section withnumerical examples which shall highlight the features of the proposed Hybrid DG method.

1.1.1 State of the Art

For most elliptic problem the standard Finite Element method which we will call CG for “continuousGalerkin” is a convenient choice. Nevertheless there are situations where discontinuous formulationshave their advantages. The CG method lacks desireable properties like conservation and the abilityto represent “nearly discontinuities” on coarse meshes which might come up at material interfacesfor example. An additional problem appears as soon as you leave the pure elliptic problems to con-vection diffusion problems. In the pure convective limit Discontinuous Galerkin methods which wecall DG methods seem to be the more convenient choice, then. So you probably want to use amethods which can naturally be extended by DG formulations. To circumvent problems of the stan-dard Finite Element method with convection diffusion equations there are conforming approaches likestreamline upwind Petrov Galerkin, Galerkin least squares and adjoint Galerkin least squares methods(see [ESW05],[DH03]) which base upon the idea of adding diffusion numerically only in streamlinedirection and only in a consistent manner by formally exchanging the test functions of the variationalformulation. Those methods certainly solve the stability problems of CG, but are often criticized to beoverdiffusive.A more general approach is the use of Standard DG formulations for elliptic operators like they arediscussed in [ABCM02]. As hyperbolic equations can be included more naturally to DG methods thanto conforming methods, the way to convection diffusion equations for DG methods is easy. Neverthe-less this convenience comes with the prize of more unknowns and larger stencils.Another approach was proposed by Egger and Schöberl in [ES09] where a mixed method for the dif-fusive part in hybridized form is coupled with a hybridized DG upwind method for the convective part.Here, even more unknowns have to be considered as in the DG case, but the stencil and so the nonzeroentries in the matrices are dramatically reduced.

1.2. HYBRID DG FOR SECOND ORDER ELLIPTIC PROBLEMS 11

Our approach is strongly connected to the last two approaches. We use a Hybrid DiscontinuousGalerkin method for the elliptic problem like they are presented and discussed in [CGL08] and thesame hybridized DG upwind method as in [ES09] and couple both.Another also promissing idea is presented in [DG10]. Therein a Discontinuous Galerkin method withtest functions differing from the trial functions is used, where the test functions are specifically tailoredfor stability resulting in a Discontinuous Petrov-Galerkin formulation which achieves L2 best approxi-mation for the linear transport equation. Thus the method also convergences optimally in the L2 normwhich does not hold for the above mentioned methods (including ours).

1.2 The Hybrid Discontinuous Galerkin method for second or-der elliptic problems

1.2.1 Introducing the method

In this section we consider poisson’s equation:−div(ε∇u) = f in Ω

u = uD on ΓD∂u∂n

= g on ΓN∂u∂n

+ βu = h on ΓR

(1.2.1)

As we want to use a discontinuous Finite Element approximation of the weak solution u ∈ H1(Ω), welook at the space of (element-)piecewise H1(T ) functions which form the broken Sobolev space

H1(Th) := u ∈ L2(Ω), u ∈ H1(T ) ∀T ∈ Th

On each element we approximate the H1(T ) functions u with functions uh which are polynomials upto degree k ∈ N, but discontinuous. Thus there only holds

uh|T ∈ Pk(T ) ∀ T ∈ Th ⇒ uh ∈ H1(Th) ⊂ L2(Ω) (1.2.2)

When solving poisson’s equation with CG you multiply the equation by a test function and integrateby parts on the whole domain to obtain a weak formulation. Due to the fact that the functions are nolonger assumed to be H1(Ω) regular but only H1(Th) regular, partial integration on the whole domaincan no longer be applied. Thus, we integrate by parts on each element T ∈ Th, instead.Starting with the problem

Find u ∈ H1(Ω) s.t. −∫

Ωdiv(ε∇u) v dx =

∫Ωf v dx ∀v ∈ H1

0 (Ω)

with H10 (Ω) := u ∈ H1(Ω), u = 0 on ΓD

and formally doing integration by parts on each element T ∈ Th we get a new variational problem foru ∈ H1(Ω) and v ∈ H1(Th)

−∑T∈Th

∫Tdiv(ε∇u) v dx =

∑T∈Th

∫Tε ∇u ∇v dx−

∫∂Tε∂u

∂nv ds

=

∑T∈Th

∫Tf v dx (1.2.3)


s

s

JJJJJJ

s

s

s @@

@@

s(a) linear 2D CG DOF

s

s

JJJJJJ

s

s

s @@@@

s(b) linear 2D DG DOF

s

s

JJJJJJ

s

c

c

s

s @@

@@

s(c) linear 2D HDG DOF

Figure 1.2.1: Degrees of freedom for different schemes

Note that we don’t move the integral terms of the domain boundary on the right hand side. Thus wedon’t get Neumann boundary conditions as natural boundary conditions yet.At this point we introduce additional unknown functions uF ∈ L2(Fh). Those functions only exist onthe facets and for the true solution they represent the trace of u. The computational advantages ofthose additional facet functions, which are the essential ingredients of Hybrid DG, will be discussedlater. In order to illustrate the situation, the degrees of freedom(DOF) of this method compared tothose of CG and DG (for linear approximations in 2D) are schematically drawn in figure 1.2.1.The above formulation obviously requires more regularity than u ∈ H1(Th), namely u ∈ Hs(Th) withs ≥ 3

2 , in order to make the normal derivative well defined. But, for ease of presentation, we willassume u ∈ H2(Th).The assumption u ∈ H1(Ω) is now extended by an additional H2(Th) regularity assumption. Thespaces we will use on the continuous and the discrete level are V and Vh and are defined in thefollowing.

V :=

(u, uF ) : u ∈ H2(Th) ∩H1(Ω), uF ∈ L2(Fh)

Vh :=

(u, uF ) : u ∈ Pk(T ) ∀ T ∈ Th, uF ∈ Pk(E) ∀ E ∈ Fh (1.2.4)

Here, u ∈ H2(Th) implies that u is in H2(T ) for every element T of Th, but no continuity constraintsare imposed on element interfaces. Dirichlet boundary conditions are posed on the facet functionsonly. Then, the spaces with (inhomogeneous and homogeneous) Dirichlet boundary conditions for thefacet unknowns in it1 can be defined analogously:

VD := (u, uF ) ∈ V, uF = uD on ΓD V0 := (u, uF ) ∈ V, uF = 0 on ΓDVh,D := (u, uF ) ∈ Vh, uF = uD on ΓD Vh,0 = (u, uF ) ∈ Vh, uF = 0 on ΓD

(1.2.5)

With these choices for the discrete spaces another problem with the current form of the new formulationbecomes evident: functions of neighbouring elements are independent from each other, such that all(element-)piecewise constant functions liein the kernel of the left hand side bilinearform. So, we haveto put some additional work into the formulation to overcome this problem. In essence, we achievethat by adding some consistent terms to the bilinearform. They either make use of the fact, that theexact solution is continuous across element interfaces or that the flux −ε∇u is normal-continuous.

1imposing boundary conditions strong for the facet functions is like imposing them weakly in a Standard DG manner(see section 1.2.2.2)


The latter justifies the following relation:

∑T∈Th

∫∂Tε∂u

∂nvF ds =

∑E∈Fint

h

∫Eε∂u

∂n+ vF ds+∫Eε∂u

∂n−vF ds

︸︷︷︸

=0 for true solution

+∑

E∈Fexth

∫Eε∂u

∂nvF ds =

∫∂Ωε∂u

∂nvF ds

(1.2.6)

where we used F ext := E ∈ Fh; E ∩ ∂Ω 6= ∅ and F int := E ∈ Fh; E ∩ ∂Ω = ∅ for the set ofall interior and exterior facets.For inner facets we get element boundary integrals from both sides which cancel out for the truesolution. Thus if we add (1.2.6) to (1.2.3) and plug in the boundary conditions for the new integralterms and introduce [[v]] := v−vF where the abbreviation v in contrast to the expression (v, vF ) iswritten bold:

Find (u, uF ) ∈ VD, such that for all (v, vF ) ∈ V0 there holds:∑T∈Th


∫∂Tε∂u

∂n[[v]] ds

+ β

∫ΓRε u vF ds

=∑T∈Th

∫Tf v dx+

∫ΓNε g vF ds+

∫ΓRε h vF ds

(1.2.7)

Due to (1.2.6), Neumann boundaries appear as natural boundary conditions in this variational for-mulation again. We add additional terms, which now make use of the H1(Ω)-regularity of the exactsolution, which ensures that for (u, uF ) ∈ V there holds [[u]] = u − uF = u − tr∂T (u) = 0 on eachfacet. Thus by adding integrals on each facet, which contain this term, we add only consistent terms,and then element neighbours couple, even though only indireclty, with each other. We choose to add

−∫∂Tε∂v

∂n[[u]] ds for symmetry

and ∫∂Tτhε [[u]] [[v]] ds for stability

with some stabilization parameter τh, which has to be chosen appropriately.We can also make use of the consistency relation to exchange u with uF at the Robin boundaryintegral. After those manipulations we get the bilinearform

Bh((u, uF ), (v, vF )) :=∑T∈Th


∫∂Tε∂u

∂n[[v]] ds

−∫∂Tε∂v

∂n[[u]] ds +

∫∂Tε τh [[u]] [[v]] ds

+ β

∫ΓRε uF vF ds

(1.2.8)

The insertion of the last two consistent integrals are degrees of freedom to tune and control the methodand so the proposed bilinearform is certainly not the only reasonable choice you could make. In thesection 1.2.2 we will try to understand the properties and the relations of this formulation to othermethods. As this method has some obvious similarities to the Standard DG symmetric interior penaltymethod, we also want to mention that some other DG methods like the NIPG method or the DGmethod by Baumann and Oden (see also [ABCM02]) can be easily translated to a similar hybridized


version.The most important properties of the proposed formulation are symmetry, coercivity for sufficientlylarge τh (we will discuss this later) and consistency. Symmetry is obvious, consistency and coercivitywill be further discussed in section 1.2.3.1 and 1.2.3.2.For ease of notation, we will often denote the pair (u, uF ) as u and also shorten the notation for thebilinearform by Bh(u,v). Thus our discrete problem reads

Find u ∈ Vh,D, such that for all v ∈ Vh,0 :

Bh(u,v) = 〈f , v〉 =∑T∈Th

∫Tf v dx+

∫ΓNε g vF ds+

∫ΓRε h vF ds

(1.2.9)

Several advantages come up when using this kind of Hybrid DG methods:

• Conservation of a discrete flux (see 1.2.4) like Standard DG or mixed methods

• Standard FEM like elementwise matrix assembling as neighbouring degrees of freedom don’tcouple directly (see 1.5.2). Inner degrees of freedom can also be eliminated by static condensa-tion.

• Extensions to DG schemes for convection diffusion problems can be considered for quite easily

Remark 1.2.1 (Stabilization paramater):For a sufficiently large stabilization parameter τh > τ ∗h the method is stable. This parameter dependson the size and shape of the neighbouring elements of a facet. From scaling arguments it is also clearthat τh should scale with 1

hwhere h is a characteristic length of the neighbouring elements. When we

show coercivity in 1.2.3.2 we’ll see that the τ ∗h can be explicitly derived from inverse inequalities.

Remark 1.2.2 (Non-Uniqueness of the Abbreviation “HDG method”):In this work the proposed symmetric interior penalty HDG method will be adressed as the “HDGmethod” even though the correct designation would be the somewhat more lenghty expression “hy-bridized symmetric interior penalty Discontinuous Galerkin method”, which distinguishes it from otherHDG methods that could also be used (see [CGL08]) . Most of the important results presented in thiswork that hold for our particular choice, the hybridized symmetric interior penalty DG method, can betranslated to other HDG methods as well.

Remark 1.2.3 (Dirichlet boundary conditions):If Dirichlet boundary conditions are used corresponding degrees of freedom can be removed from theformulation. This can be done by shifting domain boundary integrals of the bilinearform Bh on ther.h.s.. The necessary adjustments on the linearform 〈f , v〉 then lead to a modified variational problem,which incorporates Dirichlet boundary conditions in a “Nitsche” way. Another possibility is to plug inthe boundary conditions on the matrix-vector level after discretizing the variational formulation . Ofcourse, both ways end up with the same solution. However, our analysis will be carried out includingDirichlet boundary facet functions, thus the bilinearform does not have to distinguish between exteriorand interior facets, which simplifies the notation.


1.2.2 Characterizations of the method

The HDG method is deeply connected to other Finite Element methods and can be derived in differentways. We want to briefly highlight the relation to (hybridized) mixed methods, Standard DG methodsand an alternative HDG formulation. As we will see, the proposed method is actually equal to onespecific hybridized mixed method, one specific Standard DG method and also to the alternative HDGmethod. Nevertheless the methods are derived from different point of views. At the end of the section,after the several relations to other Finite Element methods are discussed, an attemp to visualize theinterconnection between the different formulations, and thereby summarize this section, is made.

For ease of presentation we don’t draw special attention on the treatment of boundary conditions,meaning we assume homogeneous Dirichlet boundary conditions on the whole boundary. However, forall three approaches different boundary conditions can be easily considered for.

1.2.2.1 HDG as a modified hybridized mixed method

Mixed methodsAs a starting point we reformulate the second order problem −div(ε∇u) = f , equation (1.2.1) as afirst order system by substituting the flux σ = −ε∇u

1εσ + ∇u = 0

div(σ) = f+ b.c. (1.2.10)

The appropriate spaces for a weak solution are H(div,Ω) for σ and L2(Ω) for u. After test with τ andv) and applying partial integration on the gradient part, we end up with the standard mixed methodfor the poisson equation2,3

(1εσ, τ)Ω − (u, div(τ))Ω = (−uD, τn)∂Ω ∀ τ ∈ H(div,Ω)

−(div(σ), v)Ω = (−f, v)Ω ∀ v ∈ L2(Ω) (1.2.11)

The appropriate discrete spaces Σh and Qh are H(div,Ω)-conforming and discontinous Finite Elementspaces, respectively.

Σh ⊂ H(div,Ω) , Qh ⊂ L2(Ω)

If we choose Σh to be an H(div,Ω)-conforming finite element space of order4 k+ 1, which we’ll labelwith Σk+1

h and Qh to be the piecewise polynomial space up to order k, which we’ll label with Qkh,

you can show stability and optimal convergence rates for the error (‖u− uh‖2L2 + ‖σ − σh‖2

L2)12 (see

standard text books, e.g. [Brae97]). The construction of H(div,Ω)-conforming Finite Element spaceswill be discussed in more detail in section 2.2.One advantage of this method is the conservation of the flux σ. But the disadvantage compared toCG is that you have solve a saddle point problem with considerably more degrees of freedom. Thisdisadvantage can be circumvented with hybridized mixed methods, which we’ll discuss next.

2for standard mixed methods you often choose σ = ε∇u and so you end up with different signs3when looking at mixed systems we will use the inner product notation instead of integrals4The order of an H(div,Ω)-conforming space describes the maximum polynomial order of the polynomials itself not

of the divergence. RT0 and BDM1 are, by this definition, examples for Σ1h.


Hybridized mixed methodsFor mixed methods the primal unknown is chosen to be discontinuous and the flux is introduced as anew unknown with normal continuity incorporated in the space Σh, so we directly get conservation ofthe mixed method. When using hybridization techniques, you break up the normal-continuity constraintfor the Finite Element space and use lagrange multipliers to reimpose the according constraint. Weachieve that by multiplying (1.2.10) with testfunctions and integrating elementwise over the wholedomain. When doing integration by parts on every element, boundary terms for all element interfacesappear. We then define the trace of u on the boundary of each element as new unknowns uF and usethem as lagrange multipliers for the conservation-constraint [[σ]]n := σ+n+ + σ−n− = 0 (in a weaksense):Find (σ, u, uF ) ∈ Σk+1

h,disc ×Qkh × F k

h , such that∑T∈Th (1

εσ, τ)T −(u, div(τ))T +(uF , τn)∂T = 0 ∀τ ∈ Σk+1

h,disc∑T∈Th −(div(σ), v)T = (−f, v)Ω ∀v ∈ Qk

h∑T∈Th (σ · n, vF )∂T = 0 ∀vF ∈ F k

h

(1.2.12)

where F kh := v ∈ L2(Fh); v|∂T ∈ Pk is the space of (facet-) piecewise polynomials. All element

unknowns u and σ only couple with other unknowns of the same element or the facet unknowns uF .The (element-)local problems are furthermore uniquely solvable (see [BF91]) and so the linear systemcan be reduced to the facet unknowns by static condensation. A remarkable property of the final linearsystem for the facet unknowns is the symmetric positive definiteness.Until this point, you might see the hybridization as an implemenatation trick. But after solving thelinear system for the facet unknowns, you can enhance your discrete solution even more: The solutionuF , which describes Dirichlet data to local problems, can be used for (element-)local postprocessingslike they are discussed in [BF91] to gain an order k + 1 reconstruction of the primal variable u.

Hybrid DG methodsThe third equation of (1.2.12) ensures continuity of σ and so conservation of the hybrid mixed method.Another possibility of ensuring conservation without enforcing continuity of σ is to pose a conditionon a numerical flux σ instead. To get there, we integrate both div-terms by parts and replace σ · nby σ · n. Now as we did integration by parts we have to adapt our Finite Element spaces. Now wechoose the primal variable u of one order higher than σ, s.t. ∇u and σ lie in the same space. Thenew problem reads:Find (σ, u, uF ) ∈ [Qk−1

h ]d ×Qkh × F k

h , such that∑T∈Th (1

εσ, τ)T + (∇u, τ)T + (uF − u, τn)∂T = 0 ∀τ ∈ [Qk−1

h ]d∑T∈Th (σ,∇v)T − (σ · n, v)∂T = (−f, v)Ω ∀v ∈ Qk

h∑T∈Th (σ · n, vF )∂T = 0 ∀vF ∈ F k

h

(1.2.13)

Here you can use several appropriate fluxes. Choosing σ = σ would lead to a formulation which isequivalent to standard mixed methods again. If we test our primal HDG-IP formulation (1.2.9) and(1.2.8) with (0, vF ), we see that we already (implicitly) used the third equation of (1.2.13), s.t. ournumerical flux is

σ · n = −ε∂u∂n

+ ετh(u− uF ) · n (1.2.14)

With the linear and element-local lifting operator L : L2(∂T )→ H(div, T ) which fulfills∫TL(w) τ dx =

∫∂Tw τ · n ds ∀τ ∈ H(div, T ) (1.2.15)


we can formally rewrite the first equation of (1.2.13) as∫T

(1εσ +∇u

)τ dx+

∫∂T

(uF − u) τ · n ds =∫T

(1εσ +∇u+ L(uF − u)

)τ dx = 0 ∀ τ ∈ [Qk−1

h ]d

(1.2.16)Thus the strong form for σ reads:

σ = −ε∇u+ εL(u− uF ) (1.2.17)

Plugging (1.2.14) and (1.2.17) into the second equation of (1.2.13) we get ∀v ∈ Qkh:∑

T∈Th

∫Tε∇u ∇v dx−

∫TεL(u− uF ) ∇v dx−

∫∂Tε∂u

∂nv ds+

∫∂Tτh(u− uF ) v ds

=∫

Ωf v dx

(1.2.18)If we plug the numerical flux of (1.2.14) also into the third equation of (1.2.13) and multiply by −1we get: ∑

T∈Th−∫Tσ · n vF dx =

∑T∈Th

∫∂Tε∂u

∂nvF ds−

∫∂Tτh(u− uF ) vF ds

= 0 (1.2.19)

Adding equation (1.2.18) and (1.2.19) together, we end up with the primal formulation of Hybrid DGagain (compare (1.2.8))

Bh(u,v) :=∑T∈Th


∫∂Tε∂u

∂n[[v]] ds −

∫∂Tε∂v

∂n[[u]] ds +

∫∂Tετh [[u]] [[v]] ds

(1.2.20)

where u and v are taken out of Qkh × F k

h = V kh

Remark 1.2.4 (Reducing the number of facet unknowns):If we compare our approximation abilities with respect to the facet unknowns’ polynomial degree,we notice that for the hybridized mixed method we can achieve approximation of order k + 1 for u(after postprocessing) when using polynomials of order k on the facets. In (1.2.13) we just achieveapproximation of order k when using polynomials of order k on the facets, so we want to briefly discusshow to overcome this drawback.The reason we take uF and vF of polynomial degree k is that σ · n is of degree k(see third equationof (1.2.13)). If we modify σ · n s.t. it is only of polynomial degree k − 1, the test function vF andso the lagrange multiplier uF can be taken from F k−1

h , what is our goal here. So the global degreesof freedom (meaning those which can not be condensated out on the element level) could be reducedfurther. To achieve this, we could introduce an L2-orthogonal projector Pk−1, which projects into thespace of (facet-) piecewise polynomials of order k − 1 and adjust our numerical flux:

σ · n = −ε∂u∂n

+ ετhPk−1(uh − uF ) (1.2.21)

With this flux and the L2-orthogonality of the projector we get a new bilinearform for which we nowtake u and v in Qk

h × F k−1h :

Bh(u,v) :=∑T∈Th


∫∂Tε∂u

∂n[[v]] ds −

∫∂Tε∂v

∂n[[u]] ds +

∫∂Tτhε Pk−1[[u]] Pk−1[[v]] ds

This “optimization” of the bilinearform only makes sense as long as we have pure elliptic problems. Assoon as we add convection we need uF and vF to be of order k again. That’s why the a priori erroranalysis does not include this “optimization”.


1.2.2.2 HDG as a Standard DG method

For standard Discontinuous Galerkin methods, there is a lot of work going on for many years and alot of theory has been derived for them. In this small part we want to derive a representation of theproposed HDG method which fits into the framework of Standard DG methods. As we presented ourHDG method in a primal formulation, the reformulation as a Standard DG formulation in primal formis the most natural approach. From the characterization of the HDG method as a hybridized mixedmethod, we know that the equations corresponding to the facet degrees of freedom ensure conserva-tion of a numerical flux. But you can also derive an expression for the lagrange multipliers uF fromthose equations. By doing so we eliminate uF , return to Standard DG unknowns and get a primal DGformulation. Afterwards, we will conclude this excursus with the numerical trace u and flux σ of theHDG method.

Eliminating the facet unknowns uFWe eliminate the facet unknowns uF of a discrete solution of the variational problem. We can do thatfor a facet E ∈ Fh by testing equation (1.2.9), (1.2.8) with v = (0, vF ) where vF ∈ L2(E). Addingup the terms of two neighbouring elements we get (with · the arithmetic average of both elementvalues) ∫

E[[ε∇u]] vF ds− 2

∫E

(τhεu − τhεuF ) vF ds = 0 (1.2.22)

In strong form 5 this meansuF = τhεu

τhε− 1

2[[ε∇u]]τhε

(1.2.23)

DG represenation of the HDG-jump operator [[u]]Now we want to make use of this representation and plug it into the bilinearform so that we end upwith a formulation which is directly comparable to other DG methods. For ease of presentation, wewill assume ε ≡ 1 and τh|E ≡ const on each facet E, s.t. uF = u − 1

2 [[∇u]]. Plugging relation(1.2.23) into the bilinearform Bh(·, ·), we get rid of the additional unknowns uF and so the additionaltest functions on the facets are no longer required. All jump terms for the test functions [[v]] = v− vFnow just get v and [[u]] is:

[[u]] = u− uF = 12

([[u]]* + 1

τh[[∇u]]n

)where we used u− u = 1

2 [[u]]* with [[u]]* := u− uext the DG jump operator and uext the neighbourelements’ value6. This gives us the Standard DG Bilinearform

BDGh (u, v) :=∑T∈Th

∫T∇u ∇v dx︸︷︷︸

A

−∫∂T

∂u

∂nv ds︸︷︷︸

B

−∫∂T

∂v

∂n[[u]] ds︸︷︷︸

C

+∫∂Tτh[[u]]v ds︸︷︷︸D

(1.2.24)

Straight forward computation gives us the following relations for the sum of element boundary integralsand according facet integrals. Here the terms with [[u]] are devided into two parts. One with 1

2 [[u]]*

5 To come to the strong form we have to assume that the test space is sufficiently large. In our applications this isalways true.

6 in contrast to [[∇u]]n the sign of [[u]]* depends on the element from where you look at the facet. That’s what weindicate by the star.


which is marked as part 1 and the one with 12τh

[[∇u]]n is marked as part 2follows:

B : −∑T∈Th

∫∂T

∂u

∂nv ds = −

∑E∈Fh

∫E∇u∗[[v]]* + v[[∇u]]n ds

C.1 : −12∑T∈Th

∫∂T

[[u]]* ∂v

∂nds = −

∑E∈Fh

∫E

[[u]]*∇v∗ ds

C.2 : −∑T∈Th

∫∂T

12τh

∂v

∂n[[∇u]]n ds = −

∑E∈Fh

∫E

12τh

[[∇u]]n[[∇v]]n ds

D.1 : 12

∑T∈Th

∫∂Tτh[[u]]*v ds =

∑E∈Fh

∫E

τh2 [[u]]*[[v]]* ds

D.2 : 12

∑T∈Th

∫∂T

[[∇u]]nv ds =∑E∈Fh

∫E

[[∇u]]nv ds

(1.2.25)

Primal DG formulationPlugging all those relations into BDGh (u, v) gives the primal DG form

BDGh (u, v) =∑T∈Th

∫T∇u ∇v dx︸︷︷︸

A

+∑E∈Fh

−∫F∇u∗ [[v]]* ds︸︷︷︸

B+D.2

−∫F∇v∗ [[u]]* ds︸︷︷︸

C.1

+12

∫Fτh[[u]]* [[v]]* ds︸︷︷︸D.1

−12

∫F

1τh

[[∇u]]* [[∇v]]* ds︸︷︷︸C.2

(1.2.26)

If we drop off the last integral and weight the stabilization parameter by a factor of two, we’d re-cover the Standard DG (symmetric) Interior Penalty formulation for the laplace operator. Thus wesee that our proposed HDG (symmetric) Interior Penalty formulation is quiet similar but not identical.

Numerical flux and trace of the HDG methodIn [ABCM02] the Standard DG (symmetric) Interior Penalty as well as other DG methods are comparedby formulating a general Standard DG bilinearform with numerical traces and fluxes u and σ. If youreplace uF in equation (1.2.13) by a numerical trace u you have the general DG formulation in mixedform:

∑T∈Th (1

εσ, τ)T +(∇u, τ)T +(u− u, τn)∂T = 0 ∀τ ∈ [Qk−1

h ]d∑T∈Th (σ,∇v)T −(σ · n, v)∂T = ∑

T∈Th(−f, v)T ∀v ∈ Qkh∑

T∈Th (σ · n, vF )∂T = 0 ∀vF ∈ F kh

(1.2.27)

For the different choices of those fluxes σ and traces u you get different DG formulations. As wasalready pointed out in [CGL08, equation (3.27)] the proposed HDG formulation has the followingnumerical trace and flux:

u = uF = u − 12τh

[[∇u]]n see (1.2.23) with ε ≡ 1 and τh ≡ const on each facet

σ · n = −∇u∗ + τh2 [[u]]* average σ from (1.2.14)

(1.2.28)

Both values are consistent as they coincide with the trace of the solution and the trace of the normalderivative, respectively, for functions u ∈ H2(Ω).


1.2.2.3 Alternative HDG formulation

As we have seen, the proposed representation of our HDG formulation can be regained from (1.2.27)by choosing a specific numerical flux and introducing u = uF as a new unknown. You might also wantto do it the other way around, set σ · n = σn,F and choose the numerical trace u as in (1.2.28). Fornow, we again want to assume ε ≡ 1 and τh|E ≡ const on each facet E. With the new unknowns wealso pose a new condition on the solution, which is 7

∑T∈Th

∫∂Tu τn,F ds−

∫∂T

1τh

(σn,F + ∂u

∂n

)τn,F ds

= 0 (1.2.29)

With the relationsu = u− 1

2[[u]]* (1.2.30a)

12[[∇u]]n = ∂u

∂n− ∇u∗ (1.2.30b)

we can reformulate the numerical trace:

u = u − 12τh

[[∇u]]n(1.2.30)= u− 1

τh

τh2 [[u]]* − 1

τh

(∂u

∂n− ∇u∗

)

= u− 1τh

(∂u

∂n− ∇u∗ + τh

2 [[u]]*

)

= u− 1τh

(∂u

∂n+ σn,F

)(1.2.31)

where we used σn,F = σ · n.Now we can use the representation of the flux (1.2.17) when replacing uF by u. For u− u we get

[[u]] = u− uF = u− u = 1τh

(∂u

∂n+ σn,F

)(1.2.32)

And we can plug it into the representaion of the flux (1.2.17)

σ = −∇u+ 1τhL(∂u∂n

+ σn,F ) (1.2.33)

Further putting this into the second equation of (1.2.27) and adding the new equation (1.2.29), weget for all (v, τn,F ) ∈ Qk

h × F kh :

∑T

∫T∇u ∇v dx+

∫∂Tσn,F v ds−

∫∂Tu τn,F ds+

∫∂T

1τh

(σn,F + ∂u

∂n

)(τn,F + ∂v

∂n

)ds

=∫

Ωfv dx

(1.2.34)Although this formulation is equivalent to our proposed method there are some computational disad-vantages:

1. the problem cannot be reduced to σn,F only, as the local problems are singular.

2. the system is indefinite

This approach was introduced in [EWY02] as a (stabilized) hybridized primal mixed method.7we pose exactly this condition as we aim to get an equivalent formulation to the one in (1.2.8)


1.2.2.4 Summary of the Interconnection between the discussed Finite Element methods

After relations to some other Finite Element methods are discussed, an attemp to visualize the inter-connection between the different formulations can be seen in Figure 1.2.2. The formulations withinthis diagram are classified by (in order)

• the set of unknowns [green]

• the differential operator acting on the Finite Element shape functions [blue]

• the constraints imposed on the formulation (explicitly, through the space or through lagrangemultipliers)[red]

• spaces (continuous/discontinous/normal-continuous) for u (and σ) [black]

Standard

mixed

hybrid

primal dual primal DG dual DG

CG FEM u∇u [[u]]=0u is cont.

M FEM (u,σ)div(σ) [[σ]]n=0u disc.,σn cont.

HM FEM (u,σ,uF )div(σ) [[σ]]n=0u, σ discont.

DG FEM u∇u [[σ]]n=[[u]]=0u discont.

DG FEM (u,σ)∇u [[σ]]n=[[u]]=0u, σ discont.

DG FEM (u,σ)div(σ) [[σ]]n=[[u]]=0u, σ discont.

HDG FEM (u,σ,uF)div(σ) [[σ]]n=[[uF]]=0u, σ discont.

HDG FEM (u,uF)∇u [[σ]]n=[[uF ]]=0u discont.

HDG FEM (u,σn,F)∇u [[σn,F ]]=[[u]]=0u discont.

-

?

@@R

?

-

⋃ ⋃

abbreviations:C: continuousD: discontinuousM: mixedH: hybridG: Galerkin

Figure 1.2.2: Interconnection of several Finite Element methods

This diagram shows the relations of the proposed HDG method (bold) to other methods. Whereaswe went the direct way from CG FEM to our formulation [turquoise path] when introducing theformulation, we also pointed out the relation to other methods in the last section. Therefore wepresented the HDG method as

1.2.2.1 a hybridized mixed method with a modification of the conservation constraint ([[σ]]n = 0 →[[σ]] = 0). We went the way from CG FEM to HDG FEM in dual DG form over mixed FEM (MFEM) and hybridized mixed FEM (HM FEM) [purple path].

1.2.2.2 a Standard DG Formulation without additional unknowns uF . It was shown that HDG FEM isa subset of DG FEM [dark red] by eliminating the lagrange multiplier uF and presenting theequivalent Standard DG formulation.

1.2.2.3 a stabilized DG formulation. Here we showed another equivalent representation of the proposedmethod which uses additional unknowns σn,F for the flux - instead of the variable uF - on thefacets [orange].


1.2.3 A priori error analysis

For a precise analysis, we have to define appropriate spaces to work with. In contrast to CG our discretevariational formulation includes normal derivates which are not defined for all function in H1(Ω). So wehave to use additional regularity. We assume our problem has a solution u ∈ H1(Ω) ∩H2(Th) =: W .Remember that we already defined V = W × L2(Fh), s.t. our bilinearform Bh : V × V → R iswell-defined.For our a priori error estimates we will first show consistency and adjoint consistency of our method,which will later on give us Galerkin orthogonality on the one hand and will allow us for using Aubin-Nitsche arguments to get error estimates in the L2-norm on the other hand. Then we’ll show coercivityand boundedness of the bilinearform with the respect to an HDG norm. These results will then beused for the error estimates, but they already imply uniquely solvability of the discrete problem itself.The last ingredient will be the approximation property of the chosen Finite Element space in the HDGnorm. Then we can put all those results together to show optimal a priori error estimates for the HDGmethod.Before we start with the analysis lets summarize the introduced spaces:

W = H1(Ω) ∩H2(Th) (1.2.35a)V =

(u, uF ) : u ∈ H2(Th) ∩H1(Ω), uF ∈ L2(Fh)

= W × L2(Fh) (1.2.35b)

V kh =

(u, uF ) : u ∈ Pk(T ) ∀ T ∈ Th, uF ∈ Pk(E) ∀ E ∈ Fh

= Qk

h × F kh ⊂ V (1.2.35c)

Additional subscripts 0 and D indicate homogeneous and nonhomogeneous Dirichlet conditions whichare incorporated into the space in an L2-projective sense. Attention: Vh,0 ⊂ V0 but Vh,D * VD, i.e.L2-projected Dirichlet data may not coincide with exact Dirichlet data.The superscripted k which describes the used (uniform) polynomial order of the space, will often bedropped. During our analysis we always assume that the solution u of (1.2.1) lies in W and ε iselement-wise continuous. For simplicity we furthermore assume that we have only Dirichlet boundaryconditions on the whole boundary ∂Ω. Modifications for Neumann and Robin boundary conditions(even without Dirichlet boundary conditions) are trivial.In order to make u ∈ W compatible with functions in V we introduce the identity operator

idh : W → V, idhu = (u, tr|Fh(u))

The subindex h indicates that the identity operator depends on the mesh as Fh does. We will use thisidentity operator at several points implicitly by writing u for idhu. Only when additional emphasis onthe different spaces seems necessary we will use this operator explicitly.At several places we will work with a “characteristic“ size of an element, which will be denoted byh. As for well shaped elements different choices of h are possible and those differ only by a constantof order one, we won’t determine the choice of this element size h further. E.g., think of it as thediameter of the element.

Remark 1.2.5 (Variable order):Throughout the analysis in this work we assume uniform polynomial order, but nevertheless, thederivation of modified results is straight forward.

1.2.3.1 Consistency

Our aim here is to show that our Bilinearform is consistent, such that Galerkin-Orthogonality holds.So we remember that tr|E(u) is single-valued for u in W , s.t. [[u]] = 0 . Let u be the solution of


(1.2.1), then for any v = (v, vF ) ∈ V0 we get with the definition of the bilinearform (1.2.8):

Bh(u,v) =∑T∈Th


∫∂Tε∂u

∂n(v − vF ) ds −

∫∂T

[[u]]︸︷︷︸=0

ε (∂v∂n− τh [[v]]) ds

=∑T∈Th

∫T−div(ε ∇u) v dx+

∫∂Tε∂u

∂nvF ds

=∫

Ωf v dx+

∑E∈Fint

h

∫E

[[ε ∇u]]n︸︷︷︸=0

vF ds+∑

E∈Fexth

∫Eε∂u

∂nvF︸︷︷︸=0

ds

(1.2.36)

In the last step we used that the integrals from both sides of an interior facet can be summed up withinone integral. Furthermore we made use of the relation [[ε ∇u]] = 0, which holds true since u is thesolution of the PDE (1.2.1). As the problem and the bilinearform are symmetric, we also directly getconsistency of the adjoint problem, which will allow us for using Aubin-Nitsche arguments for L2-normerror estimates later:

Bh(u,v) = Bh(v,u) =∫

Ωf v dx ∀ v ∈ V0 (1.2.37)

Proposition 1.2.1 (Galerkin Orthogonality). Let uh ∈ Vh,D be the solution of (1.2.9) and u be thesolution of (1.2.1). Then there holds

Bh(uh − u,vh) = 0 ∀ vh ∈ Vh,0 (1.2.38)

Proof. Due to (1.2.36) there holds

Bh(uh − u,vh) = Bh(uh,vh)− Bh(u,vh) =∫

Ωf vh dx−

∫Ωf vh dx = 0 ∀ vh ∈ Vh,0 (1.2.39)

Note that the restriction to Dirichlet boundary condition is not at all necessary here, but was onlychosen for simplicity.

1.2.3.2 Stability and boundedness

We assume that ε and τh are constant on each element and facet respectively. Due to the non-conformity of our method, i.e. due to the fact that solutions of the discrete problem are not in H1(Ω),we can not work with the H1(Ω) (semi-) norm to show stability and boundedness of the introducedbilinearform. Instead, we have to define new norms. An appropriate discrete norm for the brokenSobolev space H1(Th), which is the natural one to show coercivity in, is

|||u|||2ε,∗ :=∑T∈Th

ε‖∇u‖2

T + ε

h[[[u]]]2∂T

(1.2.40)

where [[[·]]]∂T := ‖[[·]]‖∂T .This is actually a semi-norm similar to the H1 semi-norm for conforming Finite Elements. But as thebroken Sobolev space H1(Th) does not contain any continuity restrictions on element interfaces thejump terms appear in the norm.


The kernel of the norm operator is the space of constant functions. If we now pose Dirichlet conditionson a (nonzero) part of the boundary or constrain an integral average of u

∫G u dx = C ∈ R for G some

nonzero part of the domain the kernel is only the trivial one and the proposed operator is actually anorm8.In this norm we can show coercivity of Bh : Vh × Vh → R. But as we also have to show boundednessof Bh for v ∈ Vh and u ∈ V we have to consider another norm to bound normal derivative terms byit. Therefore we introduce the second HDG norm

|||u|||2ε :=∑T∈Th

ε‖∇u‖2

T + ε

h[[[u]]]2∂T + εh‖∂u

∂n‖2∂T

(1.2.41)

which is the natural one to boundedness of the bilinearform in.For u ∈ Vh the norms |||u|||ε and |||u|||ε,∗ are equivalent. As u ∈ Vh is piecewise polynomial we canbound the difference of both norms by |||u|||2ε,∗ with the inverse inequality (see (A.2.3b)):

|||u|||2ε − |||u|||2ε,∗ =∑T∈Th

εh‖∂u∂n‖2∂T ≤

∑T∈Th

εc2T‖∇u‖2

T ≤ c2sr|||u|||2ε,∗ (1.2.42)

with csr = maxT∈Th

cT ∈ R only depending on the shape and not on the size of the elements. Obviouslythen there holds:

1√1 + c2

sr

|||u|||ε ≤ |||u|||ε,∗ ≤ |||u|||ε ∀ u ∈ Vh (1.2.43)

We will denote each element contribution by ||| · |||T and ||| · |||T .

CoercivityWith respect to those norms we can turn over to coercivity and stability. Let’s start with coercivity:

Proposition 1.2.2 (Coercivity). For a shape regular mesh and τhh (with h the local mesh size )sufficiently large Bh(·, ·) is coercive on Vh with respect to the norm ||| · |||ε, that is

Bh(u,u) ≥ c|||u|||2ε,∗ ≥ αBh|||u|||2ε ∀ u ∈ Vh (1.2.44)

with c, αBh ∈ R independent of the mesh size.

Proof. We will make use of several elementary inequalities which are listed in section A.2 of theappendix. These are Cauchy Schwarz inequalities, Young’s inequality and inverse inequalities. With

8If neither is true and ΓD = ∅ and ΓR 6= ∅ add εh‖u‖2ΓR to |||u|||2ε,∗ and you have a norm. All following results can

also be adapted very easily to this modified norm


them we can directly prove the Proposition. For all u ∈ Vh there holds

Bh(u,u) =∑T∈Th

ε ‖∇u‖2T − 2

(∂u

∂n, [[u]]

)∂T

+ τh[[[u]]]2∂T

(A.2.1a)≥

∑T∈Th

ε ‖∇u‖2T − 2

∥∥∥∥∥∂u∂n∥∥∥∥∥∂T

[[[u]]]∂T + τh[[[u]]]2∂T

(A.2.3b)≥

∑T∈Th

ε ‖∇u‖2T − ‖∇u‖T

(2 cT√

h[[[u]]]∂T

)+ τh[[[u]]]2∂T

(A.2.2)≥

∑T∈Th

ε ‖∇u‖2T − 1

2‖∇u‖2T −

12

(2cT√h

[[[u]]]∂T)2

+ τh[[[u]]]2∂T

=∑T∈Th

ε 12‖∇u‖

2T +

(τh −

2c2T

h

)[[[u]]]2∂T

≥ min(12 , τhh− 2c2

sr) |||u|||2ε,∗(1.2.43)≥ 1√

1 + c2sr

min(12 , τhh− 2c2

sr) |||u|||2ε

≥ 12√

1 + c2sr

|||u|||2ε

(1.2.45)

where in the last step we used τh ≥2c2sr+ 1

2h

. So with c = 12 and αBh = 1

2√

1+c2sr

the claim holdstrue.

Coercivity implies uniqueness of the discrete solution as for two discrete solutions u1,u2 ∈ Vh,D thereholds

0 = Bh(u1 − u2,u1 − u2) ≥ C|||u1 − u2|||ε

and as ||| · |||ε is a norm for Vh,0 and u1−u2 is zero on dirichlet boundaries u1 and u2 have to coincide.

Remark 1.2.6 (Stabilization parameter):As we see in the appendix the constant of the inverse inequality csr or locally cT depends linearly on thepolynomial degree k. So with the condition τh >

2c2sr+ 1

2h

, it is sufficient to fulfill τh > α0hk2 for an α0

only depending on the shape. It is of course sufficient if those relations are fulfilled locally (csr → cT ).

BoundednessTo show boundedness of the bilinearform in the ||| · |||ε-norm is even easier, as we will see now:

Proposition 1.2.3 (Boundedness). For all u,v ∈ V there holds:

|Bh(u,v)| ≤ βBh |||u|||ε |||v|||ε

with βBh = supx∈Fh(1 + τhh)

Proof. We can prove the claim with a combination of different versions of the Cauchy-Schwarz in-


equality:

|Bh(u,v)|=

∣∣∣∣ ∑T∈Th

ε

(∇u,∇v)T − (∂u∂n, [[v]])∂T − (∂v

∂n, [[u]])∂T + τh([[u]], [[v]])∂T

∣∣∣∣(A.2.1a)≤

∑T∈Th

ε‖∇u‖T‖∇v‖T +

∥∥∥∥∥∂u∂n∥∥∥∥∥∂T

[[[v]]]∂T +∥∥∥∥∥∂v∂n

∥∥∥∥∥∂T

[[[u]]]∂T + τh[[[u]]]∂T [[[v]]]∂T

(A.2.1b)≤

∑T∈Th

ε(‖∇u‖2

T + h

∥∥∥∥∥∂u∂n∥∥∥∥∥

2

∂T

+ 1h

[[[u]]]2∂T + τh[[[u]]]2∂T) 1

2

·(‖∇v‖2

T + h

∥∥∥∥∥∂v∂n∥∥∥∥∥

2

∂T

+ 1h

[[[v]]]2∂T + τh[[[v]]]2∂T) 1

2

(A.2.1b)≤

∑T∈Th

ε‖∇u‖2

T + h

∥∥∥∥∥∂u∂n∥∥∥∥∥

2

∂T

+ 1h

[[[u]]]2∂T + τh[[[u]]]2∂T 1

2

·

∑T∈Th

ε‖∇v‖2

T + h

∥∥∥∥∥∂v∂n∥∥∥∥∥

2

∂T

+ 1h

[[[v]]]2∂T + τh[[[v]]]2∂T 1

2

=|||u|||2ε +

∑T∈Th

ετh[[[u]]]2∂T

12

·

|||v|||2ε +∑T∈Th

ετh[[[v]]]2∂T

12

max≤ sup

x∈Fh(1 + τhh) |||u|||ε |||v|||ε

(1.2.46)

Remark 1.2.7 (non-constant coefficients):Notice that the assumption of piecewise constant coefficients is not at all necessary for proving sta-bility or boundedness. If you interprete ε‖∇u‖2

T and εh[[[u]]]2∂T as short notations for (ε∇u,∇u)T and

( εh[[[u]]], [[[u]]])∂T you directly get both inequalities in the more general case. The only difference makes

the inverse inequality (used in the step from line two to line three in equation (1.2.45)). Here, youwould have to adjust the constant cT to the properties you request on the coefficient ε.

1.2.3.3 Approximation

Let’s restate one version of the famous Bramble Hilbert Lemma

Lemma 1.2.4 (Bramble Hilbert Lemma). Let U be some Hilbert space, and L : Hm → U , m = k+1be a continuous linear operator such that Lq = 0 for all polynomials q ∈ Pk. Then there holds

‖Lv‖U ≤ C|v|Hm ∀ v ∈ Hm

with C independent of u and the size of the domain.


Now we introduce the element local L2-projector Πh,k : Hm(T ) → Pk(T )× F kh (T ) which is defined

by the following relations

Πh,k(u) :=(Πelh,k(u), Πfac

h,k (u))

with(Πelh,k(u), v

)T

= (u, v)T ∀ v ∈ Pk(T )(Πfach,k (u), vF

)∂T

=(tr|Fh(u), vF

)∂T

∀ vF ∈ F kh (T )Pk(E)

(1.2.47)

With those preparations we can now state and prove the following Lemma.

Lemma 1.2.5. For u ∈ Hm(T ), m = k + 1, k ≥ 1 and the L2-projector Πh,k : Hm(T ) →Pk(T )× F k

h (T ) of (1.2.47), there holds

|||(idh −Πh,k)u|||Tε ≤√‖ε‖∞ C |u|Hm(T ) ∀ u ∈ Hm(T ) (1.2.48)

with C independent of u and element size.

Proof. We use the Bramble Hilbert Lemma with U = H2(T )× L2(∂T ) and the scalar product

(u,v)U :=∫T∇u ∇v dx+

∫∂T

1h

[[u]][[v]] ds+∫∂Th∂u

∂n

∂v

∂nds+ 1

h2

(∫∂Tu ds

)(∫∂Tv ds

)︸︷︷︸

=:R(u,v)

such that (u,u)U = ‖u‖2U = (||| 1√

εu|||Tε )2 +R(u,v). The operator L which is chosen as

L := idh −Πh,k

is clearly linear and for all polynomials q ∈ Pk there holds Lq = 0. We still have to show that L is acontinuous operator. Note that for u ∈ H2(T ) the trace of u and the trace of the normal derivative arecontinuous linear operators such that ‖u‖U = ||| 1√

εu|||Tε +

√R(u,u) ≤ C‖u‖H2 holds. Furthermore the

introduced L2-projector (w,wF ) = Πh,ku of a function u onto the polynomial space Pk(T )× F kh (T )

can be bounded in the ‖ · ‖U -norm by the H1(T )-norm of w and the L2(∂T )-norm of wF .As a consequence of both, the U -norm of Lu can be bounded in the H2(T )-norm of u as well

‖Lu‖U ≤ ‖idhu‖U + ‖Πh,ku‖U ≤ C‖u‖H2(T ) (1.2.49)

So all requirements of the Bramble Hilbert Lemma are fulfilled and we obtain

|||(idh −Πh,k)u|||Tε ≤ ‖√εLu‖U ≤

√‖ε‖∞ C |u|Hm(T ) ∀ u ∈ Hm(T ) (1.2.50)

Now we will look at the global interpolation error to bound the approximation error infvh∈Vh |||u−vh|||ε.Therefore we will transform all element contributions on a reference element. There we will use Lemma1.2.5 and by transforming back we will get the correct h scaling of the interpolation error.


Proposition 1.2.6. On a shape regular mesh Th which consists of affine transformed elements and afunction u ∈ H1(Ω) ∩Hm(Th), m = k + 1 ≥ 2 there holds the following approximation result for V k

h

as introduced in (1.2.35c)

infvh∈V kh

|||u− vh|||2ε ≤ C∑T∈Th‖ε‖∞ h2k

T |u|2Hm(T ) (1.2.51)

with C independent of the mesh size. For quasi-uniform meshes the direct conclusion of it is

infvh∈V kh

|||u− vh|||ε ≤ C√‖ε‖∞ hk |u|Hm(Th) (1.2.52)

Proof. Let’s first note that as u ∈ H1(Ω) the trace of u onto facets is single-valued, such that the L2-projected approximation is also single-valued. So we are able to define the global projection-operatorΠΩh,k : Hm(Th)→ Vh as an element-wise projection with ΠT

h,k the L2-projector of (1.2.47) with respectto a single element T . So we start with

infvh∈V kh

|||u− vh|||2ε ≤ |||u−ΠΩh,ku|||2ε =

∑T∈Th

(|||(idTh −ΠT

h,k)u|||Tε)2

(1.2.53a)

where idThu is the restriction of idhu to the element T . Now we can transform each element T to anaccording reference element T of size O(1), where the transformation from reference element T tophysical element T is defined by FT : T → T .

T

T

TF

Figure 1.2.3: Affine transformation FT from reference domain T to the physical domain T

On each element we then get(|||(idTh −ΠT

h,k)u|||Tε)2 hdTh

−2T

(|||(idTh −ΠT

h,k)u FT |||Tε)2 (1.2.53b)

where we used that the scaling of |||v|||ε is hd−2T , which becomes clear when we look at the three

contributions of the norm separately:

part of ||| · |||2ε scaling factor integral trafo. diff. op. scaling overall scaling∫Tε(∇v)2 dx 1 ' hdT ' h−2

T ' hd−2T

1h

∫∂Tε[[v]]2 ds h−1

T ' hd−1T 1 ' hd−2

T

h∫∂Tε∂v

∂n

2ds hT ' hd−1

T ' h−2T ' hd−2

T


Here, we can use another property of the L2-projection: it commutes with the transformation, i.e. theresult of the L2-projection of a function transformed from the reference to the physical domain andthe transformation from reference to physical domain of the L2-projection of a function are the same.So there also holds:

(|||(idTh −ΠT

h,k)u|||Tε)2 hdT h

−2T

(|||(idTh −ΠT

h,k)(u FT )|||Tε)2

(1.2.53c)

We finally reached the point where Lemma 1.2.5 can be applied on each element to get(|||(idTh −ΠT

h,k)u|||Tε)2 ‖ε‖∞hdT h

−2T |u FT |2Hm(T ) (1.2.53d)

Now we transform back to physical domain and achieve(|||(idTh −ΠT

h,k)u|||Tε)2 ‖ε‖∞hdT h

−2T h−dT h2m

T |u|2Hm(T ) = ‖ε‖∞h2kT |u|2Hm(T ) (1.2.53e)

This proves the first part of the proposition. If we consider quasi-uniform meshes we get h ' hT ∀ T ∈Th and can move the factor h2k

T in front of the sum, what proves the second part of the proposition.

Remark 1.2.8 (Boundary conditions):We did not use special treatment for boundary conditions as we assume Dirichlet boundary conditionsthat are incorporated into the space are either exact or approximations in an L2-orthogonal sense. Inboth cases Vh can be exchanged by Vh,D and the results stay uneffected. If neither is true, the changesin the result are small as long as the error of approximation is of the same magnitute as the remainder.

1.2.3.4 Putting it all together

HDG norm estimatesFirst we are interested in error estimates in the HDG norm of the form

|||u− uh|||ε ≤ Chq|u|Hm(Th) (1.2.54)

with q depending on the chosen polynomial approximation and the regularity of the weak problem.Before we can show this kind of estimate we need the following Lemma:

Lemma 1.2.7 (Modified Cea’s Lemma). Let uh ∈ V kh be the solution of (1.2.9) and u ∈ H1(Ω) ∩

Hm(Th) ⊂ W , m ≥ 2 the solution of (1.2.1). Then the error is only bounded a constant away fromthe approximation error:

|||u− uh|||ε ≤ C infvh∈Vh,D

|||u− vh|||ε (1.2.55)

for C ∈ R only depending on τh and the shape regularity.

Proof. We start with a triangle inequality to separate an expression that is later on going to be thediscrete error and one that is going to be the approximation error

|||u− uh|||2ε ≤ |||u− vh|||2ε + |||vh − uh|||2ε ∀ vh ∈ Vh,D (1.2.56)


From section 1.2.3.3 we have the estimates for the approximation error. Now we need estimates forthe other part |||vh−uh|||ε, which we easily get from coercivity, Galerkin orthogonality and boundednessof the bilinearform:

|||vh − uh|||2ε(1.2.45)≤ 1

αBhBh(vh − uh,vh − uh︸︷︷︸

wh∈Vh,0

)

(1.2.38)= 1αBh

Bh(vh − u,vh − uh) + Bh(u− uh,wh)︸︷︷︸=0

(1.2.46)≤ βBh

αBh|||vh − u|||ε · |||vh − uh|||ε

(1.2.57)

Thus we obtain|||vh − uh|||ε ≤

βBhαBh|||u− vh|||ε (1.2.58)

and with (1.2.56) the discrete error depends (except for a constant factor) only on the approximationerror :

|||u− uh|||2ε ≤ (1 + βBhαBh

)|||u− vh|||2ε (1.2.59)

As the choice of vh ∈ Vh,0 is arbitrary, the same result also holds true for the infimum:

|||u− uh|||2ε ≤ (1 + βBhαBh

) infvh∈Vh,D

|||u− vh|||2ε (1.2.60)

Note that the constant C = 1 + βBhαBh

only depends on the constants of the recent propositions s.t.

1 + βBhαBh

= 1 + supx∈Fh(1+τhh2 )

√1 + c2

sr.

Now error estimates in the HDG norm are easy to get:

Proposition 1.2.8 (HDG norm estimate). Let Th be a quasi-uniform shape regular mesh, uh ∈ V kh

be the solution of (1.2.9) and u ∈ H1(Ω) ∩ Hm(Th) ⊂ W , m ≥ 2 the solution of (1.2.1). Thenthere holds the following error estimate:

|||u− uh|||ε ≤ C√‖ε‖∞h

s|u|Hm(Th) s = min(k,m− 1) (1.2.61)

with C ∈ R only depending on the shape regularity, and the choice of the stabilization parameter τh

Proof. The claim follows with Proposition 1.2.6 and Lemma 1.2.7.

L2 norm estimatesAs approximation theory predicts, an approximation of hk+1 in the L2-Norm is possible if we assumethe problem to be sufficiently smooth. That this can actually be achieved with the proposed HDGmethod, we will show next with the help of the Aubin Nitsche trick.


Lemma 1.2.9 (Aubin Nitsche trick). Let uh ∈ V kh be the solution of (1.2.9) and u ∈ H1(Ω) ∩

H2(Th) ⊂ W the solution of (1.2.1), both with homogeneous Dirichlet boundary conditions. Thenthe error converges one order faster in the L2-norm than in the HDG norm for quasi-uniform meshes(h ' hT ∀ T ∈ T )

‖u− uh‖L2 ≤ c√‖ε‖∞ h|||u− uh|||ε (1.2.62)

Proof. We consider the adjoint problem to (1.2.1)

−div(ε∇w) = u− uh in Ω (1.2.63)

with homogeneous Dirichlet boundary conditions. u−uh is in VD⊕VD,h ⊂ L2(Ω) and with a sufficientlysmooth boundary ∂Ω we obtain9 w ∈ H2(Ω) ⊂ V . As our method is also adjoint consistent (see,1.2.37 ), there holds

Bh(v,w) = (u− uh, v)Ω ∀ v ∈ V0 (1.2.64)Testing with v = u− uh gives

‖u− uh‖2L2 = (u− uh, u− uh)Ω = Bh(u− uh,w) (1.2.65)

From (1.2.38) we know that Bh(u − uh,vh) = 0 for every vh ∈ Vh,0. If we choose vh to be the(discontinuous) linear L2-projection of w, vh = Πh,1w, there holds

‖u−uh‖2L2 = Bh(u−uh,w) = Bh(u−uh,w−Πh,1w)

(1.2.46)≤ βBh|||u−uh|||ε|||w−Πh,1w|||ε (1.2.66)

Now we can use standard results for the interpolation error |||w − Πh,1w|||ε ≤ C√‖ε‖

∞h|w|H2 and

elliptic regularity (for sufficiently smooth boundaries) |w|H2 ≤ c‖u− uh‖L2 to get

‖u− uh‖L2 ≤ βBh c C︸︷︷︸=c

√‖ε‖∞h|||u− uh|||ε (1.2.67)

Remark 1.2.9 (nonhomogeneous Dirichlet boundary conditions):The results also hold true for nonhomogeneous Dirichlet boundary conditions but the proof is moreinvolved for the more general case in general u− uh /∈ V0.

And finally we obtain an error estimate in the L2-Norm:

Proposition 1.2.10 (L2 norm estimate). Let Th be a quasi-uniform shape regular mesh, uh ∈ V kh be

the solution of (1.2.9) and u ∈ H1(Ω) ∩Hm(Th) ⊂ W , m ≥ 2 the solution of (1.2.1). Then thereholds the following error estimate:

‖u− uh‖L2(Ω) ≤ C ‖ε‖∞ hs+1|u|Hm(Th) s = min(k,m− 1) (1.2.68)

with C ∈ R only depending on the shape regularity, and the choice of the stabilization parameter τh


9see for example [Evan98, 6.3]


1.2.4 Conservation property of the Hybrid DG method

One important property of most DG methods and also the proposed HDG method is the conservationof a discrete flux σh. As pointed out when deriving the Hybrid DG method as a modified mixed methodthe conservation property appears as an explicit constraint in the variational formulation. So thereholds:

Proposition 1.2.11. The proposed Hybrid DG method is conservative (locally and globally). Fur-thermore we can determine a single-valued discrete flux σh·n on each facet.

Proof. If we test the bilinearform Bh with (χT , 0), where χT is the characteristic function of oneelement, we get

Bh(u, (χT , 0)) =∫∂Tε

(−∂u∂n

+ τh [[u]])

1 ds =∫Tf dx (1.2.69)

With respect to the numerical flux σh·n := ε(−∂u∂n

+ τh[[u]])we have that the total flux which enters

the element is equal to the sum of internal sources. Using this flux we have a locally conservativeformulation.To show that we don’t lose anything when adding all elements together we have to show that thenumerical flux σh·n is single valued on each facet, s.t.

σ+h n

+ + σ−h n− = 0 (1.2.70)

for all inner facets. We show this facet by facet. Let E be an inner facet of the mesh. We choosev = (0, vF ), where vF has only support on E. Now if we plug this into the bilinearform, there holds:

Bh(u, v) =∑T∈Th

∫∂Tε

(∂u

∂n−τh [[u]]

)vF ds = 0 (1.2.71)

As vF is chosen to have only support on E we can rewrite this as

Bh(u, v) =∫E

(σ+h n

+ + σ−h n−)vF ds = 0 (1.2.72)

If the numerical flux σh·n is at most of order k where vF ∈ V kh , equation (1.2.70) holds in the strong

sense10. Otherwise it is still fulfilled in a weak sense. So the formulation is globally (strongly/weakly)conservative.

1.2.5 Numerical Example: Two dimensional incompressible potential flow

Problem description

We consider an irrotational incompressible flow field v = (vx, vy) around a circular disk. Due toirrotationality there holds curl(v) = 0 and so there exists a vector potential ψ with ∇ψ = (vy,−vx).Furthermore there holds div(v) = div(∇ψ) = −∆ψ = 0 as the fluid is incompressible. We will solvefor the flow potential, which is often called the stream function as the isolines of it coincide withthe stream lines of the flow field. The geometry is a rectangular (−1, 1)× (−1

2 ,12) without a disk of

10The assumption is true as long as ε is element-piecewise constant.


diameter 12 which is removed in the middle. We prescribe Dirichlet boundaries on ∂[(−1, 1)× (−1

2 ,12)]

and homogeneous ones on the disk (this corresponds to ∇ψ×n = v ·n = 0). Inhomogeneous Dirichletboundaries are taken from the exact solution

ψ(x, y) = y(x2 + y2 −R2)x2 + y2 R = 1

4

0.5 1

2

(a) Setup of the problem: red boundaries are inhom.Dirichlet boundaries, blue boundaries are homogeneousones

(b) Initial mesh (Level 0) and discrete solution (p = 4)

Figure 1.2.4: Setup, initial mesh and solution of Numerical Example 1.2.5

On the choice of τh

First we want to investigate the sensitivity of the error with respect to the stability parameter. Moti-vated by the estimates for the inverse inequality for simplices quoted in A.2, we choose

τh = α(p+ 1)(p+ d)

d

|∂T ||T |

, d = 2

As our mesh uses affine transformed as well as curved simplices we replace |∂T | by the absolute valueof the jacobian determinant of the transformation of the facet and |T | by the absolute value of thejacobian determinant of the transformation of the element at each integration point.

0.004

0.005

0.006

0.007

0.008

0.009

0.01

1 10 100 1000

L2-e

rror

alpha

HDGDGCG

best approx.

(a) error-dependency on α for p = 1, L = 0

3e-05

4e-05

5e-05

6e-05

7e-05

8e-05

9e-05

0.0001

1 10 100 1000

L2-e

rror

alpha

HDGDGCG

best approx.

(b) error-dependency on α for p = 4, L = 0

Figure 1.2.5: Sensitivity of the error on α

From the analysis we know that α has to be chosen sufficiently large for stability but not too large inorder not to raise the constant of the boundedness estimate. We claimed that for well shaped meshes


this parameter is of order one. In Figure 1.2.5 the error is plotted for two choices of the polynomialorder. We see that α = 2 is a good choice and that the behaviour does not differ much betweenDG11 and HDG. Furthermore the proposed scaling of τh motivated by affine transformed simplicesworks quiet well also for the curved elements. The approximation of a CG method on the same meshwith the same polynomial degree (this means a lot less degrees of freedom) is comparable with an αwhich has been chosen too large. That behaviour is expected as α penalizes the discontinuities andfor α→∞ a continuous approximation is approached.We conclude that α is a parameter which you can choose to be 2 or a bit higher for less regular meshes,but you don’t have to ”optimize“ the parameter for every new problem, unless you want to reduceeven the least few percentages of the discrete error.Notice that the solution is smooth and so the comparison to CG might not be as persuasive as in othercases. But this is not the intension of this example.

Convergence of the method

Table 1.1 confirms the L2 norm estimates of the a priori error analysis. For p = 1 we get a convergencerates of 2 and for p = 4 the rate of 5.

meshlevel unknowns ‖u− uh‖L2 unknowns ‖u− uh‖L2

(elements) p DOF error order p DOF error order0 (54) 1 354 4.601e-3 — 4 1 290 4.141e-5 —1 (216) 1 356 1.516e-3 1.602 5 010 3.790e-6 3.4502 (864) 5 304 4.028e-4 1.912 19 740 1.398e-7 4.7613 (3 456) 20 976 1.032e-4 1.965 78 360 4.996e-9 4.8064 (13 824) 83 424 2.605e-5 1.986 312 240 1.652e-10 4.9185 (55 296) 332 736 6.533e-6 1.995 1 246 560 5.253e-12 4.9756 (221 184) 1 329 024 1.636e-6 1.998

Table 1.1: L2-Error-Convergence of the HDG method for the potential flow problem

Comparison to Standard symmetric interior penalty DG method

We now want to present a comparison of DG and HDG. Both are expected to have with same orderof convergence and for this numerical example the expectations are fulfilled. However, we now wantto observe how the error decreases with respect to the total number of degrees of freedom andas a second measure with respect to the nonzero entries that are generated in the matrix. Bothmeasures discriminate the HDG formulation, as it can be made more efficient with static condensation.Nevertheless we will see, that both measures already indicate that the HDG method without any furtheroptimization is comparable and even cheaper than the DG method.On the same mesh, using the same polynomial degree, the DG and the HDG formulation should havesimilar L2 norm errors. Thus DG yields smaller errors for the same number of degrees of freedoms.But the gap between both methods is not dramatically large. Furthermore if we consider the secondmeasurement, the nonzero entries in the matrix, we see in Figure 1.2.6 that the curves of DG andHDG are very close and for higher polynomial degree HDG performs even better than DG.

11We used the symmetric interior penalty DG formulation


Thus even without taking advantage of the matrix properties (which will be discussed at the end ofthis chapter)the HDG method is comparable or even preferable to the DG method.Furthermore the DOF-plot shows the expected convergence rates O(hp+1) = O(N− p+1

2 ) for bothmethods.

1e-11

1e-10

1e-09

1e-08

1e-07

1e-06

1e-05

0.0001

0.001

1000 10000 100000 1e+06

L2-e

rror

degrees of freedom

HDG (p=1)DG (p=1)

HDG (p=4)DG (p=4)

(a) error over degrees of freedom

1e-11

1e-10

1e-09

1e-08

1e-07

1e-06

1e-05

0.0001

0.001

10000 100000 1e+06 1e+07

L2-e

rror

nonzero matrix entries

HDG (p=1)DG (p=1)

HDG (p=4)DG (p=4)

(b) error over nonzero matrix entries

Figure 1.2.6: Possible (incomplete) measurements of the computational effort


1.3 The Hybrid Discontinuous Galerkin Method for (linear)hyperbolic problems


In this section we consider the other limit of the convection diffusion equation, the linear transportequation, which is of hyperbolic type:

div(bu) = f in Ωu = uD on Γin

(1.3.1)

where b has continuous normal-component bn on element interfaces.This time we only consider two types of boundary conditions. These are inflow and outflow boundaryconditions. Boundary conditions on characteristic boundaries (these are boundaries where bn is zero)don’t influence the solution in the domain and are therefore ignored.Linear hyperbolic problems can be solved very conveniently by Discontinuous Galerkin method. Herewe use the Standard DG upwind formulation as a starting point for our Hybrid DG method. As welater on want to join the numerical scheme with the HDG formulation for the laplace operator, withoutlosing the desired features, we have to modify the Standard DG upwind formulation a little.The Standard DG upwind formulation for (1.3.1) is derived by partial integration on each element andchoosing the upwind value for the element boundary integral:

∫Ωdiv(bu)v dx =

∑T∈Th

∫Tdiv(bu)v dx DG≈

∑T∈Th

−∫Tbu∇v dx+

∫∂Tbnu

DG v ds

=: CDGh (u, v)

where uDG =u bn > 0unb bn < 0 with unb the discrete trace on the neighbour element. This means that

we always choose the value which lies in the direction from where the convection originates.Now we modify the formulation such that the access to the neighbour element functions is replacedby an access to the facet functions which we already introduced in section 1.2. So on outflow parts ofthe element boundary we choose the element value uup = u (see Fig. 1.3.1(a)). But on inflow parts ofthe element boundary we choose, in contrast to Standard DG methods, the facet value uup = uF (seeFig. 1.3.1(b)). This has the advantage that we preserve the property of the elliptic hdg formulation,where element degrees of freedom from different elements don’t couple directly.

∑T∈Th

∫Tdiv(bu)v dx HDG≈

∑T∈Th

−∫Tbu∇v dx+

∫∂Tbnu

up v ds

with uup =u bn > 0uF bn < 0

(1.3.2)

1.3. HYBRID DG FOR (LINEAR) HYPERBOLIC PROBLEMS 37

(a) at the outflow of an element (green, surroundedby dashed line) the element value (red, hatched) istaken. uup=u

(b) at the inflow of an element (green, surroundedby dashed line) the facet value (red, bold) is taken.uup=uF

Figure 1.3.1: Choice for the upwind value for the HDG formulation for one element

It is obvious that now the unknowns of different elements don’t even couple at all, because facetunknowns just couple with one neighbouring element, the downwind element, and not with upwindelement. We overcome this issue by adding a constraint which glues the facet values on the trace ofthe upwind element (in a weak sense):∑

T∈Th

∫∂Tout

bn(uF − u) vF ds = 0 with ∂Tout := x ∈ ∂T, bn(x) > 0 (1.3.3)

The “stabilization“ is illustrated in Fig. 1.3.2 and results in uF = uup

Figure 1.3.2: at the outflow the facet value and the upwind element value are glued together

Adding up both equations (1.3.3) and (1.3.2) we end up with

Ch(u,v) :=∑T∈Th

−∫Tbu∇v dx+

∫∂Tbnu

up v ds+∫∂Tout

bn(uF − u) vF ds

=∫

Ωf v dx (1.3.4a)

You can rewrite this formulation also in the following form which is sometimes more convenient for theanalysis by doing partial integration once more and noting that the first boundary term vanishes on


∂Tout as uup is exactly u there and bn is positive on ∂Tout and negative on ∂Tin := x ∈ ∂T, bn(x) < 0

Ch(u,v) =∑T∈Th

∫Tdiv(bu)v dx+

∫∂Tin|bn|(u− uF ) v ds+

∫∂Tout|bn|(uF − u) vF ds

(1.3.4b)

In this form you especially see that the chosen sign for the stabilization term results in positive entrieson the matrix of the discretized linear system.A third way of presenting the bilinearform, which will be helpful when looking at the discrete adjointof Ch, can be achieved if we move uFvF from ∂Tout to ∂Tin, or, more generally spoken (to considerthe boundary facet in a correct way), we first substract

∫∂Tout|bn|uFvF ds and add

∫∂Tin|bn|uFvF ds

on each element (∗). For inner facets both terms cancel out and the result is a simple shift from thein- to the outflow side. On boundary facets we have to add additional boundary integrals for balance.Starting from (1.3.4b), after integration by parts and the described procedure, we can formulate Chalso in the following way:

Ch(u,v) =∑T∈Th

−∫Tbu∇v dx+

∫∂Tin|bn|(−uF ) v ds+

∫∂Tout|bn|(u v + (uF − u) vF ) ds

(∗)=

∑T∈Th

−∫Tbu∇v dx+

∫∂Tbnu

up [[v]] ds+∫∂T∩Γ

bnuFvF ds

(1.3.4c)All three formulations are algebraically equivalent with Standard DG upwind, but also exhibit the samenice properties of HDG methods that we already saw for the elliptic case and will also be elaboratedlater on:

• element-wise assembly is possible

• element unknowns of different elements don’t couple directly

• inner degrees of freedom can be eliminated, i.e. the local problems are uniquely solvable

The discrete problem reads:

Find u ∈ VD, such that for all v ∈ Vh,0 :

Ch(u,v) = 〈f, v〉 =∫

Ωf v dx

(1.3.5)

Remark 1.3.1 (Discontinuities):Equation (1.3.1) could have discontinuous solutions. Those can be captured exactly by the method aslong as the mesh is aligned to the discontinuity. That’s because for a mesh aligned to the discontinuitywe have bn = 0 and so the problems on each side of the discontinuity decouple also for the numericalscheme. Notice that as we assume that at least bn is continuous and so is zero on both sides of adiscontinuity, the unknowns and testfunctions uF and vF disappear completely from the formulation.In order to make the matrix of the discretized bilinearform regular again, we could add integrals ofthe form

∫∂T ∗ |b|[[u]]vF ds on element boundaries ∂T ∗ = x ∈ ∂T, bn(x) = 0 which are aligned to

characteristics. Those integrals wouldn’t influence the solution on the neighbouring domains, but giveconvenient defining equations for uF which for equal order interpolation and single-valued b on eachfacet give uF = u. Notice that this choice for uF is arbitrary as there is no meaningful value todescribe at the discontinuity.



After the method is introduced, we want to carry out an a priori error analysis for the HDG upwindformulation. Although we already pointed out that the proposed formulation is equivalent to theStandard DG method, the analysis will be carried out completely and in a compatible way that allowsfor unification with the discretization for the elliptic operator.

Basic assumptions

We assume div(b) = 0, i.e. b a solenoidal vector. To pay tribute to the hyperbolic nature of (1.3.1) wereplace the assumption u ∈ W of section 1.2 by u ∈ W conv = v ∈ L2(Ω); b · ∇v ∈ L2(Ω) instead.That means that we, in contrast to the assumption u ∈ W , allow for discontinuities along streamlines.For ease of notation we will also call it W during the whole analysis of the pure hyperbolic problem.Further we also use V := W×L2(Fh). The discrete spaces are not exchanged. Similar to the notationof the pure elliptic case VD incorporates Dirichlet boundary conditions on the boundary. It only makessense to prescribe Dirichlet boundary conditions on Γin = x ∈ Γ, bn < 0 or Γout = x ∈ Γ, bn > 0.We assume ΓD = Γin to follow the physical causality of the equation.We further assume that discontinuities are not aligned to the mesh, meaning that bn = 0 doesn’t holdon a single facet. As mentioned in Remark 1.3.1 we could also consider this case easily, but to avoida blow-up in notation we forego that.In this section we want to devide the bilinearform into two parts

Ch(·, ·) = Celh (·, ·) + Cfach (·, ·) (1.3.6)

where Celh (·, ·) includes all integrals on the volume of each element and Cfach (·, ·) includes all elementboundary (facet) integrals.Remember also that we assumed div(b) = 0.

1.3.2.1 Consistency

We now expand the true solution of (1.3.1) by the facet values uF (u) and with the regularity alongstreamlines of V conv we are able to define the upwind trace of the true solution uF (x) = lim

ε+→0u(x−εb)

on each facet as long as bn 6= 0. 12 We use the representation (1.3.4b) to show consistency. Letu ∈ VD be the solution of (1.3.1), then for all v ∈ V

Ch(u,v) =∑T∈Th

∫Tdiv(bu)︸︷︷︸

=f

v dx+∫∂Tin|bn|(u− uF )︸︷︷︸

=0

v ds+∫∂Tout|bn|(uF − u)︸︷︷︸

=0

vF ds

=∫

Ωf v dx

(1.3.7)holds. We again directly obtain Galerkin orthogonality from that.

Proposition 1.3.1 (Galerkin Orthogonality). Let uh ∈ Vh,D be the solution of (1.3.5) and u be thesolution of (1.3.1). Then there holds

Ch(uh − u,vh) = 0 ∀ vh ∈ Vh (1.3.8)12In the case bn = 0 the definition is irrelevant because of the decoupling of the values on the neighbouring elements


Proof. Using (1.3.7) we get: For all vh ∈ Vh

Ch(u− uh,vh) =∫

Ωf vh dx−

∫Ωf vh dx = 0 (1.3.9)

We are also interested in the adjoint consistency of our formulation. Therefore we pose the adjointproblem to (1.3.1):

b · ∇u = f in Ω ; u = 0 on Γout = Γin; with b = −b (1.3.10)

Note that the boundaries without tilde still refer to the original problem. The outflow boundary of theoriginal problem Γout is the inflow boundary Γin of the adjoint problem. In the case div(b) = 0 theadjoint problem is the same as the original one except that the convective velocity is b = −b now andthe role of in- and outflow boundary are exchanged. After discretization (using representation (1.3.4c))we conclude that the discrete adjoint Ch(v,u) is consistent with the adjoint problem (1.3.10).

Ch(v,u) =∑T∈Th

∫T

(b · ∇u)︸︷︷︸f

v dx−∫∂Tbnu

up [[v]]︸︷︷︸=0

ds−∫∂T∩Γ

bn vFuF︸︷︷︸=0 on Γ

ds

(1.3.11)

where vFuF = 0 on Γ holds since uF is zero on the inflow boundary of the adjoint problem which isΓin = Γout and vF is assumed to be zero on the outflow boundary of the adjoint problem which isΓin = Γout. So we directly get

Ch(v,u) =∫

Ωfv dx ∀ v ∈ V0 (1.3.12)

for the solution of the adjoint problem.

1.3.2.2 Stability

In this paragraph we assume u = 0 on Γin = ΓD which is reasonable as we will need the stability resultfor uhu− vh with uh,uh in Vh,D and thus uhu− vh ∈ Vh,0. Here we follow the line of arguments of[ES09] to show stability of the (hybrid) DG method in a way which allows for generalization with theelliptic formulation later on (see section 1.2.3.2).After introducing the convection HDG norm

|||u|||2C,∗ :=∑T∈Th

∫T

h

|b|(b · ∇u)2 dx+

∫∂T|bn|[[u]]2 ds+

∫∂T∩Γout|bn|u2

F ds

=

∑T∈Th

(|||u|||TC,∗)2 (1.3.13)

we claim the stability statement, known as inf-sup stability

Proposition 1.3.2 (inf-sup stability). Let u be the exact solution of (1.3.1). Then, for a constantαCh ∈ R independent of the meshsize h the following is true:

supv∈Vh,0

Ch(u,v)|||v|||C,∗

≥ αCh|||u|||C,∗ ∀ u ∈ Vh,0 (1.3.14)

Furthermore, for every u ∈ Vh,0 we can find a function v ∈ Vh,0, s.t.

Ch(u, v) ≥ αCh|||u|||C,∗|||v|||C,∗ (1.3.15)


Proof. We will show that the inf-sup-condition holds with a particular choice for v ∈ Vh,0:

v = (v, vF ) = c1(u, uF ) + c2( h

|b|b · ∇u, 0) (1.3.16)

With this approach we will show the second part of the Proposition which implies the first one.First we want to present the terms coming up for the choices v = u and v = ( h|b|b · ∇u, 0) separatelyand afterwards pick one possible linear combination for v.Notice that as we have assumed that div(b) = 0 there holds:∫

Tb · ∇u u dx =

∫Tdiv(bu) u dx p.I.= −

∫Tb · ∇u u dx+

∫∂Tbnu

2 ds

⇒∫Tdiv(bu) u dx = 1

2

∫∂Tbnu

2 ds(1.3.17)

Now we use the representation (1.3.4b) for our bilinearform and test with v = u:

Ch(u,u) =∑T∈Th

∫Tdiv(bu)u dx+

∫∂Tin|bn|(u− uF ) u ds+

∫∂Tout|bn|(uF − u) uF ds

(1.3.17)=

∑T∈Th

12

∫∂Tbnu

2 ds+∫∂Tin|bn|(u− uF ) u ds+

∫∂Tout|bn|(uF − u) uF ds

=

∑T∈Th

12

∫∂Tin|bn|(u2 − 2 u uF ) ds+ 1

2

∫∂Tout|bn|(u2 − 2u uF + 2u2

F ) ds

(∗)=∑T∈Th

12

∫∂T|bn|[[u]]2 ds− 1

2

∫∂T∩Γin

|bn|2uF︸︷︷︸=0

ds+ 12


F ds

= 1

2∑T∈Th

∫∂T|bn|[[u]]2 ds+


F ds

(1.3.18)In (∗) we moved one half of u2

F from ∂Tout to ∂Tin, or more generally spoken (to consider the boundaryfacets in a correct way) we first substract 1

2∫∂Tout|bn|u2

F ds and add 12∫∂Tin|bn|u2

F ds on each element.For inner facets both terms cancel out and the result is a simple shift from the in- to the outflow side.On boundary facets we have to add additional boundary integrals for balance.So we see that we can use this for the second part of the norm ||| · |||C,∗. The volume part of that normwill be account for with v = v∗= ( h|b|b · ∇u, 0). Notice that as v∗F = 0 there holds v∗ ∈ Vh,0:

Ch(u,v∗) =∑T∈Th

∫T

h

|b|div(bu)b · ∇u dx+

∫∂Tin

h

|b||bn| [[u]] b · ∇u ds

≥∑T∈Th

∫T

h

|b|(b · ∇u)2 dx−

∫∂Tin

h

|b||bn| |[[u]]| |b · ∇u| ds

(A.2.2)≥

∑T∈Th

∫T

h

|b|(b · ∇u)2 dx− α

∫∂Tin|bn| [[u]]2 ds− 1

α

∫∂Tin

h2

|b|2|bn| (b · ∇u)2 ds

|bn|<|b|≥

∑T∈Th

∫T

h

|b|(b · ∇u)2 dx− α

∫∂Tin|bn| [[u]]2 ds− 1

α

∫∂Tin

h2

|b|(b · ∇u)2 ds

(A.2.3a)≥

∑T∈Th

∫T

h

|b|(b · ∇u)2 dx− α

∫∂Tin|bn| [[u]]2 ds− c2

T

α

∫T

h

|b|(b · ∇u)2 ds

≥∑T∈Th

(1− c2

T

α

)∫T

h

|b|(b · ∇u)2 dx− α

∫∂Tin|bn| [[u]]2 ds

(1.3.19)


Set α to 2c2sr s.t. (1− c2

T

α) is greater 1

2 and add sufficiently much of u to ensure that the total boundarycontribution has a positive factor greater 1

2 . This is achieved with the choice v = (2c2sr + 1)u + v∗.

The constant csr is the maximum on the element local constants cT ∈ R which only depend on theshape regularity of the mesh. With this linear combination we get:

Ch(u, v) ≥ 12 |||u|||

2C,∗ (1.3.20)

Now we additionally need|||v|||C,∗ ≤

12αCh

|||u|||C,∗ (1.3.21)

From scaling arguments, i.e. arguments that bound the norm of gradients of polynomials by the normof the polynomial itself and an appropriate scaling 1

h, you can bound |||v∗|||C,∗ by c|||u|||C,∗. Then with

the triangle inequality there holds

|||v|||2C,∗ ≤ (2c2sr + 1)2|||u|||2C,∗ + |||v∗|||2C,∗

≤ ((2c2sr + 1)2 + c2)︸︷︷︸=:( 1

2αCh)2

|||u|||2C,∗ (1.3.22)

And finallyCh(u, v) ≥ αCh|||u|||C,∗|||v|||C,∗ (1.3.23)

Notice that β ∈ R only depends on fixed constants and on the shape regularity of the mesh.

1.3.2.3 Boundedness

Now we want to show boundedness of the bilinearform Ch. As we have seen in chapter 1.2 there isnot only one possible choice for an HDG norm and which one might be the more natural one maydiffer from boundedness- to coercivity-estimates. When dealing with the boundedness it is convenientto work with a second convection HDG norm ||| · |||C :

|||u|||2C :=∑T∈Th

∫T

|b|hu2 dx+

∫∂T|bn|(uup)2 ds+


F ds

=

∑T∈Th

(|||u|||TC)2 (1.3.24)

We can easily show with scaling arguments that for all discrete functions u ∈ Vh

|||u|||2C,∗ ≤ C|||u|||2C (1.3.25)

holds.In the elliptic case we also introduced two different norms but there it would have been sufficient towork with the latest one only. In the hyperbolic case, here, we have to work with both norms at thesame time. That’s why both norms are involved in the following Proposition:

Proposition 1.3.3. For all u ∈ V and v ∈ V0 there holds:

|Ch(u,v)| ≤ |||u|||C |||v|||C,∗


Proof. We want to bound Ch(u,v) for u ∈ V and v ∈ V0. Therefore we start with Celh and Cfach ofthe representation (1.3.4c):

|Celh (u,v)| =

∣∣∣∣∣∣∑T∈Th−∫Tbu∇v dx

∣∣∣∣∣∣≤

∑T∈Th

∫T

∣∣∣|b| 12h− 12u∣∣∣ ∣∣∣|b|− 1

2h12 (b · ∇v)

∣∣∣ dxC.S≤

∑T∈Th

√∫T|b|h−1u2 dx

√∫T|b|−1h(b · ∇v)2 dx

(1.3.26)

|Cfach (u,v)| (1.3.4c)=

∣∣∣∣∣∣∑T∈Th

∫∂Tbnu

up [[v]] ds−∫∂T∩Γin

|bn|uF vF︸︷︷︸=0

ds+∫∂T∩Γout

|bn|uFvF ds

∣∣∣∣∣∣

C.S≤

∑T∈Th

√∫∂T|bn|(uup)2 ds

√∫∂T|bn|[[v]]2 ds+

√∫∂T∩Γout|bn|u2

F ds

√∫∂T∩Γout|bn|v2

F ds

(1.3.27)

With these preparations we can now finish the proof of boundedness:

|Ch(u,v)| ≤ |Celh (u,v)|+ |Cfach (u,v)|C.S.≤

∑T∈Th|||u|||TC · |||v|||TC,∗

≤ |||u|||C · |||v|||C,∗

(1.3.28)


As the HDG norms of the convection problem scale differently from those of the elliptic problem, werepeat and adapt the basic steps of 1.2.3.3. We will again assume some Hm(Th)-regularity which iscertainly to much to ask for in most cases. A more appropriate space would be the Hilbert spacewith L2-functions with weak directional derivatives up to mth order. Then for piecewise constant b anadapted version of the Bramble Hilbert Lemma could be stated13. The remainder would still be thesame.Later on, when looking at the convection diffusion problem for ε > 0 the full Hm(Th)-regularityassumption is again fulfilled for sufficiently smooth boundaries, but, due to nearly discontinuities (e.g.boundary layers), the upcoming constants (e.g. the Hm-semi norm) in the following estimates maydegerenate.We will stick with the more convenient although less appropriate Hm(Th)-regularity assumption anduse another application of the original Bramble Hilbert Lemma.

Lemma 1.3.4. For u ∈ Hm(T ), m = k + 1 and the L2-projector Πh,k : Hm(T )→ Pk(T )× F kh (T )

of (1.2.47), there holds

|||(idh −Πh,k)u|||TC ≤ C√‖b‖∞ |u|Hm(T ) ∀ u ∈ Hm(T ) (1.3.29)

13in this case we need piecewise constant b, s.t. the mth order seminorm has the kernel Pm−1


Proof. The proof goes analogously to the proof of Lemma 1.2.5. The scalar product is replaced by ascalar product which fulfills (u,u)U = ‖u‖2

U ≥(|||√

1|b| u|||TC

)2+ R(u,v) with R(·, ·) some symmetric

nonnegative bilinearform. Again boundedness of the operator L = idh−Πh,k is easy to get and so allrequirements for the Bramble Hilbert Lemma are fulfilled again and we get

|||(idh −Πh,k)u|||TC ≤ ‖√|b| Lu‖U ≤ C

√‖b‖∞ |u|Hm(T ) ∀ u ∈ Hm(T ) (1.3.30)

Now we will look at the global interpolation error to bound the approximation error infvh∈Vh |||u−vh|||C .Therefore we will transform all element contributions on a reference element. There we will use Lemma(1.3.4) and after transforming back to physical domain, we will get the right scaling of the interpolationerror.


h


infvh∈V kh

|||u− vh|||2C ≤ C∑T∈Th‖b‖∞h

2k+1T |u|2Hm(T ) (1.3.31)


infvh∈V kh

|||u− vh|||C ≤ C√‖b‖∞ hk+ 1

2 |u|Hm(Th) (1.3.32)

Proof. Again the proof is similar to the one of Proposition 1.2.6 except for the scaling of the normwhich is different to the diffusion HDG norms of section 1.2.3. On each element we now get

(|||(idTh −ΠT

h,k)u|||TC)2 hdTh

−1T

(|||(idTh −ΠT

h,k)u FT |||TC)2 (1.3.33)

where the scaling of the square of the norm is hd−1T , which will be explained in the following table

part of ||| · |||2C scaling factor integral trafo. diff. op. scaling overall scaling1h

∫T|b|u2 dx ' h−1

T ' hdT 1 ' hd−1T∫

∂T|bn|(uup)2 ds 1 ' hd−1

T 1 ' hd−1T∫

∂T∩Γout|bn|u2

F ds 1 ' hd−1T 1 ' hd−1

T

So the last step of the proof of Proposition 1.2.6 reads(|||(idTh −ΠT

h,k)u|||TC)2 ‖b‖∞hdT h

−1T h−dT h2m

T |u|2Hm(T ) = ‖b‖∞h2k+1T |u|2Hm(T ) (1.3.34)

which proves the claim.



HDG norm estimates

Lemma 1.3.6 (Modified Cea’s Lemma). Let uh ∈ V kh be the solution of (1.3.5) and u ∈ H1(Ω) ∩

Hm(Th) ⊂ W , m ≥ 2 the solution of (1.3.1). Then the error in the ||| · |||C,∗-norm is bounded only aconstant away from the approximation error in the ||| · |||C-norm:

|||u− uh|||C,∗ ≤ C infvh∈Vh,D

|||u− vh|||C (1.3.35)

for C ∈ R only depending on the shape regularity.

Proof. First of all we can show that the difference of a function vh ∈ Vh to the discrete solutionmeasured in the ||| · |||C,∗-norm can be bounded by the difference to the exact solution in the ||| · |||C-norm

|||vh − uh|||C,∗ ≤1αCh|||vh − u|||C (1.3.36)

This holds with inf-sup stability, consistency and boundedness of Ch (v ∈ Vh,0 is again as in Prop.1.3.2):

|||vh − uh|||C,∗ · |||v|||C,∗ ≤1αChCh(vh − uh, v)

= 1αCh

(Ch(vh − u, v) + Ch(u− uh, v)︸︷︷︸

=0

)≤ 1

αCh|||vh − u|||C · |||v|||C,∗

(1.3.37)

So after separating the discrete error and the approximation error and using (1.3.36) and (1.3.25) weget

|||u− uh|||2C,∗ ≤ |||vh − uh|||2C,∗ + |||vh − u|||2C,∗ ≤ ( 1α2Ch

+ C2)|||vh − u|||2C (1.3.38)

and so the claim holds for the constant C =√

1α2Ch

+ C2 with C2 the constant of (1.3.25).

As mentioned before, Hm(Th) is normally to much to demand from a pure hyperbolic problem. How-ever, it is the easiest way to get error estimates in the HDG norm:


be the solution of (1.3.5) and u ∈ H1(Ω)∩Hm(Th) ⊂ W , m ≥ 2 the solution of (1.3.1). Then thereholds the following error estimate:

|||u− uh|||C,∗ ≤ C√‖b‖∞ hs|u|Hm(Th) s = min(k + 1

2 ,m−12) (1.3.39)

with C ∈ R only depending on the shape regularity.



Remark 1.3.2 (L2 norm estimates):Although in principal some kind of Aubin Nitsche technique is possible to gain convergence results inthe L2 norm, the requirements on the adjoint problem are far from realistic, at least when working withHm(Th). So instead of a proof of convergence in the L2 norm, it should just be mentioned that under”good conditions“ the L2-norm is expected to convergence one half order faster than the ||| · |||C,∗-norm.However, there exist examples where the convergence rate hk+ 1

2 in the L2 norm is the best you canget even for smooth solutions.

1.3.3 Numerical Example: A steady rotating convection problem

This example is taken and modified from [NPC09]. We consider the domain Ω = (0, 1)2 \ Γ with

Γ = 12 × (0, 1

2) and the convective velocity b =(y − 1

212 − x

)which is divergence-free. On all inflow

boundaries of (0, 1)2 Dirichlet boundary conditions u = 0 and u = 1 − tanh(10(1 − 4y)) on the left(inflow) part of Γ are prescribed.

b

(a) Setup of the problem: red bound-aries are inflow boundaries, blue bound-aries are outflow boundaries

(b) Initial mesh (Level 0) and discretesolution (p = 5)

Figure 1.3.3: Setup, initial mesh and solution of Numerical Example 1.3.3

The Dirichlet boundary condition should be transported from the left to the right side of Γ.Due to non-aligned meshes the discrete solution will be diffusive on coarse meshes, i.e. the profile onthe right side of the slit will be smeared out a little bit. For different spatial and polynomial resolutionswe plotted the exact solution, the prescribed profile left of the slit and the discrete solution of theprofile right of the slit in Figure 1.3.4. As our solution is smooth, higher polynomial degree reducesthis smearing effect faster as higher grid resolution. We also see in the following table that optimalconvergence rates in the L2-norm are achieved, what indicates that for smooth functions the resultsof our analysis are often pessimistic. Notice that as gradients of the solution are steep, high resolutionis necessary before we observe the asymptotically correct convergence rate. See e.g. that for p = 0the resolution after six uniform refinement steps is still too poor to get order 1 convergence.


meshlevel unknowns ‖u− uh‖L2 unknowns ‖u− uh‖L2

(elements) p DOF error order p DOF error order0 (231) 0 603 0.3813 — 1 1 437 0.1231 —1 (924) 2 361 0.2648 0.526 5 646 0.07563 0.7032 (3 696) 9 342 0.1956 0.437 22 380 0.02667 1.5043 (14 784) 37 164 0.1352 0.533 89 112 7.560e-3 1.8194 (59 136) 148 248 0.0894 0.598 355 632 1.539e-3 2.2965 (236 544) 592 176 0.0558 0.678 1 420 896 2.861e-4 2.4286 (946 176) 2 367 072 0.0329 0.7640 (231) 3 3 798 0.0252 — 5 7 083 4.610e-3 —1 (924) 14 988 8.517e-3 1.568 28 026 8.396e-4 2.4572 (3 696) 59 544 6.981e-4 3.609 111 492 3.747e-5 4.4863 (14 784) 237 360 4.015e-5 4.120 444 744 5.985e-7 5.9684 (59 136) 947 808 2.483e-6 4.015 1 776 528 8.484e-9 6.1415 (236 544) 3 787 968 1.558e-7 3.994

Table 1.2: L2-Error-Convergence of the HDG method for the rotating convection problem

0

0.5

1

1.5

2

0.1 0.15 0.2 0.25 0.3 0.35 0.4

discrete end profilediscrete start profile

exact solution

(a) p = 0, mesh level 4

0

0.5

1

1.5

2

0.1 0.15 0.2 0.25 0.3 0.35 0.4


exact solution

(b) p = 1, mesh level 2

0

0.5

1

1.5

2

0.1 0.15 0.2 0.25 0.3 0.35 0.4


exact solution

(c) p = 3, mesh level 0

0

0.5

1

1.5

2

0.1 0.15 0.2 0.25 0.3 0.35 0.4


exact solution

(d) p = 5, mesh level 0

Figure 1.3.4: Both sides of the slit for different polynomial and grid resolutions compared to exactsolution


1.4 Adding diffusion and convection together

Subject of this section is now the steady version of the original model problem (1.1.1)div(−ε∇u+ bu) = f in Ω, div(b) = 0

u = uD on ΓD∂u∂n

= g on ΓN∂u∂n

+ βu = h on ΓR

(1.4.1)

with ε > 0 and ‖b‖ ≥ 0 everywhere.

1.4.1 Joining the limits

The discretization of the convection diffusion problem is easy as the approximation spaces of bothsubproblems are equal. The appropriate generalization of the proposed formulations is achieved byadding both bilinearforms together:

Ah(u,v) := Bh(u,v) + Ch(u,v) (1.4.2)

and the discrete problem now reads:

Find u ∈ VD, such that for all v ∈ Vh,0 :Ah(u,v) = 〈f, v〉

(1.4.3)

Remark 1.4.1 (Comparison to Standard DG):Notice that although in both limit cases the bilinearforms could be transformed to Standard DGbilinearform BDGh and CDGh respectively, the corresponding elimination of the facet unknowns in thegeneral case does not coincide with the sum of both Standard DG bilinearforms.

ADGh (u,v) 6= BDGh (u,v) + CDGh (u,v) (1.4.4)

This becomes immediately clear when you look at the result of the elimination of uF in the nextsection.

1.4.1.1 Characterization of the numerical trace

We motivated and derived convenient choices for the numerical flux in both limit cases ‖b| → 0 andε → 0. In the elliptic case we have seen that we use a central approximation which is stabilized bysome stabilization term depending on the discontinuity of the flux, whereas in the pure hyperboliccase we used the upwind choice which is known to be the natural and stable choice. When both,diffusion and convection are considered, the numerical trace lies in between those choices. Here wewant to briefly describe how the intermediate numerical trace looks like. If we again only test withfacet functions vF we can reconstruct an expression for the facet value uF :

uF =2τhεu − [[ε∂u

∂n]] + |bn|uDG

2τhε+ |bn|(1.4.5)

Thus the conserved total flux Jh is describes as

Jh := bnuup − ε∂u

∂n+ ετh(u− uF ) (1.4.6)

1.4. ADDING DIFFUSION AND CONVECTION TOGETHER 49


When adding diffusion to the convection, problems with the regularity vanish (at least from a puristpoint of view). The solution of the convection diffusion equation has the same regularity as the sameproblem without convection. Nevertheless the effect of ”nearly discontinuities“ may appear (e.g. atboundaries) which renders the apperent regularity useless in many cases.Consistency and boundedness in the a priori error analysis follow straight forward from the limitingproblems. For stability we have to put some extra technical effort into the proof. The approximationsection is in general easily reconstructed from the elliptic case, but some extra attention has to betaken as the regularity may degenerate.Again we assume to have Dirichlet data on the whole boundary.

1.4.2.1 HDG norms

The appropriate HDG norms for the convection diffusion problem also follow directly from the appro-priate norms of the subproblems.

|||u|||2∗ = |||u|||2C,∗ + |||u|||2ε,∗ |||u|||2 = |||u|||2C + |||u|||2ε (1.4.7)

In both limits the HDG norms coincide with the appropriate of the limit problem.

1.4.2.2 Consistency

Consistency follows directly from the consistency of the pure hyperbolic and the pure elliptic problem.Let u ∈ VD be the solution of (1.4.1), then for all v ∈ V0

Ah(u,v) =∑T∈Th

∫Tdiv(bu)v dx+

∫∂Tin|bn|(u− uF )︸︷︷︸

=0

v ds+∫∂Tout|bn|(uF − u)︸︷︷︸

=0

vF ds

+∫Tε ∇u ∇v dx−

∫∂Tε∂u

∂n(v − vF ) ds

=∑T∈Th

∫Tdiv(bu− ε∇u)︸︷︷︸

=f

v dx+∫∂Tε∂u

∂nvF ds

=∫

Ωf v dx+

∑E∈Fint

h

∫E

[[ε ∇u]]︸︷︷︸=0

vF ds+∑

E∈Fexth

∫Eε∂u

∂nvF︸︷︷︸=0

ds

(1.4.8)

Proposition 1.4.1 (Galerkin Orthogonality). Let uh ∈ Vh,D be the solution of (1.4.3)

Ah(uh,vh) = (f, vh)Ω ∀ vh ∈ Vh,0

and u be the solution of (1.4.1). Then there holds

Ah(uh − u,vh) = 0 ∀ vh ∈ Vh,0 (1.4.9)

Proof. Using (1.4.8) we have: For all vh ∈ Vh

Ah(u− uh,vh) =∫

Ωf vh dx−

∫Ωf vh dx = 0 (1.4.10)


1.4.2.3 Stability

Proposition 1.4.2 (inf-sup condition). Let u be the exact solution of (1.4.1). Then, for a constantαAh ∈ R independent of the meshsize h and τhh sufficiently large the following holds:

supv∈Vh,0

Ah(u,v)|||v|||∗

≥ αAh|||u|||∗ ∀ u ∈ Vh,0 (1.4.11)

Furthermore for every u ∈ Vh,0 we can find a v ∈ Vh,0, s.t.

Ah(u, v) ≥ αAh |||u|||∗|||v|||∗ (1.4.12)

Proof. Similar to the approach we used in pure convective case when showing inf-sup stability insection 1.3.2.2 we will test Ah(·, ·) with v = γu + v∗, with γ ≥ 2c2

sr + 12 . We can bound the

elliptic HDG norm of the last part v∗ = ( h|b|b · ∇u, 0) also with |||u|||ε,∗ when using scaling arguments(‖∇w‖T ≤ ck

h2‖w‖2T ∀ w ∈ Pk(T )) and an inverse inequality:

|||v∗|||2ε,∗ :=∑T∈Th

ε

∥∥∥∥∥ h|b|∇(b · ∇u)

∥∥∥∥∥2

T

+ 1h

∥∥∥∥∥ h|b|(b · ∇u)∥∥∥∥∥

2

∂T

≤

∑T∈Th

εh2 ‖∇(∇u)‖2

T + h ‖∇u‖2∂T

≤ (c+ c2

sr)∑T∈Th

ε ‖∇u‖2T ≤ (c+ c2

sr)|||u|||2ε,∗

(1.4.13)

Keeping that in mind and using (1.3.21), it follows

|||v|||∗ ≤ k|||u|||∗ (1.4.14)

where k depends linearly on c, the constant of the shape regularity, and on fixed (problem-independent)constants and γ. With the last inequality it is sufficient to show

Ah(u, v) ≥ K|||u|||2∗ for some constant K ∈ R (1.4.15)

The last ingredient to show this is to bound Bh(u,v∗) as the other arising terms were already boundedin the analysis of the pure hyperbolic and the pure elliptic case, respectively. The bound follows again


with Cauchy-Schwarz inequality and the compiled relations from above:

Bh(u,v∗) :=∑T∈Th

∫Tε ∇u ∇v∗ dx−

∫∂Tε∂u

∂n[[v∗]] ds

−∫∂Tε∂v∗

∂n[[u]] ds+

∫∂Tετh [[u]] [[v∗]] ds

C.S.≥ −

∑T∈Th

ε

‖∇u‖T‖∇v∗‖T +h 1

2

∥∥∥∥∥∂u∂n∥∥∥∥∥∂T

h−12 [[[v∗]]]∂T

+h 12

∥∥∥∥∥∂v∗∂n∥∥∥∥∥∂T

h−12 [[[u]]]∂T +τh[[[u]]]∂T [[[v∗]]]∂T

C.S.≥ −

∑T∈Th

ε

‖∇u‖2

T +h∥∥∥∥∥∂u∂n

∥∥∥∥∥2

∂T

+ 1h

[[[u]]]2∂T +τh[[[u]]]2∂T

12

‖∇v∗‖2T +h

∥∥∥∥∥∂v∗∂n∥∥∥∥∥

2

∂T

+ 1h

[[[v∗]]]2∂T +τh[[[v∗]]]2∂T

12

inv. ineq.≥ −

∑T∈Th

max(1 + c2T , 1 + τhh)|||u|||Tε,∗|||v∗|||Tε,∗

≥ −max(

1 + c2sr, sup

x∈Fh(1 + τhh)

)|||u|||ε,∗|||v∗|||ε,∗

(1.4.14)≥ −K2|||u|||2ε,∗ with K2 = max

(1 + c2

sr, supx∈Fh(1 + τhh))· k

(1.4.16)

Finally we are able to show the inf-sup condition. We choose the proposed test function but additionallyrequire that γ ≥ K2+ 1

2αBh

. Thereby we achieve the following:

Ah(u, v) = γ Bh(u,u)︸︷︷︸≥αBh |||u|||

2ε,∗

+ Bh(u,v∗)︸︷︷︸≥−K2|||u|||2ε,∗

+ Ch(u, v)︸︷︷︸≥ 1

2 |||u|||2C,∗

≥ 12 |||u|||

2∗

(1.4.14)≥ 1

2k︸︷︷︸αAh

|||u|||∗|||v|||∗ (1.4.17)

where we used coercivity of Bh and equations (1.3.20) and (1.4.16).

1.4.2.4 Boundedness

Boundedness (Continuity) follows directly from the continuity of the subproblems

Proposition 1.4.3. For all u,v ∈ V there holds:

|Ah(u,v)| ≤ βAh|||u|||C |||v|||C,∗

with βAh = βBh = supx∈Fh(1 + τhh)

Proof. We just have to use the boundedness results for Bh and Ch and directly get

Ah(u,v) = Bh(u,v) + Ch(u,v)≤ |||u|||ε,∗ |||v|||ε + supx∈Fh(1 + τhh) |||u|||C,∗ |||v|||CC.S.≤ supx∈Fh(1 + τhh) |||u|||∗ |||v|||

(1.4.18)



Plugging in the results of the pure convective and the pure elliptic case directly gives


h


infvh∈V kh

|||u− vh|||2 ≤∑T∈Th

(C1‖b‖∞h

2k+1T + C2‖ε‖∞h

2kT

)|u|2Hm(T ) (1.4.19)

with C1, C2 independent of the mesh size. For quasi-uniform meshes the direct conclusion of it is

infvh∈V kh

|||u− vh|||2 ≤(C1‖b‖∞h

2k+1 + C2‖ε‖∞h2k)|u|2Hm(Th) (1.4.20)

Proof. The claim follows with the definition of the HDG norm and the approximation results for thelimit cases of the convection diffusion equation.


HDG norm estimates

Lemma 1.4.5. Let uh ∈ V kh be the solution of (1.4.3) and u ∈ H1(Ω) ∩ Hm(Th) ⊂ W , m ≥ 2

the solution of (1.3.1). Then the error in the ||| · |||∗-norm is only bounded by a constant and theapproximation error in the ||| · |||-norm:

|||u− uh|||∗ ≤ C infvh∈Vh,D

|||u− vh||| (1.4.21)

for C ∈ R only depending on the shape regularity and the stabilization parameter τh.

Proof. The proof works exactly the same way as the one of Lemma 1.3.6, just exchange Ch by Ahand adapt the constants.

Now it is reasonable to assume Hm(Th), at least for a small mesh Peclet number Peh = |b|hε

and wecan work with the following error estimate:


be the solution of (1.4.3) and u ∈ H1(Ω) ∩ Hm(Th) ⊂ W , m ≥ 2 the solution of (1.4.1). Thenthere holds the following error estimate:

|||u− uh|||C,∗ ≤ (C1

√‖b‖∞

√h+ C2

√‖ε‖∞) hs|u|Hm(Th) s = min(k,m− 1) (1.4.22)

with C ∈ R only depending on the shape regularity.



Remark 1.4.2 (L2 norm estimates):In general L2 norm estimates are as difficult to get as in the pure convective limit. But if we assumesmall Peclét numbers, that is convection is not dominant anywhere, we can also get optimal orderL2 norm estimates. The easiest way to see that, is to remember (1.3.18) and thus the convectionbilinearform Ch does not destroy the coercivity of Ah in the ||| · |||ε norm. Additionally we have to adda boundedness estimate for Ah which uses ||| · |||ε norms only. This estimate has order one constantsas long as the Peclét numbers are small. In this framework we could easily adapt the results of thepure elliptic case to come to similar results and thereby to optimal order L2 norm a priori estimatesalso for the convection diffusion equation. Nevertheless, the convection dominated case is often themore interesting and more challenging case. The analysis carried out in this section is still valid in theconvection dominated case. When we turn over to Navier-Stokes equations or more specifically to theOseen problem, things get more complicated and we will use the kind of analysis which makes use ofthe coercivity, instead of the one we applied in this chapter.

Remark 1.4.3 (A posteriori error estimates):Although a posteriori error estimators are not being discussed here, due to the similarity to StandardDG methods, adapted versions of the error estimators proposed in [SZ09] would be possible.

1.4.3 Numerical Example: A steady convection-dominated problem

We take the numerical test of [ES09, Example 1] where the square (0, 1)2 with homogeneous dirichletboundary conditions is considered. The flow field is constant b = (b1, b2)T , the diffusion coefficient εis constant and the source term f is set to be:

f = b1

(y − eb2y/ε − 1

eb2/ε − 1

)+ b2

(x− eb1y/ε − 1

eb1/ε − 1

)(1.4.23)

Then, for ε > 0, the exact solution to the convection-diffusion problem (1.4.1) is

u(x, y) =(x− eb1x/ε − 1

eb1/ε − 1

)(y − eb2y/ε − 1

eb2/ε − 1

)(1.4.24)

For our numerical tests we take ε = 0.01 and b = (2, 1)T . The solution has boundary layers at theoutflow boundaries, i.e. the top and right side (cp. Figure 1.4.1(a)). We want to compare the proposedmethod with the results given in [ES09, Example 1], i.e. with the streamline diffusion method and themixed Hybrid DG method of [ES09]. For our Hybrid DG method we again choose

τh = α(p+ 1)(p+ d)

d

|∂T ||T |

with α = 2.Notice that p = 0 is not possible for our proposed method. For p = 1 and p = 2 the results for ourmethod are quite similar to the ones of the mixed HDG method. For all method we can not expectoptimal convergence rates on the first refinement steps as the boundary layers are not yet resolved.Nevertheless, we can observe some important properties: In contrast to both HDG methods, thestreamline diffusion method introduces large layers which stem from the stabilization of the convectiveterm, so most of the error is owed to the overdiffusive stabilization. Because of that the gap betweendiscrete solutions and best approximation(see Table 1.4 left) increases for higher k, whereas both HDGmethods are very close to their L2 best approximation(see Table 1.4 right).

14V ph is the H1-conforming finite element with piecewise complete polynomials of degree p


(a) Solution of the convection diffusion problem(1.4.1) for ε = 0.01 and b = (2, 1)T

(b) Initial mesh (Level L = 0)

Figure 1.4.1: Initial mesh(L = 0) and solution of Numerical Example 1.4.3

streamline diffusion mixed HDG ([ES09]) HDG (α = 2)L h (elements) p = 1 p = 2 p = 3 p = 0 p = 1 p = 2 p = 1 p = 20 1/4 (32) 0.21 0.19 0.19 0.065 0.040 0.033 0.040 0.0341 1/8 (128) 0.15 0.13 0.13 0.048 0.036 0.025 0.035 0.0252 1/16 (512) 0.097 0.088 0.088 0.040 0.026 0.014 0.025 0.0143 1/32 (2 048) 0.056 0.050 0.050 0.032 0.014 0.0052 0.014 0.0054

Table 1.3: L2 errors of u

1.5 Computational aspects

There are several advantages which arise with the additional facet degrees of freedom. One is theelimination of inner degrees of freedom which is possible as long as the local problems are uniquelysolvable. We will show that this is the case in section 1.5.1. In section 1.5.2 the sparsity patterns ofDG and HDG methods will be discussed.

1.5. COMPUTATIONAL ASPECTS 55

conforming14 discontinuousminvh∈V ph ‖u− vh‖L2 minvh∈Qph ‖u− vh‖L2

L h (elements) p = 1 p = 2 p = 3 p = 0 p = 1 p = 20 1/4 (32) 0.18 0.10 0.063 0.057 0.040 0.0331 1/8 (128) 0.12 0.060 0.032 0.045 0.034 0.0242 1/16 (512) 0.073 0.029 0.012 0.038 0.024 0.0133 1/32 (2 048) 0.038 0.011 0.0030 0.030 0.013 0.0051

Table 1.4: L2 best approximation error

1.5.1 Elimination of inner degrees of freedom

Now we want to briefly comment on the elimination of the element degrees of freedom. If we define

A( u , v ) =∑T∈Th

∫Tε∇u∇v dx+

∫Tdiv(bu)v dx

−∫∂Tε

(∂u

∂nv + ∂v

∂nu− τhuv

)ds+

∫∂Tout

bnuv ds

B( uF, v ) =∑T∈Th

∫∂Tε

(∂v

∂nuF − τhuFv

)ds+

∫∂Tin

bnuFv ds

C( u , vF) =∑T∈Th

∫∂Tε

(∂u

∂nvF − τhuvF

)ds−

∫∂Tout

bnuvF ds

D(uF, vF) =∑T∈Th

∫∂TετhuFvF ds+

∫∂Tout

bnuFvF ds

(1.5.1)

we end up with the following system for u and uF(A BC D

)(uuF

)=(fg

)(1.5.2)

Notice that A is block diagonal and so you can invert A in an element by element fashion. Building

up the schur complement (which in the elliptic case is s.p.d.15 as(A BC D

)is s.p.d.) is easy:

u = A−1(f −BuF )⇒ (D − CA−1B)uF = g − CA−1f (1.5.3)

The solution uF can then be used to solve for u with Au = f −BuF element by element.

Remark 1.5.1 (Existence of local solvers):To calculate A−1 we need the local element-wise problems to be uniquely solvable, that is, we have toshow

A(u, u) = 0 + homogeneous b.c. (in a dg sense) ⇒ u ≡ 0

We can do that by showing coercivity of the bilinearform A. Therefore techniques similar to thosewe already used to show inf-sup stability of the global bilinearform can be used. Actually solvingAu = f(uF ) for a given uF is like solving the PDE on the small domain T where dirichlet boundaryconditions are incorporated in a way, which is known as the Nitsche approach, i.e. in a weak sense.

15symmetric positive definite


Remark 1.5.2 (Alternative HDG formulation):This property does not hold for the alternative HDG formulation presented in (see section 1.2.2.3).There we’d have to additionally prescribe some additional moment (e.g. average value) for eachelement to get locally solvable subproblems. Then again we would have to add those moments to theglobal unknowns.

1.5.2 Sparsity Pattern of HDG compared to DG methods

As we have already seen in the numerical example for HDG in the elliptic case in section 1.2.5, theefficiency of the method depends on the criteria we evaluate it with. On the one hand HDG methods,in comparison to other DG methods, come with additional degrees of freedoms, but on the other handthe stencil of the formulation is reduced and so the coupling of degrees of freedom is significantlyreduced. If we also take a look at the sparsity pattern for example 1.4.3 in Figure 1.5.1(a) we observewhat should already be clear from the construction of the method: element DOF don’t couple directlywith each other and facet DOF don’t couple direclty with other facet DOF, s.t. the matrix blocks Aand D of the last section are block diagonal. We can then invert each block, which is possible forA as well as for D and build according schur complements. The natural approach would be to firsteliminate the element DOF (compare section 1.5.1) resulting in a linear system of equations for thefacet DOF only. The corresponding sparsity pattern for the schur complement of the matrix is shownin Figure 1.5.1(a) is shown in Figure 1.5.1(b). We observe that the size of the schur complement issignificantly smaller. The difference between the size of the large matrix and the size of the schurcomplement increases for higher polynomial degree.

(a) HDG Sparsity Pattern, Matrix has the size 2728×2728

(b) condensated HDG Sparsity Pattern, Matrix has thesize 616× 616

Figure 1.5.1: Sparsity Pattern of HDG and DG for the problem 1.4.3 and the initial mesh (see alsoFigure 1.4.1(b)) for polynomial degree p = 10. Dirichlet degrees of freedom are also included here.

Another possibility is to eliminate the facet DOF which results in a schur complement which hasexactly the same sparsity pattern as standard dg method. This is not surprising as we already saw that

1.5. COMPUTATIONAL ASPECTS 57

a Standard DG formulation of the hdg method can be achieved by eliminating the facet functions.The corresponding sparsity pattern can be seen in Figure 1.5.2.We conclude that, concerning the facet functions, direct solvers as well as iterative solvers benefitmore from them than they suffer from a few more degrees of freedom.As long as we have to solve linear systems, the hybridized version of a DG methods should be prefered.Nevertheless for time-dependent problems which are solved explicitly, s.t. no linear system has to besolved, the additional degrees of freedoms on facet are an unneccessary burden.

Figure 1.5.2: Sparsity Pattern of DG for the problem 1.4.3 and the initial mesh (see also Figure1.4.1(b)) for polynomial degree p = 10 (Matrix has the size 2112× 2112 ).

58 CHAPTER 2. DIVERGENCE-FREE HYBRID DG FOR NAVIER-STOKES EQUATIONS

Chapter 2

Navier-Stokes equations with an exactlydivergence-free Hybrid DG method

The (incompressible) Navier-Stokes equations describe the dynamics of an incompressible flow. Theyconsist of one vector-valued equation for the transport of momentum and a mass conservation equationwhich comes as a constraint to the momentum equation. We consider two types of boundary condi-tions. Either the velocity is prescribed at a boundary (e.g. rigid walls or walls moving with a knownvelocity or inflow boundaries) which is a Dirichlet boundary condition or the momentum transport overa boundary (e.g. outflow bondaries) is known which is a Neumann boundary.According to Newton’s second law the change of momentum on a control volume balances with thesum of all forces acting on it and the in- and outflow of momentum. The outer forces sum up in f ,whereas the in- and outflow of momentum is decomposed into three parts:

• convective transport due to the velocity [ ρu⊗ u ]

• viscous forces [ −ηε(u) ]

• pressure forces [ p ]

Notice that we just consider newtonian fluids for the viscous force model, where ε(u) := 12(∇ u+∇T u).

Often the Navier-Stokes equations are rewritten in a form which replaces ε(u) simply by ∇ u. Theonly difference of both formulations is actually the natural boundary conditions. The natural boundarycondition for ε(u) = ∇ u seems to be the more appropriate in most cases.As we consider only incompressible flows, there holds ρ ≡ const. So we can devide by ρ and introducethe kinematic viscosity ν = η

ρ. The volume-specific pressure p and forces f will be reinterpreted as

mass-specific pressure and volume forces without change of notation from now on. Then we end upwith the incompressible Navier-Stokes equations

∂u∂t

+ div(−ν∇ u+ u⊗ u+ pI) = f in Ωdiv(u) = 0 in Ω

u = uD on ΓD(ν∇ u− pI) · n = 0 on Γout

(2.0.1)

Our objective is to solve this set of equations numerically. There are four properties of this equationswhich make the numerical treatment for Navier-Stokes equations challenging: This set of equations is

2.1. A SEMI-CONFORMING FINITE ELEMENT SPACE FOR NAVIER-STOKES 59

1. a saddle point problem

2. vector-valued

3. nonlinear

4. unsteady

The source of the saddle point structure becomes clear if we remember that the pressure can beseen as a lagrange multiplier for the incompressibility constraint. The fact that the Navier-Stokesequations are vector-valued is not a problem of it’s own but allows for distinguishing between normaland tangential component when formulating a numerical discretization. The nonlinearity of the Navier-Stokes equations is what makes them interesting: The steady problem may no longer have a uniquesolution and the dynamic of flow problems is often dominated by this nonlinearity. To be able toconcentrate on one problem after another we will neglect the time derivative in this chapter, i.e. wejust consider steady state problems. Instationarity will be the subject of chapter 3. Furthermore, we willstart without the nonlinear convective term and derive a numerical method as well as the analysis forthe so called Stokes problem in section 2.3. Then we will introduce a given linear convective velocity,which leads to the Oseen equations in section 2.4. Finally we will replace the convective velocity withthe fluid velocity and thereby recover the Navier-Stokes equations and discuss the numerical treatmentof it in section 2.5.Before we can introduce the method, we want to propose here, for the Stokes problem, we have tomake some preparations. Section 2.1 and 2.2 will motivate and introduce the use of H(div)-conformingFinite Element methods, which will be needed for the exactly divergence-free Hybrid DiscontinuousGalekin method later on.

2.1 A semi-conforming Finite Element space for Navier-Stokes

Solutions of the Navier-Stokes equations are [H1(Ω)]d-regular under reasonable assumptions on thedomain and the data. For ongoing discussions we will write [H1(Ω)]d in a non-standard and mesh-depend way:

[H1(Ω)]d = v ∈ [H1(Th)]d; [[v · n]]* = 0 on Fh︸︷︷︸~N

∩v ∈ [H1(Th)]d; [[v × n]]* = 0 on Fh︸︷︷︸~T

Here we divided the continuity on element interfaces of all components into tangential- and a normal-continuity.Our goal is to derive a method which fulfills the following properties:

1. The method should be locally and globally conservative, meaning a discrete (momentum) fluxshould be conserved

2. The method should be stable

3. The method should possess optimal convergence properties


To account for the first point we again consider Discontinuous Galerkin methods. Following [CKS05]a DG method which is locally conservative and energy-stable has to provide H(div,Ω)-conformingsolutions, that is, the discrete solution u has to be in H(div,Ω) = v ∈ [L2(Ω)]d, div(v) ∈ L2(Ω).In this work we will use Finite Element spaces which imply this H(div,Ω)-conformity for all ansatzfunctions. Those are taken to be piecewise polynomial and normal-continuous on element interfacesand thus are a subset of ~N .Tangential continuity is not included in the space so we have to account for it in another way. Wewill do this only in a (Hybrid) DG fashion, i.e. by adapting the formulation instead of the FiniteElement space. Standard DG formulations for the tangential component of the solution are of coursepossible here. But, as in the scalar case, all degrees of freedom of neighbouring elements would coupledirectly. As well as in the scalar case, in the vector-valued case this issue can easily be overcomedwith a hybridized formulation, which introduces utF , the approximation for the tangential trace of thesolution u.

H(div)-conforming DG DG

(conforming) CG H(div)-conforming HDG HDG

Figure 2.1.1: tangential and normal continuity for different methods

In Figure 2.1.1 a sketch of several possible discretizations of the velocity field is drawn. One verycommon approach is to consider the conforming continuous Galerkin approach. Here, all componentsare restricted to be continuous, which is owed to the [H1(Ω)]d-regularity of the problem. As a rapidchange in velocity and velocity gradients appears in many flows the standard Discontinuous Galerkinapproach is also often used. No conformity condition is imposed on the Finite Element space. Ofcourse a hybrid version is most often possible and introduces additional unknowns on the elementinterfaces for all components.The approach we are interested in is drawn in the middle. Tangential components are allowed to bediscontinuous across element interfaces and are considered for only by means of consistent penalization,i.e. through the formulation. A divergence-free DG version was already proposed in [CKS05]. It usesa local DG formulation. Our approach uses a hybrid version of the (symmetric) interior penaltyformulation which was used in [CKS07]. In our formulation only the tangential component of theadditional facet functions have to be introduced as normal continuity is enforced on the Finite Elementspace. The global linear systems that arise after discretization can be reduced to degrees of freedomsassociated with the facets, which are:

• normal flow degrees of freedom of the H(div)-conforming Finite Element space

• tangential facet degrees of freedom

2.2. H(DIV)-CONFORMING FINITE ELEMENTS 61

These are as much degrees of freedom on each facet as we would get for a hybridized version ofthe completely discontinuous approach. But we still have less overall degrees of freedoms due to theH(div)-conforming Finite Element space. The construction of H(div)-conforming Finite Elements,which we need for our formulation, is the subject of the next section.

2.2 H(div)-conforming Finite Elements

This section is about H(div)-conforming discretizations with Finite Elements. We want to present theconcepts and main ideas which are necessary for the comprehension and construction of higher orderH(div)-conforming Finite Elements.H(div)-conforming Finite Elements are supposed to provide that every discrete function uh of thediscretization lies in H(div,Ω), that is, that

uh ∈ H(div,Ω) = uh ∈ [L2(Ω)]d : div(uh) ∈ L2(Ω) (2.2.1)

We work with piecewise polynomials on each element. Thus locally the functions are automaticallyin H(div, T ). But similar to H1-conforming Finite Elements we have to enforce another compatibilitycondition on the global functions to satisfy uh ∈ H(div,Ω). For H1-conformity the requirement is thecontinuity over element interfaces, for H(curl)-conformity1 it would be the continuity of the tangentialcomponent over element interfaces and in our case, we need continuity of the normal-component.

Lemma 2.2.1. Let Th be a triangulation of a domain Ω. A function v is in H(div,Ω) if and only ifv ∈ H(div, T ) ∀ T ∈ Th and [[v]]n = 0 on all interior facets E ∈ F inth .

We will speak of element contributions to the divergence, which is simply div(·) on each elementT (which may represent sources and sinks) and facet contributions, which come from the normalcomponents on the element boundaries (think of the sum of in- and outflows).H(div)-conforming discretizations were first used in the context of mixed formulations of Poisson’sequation, as we have already seen in the last chapter. Typical examples for such spaces are the(lowest order) Raviart-Thomas and the (lowest order) Brezzi Douglas Marini Elements which we wantto discuss in the next section. Afterwards the Piola transformation will be discussed, which allowsfor reducing the problem of construction Finite Element shape functions on arbitrary elements tothe same problem on a reference element. Further, we present the exact sequence property of H1,H(curl), H(div) and L2 in section 2.2.3 which is the foundation for the chosen realization of higherorder H(div)-conforming methods of section 2.2.4. We conclude with interpolation operators on anH(div)-conforming space.All results on H(div)-conformity will be given without proof as we want to focus on the conceptsin this section. Additionally we will focus on the two dimensional case involving triangles only. Forextensions to three dimensions and proofs we refer to [Zagl06], [SZ05], [DB05] and [DGS08].

2.2.1 Lowest Order H(div)-conforming Finite Elements

Let’s start looking at lowest order realizations of those Finite Elements.1H(curl,Ω) := uh ∈ [L2(Ω)]d : curl(uh) ∈ [L2(Ω)]d


Lowest order Raviart Thomas Finite ElementsThe most famousH(div)-conforming Finite Element is the lowest order Raviart Thomas Finite Element,abbreviated by RT0. An RT0 Finite Element has m shape functions, where m is the number of facets(faces (in 3D) or edges (in 2D)). Let’s focus on two dimensional, triangular meshes here. Those trialfunctions have non-zero normal component one on one facet and zero on all others. Accordingly everyfacet possesses one degree of freedom, which is the total flow over a facet:

N0E(ψ) =

∫Eψ · n ds (2.2.2)

With the help of the barycentric coordinates λi, the trial functions can be easily constructed. The trialfunction associated with the edge E connecting the nodes (w.l.o.g.) 1 and 2 is

ψ0E

(λ1, λ2) = λ2 curl(λ1)− λ1 curl(λ2) with curl =(

∂x2

−∂x1

)

Those trial functions’ element contribution to the divergence are constant by construction as is thenormal flow on each edge. We note further that the local space of the RT0 Finite Element is

RT0(T ) = v(x) = a+ x · b, a ∈ Rd, b ∈ R (2.2.3)

and as div(x) = d = const the relation

div(RT0(Th)) = Q0h(Th) , with RT0(Th) :=

⊕E∈Fh

spanψ0E ⊂ H(div,Ω) (2.2.4)

holds.

Lowest order Brezzi Douglas Marini Finite ElementsWe can enrich the RT0 element by adding one additional shape function on each facet. We thereforeintroduce the additional DOF which we mark with a star:

N∗E(ψ) =∫Eψ · n v ds for v ∈ P1(E) ∩ P0(E)⊥

NE N∗E

Figure 2.2.1: Degrees of freedom for the lowest order BDM element

The associated shape functions are

ψ∗E

(λ1, λ2) = curl(λ1λ2)


There holds div(ψ∗E

) = 0 and NE(ψ∗E

) = 0 ∀E ∈ ∂T . The element consisting of both DOF on eachfacet is called Brezzi Douglas Marini element of first order, or short BDM1. The number of DOF ona simplex exactly matches the dimension of the polynomial space [P1(T ))]d and indeed the proposedshape functions are linearly independent and span the local space [P1(T ))]d.If we identify the degrees of freedom associated with the faces of the mesh we obtain the global BDM1Finite Element Space

BDM1(Th) :=⊕E∈Fh

spanψ0E, ψ∗

E ⊂ H(div,Ω) (2.2.5)

Both elements, RT0 as well as BDM1 fulfill

div(RT0(Th)) = div(BDM1(Th)) = Q0h(Th) (2.2.6)

but the approximation of the function itself is better for BDM1 as

[P0(T )]d ⊂ RT0(T ) ⊂ BDM1(T ) = [P1(T )]d (2.2.7)

holds.For both approaches there exist higher order versions. The local spaces RTk(T ) and BDMk(T ) are

RTk(T ) = [Pk(T )]d ⊕ x · Pk,∗(T ) (2.2.8)BDMk(T ) = [Pk(T )]d (2.2.9)

where Pk,∗(T ) are all polynomials on T of exact degree k. So the index k of RT and BDM indicateswhich polynomial space is completely contained within the respective local spaces. The relationsbetween the higher order spaces is similar to the lower order ones:

div(RTk(Th)) = div(BDMk+1(Th)) = Qkh(Th) (2.2.10)

and[Pk(T )]d = BDMk(T ) ⊂ RTk(T ) ⊂ BDMk+1(T ) = [Pk+1(T )]d (2.2.11)

We see that the approximation of the divergence is, w.r.t. the total degrees of freedom, faster achievedwith Raviart Thomas elements. But if the approximation of the function itself is important the methodof choice is the BDM element. In our case we use the BDM element as we are especially interestedin the velocity-field itself. How to choose the functionals N and how to construct appropriate trialfunctions for the high order BDM-case is explained in section 2.2.4. Before we turn over to that, wewill show how H(div)-conforming Finite Elements on the physical domain can be achieved by meansof transforming Finite Elements from the reference domain to the physical domain in a way whichpreserves the normal component at the boundary and thereby the H(div)-conformity.

2.2.2 H(div)-conforming Transformation (Piola Transformation)

To obtain globally H(div)-conforming Finite Elements, i.e. shape functions whose normal componentsat element interfaces coincide also on the physical domain, we make use of the H(div)-conformingtransformation also known as Piola transformation. It further ensures that the introduced degrees offreedom are meaningful also on the physical domain. The following statements on the Piola transfor-mation are given without proof. Those can be found in [Zagl06] or [Monk03, section 3.9]


Lemma 2.2.2. (H(div)-conforming transformation) Let ΦT : T → T be a diffeomorphism andu ∈ H(div, T ). Then the Piola transformation

u := J−1T F

Tu Φ−1

T

with FT

= ∇ ΦT , JT = det(FT

), implies u ∈ H(div, T ). The divergence of u is transformed evensimpler:

divx(u) := J−1T divx(u Φ−1

T )

This transformation results in some further desirable properties:

1. The normal components transform like

u(x) · n(x) = u(x) · n(x) · 1JT‖F T

T n‖

The transformation factor cancels out (except for the sign) with the transformation factor ofsurface integrals ∫

Eds =

∫E|JT |‖F−TT n‖ ds

which leads to (except for the sign) invariant surface integrals∫Eu · nE · q ds = sign(JT )

∫Eu · nE · q ds (2.2.12)

2. Integrals of the form∫T div(u) · q dx are invariant (except for the sign) w.r.t. the Piola trans-

formation: ∫Tdivx(u) · q dx = sign(JT )

∫Tdivx(u) · q dx (2.2.13)

3. In the case of affine linear transformation, the Raviart Thomas space RTk(T ) is invariant withrespect to the Piola transformation2.

With the help of the Piola transformation shape functions have to be defined only on the referencedomain. Additionally some integral terms are geometry-independent what can be exploited for a fastermatrix assembly. Using the Piola transformation on curved elements may lead to non-polynomial localspaces on the physical domain as the local space on the reference domain is set to be polynomial and thetransformation doesn’t have to be. In this case the Piola transformation still ensures H(div)-conformityas the transformation of the normal component coincides for two (curved) adjacent elements.

2.2.3 The Exact sequence and the DeRham Complex

The exact sequence of the Sobolev spaces H1, H(div), H(curl) and L2 characterizes the relationshipbetween spaces and the according differential operators and provides a foundation for the constructionof the higher order H(div)-conforming Finite Elements we want to use.

2The same holds trivially for the Brezzi Douglas Marini Element space BDMk


2.2.3.1 Exact sequence of H1, H(div), H(curl) and L2

Let’s first state the exact sequence property on the continuous level in two and three dimensions. Itis called the de Rham Complex and in three dimensions it is written down as

R id−→ H1(Ω) ∇−→ H(curl,Ω) curl−→ H(div,Ω) div−→ L2(Ω) 0−→ 0 (2.2.14)

and should be read in the following way: the kernel of the operator standing on the right side of thespace is exactly the image of the operator standing on the left side of the same space. For examplethe kernel of the curl operator in H(curl,Ω) coincides with the gradients of H1(Ω).The exact sequence implies the well known identities curl(∇(·)) = 0 and div(curl(·)) = 0. It actuallytells us more: Assume div(u) = 0 for a function u ∈ H(div,Ω), then there exists a function φ ∈H(curl,Ω), s.t. u = curl φ or assume curl u = 0 for a function u ∈ H(curl,Ω), then there exists afunction φ ∈ H1(Ω), s.t. u = ∇φ.

In two dimensions, the sequence is reduced by one stage, which is either H(curl) and H(div). Thuswe have two different versions of the de Rham Complex :

R id−→ H1(Ω) ∇−→ H(curl,Ω) curl−→ L2(Ω) 0−→ 0 (2.2.15a)

R id−→ H1(Ω) curl−→ H(div,Ω) div−→ L2(Ω) 0−→ 0 (2.2.15b)

The latter is the more important one for us as it involves the space H(div,Ω).

2.2.3.2 The exact sequence for Finite Element spaces

Conforming Finite Elements which also fulfill the exact sequence property are desirable as the separationinto the kernel of the according differential operator and the remainder of the space can be exploited.The following diagram includes the finite dimensional analog to the de Rham Complex of the continuousspaces and presumes conforming discretizations (we just consider the 2D case (2.2.15b)).

R id−→ H1(Ω) curl−→ H(div,Ω) div−→ L2(Ω) 0−→ 0⋃ ⋃ ⋃R id−→ Wh

curl−→ Σhdiv−→ Qh

0−→ 0

(2.2.16)

2.2.3.3 Exact sequence for the lowest order Finite Element spaces

The group of Finite Elements consisting of first order nodal Finite Elements, lowest order RaviartThomas Finite Elements and piecewise constant Finite Elements possesses the exact sequence prop-erty3, which is drawn in Figure 2.2.2. The two numbers below the drawn Finite Element indicate thedimension of the image of the left hand side operator which coincides with the dimension of the kernelof the right hand side operator (left number) and the dimension of the rest of the Finite Element space(right number). The sum of both is the number of degrees of freedom on each element.

3in 3D, the lowest order Nedelec Finite Elements of first kind would be the appropriate representative for H(curl)Finite Elements


R id−→ curl−→ div−→ 0−→ 0

1 + 2 2 + 1 1 + 0

Figure 2.2.2: exact sequence involving the lowest order Raviart Thomas element

Another exact sequence involves the BDM1 element and is the following.

R id−→ curl−→ div−→ 0−→ 0

1 + 5 5 + 1 1 + 0

Figure 2.2.3: exact sequence involving the lowest order BDM element

We see that the additional edge DOF of the H1-conforming element lead to three additional DOFfor the H(div)-conforming Finite Element space. But as we already saw, possible additional shapefunctions are ψ∗

E(λ1, λ2) = curl(λ1λ2) which are curl-fields (of H1-functions)4 and thereby divergence-

free. The last step of the sequence is thus uneffected. This exact sequence reduces the polynomialorder by one for each differential operator in contrast to the one involving the Raviart Thomas element.

2.2.4 Construction of Higher Order H(div)-conforming Finite Elements

After having discussed the lowest order cases, we want to turn to the construction ofH(div)-conformingFinite Elements with arbitrary polynomial degree > 1. Several realizations are possible. Our approachof choice is the one proposed and extensively discussed in [SZ05, Zagl06], where the construction ofvariable order H1-, H(curl)-, H(div)- and L2-conforming Finite Element spaces on hybrid meshes withtetrahedrals, hexahedrals, prisms, triangles and quadrilaterals is presented. The construction of thehigher order shape functions mimics the exact sequence property of the spaces H1, H(curl), H(div)and L2 and thereby facilitates further approvements.We proceed as follows: For the H(div)-conforming case only, we will show which face- and cell degreesof freedom should be used to enable an H(div)-conforming Finite Element space which inherits theexact sequence property from the continuous level. For completeness we will also give the set of basisfunctions which were derived in [SZ05, Zagl06]. These use the separation of the element-local FiniteElement space into low order RT0, face-associated, divergence-free element-associated and polynomialspace completing basis functions. This classification allows for reducing the degrees of freedom withoutloss of accuracy, which will be discussed in the last part of this section. As before, we will limit ourdiscussion on the concepts and the two-dimensional case with a p-uniform triangular meshes. Forthe more involved cases and further details we recommend to consult [Zagl06], where all details areincluded and comprehensibly explained.

4λ1λ2 is actually the standard second order edge bubble function


2.2.4.1 Transfering the exact sequence property to the discrete level

We now want to obtain a sequence of higher order Finite Element spaces (Wh,Σh, Qh) which fulfill(2.2.16). All Finite Element spaces should be approximated with piecewise complete polynomials,which is the generalization of the lowest order case of the exact sequence which involves the BDM1element. We devide the degrees of freedom into three classes:

1. lowest order degrees of freedom (vertex-, edge- and cell-based5)

2. higher order edge6 degrees of freedom

3. higher order element degrees of freedom

The lowest order degrees of freedom are exactly those we already used in Figure 2.2.2, i.e. vertexshape functions, RT0 and piecewise constants. For each polynomial degree an additional degree offreedom is associated with an edge for Wh and Σh. These are called higher order edge degrees offreedom and ensure continuity and normal-continuity, respectively. All other degrees of freedom arehigher order element degrees of freedom and have to provide shape functions which are zero at theelement boundary for Wh or which have zero normal component for Σh. For Qh no such classificationhas to be made as no continuity restriction is necessary.In the next section we will give an appropriate choice of moments forH(div)-conforming Finite Elementsonly, before we discuss the separation of the Finite Element space into four classes.

Moments / Degrees of freedomWe already discussed the lowest order degrees of freedom related to the RT0 and BDM1 Finite Elementsin section 2.2.1. In the high order case we use the following degrees of freedom for a Finite Elementof uniform order k:

• lowest order edge-based DOF:NE(φ) =

∫Eφ · n ds

• higher order edge-based DOF:

N∗E(φ) =∫Eφ · n v ds for v ∈ Pk(E) ∩ P0(E)⊥

• higher order nonzero divergence cell-based DOF:

NAC (φ) =

∫Tφ · vi dx for via basis of curl([Pk+1]dt,0)

where [Pk]dt,0 = [Pk]d ∩H0(curl, T )

• higher order divergence-free cell-based DOF:

NBC (φ) =

∫Tdiv(φ) div(vi) dx for vi s.t div(vi)is a basis of div([Pk]dn,0)

where [Pk]dn,0 = [Pk]d ∩H0(div, T )5in 3D vertex-, face-, edge- and cell-based6in 3D face and edge


The last separation of cell-based DOF is owed to the exact sequence as it naturally divides them intoshape functions which have nonzero divergence (↔ NB

C ) and those who are divergence-free (↔ NAC ).

Additionally the higher order edge-based shape functions can be taken as the curl of H1 edge-basedshape functions and are thereby also divergence-free.We consider the following shape functions related to the degrees of freedom we just defined:

1. edge-based shape functions for edge E connecting vertices (w.l.o.g.) 1 and 2:

(a) lowest order:

ψ0E

= curl(λ1)λ2 − curl(λ2)λ1

(b) higher order:

ψiE

= curl(LSi+1(λ1 − λ2, λ1 + λ2)) [= curl(ϕi+1E )], for 1 ≤ i ≤ k

2. cell-based shape functions for element T sharing vertices 1, 2 and 3:

(a) divergence-free:

ψA;i,jC

= curl(ui)vj − curl(vj)ui = curl(ui vj) [= curl(ϕi,jC )] for 0 ≤ i+ j ≤ k − 2

(b) nonzero divergence:

ψB1;i,jC

= curl(ui)vj + curl(vj)ui for 0 ≤ i+ j ≤ k − 2ψB2;jC

= (curl(λ1)λ2 − curl(λ2)λ1)vj for 0 ≤ j ≤ k − 2 (2.2.17)

where ui := LSi+2(λ1 − λ2, λ1 + λ2) and vj := λ3lj(2λ3 − 1). Here we made use of the ScaledIntegrated Legendre Polynomials of degree i+2 LSi+2 and the Legendre polynomials lj of degreej . Both are defined in the appendix A.3.

The functions ϕi+1E and ϕi,jC are the edge- and cell-based shape functions used in [SZ05, Zagl06] for the

H1-conforming Finite Element space. So we see curl(Wh) is completely included in Σh and furthermorethe subspace curl(Wh \ φV ), where φV , the set of all lowest order vertex-based shape functions,is explicitly separable, i.e. we can, except for the lowest order RT0 functions, divide the space intodivergence and divergence-free functions.Before we will make use of this, let’s summarize the properties of the separation of the Finite Elementspace in the next Table (2.1). Note that now each symbol (e.g. an arrow) corresponds not to exactlyone degree of freedom but to a class of degrees of freedom.


RT0 DOF higher orderedge DOF

higher order div.-free DOF bubble

higher orderDOF with nonz.div

Σkh =

⊕E∈Fh

spanψ0E ⊕

⊕E∈Fh

spanψiE ⊕

⊕T∈Th

spanψA;i,jC ⊕

⊕T∈Th

spanψB1;i,jC

, ψB2;jC

div(Σkh) =

⊕T∈ThP0(T ) ⊕ 0 ⊕ 0 ⊕

⊕T∈Th

[Pk−1 ∩ P0⊥](T )

#DOF = 3 + 3k + 12k(k − 1) + 1

2k(k + 1)− 1

Table 2.1: Separation of the H(div)-conforming Finite Element space

Let’s also try a more physical interpretation of this space separation. Every flow field can be interpretedas a superposition of a divergence- and vortex-free mean flow, (smaller and larger) eddies, sources andsinks. The eddies can again be subdivided into three kind of eddies:

1. eddies aroung one or more vertices

2. eddies on an edge, i.e. eddies which have zero normal velocity on all edges except for one

3. eddies which lie within one element

To resolve the mean flow and the first kind of eddies we need the lowest order RT0 shape functions.Now the eddies we need to represent in order to increase the resolution can either be located at oneedge or within an element.

2.2.4.2 Excluding higher order divergence functions - Reducing the basis

With the introduced separation, we can neglect all higher order nonzero divergence shape functions ifwe consider an incompressible flow explicitely. For u ∈ Σk

h there holds

• the divergence of the higher order divergence shape functions are linear independent (due to thechoice of the functionals N)

• the divergence of the space of all higher order divergence shape functions is 12p(p + 1) − 1

dimensional on each element

• the divergence of the higher order divergence shape functions is L2-orthogonal to (element-piecewise) constant functions (due to

∫∂T ψ

B

C· n ds =

∫T div(ψBC) = 0)


• we can uniquely decompose u into uRT0E , uhoE , uAC and uBC due to the separation of the space.

• div(u) = 0 now requires div(uRT0E ) = 0, where div(uRT0

E ) ∈ Q0h(Th) and div(uB) = 0, where

div(uB) ∈ Pk−1(T ) ∩ P0(T )⊥ on each element T

• This gives 12p(p+ 1)− 1 equations for uB which is a linear combination of 1

2p(p+ 1)− 1 shapefunctions.

Thus div(u) = 0 always implies uB = 0, s.t. those DOF can be neglected.

This reduces the number of unknowns drastically. Additionally for the incompressible Navier-Stokesequations the pressure unknowns, which are lagrange multipliers for the incompressibility constraint,can be reduced to the lowest order case, i.e. piecewise constant functions, because we only have tocontrol the RT0 degrees of freedom. The overall unknowns can thereby approximately halfed withoutloss of accuracy for higher order approximation (p ≥ 2, see section 2.5.2). The neglected velocityDOF would be zero the one way or the other. Only the pressure approximation has to be recoveredwith a post-processing step. This is not discussed in the proceeding. So instead of Σk

h we can workwith the space

Σk,∗h :=

⊕E∈Fh

spanψ0E ⊕

⊕E∈Fh

spanψiE ⊕⊕T∈Th

spanψAC (2.2.18)

2.2.5 H(div)-conforming Interpolation

At several places during the analysis later on we will need an H(div)-conforming interpolation ΠΣh,k :

[H1]d → Σkh , which possesses the following properties:

• H(div)-conformity: ΠΣk (u) ∈ Σk

h

• L2-optimal approximation of the divergence: (div(ΠΣk (u)), qh) = (div(u), qh) ∀ qh ∈ Qk−1

h

• boundedness of the interpolation in the H1 norm: ‖ΠΣk (u)‖H1 ≤ C‖(u)‖H1

• the interpolator is a projector, s.t. ΠΣk (v) = v ∀ v ∈ Σk

h

Those properties can be achieved with the interpolation operator proposed in [BF91, III.3.3]. Also theprojection-based interpolation proposed in [DB05] fulfills all the above requirements and is furthermorep-robust. Proofs and further discussion can be found in [DB05] and [DGS08].

2.3. DIVERGENCE-FREE HYBRID DG FOR STOKES EQUATIONS 71

2.3 Exactly divergence-free Hybrid DG method for the Stokesproblem

In this section we want to present the variational formulation as well the a priori error analysis for thesteady Stokes problem, which reads as

div(−ν∇ u+ pI) = f in Ωdiv(u) = 0 in Ω


(2.3.1)

On Γout the natural boundary conditions are used. Those boundaries prescribe zero momentum fluxon the outflow region. It naturally comes up like Neumann boundaries (∂u

∂n[= 0] ) in the scalar case.


As we already pointed out in section 2.1 our approach to a good (H)DG method is to decompose theH1-consistency into normal and tangential consistency. Both components will be treated differently.The normal consistency (i.e. normal continuity) will be enforced by the choice of the discrete space.The tangential consistency will be included with a HDG formulation. The procedure to reduce trial andtest functions will not be adressed explicitely as the results stay the same and only make a differencefor the implementation.

viscous part

The procedure of deriving our Hybrid DG formulation for the viscous part is similar to the one usedfor scalar diffusion. First we integrate by parts on each element.∫

Ωdiv(−ν∇u)v dx =

∑T∈Th

∫Tdiv(−ν∇u)v dx =

∑T∈Th

∫Tν∇u : ∇ v dx−

∫∂Tν∂u

∂n· v ds (2.3.2)

Adding a consistent term ∑T∈Th

∫∂Tν∂u

∂n· v ds =

∫ΓNν∂u

∂n· v ds (2.3.3)

with the function v = vn + vtF , which is single-valued, results in∑T∈Th


∫∂Tν∂u

∂n· (v − v)︸︷︷︸

=[[vt]]

ds−∫

ΓNν∂u

∂n· v ds (2.3.4)

where we used v − v = vt − vtF =: [[vt]]. The boundary integral allows for naturally plugging in theoutflow boundary condition. On the other parts the boundary term will vanish as v is zero there.Similar to the procedure in the scalar case we can now use that [[ut]] is zero for the true solution toadd consistency terms for symmetry and stability. Finally we end up with

Bh(u, v) :=∑T∈Th

∫Tν∇u : ∇v dx−

∫∂Tν∂u

∂n·[[vt]] ds−

∫∂Tν∂v

∂n·[[ut]] ds+

∫∂Tν τh [[ut]]·[[vt]] ds (2.3.5)

This actually is the analog to the scalar case. Instead of [[ut]] and [[vt]] you could also write [[u]] and[[v]] as [[un]] and [[vn]] is always zero due to H(div)-conformity. However, we stick with [[ut]] and [[vt]]to highlight the decomposition of the continuity.


pressure part

For the pressure we again integrate by parts on each element and replace the boundary value for pwith a numerical flux p.∫

Ω∇p · v dx =

∑T∈Th

∫T∇p · v dx =

∑T∈Th−∫Tp div(v) dx+

∫∂Tp vn ds (2.3.6)

We choose p =p on inner facetp on boundary facets , s.t.

∑T∈Th

∫∂Tp vn ds =

∫Γout

p vn ds and results in

∑T∈Th−∫Tp div(v) dx+

∫Γout

p vn ds (2.3.7)

The boundary integrals are part of the natural boundary condition and so won’t be considered for inthe bilinearform:

Dh(v, p) :=∑T∈Th−∫Tp div(v) dx (2.3.8)

The discrete problem

Now we have described the formulation, but didn’t specify completely which spaces we use. First wehave to introduce a new facet Finite Element space for the tangential facet functions:

F kh :=

utF : utF ∈ [Pk(E)]d × n ∀ E ∈ Fh

(2.3.9)

Then our velocity approximation u is again a composition of an element function u ∈ Σkh and uF ∈ F k

h

which we again indicate by the bold writing.The compound space is called Skh := Σk

h × F kh and Skh,D = Σk

h,D × F kh,D (or D = 0) denotes that

(homogeneous) dirichlet boundary conditions are incorporated into the space. Normal components ofboundary values have to be imposed on Σk

h and tangential ones on F kh . That means that the normal

component is imposed in a strong sense, whereas the tangential component is only enforced in a weak(Hybrid) DG sense.Equally to the scalar case, the pressure space is the space Qk

h which consists of discontinuous, piecewisepolynomials of order k. It is chosen to be one order less than the velocity approximation.Then the discrete problem reads:Find u ∈ Sk+1

h,D and p ∈ Qkh, s.t.Bh(u,v) +Dh(v, p) = 〈f, v〉 ∀ v ∈ Sk+1

h,0Dh(u, q) = 0 ∀ q ∈ Qk

h

(2.3.10)

We also introduce the larger Bilinearform to shorten the notation at some parts of the analysis:

Kh ((u, p), (v, q)) := Bh(u, v) +Dh(v, p) +Dh(u, q) (2.3.11)

Thus (2.3.10) can also be written as (with U = (u, p) the velocity-pressure pair)

Kh (U, V ) = 〈f, v〉 ∀ V ∈ Sk+1h ×Qk

h =: Zkh (2.3.12)


We additionally introduce an appropriate norm:

‖U‖Kh = ‖(u, p)‖Kh :=√ν|||u|||1 + ‖p‖L2 (2.3.13)

where ||| · |||1 is the hybrid version of the modified broken H1-norm, similar to the one we saw duringthe analysis of the scalar poisson problem:

|||u|||21 :=∑T∈Th

‖∇ u‖2

T + 1h

[[[ut]]]2∂T + h‖∂u∂n‖2∂T

(2.3.14)

Remark 2.3.1 (Weak incompressibility ⇒ strong incompressibility):Due to the choice of the velocity-pressure pair, there holds div(Σk+1

h ) ⊆ Qkh. Then Dh(u, q) = 0 ∀ q

results in div(u) = 0 in a strong sense on each element. Additionally taking normal continuity, i.e.H(div)-conformity, into account gives div(u) = 0 globally, thus the approximate solution is actuallyexactly divergence-free. In other words the formulation imposes a strong incompressibility conditionon u.

Remark 2.3.2 (Reducing the facet unknowns):As was pointed out in Remark 1.2.4 for the scalar case, in the pure diffusive limit the jump operatorcould be replaced by a projection on the polynomial space of a degree lower and so the facet functionscould be chosen of an order less. As soon as we also consider convection this is not possible withoutreducing the convergence rate.

Interpolation on the hybrid space ShOwed to the fact that tangential facet functions where added to Σk

h, resulting in Skh, we have tointroduce an appropriate interpolator / projector for the complete Finite Element space Sh:

ΠSh,k(u) :=

(ΠΣh,k(u), Πf,t

h,k(u))

with Πf,th,k(u) ∈ F t,k

h and(Πf,th,k(u), vtF

)∂T

=(tr|Fh(u), vtF

)∂T

∀ vtF ∈ F thk(T )

(2.3.15)

Notice that both projections w = ΠΣh,k(u), wtF = Πf,t

h,k(u) applied on adjacent elements result in thesame facet function w = wn+wtF , so the global projector can be constructed in an element by elementfashion. Furthermore the projection ΠS

h,k(u) is bounded in the ||| · |||1 norm as ΠΣh,k(u) and Πf,t

h,k(u)are bounded in weaker norms. Of course, as ΠS

h,k(u) is actually a projection onto Skh there holdsΠSh,k(uh) = uh ∀ uh ∈ Skh.

Let’s state the following Lemma, which will be made use of shortly after:

Lemma 2.3.1. The interpolation operator ΠSh,k : [H1(Ω)]d → Skh for k ≥ 1 fulfills the following

relations with c independent of the mesh size:

Dh(u, q) = Dh(ΠSh,ku, q) ∀ q ∈ Qk−1

h (2.3.16)|||ΠS

h,ku|||1 ≤ c‖u‖H1 (2.3.17)

Proof. See the appendix, chapter B.



To show well-posedness and later on optimal convergence of the method we will make use of Brezzi’sTheorem. Essential conditions for this theorem are the coercivity of Bh, the boundedness of bothbilinearforms Bh and Dh and an inf-sup condition on Dh.For the continuous problem with the canonical bilinearforms, i.e. replacing Bh with

∫Ω∇ u ∇ v dx,

coercivity, boundedness and inf-sup condition hold, and due to the subsequent Theorem the problemis well-posed and thus has a unique solution u ∈ [H1(Ω)]d and p ∈ L2(Ω). The task now is to showthat the same properties translate to our discretization. Therefore many of the scalar results can beovertaken for Stokes.

The boundary is splitted into two parts ΓD and Γout where either Dirichlet boundary conditions orNeumann boundary conditions are prescribed. By Neumann boundary conditions we mean bound-aries where

(ν∇ u− pI

)· n is given. We call this part of the boundary Γout as we just consider

homogeneous Neumann boundary conditions which are normally the outflow conditions of a domain.Further we assume that wherever Dirichlet boundary conditions are prescribed, they are prescribed forall components (just for ease of notation) and that ΓD is a nonempty part of the boundary.

We will use the following discrete and continuous spaces during the analysis:

BDMk = Σkh ⊂ Σ = H(div,Ω)

utF : utF ∈ [Pk(E)]d × n ∀ E ∈ Fh

= F t,kh ⊂ F t =

utF : utF ∈ [L2(Fh)]d × n

Σkh × F

t,kh = Skh ⊂ S = Σ× F t

p : p ∈ Pk(T ) ∀ T ∈ Th

= Qkh ⊂ Q = L2(Ω)

Skh ×Qk−1h = Zk

h ⊂ Z = [H1(Ω)]d ∩ [H2(Th)]d × L2(Ω)

(2.3.18)

We introduce further the spaces Σ∗,kh and S∗,kh which in contrast to Σk,∗h (which was the Finite Element

space without the higher order nonzero divergence functions) indicate the divergence-free subspacesof Σk

h and Skh:

S∗,kh := u ∈ Skh : div(u) = 0 Rem.2.3.1⇐⇒ u ∈ Skh : Dh(u, q) = 0 ∀ q ∈ Qk−1h (2.3.19)

As we need to mimic H1 regularity the BDM approach fits better than Raviart Thomas element aswith the full polynomial space the approximation of ∇ u and u is more natural.The identity operator idh : [H1(Ω)]d ∩ [H2(Th)]d → S, which uses the well-defined tangential tracefor the tangential facet functions, will again be used at several occasions implicitely. It provides anextension of the bilinearforms, s.t. Bh : S × S → R and Dh : Q× S → R are well defined.For ease of presentation we only consider constant viscosity ν.

Brezzi’s TheoremWe restate a rather extensive version of Brezzi’s Theorem:


Theorem 2.3.2 (Brezzi’s Theorem). Assume that we have two Hilbert spaces V and Q and twocontinuous bilinearforms a : V × V → R and b : Q× V → R

a(u, v) ≤ βa‖u‖V ‖v‖V ∀ u, v ∈ Vb(v, q) ≤ βb‖v‖V ‖q‖Q ∀ u ∈ V, ∀ q ∈ Q

and associated linearforms f : V → R and g : Q → R, which are bounded in the operator norms‖ · ‖V ′ and ‖ · ‖Q′ . Assume further that coercivity of a(·, ·) on the kernel of b, i.e.

a(u, u) ≥ αa‖u‖2V ∀ u ∈ V ∗ := v ∈ V : b(v, q) = 0 ∀ q ∈ Q

and the inf-sup condition holds

inf06=q∈Q

supu∈V

b(u, q)‖u‖V ‖q‖Q

≥ αb > 0 ⇔ supu∈V

b(u, q)‖u‖V

≥ αb‖q‖Q ∀ q ∈ Q

Then, the compound bilinearform

k ((u, p), (v, q)) := a(u, v) + b(u, q) + b(v, p)

is inf-sup-stable, i.e. with the definition of the compound functions, spaces and norms W := V ×Q,U := (u, p), V = (v, q) and ‖U‖W := ‖u‖V + ‖p‖Q there holds:

supV ∈W

k(U, V )‖V ‖W

≥ αk‖U‖W ∀ U ∈ W

So the mixed problema(u, v) + b(v, p) = f(v) ∀ v ∈ Vb(u, q) = g(q) ∀ q ∈ Q ⇔ k(U, V ) = f(v) + g(q) ∀ V ∈ W

is uniquely solvable and the solution fulfills the stability estimate

‖U‖W = ‖u‖V + ‖p‖Q ≤ c (‖f‖V ′ + ‖g‖Q′)

with the constant c depending only on αa, αb, βa, βb.

Our task now is to show stability, consistency and boundedness for the discrete problem to get uniquelysolvability and the stability estimates. One essential ingredient for this is the treatment of the incom-pressibility constraint. The inf-sup condition on the discrete level is often the most important criteriafor stability of numerical methods for incompressible flows. If αb of Brezzi’s Theorem or in our caseαDh is required to be independent of the mesh size the inf-sup criteria is most often called the LBB(Ladyshenskaja-Babuška-Brezzi) condition. We will again show consistency, stability, boundedness andapproximation one after another before we conclude with a priori error estimates in a discrete brokenH1 norm.

2.3.2.1 Consistency

One important property of the discretization is the consistency, i.e. that the discrete equation is fulfilledif we replace the discrete solution with the exact one. The incompressibility condition is trivially fulfilled


as the exact solution is exactly divergence-free.Let (u, p) ∈ Z be the solution of (2.3.1). Then for all q ∈ L2(Ω)

Dh(u, q) =∑T∈Th

∫T− div(u)︸︷︷︸

=0

q dx = 0 (2.3.20)

And for Bh +Dh we get for all v ∈ S0:

Bh(u,v) +Dh(v, p)

=∑T∈Th


∫∂Tν∂u

∂n[[vt]] ds−

∫∂Tν∂v

∂n[[ut]]︸︷︷︸=0

ds

+∫∂Tντh [[ut]]︸︷︷︸

=0

[[vt]] ds−∫Tdiv(v) p dx

=∑T∈Th

∫Tdiv(−ν∇u) v dx+

∫∂Tν∂u

∂n

(v − [[vt]]

)︸︷︷︸

=vn+vtF

ds

+∫T∇p v dx

−

∑E∈Fint

h

∫E

[[p]]*︸︷︷︸=0

vn ds−∑

E∈Fexth

∫Ep vn ds

(∗)=∑T∈Th

∫Tdiv(−ν∇u+ pI)︸︷︷︸

= f

v dx+∑

E∈Finth

∫E

[[ν ∂u∂n

]]︸︷︷︸=0

(vn + vtF

)ds

+∑

E∈Fexth

∫E

(ν∇ u− p I

)· n(vn + vtF

)ds

=∫

Ωf v dx+

∫ΓD

(ν∇ u− p I

)· n(vn + vtF

)︸︷︷︸

=0

ds

+∫

Γout

(ν∇ u− p I

)· n︸︷︷︸

=0

(vn + vtF

)ds

=∫

Ωf v dx

(2.3.21)

With those preparations we can conclude the following:

Proposition 2.3.3 (Galerkin Orthogonality). Let Uh = (uh, ph) ∈ Zh,D be the solution of (2.3.10)and U = (u, p) ∈ Z be the solution of (2.3.1). Then there holds

Bh(uh − u,vh) +Dh(uh − u, qh) +Dh(vh, ph − p) = 0 ∀ (vh, qh) ∈ Sh,0 ×Qh (2.3.22a)⇐⇒ Kh(Uh − U, Vh) = 0 ∀ Vh ∈ Zh,0 (2.3.22b)

Proof. There holds Dh(u, qh) = Dh(uh, qh) = 0 ∀qh ∈ Qh (see eq. (2.3.20)). And with (2.3.21) weobtain

Bh(uh,vh) +Dh(vh, ph)− Bh(u,vh)−Dh(v, ph) =∫

Ωf v dx−

∫Ωf v dx = 0 ∀ vh ∈ Sh,0

(2.3.23)which completes the proof.


2.3.2.2 Stability

The stability of the mixed problem consists of two parts, the (kernel-) coercivity of Bh and the LBBcondition of Dh. For Bh we can reuse the scalar results. For the LBB condition we will essentiallybenefit from the H(div)-conformity and the interpolation properties presented in 2.2.5.First, let us define discrete HDG norms analogously to the scalar case (see section 1.2.3.2):

|||u|||21,∗ :=∑T∈Th

‖∇ u‖2

T + 1h

[[[ut]]]2∂T

(2.3.24)

|||u|||21 :=∑T∈Th

‖∇ u‖2

T + 1h

[[[ut]]]2∂T + h‖∂u∂n‖2∂T

(2.3.25)

They are again equivalent on the discrete Space Sh which follows with the same line of argument asin the scalar case. In contrast to the scalar case we don’t include a problem-dependent scaling like εor ν which is owed to the different type of analysis we use in this chapter. Scalings with ν now appearexplicitely in the estimates.

Proposition 2.3.4 (Coercivity). For a shape regular mesh and τhh (with h the local mesh size)sufficiently large Bh(·, ·) is coercive on Sh, that is

Bh(u,u) ≥ cν|||u|||21,∗ ≥ αBhν|||u|||21 ∀ u ∈ Sh (2.3.26)

with c, αBh ∈ R independent of the mesh size.

Proof. The proof can be straightforwardly taken from the proof of the scalar equivalent (Proposition1.2.2). As the norms are accordingly defined all inequalities hold even with the same constants. Eventhe inverse inequality can be applied component-wise which leads to the same overall constant as inthe scalar case.

As the claim and the proof follow directly from the scalar case, even the stabilization parameter τhcan be chosen in the same way as already discussed in chapter 1.

Proposition 2.3.5 (discrete LBB-condition for Dh). For Dh(·, ·) there holds

supu∈Sk+1

h

Dh(u, q)|||u|||1

≥ αDh‖q‖L2 ∀ q ∈ Qkh (2.3.27)

(as long as ΓD 6= ∂Ω (see also the subsequent Remark 2.3.3)).

Proof. In short we use several relations to come to

supu∈Sk+1

h

Dh(u, q)|||u|||1

≥ Dh(v, q)|||v|||1

2.≥ Dh(v

∗, q)c2‖v∗‖H1

1.≥ 1c1c2

‖q‖2L2

‖q‖L2≥ 1c1c2‖q‖L2 (2.3.28)

The idea of the inequality labeled with ”1.“ is to construct an [H1]d-regular function v∗ which fulfillsdiv(v∗) = q. In 2. we find a discrete function v ∈ Sk+1

h which fulfills Dh(v, q) = Dh(v∗, q) and canbe bounded in the HDG norm by the H1(Ω) norm of v∗. This procedure can also be found in [HL02].We will explain the steps 1. and 2. in more detail now:


1. Choose v∗, s.t. v∗ = ∇φ with φ ∈ H2 the solution of

−∆φ = q on Ω and ∂φ

∂n= 0 on ∂Ω (2.3.29)

The solution exists under appropriate regularity assumptions.Then we use

(a) Dh(v∗, q) =∑T∈Th

∫T−div(v∗)︸︷︷︸−∆φ=q

q dx =∫

Ωq2 dx = ‖q‖2

L2

(b) ‖v∗‖H1 ≤ ‖φ‖H2 ≤ c1‖q‖2L2 due to elliptic regularity

2. Now we use Lemma 2.3.1 which garantuees the existence of a discrete function v ∈ Sk+1h , s.t.

(a)Dh(v∗, q) = Dh(v, q) ∀ q ∈ Qk

h (2.3.30)

(b)|||v|||1 ≤ c2‖v∗‖H1 (2.3.31)

Remark 2.3.3 (Non-Uniqueness of the pressure if ΓD = ∂Ω):Note, that if Dirichlet boundary conditions are prescribed on the whole boundary, Dh(v, p) = 0 forany constant pressure p. Thus the inf-sup condition does not work in this case. Then, the proofgoes wrong as the Neumann problem (2.3.29) posed therein does not fulfill the compability condition∫Ω q dx =

∫∂Ω

∂φ∂nds = 0 and thus has no solution.

The problem we observe here, is owed to the lack of uniqueness of the pressure of the Stokes problem.This problem obviously also exists on the continuous level (as p only appears as a ∇p in the strongform).We can solve this issue by forbidding constant-pressure contributions to p, i.e. replacing L2(Ω) byL2

0(Ω) := p ∈ L2Ω :∫Ω p dx = 0 in the weak formulation. Then, the pressure is made unique again.

This also solves the problem with the inf-sup-condition, because if we replace Qkh by Qk

h \ R in thediscrete formulation, the compability condition for the Neumann problem (2.3.29) is fulfilled again andthe proof works.

2.3.2.3 Boundedness

Proposition 2.3.6 (Boundedness). For all u,v ∈ S and q ∈ Q there holds:

|Bh(u,v)| ≤ βBh ν |||u|||1 |||v|||1 (2.3.32)

with βBh = supx∈Fh(1 + τhh) and

|Dh(u, q)| ≤√d︸︷︷︸

= βDh

‖q‖L2 |||u|||1 (2.3.33)

Proof. (2.3.32) follows from the scalar poisson equation, where we used Proposition 1.2.3.With |div(u)| = |tr(∇ u)| ≤

√d‖∇ u‖F and Cauchy Schwarz inequalities the second equation also

holds true.



Proposition 2.3.7 (Approximation of the divergence-free subspace). Let u ∈ [H1(Ω)]d ∩ [Hm(Th)]dbe divergence-free, then the spaces S∗h and Sh fulfill the following approximation relation

infwh∈S∗h

|||wh − u|||1 ≤ C infvh∈Sh

|||vh − u|||1 (2.3.34)

with C independent of the mesh size. This means that the approximation error on the divergence-freesubspace is (up to a constant) as small as the approximation error on the whole space.

Proof. This proof is taken and adapted from [Brae97, remark 4.10]. As vh ∈ Sh there holds (seeLemma 2.3.1) ΠS

h(vh) = vh and for the solution u there holds ΠSh(u) ∈ S∗h. If we additionally make

use of the boundedness of the projection ΠSh we obtain:

|||ΠSh(u)−u|||1 = |||ΠS

h(u−vh)−(u−vh)|||1 ≤ |||u−vh|||1+c‖u−vh‖H(div) ≤ (1+c)|||u−vh|||1 (2.3.35)

So for every arbitrary vh ∈ Sh we can construct a function wh = ΠSh(u) ∈ S∗h, s.t. the claim holds.

Proposition 2.3.8. On a shape regular mesh Th which consists of affine transformed elements and afunction u ∈ [H1(Ω)]d ∩ [Hm(Th)]d, m = k+ 1 ≥ 2 there holds the following approximation result forSkh

infvh∈Skh

|||u− vh|||21 ≤ C∑T∈Th

h2kT |u|2Hm(T ) (2.3.36)


infvh∈Skh

|||u− vh|||1 ≤ C hk |u|Hm(Th) (2.3.37)

Proof. Use Bramble-Hilbert Lemma with L = idh−ΠSh,k and a scalar product similar to the one used

in the scalar case. Then, as in the scalar case we directly get the local estimate:

|||(idh −ΠSh,k)u|||T1 ≤ C |u|Hm(T ) ∀ u ∈ [Hm(T )]d (2.3.38)

As the scaling of the norm is exactly the same as in the scalar case an adapted version of the proof ofProposition 1.2.6 also works in the vector-valued case.

As already pointed out we could replace Σkh = BDMk with Σk

h = RTk−1 to get a stable solution. Butwe would loose one order of accuracy in the last Proposition if we’d do so.

Proposition 2.3.9. On a shape regular mesh Th which consists of affine transformed elements anda function p ∈ Hm(Th), m = k + 1 ≥ 1 there holds the following approximation result for Qk

h asintroduced in section 1.2.2.1.

infqh∈Qkh

‖p− qh‖2L2 ≤ C

∑T∈Th

h2k+2T |p|2Hm(T ) (2.3.39)


infqh∈Qkh

‖p− qh‖L2 ≤ C hk+1 |p|Hm(Th) (2.3.40)


Proposition 2.3.10. On a shape regular mesh Th which consists of affine transformed elements anda function (u, p) ∈ [H1(Th)]d ∩ [Hm+1(Th)]d × Hm(Th), m = k + 1 ≥ 1 there holds the followingapproximation result

infqh∈Qkh

‖p− qh‖2L2 + ν inf

vh∈Sk+1h

|||u− vh|||21 ≤ Chk+1(|p|2Hm(Th) + ν |u|2Hm(Th)

)(2.3.41)


Lemma 2.3.11 (Cea’s Lemma for Stokes). Let (u, p) ∈ Z be the solution of (2.3.1) and (uh, ph) ∈Zk+1h the solution of (2.3.10). Then there holds

‖p− ph‖L2 +√ν|||u− uh|||1 ≤ c

inf

qh∈Qkh‖p− qh‖L2 +

√ν inf

vh∈Sk+1h

|||u− vh|||1

(2.3.42)

for c independent of the mesh size.

Proof. We basically use Cea’s Lemma for Kh and again substitute U := (u, p) and V := (v, q) anddefine Uh, Vh,Wh ∈ Zk+1

h accordingly . Due to Brezzi’s Theorem, Lemma 2.3.2, we have inf-supstability for Kh with a constant αKh . Boundedness is easily achieved: For all Uh, Vh ∈ Zk+1

h :

|Kh(Uh, Vh)| ≤ |Bh(uh,vh)| + |Dh(uh, qh)| + |Dh(vh, ph)|≤ βBhν|||uh|||1|||vh|||1 + βDh‖ph‖L2 |||vh|||1 + βDh‖qh‖L2|||uh|||1≤ max(βBh , βDh)︸︷︷︸

=βKh

· (√ν|||uh|||1 + ‖ph‖L2) · (

√ν|||vh|||1 + ‖qh‖L2)

= βKh ‖Uh‖Kh ‖Vh‖Kh(2.3.43)

Now again we use the triangle inequality for the error U − Uh

‖U − Uh‖Kh ≤ ‖U − Vh‖Kh + ‖Vh − Uh‖Kh (2.3.44)

The second term can be bounded by the first term due to inf-sup stability, boundedness and consistencyof Kh. We take Wh such that inf-sup stability holds and Vh is an arbitrary element of Zk+1

h :

‖Uh − Vh‖Kh‖Wh‖Kh ≤ 1αKhKh(Uh − Vh,Wh)

≤ 1αKh

(Kh(Uh − U,Wh)︸︷︷︸

=0

+ Kh(U − Vh,Wh))

≤ βKhαKh‖U − Vh‖Kh‖Wh‖Kh

(2.3.45)

So we have‖U − Uh‖Kh ≤ (1 + βKh

αKh)‖U − Vh‖Kh (2.3.46)

which holds for all discrete functions Vh. So the limit is also true and if we keep in mind the definitionof the ‖ · ‖Kh norm the claim holds true with c = (1 + βKh

αKh). Note that the constant scales in the

same manner as for the scalar HDG discretization.


Lemma 2.3.12 (An estimate in the ||| · |||1 norm). Let Th be a quasi-uniform shape regular mesh,(uh, ph) ∈ Zk

h be the solution of (2.3.10) and (u, p) ∈ [H1(Ω)]d ∩ [Hm(Th)]d × Hm−1(Th), m ≥ 2the solution of (2.3.1). Then there holds the following error estimate:

‖p− ph‖L2 +√ν|||u− uh|||1 ≤ Chs

(√ν|u|Hs(Th) + |p|Hs−1(Th)

)s = min(k,m− 1) (2.3.47)

Proof. Marriage Lemma 2.3.11 and the approximation results, i.e. Propositions 2.3.8 and 2.3.9.

Lemma 2.3.13 (Ceá-like Lemma for the velocity). Let (u, p) ∈ Z be the solution of (2.3.1) and(uh, ph) ∈ Zh the solution of (2.3.10). Then there holds the following estimate not including thepressure field:

|||u− uh|||1 ≤ c1 infwh∈S∗h

|||u−wh|||1 ≤ c2 infwh∈Sh

|||u−wh|||1 (2.3.48)

with c1 and c2 independent of the mesh size.

Proof. Use coercivity, consistency (Bh(uh − u,vh) = −Dh(vh, p− ph)) and boundedness to get

αBh ν |||uh −wh︸︷︷︸=vh

|||21 ≤ Bh(uh−wh,vh) = Bh(uh − u,vh)︸︷︷︸=−Dh(vh,p−ph)=0

+Bh(u−wh,vh) ≤ βBhν|||u−wh|||1|||vh|||1

(2.3.49)where Dh(vh, p − ph) = 0 just holds true on the (exactly) divergence-free subspace. Divide (2.3.49)by αBh |||uh −wh|||1 and use the triangle inequality to get:

|||u− uh|||1 ≤ |||u−wh|||1 + |||wh − uh|||1 ≤(

1 + βBhαBh

)|||u−wh|||1 (2.3.50)

which completes the proof.

With the last estimate we see that the approximation of the pressure field p is not at all required tofulfill any approximation condition explicitly to get optimal convergence results in the |||·|||1 norm for thevelocity field u. Nevertheless exactly divergence-free velocities are essential for the last estimate andso the pressure field has to be chosen rich enough. For the reduced basis presented before (see section2.2.4.2), this means that only a constant pressure approximation (corresponding to the RT0 DOF)works without destroying higher order convergence of the velocity field as long as the incompressibilityconstraint is ensured exactly (e.g. with the Finite Element space presented before).

We can show better results for the velocity field in the L2 norm, if we use standard duality techniques.Therefore we have to assume regularity of the adjoint problem, which coincides with the primal problemexcept for the different force term f ∈ L2. So, those assumptions are as realistic as in the estimatesbefore.

Lemma 2.3.14 (Aubin Nitsche for the velocity). Let (u, p) ∈ Z be the solution of (2.3.1) and(uh, ph) ∈ Zh the solution of (2.3.10). Then there holds

‖u− uh‖L2 ≤ c h|||u− uh|||1 (2.3.51)

for c independent of the mesh size.


Proof. We pose the adjoint stokes problem with f(v) := (u− uh, v)Ω. Find W ∈ Z, s.t.

Kh(V,W ) = (u− uh, v)Ω ∀ V ∈ Z (2.3.52)

Under regularity assumptions on the domain we may assumeW ∈ [H2(Ω)]d×H1(Ω). Then we chooseV = U − Uh = (u− uh, p− ph) and so get

‖u− uh‖2L2 = Kh(U − Uh,W ) (∗)= Kh(U − Uh,W −Wh)

= Bh(u− uh, w − wh) +Dh(w − wh, p− ph)︸︷︷︸=0

+Dh(u− uh, q − qh)︸︷︷︸=0

≤ βBhν|||u− uh|||1|||w−wh|||1

(2.3.53)

where we shifted an arbitrary function Wh ∈ Z∗ into the bilinearform due to consistency in (∗).Dh(·, ·) = 0 holds as u, uh, w and wh are divergence-free. Now we can choose Wh as the linearprojection of the velocity (which is divergence-free) and the constant projection of the pressure. Aswe have assumed sufficient regularity we have

|||w−wh|||1 ≤ Cappr h |w|H2 (2.3.54)

Furthermore the | · |H2 semi-norm of u can be bounded by u − uh in the L2 norm since we assumedsufficient regularity:

ν|w|H2 ≤ Creg‖f‖L2 = Creg‖u− uh‖L2 (2.3.55)

Putting together (2.3.54) and (2.3.55) and dividing by ‖u − uh‖L2 the statement holds true withc = βBh Creg Cappr.

We conclude the a priori error analysis section with the statement of optimal convergence rate for thevelocity L2 error:

Lemma 2.3.15 (L2 norm estimate). Let Th be a quasi-uniform shape regular mesh, (uh, ph) ∈ Zkh

be the solution of (2.3.10) and (u, p) ∈ [H1(Ω)]d ∩ [Hm(Th)]d × Hm−1(Th), m ≥ 2 the solution of(2.3.1). Then there holds the following error estimate:

‖u− uh‖L2 ≤ Chs+1(√

ν|u|Hs(Th) + |p|Hs−1(Th))

s = min(k,m− 1) (2.3.56)

Proof. Use Lemma 2.3.13, Lemma 2.3.14 and Proposition 2.3.8

2.4. DIVERGENCE-FREE HYBRID DG FOR OSEEN EQUATIONS 83

2.4 Exactly divergence-free Hybrid DG method for the Oseenproblem

In this section we want to add a vector-valued version of the scalar HDG upwind method to the StokesEquations to include convective effects. We therefore start by replacing the nonlinear convectivevelocity u of the Navier-Stokes equations by a known velocity w. Normally w is an old value of uduring an iterative procedure. So the Oseen problem is

div(−ν∇ u+ u⊗ w + pI) = f in Ωdiv(u) = 0 in Ω


(2.4.1)


Except for the (for now linear) convective term div(u⊗w) the numerical treatment of the set of equa-tions is the same as for the Stokes Equations. We additionally introduce the bilinearform Ch(w; u,v)where w originates from a function w ∈ Z. We use the same procedure to derive the Hybrid DG up-wind bilinearform as in the scalar convective case by looking at each component one by one. The onlydifference appears when looking at the upwind-part. As un is continuous for each discrete function,we just have to treat the tangential component in an upwind fashion.We directly give all three versions of the bilinearform as the derivation is straight-forward when knowingthe scalar version and the notation used for the viscous bilinearform Bh.With uup = un +

utF wn < 0ut wn > 0

Ch(w,u,v)

=∑T∈Th

−∫Tu⊗ w : ∇ v dx+

∫∂Twnu

up v ds+∫∂Toutwn(utF − ut) vtF ds

(2.4.2a)

=∑T∈Th

∫Tdiv(u⊗ w) v dx+

∫∂Tin|wn|(ut − utF ) vt ds+

∫∂Tout|wn|(utF − ut) vtF ds

(2.4.2b)

=∑T∈Th

−∫Tu⊗ w : ∇ v dx+

∫∂Twnu

up [[vt]] ds+∫∂T∩Γwnu

tF vtF ds

(2.4.2c)

So the discrete version of (2.4.1) is

Find U = (u, p) ∈ Zh,D, s.t. Kh(U, V ) + Ch(w,u,v) = 〈f, v〉 ∀ V = (v, q) ∈ Zh,0 (2.4.3)

The bilinearform consisting of convection and viscosity(diffusion) is again called Ah : S × S → R:

Ah(u,v) := Bh(u,v) + Ch(w,u,v) (2.4.4)


The analysis used in chapter 1 can, for the most part, be translated directly to the Oseen problem toshow that Ah is consistent and coercive. Approximation results can be adopted from the last section


because we won’t introduce another convection norm, but work with the discrete broken Sobolev norm||| · |||1. So the analysis here will be most meaningful in the diffusion dominant regions, i.e. where thePeclét number is small. The only ingredient we have to put some effort into is the boundedness of thebilinearform Ch.

ConsistencyThe discretization of the vector-valued convection is also consistent:

Proposition 2.4.1 (Galerkin Orthogonality). Let (uh, ph) ∈ Zh,D be the solution of (2.4.3) and(u, p) be the solution of (2.4.1). Then there holds

Ch(w; uh − u,vh) = 0 ∀ vh ∈ Sh,0 (2.4.5)

Proof. See the proof of the scalar analogous Proposition 1.2.1.

StabilityIn this section we won’t consider purely convective equations and so we always consider just thebilinearform Ah directly. That’s why the used norms are taken from the diffusive limit (withoutviscous scaling), i.e. the discrete H1-seminorm ||| · |||1. To get stability, that is coercivity for this kind ofanalysis, is easy. From equation (1.3.18) we know that the contribution of the convective bilinearformis positive and so

Ah(u,u) = Bh(u,u) + Ch(w; u,u)︸︷︷︸≥0

≥ αBh|||u|||21 (2.4.6)

and we conclude:

Proposition 2.4.2 (Coercivity). For a shape regular mesh and τhh (with h the local mesh size)sufficiently large Ah(·, ·) is coercive on Sh, that is

Ah(u,u) ≥ αAh|||u|||21 ∀ u ∈ Sh (2.4.7)

with αAh = αBhν ∈ R only depending on the shape regularity of the mesh, the polynomial degree andthe kinematic viscosity ν.

BoundednessWe now want the trilinearform Ch to be bounded (in all three arguments) w.r.t. the discrete norm||| · |||1. Unfortunately we can show this only in two dimensions:

Proposition 2.4.3. For u,w ∈ (Σ ∩ [H1(Th)]2)× F t, div(w) = 0 and v ∈ Sh there holds

|Ch(w,u,v)| ≤ β∗Ch|||w|||1|||u|||1|||v|||1 = βCh|||u|||1|||v|||1 (2.4.8)

for β∗Ch ∈ R and βCh = β∗Ch|||w|||1independent of the mesh size.

Proof. See appendix, B

2.4. DIVERGENCE-FREE HYBRID DG FOR OSEEN EQUATIONS 85

Remark 2.4.1:The assumption u,w ∈ (Σ∩ [H1(Th)]d)× F t is satisfied for discrete functions of Sh as well as for allfunctions of [H1(Ω)]d if we again use the tangential trace operator to identify a function of [H1(Ω)]dwith one of F t. Thus the discrete and the continuous solution can be plugged in here for u and w.Actually for w the tangential component is not important at all, s.t. the regularity and the used normfor w could be optimized further and a wider class of functions (e.g. functions without facet-values)could be used. We indicate that by not using the bold writing for the first argument of Ch. However,in our case w is going to be an approximation to u, s.t. using the same regularity assumptions andnorms is reasonable.

Putting it all togetherAs we now still have coercivity, boundedness and consistency with respect to the ”old“ norm ||| · |||1, wecan use the results we got for the Stokes problem. Thus optimal convergence in the broken H1 normas well as L2 norm for the velocity field are achieved. Nevertheless there remains a disadvantage ofthe analysis we used here:If we look at the proof of Lemma 2.3.11 or Lemma 2.3.13 we have to replace the constants for Bhwith the constants for Ah. But this results in(

1 + βAhαAh

)=(

1 +νβBh + |||w|||1β∗Ch

ναBh

)= c1 + c2 ·

|||w|||1ν

for c1, c2 ∈ R (2.4.9)

Thus for convection dominating over viscosity (i.e. diffusion) the constant increases. That’s why thiskind of analysis, which is based upon ellipticity of the bilinearform Ah is less appropriate for largeconvection. Nevertheless existence and uniqueness of the discrete solution is still ensured in that case!Actually the method itself works fine also in the convection dominated case. As long as viscositydominates the behaviour of the system the constants which appear have the same magnitude as inthe Stokes case and the results are comparable also w.r.t. the constant.


2.5 Exactly divergence-free Hybrid DG method for the Navier-Stokes problem

If we replace the convective velocity of the Oseen problem we recover the original steady Navier-Stokesproblem:

div(−ν∇ u+ u⊗ u+ pI) = f in Ωdiv(u) = 0 in Ω


(2.5.1)

This set of nonlinear equation can then be solved iteratively by means of a Picard iteration, whichsolves an Oseen problem with the convective velocity of the last iteration at each step. Here thesolution of the Stokes problem is a good initial guess as it garantues a divergence-free velocity.The discrete problem reads

Find U = (u, p) ∈ Zh,D, s.t. Kh(U, V ) + Ch(u,u,v) = 〈f, v〉 ∀ V = (v, q) ∈ Zh,0 (2.5.2)


The only step we have to take to get from the discretization of the Oseen problem to the nonlinearproblem is to describe an iterative procedure converging to the discrete solution. E.g. this iterativeprocedure could be a Newton-Raphson method. Here, we consider the simpler approach, the Picard-iteration:

Algorithm 2.5.1 Picard-iteration for the steady Navier-Stokes equationsSet i = 0Solve the discrete steady Stokes problem (2.3.10) for u0

repeatincrease i by 1Solve the discrete steady Oseen problem (2.4.3) for un with w = un−1

until desired accuray is achieved (‖un − un−1‖L2 ≤ ε)

Now we hope that this fixpoint iteration converges. If it does, it obviously converges to a discretesolution of (2.5.2). Convergence of the fixpoint iteration, which would imply existence and uniquenessof (2.5.2), is only ensured under so called ”smallness assumptions”, i.e. volume forces and boundaryconditions are assumed to be sufficiently small. To use those smallness assumptions is reasonable asexistence and uniqueness of the continuous problem also depends on them. In [CKS05, Theorem 4.7]and [CKS05, Theorem 4.8] those smallness assumptions where used to proof existence and uniquenessof the discrete solution for a similar exactly divergence-free DG method. They also provide a priorierror estimates which ensure the optimal convergence properties we achieved for the Stokes problemand the Oseen problem. The same technique could be applied here. A detailed look at those proofsshows that except for adapting the bilinearforms and norms used therein the proof can actually bedirectly applied to our discretization.

2.5. DIVERGENCE-FREE HYBRID DG FOR NAVIER-STOKES EQUATIONS 87

2.5.2 Numerical Example: Kovasznay flow

We consider an exact solution to the Navier-Stokes equations obtained by Kovasznay [Kova48] whichwe will use as a reference solution. The Kovasznay flow is two-dimensional and describes the fluid flowbehind a grid. The solution is given by

u(x, y) =(

1− eλx cos(2πy)λ2πe

λx sin(2πy)

)

p(x, y) = −12e

2λx + p with p ∈ R

withλ = −8π2

ν−1 +√ν−2 + 64π2

On the boundary of the domain Ω = [−12 ,

32 ]× [0, 2] we prescribe inhomogeneuos Dirichlet data given

by the exact solution. If we additionally demand the average of p ( 1|Ω|∫

Ω p dx) to be zero the solutionof the steady Navier-Stokes equations is unique and is exactly the solution proposed by Kovasznaywith a suitable p.Again we choose τh as in the numerical examples before with α to be 2 and consider the case whereν = 1. The chosen meshes we used and a view on the solution can be found in Figure 2.5.1.

(a) Initial and finest mesh used for the discretiza-tion and absolute value of the velocity solution

(b) Absolute value of the velocity solution andstreamlines

Figure 2.5.1: Solution and meshes for the Kovasznay problem

Convergence of the method

In the next table we see the convergence of the method for the chosen example. We observe thatthe pressure error in the L2-norm and velocity error eh in the broken H1 seminorn converge with theexpected rate hp whereas the L2-norm of the velocity error eh converges even an order faster. That iswhat we expected from the analysis carried out before. We can furthermore observe how the reductionof the Finite Element space reflects on the number of unknowns (compare DOF with DOF∗) and the


number of nonzero elements in the matrices (compare nze, which stands for ”nonzero entries“, withnze∗). For p = 2 the unknowns are reduced roughly by a factor of 1

4 and a factor of nearly 12 is achieved

for p = 4. We see that the asymptotial reduction factor of 12 is already achieved for moderately high

p. For all meshes and polynomial degrees the iteration counts are bounded by 10.

meshlevel ‖eh‖L2 ‖div(uh)‖L2 ‖∇(eh)‖L2 ‖p− ph‖L2

(elements) p DOF DOF∗ error order error error order error order0 (18) 2 306 234 3.21 — 6.77e-9 5.76e1 — 4.19e1 —1 (72) 1 152 864 6.25e-1 2.36 5.48e-9 1.91e1 1.60 2.16e1 0.962 (288) 4 464 3 312 8.62e-2 2.86 5.20e-9 5.22 1.87 6.52 1.733 (1 152) 17 568 12 960 1.00e-2 3.11 5.69e-9 1.31 1.99 1.75 1.904 (4 608) 69 696 51 264 1.17e-3 3.10 7.83e-9 3.26e-1 2.01 4.49e-1 1.965 (18 432) 277 632 203 904 1.42e-4 3.04 1.28e-8 8.11e-2 2.01 1.13e-1 1.99

Table 2.2: Errors and order of convergence of pressure and velocity in different norms as well as acomparison between DOF of the hybrid divergence-free DG method with and without the reducedbasis for the kovasznay flow problem using polynomial spaces up to degree 2.

meshlevel ‖eh‖L2 ‖div(uh)‖L2 ‖∇(eh)‖L2 ‖p− ph‖L2

(elements) p DOF DOF∗ error order error error order error order0 (18) 4 780 456 3.24e-1 — 2.70e-9 9.36 — 3.8e1 —1 (72) 3 000 1 704 1.56e-2 4.38 4.86e-9 8.59e-1 3.45 1.22 4.962 (288) 11 760 6 576 5.62e-4 4.79 5.03e-9 6.06e-2 3.83 1.09e-1 3.483 (1 152) 46 560 25 824 1.84e-5 4.93 5.80e-9 3.89e-3 3.96 7.83e-3 3.804 (4 608) 185 280 102 336 5.73e-7 5.01 1.88e-9 2.43e-4 4.00 5.08e-4 3.955 (18 432) 739 200 407 424 1.80e-8 4.99 7.81e-8 1.52e-5 4.00 3.19e-5 3.99

Table 2.3: Errors and order of convergence of pressure and velocity in different norms as well as acomparison between DOF of the hybrid divergence-free DG method with and without the reducedbasis for the kovasznay flow problem using polynomial spaces up to degree 4.

89

Chapter 3

Time Integration for incompressible flowproblems

3.1 Semi-discrete form of convection diffusion type problems

So far we only considered steady flow problems in the last chapters in order to concentrate on theproposed spatial discretization. Now we to take time discretization into consideration. We will usethe method of lines here, i.e. we use spatial discretization and apply time discretization schemesafterwards. As we discussed convection diffusion and the Navier-Stokes Equations we want to castboth problems into a unified form. The discretization of the unsteady problems (1.4.1) and (2.5.1)then read

M

(∂u

∂t, v

)+B(u, v) + C(u;u, v) = f(v) ∀ v ∈ V (3.1.1)

or in operator notation:

M∂u

∂t+Bu+ C(u)u = f (3.1.2)

In the following table we summarize the interpretations of this notation which result either in theunsteady convection diffusion equation or the Navier-Stokes equations.

unified notation scalar convection diffusion counterpart Navier-Stokes counterpartu, v ∈ V scalar quantity u = (u, uF ) ∈ Vh velocity-pressure pairU=(u, uF , p)∈ZhM(·, ·) Mass matrix M(u, v) :=

∫Ωu v dx Mass matrix M(u, v) :=

∫Ωu v dx

matrix properties block diagonal not block diagonalB(·, ·) diffusion bilinearform Bh(·, ·) (1.2.8) stokes bilinearform, i.e. viscosity, pres-

sure and incompressibility constraint:Kh(·, ·) (2.3.11)

matrix properties symmetric positive definite symmetric indefiniteC(·; ·, ·) convection bilinearform Ch(·, ·) (1.3.4) nonlinear vector-valued convection bi-

linearform Ch(·; ·, ·) (2.4.2)properties linear nonlinearf(·) source terms force terms

90 CHAPTER 3. TIME INTEGRATION FOR INCOMPRESSIBLE FLOW PROBLEMS

We see that the bilinearforms for the Navier-Stokes equations are less comfortable as those of theconvection diffusion equation. Actually due to the incompressibility constraint the system of equationsarising from spatial discretization is not even an ordinary system of equations, but a differential algebraic(DAE) system. That is why we will concentrate on this more involved case, transferring the results toconvection diffusion is then straight forward.

Remark 3.1.1 (On facet unknowns):Time derivatives of facet unknowns do not appear in either of both formulations and thus even theconvection diffusion discretization is not an ordinary system of equations (the equations for vF arealgebraic). But as facet unknowns can be eliminated easily, see section 1.2.2.2, 1.4.1.1 and 1.5.1, thisproblem is only a formal one and can be overcomed easily. In the proceeding we will neglect to adressthis issue.

3.2 Approaches for Time Integration: explicit vs. implicit

We now want to discuss the numerical integration of (3.1.2). Therefore let’s have a look on the twomost famous lowest order time integration methods, the explicit (forward) and the implicit (backward)Euler method. Denoting with un, un+1 the values of the discrete solutions at time tn and tn+1,respectively, and furthermore defining the increment ∆u := un+1 − un, the explicit Euler methodreads:

M∆u∆t = fn −Bun − C(un)un (3.2.1)

whereas the implicit Euler method is

M∆u∆t = fn+1 −Bun+1 − C(un+1)un+1 (3.2.2)

Those methods are the simplest representatives for explicit and implicit time integration methodswhose advantages and disadvantages will be discussed in the following.

3.2.1 Explicit methods

For the Navier-Stokes equations fully explicit methods are not directly possible as the incompressibilityconstraint is algebraic and thus has to be treated implicitely. At least you need to treat the pressureand the incompressibility constraint implicitely and the rest explicitely.A more elegant way to derive explicit methods for Navier-Stokes equations is the projection onto thedivergence-free subspace. On the divergence-free subspace, the Navier-Stokes equations are a systemof ordinary differential equations. For this we can apply every explicit method directly. Afterwards wehave to go back on the full space by reintroducing the incompressibility constraint for the new timestep, i.e. an implicit term. As this is nevertheless the minimalistic approach for time integration ofNavier-Stokes equations we may still call it explicit.Explicit time integration comes with basically two nice properties: Discrete systems arising from explicittime integration are linear and normally involve only the mass matrix. For Navier-Stokes equations

3.2. APPROACHES FOR TIME INTEGRATION: EXPLICIT VS. IMPLICIT 91

the situation is - as we have seen - more involved and the mass matrix has to be extended with theincompressibility constraint:(

M 00 0

)∆u∆t = ... −→

(M DT

D 0

)∆u∆t = ...

where M is the mass matrix and DT , D the discrete gradient of p and the discrete incompressibilityconstraint, respectively.

Remark 3.2.1 (Modifications of the spatial discretization for explicit time integration):Solving linear systems of equations with the extended mass matrix(

M DT

D 0

)

for the discretization proposed in chapter 2 is not at all cheap (mainly due to the fact that M−1 isnot sparse), at least we can do better. Therefore let’s shortly explain possible modifications of theformulation which make the use of the spatial discretization more efficient for explicit time integration:If we use completely discontinuous finite elements instead of H(div)-conforming ones and introducelagrange multipliers pF , qF to ensure normal-continuity, the mass matrix M is block diagonal and theschur complement DM−1DT is a sparse symmetric positive definite matrix which is comparable toa discrete laplacian for the pressure. As H(div)-conformity is then ensured through the additionallagrange multipliers the spatial discretization coincides with the one proposed in chapter 2.

The major disadvantage of explicit methods is that they are only conditionally stable, s.t. a certaintime step restriction has to be fulfilled. Due to the time step restriction this approach renders inefficientfor non-negligibly large viscosity. Those time step restrictions behave in the following way:Convection owes time step restrictions which scale linearly with the spatial resolution and the magnituteof the convective term |u|

∆t ≤ cCh

p21|u|

(3.2.3)

for some constant cC ∈ R whereas the viscosity demands more:

∆t ≤ cBh2

p41ν

(3.2.4)

for some constant cB ∈ R. The latter condition is often the more disturbing one as it scales quadrat-ically with h

p2 and dominates over the first condition for moderate values of ν.

3.2.2 Implicit methods

In contrast to explicit methods most practically relevant implicit methods are unconditionally stable.This of course is desireable, but it doesn’t come for free. The system of equation you have to solvein each time step may be considerably expensive: In our case it is nonlinear and each linearizationresults in a linear system of equations with a nonsymmetric and indefinite matrix. This implies thatthe matrices are, in contrast to explicit methods, varying at each time step.For stiff systems implicit methods are indispensable and the viscosity matrix of the Navier-Stokesdiscretization is considerably stiff for moderate values of ν which results in the inconvenient time steprestriction (3.2.4) for explicit methods. Thus for problems with non-small diffusity (viscosity) implicitmethods should be prefered.


3.2.3 Semi-Implicit methods

Let’s also discuss a mixture which combine the features of explicit and implicit methods with thefollowing aim:

• Avoid solving nonlinear, nonsymmetric equations (too expensive)

• Avoid quadratical dependence of the time step on the mesh size (see (3.2.4))

or, to put the second point in another way:

• Avoid integration of the stiff part with an explicit scheme

We achieve this by splitting the r.h.s. of (3.1.2) into an implicit and an explicit part. The diffusiveterm, i.e. B is treated implicitely and the (nonlinear) convective term C is treated explicitely resultingin the first order semi-implicit Euler :

M∆u∆t = fn −Bun+1 − C(un)un

⇔ (M + ∆tB)∆u = ∆t(fn −Bun − C(un)un)(3.2.5)

As f is (assumed to be) independent of u the choice of the time level for f is arbitrary, i.e. the costsfor evaluation are the same for tn and tn+1. In the following we will treat f explicitely, i.e. in the sameway as the convective term. Nevertheless other choices are possible.

3.2.4 Stability vs. Efficiency

In practice all three approaches discussed above are use. Which of those methods is the most appro-priate, i.e. the more efficient, normally depends on the spatial discretization, the flow parameters andperhaps also on the computer you are working on1. Implicit methods are most robust and could beused for all practically relevant flow configurations, but if the use of explicit or semi-implicit methodsis possible without dramatical restrictions of the time step, the implicit method is often (far) moreexpensive. In general there can’t be stated a best choice. Nevertheless as (fully) explicit and (fully)implicit approaches are well-known we will concentrate on the discussion of semi-implicit approaches.

3.3 Higher Order methods

All three time discretizations mentioned above ((3.2.1), (3.2.2) and (3.2.5)) are only of first order.For fully explicit (and with this respect the method with implicit treatment of the incompressibilitystill applies as fully explicit) and fully implicit schemes extensions to higher order methods by e.g.Runge-Kutta formulas are well-known. For semi-implicit schemes the extension to higher order seemsto be less famous. We want to discuss the so called implicit-explicit (IMEX) schemes which allow forsuch an extension in the next section.

1e.g. some methods may be easier to parallalize then other

3.4. IMEX TIME INTEGRATION 93

3.4 IMEX Time Integration

Starting point is the equation (3.1.2)

M∂u

∂t= f −Bu− C(u)u

where B is a linear operator describing a diffusive process and thereby is stiff, whereas C(u) is apossibly nonlinear convective operator. In the following we want to treat the diffusive operator alwaysimplicitely and the convective operator always explicitely. For convenience we thus write (rememberthat f is also treated explicitely)

M∂u

∂t+Bu = f − C(u)u

Terms on the right hand side are treated explicitely and those on the left hand side implicitely. Wealready saw the first order version of these IMEX schemes. It is obtained by replacing the timederivative by a finite difference and by exchanging u by un+1 on the l.h.s. and u by un on the r.h.s.resulting in

M∆u∆t +Bun+1 = fn − C(un)un

We will discuss two different realizations of higher order extensions in the sequel. These are eitherachieved with special partitioned Runge-Kutta methods or Multistep methods. Further we will focuson methods which are stiffly accurate as they garantuee that the solution at the new time level isdivergence-free. DAE systems can be seen, with this regard, as a limit case of stiff systems and sostiffly accurate schemes are appropriate for those kind of problems.

3.4.1 IMEX Runge-Kutta methods

The basic idea of IMEX Runge-Kutta methods is simple: combine two compatible Runge-Kuttaschemes from which one is implicit and one is explicit, in a way such that a scheme is achievedwhich has the overall order of accuracy which coincides with the minimal order of accuracy of theindividual Runge-Kutta schemes.To find the compatibility condition for two Runge-Kutta schemes to match in the desired sense isactually quite simple: They just have to work on the same time levels. Methods fulfilling these criteriawere introduced and discussed in [ARS97].In the following we will only consider the pair of Runge-Kutta schemes which consists of one explicitRunge-Kutta scheme and one diagonally implicit Runge-Kutta scheme which is stiffly accurate. Fur-thermore we prefer methods which are singly diagonally implicit and L-stable2. A method is calleddiagonally implicit if the according butcher tableau has only zero entries above the diagonal, i.e. eachstage value does not depend on future stage values. It is furthermore called singly diagonally implicitif the diagonal entries of the butcher tableau are all the same. If we have singly diagonally implicitmethods and the stiff part, i.e. the part we want to treat implicitely, is linear, the linear equations wehave to solve within one step use the same matrix. This allows for reusing factorizations and precondi-tioners. If the time step is additionally used for multiple time steps the l.h.s. matrix can be used evenmore often. We require stiffly accurate schemes to deal with the incompressibility-constraint. Eachnew time step has to fulfill it exactly and not only in an interpolated sense. Stiffly accurate schemes

2Thus all examples we present here are singly diagonally implicit and L-stable what does not hold in general for allmethods fitting into the presented framework.


can always be written, s.t. the last stage value of the Runge-Kutta scheme can be taken as the newinstance for the new time step tn+1.The described IMEX Runge-Kutta methods (diagonally implicit & stiffly accurate) work as showed inalgorithm 3.4.1.

Algorithm 3.4.1 Stiffly accurate s stage IMEX Runge-Kutta Time Integration1: Use Startvalue un2: Evaluate C1 = C(un)un − f(tn)3: for i = 1 to s do :4: Solve (M + ∆t ←a i,iB)ui = Mun −∆t

i∑j=1

→a i+1,jCj −∆t

i−1∑j=1

←a i,jBj for ui

5: Evaluate Bi = Bui

6: Evaluate Ci+1 = C(ui)ui − f(tn + ci∆t)7: end for8: Set un+1 = us

Coefficients with left-sided arrows refer to a Butcher tableau of a (diagonally) implicit Runge-Kuttascheme and those with a right-sided arrow are taken from a Butcher tableau of a compatible explicits-stage Runge-Kutta scheme.

Normally Butcher tableaus also include coefficients bj which describe which linear combination of stagevalues should be used for the new instance un+1. Here, the last stage value is automatically the newvalue at time level tn+1, s.t. the coefficients →as,j and ←as,j coincide with

→bj and

←bj, respectively. In

table 3.1 the general pair of Butcher tableaus for an s-stage IMEX Runge-Kutta schemes is shown.

explicit implicit0 0 0 0 . . . 0 0c1

→a2,1 0 0 . . . 0 0 ←a1,1 0 . . . 0 0c2

→a3,1→a3,2 0 . . . 0 0 ←a2,1

←a2,2. . . 0 0

... ... ... ... . . . ... ... ... ... . . . . . . ...cs−1

→as,1→as,2

→as,3 . . . 0 0 ←as−1,1←as−1,2 . . .

←as−1,s−1 0cs

→as+1,1→as+1,2

→as+1,3 . . .→as+1,s 0 ←as,1

←as,2 . . .←as,s−1

←as,sq q q . . . q q q q . . . q q→b1

→b2

→b3 . . .

→bs

→bs+1

←b1

←b2 . . .

←bs−1

←bs

Table 3.1: IMEX Runge-Kutta Butcher tableau consisting of an embedded pair of Butcher tableausfor compatible Runge-Kutta methods

The coefficients ci coincide for both single Butcher tableau, s.t. both Runge-Kutta methods workon the same time levels, thus both methods are compatible. We see on the left side the Butchertableau of the explicit Runge-Kutta method embedded in the IMEX Runge-Kutta Butcher tableau.The coefficients of the original scheme are extended by →as+1,j =

→bj. The embedding of the implicit

Runge-Kutta method is straight forward. Now the new Butcher tableau reflects the procedure ofalgorithm 3.4.1:In the first row we have only zeros, i.e. only the convective part at tn is evaluated. Then the next row


describes which linear combination of convection terms and diffusion terms (including the new stagevalue) should be used to calculate the new stage value. Continuing until stage s is reached gives usan approximation for u at tn+1. Therefore no additional linear combination of old stage values has tobe computed as ←as,j =

←bj due to the stiffly accuracy of the implicit method.

Note also that lines 5 and 6 of algorithm 3.4.1 don’t have to be executed for i = s.If the operator B is zero we recover the explicit Runge-Kutta method and for C = 0 we recover the(diagonally) implicit Runge-Kutta method.Let’s discuss the lowest order case s = 1. Here, we have exactly four free parameters, these are c1,→a1,1,

←a1,1 and b1. For consistency (of the explicit method) we need b1 = →a1,1 = 1. This correspondsto the use of the explicit euler for the explicit part. For our problems (f is not treated implicitely andno other term depends directly on t) the choice of cs is arbitrary (see algorithm 3.4.1). So the onlyimportant parameter we have is to choose is ←a1,1 =: θ. For θ ∈ (0, 1] we recover the implicit θ-schemefor the implicit part. For θ ≥ 1

2 the implicit method is unconditionally stable and at least of first order(for θ = 1

2 we recover the second order accurate Crank-Nicholson method). So both methods are offirst order and compatible. For θ = 1 we recover the semi-implicit Euler method which is completelydescribed with the following Butcher tableau:

Implicit Euler → Semi-Implicit Euler ← Explicit Euler

1 11 →

0 0 01 1 0 1

1 0 1← 0 0

1

Table 3.2: The semi-implicit Euler as a 1-stage IMEX Runge-Kutta scheme

We want to complete the presentation of IMEX Runge-Kutta methods by stating second and thirdorder accurate IMEX Runge-Kutta methods which use SDIRK (singly diagonally implicit Runge-Kutta)schemes. Both are taken from [ARS97].For second order accuracy it is sufficient to use a 2-stage IMEX Runge-Kutta scheme with the followingbutcher tableau:

((SA)SDIRK2): → (SA)RKIMEX2 ← ERK2

γ γ 01 1− γ γ

1− γ γ→

0 0 0 0γ γ 0 0 γ 01 δ 1− δ 0 1− γ γ

δ 1− δ 0 1− γ γ

←0 0 0γ γ 0

δ 1− δ

Table 3.3: 2-stage second order accurate IMEX Runge-Kutta scheme

To achieve third order we have to use at least s = 4 stages. This is due to the combination ofconsistency and compatibility conditions that have to be fulfilled. One possibility of a 4-stage thirdorder accurate IMEX Runge-Kutta scheme is given next.


0 0 0 0 0 01/2 1/2 0 0 0 0 1/2 0 0 02/3 11/18 1/18 0 0 0 1/6 1/2 0 01/2 5/6 −5/6 1/2 0 0 1/2 1/2 1/2 01 1/4 7/4 3/4 −7/4 0 3/2 −3/2 1/2 1/2

1/4 7/4 3/4 −7/4 0 3/2 −3/2 1/2 1/2

Table 3.4: 4-stage third order accurate Runge-Kutta IMEX scheme

In [KC03] IMEX Runge-Kutta methods with 6 stages which has overall accuracy order 4 is given aswell as an alternative for the third order accurate IMEX Runge-Kutta method. Both use an ESDIRKmethod which in contrast to the SDIRK methods we stated has one explicit first stage even for theimplicit method. This reduces the computational effort one step further3. Furthermore those methodsuse embedded Runge-Kutta formulas which allow for error estimation and adaptation as well as denseoutput. For readers with solid background in numerical methods for ordinary differential equations[KC03] also provides deeper insight into generalizations of IMEX schemes (to splitting of the operatorinto more than two parts) and order and coupling conditions and their derivations. In [KCGH07] thesemethods are used to solve compressible flow problems.

3.4.2 Multistep IMEX methods

Another possibility to generalize the first order semi-implicit Euler scheme to higher order accuracyis the use of multistep methods. The basic idea is to extrapolate information from old time stepsto achieve higher order accuracy. Thus the number of function evaluations and linear systems whichhave to be solved is reduced to one for both, which makes the implementation efficient and straightforward: Extrapolated values for the stiff, the nonstiff and the mass term replace the evaluations attime level tn of the semi-implicit Euler method and with the coefficient c∗ denoting the weight of thestiff term at time level tn+1 we obtain:

(M + ∆tc∗B)un+1 = Ms−1∑j=0

ajun−j + ∆t

s−1∑j=0

bjC(un−j) + ∆ts−1∑j=0

cjB(un−j)

The following table shows possible choices up to order 4 which are taken from [ARW95]. Notice thatSBDF methods can be generated generically of arbitrary order by left-sided finite differences, but asthe BDF method is only L-stable up to order 4, the SBDF method is less desirable for p > 4.

3Thus the alternative third order IMEX Runge-Kutta method only has to solve three systems of equations


Name s c∗ a0 a1 a2 a3 b0 b1 b2 b3 c0 c1 c2 c3TFE(θ) 1 θ 1 1 1− θ

BEFE, [1],SBDF1 1 1 1 1 0CNFE[1/2] 1 1/2 1 1 1/2

[[γ, c]] 2 γ+c/2γ+1/2

2γγ+1/2

−γ+1/2γ+1/2

γ+1γ+1/2 −

γγ+1/2

1−γ−cγ+1/2

c2γ+1

CNAB [[1/2, 0]] 2 1/2 1 0 3/2 −1/2 1/2 0MCNAB [[1/2, 1/8]] 2 9/16 1 0 3/2 −1/2 3/8 1/16

CNLF [[0, 1]] 2 1 0 1 2 0 0 1SBDF2 [[1, 0]] 2 2/3 8/3 −2/3 2 −1 0 0

SBDF3 3 6/11 18/11 −9/11 2/11 18/11 −18/11 6/11 0 0 0

SBDF4 4 12/25 48/25 −36/25 16/25 −3/25 48/25 −72/25 48/25 −12/25 0 0 0 0

Table 3.5: Some Multistep IMEX schemes (source: [ARW95])

Some explanations to the table: The first order 1-stage IMEX multistep methods can be parametrizedby one parameter θ. The first method is the ”Theta-Forward-Euler” (TFE) method for general θ.The ”Semi-implicit Backward Differentiation Formula“ of first order (SDBF1) which is exactly the”Backward-Euler-Forward-Euler“ (BEFE) is then again the ”Theta-Forward-Euler” method for θ = 1,whereas the “Crank-Nicholson-Forward-Euler”(CNFE) is the TFE method for θ = 1/2.The set of all second order 2-stage IMEX multistep methods can be parametrized by 2 parametersγ and c. Again, the general case ist presented first and then four special cases, namely the “Crank-Nicholson-Adams-Bashforth” (CNAB), the “Modified-Crank-Nicholson-Adams-Bashforth” (MCNAB),the “Crank-Nicholson-Leap-Frog” (CNLF) and the ”Semi-implicit Backward Differentiation Formula“of second order (SBDF2) follow. As mentioned before the SBDF method can be derived for arbitraryorder, but only the case s ≤ 4 is relevant for practical applications, as L-stability is lost for larger s.Although multistep methods are obviously computationally cheap (at the cost of a bit more memory tokeep old time values), you can not, in contrast to the Runge-Kutta methods presented before, expect again in the time step restriction. On the contrary: the time step restriction gets worse for higher ordermultistep methods. So the apparent efficiency is compensated. In our calculations IMEX Runge-Kuttamethods have been more robust and also a bit cheaper than IMEX Multistep methods.

98 APPENDIX A. NOTATION, ELEMENTARY INEQ. AND ORTHOGONAL POLYNOMIALS

Appendix A

Notation, elementary inequalities andorthogonal polynomials

A.1 Notation

In this section the notation used through out the diploma thesis is summarized.

General Abbreviationscharacteristic elementlength, e.g. diameter

: h

vector : vouter normal : ncompound function : v = (v, vF )compound vectorfunction

: v = (v, vtF )

velocity pressure pair : U = (u, utF , p)dyadic product : ⊗ (a⊗ b)i,j = aibjsmaller up to const. : greater up to const. : smaller and greater upto constant

: '

Degree(s) of freedom : DOFFinite Element Method : FEMContinuous Galerkin : CGDiscontinuous Galerkin : DGHybrid Disc. Galerkin : HDGHybrid mixed : HMCauchy Schwarz (ineq.) : C.S.

Constants Triangulationstability constant for bilinearform Ah : αAhboundedness constant for bilinearform Ah : βAhshape regularity constant on one element : cTshape regularity for the whole mesh : csr

Triangulation : Thset of all element boundaries : Fhset of interior el. interfaces : F inth

set of exterior el. interfaces : F exth

A.1. NOTATION 99

Norms and scalar products

For f, g ∈ L2(Ω) we use the following notation for the L2 inner product and the L2 norm:

(f, g)Ω :=∫

Ωf(x) g(x) dx

‖f‖2Ω :=

∫Ωf(x)2 dx

DG and HDG operators

As we consider having double and even triple-valued functions on interfaces, we have to distinguishbetween values which come from functions of different sides and functions which live only on theinterface. The interface functions are always indicated by the index F . When we work with interfaceunknowns, neighbouring unknowns don’t interact directly, so that for this case the following jumpoperator which is single-sided and an according norm is sufficient:

[[f ]] := f − fF[[f t]] := f t − f t

F

[[[f ]]]∂Ω := ‖f − fF‖∂Ω

If we are interested in average values and jump operators between neighbouring values without incor-perating interface unknowns 1, we define for a scalar function f and a vector function h - with n+,n− the outwards pointing normals of each side - the following average and jump operators:

f := 12(f+ + f−)

[[f ]]* := f+ − f−

h∗ := h+ n+ − h− n−

[[h]]n := h+ n+ + h− n−

All those operators are scalar values. The stars indicate that the signs of those values are direction-dependent2.

1like it is the case for Standard DG2depending on which side is marked with + and which with −

100APPENDIX A. NOTATION, ELEMENTARY INEQ. AND ORTHOGONAL POLYNOMIALS

Infinite dimensional spaces Finite Element spacesPk(T ) =

k∑i=0

∑l,m,n≥0l+m+n=k

al,m,nxlymzn, al,m,n ∈ R

Q = L2(Ω) Qkh =

u : u ∈ Pk(T ) ∀ T ∈ Th

F = L2(Fh) F k

h =uF : uF ∈ Pk(E) ∀ E ∈ Fh

F t =

utF : utF ∈ [L2(Fh)]d × n

F k,th =

uF : uF ∈ [Pk(E)]d × n ∀ E ∈ Fh

V =

(u, uF ) : u ∈ H2(Th) ∩H1(Ω), V k

h = Qkh × F k

h

uF ∈ L2(Fh)

VD = (u, uF ) ∈ V, uF = uD on ΓD Vh,D = (u, uF ) ∈ Vh, uF = uD on ΓDV0 = (u, uF ) ∈ V, uF = 0 on ΓD Vh,0 = (u, uF ) ∈ Vh, uF = 0 on ΓDΣ = H(div,Ω) Σk

h =u : u ∈ [Pk(T )]d ∀ T ∈ Th,

[[u]]n = 0 ∀ E ∈ Fh

= BDMk

S = Σ× F t Skh = Σkh × F

k,th

Z = [H1(Ω)]d ∩ [H2(Th)]d × L2(Ω) Zkh = Skh ×Qk−1

h

W = H1(Ω) W kh = u;u ∈ Pk(T ) ∀ T ∈ Th, [[u]]* = 0 ∀ E ∈ Fh

Projections

L2 projection on all facets of one element ∂T : Πfach,k

L2 projection on one element T : Πelh,k

L2 projection on Vh : Πh,k

H(div)-conforming projection on Σh : ΠΣh,k

L2 projection on F th : Πf,t

h,k

H(div)-conforming projection on Sh : ΠSh,k

A.2. ELEMENTARY INEQUALITIES 101

A.2 Elementary Inequalities

Cauchy Schwarz inequality

Let V be a vector space over R with the inner product (a, b)N for a, b ∈ V . Then

(a, b)N ≤ ‖a‖N‖b‖N

holds with ‖x‖N =√

(x, x)N . This inequality is useful for different norms on different vector spaces.In many proofs (especially boundedness proofs) we make use of two special cases of it:

a) In the first version we use the Cauchy Schwarz inequality for the L2 scalar product on Ω

(f, g)Ω :=∫

Ωf(x)g(x) dx

Thus for all f, g ∈ V = L2(Ω) there holds:

∫Ωf(x)g(x) dx = (f, g)Ω ≤ ‖f‖Ω‖g‖Ω =

√∫Ωf 2(x) dx

√∫Ωg2(x) dx (A.2.1a)

b) The second version appears often when sums (over 3 or 4 components or over the triangulation)are involved. The l2 scalar product is defined as

(a, b)2 =m∑i=1

ai bi

for a, b ∈ Rm. Thus for sums there holds

m∑i=1

ai bi = (a, b)2 ≤ ‖a‖2‖b‖2 =√√√√ m∑i=1

a2i

√√√√ m∑i=1

b2i (A.2.1b)

Young’s Inequality

Let a, b, γ ∈ R. Then there holds

0 ≤ (a− γb)2 = a2 − 2aγb+ (γb)2

s.t. with straight forward manipulations we get Young’s inequality:

⇒ ab ≤ a2

2γ + γb2

2 (A.2.2)

102APPENDIX A. NOTATION, ELEMENTARY INEQ. AND ORTHOGONAL POLYNOMIALS

Inverse Inequalities

At several places we will make use of the inverse inequality which bounds the L2 norm on the boundaryby the L2 norm at the domain for polynomials. Let u be a polynomial in x ∈ Rd, d = 1, 2, 3 with orderp on a polyhedra T . Then there holds

‖u‖∂T ≤c√h‖u‖T (A.2.3a)

with c ∈ R independent of the diameter of the polyhedra T .Especially for the normal derivatives of polynomials we will use this inequality. That’s why we alsogive the following direct conclusion:∥∥∥∥∥∂u∂n

∥∥∥∥∥∂T

≤ ‖∇u‖∂T ≤c√h‖∇u‖T (A.2.3b)

The constant c does not depend on the mesh size but it clearly depends on p, as the variation fromboundary surfaces to the inner volume may increase with rising polynomial degree. For clarity let’salso cite an expression for c√

hon a simplex from [WH03]. The equation (A.2.3a) holds true for

c√h

=√

(p+ 1)(p+ d)d

√√√√ |∂T ||T |

A.3 Orthogonal polynomials

Legrendre polynomials denoted by li are orthogonal polynomials spanning Pp([−1, 1]). They satisfy

(li, lj)L2([−1,1]) = 22i+ 1δi,j

with δi,j the kronecker symbol. They can be computed with the follow recurrence relations:

l0(x) = 1,l1(x) = x,

(n+ 1)ln+1(x) = (2n+ 1)ln(x)x− nln−1(x), n ≥ 1

Integrated Legrendre polynomials denoted by Li are defined as follows:

Ln(x) :=∫ x

−1ln−1(ξ) dξ for x ∈ [−1, 1] and n ≥ 2

They are mutually orthogonal with respect to the H1-seminorm, i.e.∫ 1

−1L′i(x)L′j(x) dx = 0 for i 6= j

The polynomials Ln, n ≥ 2 vanish at the interval bounds, i.e. Ln(−1) = Ln(1) = 0 and spanPp0 ([−1, 1]), that is the space of all polynomials up to degree p which have zero values on the bounds.For those polynomials there also holds a three-term-recurrence relation:

L1(x) = x,

L2(x) = 12(x2 − 1),

(n+ 1)Ln+1(x) = (2n− 1)Ln(x)x− (n− 2)Ln−1(x), n ≥ 2

A.3. ORTHOGONAL POLYNOMIALS 103

Note that L1 was added to allow for a general recursion relation although there does hold L1(x) 6=∫ x−1 l0(ξ) dξ An interesting property of the integrated Legendre polynomials is that they are almostorthogonal also with respect to the L2 inner product:

(Li, Lj)L2([−1,1]) = 0 for |i− j| > 2

Scaled Legrendre polynomials denoted by lSi are defined as follows:

lSn(x, t) := tnln(xt

) for x ∈ [−1, 1], t ∈ (0, 1]

Thanks to the multiplication with the monomial tn, the functions are free of fractions and stay poly-nomial of order n, i.e. lSn ∈ Pn([−1, 1] × [0, 1]) and the limit t → 0 is well defined. We again getthree-term recurrences:

lS0 (x, t) = 1,lS1 (x, t) = x,

(n+ 1)lSn−1(x, t) = (2n+ 1)lSn(x, t)x− nt2lSn−1(x, t), n ≥ 2

Scaled Integrated Legrendre polynomials denoted by LSi are defined as

LSn(x, t) := tnLn(xt

) for x ∈ [−1, 1], t ∈ (0, 1]

=∫ x

−tlSn(x, t) dξ for x ∈ [−1, 1], t ∈ (0, 1] and n ≥ 2

and they satisfy LSn ∈ Pn0 ([−1, 1]× [0, 1]) and the following three-term-recurrence relations:

L1(x)S = x,

L2(xS) = 12(x2 − t2),

(n+ 1)LSn+1(x) = (2n− 1)LSn(x)x− (n− 2)t2LSn−1(x), n ≥ 2

104 APPENDIX B. SELECTED PROOFS

Appendix B

Selected Proofs

Proof of Lemma 2.3.1

Proof. The proof is similar to the LBB-proof used in [HL02], but has some modifications to accountfor the additional facet functions.Let’s define w := ΠS

ku. Now start with the definition of the norm and work on each term one afteranother.

|||ΠSku|||1 = |||w|||1 =

∑T∈Th‖∇ w‖2

T + 1h

[[[wt]]]2∂T + h‖∂w∂n‖2∂T = ...

As w is piecewise polynomial we can use inverse estimates to absorb the last term by the first one

... ≤∑T∈Th

c1‖∇ w‖2T + 1

h[[[wt]]]2∂T = ...

For the first term we can use boundedness results of the H(div)-conforming interpolation. The secondpart is more tricky. To separate the element functions and the facet functions we shift the function uin by using the triangle inequality

... ≤∑T∈Th

c1‖∇ w‖2T + 1

h‖wt − ut‖2

∂T + 1h‖ut − wtF‖2

∂T = ...

For the last part we can use the L2-projection approximation property (k ≥ 1)

‖ut − wtF‖2∂T ≤ c2h

2‖u‖2∂T

and with the trace inequality [BS94, Theorem 1.6.6.] we additionally get

‖u‖2∂T ≤

c3

h‖u‖T‖u‖H1(T ) ≤

c3

h‖u‖2

H1(T )

This leaves us with (c4 = c2c3)

... ≤∑T∈Th

c1‖∇ w‖2T + c4‖u‖2

H1(T ) + 1h‖wt − ut‖2

∂T = ...

We now concentrate on the last problem, i.e. to bound ‖wt − ut‖2∂T . Let’s therefore substitute

v := wt − ut. Now we use the same trace inequality that we used before, but transform to thereference element T before and back afterwards:

‖v‖2∂T ≤ c5h

d−1‖v‖2∂T≤ c6

(hd2‖v‖T

) (hd−2

2 ‖v‖H1(T )

)≤ c9‖v‖T ( 1

h‖v‖T + |v|H1(T ))

105

due to

‖v‖2T≤ c7h

−d‖v‖2T and ‖v‖2

H1(T ) = ‖v‖2T

+ ‖∇ v‖2T≤ c7h

−d‖v‖2T + c8h

−(d−2)‖∇ v‖2T

Now with the interpolation results for the H(div)-conforming interpolation(see e.g. [BF91, III.3.6])we have

‖v‖T = ‖ΠΣk u− u‖T ≤ c10h‖u‖H1(T )

So finally we get (with c11 = c5c10 and C = c11 + c4 + c1)

... ≤ C∑T∈Th‖u‖2

H1(T ) ≤ C‖u‖2H1(Ω)x

Proof of Proposition 2.4.3

Proof. We won’t proof the Proposition down to all details. We will, at certain points, refer to [KJ98],where the proof for a similar situation has been performed in detail, instead.We start with representation (2.4.2b) to proof the result. Let’s divide volume and facet integrals intotwo separate expressions for which we will show boundedness:

Ch = (A) + (B) with

(A) =∑T∈Th

∫Tdiv(u⊗ w) v dx =

∑T∈Th

∫T

(w · ∇)u v dx (B.0.1a)

(B) =∑T∈Th

∫∂Tin|wn|(ut − utF ) vt ds+

∫∂Tout|wn|(utF − ut) vtF ds

(B.0.1b)

1. There holds(A) ≤ ‖w‖L4‖∇ u‖L2‖v‖L4 (B.0.2)

by applying Cauchy-Schwarz two times.

2. We now show that ‖w‖L4 ≤ C|||w|||1,∗∀ w ∈ [H1(Th)]d:

(a) We pose the problem

− div(∇φ) = w3, φ = 0 on ∂Ω where (w3)i = w3i (B.0.3)

and as w3 ∈ L4(Ω) we have φ ∈ W 2,p(Ω)∩H10 (Ω) ∀ p ≥ 1 if we assume smooth boundaries.

(b) The solution operator is continuous onto W 2, 43 (Ω) and so we get

‖φ‖2, 43≤ c1‖w3‖

L43

= c‖v‖3L4 (B.0.4)

(c) Now we pose the discrete problem to (B.0.3). Therefore we use a consistent and boundedbilinearform to discretize the vector-valued laplace operator. We use the results of [KJ98,Prop.4.5] where a Standard DG interior penalty bilinearform Xh is used. Due to consistencythere holds

‖w‖4L4 = (w3, w) = Xh(φ,w) (B.0.5)

106 APPENDIX B. SELECTED PROOFS

In the proof of [KJ98, Prop.4.5] it was shown that the bilinearform can be bounded in the‖ · ‖2, 4

3norm for φ and in a discrete H1-seminorm ‖ · ‖1,k similar to one used here for w:

Xh(φ,w) ≤ c2‖φ‖2, 43‖w‖1,k (B.0.6)

With [[w]]*2 ≤ [[w−]]2 + [[w+]]2 (where + and − indicate the different sides of a facet) thediscrete norm ‖ · ‖1,k is smaller than the one used here.

Xh(φ,w) ≤ c2‖φ‖2, 43|||w|||1,∗ (B.0.7)

(d) Putting together (B.0.5), (B.0.7) and (B.0.4) we get

‖w‖4L4 ≤ c1c2‖w‖3

L4|||w|||1,∗ ⇒ ‖w‖L4 ≤ c1c2|||w|||1,∗ (B.0.8)

3. So we conclude (A) ≤ |||w|||1,∗|||u|||1,∗|||v|||1,∗

4. We now turn over to the facet integrals.

(a) We get ( with ‖vtF‖ ≤ ‖v‖+ |[[vt]]|)

(B) ≤∑T∈Th

∫∂T‖w‖ ‖[[ut]]‖ (‖v‖+ |[[vt]]|) ds

≤∑T∈Th

h12 ‖w‖L4(∂T )

(‖v‖L4(∂T ) + ‖[[vt]]‖L4(∂T )

)|||u|||T1,∗

≤ |||u|||1,∗ ·( ∑T∈Th

h‖w‖2L4(∂T )

(‖v‖L4(∂T ) + ‖[[vt]]‖L4(∂T )

)2) 1

2

︸︷︷︸(C)

(B.0.9)

(b) And with a, b ≥ 0 and (a+ b)2 ≤ 2(a2 + b2) and√a+ b ≤

√a+√b we get

(C) ≤ 2(( ∑

T∈Thh‖w‖2

L4(∂T )‖v‖2L4(∂T )

) 12

︸︷︷︸(D)

+( ∑T∈Th

h‖w‖2L4(∂T )‖[[vt]]‖2

L4(∂T )

) 12

︸︷︷︸(E)

)(B.0.10)

(c) For (D) the estimates are clear if we use [KJ98, eq. 4.19] in (∗):

(D) ≤( ∑T∈Th

h ‖w‖4L4(∂T )

) 14( ∑T∈Th

h ‖v‖4L4(∂T )

) 14

(∗)≤ C(D)

(‖w‖4

L4(Ω) + ‖w‖41,k

) 14(‖v‖4

L4(Ω) + ‖v‖41,k

) 14

≤ C(D)|||w|||1,∗|||v|||1,∗

(B.0.11)

(d) For (E) the estimate is different and as the hybrid version is less common we have to putsome additional effort into this part. Remember that we just consider the two-dimensionalcase here! We use a scaled version of the Bernstein inequality (see [LS03, p.7]) in (∗)

107

which we are allowed to do as v is piecewise polynomial. We also use the rough estimateAT ≥ 0 :

√∑iAi ≤

∑i

√Ai:

(E) ≤∑T∈Th

h12‖w‖L4(∂T )‖[[vt]]‖L4(∂T )

(∗)≤

∑T∈Th

cph14‖w‖L4(∂T )‖[[vt]]‖L2(∂T )

≤ C

( ∑T∈Th

c2p h

32 ‖w‖2

L4(∂T )

) 12

|||v|||1,∗

≤ C

( ∑T∈Th

h ‖w‖4L4(∂T )

) 14( ∑T∈Th

h2 c4p︸︷︷︸

' |Ω|

) 14

|||v|||1,∗

≤ C(E)(|Ω|) |||w|||1,∗|||v|||1,∗

(B.0.12)

(e) with C(B) = 2C(D)C(E)(|Ω|) the estimate holds true for the facet integrals

5. Keep in mind ||| · |||1,∗ ≤ ||| · |||1

6. with β∗Ch = 1 + C(B) the proposition holds true.

108 BIBLIOGRAPHY

Bibliography

[ABCM02] Douglas N. Arnold, Franco Brezzi, Bernardo Cockburn, L.Donatella MariniUnified Analysis of Discontinuous Galerkin Methods for Elliptic ProblemsSIAM Journal for Numerical Analysis, Vol. 39, No. 5, pp. 1749-1779, 2002

[ARW95] Uri M. Ascher, Steven J. Ruuth, Brian T.R. WettonImplicit-explicit methods for time-dependent partial differential equationsSIAM Journal for Numerical Analysis, Vol. 32, No. 3, pp. 797-823, 1995

[ARS97] Uri M. Ascher, Steven J. Ruuth, Raymond J. SpiteriImplicit-explicit Runge-Kutta methods for time-dependent partial differential equationsApplied Numerical Mathematics, Vol. 25, pp. 151-167, 1997

[Brae97] Dietrich BraessFinite Elements: theory, fast solvers, and applications in solid mechanicsCambridge University Press, Cambridge, UK, 1997

[BF91] F. Brezzi, M. Fortin,Mixed and Hybrid Finite Element MethodsSpringer-Verlag, New York, 1991

[BS94] Susanne C. Brenner, L. Ridgway Scott,The Mathematical Theory of Finite Element MethodsSpringer-Verlag, Berlin, 1994

[CGL08] Bernardo Cockburn, Jayadeep Gopalakrishnan, Raytcho LazarovUnified Hybridization of Discontinuous Galerkin, Mixed and Continuous Galerkin Methods for SecondOrder Elliptic ProblemsPreprints (http://www.math.ufl.edu/ jayg/research/papers.html, November 2008)

[CKS07] Bernardo Cockburn, Guido Kanschat, Dominik SchötzauA Note on Discontinuous Galerkin Divergence-free Solutions of the Navier-Stokes EquationsJournal of Scientific Computing Vol. 31, p. 61 ff, May 2007; doi: 10.1007/s10715-006-9107-7

[CKS05] Bernardo Cockburn, Guido Kanschat, Dominik Schötzau,A locally conservative LDG method for incompressible Navier-Stokes equationsMathematics of Computations Vol. 74, pp. 1067-1095, 2005

[DB05] Leszek Demkowicz, Annalisa BuffaH1,H(curl) and H(div)-conforming projection-based interpolation in three dimensions, Quasi-optimal p-interpolation estimatesComput. Methods Appl. Mech. Eng. 194, pp. 267-296, 2005

BIBLIOGRAPHY 109

[DG10] Leszek Demkowicz, Jayadeep Gopalakrishnan,A class of Discontinuous Petrov-Galerkin Methods. Part I: The Transport EquationComputer Methods in Applied Mechanics and Engineering Volume 199, Issues 23-24, pp. 1558-1572,15 April 2010

[DGS08] Leszek Demkowicz, Jayadeep Gopalakrishnan, Joachim SchöberlPolynomial Extension Operators. Part ISIAM J Numerical Analysis, Vol46(6), pp. 3006-3031, 2008

[DH03] Jean Donea and Antonio HuertyFinite Element Methods for Flow ProblemsJohn Wiley & Sons Ltd, 2003

[ES09] Herbert Egger, Joachim SchöberlA Mixed-Hybrid-Discontinuous Galerkin Finite Element Method for Convection-Diffusion ProblemsIMA Journal of Numerical Analysis 2009; doi: 10.1093/imanum/drn083

[ESW05] Howard Elman, David Silvester, Andy WathenFinite Elements and Fast Iterative Solvers with applications in incompressible fluid dynamicsNumerical Mathematics and Scientific Computation, Oxford Science Publications, 2005

[Evan98] Lawrence C. EvansPartial Differential EquationsGraduate Studies in Mathematics, American Mathematical Society, Providence, Rhode Island, 1998

[EWY02] R. Ewing, J. Wang, and Y. YangA Stabilized Discontinuous Finite Element Method for Elliptic ProblemsNumerical Linear Algebra with Applications, 2002

[FHV97] J. Frank, W. Hundsdorfer, J.G. VerwerOn the stability of implicit-explicit linear multistep methodsApplied Numerical Mathematics, Vol. 25, pp. 193-205, 1997

[HL02] Peter Hansbo, Mats G. LarsonDiscontinuous Galerkin methods for incompressible and nearly incompressible elasticity by Nitsche’smethodComput. Methods Appl. Mech. Eng. 191, pp. 1895-1908, 2002

[KCGH07] Alex Kanevsky, Mark H. Carpenter, David Gottlieb, Jan S. HesthavenApplication of implicit-explicit high order Runge-Kutta methods to discontinuous-Galerkin schemesJournal of Computational Physics, Vol. 225, pp. 1753-1781, 2007

[KC03] C.A. Kennedy, Mark H. CarpenterAdditive Runge-Kutta schemes for convection-diffusion-reaction equationsApplied Numerical Mathematics, Vol. 44, pp. 139-181, 2003

[KJ98] Ohannes A. Karakashian, Wadi N. JureidiniA Nonconforming Finite Element Method for the Stationary Navier-Stokes EquationsSIAM Journal on Numerical Analysis, Vol. 35, pp. 93-120, 1998

110 BIBLIOGRAPHY

[Kova48] L.I.G. KovasznayLaminar flow behind a two-dimensional gridMathematical Proceedings of the Cambridge Philosophical Society, Vol. 44, pp. 58-62, 1948

[LS03] Andris Lasis, Endre SüliPoincaré-type inequalities for broken Sobolev SpacesTechnical Report 03/10, Oxford University Computing Laboratory, 2003

[Monk03] Peter MonkFinite Element Methods for Maxwell’s EquationsNumerical Mathematics and Scientific Computation. The Clarendon Press Oxford University Press,New York, 2003

[NPC09] N.C. Nguyen, J.Peraire, B.CockburnAn implicit high-order hybridizable discontinuous Galerkin method for linear convection-diffusionequationsJournal of Computational Physics, Vol. 228, pp. 3232-3254, 2009

[Schö08] Joachim SchöberlDivergence-free Hybrid Discontinuous Galerkin Finite Elements for Incompressible Navier StokesEquations, Augmented Lagrangians, and Robust PreconditionersTalk, Durham Symposium on CLAPDE, July 2008

[SZ05] Joachim Schöberl, Sabine ZaglmayrHigh order Nedelec elements with local complete sequence propertiesInternational Journal for Computation and Mathematics in Electrical and Electronic Engineering,Vol. 24(2), pp. 374-384, 2005

[SZ09] Dominik Schötzau, Liang ZhuA robust a-posteriori error estimator for discontinuous Galerkin methods for convection-diffusionequationsApplied Numerical Mathematics, Vol. 59, pp. 2236-2255, 2009

[WH03] T. Warburton, J.S. HesthavenOn the constants in hp-finite element trace inverse inequalitiesComput. Methods Appl. Mech. Eng., Vol. 192, pp. 2765-2773, 2003

[Zagl06] Sabine ZaglmayrHigh Order Finite Element Methods for Electromagnetic Field ComputationDissertation, Institut für Numerische Mathematik, July 2006

HDG for incompressible flow - RWTH Aachen …...Chapter 2translates the discretization discussed in...

Documents

Transcript of HDG for incompressible flow - RWTH Aachen …...Chapter 2translates the discretization discussed in...