Solving the Heat Equation with a Multistep Method in Irksome

Consider the heat equation on \(\Omega = [0,10] \times [0,10]\), with boundary \(\Gamma\):

\[ \begin{align}\begin{aligned}u_t - \Delta u &= f\\u & = 0 \quad \textrm{on}\ \Gamma\end{aligned}\end{align} \]

for some known function \(f\). At each time \(t\), the solution to this equation will be some function \(u\in V\), for a suitable function space \(V\).

We transform this into weak form by multiplying by an arbitrary test function \(v\in V\) and integrating over \(\Omega\). We now have the variational problem of finding \(u:[0,T]\rightarrow V\) such that

\[(u_t, v) + (\nabla u, \nabla v) = (f, v)\]

This demo implements an example used by Solin with a particular choice of \(f\) given below

As usual, we need to import firedrake:

from firedrake import *

We will also need to import certain items from irksome:

from irksome import Dt, TimeStepper, BDF

We will use the 3-step Backward Difference Formula:

method = BDF(3)

Now we define the mesh and piecewise linear approximating space in standard Firedrake fashion:

N = 100
x0 = 0.0
x1 = 10.0
y0 = 0.0
y1 = 10.0

msh = RectangleMesh(N, N, x1, y1)
V = FunctionSpace(msh, "CG", 1)

We define variables to store the time step and current time value:

dt = Constant(5.0 / N)
t = Constant(0.0)

This defines the right-hand side using the method of manufactured solutions:

x, y = SpatialCoordinate(msh)
S = Constant(2.0)
C = Constant(1000.0)
B = (x-Constant(x0))*(x-Constant(x1))*(y-Constant(y0))*(y-Constant(y1))/C
R = (x * x + y * y) ** 0.5
uexact = B * atan(t)*(pi / 2.0 - atan(S * (R - t)))
rhs = Dt(uexact) - div(grad(uexact))

We define the initial condition for the fully discrete problem, which will get overwritten at each time step:

u = Function(V)
u.interpolate(uexact)

Now, we will define the semidiscrete variational problem using standard UFL notation, augmented by the Dt operator from Irksome:

v = TestFunction(V)
F = inner(Dt(u), v)*dx + inner(grad(u), grad(v))*dx - inner(rhs, v)*dx
bc = DirichletBC(V, 0, "on_boundary")

We’ll use a basic solver for this demo:

luparams = {"mat_type": "aij",
            "ksp_type": "preonly",
            "pc_type": "lu"}

Irksome’s TimeStepper automates the transformation of our semidiscrete form F into a fully discrete form for the next approximate value and sets up a variational problem to solve at each time step.

An s-step multistep method requires starting values \(u^0,\dots, u^{s - 1}\) in order to begin the approximation. Irksome provides two ways to set these values. If you wish to set the starting values manually, you can create a TimeStepper as shown below, but with no startup_parameters. You can then access the list TimeStepper.us which contains the starting approximations, assign the desired starting values, and advance t by hand. On the other hand, if you wish to use a method which Irksome supports to obtain these starting values, then Irksome allows this process to be completed automatically.

We’ll showcase both methods. The automatic process requires defining a dict of keyword arguments used to setup a TimeStepper representing a single-step method. This TimeStepper is then used to compute the initial approximations needed to start the method. We’ll import RadauIIA and use the backward Euler method to obtain our starting values. Formally, the backward Euler method is only first order accurate, and we wish to use it obtain the starting values for the third-order accurate BDF(3) method. A crude way to increase the accuracy of the starting values is to use a smaller timestep for the startup procedure. Here, we use timesteps of size \(\Delta t / 8\) for the backward Euler method which is accessible through the keyword num_startup_steps. Setting num_startup_steps = 8 will instruct the startup TimeStepper to use 8 substeps to compute each starting value. The keyword stepper_kwargs allows for easy customization of the startup TimeStepper:

from irksome import RadauIIA

startup_stepper_kwargs = {'stage_type': 'value',
                          'solver_parameters': luparams}

startup_parameters = {'tableau': RadauIIA(1),
                      'num_startup_steps': 8,
                      'stepper_kwargs': startup_stepper_kwargs}

We can then set up the TimeStepper as follows:

stepper = TimeStepper(F, method, t, dt, u, bcs=bc,
                               solver_parameters=luparams,
                               startup_parameters=startup_parameters)

Note that the creation of a TimeStepper configured for a multistep method with parameters for the startup procedure will not automatically solve for the required starting values. One must call the TimeSteppers’s startup method, which will internally construct the single-step TimeStepper and solve for the required starting values. It is the user’s responsibility to advance t to t + (s-1)*dt. Unless more refined control is needed, this value is available as TimeStepper.startup_t:

stepper.startup()
t.assign(stepper.startup_t)

If you would rather set the starting values by hand (for instance, in this case the exact solution is available) the process is straightforward. We’ll first reset t so that everything is consistent, and then set the starting values by interpolating the exact solution.:

t.assign(0.0)
stepper.us[0].interpolate(uexact)
t.assign(t + dt)
stepper.us[1].interpolate(uexact)
t.assign(t + dt)
stepper.us[2].interpolate(uexact)

This remaining logic is pretty self-explanatory. We use the TimeStepper’s advance method, which solves the variational problem to compute the next approximate value and updates the solution:

while (float(t) < 1.0):
  if (float(t) + float(dt) > 1.0):
      dt.assign(1.0 - float(t))
  stepper.advance()
  t.assign(float(t) + float(dt))
  print(float(t))

Finally, check the relative \(L^2\) error:

print()
print(norm(u-uexact)/norm(uexact))