General Linear Methods

The implementation of time integration schemes in Nektar++

Introduction

General linear methods can be considered as the generalisation 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. From an implementational 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 the following references:
[1] Butcher, J.C. (2006) General linear methods, Acta Numerica 15, 157-256
[2] http://www.math.auckland.ac.nz/~butcher/conferences.html

General linear methods

The standard initial value problem can written in the form

\[ \frac{d\boldsymbol{y}}{dt}=\boldsymbol{f}(t,\boldsymbol{y}),\quad \boldsymbol{y}(t_0)=\boldsymbol{y}_0 \]

where $\boldsymbol{y}$ is a vector containing the variable (or an array of array contianing the variables).
In the formulation of general linear methods, it is more convenient to consider the ODE in autonomous form, i.e.

\[ \frac{d\boldsymbol{\hat{y}}}{dt}=\boldsymbol{\hat{f}}(\boldsymbol{\hat{y}}),\quad \boldsymbol{\hat{y}}(t_0)= \boldsymbol{\hat{y}}_0. \]

Formulation

Suppose the governing differential equation is given in autonomous form, the $n^{th}$ step of the general linear method comprising

is formulated as:

\begin{eqnarray*} \boldsymbol{Y}_i = \Delta t\sum_{j=0}^{s-1}a_{ij}\boldsymbol{F}_j+\sum_{j=0}^{r-1}u_{ij} \boldsymbol{\hat{y}}_{j}^{[n-1]}, \qquad i=0,1,\ldots,s-1 \\ \boldsymbol{\hat{y}}_{i}^{[n]}=\Delta t\sum_{j=0}^{s-1}b_{ij}\boldsymbol{F}_j+\sum_{j=0}^{r-1}v_{ij} \boldsymbol{\hat{y}}_{j}^{[n-1]}, \qquad i=0,1,\ldots,r-1 \end{eqnarray*}

where $\boldsymbol{Y}_i$ are referred to as the stage values and $\boldsymbol{F}_i$ as the stage derivatives. Both quantities are related by the differential equation:

\begin{displaymath} \boldsymbol{F}_i = \boldsymbol{\hat{f}}(\boldsymbol{Y}_i). \end{displaymath}

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 partioned $(s+r)\times(s+r)$ matrix

\begin{displaymath} \left[ \begin{array}{cc} A & U\\ B & V \end{array}\right] \end{displaymath}

Matrix notation

Adopting the notation:

\begin{displaymath} \boldsymbol{\hat{y}}^{[n-1]}= \left[\begin{array}{c} \boldsymbol{\hat{y}}^{[n-1]}_0\\ \boldsymbol{\hat{y}}^{[n-1]}_1\\ \vdots\\ \boldsymbol{\hat{y}}^{[n-1]}_{r-1} \end{array}\right],\quad \boldsymbol{\hat{y}}^{[n]}= \left[\begin{array}{c} \boldsymbol{\hat{y}}^{[n]}_0\\ \boldsymbol{\hat{y}}^{[n]}_1\\ \vdots\\ \boldsymbol{\hat{y}}^{[n]}_{r-1} \end{array}\right],\quad \boldsymbol{Y}= \left[\begin{array}{c} \boldsymbol{Y}_0\\ \boldsymbol{Y}_1\\ \vdots\\ \boldsymbol{Y}_{s-1} \end{array}\right],\quad \boldsymbol{F}= \left[\begin{array}{c} \boldsymbol{F}_0\\ \boldsymbol{F}_1\\ \vdots\\ \boldsymbol{F}_{s-1} \end{array}\right]\quad \end{displaymath}

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

\begin{displaymath} \left[ \begin{array}{c} \boldsymbol{Y}\\ \boldsymbol{\hat{y}}^{[n]} \end{array}\right] = \left[ \begin{array}{cc} A\otimes I_N & U\otimes I_N \\ B\otimes I_N & V\otimes I_N \end{array}\right] \left[ \begin{array}{c} \Delta t\boldsymbol{F}\\ \boldsymbol{\hat{y}}^{[n-1]} \end{array}\right] \end{displaymath}

where $I_N$ is the identity matrix of dimension $N\times N$.

Evaluation of an General Linear Method

Although the General linear method 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,

\[ \frac{d\boldsymbol{y}}{dt}=\boldsymbol{f}(t,\boldsymbol{y}),\quad \boldsymbol{y}(t_0)=\boldsymbol{y}_0 \]

a single step of General Linear Method can then be evaluated in the following way:

General Linear Methods in Nektar++

In Nektar++, we do not use the standard General Linear Methods formulation. As mentioned above, we will use the non-autonomous form as the explicit treatment of the time variable $t$ allows for more flexibility.

For a detailed describtion of the formulation and a deeper insight of the numerical method see:

Peter E.J. Vos, Claes Eskilsson, Alessandro Bolis, Sehun Chun,Robert M. Kirby & Spencer J. Sherwin, (2011),
A generic framework for time-stepping partial differential equations (PDEs):
general linear methods, object-oriented implementation and application to fluid problems,
International Journal of Computational Fluid Dynamics, 25:3, 107-125

Type of Time Integration Schemes

Nektar++ contains various classes and methods which implement the concept of the General Linear Methods. 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:

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

How to use

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 TimeIntegrationScheme class allow the user to numerically integrate ODE's using high-order complex schemes, as if it were done using the forward euler method. Switching between time-stepping schemes should be as easy as changing a parameter in an input file. The only thing the user should provide, is an implementation of the left and right hand side operation of the generalised ODE to be solved.

To introduce the implementation of time stepping schemes in Nektar++, consider the following example: 

 NekDouble timestepsize = 0.1;
 NekDouble time         = 0.0;
 
 Array<OneD, Array<OneD, NekDouble> >  y;
 
 LibUtilities::TimeIntegrationSchemeOperators ode;

 ode.DefineImplicitSolve (&function1, &object1);
 ode.DefineOdeRhs        (&function2, &object2);


 LibUtilities::TimeIntegrationSchemeKey       IntKey(LibUtilities::eClassicalRungeKutta4);
 LibUtilities::TimeIntegrationSchemeSharedPtr IntScheme = LibUtilities::TimeIntegrationSchemeManager()[IntKey];

 LibUtilities::TimeIntegrationSolutionSharedPtr y_0;

 y_0 = IntScheme->InitializeScheme(timestepsize,y,time,ode);

 for(int n = 0; n < nsteps; ++n)
 {
  y = IntScheme->TimeIntegrate(timestep,y_0,ode);
  time += timestepsize;
 }


We can distinguish four different sections in the code above:

Required Methods

An object bar of the class Foo should be passed to the TimeIntegrationScheme class routines. This class Foo should have the following public member functions:

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 lifted function, $y^{\mathcal{D}}$, which satisfies the Dirichlet boundary conditions, and an unknown homogeneous function, $y^{\mathcal{D}}$, which is zero on the Dirichlet boundaries, i.e.

\[ y = y^{\mathcal{D}} + y^{\mathcal{H}} \]

In a Finite Element discretisation, this corresponds to splitting the solution vector of coefficients $\boldsymbol{y}$ into the known Dirichlet degrees of freedom $\boldsymbol{y}^{\mathcal{D}}$ and the unknown homogeneous degrees of freedom $\boldsymbol{y}^{\mathcal{H}}$. If ordering the known coefficients first, this corresponds to:

\[ \boldsymbol{y} = \left[ \begin{array}{c} \boldsymbol{y}^{\mathcal{D}} \\ \boldsymbol{y}^{\mathcal{H}} \end{array} \right] \]

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.

\[ \boldsymbol{M}\boldsymbol{Y}_i = \Delta t\sum_{j=0}^{i-1}a_{ij}\boldsymbol{L}\boldsymbol{Y}_j+\sum_{j=0}^{r-1}u_{ij}\boldsymbol{y}_{j}^{[n-1]}, \qquad i=0,1,\ldots,s-1 \]

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

\[ \left[ \begin{array}{cc} \boldsymbol{M}^{\mathcal{DD}} & \boldsymbol{M}^{\mathcal{DH}} \\ \boldsymbol{M}^{\mathcal{HD}} & \boldsymbol{M}^{\mathcal{HH}} \end{array} \right] \left[ \begin{array}{c} \boldsymbol{Y}^{\mathcal{D}}_i \\ \boldsymbol{Y}^{\mathcal{H}}_i \end{array} \right] = \Delta t\sum_{j=0}^{i-1}a_{ij} \left[ \begin{array}{cc} \boldsymbol{L}^{\mathcal{DD}} & \boldsymbol{L}^{\mathcal{DH}} \\ \boldsymbol{L}^{\mathcal{HD}} & \boldsymbol{L}^{\mathcal{HH}} \end{array} \right] \left[ \begin{array}{c} \boldsymbol{Y}^{\mathcal{D}}_j \\ \boldsymbol{Y}^{\mathcal{H}}_j \end{array} \right] +\sum_{j=0}^{r-1}u_{ij} \left[ \begin{array}{c} \boldsymbol{y}^{\mathcal{D}[n-1]}_j \\ \boldsymbol{y}^{\mathcal{H}[n-1]}_j \end{array} \right], \qquad i=0,1,\ldots,s-1 \]

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

Implemented integration schemes

Currently the time integration schemes below are implemented in Nektar++. We will list their coefficients here and we will also indicate what the auxiliary parameters $\boldsymbol{\hat{y}}_{j}^{[n]}$ ( $j=1,\ldots,r-1$) of the multistep methods represent:

How to add a method

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

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Thu Mar 29 2012 17:00:28 for Nektar++ by  doxygen 1.7.1