UHasselt - A Hybridized Discontinuous Galerkin Method for Time … · 2017. 2. 3. · Dezember 2013...

A Hybridized Discontinuous GalerkinMethod for Time-Dependent

Compressible Flows

Master’s Thesis

by

Alexander Jaust

Institut für Geometrie und Praktische Mathematik

Rheinisch-Westfälische Technische Hochschule Aachen

December 9th, 2013

Supervisor:Prof. Dr. Sebastian Noelle

Co-Examiner:Dr. Jochen Schütz

Eidesstattliche Versicherung

Hiermit erkläre ich, dass ich die vorliegende Masterarbeit selbstständig verfasst habe. Eswurden keine anderen als die in der Arbeit angegebenen Quellen und Hilfsmittel benutzt.Die wörtlichen oder sinngemäß übernommenen Zitate habe ich als solche kenntlich gemacht.

Aachen, den 9. Dezember 2013 Alexander Jaust

Acknowledgments

I want to thank my supervisor Prof. Dr. Sebastian Noelle for giving me the opportunity to work atthis thesis in an interesting field of computational fluid dynamics. I thank my advisor Dr. JochenSchütz who introduced me to hybridized discontinuous Galerkin methods and supported me duringmy thesis. Thanks to Michael Woopen who was always giving helpful advises about the code.Further thanks go to Tobias Leicht from the DLR for providing me the opportunity to participateat the 2nd International Workshop on High-Order CFD Methods in Cologne.

I also want to thank my family for their support. Without them I could not have pursued mystudies.

Contents

1. Introduction 1

2. Governing Equations 32.1. Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2. Scalar Convection-Diffusion Equation . . . . . . . . . . . . . . . . . . . . . . . . . 42.3. Euler Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.4. Navier-Stokes Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.5. Dimensionless Numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3. Spatial Discretization 83.1. The Hybridized Discontinuous Galerkin Method . . . . . . . . . . . . . . . . . . . 8

3.1.1. Hybridization of the Nonlinear Convection-Diffusion Equation . . . . . . . . 93.1.2. Discretization of the Nonlinear Convection-Diffusion Equation . . . . . . . 10

4. Time Integration Methods 134.1. General Setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

4.1.1. Consistency and Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . 144.1.2. Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154.1.3. A- and A(α)-Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174.1.4. L-Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

4.2. Diagonally Implicit Runge-Kutta Methods . . . . . . . . . . . . . . . . . . . . . . . 184.2.1. Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214.2.2. Time Step Adaption . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

4.3. Backward Differentiation Formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . 254.3.1. Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

5. Time Integration for the Hybridized Discontinuous Galerkin Method 315.1. Notes on Differential-Algebraic Equations . . . . . . . . . . . . . . . . . . . . . . . 315.2. Backward Differentiation Formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . 335.3. Singly Diagonally Implicit Runge-Kutta Methods . . . . . . . . . . . . . . . . . . . 33

6. Shock-Capturing 356.1. An Artificial Viscosity Model for the Hybridized Discontinuous Galerkin Method . 35

7. Numerical Results 387.1. Rotating Gaussian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

7.1.1. Convergence Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 397.1.2. Time Step Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

7.2. Radial Expansion Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 437.3. Kármán Vortex Street . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 477.4. Sod’s Shock Tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

7.4.1. Constant Viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 517.4.2. Adaptive Artificial Viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

8. Conclusions and Outlook 58

A. Additional and Tabulated Results 60A.1. Rotating Gaussian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60A.2. Radial Expansion Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63A.3. Sod’s Shock Tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

Contents ii

Bibliography 68

List of Figures

3.1. Numerical method: Domain Ω split into two parts by an edge Γ. . . . . . . . . . . 9

4.1. Stability region: Examples of the stability regions of an A-stable method (left) andan A(α)-stable method with α = π4 (right) are shown. . . . . . . . . . . . . . . . . 18

4.2. Stability region: The stability regions of all SDIRK methods. Stable regions withR(z) ≤ 1 are indicated white while unstable regions are red. . . . . . . . . . . . . . 22

4.3. Initialization of the time integration for the BDF3 method. . . . . . . . . . . . . . 264.4. Stability region: The stability regions of BDFk methods for k = 1, . . . , 4. Stable

regions with R(z) ≤ 1 are indicated white while unstable regions are red. . . . . . . 284.5. Solution strategy with error estimator. . . . . . . . . . . . . . . . . . . . . . . . . . 30

6.1. Artificial Viscosity: Approximation of the ramp function (left) and sonic speed (right). 36

7.1. Rotating Gaussian: The initial distribution (left) and final distribution at t = π4(right) of the Gaussian. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

7.2. Rotating Gaussian: Convergence of ‖w−wh‖L2 under uniform spatial and temporalrefinement with p = 2 (left) and p = 3 (right). . . . . . . . . . . . . . . . . . . . . . 39

7.3. Rotating Gaussian: Error in dependence of time step size ∆t (left) and tolerance tol(right) on a fixed mesh with Ne = 512 elements and cubic ansatz functions. . . . . 40

7.4. Rotating Gaussian: time step evolution for different tolerances tol on a fixed meshwith Ne = 512 elements and cubic ansatz functions. . . . . . . . . . . . . . . . . . 41

7.5. Rotating Gaussian: Error evolution on different grids with cubic ansatz functions(left). The tolerance is divided by a factor of eight or sixteen depending on the orderof the method (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

7.6. Radial Expansion Wave: Density ρ on the domain at t = t0 (left) and t = T (right). 437.7. Radial Expansion Wave: Quadratic ansatz functions and Ne = 2048 elements . . . 457.8. Radial Expansion Wave: Cubic ansatz functions and Ne = 8192 elements . . . . . 467.9. Kármán vortex street: Section of the mesh. . . . . . . . . . . . . . . . . . . . . . . 477.10. Kármán vortex street: Mach number distribution at two instances approximately a

half period apart showing the periodic behavior of the flow. . . . . . . . . . . . . . 477.11. Von Kármán vortex street: Results obtained for tol = 0.05 and 10−2 ≤ ∆t ≤ 8. . . 507.12. Sod’s shock tube: Initial data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 517.13. Sod’s shock tube: Exact solution at T = 0.2. . . . . . . . . . . . . . . . . . . . . . 527.14. Sod’s shock tube: Sketch of the mesh. . . . . . . . . . . . . . . . . . . . . . . . . . 527.15. Sod’s shock tube: Errors in the density and velocity for constant viscosity ε̃. . . . . 537.16. Sod’s shock tube: Exact and approximated solution of density and velocity for

different bulk viscosities ε̃0 at the discontinuity and shock. A mesh with Ne = 400elements and polynomials of degree p = 2 have been used. . . . . . . . . . . . . . . 54

7.17. Sod’s shock tube: Error in density ρ and velocity u for different polynomial degreesp and values for bulk viscosity ε̃0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

7.18. Sod’s shock tube: Exact and approximated solution of density and velocity fordifferent bulk viscosities ε̃0 at the discontinuity and shock. A mesh with Ne = 400elements and polynomials of degree p = 2 have been used. . . . . . . . . . . . . . . 57

A.1. Radial Expansion Wave: Quadratic ansatz functions and Ne = 8192 elements. . . . 64A.2. Radial Expansion Wave: Cubic ansatz functions and Ne = 2048 elements. . . . . . 65

List of Tables

4.1. Coefficients of Cash’s SDIRK method. . . . . . . . . . . . . . . . . . . . . . . . . . 204.2. Coefficients of Al-Rabeh’s SDIRK method. . . . . . . . . . . . . . . . . . . . . . . 214.3. Coefficients of Hairer’s and Wanner’s SDIRK method. . . . . . . . . . . . . . . . . 214.4. Range of values for γ to retain an A-stable method. . . . . . . . . . . . . . . . . . 234.5. Range of values for γ to retain an L-stable method. . . . . . . . . . . . . . . . . . . 234.6. Coefficients of BDFk-methods up to order 4. . . . . . . . . . . . . . . . . . . . . . 254.7. Angle α of A(α)-stable BDFk-methods up to order 4. . . . . . . . . . . . . . . . . 28

7.1. Kármán vortex street: Mean drag coefficients cD and Strouhal numbers Sr fromexperiments and previous publications. . . . . . . . . . . . . . . . . . . . . . . . . . 48

7.2. Kármán vortex street: Mean drag coefficients cD and Strouhal numbers Sr determinedby Cash’s SDIRK method using different tolerances and bounds for the time step. 48

7.3. Kármán vortex street: Mean drag coefficients cD and Strouhal numbers Sr determinedwith Al-Rabeh’s SDIRK method using different tolerances and bounds for the timestep. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

7.4. Kármán vortex street: Mean drag coefficients cD and Strouhal numbers Sr determinedwith Hairer’s and Wanner’s SDIRK method using different tolerances and boundsfor the time step. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

A.1. Rotating Gaussian: Convergence of BDF methods on fixed mesh Ne = 512, cubicansatz functions and fixed time step. . . . . . . . . . . . . . . . . . . . . . . . . . . 60

A.2. Rotating Gaussian: Convergence of SDIRK methods on fixed mesh Ne = 512, cubicansatz functions and fixed time step. . . . . . . . . . . . . . . . . . . . . . . . . . . 60

A.3. Rotating Gaussian: Convergence of BDF2 for quadratic and cubic ansatz functions. 61A.4. Rotating Gaussian: Convergence of BDF3 method for quadratic and cubic ansatz

functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61A.5. Rotating Gaussian: Convergence of Cash’s method for quadratic and cubic ansatz

functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61A.6. Rotating Gaussian: Convergence of Al-Rabeh’s method for quadratic and cubic

ansatz functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61A.7. Rotating Gaussian: Convergence of Hairer’s and Wanner’s method for quadratic and

cubic ansatz functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61A.8. Rotating Gaussian: Convergence of SDIRK methods on fixed mesh Ne = 512, cubic

ansatz functions and variable time step. . . . . . . . . . . . . . . . . . . . . . . . . 62A.9. Rotating Gaussian: Convergence of SDIRK methods when mesh and tolerance are

refined using cubic ansatz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62A.10.Radial Expansion Wave: Entropy error and total number of Newton steps on the

coarse mesh with quadratic ansatz functions. . . . . . . . . . . . . . . . . . . . . . 63A.11.Radial Expansion Wave: Entropy error and total number of Newton steps on the

fine mesh with quadratic ansatz functions. . . . . . . . . . . . . . . . . . . . . . . . 63A.12.Radial Expansion Wave: Entropy error and total number of Newton steps on the

coarse mesh with cubic ansatz functions. . . . . . . . . . . . . . . . . . . . . . . . . 63A.13.Radial Expansion Wave: Entropy error and total number of Newton steps on the

fine mesh with cubic ansatz functions. . . . . . . . . . . . . . . . . . . . . . . . . . 64A.14.Sod’s Shock Tube: Density errors ‖ρ(T )− ρh(T )‖L2 for constant viscosity. . . . . . 66A.15.Sod’s Shock Tube: Velocity errors ‖u(T )− uh(T )‖L2 for constant viscosity. . . . . 66A.16.Sod’s Shock Tube: Density errors ‖ρ(T )− ρh(T )‖L2 for adaptive viscosity. . . . . . 67A.17.Sod’s Shock Tube: Velocity errors ‖u(T )− uh(T )‖L2 for adaptive viscosity. . . . . 67

1. Introduction

For a long time experiments have been one of the most important tools in science and industry.However, experiments are usually expensive, time consuming and sometimes hardly executable.With the occurrence of computers and their performance increasing, simulations become a moreand more interesting alternative to experiments. They offer a cheaper environment where certainexperiments can be modeled. It is possible to easily change key quantities to study their influenceon the experiment.

In this work we focus on computational fluid dynamics (CFD), i.e., the modeling of fluid flowssuch as the flow of air around an airfoil, for example. Such flows are usually described by partialdifferential equations (PDEs). These are differential equations of multivariate functions in which theirpartial derivatives occur. In many cases it is only possible to solve these equations approximatelyusing simulations. A common simulation technique for PDEs arising in CFD are finite-volumemethods (FV methods) [41–43, 70–73]. Current codes employing finite-volume schemes are usuallysecond order accurate. With order of accuracy q we denote that the error e is proportional tohq with h being the mesh size. One often speaks of high-order methods when the order of themethod is higher than two. By now, finite-volume methods are well-accepted in engineering. Theyshow stable behavior in particular for steady-state problems due to strong numerical dissipation.However, for unsteady problems that need to preserve complex structures of a flow over a longertime, e.g. vortices, the dissipation becomes disadvantageous because the structures are smoothenedout. Moreover, for increased accuracy the mesh size needs to be decreased much more than forhigher order methods. In general, this leads to a much higher number of unknowns such that alsothe time needed for a simulation is larger.

A promising method for high-order CFD are the discontinuous Galerkin methods (DG) [5, 11,19, 28, 29, 34, 35, 37, 38, 53, 63]. The computational domain is triangulated and then the solutionis represented by a polynomial of degree p on each of the elements. A drawback is the highernumber of degrees of freedom compared to continuous Galerkin (CG) methods [37]. Similar toDG methods the solution is approximated using polynomials, but the approximated solution isassumed to be continuous over the whole triangulation. For the DG method continuity is onlyenforced on each element of the triangulation such that jumps in the solution between elementsare possible. Problems in CFD are usually described by hyperbolic or nearly hyperbolic partialdifferential equations for which classical CG methods are usually not stable without stabilization[14, 39]. However, DG methods are stable for such equations and therefore of greater interest inCFD.

The degrees of freedom that are coupled globally for a DG, but also other methods, can bereduced using hybridization [26]. A new characterization for elliptic problems has been given byCockburn and Gopalakrishnan [20] and later unified to be applied to different numerical methodsby Cockburn et al. [21]. Based on this the hybridized discontinuous Galerkin method (HDG) hasbeen developed for different partial differential equations including the equations of fluid dynamics,especially the compressible Euler and Navier-Stokes equations [45–47, 49, 59, 61]. The methodintroduces additional unknowns on the edges between elements increasing the total number ofdegrees of freedom. However, this allows the reformulation of the problem such that the numberof globally coupled unknowns is much lower than for the initial problem at the cost of a morecomplicated assembling of the system of equations. Since the computational costs and the memoryrequirements are mainly dependent on the system size of globally coupled unknowns the HDGmethod still reduces the overall computational work.

We focus on the HDG method for time dependent problems in this work. An implementationbased on Netgen, NGSolve [57] and PETSc [7–9] has been done in an earlier work [59, 77]. It hasbeen extended by time integration methods based on backward differentiation formulas (BDF) bySchütz et al. [62]. Other experiments with BDF and diagonally implicit Runge-Kutta (DIRK)

1. Introduction 2

methods have been executed by Nguyen and Peraire [45, 47].We examine time integration based on singly diagonally implicit Runge-Kutta (SDIRK) methods

with an embedded error estimator [1, 15, 18, 33]. These methods are quite popular for solving stiffordinary differential equations (ODEs) that arise from many different real world applications aschemical reactions, for example. With the error estimator it is possible to adapt the time stepduring the simulation in order to control the temporal error. It is investigated how well the errorestimator works to obtain simulation results with a desired accuracy. Furthermore, these methodsare compared to the BDF methods for which it is much more complicated to have a time stepadaptation because of the nature of the method and the absence of an integrated error estimator.

In order to apply the HDG method to a wide range range of problems the solver has to be able tohandle discontinuities. These often occur in supersonic fluid flow in the shape of shocks or contactdiscontinuities. One way to handle these discontinuities is to use of an artificial viscosity model. Inthe work we study an artificial viscosity model that has been introduced by Nguyen and Peraire[44] for HDG methods. It has been originally used for steady state problems only. We investigatewhether it can be applied to time-dependent problems as well.

This work is structured in the following way. First, the governing equations used in this work areintroduced, namely the convection-diffusion equation, the Euler equations and the Navier-Stokesequations. In the following chapter we introduce the spatial discretization using the hybridizeddiscontinuous Galerkin method by applying it to the convection-diffusion equation. Then, inChapter 4, different time integration methods for ordinary differential equations are presented.Especially their stability properties and the possibility of time step adaptation for these methods areexamined. The application of these methods to the hybridized DG method is shown afterwards inChapter 5. Following, we describe the artificial viscosity model in Chapter 6. The implementationis validated in Chapter 7. Moreover, some results for different test cases are discussed to quantifyhow good the time integration methods work. Finally, a short conclusion and an outlook of futurework is given in Chapter 8.

2. Governing Equations

The problems discussed in this work are mainly motivated by computational fluid dynamics (CFD).This means in most cases the flow field of a certain fluid, mostly air, has to be determined.Additionally to the flow field itself and occuring phenomena as vortices or shocks e.g., the transportof some scalars within the flow field may be of interest.

For a simulation a model of the problem in form of a set of equations must exist. These equationsare presented in this section and are all partial differential equations (PDE). First we introducethe general notation. Then we present the (nonlinear) convection-diffusion equation. It is oftenused for developing numerical methods in computational fluid dynamics as it describes the twomain types of transport that occur, namely convection and diffusion. Real-world flows are oftenmodeled using the Euler equations or the Navier-Stokes equations that are introduced afterwards.The Euler equations describe inviscid flows which can be seen as a special case of the Navier-Stokesequations under the assumption of vanishing viscosity.

2.1. Notation

To apply the discretization method to a wide range of problems in a satisfactory manner, it isbeneficial to introduce generalized notation for problems. For the hybridized DG method we writeproblems in the following form:

wt +∇ · (f(w)− fv(w,∇w)) = h(w,∇w). (2.1)

Here, w is unknown and may be vector valued. With wt we denote the time derivative wt :=∂w∂t .

Moreover, there are spatial derivatives, denoted by the gradient operator ∇ and divergence operator∇·, of the convective flux f(w) and the viscous flux fv(w,∇w). On the right hand side h(w,∇w) isa source or sink term. Depending on the problem certain terms may vanish.

In this work we focus on two dimensional (d = 2) problems. Therefore two spatial directions x1and x2 exist. This is no restriction of the generality of the method. It can also be applied to one(d = 1) or three dimensional (d = 3) problems.

For the sake of completeness we define the used operators for spatial derivatives of vectors andmatrices. We introduce these notations under the assumption of exactly two spatial dimensions(d = 2). The gradient operator ∇ applied to some vector v ∈ Rm is defined as

∇v =

∂v1∂x1

∂v1∂x2

......

∂vm∂x1

∂vm∂x2

(2.2)and is a matrix of size m× 2. The divergence of a vector v ∈ R2 and of a matrix A ∈ Rm×2 aredefined as

div v = ∇ · v =2∑i=1

∂vi∂xi

=∂v1∂x1

+∂v2∂x2

(2.3)

and

∇ ·A := ∇ ·

a11 a12... ...am1 am2

= ∇ · (a11, a12)

T

...

∇ · (am1, am2)T

=

∂a11∂x1

+ ∂a12∂x2...

∂am1∂x1

+ ∂am2∂x2

. (2.4)

2. Governing Equations 4

2.2. Scalar Convection-Diffusion Equation

Convective transport and diffusion processes can be described using the convection-diffusion equation.In general, the equation can be written as

wt +∇ · (f(w)− fv(w,∇w)) = 0, (2.5)

In our case w is an unknown scalar. Nevertheless, the convection-diffusion equation can also beformulated in case w is a vector. The fluxes are given by

f(w) = uw, fv(w,∇w) = ε∇w. (2.6)

The term convective flux f(w) describes the convective transport of quantity w by a flow field withvelocity u = (u1, u2)

T . In case w contains components of the velocity the problem becomes nonlinear.The viscous flux fv(w,∇w) accounts for diffusive transport with ε > 0 being the diffusivity. Werequire positive values for ε because otherwise the problem is ill-conditioned. The right hand sideis zero, because no sources and sinks h(w,∇w) are present. All terms may be functions of space(x1, x2) and time t.

As already mentioned in the introduction to this chapter, equation (2.5) is often used as a modelequation when developing numerical methods. The reason for this is that flow problems are describedby coupled convective and diffusive transport processes. While yielding a simpler equation thanthe full equations of fluid dynamics the main features are still covered by the convection-diffusionequation. It also covers different types of partial differential equations. For u = (0, 0)T and wt = 0one has a so-called elliptic problem. A common property is that this equation leads to smoothsolutions even if the initial conditions are non-smooth. On the other side one has a hyperbolicproblem for vanishing diffusivity ε = 0. Then in the nonlinear case non-smooth solutions maydevelop even if the initial data is smooth. In the setting of fluid dynamics the discontinuities areoften referred to as shocks.

Since the equation is a partial differential equation of mixed type many different effects mayoccur depending on the initial conditions and diffusivity ε. This makes it challenging to developnumerical methods because the discretizetion method must be able to handle a large number ofeffects that may occur. On the other hand, this makes the convection-diffusion a very powerful testequation.

2.3. Euler Equations

A very important set of equations in fluid dynamics are the Euler equations. They describe thespecial case of an inviscid flow and are often used in gas dynamics to model high speed flows of air.Using our notation the equations can be written as

wt +∇ · f(w) = 0. (2.7)

Neither viscous fluxes fv(w,∇w) nor source terms h(w,∇w) appear in the Euler equations. Theconvective flux is split into two parts f(w) := (f1(w), f2(w)) with

f1(w) =(ρu1, P + ρu

21, ρu1u2, u1(E + P )

)T(2.8)

f2(w) =(ρu2, ρu1u2, P + ρu

22, u2(E + P )

)T. (2.9)

The vector of unknownsw = (ρ, ρu1, ρu2, E)

T(2.10)

contains the density ρ, momentum ρu1 and ρu2 in the two spatial directions and the energy E.Therefore, this is a system of four equations. The first equation is called continuity equation and itensures mass conservation. The second and third equation describe the momentum conservation ineach spatial direction. There occurs a term containing the gradient of the pressure P . It describesthe change of momentum due to forces acting on the surface of the fluid. The last equation describesthe energy conservation with E being the energy. In general, the gravitational acceleration g may


have an influence on the fluid, but as we focus on flows of gases in this work, the density has rathersmall values. Therefore, the influence of the gravitational acceleration on the fluid is neglected suchthat the terms containing g are assumed to be zero.

In two dimensions, there are five unknowns (density ρ, velocities u1 and u2, energy E and thepressure P ), but only four equations, leading to an under-determined system of equation. To resolvethis issue a closure is needed. A common choice is the equation of state

P = (γ − 1)(E − 1

2ρ‖u‖22

)(2.11)

where γ is the ratio of specific heats. It depends on the fluid and is approximately γ = 1.4 for air.The Euler equations are a nonlinear, hyperbolic system of PDEs. As already mentioned, problems

of this type may lead to non-smooth solutions containing shocks. This coincides with observationsof so-called transonic and supersonic flows in gas dynamics where shocks appear.

2.4. Navier-Stokes Equations

A wider set of flows than described by the Euler equations are covered by the Navier-Stokesequations. In contrast to the Euler equations the fluid is assumed to be viscous. Therefore, also aviscous flux appears in the equations

wt +∇ · (f(w)− fv(w,∇w)) = 0. (2.12)

The convective flux f(w) = (f1(w), f2(w)) equals the convective flux of the Euler equations (cf. eq.(2.8) and (2.9)). The viscous flux fv(w,∇w) is split in two parts similar to the convective flux withcomponents

fv,1(w,∇w) =(

0, τ11, τ21, τ11u1 + τ12u2 + κ∂θ

∂x1

)T, (2.13)

fv,2(w,∇w) =(

0, τ21, τ22, τ21u1 + τ22u2 + κ∂θ

∂x2

)T. (2.14)

The Navier-Stokes equations describe conservation of mass, momentum and energy as the Eulerequations. Therefore, the vector of unknowns w is the same as for the Euler equations (cf. eq.(2.10)). However, additional terms occur for the momentum and energy equations because of viscouseffects.

The entries τij in the viscous flux refer to the entries of the stress tensor S. For a Newtonianfluid it can be written as

S =

(τ11 τ12τ21 τ22

)= µ

(∇u+ (∇u)T − 2

3(∇ · u)I

). (2.15)

with the dynamic viscosity µ and the identity matrix I. It causes change of momentum and energydue to friction. The second additional term is the heat flux q which is approximated using Fourier’slaw [30]

q = κ∇θ (2.16)with the temperature θ and heat conductivity κ. In the case of an inviscid flow µ = 0 also the heatflux q = 0 vanishes and one obtains the Euler equations again. The viscosity is approximated usingSutherland’s law [68]

µ = µ(θ) =C1θ

32

θ + C2(2.17)

with C1 and C2 being constants depending on the fluid. Yet another closure is needed for thetemperature θ. It can be computed using the ideal gas law

θ =µγ

κ · Pr

(E

ρ− 1

2

(u21 + u

22

))=

1

(γ − 1)cvP

ρ. (2.18)


Here, the Prandtl number Pr =µcpκ and the specific heats of constant pressure cp and of constant

volume cv occur.Compared to the Euler equations the Navier-Stokes equations are more complex especially due

to the stress tensor. This allows the description of phenomena as turbulence [51] and boundarylayers which can not be described by the Euler equations because these are viscous effects.

2.5. Dimensionless Numbers

An important tool in fluid dynamics are dimensionless numbers. They are used in similitude andcan be determined by doing a dimensional analysis of the equations describing a problem. Webriefly introduce the dimensionless numbers we are using in this work.

Reynolds Number

The Reynolds number Re quantifies the ratio of inertial to viscous forces. It is defined as

Re =ρ0UL

µ0=UL

ν0. (2.19)

For the density ρ0 and viscosity µ0 a specific value at a reference point is taken. This could be inthe free stream, for example. The viscosity ν is the kinematic viscosity ν = µρ that allows a shorterrepresentation since the density cancels out. The velocity U is some representative velocity of thefluid and L is a characteristic length influencing the fluid. This may be the height or diameter of achannel or the the diameter of an obstacle. The Reynolds number is useful when determining ifa flow contains instationary effects as vortex shedding or turbulence. However, it is usually notpossible to specify sharp values of Re for which these effects occur because they also depend onquantities not covered by the Reynolds number [51, 58].

Prandtl Number

The Prandtl number Pr describes the ratio of kinematic viscosity to the thermal diffusivity:

Pr =ν

a=µcpκ. (2.20)

It only depends on the fluid properties, such as the thermal diffusivity a, but not on the length scaleas the Reynolds number does. We assume the Prandtl number to be constant and set Pr = 0.72which is a common choice for air at moderate temperatures. It is needed to compute the temperatureθ from the ideal gas law (cf. equation (2.18)).

Mach Number

For flows described by the Euler or Navier-Stokes equations the Mach number Ma is important tomeasure the influence of compressible effects. It is defined as the ratio of the characteristic velocityU and the speed of sound c:

Ma =U

c. (2.21)

In case of an ideal gas the speed of sound can be determined from c =√γ Pρ . The speed of sound

actually defines the speed with which information is transported within the fluid. For Ma < 0.3 thedensity variations are less than 5% such that simplifications as incompressibility, i.e., ρ = const,are often applied. For Mach numbers greater than one, i.e., Ma > 1, information can only betransported downstream.


Strouhal Number

In the special case of a flow around an obstacle, vortices shed for certain Reynolds numbers. Usingthe frequency f of the vortex shedding the Strouhal number Sr is defined as:

Sr =fL

U. (2.22)

The characteristic length L is the diameter of the obstacle and U is the flow velocity.

3. Spatial Discretization

In this section the hybridized DG method (HDG) is introduced. It is explained using the steady-stateconvection-diffusion equation. First we present some general notes on the hybridized DG method.Then we illustrate how to hybridize the convection-diffusion. Finally, the actual discretization ofthe convection-diffusion equations is described.

For this, some function spaces are needed. The space of square integrable functions on a domainΩ is

L2(Ω) :=

{f :

∫Ω

|f |2dx =: ‖f‖L2(Ω)

3. Spatial Discretization 9

Ω1

Ω2

Γ

n2

n1

Figure 3.1.: Numerical method: Domain Ω split into two parts by an edge Γ.

dimensions. Therefore, the hybridized DG method becomes more beneficial as the polynomialdegree is increased [59].

All in all, the hybridized DG method is a promising high-order method. Besides the reducednumber of globally coupled unknowns, it allows many local computations which can be exploited byparallel computations. Furthermore, spatial adaptation is not only possible by refining the mesh,but also by choosing ansatz functions of different polynomial degree within a mesh.

3.1.1. Hybridization of the Nonlinear Convection-Diffusion Equation

In Section 2.2 the general case of the convection-diffusion equation was introduced. Here, wefocus on the problem for a given open bounded domain Ω with homogeneous Dirichlet boundaryconditions. Other boundary conditions can also be applied as von Neumann boundary conditions,for example. A more detailed discussion on the application of boundary conditions to the hybridizedmethods is given in [45]. Here, we investigate the stationary case, i.e., wt = 0; the time-dependentcase is treated in Section 5. The nonlinear convection-diffusion equation with boundary conditionsis given by:

∇ · f(w)− ε∆w = h, ∀x ∈ Ωw = 0, ∀x ∈ ∂Ω. (3.4)

We assume that the solution w has certain regularity such that it is at least in H2(Ω). Then, theweak formulation of the problem is:

N(w,ϕ) := −∫Ω

f(w) · ∇ϕdx+ ε∫Ω

∇w · ∇ϕdx !=∫Ω

ϕhdx, ∀ϕ ∈ H10 (Ω). (3.5)

with a test function ϕ. Further, we assume that the domain Ω can be split into two parts Ω = Ω1∪Ω2by an edge Γ := Ω1 ∩ Ω2 with normal vectors n1,2 as in Figure 3.1. Then, the initial problem (3.4)can be formulated separately on each domain Ω1 and Ω2:

∇ · f(w1)− ε∆w1 = h ∀x ∈ Ω1, w1 = 0 ∀x ∈ ∂Ω ∩ ∂Ω1, w1 = λ ∀x ∈ Γ∇ · f(w2)− ε∆w2 = h ∀x ∈ Ω2, w2 = 0 ∀x ∈ ∂Ω ∩ ∂Ω2, w2 = λ ∀x ∈ Γ .

(3.6)

It has the same structure as the convection-diffusion equation on the initial domain Ω, but with anadditional boundary condition on the edge Γ. This condition λ is introduced to retain continuity ofthe solution over the edge. Moreover, the solutions w1 and w2 of Problem (3.6) inherit the regularityassumption of the global solution w such that they have to be at least in H2(Ωi), i = {1, 2}. Now,we introduce a new global solution ŵ on Ω that is constructed by the partial solutions w1 and w2:

ŵ =

{w1, ∀x ∈ Ω1w2, ∀x ∈ Ω2

, wi ∈ H2(Ωi), i = {1, 2} . (3.7)

The initial problem (3.4) and the split problem (3.6) must have the same solution ŵ!= w since the

problems shall be equivalent. To ensure that ŵ coincides with w the variable λ has to be chosensuch that continuity is retained and that the initial problem is solved. This means we require

N(w,ϕ)!= N(ŵ, ϕ) =

∫Ω

ϕhdx. (3.8)


From this, an additional equation, that λ has to fulfill, can be derived by partially integratingN(ŵ, ϕ) such that

N(ŵ, ϕ)=−∫

Ω

f(ŵ) · ∇ϕdx+ ε∫

Ω

∇ŵ · ∇ϕdx (3.9)

=

∫Ω1

∇ · f(w1)ϕdx+∫

Ω2

∇ · f(w2)ϕdx+∫

Γ

ϕε(∇ŵ1 −∇ŵ2) · n2dσ (3.10)

−∫

Γ

ϕ(f(ŵ1) · n2 + f(ŵ2) · n1)dσ − ε∫

Ω1

∆w1ϕdx− ε∫

Ω2

∆w2ϕdx

=

∫Ω

ϕhdx−∫

Γ

ϕ(f(ŵ1) · n2 + f(ŵ2) · n1)dσ + ε∫

Γ

ϕ(∇w1 −∇w2) · n2dσ (3.11)

with ŵi = limτ→0+

w(x+ τni). In comparison to the weak formulation of the initial problem (3.5) two

additional integrals over the edge Γ occur. To ensure the equivalence of the problems the additionalterms must vanish such that equation (3.8) holds true. Therefore we require:∫

Γ

ϕ (ε(∇ŵ1 −∇ŵ2) · n2 − f(ŵ1) · n2 − f(ŵ2) · n1) dσ = 0. (3.12)

In our case ŵ1,2 are actually ŵ2 = w1 and ŵ1 = w2 on the edge Γ. The variable λ is used to enforcethat the integral (3.12) vanishes. With this additional equation the split problem (3.6) can beformulated as

∇ · f(w1)− ε∆w1 = h ∀x ∈ Ω1, w1 = 0 ∀x ∈ ∂Ω ∩ ∂Ω1, w1 = λ ∀x ∈ Γ (3.13)∇ · f(w2)− ε∆w2 = h ∀x ∈ Ω2, w2 = 0 ∀x ∈ ∂Ω ∩ ∂Ω2, w2 = λ ∀x ∈ Γ (3.14)

Jε∇w − f(w)K = 0,∀x ∈ Γ (3.15)

such that it is finally equivalent to the original problem (3.4). Here J·K denotes the jump-operatorthat is defined as

JwK := w1n1 − w2n1. (3.16)It can be seen that it is possible to divide the initial problem into two problems. These subproblemsonly depend on given boundary conditions and the unknown λ on the edge Γ. From this it followsthat the problems can be solved independently once λ is known. The idea of HDG methods isto reformulate the problem in a way such that a problem arises which depends only on λ. Thisrepresents all globally coupled unknowns and then allows to solve locally for the solution withineach subdomain Ωi. For more than two subdomains, more edges respectively, the problem can berewritten accordingly.

3.1.2. Discretization of the Nonlinear Convection-Diffusion Equation

Now, we rewrite the steady-state convection-diffusion equation (3.4) into the mixed form. Thismeans we introduce an additional unknown σ = ∇w such that only first derivatives occur in theproblem:

σ = ∇w ∀x ∈ Ω (3.17)∇ · (f(w)− εσ) = h ∀x ∈ Ω (3.18)

w = g ∀x ∈ ∂Ω . (3.19)

The mixed form is often used for DG methods because second derivatives have strong continuityrequirements. These requirements cannot be guaranteed due to the discontinuities that are allowedin the solution.

Furthermore, the domain Ω is partitioned into N disjoint partitions such that Ω =⋃Ni=1 Ωk with

∂Ωk being the set of edges of the domain Ωk. The set Γ denotes all edges ek that occur in thepartition. This means edges that intersect with the domain boundary ∂Ω and edges of intersectingsub-domains Ωi ∩Ωj 6= ∅, i 6= j. Edges ek between intersecting sub-domains occur only once in the


set Γ. The total number of edges ek is given by N̂ := |Γ|. Additionally we define the set Γ0 ⊂ Γthat contains all edges ek which do not intersect with the domain boundary ∂Ω. All elements andedges have a normal vector n pointing outwards the corresponding domain Ωk. Based on this w

±

is defined as w± = limτ→0+

w(x± τn).For the actual discretization, we need ansatz spaces for the unknowns. We choose

Hh := {f ∈ L2(Ω)|f|Ωk ∈ Πp(Ωk) ∀k = 1, . . . , N}d·m (3.20)Vh := {f ∈ L2(Ω)|f|Ωk ∈ Πp(Ωk) ∀k = 1, . . . , N}m (3.21)Mh := {f ∈ L2(Γ)|f|ek ∈ Πp(ek) ∀k = 1, . . . , N̂ , ek ∈ Γ}m . (3.22)

with d = 2 being the spatial dimension. The space of polynomials up to a degree p on a domain Ωis Πp(Ω) while m is the number of unknowns. In case of the scalar convection-diffusion equationthe number of unknowns is m = 1. The Euler and Navier-Stokes equations in two dimensionshave m = 4 unknowns (density, momentum in both spatial directions and energy). For a moreconvenient notation we introduce a set of abbreviations for integration and summation:

(g, h) :=

N∑k=1

∫Ωk

g · h dx, 〈g, h〉∂Ωk :=N∑k=1

∫∂Ωk

g · h dσ, 〈g, h〉Γ :=N̂∑k=1

∫ek

g · h dσ (3.23)

Then the weak formulation of the hybridized method is:

Method 1 (Hybridized DG-Method) Find a solution (σh, wh, λh) ∈ Hh × Vh ×Mh such that forall test functions (τh, ϕh, µh) ∈ Hh × Vh ×Mh the following equations hold true:

(σh, τh) + (wh, ∇ · τh)−〈g, τ−h · n

〉∂Ωk∩∂Ω −

〈λh, τ

−h · n

〉∂Ωk\∂Ω = 0 (3.24)

− (f(wh), ∇ϕh)− (ε∇σh, ϕh) +〈f(g) · n, ϕ−h

〉∂Ωk∩∂Ω

+〈f(λh) · n− α(λh − w−h ), ϕ−h

〉∂Ωk\∂Ω = (h, ϕh) (3.25)〈

σ−h · n− σ+h · n+ α(2λh − w−h − w+h ), µh〉

Γ0+ 〈λh − g, µh〉Γ\Γ0 = 0 . (3.26)

The analytical flux f(w) has been replaced on edges by a numerical flux, a modified Lax-Friedrichsflux [43]

f(u) · n ≈ fnum(λ, v, n) := f(λ) · n− α(λ− v) (3.27)with a parameter α that is problem dependent.

Since the solution is only required to be continuous on an element Ωk, there may occur two differentsolutions on the edge between two neighboring elements Ωk and Ωk′ . To retain conservation thenumerical flux is introduced. It ensures that the in- or outgoing flux is the same for both elements.It is a common technique borrowed from finite volume methods.

This method does not yet look too interesting because a lot of new unknowns for σ and thefunction λ on the edges have been introduced. However, this formulation can be rewritten in amore convenient way. As already indicated in a previous section, the unknowns (wh, σh) can beexpressed as a function of λ such that (wh(λ), σh(λ)). Then wh(λ) and σh(λ) can be computedonce λ is known.

One way to interpret this, is given by the system of equations. After linearization of the possiblynonlinear problem, one retains a linear system of equations that can be written in the followingform: A B RC D S

L M N

ΣWΛ

=FGH

. (3.28)This is a very general way to express the problem with block-matrices. The vectors Σ,W and Λcorrespond to the unknown basis coefficients of (σh, wh, λh) and the vectors F ,G and H correspond


to the residual, i.e., equations (3.24)-(3.26). The system of equations can be divided into two partsthat are given as (

A BC D

)(ΣW

)=

(FG

)−(RS

)Λ (3.29)

⇒(

ΣW

)=

(A BC D

)−1((FG

)−(RS

)Λ

)(3.30)

and (L M

)(ΣW

)+NΛ = H. (3.31)

Then equation (3.30) can be inserted into equation (3.31) yielding a linear system of equations thatonly depends on the hybrid variable λh:(

N −(L M

)(A BC D

)−1(RS

))Λ = H −

(L M

)(A BC D

)−1(FG

). (3.32)

The block structure of the matrices allows to invert the matrices locally in an efficient manner. Ithas to be noted that the construction of the matrices before the system of equations is solved for Λis more involved than for other methods. It requires solving a problem on each element. Theselocal solves are needed to construct the system of equations and reconstruct the solution once Λ isknown. Depending on the memory available for a simulation, these local solves can be stored for afaster reconstruction. More specific information about the implementation of the used frameworkcan be found in the Master’s thesis of Woopen [77].

After obtaining the vector Λ the remaining unknowns can be computed by solving equation (3.30)in an elementwise fashion. It leads to a larger number of total computations, but the requiredpeak memory is lowered due to the smaller system when solving for the number of globally coupledunknowns.

In some cases the hybridized method actually increases the number of globally coupled unknownscompared to a standard DG method. This happens if the polynomial degree p is small. Then, thereare only few degrees of freedom on each element Ωk. Since the degrees of freedom of the hybridvariable λh exist on the edges and there are usually more edges than elements, this may lead to anincreased number of unknowns in total [59]. Nevertheless, hybridization is usually beneficial. Thecondition of the matrix resulting from the system of equations is often better than for standard DGmethods leading to faster convergence. At the “2nd International Workshop on High-Order CFDMethods” in Cologne, the results have shown that even for ansatz functions with low polynomialorder p the hybridized method can be faster than standard DG methods.

4. Time Integration Methods

In this section the integration methods are explained and a way of time step adaptation is discussed.These methods have been developed for ordinary differential equations (ODE) and we also stay withthese type of problems in this chapter. How they can be applied to the HDG method is describedafterwards in Chapter 5.

First, a short introduction to the general problem setting is given. After that, the key properties– convergence, consistency and stability – of time integration methods are introduced briefly. Thentwo types of schemes, diagonally implicit Runge-Kutta methods (DIRK methods) and backwarddifferentiation formulas (BDF), are presented and especially their stability properties are discussed.For DIRK methods a strategy for time step adaptation is shown.

The information given in this section are taken from the books of Dahmen and Reusken [25],Hairer and Wanner [33], Hairer et al. [32] and Strehmel et al. [67] if nothing else is stated.

4.1. General Setting

The methods that are used in this work were actually developed to solve a first order ordinarydifferential equation with given initial conditions. Such a problem can be formulated as:

Problem 1 Find a function y = y(t), with t being the time variable, satisfying

y′(t) = g(t, y(t)), t ∈ [t0, T ], y ∈ Rn (4.1)

with given initial conditiony(t0) = y

0. (4.2)

The idea of time integration methods is to use the slope of the solution y(t) at certain points thatis given by g(t, y(t)), that has to be Lipschitz continuous, to determine a solution. Although themethods have been developed for ODEs, they are also useful in the case of partial differentialequations. They are often applicable to the system of equations that occurs after the spatialdiscretization of a PDE, not only in the case of hybridized DG methods.

First we consider so-called one step methods. A very famous and also simple method is theexplicit Euler method:

Method 2 (Explicit Euler Method) Given a step size ∆t = T−t0N with N ∈ N. Compute forj = 0, . . . , N − 1

tj+1 = tj + ∆t

yj+1 = yj + ∆tg(tj , yj).

(4.3)

The idea behind the method is to use the slope of the solution given by g(t, y) and construct thetangent at point (tj , y

j). Another important scheme is the implicit or backward Euler method.

Method 3 (Implicit Euler Method) Given a step size ∆t = T−t0N with N ∈ N. Compute forj = 0, . . . , N − 1

tj+1 = tj + ∆t

yj+1 = yj + ∆tg(tj+1, yj+1).

(4.4)

In this case the slope at point (tj+1, yj+1) is used. However, this point is not known beforehand,

as it is actually the point that shall be computed. Therefore, the method is implicit and involvessolving a, possibly nonlinear, system of equations. Although this renders the method much moreexpensive than the explicit method, where the solution can be directly determined from givendata, implicit methods are used for the HDG method. The reason for this is on the one hand side

4. Time Integration Methods 14

the stability properties of implicit methods (cf. Section 4.1.2) which allow larger time steps thanexplicit methods. On the other hand side the HDG method itself is implicit and leads to a set ofdifferential algebraic equations, as already mentioned, to which explicit time integration methodscannot be applied.

4.1.1. Consistency and Convergence

For the analysis of one step methods we introduce an alternative notation, i.e.,

yj+1 := yj + ∆tΦg(tj , yj ,∆t) . (4.5)

The function Φg is called increment function and gives a convenient way to encode the wholemethod. All methods, no matter if explicit or implicit, can be written this way. We assume thatthe time step size ∆t is constant with ∆t = T−t0N , N ∈ N. It is possible to use a varying time step,what is actually our aim in the end. However, a constant time step size makes the notation mucheasier without loss of generality. We come back to variable time steps when we discuss the timestep adaptation (cf. Section 4.2.2). For the explicit Euler method the increment function is

Φg(tj , yj ,∆t) = g(tj , y

j)

while for the implicit Euler method it is

Φg(tj , yj ,∆t) = g(tj + ∆t, y

j + ∆tΦg(tj , yj ,∆t)) .

The time interval [t0, T ] is discretized into a set of discrete points in time t0 < t1 < . . . < tN = Tdefining the temporal grid Gh = {t0, . . . , tN}. We define the grid function yh as the solution thatis computed for a given point in time tj that means yh(tj) = y

j , j = 0, . . . , N . From this we canintroduce the local error in two different ways as in Dahmen and Reusken [25]. It quantifies theerror that occurs when computing a single time step. First, we define the general form of the localerror, also called local error per unit time step,

δ(ta, ya,∆t) = y(ta + ∆t; ta, y

a)− yh(ta + ∆t; ta, ya) (4.6)

with (ta, ya) being initial values in a neighborhood of the global solution and y(t; ta, y

a) is thesolution at time t for the initial value problem

y′(t) = g(t, y), t ∈ [ta, ta + ∆t]y(ta) = y

a.

It is assumed that there exists a unique solution to this problem. Since we determine a solution fora certain set of times tj we take certain points

(ta, ya) = (tj , y(tj)).

from our grid Gh. This leads to a local error of

δj,∆t := y(tj+1; tj , y(tj))− yh(tj+1; tj , y(tj))= y(tj+1)− y(tj)−∆tΦg(tj , y(tj),∆t)

(4.7)

which is actually the difference between the ’true’ solution y(tj+1) and the solution determinedwith the one step method assuming one has the exact value of y(tj) at time tj . It is also possibleto define the local error when doing one time step based on the approximated solution

(ta, ya) = (tj , y

j).

This leads to a second definition of the local error

δ̃j,∆t := y(tj+1; tj , yj)− yh(tj+1; tj , yj)

= y(tj+1; tj , yj)− yj −∆tΦg(tj , yj ,∆t).

(4.8)


Both definitions can be used for analysis of the time integration methods. However, the firstdefinition is used in this section since it allows easier notation. The second one is important whenintroducing error estimators for time step control. With the local error we can establish the notionof consistency. A one step method is consistent iff

lim∆t→0

‖δ(ta, ya,∆t)‖∆t

= 0, ∀ta ∈ [t0, T ]

holds true for all initial value problems (4.1) with a vector norm ‖·‖ on Rn. This can additionallybe used to quantify the quality of the approximation.

Definition 1 (Consistency of one step methods) A one step method is consistent of orderq ∈ N with the initial value problem (4.1), if the relation

‖δ(ta, ya,∆t)‖ ≤ C∆tq+1, ∆t→ 0 (4.9)

holds true for all points (ta, ya) in a neighborhood of the solution {(t, y(t)) | t ∈ [t0, T ]} with a

constant C independent of ∆t.

Besides the local error there exists the global discretization error eh(tj) = y(tj) − yh(tj), j =0, . . . , N . It takes into account that the local errors sum up from t0, . . . , tN . Convergence is definedas

‖eh(tj)‖∞ := maxj‖eh(tj)‖ = 0, ∆t→ 0. (4.10)

The order of convergence can also be determined and is an important quantity that describes howfast the solution converges.

Definition 2 (Convergence of one step methods) A one step method is of order q ∈ N con-vergent if the relation

‖eh‖∞ ≤ C∆tq, ∆t→ 0 (4.11)holds true with a constant C independent of ∆t.

In the end the global error is of most interest because it determines the temporal error added to thesolution. Nevertheless, the knowledge of the local errors is important for the time step size controland there is a useful connection between the consistency and convergence of a one step method.

Theorem 1 Let g(t, y) and Φg(t, y,∆t) be Lipschitz continuous in y. Then from consistency oforder q follows convergence of order q.

This theorem is very useful as it is sufficient to do an analysis of the error in one time step to get aconclusion about the global error.

4.1.2. Stability

As already mentioned only implicit methods are used in this work. This comes from the factthat the discretized equations impose restrictions on the solution technique. One constraint isthe CFL-condition introduced by Courant, Friedrichs and Lewy [23]. For DG methods a modifiedfashion has been introduced by Cockburn et al. [22]

∆t|el|

min (|Ωk|, |Ωk′ |)‖f ′(w) · n‖L∞ ≤ CFL, ∀l = 1, . . . , N̂ .

Here el denotes an edge of the triangulation. It restricts the largest possible time step size forhyperbolic equations. This restriction can be loosened by using of implicit methods. On the otherside, the discretized equations usually form a stiff problem. It is hard to give an exact definitionwhat stiff actually means. Most textbooks refer to Curtiss and Hirschfelder [24] who observed thatexplicit methods had trouble solving a system of ordinary differential equations describing chemicalreactions. The reaction rates of this system differed in several orders of magnitude. Implicitmethods, in their case backward differentiation formulas, where much more successful. Therefore,


they identified stiff equations as ”equations where certain implicit methods [...] perform better [...]than explicit ones” [24].

Further analysis of stiff problems has shown that often the Jacobian J = ∂g∂y of a problem andespecially its eigenvalues λ determine the stiffness of a problem. To show this, we linearize theinitial value problem (4.1) around a smooth solution h(t) such that it reads

y′(t) = g(t, h(t)) +∂g

∂y(t, h(t))(y(t)− h(t)) + . . . . (4.12)

By substituting y(t)− h(t) =: ȳ(t) we obtain

ȳ′(t) =∂g

∂y(t, h(t))ȳ(t) + . . . = J(t)ȳ(t) + . . . . (4.13)

We assume that the Jacobian is constant J(t) = J = const and truncate terms of higher order. Forimproved readability we additionally drop the bars. Then the problem becomes

y′(t) = Jy(t), y(0) = y0. (4.14)

It is obvious that the solution of this problem is given by the exponential function as

y(t) = y0eJt. (4.15)

Now, we can apply one step methods to the problem. In the case of the explicit Euler method (4.3)the solution can be determined by

yj+1 = (1 + hJ)yj =: R(hJ)yj . (4.16)

With R(hJ) we denote the so-called stability function. It plays an important role when analyzingthe time integration methods. Under the assumption that the Jacobian is diagonalizable witheigenvectors v1, . . . , vn the initial data can be written in the eigenvector basis as

y0 =

n∑i=1

civi . (4.17)

By using (4.16) recursively the solution yj+1 can be written in terms of the initial data and thereforecan also be expressed in an eigenvector basis

yj = R(hJ)jy0 =

n∑i=1

R(∆tλi)jcivi (4.18)

with the eigenvalues λi of J . The problem itself now is called stiff in the case that all componentsdecay while the speed of decay varies strongly. In the linear case this means the problem is stiff if

Re(λi) < 0, maxi,j

|λi||λj |� 1, λi ∈ C, i = 1, . . . , n (4.19)

It can be seen that yj stays bounded if z = ∆tλi lies in the stability region

S := {z ∈ C | |R(z)| ≤ 1} . (4.20)

In this case, stiff components of the solution are damped by the stability function R(z). Thisis hardly possible for the explicit Euler method. From equation (4.18) follows that the stabilityfunction is R(z) = 1 + z. Therefore extremely small time steps are needed such that |R(z)| is atleast near to one. This is a problem for all explicit methods making it very hard to apply themto stiff problems. In this case implicit methods are needed. The stability function of the implicitEuler method is R(z) = 11+z and therefore there exist z ∈ C for which the method is stable.


From the fact that we assumed the Jacobian to be diagonalizable it is also possible to decouplethe linearized problem (4.14) into a set of n scalar problems. Based on this results we define ascalar test equation. It is called the Dahlquist test equation given by

y′(t) = λy, y0 = 1. (4.21)

We use this problem for the stability analysis of the methods. The solution to this problem isy = eλt.

The eigenvalues of the linearized problem (4.14) are an important tool when investigating thestiffness of a problem. However, this is not always a sufficient condition. A lot of assumptionshave been made to get to the linearized problem. Furthermore, smoothness of the solution, theintegration interval and the size of the problem have also shown to influence on the stiffness of aproblem [33].

4.1.3. A- and A(α)-Stability

A very important kind of stability, when dealing with stiff problems, is A-stability. A method iscalled A-stable if the whole negative half-plane C− = {z ∈ C | Re(z) ≤ 0} is part of the stabilityregion S. This is the region where also the solution of Dahlquist’s problem (4.21) is stable, meaningeλt decays. Using the definition of the stability function from the previous section a method isA-stable iff

|R(z)| ≤ 1, ∀ z ∈ C− (4.22)holds true. The explicit Euler method with stability function R(z) = 1 + z is therefore not A-stable.In fact, there does not exist any explicit one step methodwhich is A-stable as already mentionedbefore. Only implicit methods can be A-stable. The implicit Euler method, for example, with thestability function R(z) = 11+z is obviously A-stable.

In some cases A-stability is a too strong requirement for a solution technique. It may be sufficientif the method is stable for

”most“, but not all, z of the negative half plane C−. Therefore, also a

weaker form of stability is introduced called A(α)-stability. A method is A(α)-stable if its stabilityregion contains a symmetric domain on the left half-plane

S ⊃ {z ∈ C | | arg(z)− π| < α}. (4.23)

The domain is limited by two axis-symmetric lines with an inner angle of 2α. In the case of α = π2the stability region is the whole negative half-plane z ∈ C−, meaning it is the same as A-stability.This kind of stability is especially interesting for linear multistep methods as is explained in Section4.3.

The stability regions are often displayed using a two dimensional plot. In Figure 4.1 we showsuch plots with the regions that have to be at least in the stability region S for A-stability (left)and A(α)-stability with α = π4 (right).

4.1.4. L-Stability

It has been mentioned that the damping of stiff components is a desired property. A method isdamping if the norm of the stability function is strictly smaller than one, |R(z)| < 1. It may appearthat a method loses its damping property. For example, the trapezoidal rule

yj+1 = yj +∆t

2

(g(tj , y

j) + g(tj+1, yj+1)

)(4.24)

is an implicit one step method with stability function

R(z) =1 + 0.5z

1− 0.5z . (4.25)

It is A-stable since it fulfills |R(z)| ≤ 1,∀z ∈ C−, but it looses its damping property for largeamplitudes because the stability function reaches the stability bound for large arguments z, i.e.,

limz→−∞

R(z) = 1. (4.26)


−3 −2 −1 0 1 2 3−3

−2

−1

0

1

2

3

Re(z)

Im(z)

(a) Stability region: A-stable

α

α

−3 −2 −1 0 1 2 3−3

−2

−1

0

1

2

3

Re(z)

Im(z)

(b) Stability region: A(α)-stable

Figure 4.1.: Stability region: Examples of the stability regions of an A-stable method (left) and anA(α)-stable method with α = π4 (right) are shown.

This may become a problem depending on the problem to be solved because stiff components maybe damped only very slowly. Therefore an additional kind of stability is introduced.

Definition 3 (L-Stability) A one step method is L-stable if it is A-stable and for the stabilityfunction the following holds true:

limz→∞

R(z) = 0. (4.27)

L-stability is not a requirement as important as A-stability for solving stiff problems becauseA-stable methods that are not L-stable can solve many stiff problems successfully. Nevertheless,it may be useful because our solver should be applicable to a wide range of different problems.Therefore, it is possible to have problems for which L-stability is beneficial.

In contrast to one step methods, multistep methods may be L-stable even if they are onlyA(α)-stable. A more detailed discussion on the stability of multistep methods, especially BDFmethods, is given in Section 4.3.1.

4.2. Diagonally Implicit Runge-Kutta Methods

A famous class of one step methods are the so-called Runge-Kutta methods (RK methods). Theidea is that the initial value problem (4.1) can be solved by integrating g(t, y(t)) over time from aknown point (tj , y

j) to a new time tj+1

yj+1 = yj +

∫ tj+1tj

g(x, y(x))dx. (4.28)

The integral is approximated using a quadrature rule∫ tj+1tj

g(x, y(x))dx ≈ ∆ts∑i=1

biki (4.29)

with some weights bi and certain nodes

ki ≈ g(xi, y(xi)), i = 1, . . . , q. (4.30)

All Runge-Kutta methods can be written in the following form:


Method 4 (Runge-Kutta method with s stages) Given weights ci, bi, 1 ≤ i ≤ s and ai,l,1 ≤ i, l ≤ s with step size ∆t = T−t0N . Compute for j = 0, . . . , N − 1:

tj+1 = tj + ∆t (4.31)

ki = g(tj + ci∆t, yj + ∆t

s∑l=1

ai,lkl), i = 1, . . . , s (4.32)

yj+1 = yj + ∆t

s∑l=1

blkl (4.33)

We call s the number of stages. The weights have to be chosen such that a desired stability andconsistency is achieved. A common way to write the weights is the Butcher-Tableau

c1 a1,1 a1,2 . . . a1,s

c2 a2,1. . .

. . ....

......

. . .. . .

...cs as,1 as,2 . . . as,s

b1 b2 . . . bs

(4.34)

where all coefficients are displayed. This can also be viewed as an array of two vectors b, c and amatrix A with

b = (b1, . . . , bs)T , c = (c1, . . . , cs)

T , A = (ai,l)si,l=1 (4.35)

such that the Butcher-Tableau can also be displayed as

c AbT

. (4.36)

For explicit methods the matrix A has a lower triangular structure with all diagonal elementsbeing zero while implicit methods have at least one non-zero element on the diagonal or uppertriangular part of the matrix. In this work we consider of the special case where A has a lowertriangular structure with diagonal entries being non-zero. We call these methods diagonally implicitRunge-Kutta (DIRK) methods. In fact, we focus on methods where all diagonal elements are alike,i.e., a11 = aii = ass = γ 6= 0. Then the Butcher-Tableau becomes

c1 γ 0 . . . 0

c2 a2,1 γ. . .

......

.... . .

. . . 0cs as,1 as,2 . . . γ

b1 b2 . . . bs

(4.37)

and the method is called singly diagonally implicit Runge-Kutta method (SDIRK method).An important work about these methods is from Alexander [2]. The actual naming of these

methods varies depending on the sources and authors. Alexander has called his methods DIRKmethods although they had the same element on the diagonal. Other authors have called thesemethods semi-explicit [16, 48] or semi-implicit [15, 17] Runge-Kutta methods because these methodsare somewhere between the full-implicit and the explicit case.

One advantage of (S)DIRK methods compared to the fully implicit methods is that they requiresolving one system of equations of size n×n in each stage instead of solving one system of equationsof size sn× sn where n is the number of unknowns. Therefore, the required memory to solve thesystem of equations is much lower for (S)DIRK methods compared to full-implicit Runge-Kuttamethods. This especially becomes important for growing problem sizes n and number of stages swhere the memory needed for solving such a system of equations becomes a bottleneck. Moreover,there is the possibility to freeze the Jacobian arising in the linearized system of equations for SDIRKmethods to lower the computational effort because of the constant diagonal term γ. However, we do


0.435866521508 0.435866521508 0 00.717933260755 0.2820667320 0.435866521508 0

1.0 1.208496649 -0.6443632015 0.435866521508bi 1.208496649 -0.6443632015 0.435866521508

b̂i 0.77263013745746 0.22736986254254 0.0

Table 4.1.: Coefficients of Cash’s SDIRK method.

not exploit on this feature in this thesis. This is left open for future work. The actual applicationof the SDIRK methods to the hybridized DG method is discussed in detail in Chapter 5.

A drawback compared to the fully implicit schemes is the lower order of the SDIRK methods.The Runge-Kutta-Gauss methods have order q = 2s while being A-stable, for example, and arefully implicit. An A-stable SDIRK method, however, has maximum order q ≤ s + 1 and for anL-stable method the order is limited to q ≤ s [33].

A huge advantage of all one step methods, that already has been mentioned, is the possibility toadapt the time step size. For this an error estimator is needed. In the case of (S)DIRK methods itis possible to embed an additional (S)DIRK method of different order into an existing method forerror estimation. This way it is possible to get an error estimate with only few extra work.

To refer to the stages and the order of the SDIRK methods we use the notation (s, q) whichrefers to a method of order q with s stages. For the case with an embedded method we write(s1/s2, q1/q2) meaning the method with s1 stages has order q1 and second scheme has s2 stages andorder q2. For example, (2/3, 4/5) would refer to a SDIRK scheme with 3 stages of order 5 where anSDIRK scheme with 2 stages of order 4 is embedded. The Butcher-Tableau of such a scheme canbe displayed as

c1 γ 0 . . . 0

c2 a2,1 γ. . .

......

.... . .

. . . 0cs as,1 as,2 . . . γ

b1 b2 . . . bsb̂1 b̂2 . . . b̂s.

. (4.38)

We denote the coefficients of the higher order method with bi and the coefficients of the embeddedmethod with b̂i. From this, two different solutions can be determined

yj+1 = yj + ∆t

s∑l=1

blkl, (4.39)

ŷj+1 = yj + ∆t

s∑l=1

b̂lkl. (4.40)

The local error then can be estimated using the two solutions to adjust the time step size. Furtherremarks regarding the estimated error and time step size adaptation are given later in Section 4.2.2.

In this work we restrict ourselves to three different SDIRK methods. Cash has derived differentformulas for embedding an error estimating method into a DIRK method [18]. We use Cash’s (2/3,2/3) SDIRK method. It is the method of lowest order in our analysis, but also has the lowestnumber of stages meaning it is possibly less costly for a time step than the other methods. Itscoefficients are given in Table 4.1. Another approach has been used by Al-Rabeh [1]. He hasconstructed a (3,3) method very similar to the one of Cash and then developed a (4,4) SDIRKmethod enveloping the (3,3) formula. The coefficients are given in Table 4.2. The last method is a(5,4) method of Hairer and Wanner [33]. It has a (4,3) method embedded. In contrast to the othermethods it has one additional stage compared to its order q meaning q = s− 1. The other methodshave the same order as number of stages q = s. This may lead to a method with higher amount ofcomputational work due to the additional stage, but also may be covered if the error estimatorworks better than for the other methods.


0.4358665 0.4358665 0 0 00.0323722 -0.4034943 0.4358665 0 00.9676278 -0.3298751 0.8616364 0.4358665 00.5641335 0.5575315 -0.1930865 -0.2361781 0.4358665

bi 0.3153914 0.1846086 0.1846086 0.3153914

b̂i 0.6307827 0.1413538 0.2278634 0

Table 4.2.: Coefficients of Al-Rabeh’s SDIRK method.

14

14 0 0 0 0

34

12

14 0 0 0

1120

1750 -

125

14 0 0

12

3711360 -

1372720

15544

14 0

1 2524 -4948

12516

−8512

14

bi2524 -

4948

12516 -

8512

14

b̂i5948 -

1796

22532 -

8512 0

Table 4.3.: Coefficients of Hairer’s and Wanner’s SDIRK method.

4.2.1. Stability

As already mentioned in Section 4.1.2, stability of the time integration method is a crucial propertyto solve stiff problems. In general, the stability function of a Runge-Kutta method can be computedusing the Butcher-Tableau in matrix vector notation (4.35). The stability function then becomes

R(z) =det(I − zA+ zobT

)det (I − zA) (4.41)

with I being the identity matrix and o := (1, . . . , 1)T is a vector of ones [33]. For SDIRK-methodsthe stability function only depends on the diagonal element aii = γ. Using this knowledge, Hairerand Wanner [33] computed valid values for γ to retain a certain order and stability depending onthe number of stages s. The values for γ to have an A-stable method are given in Table 4.4 and forL-stable methods are given in Table 4.5.

The method of Cash has three stages (s = 3) and γ = 0.435866521508. As it is third-orderconsistent, it is A- and L-stable according to the Tables. The embedded formula with two stagess = 2 is only A-stable.

Al-Rabeh’s embedded method with three stages is A- and L-stable similar to the method ofCash. However, the higher-order method with four stages is only A-stable. Whether this may be adisadvantage for stiff problems is examined for different problems in section 7.

The (5,4) SDIRK method of Hairer and Wanner with γ = 0.25 is A- and L-stable. The embedded(4,3) method, however, is not A-stable since γ < 13 .

For a better understanding, the stability regions are plotted in Figure 4.2. The x-axis representsthe real part of z while the y-axis represents the imaginary part of z. White regions in the plotsshow where R(z) ≤ 1 holds true, meaning the SDIRK method is stable. One can see that the (3,3)method of Cash and Al-Rabeh coincide due to the same γ. The (4,3) method of Hairer and Wannerhas only a small stability region. This, however, is acceptable if this method is only used as anerror estimator. For the other methods it is possible to choose between the higher and lower ordermethod as the error estimator. It depends on how important L-stability is for the problem. In thiswork we always use the lower order method for error estimating.

4.2.2. Time Step Adaption

A key feature of one step methods we utilize in this work, is the possibility to adapt the time step.The structures of flows that have to be resolved during a simulation may change over time. In this


−4 −3 −2 −1 0 1 2 3 4−8

−6

−4

−2

0

2

4

6

8

Re(z)

Im(z)

(a) Stability region: Cash’s (2,2) method

−4 −3 −2 −1 0 1 2 3 4−8

−6

−4

−2

0

2

4

6

8

Re(z)

Im(z)

(b) Stability region: Cash’s (3,3) method

−4 −3 −2 −1 0 1 2 3 4−8

−6

−4

−2

0

2

4

6

8

Re(z)

Im(z)

(c) Stability region: Al-Rabeh’s (3,3) method

−4 −3 −2 −1 0 1 2 3 4−8

−6

−4

−2

0

2

4

6

8

Re(z)

Im(z)

(d) Stability region: Al-Rabeh’s (4,4) method

−4 −3 −2 −1 0 1 2 3 4−8

−6

−4

−2

0

2

4

6

8

Re(z)

Im(z)

(e) Stability region: Hairer’s and Wanner’s (4,3)method

−4 −3 −2 −1 0 1 2 3 4−8

−6

−4

−2

0

2

4

6

8

Re(z)

Im(z)

(f) Stability region: Hairer’s and Wanner’s (5,4)method

Figure 4.2.: Stability region: The stability regions of all SDIRK methods. Stable regions withR(z) ≤ 1 are indicated white while unstable regions are red.


q ≤ s q = s+ 1

s = 1 12 ≤ γ


For the analysis, we define an error estimator s(∆tj) for δ̃j,∆tj such that

δ̃j,∆tj ≈ s(∆tj) ≈ c∆tq+1j . (4.46)

Here q is the order of accuracy of the estimated error. To test if a time step is valid or not, the ratio

r(∆tj) :=s(∆t)

tol ·∆tj(4.47)

is very important. It describes the ratio of the estimated error to the error that is allowed for atime step. In the case r(∆tj) > 1 the error in the current time step is too large and it has to berecomputed with a smaller time step size ∆tj . If r(∆tj) ≤ 1 the time step is accepted. Dependingon the actual value of r(∆tj) the time step size ∆tj+1 may be increased for the next time step as itis preferable to have r(∆tj) close to one. Therefore, we require for a time step

1!≈ c∆t

q+1new

tol ·∆tnew. (4.48)

When we cancel ∆tnew from the denominator and multiply with the old time step size ∆t we obtain

c∆tq+1newtol ·∆tnew

=c∆tq+1

tol ·∆t

(∆tnew

∆t

)q. (4.49)

Since our estimator shall be s(∆tj) ≈ c∆tq+1j (cf. eq. (4.46)), we finally get

s(∆t)

tol ·∆t

(∆tnew

∆t

)q= r(∆t)

(∆tnew

∆t

)q. (4.50)

This expression still is required to be close to one. Therefore, the new time step size ∆tnew iscomputed as

r(∆t)

(∆tnew

∆t

)q≈ 1⇒ ∆tnew ≈ r(∆t)−

1q ∆t. (4.51)

In our case we use the embedded formula to compute two solutions yj+1 (cf. eq. (4.39)) and ŷj+1

(cf. eq. (4.40)). The local error is estimated as

‖δ̃j,∆tj‖ ≈ s(∆tj) :=‖ŷj − yj‖

=‖yj + ∆tjs∑l=1

b̂lkl −(yj + ∆tj

s∑l=1

blkl

)‖

=‖∆tjs∑l=1

(b̂l − bl)kl‖.

(4.52)

Additionally, one can use information from the solution process of the nonlinear system of equations.The number of iterations nit needed to solve the system of equations using Newton’s methoddecreases, the better the starting value and the condition of the system of equations is. If thesolution changes only slowly, then Newton’s method is more likely to converge after a few steps.The opposite happens in the case of too large time steps or a rapidly changing solution. Therefore,we introduce information about the solution process into the time step adaptation. For this wedefine an additional safety factor α that depends on the highest number of iterations of all stagesnit and the maximum allowed number of iterations for the Newton solver nit,max as:

α(nit) = 0.92nit,max + 1

2nit,max + nit. (4.53)

If nothing else is stated we set nit,max := 10 such that 0.63 ≤ α(nit) ≤ 0.9. Combining equation(4.51) and the safety factor (4.53) we set the new time step size using the following formula:

∆tnew := α(nit)r(∆tj)− 1q ∆tj . (4.54)


a0 a1 a2 a3 a4

BDF1 −1.0 1.0 0.0 0.0 0.0BDF2 12 −2.0 32 0.0 0.0BDF3 − 13 32 −3.0 116 0.0BDF4 14 − 43 3.0 −4.0 2512

Source: Textbook of Dahmen and Reusken [25]

Table 4.6.: Coefficients of BDFk-methods up to order 4.

Even when setting the tolerance tol carefully, it may occur that the time step size grows too largeto resolve all interesting features of a solution. On the other hand it may happen that the timestep size gets ridiculously small, e.g. when a boundary layer evolves. To maintain some certaintyabout the solution and the time needed to run a simulation we specify an upper ∆tmax and lowerbound ∆tmin for the time step. In each time step it is checked if the time step size stays withinthese bounds ∆tmin ≤ ∆tj ≤ ∆tmax. If this is not the case the time step size ∆tnew is rejected andset to ∆tmin or ∆tmax corresponding to the bound reached. The full time step adaptation processis also displayed in a flow chart, see Fig. 4.5.

4.3. Backward Differentiation Formulas

Another popular class of time integration methods are linear multistep methods. As the namesuggests these methods approximate a solution of (4.1) based on the solution of one or moreprevious time steps. It is possible to use an adaptive time step for BDF methods like for one stepmethods [32]. However, BDF methods have no integrated error estimator such that one cannotdecide whether a time step is adequate. Therefore, we only use multistep methods with a fixedtime step size ∆t = T−t0N . The general form of a k-step method can be written as

yj+k = Φh(tj+k−1, yj , yj+1, . . . , yj+k), j = 0, . . . ,m− k. (4.55)

Then the method is given as

k∑l=0

alyj+l = ∆t

k∑l=0

bl · g(tj+l, yj+l) (4.56)

where yj+k is the value that has to determined. From this follows that ak 6= 0 is required. Thecoefficients al and bl are fixed and depend on the method. Multistep methods can be implicit(bk 6= 0) or explicit (bk = 0). Similar to the one step methods, explicit methods have poor stabilityproperties. Therefore, we focus on a class of implicit methods called backward differentiationformula. It is a very popular class of multistep methods for stiff problems. A short form to refer tothese methods is BDFk-method where k is the number of previous steps needed for method. Atthe same time k is the order of consistency and convergence of the method.

Method 5 (BDFk-method) Given a step size ∆t = T−t0N with N ∈ N, coefficients al (0 ≤ l ≤ k)and initial values y0, . . . , yk−1. Compute for j = 0, . . . , N − k

tj+1 = tj + ∆t

k∑l=0

alyj+l = ∆t · g(tj+k, yj+k)

(4.57)

The coefficients of the BDFk-methods up to fourth order are given in Table 4.6. They are determinedby polynomial interpolation using the Lagrange polynomials with the data points

(tj , yj), (tj+1, y

j+1), . . . , (tj+k, yj+k). (4.58)


t0 t0 +∆t t0 + 2∆t t0 + 3∆t

t

∆t

∆t̃

impl. Euler BDF2 BDF3

y0 y1 y2

Figure 4.3.: Initialization of the time integration for the BDF3 method.

The resulting polynomial pk(t) then is differentiated and required to equal the right hand side ofthe initial value problem (4.1) such that

p′k(tj+k) = g(tj+k, yj+k) (4.59)

holds true.A problem occurs for k > 1 since only initial data y0 are given, but for the computations data

at the points yj , j = 0, . . . , k − 1 are needed. The idea is to compute the missing values usingalternative methods. They have to be computed with sufficient accuracy because the order ofconsistency and convergence misbehaves if the error in the initial data is too large.

For the BDF2 method one implicit Euler step is computed with time step size ∆t to compute y1.The implicit Euler method is of second order consistent for one step and therefore preserves theneeded accuracy. For higher order methods another initialization procedure is needed. We sketch itfor the BDF3 method (cf. Figure 4.3):

1. Determine a new time step size ∆t̃ that fulfills ∆t̃2 < ∆t3 and ∆t = l ·∆t̃ for an l ∈ N.

2. Compute one implicit Euler step using ∆t̃.

3. Compute BDF2 steps until t1 = t0 + ∆t and t2 = t0 + 2∆t to obtain y1 and y2.

After the starting values y0, y1 and y2 are computed the BDF3 method is used.An alternative approach is to use one step methods to compute starting values. Depending

on the order of the one step method the time step size has to be adjusted, but it may yield lesscomputational costs than initializing with lower order BDF methods. This is not investigated inthis work as it is not directly obvious whether the computational work is actually reduced. Theamount of computational work depends on the time step size that has to be used for initializationand error constants of the methods, for example.

4.3.1. Stability

In contrast to one step methods, an additional kind of stability is needed for multistep methodssuch that stability and consistency lead to a convergent method. For this the test equation

y′(t) = 0, y(0) = 1 (4.60)

is defined. Obviously, the solution is y(t) = 1. Applying a multistep method to this problem leadsto the equation

k∑l=0

alyj+l = 0, j = 0, 1, . . . . (4.61)

meaning that the solution yj fulfills a homogeneous difference equation. In order to determine thesolution of this equation we use the ansatz yj = ξj with j being the jth power of ξ ∈ C, ξ 6= 0.Then, the so-called characteristic polynomial ρ(ξ) is

ρ(ξ) =

k∑l=0

alξl = 0 (4.62)


where we have already divided the equation by ξj . The equation is fulfilled if ξ is a root of theequation. Using all roots ξ1, . . . , ξk of ρ(ξ), the general solution of equation 4.61 can be written as

yj = p1(j)ξj1 + . . .+ pk(j)ξ

jk (4.63)

where p1(j), . . . , pk(j) are polynomials of degree m− 1 with m being the multiplicity of the root. Inorder to obtain a stable solution we require yj to be bounded for j →∞. This is fulfilled for |ξl| ≤ 1if all ξl are single roots. For a multiple root |ξl| the series jiξjl with i = 1, . . . ,m is a solution, too.Since ji grows unbounded for i > 1 the root must fulfill |ξl| < 0 in this case to obtain a boundedsolution. This requirement leads to zero-stability.

Definition 4 (Zero-Stability) A method is zero-stable if the root condition holds true: In caseρ(ξl) = 0 for a ξl ∈ C, then

|ξl| ≤ 1 (4.64)and if ξl is a multiple root then

|ξl| < 1 (4.65)is fulfilled.

Linear one step methods are always zero-stable. Applying a one step problem to equation 4.60 andusing the ansatz yj = ξj again gives the characteristic polynomial

ρ(ξ) = 1.0ξ0 − 1.0ξ1 = 0 (4.66)

with only one root ξ1 = 1.Now, with consistency and zero-stability it is possible to formulate the following theorem:

Theorem 2 If a multistep method is zero-stable and consistent of order q then it is convergent.The order of convergence is also q.

For the stability concerning stiff problems we use again Dahlquist’s test equation (4.21). Applyinga multistep method on this problem gives

k∑l=0

alyj+l = ∆tλ

(k∑l=0

blyj+l

). (4.67)

In contrast to zero-stability the right hand side does not vanish. The problem can be rewritten as

k∑l=0

((al − blz)yj+l

)= 0, z = λ∆t. (4.68)

Now, we use the ansatz yj = ξj again, where ξj is still the jth power of ξ, and divide by ξj suchthat we obtain the characteristic equation

k∑l=0

((al − blz)ξl

)= ρ(ξ)− zµ(ξ) = 0 (4.69)

with ρ(ξ) as defined in equation 4.62 and

µ(ξ) =

k∑l=0

blξl. (4.70)

The polynomials ρ(ξ) and µ(ξ) are also called the generating polynomials of the method. Asfor the zero-stability, the solution can be expressed as a linear combination (4.63) of the rootsξl, l = 1, . . . , k of equation (4.69). In this case the roots also depend on the complex number z suchthat ξl = ξl(z). So, the solution is, again, stable if the solution y

j stays bounded for j →∞. Thisleads to the same conditions for a stable method as for the zero stability. The solution is stable ifthe norm of all roots are less or equal than one (|ξl(z)| ≤ 1) or strictly less than one (|ξl(z)| < 1)in the case of multiple roots. From this a more general stability concept is derived than includeszero-stability:


−4 −3 −2 −1 0 1 2 3 4−8

−6

−4

−2

0

2

4

6

8

Re(z)

Im(z)

(a) Stability region: BDF1 method

−4 −3 −2 −1 0 1 2 3 4−8

−6

−4

−2

0

2

4

6

8

Re(z)

Im(z)

(b) Stability region: BDF2 method

−4 −3 −2 −1 0 1 2 3 4−8

−6

−4

−2

0

2

4

6

8

Re(z)

Im(z)

(c) Stability region: BDF3 method

−4 −3 −2 −1 0 1 2 3 4−8

−6

−4

−2

0

2

4

6

8

Re(z)

Im(z)

(d) Stability region: BDF4 method

Figure 4.4.: Stability region: The stability regions of BDFk methods for k = 1, . . . , 4. Stable regionswith R(z) ≤ 1 are indicated white while unstable regions are red.

BDF1 BDF2 BDF3 BDF4

α 90.00◦ 90.00◦ 86.03◦ 73.35◦

Source: Textbook of Hairer and Wanner [33]

Table 4.7.: Angle α of A(α)-stable BDFk-methods up to order 4.

Definition 5 (Stability of Multistep Methods) The characteristic equation of a linear multi-step method

ρ(ξ)− zµ(ξ) = 0 (4.71)has roots ξl depending on z ∈ C. The stability domain of the method is defined asS = {z ∈ C; for all roots |ξl| ≤ 1 holds true, for all multiple roots |ξl| < 1 holds true } . (4.72)

For z = 0 characteristic equation of the method becomes ρ(ξ) = 0 which is actually the conditionfor zero-stability. Therefore, a multistep method is zero-stable if 0 ∈ S.

Using this definition of stability, one can transfer the A- and A(α)-stability conditions that havebeen already introduced for one step methods (cf. eq. (4.22) and (4.23)) to multistep methods. Alinear multistep method is A-stable if the negative complex half plane is part of the stability region

C− ⊂ S. (4.73)


The method is A(α)-stable if

{z ∈ C−| |arg(z)− π| ≤ α} ⊂ S (4.74)

holds true. The stability regions for multistep methods can be displayed in the same way as for theone step methods. In Figure 4.4 the stability regions of BDFk methods are shown for k = 1, . . . , 4where red regions denote unstable values for z. Therefore, the BDF methods are only A-stableup to order q = 2 and A(α)-stable for higher orders with angles α given in Table 4.7. The angleα decreases with increasing order of the method. Therefore, BDF methods of order larger thanq = 3 are rarely used. Actually, for k ≥ 7 BDF methods are not zero-stable and therefore are notconvergent.

The observation that multistep methods of order q > 2 are not A-stable holds true for allmultistep methods, not only BDF methods. It is also known as the second Dahlquist barrier.Nevertheless, it is possible to construct multistep methods of hi

UHasselt - A Hybridized Discontinuous Galerkin Method for Time … · 2017. 2. 3. · Dezember 2013...

Documents

Transcript of UHasselt - A Hybridized Discontinuous Galerkin Method for Time … · 2017. 2. 3. · Dezember 2013...