Kaskade 7 > Time Integration > LIMEX

LIMEX

Introduction

The name LIMEX originates from L inearly IM plicit EX trapolation method. LIMEX is a one-step extrapolation method with linearly implicit Euler as base method for differential-algebraic equations (DAE) combined with an advanced order and stepsize control. LIMEX is also suitable for use in parallel programming. It is used for numerically solving linearly-implicit differential systems of the form $B \dot u = L(u)$ with initial conditions. The matrix $B$ may be singular or non-constant. Extrapolation methods allow the construction of higher order approximations.

LIMEX

Extrapolation method for linearly implicit Euler

The algorithm for one basic integration step is depicted in the following pseudocode from EhrigNowakDeuflhard1996:

$\begin{aligned} & \text{for $j=1,\ldots,j_{max}$ while convergence criterion not met} \\ & \qquad \text{set stepsize $\tau_j$} \\ & \qquad \text{for $k=0,\ldots,j-1$} \\ & \qquad \qquad \delta u^k = (B-\tau_j L’(u_0))^{-1}\tau_j L(u^k), \quad \text{where } \delta u^k = u^{k+1}-u^k \\ & \qquad T_{j,1} = u^j \\ & \qquad \text{if $j>1$ compute $T_{j,j}$ and check convergence} \\ & u^{new} = T_{j,j} \end{aligned}$

With basic stepsize $H$ and definition of stepsizes $\tau_j = \frac{H}{n_j}$ for the harmonic sequence $n_j = 1,2,3, …$ we obtain different numerical solutions $T_{j,1}$ at points $(t_0 + H)$ by applying $n_j$ steps of the linearly implicit Euler method $\delta u^k = (B-\tau_j L’(u^k))^{-1}\tau_j L(u^k)$ . The monotonically increasing sequence of positive integers $n_j$ can be chosen alternatively as Romberg sequence, Bulirsch sequence, etc.

After the approximations $T_{j,1}$ for $u^j = u(t_0+\tau_j n_j)$ are calculated, the Aitken-Neville formula defines the higher order approximations $T_{j,k}$ for $\gamma = 1, 2$ recursively by $T_{j,k} = T_{j,k-1} + \frac{T_{j,k-1}-T_{j-1,k-1}}{(\frac{n_j}{n_{j-k+1}})^\gamma-1}, \qquad k=1,\ldots,j.$ Parameter $\gamma$ occurs in the asymptotic expansion of the error of one basic approximation step as exponent of stepsize h. The extrapolation method is designed in a way to cancel those terms out, thus yielding a higher order approximation scheme.

Afterwards, the integration proceeds with a new basic stepsize. Such a numerical approach with a base method (here: linearly implicit Euler) is called extrapolation method in general. A representation as extrapolation tableau for this approach for the linearly implicit Euler method of order p is given by $\begin{aligned} &T_{1,1} \\ &T_{2,1} && T_{2,2} \\ &T_{3,1} && T_{3,2} && T_{3,3} \\ & \quad \vdots && && \qquad \ddots \end{aligned}$ where $T_{1,1}, T_{2,1}, T_{3,1}, \ldots$ correspond to order p from our base method and $T_{j_{max},j_{max}}$ is the desired final value. The value, $j_{max}$ from the algorithmic pseudocode defines us the maximum number of rows and is called the extrapolation order. For each $j=1,\ldots,j_{max}$ , a linearly implicit Euler solution is computed on an equidistant time grid with shorter and shorter time sub-steps. At the joint final time, the value is extrapolated. Setting the extrapolation order to 0, only one step will be performed. Then, the algorithm corresponds simply to the base method which is the linearly implicit Euler.

Error control and stepsize adaption

In each step of the time integration, a new basic stepsize has to be chosen. The temporal discretization error $\epsilon := ||T_{j,j}-T_{j,j-1}|| \leq \mathrm{eps}$ for the diagonal entries $T_{j,j}$ in a suitable norm is estimated by taking the difference of the diagonal and subdiagonal extrapolation values of maximal order. This error criterion with $\mathrm{eps}$ as a user defined relative accuracy is then also called subdiagonal error criterion and used as basis for joint order and stepsize control after each time integration step.

To define an appropriate accuracy criterion, two important parameters for controlling the adaptive process of the stepsizes are the absolute tolerance atolT and the relative tolerance rTolT. Besides the subdiagonal error criterion, the local temporal error $\mathrm{eps}$ of the current approximate solution $u$ satisfies $|\mathrm{eps}| \leq \mathrm{aTolT} + \mathrm{rTolT} |u|$ , where the accuracy is chosen as $\mathrm{aTolT} + \mathrm{rTolT} |u|$ . Here, the relative tolerance becomes more important for large solutions and the absolute tolerance becomes more important for small solutions. (Rule of thumb: If m is the number of digits required for solution $u$ , set relative tolerance rTol to $10^{-(m+1)}$ and set absolute tolerance aTol to a value at which $|u|$ is insignificant. In many well scaled problems a reasonable setting is aTol = rTol.)

For fixed index $j$ and current basic stepsize $H$ , a new stepsize estimate based on the accuracy will be chosen by $H_j = H \bigl(\rho \frac{\mathrm{eps}}{\epsilon}\bigr)^{\frac{1}{p_j}}$ with $p_j=\gamma j+1$ and $\rho \in (0,1]$ a constant safety factor, describing how much the error can be underestimated. To prevent the stepsize factor $\bigl(\rho \frac{\mathrm{eps}}{\epsilon}\bigr)^{\frac{1}{p_j}}$ from growing or shrinking too much, the factor may be bounded to a value in $(0,1]$ from below and a value in $[1,q)$ from above, $q$ depending on how much the factor should be allowed to grow.

Implementation in Kaskade

Have a look at the time stepping tutorial and the limex header files in Kaskade’s timestepping/ folder for implementation and usage of the Limex method. The best choice is to use limexWithoutJens.hh for calculations. It works as well for non-constant matrix $B$ and is based on the ParabolicEquation interface and the class SemiImplicitEulerStep. In addition to limex.hh, it contains a mesh adaption loop for the first iteration and allows the choice of a direct or iterative solver as well as solver type and preconditioner type.

The procedure is as follows:

construct LIMEX integrator for time integration,
loop over number of desired steps to perform time stepping with LIMEX and
adapt basic stepsize by estimated errors after each step.

First, the LIMEX integrator has to be constructed. Then, the time integration steps are performed. The limexWithoutJens.hh method for executing one step of the time integration consists of a loop for the solution of $(B-\tau_j L’(u_0))^{-1}\tau_j L(u^k)$ and a nested propagation loop $(B-\tau_j L’(u_0))^{-1}\tau_j L(u^k+\delta u^k)$ , followed by the insertion of the intermediate results and the current step sizes into the extrapolation tableau. It also provides a method for error estimation to adapt the basic stepsize after each step of the LIMEX.

The method limexWithoutJensHierarchicEst.hh is used with an hierarchic error estimator for mesh adaption, since the embedded error estimator cannot be used for linear finite elements.

Page last modified: 2022-01-26 15:49:29 +0100 CET