Solving the Mixed form of the Heat Equation

This shows a different variational formulation of the heat equation. As before, we put \(\Omega = [0,10] \times [0,10]\), with boundary \(\Gamma\): but now we write the heat equation in mixed form:

\[ \begin{align}\begin{aligned}\sigma + \nabla u & = 0\\u_t + \nabla \cdot \sigma & = f\end{aligned}\end{align} \]

which gives rise to the weak form

\[ \begin{align}\begin{aligned}(\sigma, v) - (u, \nabla \cdot v) &= 0\\(u_t, w) + (\nabla \cdot \sigma, w) &= (f, w)\end{aligned}\end{align} \]

Here, \(\sigma\) is a vector-valued variable and will be discretized in a suitable \(H(\mathrm{div})\) -conforming space. The variable \(u\) actually only requires \(L^2\) regularity and will be discretized with discontinuous piecewise polynomials.

Note that this gives us a differential-algebraic system at the fully discrete level as there is no time derivative on \(\sigma\).

Standard imports, although we’re using a different RK scheme this time:

from firedrake import *
from irksome import LobattoIIIC, Dt, MeshConstant, TimeStepper
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

Create time and time-step variables:

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

As in the first heat demo, we build the RHS via 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)))
sigexact = -grad(uexact)

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

Set up the initial condition:

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

And define the variational form:

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)

As before, we use a sparse direct method:

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

We set the time stepper as before, except there are no strongly-enforced boundary conditiuons for the mixed method:

stepper = TimeStepper(F, butcher_tableau, t, dt, sigu,
                      solver_parameters=params)

And we advance the solution in time:

while (float(t) < 1.0):
    if (float(t) + float(dt) > 1.0):
        dt.assign(1.0 - float(t))
    stepper.advance()
    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'))