Getting Started with Time Integrators

In this example, we go through the basic setup needed when using the TimeIntegrators module in Mantis. We will use a very simple ODE and show how to setup different time integrators for this problem.

Background knowledge

We will consider the following ODE:

$$\frac{dy}{dt}=\lambda y,y(t=0)=1.0,$$

which has exact solution $y(t)=e^{\lambda t}$. We will set $\lambda$ as

const lambda = -4

-4

Packages

As the goal of this example is to setup basic time integration problems in Mantis, we will only use Mantis for now.

using Mantis

Explicit Integrators

In this section, we will treat the previously introduced ODE with an TimeIntegrators.Explicit integrator. We will go through the setup step-by-step.

Step 1: Defining an explicit ODE

For an explicit integrator, we need to provide a function that provides the explicit part of the ODE. In the terminology of the TimeIntegrators-module, this is the function $F$. In our case, we simply have $F(y,t)=\lambda y$.

This is implemented as follows.

Make sure to overwrite the first argument (here called output).

As explained in the TimeIntegrators.TimeIntegrationOperators docstring, the evaluate function for an explicit ODE should have three input arguments, and the first one must be overwritten. Not overwriting the first argument will lead to incorrect results.

function ode_explicit_function!(output, yn, t)
    for n in eachindex(output)
        output[n] = lambda * yn[n]
    end
    return nothing
end

ode_explicit_function! (generic function with 1 method)

We can use this newly defined function and pass it to the TimeIntegrators.define_explicit_ode-function, which will create the required TimeIntegrators.TimeIntegrationOperators object.

ode_explicit = TimeIntegrators.define_explicit_ode(ode_explicit_function!)

Mantis.TimeIntegrators.TimeIntegrationOperators{typeof(Main.ode_explicit_function!), Nothing, Nothing}(Main.ode_explicit_function!, nothing, nothing)

Step 2: Picking and initialising a scheme

Now we can pick an TimeIntegrators.Explicit integrator and initialise it. We will pick the provided TimeIntegrators.RK4 method. Since this is a multi-stage but single-step method, the initialisation procedure for this integrator is simple.

const method_ex = TimeIntegrators.RK4

const y₀_ex = [1.0]
const y_n_ex = TimeIntegrators.initialise_scheme(y₀_ex, method_ex)

Mantis.TimeIntegrators.TimeIntegrationSolution{Mantis.TimeIntegrators.Explicit{4, 1, Float64, 16, 4, 1}, Nothing, Float64, Nothing}(1, [1.0;;], Mantis.TimeIntegrators.Explicit{4, 1, Float64, 16, 4, 1}([0.0 0.0 0.0 0.0; 0.5 0.0 0.0 0.0; 0.0 0.5 0.0 0.0; 0.0 0.0 1.0 0.0], [0.16666666666666666 0.3333333333333333 0.3333333333333333 0.16666666666666666], [1.0; 1.0; 1.0; 1.0;;], [1.0;;], [0.0, 0.5, 0.5, 1.0], Mantis.TimeIntegrators.TimeLevels([0], Int64[], Int64[]), 4), nothing, 0, [6.92898168621755e-310;;], [0.0 0.0 0.0 0.0], [0.0 0.0 0.0 0.0], nothing, [0.0], [0.0])

Step 3: Integrating the ODE

Now we can integrate our ODE. Say we want to integrate for 100 time steps with a time step size of 0.1, then our code would look like this.

const n_steps_ex = 100
const dt_ex = 0.1
t_ex = 0.0
for step in 1:100
    global t_ex
    # Call this function to advance the solution. Note that the input time is the current
    # time.
    TimeIntegrators.time_integrate!(y_n_ex, ode_explicit, t_ex, dt_ex)
    t_ex += dt_ex
end

Step 2 take 2: Picking and initialising a scheme

If we instead want to use a multi-step method, we will need to be a little more careful when initialising the scheme. Say we pick the provided TimeIntegrators.AB2 method. We will also have to pick a startup scheme. In this case, we pick TimeIntegrators.HEUN2. Note that a startup scheme cannot be a multi-step scheme itself. Also, make sure that the startup scheme is accurate enough to initialise your multi-step scheme, or you may lose accuracy.

const method_ex2 = TimeIntegrators.AB2
const startup_method_ex2 = TimeIntegrators.HEUN2

const y₀_ex2 = [1.0]
const y_n_ex2 = TimeIntegrators.initialise_scheme(y₀_ex2, method_ex2, startup_method_ex2)

Mantis.TimeIntegrators.TimeIntegrationSolution{Mantis.TimeIntegrators.Explicit{1, 2, Float64, 1, 2, 4}, Mantis.TimeIntegrators.Explicit{2, 1, Float64, 4, 2, 1}, Float64, Mantis.TimeIntegrators.TimeIntegrationSolution{Mantis.TimeIntegrators.Explicit{2, 1, Float64, 4, 2, 1}, Nothing, Float64, Nothing}}(1, [1.0 0.0], Mantis.TimeIntegrators.Explicit{1, 2, Float64, 1, 2, 4}([0.0;;], [1.5; 1.0;;], [1.0 0.0], [1.0 -0.5; 0.0 0.0], [0.0], Mantis.TimeIntegrators.TimeLevels([0], Int64[], [1]), 2), Mantis.TimeIntegrators.Explicit{2, 1, Float64, 4, 2, 1}([0.0 0.0; 1.0 0.0], [0.5 0.5], [1.0; 1.0;;], [1.0;;], [0.0, 1.0], Mantis.TimeIntegrators.TimeLevels([0], Int64[], Int64[]), 2), 1, [5.0e-324 2.87e-322], [0.0;;], [0.0;;], Mantis.TimeIntegrators.TimeIntegrationSolution{Mantis.TimeIntegrators.Explicit{2, 1, Float64, 4, 2, 1}, Nothing, Float64, Nothing}(1, [1.0;;], Mantis.TimeIntegrators.Explicit{2, 1, Float64, 4, 2, 1}([0.0 0.0; 1.0 0.0], [0.5 0.5], [1.0; 1.0;;], [1.0;;], [0.0, 1.0], Mantis.TimeIntegrators.TimeLevels([0], Int64[], Int64[]), 2), nothing, 0, [5.0e-324;;], [0.0 0.0], [0.0 0.0], nothing, [0.0], [0.0]), [0.0], [0.0])

Diagonally Implicit Integrators

In this section, we will treat the previously introduced ODE with an TimeIntegrators.DiagonallyImplicit integrator. We will go through the setup step-by-step.

Step 1: Defining an diagonally Implicit ODE

For a diagonally implicit integrator, we need to provide a solver that solves the implicit part of the ODE. That is, in the terminology of the TimeIntegrators-module, we need to solve $\mathbf{Y}-h\mathbf{g}(\mathbf{Y})=\mathbf{x}$. If $\mathbf{g}(\mathbf{Y})$ is linear, as it is in this example, we can use a direct solver. This is implemented as follows. Note that, for efficiency reasons, you may want to use solver that directly overwrite the output argument. For simplicity, that is not done here.

Make sure to overwrite the first argument (here called output).

As explained in the TimeIntegrators.TimeIntegrationOperators docstring, the solver function for a diagonally implicit ODE should have four input arguments, and the first one must be overwritten. Not overwriting the first argument will lead to incorrect results.

An implicit evaluate function may be needed.

As explained in the TimeIntegrators.define_diagonally_implicit_ode docstring, if you use a multi-step scheme or a scheme with zeros on the diagonal, you will also have to provide an implicit evaluate function.

import LinearAlgebra
function implicit_solve!(output, x, h, t)
    output .= (LinearAlgebra.I - h * lambda) \ x
    return nothing
end

implicit_solve! (generic function with 1 method)

We can use this newly defined function and pass it to the TimeIntegrators.define_diagonally_implicit_ode-function, which will create the required TimeIntegrators.TimeIntegrationOperators object.

ode_diag_implicit = TimeIntegrators.define_diagonally_implicit_ode(implicit_solve!)

Mantis.TimeIntegrators.TimeIntegrationOperators{Nothing, typeof(Main.implicit_solve!), Nothing}(nothing, Main.implicit_solve!, nothing)

Step 2: Picking and initialising a scheme

Now we can pick a TimeIntegrators.DiagonallyImplicit integrator and initialise it. We will pick the provided TimeIntegrators.DIRK3 method. Since this is a multi-stage but single-step method, the initialisation procedure for this integrator is simple.

const method_di = TimeIntegrators.DIRK3

const y₀_di = [1.0]
const y_n_di = TimeIntegrators.initialise_scheme(y₀_di, method_di)

Mantis.TimeIntegrators.TimeIntegrationSolution{Mantis.TimeIntegrators.DiagonallyImplicit{2, 1, Float64, 4, 2, 1}, Nothing, Float64, Nothing}(1, [1.0;;], Mantis.TimeIntegrators.DiagonallyImplicit{2, 1, Float64, 4, 2, 1}([0.7886751345948129 0.0; -0.5773502691896258 0.7886751345948129], [0.5 0.5], [1.0; 1.0;;], [1.0;;], [0.7886751345948129, 0.21132486540518708], Mantis.TimeIntegrators.TimeLevels([0], Int64[], Int64[]), 3), nothing, 0, [4.0e-323;;], [0.0 0.0], [0.0 0.0], nothing, [0.0], [0.0])

Step 3: Integrating the ODE

Now we can integrate our ODE. Say we want to integrate for 100 time steps with a time step size of 0.1, then our code would look like this. Note that this has the same structure as the explicit case.

const n_steps_di = 100
const dt_di = 0.1
t_di = 0.0
for step in 1:100
    global t_di
    # Call this function to advance the solution. Note that the input time is the current
    # time.
    TimeIntegrators.time_integrate!(y_n_di, ode_diag_implicit, t_di, dt_di)
    t_di += dt_di
end

Fully Implicit Integrators

In this section, we will treat the previously introduced ODE with an TimeIntegrators.Implicit integrator. We will go through the setup step-by-step.

Step 1: Defining an Implicit ODE

For an implicit integrator, we need to provide both a solver that solves the implicit part of the ODE as well as an evaluate function for this implicit part. In the terminology of the TimeIntegrators-module, we need to solve $\mathbf{Y}-h\mathbf{g}(\mathbf{Y})=\mathbf{x}$. If $\mathbf{g}(\mathbf{Y})$ is linear, as it is in this example, we can use a direct solver. Note that, for efficiency reasons, you may want to use a solver that directly overwrites the output argument. For simplicity, that is not done here. Additionally, we need to define the function $G(y,t)=\lambda y$. This is implemented as follows.

Make sure to overwrite the first argument (here called output).

As explained in the TimeIntegrators.TimeIntegrationOperators docstring, the solver function for an implicit ODE should have four input arguments, and the first one must be overwritten. Not overwriting the first argument will not lead to correct results. The evaluate function should have three argument, and the first one must also be overwritten.

function implicit_solve!(output, x, h, t)
    output .= (LinearAlgebra.I - h * lambda) \ x
    return nothing
end

function implicit_evaluate!(output, y, t)
    for n in eachindex(output, y)
        output[n] = lambda * y[n]
    end
    return nothing
end

implicit_evaluate! (generic function with 1 method)

We can use this newly defined function and pass it to the TimeIntegrators.define_implicit_ode-function, which will create the required TimeIntegrators.TimeIntegrationOperators object.

ode_implicit = TimeIntegrators.define_implicit_ode(implicit_solve!, implicit_evaluate!)

Mantis.TimeIntegrators.TimeIntegrationOperators{Nothing, typeof(Main.implicit_solve!), typeof(Main.implicit_evaluate!)}(nothing, Main.implicit_solve!, Main.implicit_evaluate!)

Step 2: Picking and initialising a scheme

Now we can pick an TimeIntegrators.Implicit integrator and initialise it. We will pick the provided TimeIntegrators.GAUSS_LEGENDRE_4 method. Since this is a multi-stage but single-step method, the initialisation procedure for this integrator is simple.

const method_impl = TimeIntegrators.GAUSS_LEGENDRE_4

const y₀_impl = [1.0]
const y_n_impl = TimeIntegrators.initialise_scheme(y₀_impl, method_impl)

Mantis.TimeIntegrators.TimeIntegrationSolution{Mantis.TimeIntegrators.Implicit{2, 1, Float64, 4, 2, 1}, Nothing, Float64, Nothing}(1, [1.0;;], Mantis.TimeIntegrators.Implicit{2, 1, Float64, 4, 2, 1}([0.25 -0.038675134594812866; 0.5386751345948129 0.25], [0.5 0.5], [1.0; 1.0;;], [1.0;;], [0.21132486540518713, 0.7886751345948129], Mantis.TimeIntegrators.TimeLevels([0], Int64[], Int64[]), 4), nothing, 0, [6.92898056555334e-310;;], [0.0 0.0], [0.0 0.0], nothing, [0.0], [0.0])

Step 3: Integrating the ODE

Now we can integrate our ODE. Say we want to integrate for 100 time steps with a time step size of 0.1, then our code would look like this. Note that this again has the same structure as in the previous cases.

const n_steps_impl = 100
const dt_impl = 0.1
t_impl = 0.0
for step in 1:100
    global t_impl
    # Call this function to advance the solution. Note that the input time is the current
    # time.
    TimeIntegrators.time_integrate!(y_n_impl, ode_implicit, t_impl, dt_impl)
    t_impl += dt_impl
end

IMEX Integrators

In this section, we will treat the previously introduced ODE with an TimeIntegrators.IMEX integrator. We will go through the setup step-by-step.

Step 1: Defining an IMEX ODE

IMEX integrators can treat one part of the ODE with an implicit step, while treating the rest with an explicit step. The stiffness of the terms in your ODE often influences this split. In this simple example, we introduce the following split:

$$\frac{dy}{dt}=0.5\lambda y+0.5\lambda y.$$

That is, the implicit and explicit parts are the same.

For an IMEX integrator, we need to provide a function for evaluating the explicit part, here $F(y,t)=0.5\lambda y$, and a solver for the implicit part, just like in the diagonally implicit case. This is implemented as follows.

Make sure to overwrite the first argument (here called output).

As explained in the TimeIntegrators.TimeIntegrationOperators docstring, the evaluate function for an explicit ODE should have three input arguments, and the first one must be overwritten. Not overwriting the first argument will not lead to correct results.

An implicit evaluate function may be needed.

As explained in the TimeIntegrators.define_imex_ode docstring, if you use a multi-step scheme or a scheme with zeros on the diagonal (as done here), you will also have to provide an implicit evaluate function.

function explicit_imex_function!(output, yn, t)
    for n in axes(output, 1)
        output[n] = 0.5 * lambda * yn[n]
    end
    return nothing
end

function implicit_imex_solver!(output, x, h, t)
    LinearAlgebra.ldiv!(output, LinearAlgebra.lu(LinearAlgebra.I - 0.5 * h * lambda), x)
    return nothing
end

function implicit_imex_function!(output, yn, t)
    for n in axes(output, 1)
        output[n] = 0.5 * lambda * yn[n]
    end
    return nothing
end

implicit_imex_function! (generic function with 1 method)

We can use this newly defined function and pass it to the TimeIntegrators.define_imex_ode-function, which will create the required TimeIntegrators.TimeIntegrationOperators object.

ode_imex = TimeIntegrators.define_imex_ode(
    explicit_imex_function!, implicit_imex_solver!, implicit_imex_function!
)

Mantis.TimeIntegrators.TimeIntegrationOperators{typeof(Main.explicit_imex_function!), typeof(Main.implicit_imex_solver!), typeof(Main.implicit_imex_function!)}(Main.explicit_imex_function!, Main.implicit_imex_solver!, Main.implicit_imex_function!)

Step 2: Picking and initialising a scheme

Now we can pick an TimeIntegrators.IMEX integrator and initialise it. We will pick the provided TimeIntegrators.RK3_IMEX method. Since this is a multi-stage but single-step method, the initialisation procedure for this integrator is simple.

const method_imex = TimeIntegrators.RK3_IMEX

const y₀_imex = [1.0]
const y_n_imex = TimeIntegrators.initialise_scheme(y₀_imex, method_imex)

Mantis.TimeIntegrators.TimeIntegrationSolution{Mantis.TimeIntegrators.IMEX{3, 1, Float64, 9, 3, 1}, Nothing, Float64, Nothing}(1, [1.0;;], Mantis.TimeIntegrators.IMEX{3, 1, Float64, 9, 3, 1}([0.0 0.0 0.0; 0.0 0.7886751345948128 0.0; 0.0 -0.5773502691896255 0.7886751345948128], [0.0 0.0 0.0; 0.7886751345948128 0.0 0.0; -0.21132486540518725 0.4226497308103745 0.0], [0.0 0.5 0.5], [0.0 0.5 0.5], [1.0; 1.0; 1.0;;], [1.0;;], [0.0, 0.7886751345948128, 0.21132486540518725], [0.0, 0.7886751345948128, 0.21132486540518725], Mantis.TimeIntegrators.TimeLevels([0], Int64[], Int64[]), 3), nothing, 0, [1.0e-323;;], [0.0 0.0 0.0], [0.0 0.0 0.0], nothing, [0.0], [0.0])

Step 3: Integrating the ODE

Now we can integrate our ODE. Say we want to integrate for 100 time steps with a time step size of 0.1, then our code would look like this. Note that this again has the same structure as in the previous cases.

const n_steps_imex = 100
const dt_imex = 0.1
t_imex = 0.0
for step in 1:100
    global t_imex
    # Call this function to advance the solution. Note that the input time is the current
    # time.
    TimeIntegrators.time_integrate!(y_n_imex, ode_imex, t_imex, dt_imex)
    t_imex += dt_imex
end

---

This page was generated using Literate.jl.