Mixed problems with the low-level interface

Updating the solution at each time step is more complicated when we have problems on mixed function spaces. This demo peels back the TimeStepper abstraction in the mixed heat equation demo. In this case, the Dirichlet boundary conditions are weakly enforced. However, mixed problems do not change the way strongly-enforced BC are handled, just how the solution is updated.

Imports:

from firedrake import *
from irksome import LobattoIIIC, Dt, getForm, MeshConstant
from ufl.algorithms.ad import expand_derivatives

butcher_tableau = LobattoIIIC(2)

Build the mesh and approximating spaces:

N = 32
x0 = 0.0
x1 = 10.0
y0 = 0.0
y1 = 10.0
msh = RectangleMesh(N, N, x1, y1)

V = FunctionSpace(msh, "RT", 2)
W = FunctionSpace(msh, "DG", 1)
Z = V * W

MC = MeshConstant(msh)
dt = MC.Constant(10.0 / N)
t = MC.Constant(0.0)

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)))
sigexact = -grad(uexact)

rhs = expand_derivatives(diff(uexact, t) + div(sigexact))

sigu = project(as_vector([0, 0, uexact]), Z)
sigma, u = split(sigu)

v, w = TestFunctions(Z)

F = (inner(Dt(u), w) * dx + inner(div(sigma), w) * dx - inner(rhs, w) * dx
     + inner(sigma, v) * dx - inner(u, div(v)) * dx)

Because we aren’t concerned with any strongly-enforced boundary conditions, we drop that information in calling get_form:

Fnew, k, _, _ = getForm(F, butcher_tableau, t, dt, sigu)

We set up the variational problem and solver using a sparse direct method:

params = {"mat_type": "aij",
          "snes_type": "ksponly",
          "ksp_type": "preonly",
          "pc_type": "lu"}

prob = NonlinearVariationalProblem(Fnew, k)
solver = NonlinearVariationalSolver(prob, solver_parameters=params)

Advancing the solution in time is a bit more complicated now:

num_fields = len(Z)
b = butcher_tableau.b

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

    solver.solve()

    for s in range(butcher_tableau.num_stages):
        for i in range(num_fields):
            sigu.dat.data[i][:] += float(dt) * b[s] * k.dat.data[num_fields * s + i][:]

    print(float(t))
    t.assign(float(t) + float(dt))

Finally, we check the accuracy of the solution:

sigma, u = sigu.split()
print("U error      : ", errornorm(uexact, u) / norm(uexact))
print("Sig error    : ", errornorm(sigexact, sigma) / norm(sigexact))
print("Div Sig error: ",
      errornorm(sigexact, sigma, norm_type='Hdiv')
      / norm(sigexact, norm_type='Hdiv'))