Gryphon - a Module for Time Integration of Partial Differential Equations in FEniCS

Gryphon - a Module for Time Integration of Partial Differential Equations in FEniCS

Knut Erik Skare

Master of Science in Physics and Mathematics Supervisor: Anne Kværnø, MATH Submission date: June 2012

Norwegian University of Science and Technology

Sammendrag

Denne oppgaven tar sikte på å implementere tidsintegratorer i FEniCS-rammeverket.

Mer spesifikt går oppgaven ut på å velge egnede tidsintegratorer, implementere disse og verifisere at de virker ved å anvende dem på et utvalg relevante testproblemer. Dette arbeidet resulterte i en modul til FEniCS som fikk navnet Gryphon. Oppgaven er delt inn i fire deler.

Del I bygger et teoretisk rammeverk som motiverer hvorfor ESDIRK-metoder (Singly Diagonally Implicit Runge-Kutta method with an Explicit first stage) er gode løsere for systemer av stive ordinære differensialligninger (ODEer). Det vil også bli vist hvordan en ESDIRK metode kan brukes til å løse tidsavhengige partielle differensialligninger (PDEer) ved å løse det semi-diskretiserte systemet som oppnås ved å først anvende en endelig elementmetode. Vi vil begrense oss til PDEer som enten semi-diskretiseres til et rent ODE-system eller et differensial-algebraisk system (DAE) av indeks 1.

Del II tar for seg implementasjonen av Gryphon, med fokus på nytteverdi og kodestruk- tur.

Del III tar for seg numeriske eksperimenter på ESDIRK-metodene som er implementert i Gryphon. Eksperimentene vil etablere konvergens og gi kjøretidsresultater for ulike ESDIRK-metoder. Vi vil også se atL-stabilitet er en nyttig egenskap når en jobber med stive ligninger, ved å sammenligne en ESDIRK metode med trapesmetoden. Det blir også verifisert at skrittlengde-kontrollerne implementert i Gryphon oppfører seg som for- ventet. Som testproblemer vil vi se på varmeligningen, Fisher-Kolmogorov-ligningen, Gray-Scott-ligningene, Fitzhugh-Nagumo-ligningene og Cahn-Hilliard-ligningen.

Del IV er en brukermanual for Gryphon hvor alle parameterne brukeren kan endre vil bli forklart. Manualen inneholder også kode for å løse varmeligningen, Gray-Scott- ligningene og Cahn-Hilliard-ligningen, for å hjelpe leseren i gang med å løse egne prob- lemer.

Abstract

This thesis aims to implement time integrators in the FEniCS framework. More specifi- cally, the thesis focuses on selecting suitable time integrators, implement these and verify that the implementation works by applying them to various relevant test problems. This work resulted in a module for FEniCS, named Gryphon. The thesis is divided into four parts.

The first part builds a theoretical framework which will motivate why singly diagonally implicit Runge-Kutta methods with an explicit first stage (ESDIRKs) should be consid- ered for solving stiff ordinary differential equations (ODEs). It will also be shown how an ESDIRK method can be utilized to solve time dependent partial differential equa- tions (PDEs) by solving the semidiscretized system arising from first applying a finite element method. We will restrict our attention to PDEs which either give rise to a pure ODE system or a DAE (differential-algebraic equation) system of index 1.

The second part discusses the implementation of Gryphon, focusing on why such a mod- ule is useful and how the source code is structured.

The third part is devoted to numerical experiments on the ESDIRK solvers implemented in Gryphon. The experiments will establish convergence and give some run-time statis- tics for various ESDIRK schemes. We will also see thatL-stability is a favorable trait when working with stiff equations, by comparing an ESDIRK method to the trapezoidal rule. It will also be verified that the step size selectors implemented in Gryphon behaves as expected. As test problems we consider the heat equation, the Fisher-Kolmogorov equation, the Gray-Scott equations, the Fitzhugh-Nagumo equations and the Cahn-Hilliard equations.

The fourth part is a user manual for Gryphon. All the parameters which can be changed by the user are explained. The manual also includes example code for solving the heat equation, the Gray-Scott equations and the Cahn-Hilliard equation, to get the reader starting on solving their own problems.

Preface

This thesis completes my master’s degree in Industrial Mathematics at the Norwegian University of Science and Technology (NTNU). The work has been carried out at the Department of Mathematical Sciences during the spring semester of 2012 under the supervision of Anne Kværnø.

I would like to thank Anne for giving me an interesting assignment which allowed me to combine my passion for programming with a mathematical application, and for always having an open door to help and discussions.

I would also like to thank Anders Logg, Marie Rognes and Benjamin Kehlet at Simula Research Laboratory in Oslo for valuable input on my work and for giving me insight in how to use FEniCS.

Finally I would like to thank my friends and family for providing me with both moral support and distractions when needed. A special thanks goes to my class mate Olav Møyner for many a useful discussion on programming and numerics.

Knut Erik Skare June 11, 2012 Trondheim

Contents

I Theoretical Background 6

1.1 Introduction . . . 7

1.2 Applying a Runge-Kutta method . . . 7

1.3 Stiff Equations . . . 8

1.4 Stability Properties . . . 9

1.5 Classification of Runge-Kutta Methods . . . 11

1.6 Runge-Kutta pairs . . . 14

1.7 Adaptive Step Size Selection . . . 15

1.8 Differential Algebraic Equations . . . 19

1.9 Applying a Finite Element Method . . . 20

II Implementation 23

2.2 The FEniCS Project . . . 24

2.3 The Gryphon Module . . . 26

2.3.1 The ESDIRK class . . . 29

2.4 A Code Example . . . 31

III Experiments 34

3.2 Test Cases . . . 35

3.2.1 Case 1: The Heat Equation . . . 36

3.2.2 Case 2: The Fisher-Kolmogorov Equation . . . 36

3.2.3 Case 3: The Gray-Scott Model . . . 38

3.2.4 Case 4: A FitzHugh-Nagumo Reaction-Diffusion model . . . . 39

3.2.5 Case 5: The Cahn-Hilliard equation . . . 42

3.3 Verification of Convergence . . . 42

3.4 Verification of Constructive Step Size Selection . . . 46

3.5 Run Time Statistics . . . 46

3.6 Comparison to Trapezoidal Rule . . . 47

IV Gryphon User Manual 55

4.2 Handling Explicit Time Dependency . . . 56

4.3 Solver: ESDIRK . . . 56

4.3.1 Parameters . . . 57

4.3.2 Example Output . . . 62

4.4 Example Problems . . . 63

A Appended Code and Documentation 71

Part I

Theoretical Background

1.1 Introduction

The purpose of this part is to give an introduction to what a Runge-Kutta (RK) method is and how it can be applied to solve ordinary differential equations (ODEs) and index 1 differential algebraic equations (DAEs). It will also be shown how Runge-Kutta meth- ods can be extended to solve partial differential equations (PDEs) by applying them to a semidiscretized system. Different classes of Runge-Kutta methods are presented along- side with their advantages and disadvantages related to computational cost, order and stability properties. This will build the theoretical framework which in turn will mo- tivate why singly diagonally implicit Runge-Kutta methods with an explicit first stage (ESDIRKs) should be considered when working with stiff equations.

1.2 Applying a Runge-Kutta method

A Runge-Kutta method is a one-step time integration scheme which can be applied to approximate a system of ODEs. Given the following system,

d

dty(t) = f(t,y), y(0) =y₀, y∈R^m, ans-stage Runge-Kutta method is applied by setting

yn+1=yn+∆t

∑

biY˙i, Y˙i= f(t_n+ci∆t,Y_i) Y_i=y_n+∆t

∑

a_{i j}Y˙_j, i=1, ...,s,

wherea_{i j},b_iandc_iare coefficients which define the method applied and∆t is the step size. It is common to characterize a Runge-Kutta method by a table called the Butcher tableau. The structure of such a tableau can be seen in table 1.1.

c₁ a₁₁ a₁₂ . . . a_1s

c₂ a₂₁ a₂₂ . . . a_2s

... ... ... ...

c_s a_s1 a_s2 . . . a_ss

b₁ b₂ . . . b_s

Table 1.1: A Butcher tableau used to characterize Runge-Kutta methods.

When deciding upon a Runge-Kutta method to solve a problem, we are first and foremost interested in efficiency. We want to solve our problem within reasonable time to a certain accuracy without using too much computing power. Depending on the nature of our problem, we may be forced to use computationally demanding implicit methods or we may get satisfactory results using relatively cheap, explicit methods.

1.3 Stiff Equations

To give a precise mathematical definition to the phenomena of stiff equations, has shown to be a difficult task. Instead it is helpful to talk about qualitative features which can help us decide whether or not our problem is stiff. In general, stiff problems are known to cause poor performance in explicit Runge-Kutta solvers, meaning that different solvers may give different answers or that the solution grows exponentially. The remedy for this is to use implicit solvers which tend to handle such problems a lot better. This behavior can be linked to stability properties which will be discussed in the next section.

Another source for stiffness comes from semidiscretizing a PDE with high order spa- tial derivatives. The stiffness of the problem is then related to the eigenvalues of the corresponding ODE/DAE system differing in several orders of magnitude as the spatial discretization becomes more refined.

(a) Explicit solver (b) Implicit solver

Figure 1.1: Solution of the heat equation.

A visual example of how a stiff problem behaves when subjected to an implicit method versus an explicit method, can be seen in figure 1.1. The problem solved is the heat equation with homogeneous Dirichlet boundary conditions and a simple Gauss pulse as initial condition. The implicit solver is able to maintain the smooth solution while the explicit one gets oscillations which eventually causes the solution to grow unbounded.

It should be noted that the behavior in figure 1.1a takes place after just five time steps.

After five more, the solution have diverged into meaningless data. For larger step sizes, the effect manifests quicker.

1.4 Stability Properties

The material for this section was found in [HW10, IV.3], where a more elaborate pre- sentation on the following topics can be found.

When discussing stability properties of a Runge-Kutta method, it is useful to consider how it behaves applied to a linear test problem given as

d

dty=λ^{y (1.1)}

whereλ is often called the stiffness parameter. This equation is often referred to as the Dahlquist test equation, and any Runge-Kutta method (explicit or implicit) applied to it, can be written as

y_n+1=R(λ∆t)y_n, R(λ∆t) =P(λ∆t)

Q(λ∆t), P,Q∈P.

The functionR(λ∆t), commonly referred to as the stability function, can tell us a great deal about the stability properties of the method in question. It takes the form of a rational function defined in the complex plane and can be written as

R(z) = P(z)

Q(z)=1+zb^T(I−zA)⁻¹1, z=λ∆t, (1.2) whereAandbconstitute the Butcher tableau and1is a vector of just ones. This stability function also satisfies the following expression

R(z) =det(I−zA+z1b^T)

det(I−zA) (1.3)

which is an easier expression to analyze than (1.2). The stability region, which is where our method is expected to produce stable solutions, is defined to be the areas where

|R(z)|<1. If this is a very narrow region and our problem is very stiff due toλ 0, we may have to choose unreasonably small step sizes in order for our solution to remain stable. An example of this behavior can be seen by considering the explicit Euler method applied to (1.1). The stability function will in this case amount to

R(z) =1+z.

The stability region, where|R(z)|<1, will amount to a circle of radius 1 centered in the point(−1,0)in the negative complex half plane. We now see that if we have a very large, negative value forλ, we have to choose a correspondingly small value for∆t in order to stay within the stability region since we get the bound

|R_euler(λ∆t)|<1=⇒0<∆t<−2

λ^, λ^<0,

which is an increasingly small interval as|λ|grows. Methods which do not suffer from this behavior are said to beA-stable, which are methods whereC⁻is contained in the stability region. This criterion is satisfied if and only if

|R(iy)|<1, y∈R and

R(z)is analytic for Re(z)<0.

A slightly weaker criterion thanA-stability isA(α⁾stability. This concept arose from the observation that methods which are not completelyA-stable, still may be decent methods. A method is said to beA(α^{)stable if}

S_α={zsuch that|arg(z)| ≤α^,z6=0} is contained in the stability region.

We may however still run into problems if we letzapproach the real axis with a very large negative real value (very stiff problem). In this situation,|R(z)|is very close to one which will cause bad convergence for our method due to the stiff parts being damped out very slowly. This motivates the concept ofL-stability which is defined to beA-stable methods which also satisfy

z→∞limR(z) =0.

To further understand whyL-stable methods are beneficial, we will quote some results from [HW10, IV.15]. A more sophisticated test equation than the Dahlquist equation (1.1), is the Prothero-Robinson equation given as

y⁰=λ^(y−γ^{(x)) +}γ^{0(x), y(x}0) =γ^(x0), Re(λ^{)<0, (1.4)} with analytical solutiony(x) =γ(x). The constantλ is still the stiffness parameter.

Applying a Runge-Kutta to (1.4) yields Y_i=y₀+∆t

∑

a_{i j}[λ^(Yi−γ^(x0+c_i∆t)) +γ^0(x0+c_i∆t)], (1.5) y₁=y₀+∆t

∑

b_i[λ^(Yi−γ^(x0+c_i∆t)) +γ^0(x0+c_i∆t)]. (1.6) By inserting the analytical solution fory, we get

γ^(x0+c_i∆t) =γ^(x0) +∆t

∑

a_{i j}γ^0(x0+c_i∆t) +∆i,∆t(x₀), (1.7) γ^(x0+∆t) =γ^(x0) +∆t

∑

b_iγ^0(x0+c_i∆t) +∆0,∆t(x₀), (1.8)

where∆i,∆t(x₀),∆0,∆t(x₀)are known as the numerical defect. By doing Taylor expansion of the above statements it can be shown that

∆0,∆t(x₀) =O(∆t^p+1), ∆i,∆t(x₀) =O(∆t^q+1),

wherepis the order andqis the stage order of the Runge-Kutta method in question. We can now express the error by taking the difference between (1.5)-(1.6) and (1.7)-(1.8):

y₁−γ^(x0+∆t) =R(z)(y₀−γ^(x0))−zb^T(I−zA)⁻¹∆∆t(x₀)−∆0,∆t(x₀)

wherez=λ∆t,R(z)is the stability function,∆∆t(x₀) = (∆1,∆t(x₀), . . . ,∆s,∆t(x₀)). If we replacex₀withx_nwe get the recursion

y_n+1−γ^(xn+1) =R(z)(y_n−γ^(xn)) +δ∆t(x_n) (1.9) where

δ∆t(x_n) =−zb^T(I−zA)⁻¹∆∆t(x_n)−∆0,∆t(x_n).

While we can not control the latter term of (1.9), the first term can be controlled by imposingL-stability (R(∞) =0), making the term vanish asymptotically.

1.5 Classification of Runge-Kutta Methods

This section will present some classes of Runge-Kutta methods. We will focus on which problems the methods can be applied to versus computational complexity. The material regarding DIRK/SDIRK methods was found in [KNO96].

Explicit Runge-Kutta methods (ERKs)

An explicit Runge-Kutta method is best characterized by the Butcher tableau being lower triangular and thus take the form found in table 1.2. These methods have very low computational cost since all the stage values can be expressed explicitly. To perform a time step we only have to evaluate the right hand side in different points, rather than having to solve linear/nonlinear equations. The drawback of using explicit RK methods can be seen by considering the stability function. It takes the form of a polynomial,

R(z) =P(z) =1+O(z),

implying that these methods can never beA-stable. From a practical viewpoint, this is enough to classify explicit Runge-Kutta methods as poor candidates when working with stiff equations, and the methods should not be expected to perform well if extended to PDE solvers. There do however exist techniques where explicit methods can be used to give satisfactory results. See [EJL03] for an example.

0 0 c₂ a₂₁ 0 c₃ a₃₁ a₃₂ 0

... ... ... . ..

cs as1 as2 a_s,s−1 0

b₁ b₂ . . . b_s−1 b_s

Table 1.2: Butcher tableau for an explicit Runge-Kutta method.

Implicit Runge-Kutta methods (IRKs)

An implicit Runge-Kutta method has a Butcher tableau for which one or more stage value is implicitly dependent on any of the other. In the extreme case of a full implicit Runge-Kutta method (FIRK), all the stage values are implicitly dependent on all the other, meaning that performing one time step will involve solvingscoupled systems of equations.

These methods can however be constructed to incorporate stability properties likeA- and L-stability, making them able to handle stiff problems.

Diagonally Implicit Runge-Kutta methods

A subclass of the implicit Runge-Kutta methods are the diagonally implicit Runge-Kutta methods, or DIRKs. The Butcher tableau for this class of methods takes the form found in table 1.3, where at least one of the diagonal elementsa_iimust be nonzero. Compu-

c₁ a₁₁ c₂ a₂₁ a₂₂

... ... ... . ..

c_s a_s1 a_s2 . . . a_ss

b1 b2 . . . bs

Table 1.3: Butcher tableau for a diagonally implicit Runge-Kutta method.

tationally speaking, these methods are less demanding than a full implicit Runge-Kutta method since each stage is only dependent on itself and previous stages, if we start by solving from the top.

Singly Diagonally Implicit Runge-Kutta methods

A subclass of DIRK methods known as SDIRK methods (Singly Diagonally Implicit Runge-Kutta methods) have the property that

a_ii=γ^{, i=1, . . . ,s.}

These methods have a significant computational advantage over DIRK methods since the Newton matrix will be the same for all stages. This can be seen by realizing that solving for each stage value will be a problem on the form

Y_i−∆tγ^Y˙i=y_n+∆t

i−1

∑

i=1

a_{i j}Y˙_i, for which the Newton matrix is

I−∆tγ^fy,

wheref is the right hand side of your problem. Successful realizations of SDIRK meth- ods include SIMPLE by Nørsett and Thomsen [NT84] as well as SDIRK4 by Hairer and Wanner [HW10].

The main restrictions on SDIRK methods is their relatively low order (at most order s+1 forA-stable methods) as well as suffering from order reduction when applied to very stiff problems.

Singly Diagonally Implicit Runge-Kutta methods with an Explicit First Stage

A singly diagonally implicit Runge-Kutta method with an explicit first stage (ESDIRK) has the same format as an SDIRK method with the exception that the first row in the Butcher tableau is equal to zero. The structure of this table can be found in table 1.4.

This class of methods have been studied by Anne Kværnø [Kvæ04], and in her article several ESDIRK methods are developed and presented as adaptive Runge-Kutta pairs of orderpandp−1. Each pair gives rise to two different methods depending on whether or not local extrapolation is used.

To differentiate between the advancing method and the error estimating method, the hat-symbol will be used to mark the error estimating methods, i.e. R(x)is the stabil- ity function of the advancing method while ˆR(x)is the stability function of the error estimating method.

The methods presented in [Kvæ04] are constructed according to the following criteria:

0 0

c₂ a₂₁ γ

c₃ a₃₁ a₃₂ γ

... ... . ..

... ... . ..

c_s−2 a_s−2,1 a_s−2,2 a_s−2,3 . . . . . . γ

1 a_s−1,1 a_s−1,2 a_s−1,3 . . . . . . a_s−1,s−2 γ

1 a_s1 a_s2 a_s3 . . . . . . a_s,s−2 a_s,s−1 γ^.

Table 1.4: Butcher tableau for a singly diagonally implicit Runge-Kutta method with an explicit first stage.

1. Stiff accuracy in both the advancing and the error estimating methods.

2. R(∞) =0, and|R(ˆ ∞)|as small as possible, at least less than one.

3. A-stability, or at leastA(α⁾stability for both methods.

4. As high stage order as possible.

Stiff accuracy in the advancing method is a well known remedy for the problem of order reduction SDIRK methods may experience. It is however not as common to incorporate stiff accuracy in the error estimator, causing the order of the error estimate to be lower than expected. In this sense, the methods developed by Kværnø stands out by requiring stiff accuracy in the error estimator as well.

A-stability together with the requirementR(∞) =0 makes the methodsL-stable which is beneficial when working with very stiff problems. Computationally speaking we nat- urally inherit all the benefits of SDIRK methods in addition to the local error estimate being reduced to

le=Y_s−1−Y_s.

1.6 Runge-Kutta pairs

We will now discuss two strategies for selecting step size when integrating ODE/DAE systems over a given time interval. The simplest strategy is to specify a fixed step size which the program will use until it reaches the end of the domain. This can be said to be a reasonable approach if the behavior of your problem in the desired time interval is

relatively uniform. If this is not the case, allowing your method to change time steps adaptively becomes advantageous. We can save computational power by doing large time steps when the solution is slowly varying and use that to take smaller time steps when the solution is varying more rapidly.

In order to determine the size of the next time step, the program needs an estimate for the local error in the current time step. The user can then specify a tolerance which the new step size is adjusted according to. To allow for such functionality, Runge-Kutta methods are often presented in pairs of orderp andp−1. The idea is that one method is used as an advancing method while the other one is used as a reference to measure the error against. If the method of highest order is used to advance the solution, we do what is known as local extrapolation. The Butcher tableau for a Runge-Kutta pair is written as

c₁ a₁₁ a₁₂ . . . a_1s

... ... ... ...

c_s a_s1 a_s2 . . . a_ss

b₁ b₂ . . . b_s

b^∗₁ b^∗₂ . . . b^∗_s ,

wherebi refers to the coefficients of the advancing method andb^∗_i refers to the coef- ficients of the error measuring method. The local error (denoted byle) can then be estimated by

le=y_n−y^∗_n=∆t

∑

(b_i−b^∗_i)f(Y_i). (1.10) A special case for this estimate occurs when both the advancing method and the error estimating method are stiffly accurate. The estimate for the local error (1.10) simply reduce to

le=y_n−y^∗_n=∆t

∑

(a_s,i−a_s−1,i)f(Y_i) =Y_s−Y_s−1.

1.7 Adaptive Step Size Selection

This section will consider adaptivity in the framework of adaptive Runge-Kutta pairs.

Two different step size selectors will be discussed.

Consider a Runge-Kutta method with advancing method of orderpapplied to a problem with step size∆t. Asymptotic theory sates that the local error in that step is approxi-

mately equal to

le_n≈ϕn·∆t_n^p (1.11)

whereϕnis some unknown quantity. It is then common to make the prediction¹that ϕˆn=ϕn−1,

that is, the distribution inϕnis varying slowly over the course of the time stepping. If we now want to select a new step size corresponding to a user specified tolerancetol, we get the same relation,

tol≈ϕn+1·∆t_n+1^p . (1.12)

Since we assume thatϕnis constant, we can divide (1.12) by (1.11) and get tol

le_n ≈ (∆t_n+1

∆t_n )p

⇒∆t_n+1≈∆t_n· (tol

le_n )1/p

. (1.13)

To compensate for this oversimplified model, it is common to include what is called a pessimistic factor(denotedP) in (1.13) in the following way

∆t_n+1=P·∆t_n (tol

le_n )1/p

, P∈[0,1]. (1.14)

This factor, which can be viewed as a measure on how trustworthy we consider the step size selector to be, must be adjusted according to experimental results on a particular problem. It is common to start withP=0.8, but if the program still rejects/accepts a lot of steps it should be decreased/increased accordingly. A nice feature by selecting step sizes in this way, is that the step size is automatically increased/decreased if the local error is smaller/greater than the specified tolerance.

There are however cases where the assumptionϕn≈ϕn−1can be said to be a poor ap- proximation. In his article, Kjell Gustafsson [Gus94], who has studied control-theoretic techniques for step size selection, lists the following examples where the assumption may fail:

• The properties of the differential equation change along the solution.

• The step size is nonzero during the integration and is not necessarily the lowest- order term that dominates in the error expression. Consequently the error may behave as ifpis larger than expected in (1.11).

1Note that predicted values are denoted with a hat-symbol, so that ˆϕmeans the predicted value forϕ^.

• Some implicit methods lose convergence order when applied to stiff problems (order reduction), causingpin (1.11) to be smaller than expected.

To capture the variation inϕ, Gustafsson uses a model which assumes a linear trend in logϕ^:

log ˆϕn=logϕn−1+∇logϕn−1, ∇logϕn−1=logϕn−1−logϕn−2. The expression for log ˆϕncan be rewritten in the following way

log ˆϕn= (2−q⁻¹)logϕn−1

whereqis the forward shift operator². By taking the logarithm of (1.11) and considering a step size satisfyingle_n=tol, we get

log∆t_n=1

p(logtol−log ˆϕn) By inserting the expression for log ˆϕnwe get

log∆t_n= 1

p(2−q⁻¹)(logtol−logϕn−1), (1.15) where we have used thattol is constant (and thus unaffected by the shift operator).

Finally we have that

logϕn−1=logle_n−1−plog∆t_n−1, which inserted in (1.15) gives

log∆t_n=1

p(2−q⁻¹)(logtol−logle_n−1+plog∆t_n−1)

By applying the shift operators and multiplying out the parenthesis, we arrive at log∆t_n= 1

p(logtol−2 logle_n−1+logle_n−2) +2 log∆t_n−2−log∆t_n−1, which, by removing the logarithms, transforms into

∆t_n=∆t_n−1

∆t_n−2 (

tol·le_n−2 le²_n−1

)1/p

∆t_n−1,

2The forward shift operator is defined asq(ϕn−1) =ϕn.

which is a step size controller based on the two most recent values ofleand∆t. Before describing the actual implementation of the step size controller, we first have to consider what to do when we get rejected steps. If a step is rejected, we have to restart the step size selector since it requires information from two consecutive successful steps. This is done by using the standard asymptotic step size selector (1.14). In the case of two or more consecutive rejects, it is possible that the value ofpis lower than expected (so called order reduction). When this occurs, we can exploit the fact that we have several measurements in the same timestep, which satisfyϕn=ϕn−1, to form the following prediction forp:

le_n le_n−1=

( ∆t_n

∆t_n−1 )_p

⇒pˆ= log(le_n/le_n−1) log(∆t_n/∆t_n−1).

For robustness, Gustafsson suggests that predictions outside the range[0.1,p]should be rejected. We now have all the components needed to present an outline of the complete step size algorithm which can be found in Algorithm 1.

Algorithm 1Gustafsson step size selector

1: ifcurrent step is acceptedthen

2: iffirst steporfist after consecutive rejectsorrestricted timestepthen

3: ∆t_n+1←(

tol len

)1/p

∆t_n .Asymptotic step size selector to restart.

4: else

5: ∆t_n+1←_∆t^∆t_n−1ⁿ (_tol·le

n−1

le²_n

)∆t_n .Gustafsson step size selector.

6: end if

7: else

8: ifconsecutive rejectsthen

9: pˆ←_log(^log(le_∆t^n/len−1)

n/∆t_n−1) 10: pˆ←Restrict(p)ˆ

11: ∆t_n+1←(

tol len

)1/pˆ

∆t_n

12: else

13: ∆t_n+1←(

tol len

)1/p

∆t_n

14: end if

15: end if

16: ∆t_n+1←Restrict(∆t_n+1)

TheRestrict-method should be designed to cap predictions of ˆpto the interval[0.1,p]

and to make sure that the timestep ∆t_n+1 is not increased/decreased too much. If a

timestep is restricted, the algorithm should regard this as a restart.

1.8 Differential Algebraic Equations

Consider the following system of implicit differential equations:

F (d

dty,y,t )

=F(y⁰,y,t) =0. (1.16)

If the Jacobian^∂_∂y^F₀ is nonsingular, we can, by the implicit function theorem, solve (1.16) fory⁰and get a system of ordinary differential equations on the form

y⁰=F(y,t).

If this is not the case, the functionywill have to satisfy some algebraic constraints and we have what is called adifferential algebraic equation(DAE). In order to measure how far a DAE system is from being an ODE system, we can assign the system anindex. It is however not a trivial task to come up with a single definition on what this index should be, different definitions are suitable for different problems. For our purpose, we will consider the definition given by Brenan, Campbell and Petzold [BCP89, Chapter 2.2]

known as thedifferentiation index:

Definition 1 The minimum number of times that all or the part of the system F(y⁰,y,t) =0

must be differentiated with respect to t in order to determine y⁰as a continuous function of y,t, is the differentiation index of F(y⁰,y,t).

When solving a DAE system, selecting initial conditions is not quite as straightforward as when solving an ODE system. For a well posed system of ODEs, a set of initial conditions uniquely determines a solution. For a general high index DAE, on the other hand, finding a set of consistent initial conditions may be very hard due to (potentially numerous) complicated algebraic constraints.

In her article, Kværnø [Kvæ04] states that the ESDIRK methods (presented in section 1.5 of this document) can be directly applied to semi-explicit systems of index 1. This will amount to the following system,

y⁰=f(y,z),

0=g(y,z), (1.17)

where we assumeg_zto be nonsingular. We now show that this system indeed has index 1 according to Definition 1. Differentiating the second component of (1.17), with respect tot, yields

∂^g(y,z)

∂^{t =gyy0(y,z) +gzz0(y,z) =gyf(y,z) +gzz0(y,z).}

We can now solve forz⁰and formulate the following pure ODE problem:

y⁰=f(y,z),

z⁰=−g⁻¹_z (g_yf)(x,y).

Since we needed to carry out one differentiation, we say that this system has index 1. By the same definition, a pure ODE system has index 0.

When solving a semi-explicit system on the form (1.17), we could, sinceg⁻¹_z is assumed to be nonsingular, solveg(y,z) for zand insert into the first equation to get the pure ODE-system

y⁰=f(y,G(y)), z=G(y).

This is mathematically equivalent to applying a Runge-Kutta method in the following way:

Y_i=y_n+∆t

∑

a_{i j}f(Y_i,Z_i), 0=g(Y_i,Z_i), y_n+1=y_n+∆t

∑

b_if(Y_i,Z_i), 0=g(y_n+1,z_n+1).

fori=1, . . . ,s. It is possible to avoid the step of solving forz_n+1by applying a stiffly accurate Runge-Kutta method for which

y_n+1=Y_s⇒z_n+1=Z_s.

1.9 Applying a Finite Element Method

We will now show how ESDIRK methods can be extended to solve PDEs. Our strategy will be to use the ESDIRK method to handle time dependencies, and a finite element method to handle the spatial dependencies of our problem.

Say that we want to solve a PDE given as

u(x,y)−∇²u(x,y)−f(u(x,y)) =0, (x,y)∈Ω, (1.18) u(x,y) =0, (x,y)∈∂Ω,

where the functionf is assumed to be nonlinear. This way of formulating a PDE problem is known as a strong formulation, where we require that our solutionusatisfies (1.18) in every point(x,y), and that it belongs to the space

H₀²(Ω) = {

v:

∫

Ωv²dΩ<∞,

∫

Ω|∇v²|dΩ<∞,

∫

Ω|∇²v²|dΩ<∞, v|_∂Ω=0 }

. We can relax the regularity conditions by multiplying (1.18) by a functionvcoming from some function space ˆVto get

uv−∇²uv−f(u)v=0.

If we now integrate this equation over the domainΩ, and do integration by parts on the diffusion term, we end up with

∫

ΩuvdΩ+

∫

Ω∇u∇vdΩ−^∫

Ωf(u)vdΩ=0. (1.19) We now defineu∈V to be a weak solution of (1.18) if it satisfies (1.19) for allv∈Vˆ. This is imposed for allv∈Vˆ since our choice ofvis completely arbitrary. The space V is known as the trial space and is the space where the weak solution is located, while the space ˆV is known as the test space. The functions u and v are thus referred to as trial and test functions. In our example the two spaces coincides with the space H₀¹(Ω)due to the homogeneous Dirichlet boundary conditions. In general, the spaceV is the space of functions with the required regularity which also satisfies the Dirichlet boundary conditions of our problem. A finite element method seeks to approximate the weak solution of (1.18) by constructing the finite dimensional subspacesV_h⊂V and Vˆ_h⊂Vˆ and use those to construct an approximation to (1.19). In our example, the two spaces are the same and will be referred to as justV_h. Let nowV_hbe defined as

V_h=span{ϕ1,ϕ2, . . . ,ϕm},

where all the functionsϕiare linearly independent functions selected fromV. Which functions we choose to span this subspace depends on which kind of element we want to use. Our numerical approximation to the solutionu, denoted byu_h, can can be written as a linear combination of the functions spanningV_has such

u_h=

∑

U_iϕi, ϕi∈V_h,

whereU= [U₁,U₂, . . . ,U_m]are the degrees of freedom we need to solve for. If we insert this approximation into (1.19) we get

∑

( U_i

∫

Ω[ϕiϕˆj+∇ϕi·∇ϕˆj]dΩ )

−^∫

Ωf (m

i=1

∑

U_iϕi

)

ϕˆjdΩ=0, j=1, . . . ,m, which is a system of nonlinear equations inUsince f is assumed to be nonlinear. By solving this system, we can construct our approximationu_h=∑^mi=1U_iϕi.

We now show how applying a finite element method to a PDE can give rise to a system of DAEs. Consider the Cahn-Hilliard system (which will be studied in some detail in the experiments chapter) withw=w(x,y)andz=z(x,y):

dw dt =∇²z,

0=z−f(w)−∇²w.

The weak formulation of this problem amounts to findingw×z∈V×Vsuch that

∫ Ω

dw

dt vdΩ=−^∫

Ω∇z∇vdΩ, 0=

∫

Ω(z−f(w))qdΩ+

∫

Ω∇w∇qdΩ,

for allv×q∈Vˆ×Vˆ whereV,Vˆ are appropriate finite element spaces. By applying a finite element method we seek the approximationsw_h×z_hsuch that

∫ Ω

dw_h

dt vdΩ=−

∫

Ω∇z_h∇vdΩ, 0=

∫

Ω(z_h−f(w_h))qdΩ+

∫

Ω∇w_h∇qdΩ, for allv×q∈Vˆ_h×Vˆ_hwhere ˆV_h⊂V_handV_h⊂V. By inserting

w_h=

∑

W_iϕi, z_h=

∑

Z_iϕi, q=v=ϕˆ ∈Vˆ, we get the following system of equations:

∑

∫ Ω

dW_i

dt ϕiϕˆjdΩ=

∑

Z_i

∫

Ω∇ϕi·∇ϕˆjdΩ, 0=

∑

( Z_i

∫

ΩϕiϕˆjdΩ )

+

∫ Ωf

(m i=1

∑

W_iϕi

)

ϕˆjdΩ+

∑

( W_i

∫

Ω∇ϕi·∇ϕˆjdΩ )

, for j=1, . . . ,m. This is now a DAE system of index 1 (withW_i as ODE-components andZ_i as algebraic components) which can be solved directly by applying one of the ESDIRK methods described in section 1.5.

Part II

Implementation

2.1 Introduction

This part aims to give a presentation of my work resulting in the Gryphon module. First, a short introduction to FEniCS with special focus on the components relevant for the development of Gryphon, will be given. Next, the need for a module like Gryphon will be motivated through an example. The focus will be on how the work flow of solving a time dependent PDE with FEniCS compares to the work flow of solving the same problem with FEniCS and Gryphon. This part will not go into details regarding the source code, instead we will focus on how I wanted Gryphon to work and how I used object orientation to achieve a clean program architecture. For details on the source code, the reader is encouraged to consult the attached documentation (see appendix A for details).

2.2 The FEniCS Project

Before discussing implementation details concerning Gryphon, we will give a short in- troduction to underlying framework, namely the FEniCS project. A brief summary of the project can be found in the "about" section on the FEniCS web page³ where it is stated that:

The FEniCS Project is a collaborative project for the development of inno- vative concepts and tools for automated scientific computing, with a partic- ular focus on automated solution of differential equations by finite element methods.

We will now dive into some key features in FEniCS in order to build some termi- nology which will be useful later. For a more complete coverage of the topics pre- sented, the reader is encouraged to download the FEniCS book [LMW11] fromhttp:

//fenicsproject.org/book/.

UFL (Unified Form Language) UFL is the language used to specify problems in FEniCS. It is described in the following way by its authors on the UFL web page⁴:

The Unified Form Language (UFL) is a domain specific language for decla- ration of finite element discretizations of variational forms. More precisely,

3http://fenicsproject.org/about

4http://launchpad.net/ufl

it defines a flexible interface for choosing finite element spaces and defining expressions for weak forms in a notation close to mathematical notation.

The elegance of this language is best shown with an example. Say that we would like to solve the Laplace equation given as

∇²u=0

on some domainΩ. The weak form of this equation amounts to findingu∈V such that a(u,v) =0,

−^∫

Ω∇u∇vdΩ=0,

for allv∈Vˆ for appropriate trial/test spacesV,Vˆ. This problem can be described by the following UFL code where we have used first order Lagrange elements as basis functions on a 2D domain:

element = FiniteElement("Lagrange", triangle , 1) v =TestFunction(element)

u =TrialFunction(element) a =−inner(grad(u)∗grad(v))∗dx

Listing 2.1: UFL for code defining the Laplace equation.

As this code snippet shows, UFL indeed provides a very clean and intuitive way of defin- ing weak forms by using familiar names likeinnerfor inner product andgradfor the gradient operator. UFL also supports symbolic differentiation which is a tremendous ad- vantage when working with nonlinear PDEs, since we are no longer required to estimate the Jacobian matrix used when doing Newton iterations.

FFC (The FEniCS Form Compiler) and UFC (Unified Form-assembly Code) The FEniCS Form Compiler is responsible for interpreting UFL and translate it into UFC.

Simply put, UFC is C++ code responsible for assembling the systems of equations de- scribed in UFL. The way FFC is used depends on which FEniCS API you are working with.

If you are working with the C++-API, the common approach is to define your problem in UFL, send the UFL-file to FFC to get a UFC-object, and then include the header file of the UFC-object in your C++ program.

If you are working with the Python-API, the UFL/UFC handling is more seamless. UFL can be imported as a module in Python, allowing users to create and manipulate UFL

forms while the script is running. The UFC generation can be postponed until it is specif- ically requested by a FEniCS component, like for instance a linear/nonlinear solver. This is made possible by a module calledInstant, which make FFC support just-in-time (JIT) compilation. Instant also supports caching so that already generated UFC code can be reused if possible.

For details on the concepts presented above, the reader is encouraged to consult [LMW11]

chapter 11 for FFC, chapter 14 for Instant, chapter 16 for UFC and chapter 17 for UFL.

2.3 The Gryphon Module

We will start by motivating why a module like Gryphon is useful when working with time dependent PDEs. This will be done by showing how the work flow for solving a time dependent in FEniCS, without Gryphon, can be improved upon. As working example we will consider the Heat equation with a source term given as

∂^u

∂^{t =∇2u+ (}β−2−2α^{), u}∈Ω,

u=1+x²+y²+β^{t, u}∈∂Ω, (2.20) u=1+x²+y², t=0,

on the domainΩ= (x,y)∈[0,1]×[0,1]withα^{=3 and}β^{=1.2 withu=u(x,y). This} example was found in the FEniCS tutorial in the section "time-dependent problems"

[Pro12b]. Because of the time derivative, we can not solve the problem in FEniCS in its current form, we have to apply some sort of time integrator. In the FEniCS tutorial, the backward Euler method is applied, giving rise to the new problem:

u_n−u_n−1

∆t =∇²u_n+ (β−2−2α^{), u}n∈Ω, u_n=1+x²+y²+β^{t, u}n∈∂Ω, u₀=1+x²+y²,

where∆t is the step size in time andu_nrefers to the numerical solution in time stepn.

Rearranging the first expression yields

u_n−∆t∇²u_n=u_n−1+∆t(β−2−2α^).

The linear variational problem corresponding to each time step amounts to findingu_n∈V such that

a(u_n,v) =`(v),

∫

Ω(u_nv+∆t∇u_n∇v)dΩ=

∫

Ω(u_n−1+∆t(β−2−2α^))dΩ, for allv∈Vˆ whereV,Vˆ are appropriate test/trial spaces.

The problem is now in a form suitable for FEniCS, but we still have to write additional code for doing the actual time stepping process. In the simplest case, this can be achieved by a while-loop which solves the variational problem for each time step and updates the solutionu_n+1←u_n. If we want to utilize more advanced tools for the time stepping, like for instance adaptive step size control, the required code becomes considerably more involved.

Additional complications arise if we want to use a more advanced time integrator than the backward Euler method. Reformulating your initial problem to a form suitable for FEniCS, and then writing the code for solving each time step, may not be straightfor- ward. If you are experimenting with different equations this work flow will most likely lead to a lot of boilerplate code which is both tedious and hard to maintain.

In the end we have distanced ourselves from the actual problem we were trying to solve in the first place. Gryphon aims to bridge this gap by the following philosophy:

Just as FEniCS provides an easy and intuitive way of defining and solving differential equations, Gryphon should provide an easy and intuitive way of applying and customizing a time integrator to solve a time dependent PDE.

In other words, I wanted Gryphon to be a tool which captured the feel of working in a FEniCS environment, where the main focus is on the problem you are trying to solve.

This goal gave rise to the work flow presented in figure 2.2, where the already described work flow for solving a time dependent PDE in FEniCS is compared to solving the same problem in FEniCS with Gryphon. Gryphon adapts the already existing parameter sys- tem from FEniCS, allowing users to familiarize themselves with the available parameters by calling

info (Gryphon_object.parameters,verbose=True)

All the available parameters are explained in the user manual found in part IV.

.

. Initial value problem (strong formulation)

. Apply time integrator.

. Derive weak form.

. Formulate the modified problem in FEniCS.

. Write code for do- ing the time stepping.

. Get solution.

. Derive weak form.

. Formulate the initial problem in FEniCS.

. Create Gryphon object.

Set desired parameters.

. Call.solve on Gryphon object.

. Get solution.

.Optional output Run time statistics

Plot of selected/rejected step sizes

Plot of solution in each time step

Figure 2.2: Work flow for solving a time dependent PDE in FEniCS without Gryphon (red) and with Gryphon (green).

To give a better understanding on how Gryphon works, we will give a short summary of the program architecture. When developing Gryphon, I sought to follow the don’t re- peat yourself (DRY) principle, which is a principle formulated by Andy Hunt and Dave Thomas stating that"Every piece of knowledge must have a single, unambiguous, au- thoritative representation within a system"[HT99]. This gave rise to the class hierarchy found in figure 2.3. To elaborate on the reason behind this structure, each class will be presented with its intended usage.

gryphon_base: This class is the highest superclass in the Gryphon hierarchy. Its con- structor is responsible for assigning a variety of class variables to be used in both the gryphon_toolbox-class and thetime_integrator-class. It is also contains methods for input verification, printing of program progress and error handling. In short, it is designed to be a platform for which to build tools for doing time integration, and for the time integrators themselves.

. .

gryphon_base

. gryphon_toolbox

. time_integrator

Figure 2.3: Class hierarchy in Gryphon gryphon_toolbox: This class is in-

tended to contain tools found useful when implementing a time integrator. This res- onates well with the DRY principle in the sense that we collect tools relevant for several time integrators in one place, which makes for easy maintenance. No- table tools in this class includes step size selectors and methods for output gener- ation (run time statistics, plot of accept- ed/rejected step sizes).

time_integrator: This final class layer is intended to contain the realization of some time integrator. Each class in this layer must contain the code for augment- ing a user specified UFL form with the

code amounting to applying the desired time integrator. As of today, the only class of methods realized in Gryphon are ESDIRK methods. To show how other time integrators can be implemented, we will take a closer look at how the ESDIRK class is implemented, with the intent that the methodology can be copied.

2.3.1 The ESDIRK class

The ESDIRK class represents the realization of the ESDIRK methods developed by Anne Kværnø presented in section 1.5. It is capable of solving systems of PDEs which

either semidiscretize into a pure ODE system or a DAE system of index 1. The construc- tor expects to receive UFL forms representing the right hand side of a system of PDEs on the form

Mu⁰=rhs whereMmay be singular.

Constructor: The constructor is responsible for passing on the user given data to the gryphon_toolboxclass which then passes the data to the gryphon_baseclass for input verification and initialization.

getLinearVariationalForms: This method is responsible for creating the UFL-forms arising from applying an ESDIRK method to the user given right hand side. This method assumes the problem to be linear and creates a variational problem on the form: Find u_n∈V such that

a(u_n,v) =`(v)

for allv∈Vˆ whereais bilinear inu_nandvand`is linear inv.

getNonlinearVariationalForms: This method is responsible for creating the UFL- forms arising from applying an ESDIRK to the user given right hand side. This method assumes the problem to be nonlinear and creates a variational problem on the form: Find u_n∈V such that

F(u_n;v) =0 for allv∈Vˆ whereFis linear inv.

solve: This method is responsible for initializing and performing the time stepping loop. The process is outlined in algorithm 2, but a more detailed description can be found by inspecting the attached documentation.

Algorithm 2ESDIRK.solve

1: Get Butcher-tableau and createsFunction-objects to store the stage values

2: iflinear problemthen

3: CallgetLinearVariationalForms

4: else

5: CallgetNonlinearVariationalForms

6: end if

7: whilein the time stepping loopdo

8: Set first stage value equal to previous time step

9: Solve for the nexts−1 stage values

10: Estimate local error

11: ifLocal error in current time step<tolerancethen

12: Accept step, calculate new step size.

13: else

14: Reject step, calculate new step size.

15: end if

16: end while

17: Generate specified output and terminate program

2.4 A Code Example

We will round off this part by showing the code required to solve the problem introduced in section 2.3:

∂^u

∂^{t =∇2u+ (}β−2−2α^{), u}∈Ω, u=1+x²+y²+β^{t, u}∈∂Ω, u=1+x²+y², t=0,

on the domainΩ= (x,y)∈[0,1]×[0,1]withα^{=3 and}β^{=1.2 withu=u(x,}y). Listing 2.2 shows the code for solving (2.20) using FEniCS alone. The code was constructed by following the FEniCS tutorial. Listing 2.3 shows the code for solving (2.20) using FEniCS in combination with Gryphon. The code snippets can also be found in the attached filessol_FEniCS.py and sol_FEniCS_Gryphon.py. See appendix A for details.