4.11 Time Integration

This directory consists of the following files:

The original time-stepping classes (found in TimeIntegrationScheme) where implemented by Vos et al. [62] and a detailed explanation of the mathematics and computer science concepts are provided there. After the original implementation, Cantwell et al. extended Vos’ original work by adding the time-stepping routines into a factory pattern (found in TimeIntegrationWrapper).

4.11.1 General Linear Methods

General linear methods (GLMs) can be considered as the generalization of a broad range of different numerical methods for ordinary differential equations. They were introduced by Butcher and they provide a unified formulation for traditional methods such as the Runge-Kutta methods and the linear multi-step methods such as Adams-Bashforth. From an implementation point of view, this means that all these numerical methods can be abstracted in a similar way. As this allows a high level of generality, it is chosen in Nektar++ to cast all time integration schemes in the framework of general linear methods.

For background information about general linear methods, please consult [14].

4.11.2 Introduction

The standard initial value problem can written in the form

= f(t,y), y(t0) = y0

where y is a vector containing the variable (or an array of array containing the variables).

In the formulation of general linear methods, it is more convenient to consider the ODE in autonomous form, i.e.

= ^ f (^y ), ^ y (t0) = ^ y 0.

4.11.3 Formulation

Suppose the governing differential equation is given in autonomous form, the n^th step of the GLM comprising

• r steps (as in a multi-step method)
• s stages (as in a Runge-Kutta method)

is formulated as:

where Y _i are referred to as the stage values and F_j as the stage derivatives. Both quantities are related by the differential equation:

The matrices A = [a_ij], U = [u_ij], B = [b_ij], V = [v_ij] are characteristic of a specific method. Each scheme can then be uniquely defined by a partitioned (s + r) × (s + r) matrix

4.11.4 Matrix notation

Adopting the notation:

^y [n-1] = , ^ y [n] = , Y = , F =

the general linear method can be written more compactly in the following form:

=

where I_N is the identity matrix of dimension N × N.

4.11.5 General Linear Methods in Nektar++

Although the GLM is essentially presented for ODE’s in its autonomous form, in Nektar++ it will be used to solve ODE’s formulated in non-autonomous form. Given the ODE,

= f(t,y), y(t0) = y0

a single step of GLM can then be evaluated in the following way:

• input

  y^[n-1], i.e. the r subvectors comprising y_i^[n-1] - t^[n-1], i.e. the equivalent of y^[n-1] for the time variable t.

• step 1: The stage values Y _i, T_i and the stage derivatives F_i are calculated through the relations:
  1. Y _i = Δt∑ _j=0^s-1a_ijF_j + ∑ _j=0^r-1u_ijy_j^[n-1], i = 0,1,…,s - 1
  2. T_i  =  Δt∑ _j=0^s-1a_ij + ∑ _j=0^r-1u_ijt_j^[n-1], i = 0,1,…,s - 1
  3. F_i  =  f(T_i,Y _i), i = 0,1,…,s - 1
• step 2: The approximation at the new time level y^[n] is calculated as a linear combination of the stage derivatives F_i and the input vector y^[n-1]. In addition, the time vector t^[n] is also updated
  1. y_i^[n] = Δt∑ _j=0^s-1b_ijF_j + ∑ _j=0^r-1v_ijy_j^[n-1], i = 0,1,…,r - 1
  2. t_i^[n]  =  Δt∑ _j=0^s-1b_ij + ∑ _j=0^r-1v_ijt_j^[n-1], i = 0,1,…,r - 1
• output
  1. y^[n], i.e. the r subvectors comprising y_i^[n]. y₀^[n] corresponds to the actual approximation at the new time level.
  2. t^[n] where t₀^[n] is equal to the new time level t + Δt.

For a detailed description of the formulation and a deeper insight of the numerical method see [62].

4.11.6 Types of time Integration Schemes

Nektar++ contains various classes and methods which implement the concept of GLMs. This toolbox is capable of numerically solving the generalised ODE using a broad range of different time-stepping methods. We distinguish the following types of general linear methods:

• Formally Explicit Methods: These types of methods are considered explicit from an ODE point of view. They are characterised by a lower triangular coefficient matrix A, ”i.e.” a_ij = 0 for j ≥ i. To avoid confusion, we make a further distinction:
  • direct explicit method: Only forward operators are required.
  • indirect explicit method: The inverse operator is required.
• Diagonally Implicit Methods: Compared to explicit methods, the coefficient matrix A has now non-zero entries on the diagonal. This means that each stage value depend on the stage derivative at the same stage, requiring an implicit step. However, the calculation of the different stage values is still uncoupled. Best known are the DIRK schemes.
• IMEX schemes: These schemes support the concept of being able to split right hand forcing term into an explicit and implicit component. This is useful in advection diffusion type problems where the advection is handled explicitly and the diffusion is handled implicit.
• Fully Implicit Methods Methods: The coefficient matrix has a non-zero upper triangular part. The calculation of all stages values is fully coupled.

The aim in Nektar++ is to fully support the first three types of GLMs. Fully implicit methods are currently not implemented.

4.11.7 Usage

The goal of abstracting the concept of general linear methods is to provide users with a single interface for time-stepping, independent of the chosen method. The classes tree allows the user to numerically integrate ODE’s using high-order complex schemes, as if it was done using the Forward-Euler method. Switching between time-stepping schemes is as easy as changing a parameter in an input file.

In the already implemented solvers the time-integration schemes have been set up according to the nature of the equations. For example the incompressible Navier-Stokes equations solver allows the use of three different Implicit-Explicit time-schemes if solving the equations with a splitting-scheme. This is because this kind of scheme has an explicit and an implicit operator that combined solve the ODE’s system.

Once aware of the problem’s nature and implementation, the user can easily switch between some (depending on the problem) of the following time-integration schemes:

Nektar++ input file for your problem will ask you just the string corresponding the time-stepping scheme you want to use (between quotation marks in the previous list), and few parameters to define your integration in time (time-step and number of steps or final time). For example:

4.11.8 Implementation of a time-dependent problem

In order to implement a new solver which takes advantage of the time-integration class in Nektar++, two main ingredients are required:

• A main files in which the time-integration of you ODE’s system is initialized and performed.
• A class representing the spatial discretization of your problem, which reduces your system of PDE’s to a system of ODE’s.

Your pseudo-main file, where you actually loop over the time steps, will look like

We can distinguish three different sections in the code above

Definitions

In this section you define the basic parameters (like time-step, initial time, etc.) and the time-integration objects. The operators are not all required, it depends on the nature of your problem and on the type of time integration schemes you want to use. In this case, the problem has been set up to work just with Forward-Euler, then for sure you will not need the implicit operator. An object named equation has been initialized, is an object of type YourClass, where your spatial discretization and the functions which actually represent your operators are implemented. An example of this class will be shown later in this page.

Initialisations

The second part consists in the scheme initialization. In this example we set up just Forward-Euler, but we can set up more then one time-integration scheme and quickly switch between them from the input file. Forward-Euler does not require any other scheme for the start-up procedure. High order multi-step schemes may need lower-order schemes for the start up.

Integration

The last step is the typical time-loop, where you iterate in time to get your new solution at each time-level. The solution at time tⁿ⁺¹ is stored into vector U (you need to properly initialize this vector). U is an Array of Arrays, where the first dimension corresponds to the number of variables (eg. u,v,w) and the second dimension corresponds to the variables size (e.g. the number of modes or the number of physical points).

The variable ODE is an object which contains the methods. A class representing a PDE equation (or a system of equations) must have a series of functions representing the implicit/explicit part of the method, which represents the reduction of the PDE’s to a system of ODE’s. The spatial discretization and the definition of this method should be implemented in YourClass. &YourClass::YourExplicitOperatorFunction is a functor, i.e. a pointer to a function where the method is implemented. equation is a pointer to the object, i.e. the class, where the function/method is implemented. Here a pseudo-example of the .h file of your hypothetical class representing the set of equations. The implementation of the functions is meant to be in the related .cpp file.

4.11.9 Strongly imposed essential boundary conditions

Dirichlet boundary conditions can be strongly imposed by lifting the known Dirichlet solution. This is equivalent to decompose the approximate solution y into an known part, y^D, which satisfies the Dirichlet boundary conditions, and an unknown part, y^H, which is zero on the Dirichlet boundaries, i.e.

In a Finite Element discretisation, this corresponds to splitting the solution vector of coefficients y into the known Dirichlet degrees of freedom y^D and the unknown homogeneous degrees of freedom y^H. Ordering the known coefficients first, this corresponds to:

y =

The generalised formulation of the general linear method (i.e. the introduction of a left hand side operator) allows for an easier treatment of these types of boundary conditions. To better appreciate this, consider the equation for the stage values for an explicit general linear method where both the left and right hand side operator are linear operators, i.e. they can be represented by a matrix.

In case of a lifted known solution, this can be written as:

= Δt∑ j=0i-1a ij + ∑ j=0r-1u ij,
i = 0,1,…,s - 1

In order to calculate the stage values correctly, the explicit operator should now be implemented to do the following:

=

Note that only the homogeneous part b^H will be used to calculate the stage values. This means essentially that only the bottom part of the operation above, i.e. L^HDy^D + L^HHy^H is required. However, sometimes it might be more convenient to use/implement routines for the explicit operator that also calculate b^D.

An implicit method should solve the system:

= =

for the unknown vector y. This can be done in three steps:

• Set the known solution y^D
• Calculate the modified right hand side term b^H-H^HDy^D
• Solve the system below for the unknown y^H, i.e. H^HHy^H = b^H-H^HDy^D

4.11.10 How to add a new time-stepping method

To add a new time integration scheme, follow the steps below:

• Choose a name for the method and add it to the TimeIntegrationMethod enum list.
• Populate the switch statement in the TimeIntegrationScheme constructor with the coefficients of the new method.
• Use ( or modify) the function InitializeScheme to select (or implement) a proper initialisation strategy for the method.
• Add documentation for the method (especially indicating what the auxiliary parameters of the input and output vectors of the multi-step method represent)

4.11.11 Examples of already implemented time stepping schemes

Here we show some examples time-stepping schemes implemented in Nektar++, to give an idea of what is required to add one of them.

Forward Euler

= , y[n] = =

Backward Euler

= , y[n] = =

2^nd order Adams-Bashforth

= , y[n] = =

1^st order IMEX Euler-Backward/ Euler-Forward

= with y[n] = =

2^nd order IMEX Backward Different Formula & Extrapolation

=

with

y[n] = =

3rdorder IMEX Backward Different Formula & Extrapolation

=

with

y[n] = =

2^nd order Adams-Moulton

= , y[n] = =

Midpoint Method

= , y[n] = =

RK4: the standard fourth-order Runge-Kutta scheme

= , y[n] = =

2^nd order Diagonally Implicit Runge-Kutta (DIRK)

= with λ = , y[n] = =

3^rd order Diagonally Implicit Runge-Kutta (DIRK)

=

with

λ = 0.4358665215, y[n] = =

3^rd order L-stable, three stage IMEX DIRK(3,4,3)

	=

with

λ = 0.4358665215, y[n] = =