# Heat equation¶

We will be interested in solving heat equation:

\begin{align}\begin{aligned}u_t - \Delta u &= f &&\quad\text{ in }\Omega\times(0, T),\\\tfrac{\partial u}{\partial\mathbf{n}} &= g &&\quad\text{ on }\partial\Omega\times(0, T),\\u &= u_0 &&\quad\text{ on }\Omega\times\{0\}\end{aligned}\end{align}

using $$\theta$$-scheme discretization in time and arbitrary FE discretization in space with given data $$f$$, $$g$$, $$u_0$$. $$\theta$$-scheme time-discrete heat equation reads:

(1)\begin{align}\begin{aligned}\frac{1}{\Delta t} \Bigl(u^{n+1} - u^n\Bigr) - \theta\Delta u^{n+1} - (1-\theta)\Delta u^n &= \theta f(t_{n+1}) + (1-\theta) f(t_n) &&\quad\text{ in }\Omega, \; n=0,1,2,\ldots\\\tfrac{\partial u^n}{\partial\mathbf{n}} &= g(t_n) &&\quad\text{ on }\partial\Omega, \; n=0,1,2,\ldots\\u^0 &= u_0 &&\quad\text{ in }\Omega\end{aligned}\end{align}

for a certain sequence $$0=t_0 < t_1 < t_2 < ... \leq T$$. Special cases are:

 $$\theta=0$$ explicit Euler scheme, $$\theta=\frac12$$ Crank-Nicolson scheme, $$\theta=1$$ implicit Euler scheme.

Test (1) by functions from $$H^1(\Omega)$$ and derive a weak formulation of $$\theta$$-scheme for heat equation.

## First steps¶

Consider data

(2)\begin{align}\begin{aligned}\Omega &= (0, 1)^2,\\T &= 2,\\f &= 0,\\g &= 0,\\u_0(x,y) &= x.\end{aligned}\end{align}

Write FEniCS code implementing problem (1), (2), assuming general $$\theta$$, and arbitrary but fixed $$\Delta t$$. In particular assume:

from dolfin import *

mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

theta = Constant(0.5)
dt = Constant(0.1)


Proceed step-by-step.

1. Define all relevant data from (2). Use Constant or Expression classes to define $$f$$, $$g$$, $$u_0$$.

2. Define a finite element function for holding solution at a particular time step:

u_n = Function(V)


and arguments of linear and bilinear forms:

u, v = TrialFunction(V), TestFunction(V)

3. Define bilinear and linear forms describing Galerkin descretization of the weak formulation derived in Task 1 on the space V.

You can conveniently mix bilinear and linear terms into a single expression:

F = 1/dt*(u - u_n)*v*dx + ...


and separate bilinear and linear part using lhs, rhs:

a, L = lhs(F), rhs(F)


Tip

It is good to execute your code every once in a while, even when it is not doing anything useful so far, e.g., does not have time-stepping yet. You will catch the bugs early and fix them easily.

4. Prepare for the beggining of time-stepping. Assume u0 is an Expression or Constant. You can use Function.interpolate() or interpolate():

u_n.interpolate(u0)
# or
u_n = interpolate(u0, V)

5. Implement time-stepping. Write a control flow statement (for example a while loop) which executes the solver for problem a == L repeatedly while updating what needed.

Hint

Note that a single Function object is needed to implement the time-stepping. The function can be used to hold the value of $$u_n$$ and then be updated by calling solve(...).

6. Run with different values of $$\theta=1,\frac12,0$$.

As a first indicator of correctness of the implementation you can drop into the loop lines like:

energy = assemble(u_n*dx)
print("Energy =", energy)


Are you observing expected value?

## Data IO, plotting¶

There are several possibilities for visualization of data.

XDMF output and Paraview

One possibility is to use IO facilities of FEniCS and visualize using external software, for example Paraview.

Note

This approach allows to separate

• actual computation, which can happen in headless HPC environment, for example big parallel clusters of thousands of CPU cores,
• and visualization, which many times needs human interaction.

One can used XDMFFile to store data:

# Open file for XDMF IO
f = XDMFFile('solution.xdmf')

while t < T:

# Compute time step
perform_timestep(u_n, t, dt)
t += dt

# Save the result to file at time t
f.write(u_n, t)


Then you can open Paraview by shell command

paraview &


and visualize the file solution.xdmf.

Matplotlib – native plotting in Python

Another possibility is to use Python plotting library Matplotlib.

Note

Matplotlib is Python native plotting library, which is programmable and supports

• interactive use from Python interpreters, including popular shells like Jupyter,
• high-quality vector output suitable for scientific publishing.

FEniCS plot(obj, **kwargs) function implements plotting using Matplotlib for several different types of obj, for instance Function, Expression, Mesh, MeshFunction. As Matplotlib is highly programmable and customizable, FEniCS plot() is typically accompanied by some native matplotlib commands. Mimimal example of interaction of FEniCS and matplotlib:

from dolfin import *
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(64, 64)
plot(mesh)
plt.savefig('mesh_64_64.pdf')  # Render to PDF
plt.show()  # Render into interactive window


Add something along the lines of:

import matplotlib.pyplot as plt

# Open a plot window
fig = plt.figure()
fig.show()

while t < T:

# Compute time step
perform_timestep(u_n, t, dt)
t += dt

# Update plot to current time step
fig.clear()
p = plot(u_n, mode="warp")
fig.colorbar(p)
fig.gca().set_zlim((0, 2))
fig.canvas.draw()


Warning

Matplotlib’s interactive capabalities aparently depend on used Matplotlib backend. In particular updating the contents of the plot window seems to work fine with TkAgg backend. Issue shell command

export MPLBACKEND=tkagg


to choose TkAgg in the current shell session.

Implement at least one of the aforementioned ways to plot your solutions in time. Check that your solution of Task 2 looks reasonable.

## Nonhomogeneous Neumann BC¶

Consider (1), (2) but now with nonhomogeneous Neumann data

(3)\begin{align}\begin{aligned}g &= 1 \text{ on } \{ x = 0 \},\\g &= 0 \text{ elsewhere}.\end{aligned}\end{align}

1. Derive weak formulation describing (1), (2), (3).

2. Define surface measure supported on the left boundary of the unit square mesh by following steps:

1. subclass SubDomain,
2. define MeshFunction,
3. mark the mesh function using SubDomain.mark method,
4. define integration Measure.

Hint

Show/Hide Code

# Define instance of SubDomain class
class Left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0)
left = Left()

# Define and mark mesh function on facets
facets = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
left.mark(facets, 1)

# Define exterior facet measure where facets==1
ds_left = Measure("ds", mesh, subdomain_data=facets, subdomain_id=1)

1. Using the surface measure, modify the implementation from Task 2 to incorporate boundary condition (3).
2. Run the code with $$\theta=1$$ and check that the results look as expected.

## Time-dependent BC¶

Consider time-dependent data

(4)\begin{align}\begin{aligned}f(x, t) &= 2-t,\\g(x, t) &= \left\{\begin{array}{ll} t & x = 0, \newline 0 & \text{otherwise}. \end{array}\right.\end{aligned}\end{align}

Modify solution of the previous task to use data (4).

Hint

You can use Constant.assign() or Expression.<param> = <value> to change existing Constant or Expression. Look for User defined parameters in Expression documentation.

Now consider different time-dependent data

(5)\begin{align}\begin{aligned}f(x, t) &= 0,\\g(x, t) &= \left\{\begin{array}{ll} \max\bigl(0, \tfrac{1-t}{2}\bigr) & x = 0, \newline 0 & \text{otherwise}. \end{array}\right.\end{aligned}\end{align}

Modify solution of the previous task to use data (5).

Consider solution of low precision generated by timestep $$\Delta t$$:

(6)$\frac{1}{\Delta t} \Bigl(u^{n+1}_\mathrm{low} - u^n\Bigr) - \theta\Delta u^{n+1}_\mathrm{low} - (1-\theta)\Delta u^n = \theta f(t_{n+1}) + (1-\theta) f(t_n)$

and solution of high precision computed by two timesteps of a half size:

(7)\begin{align}\begin{aligned}\frac{1}{\Delta t/2} \Bigl(u^{n+1/2}_\mathrm{high} - u^n\Bigr) - \theta\Delta u^{n+1/2}_\mathrm{high} - (1-\theta)\Delta u^n &= \theta f(t_{n+1/2}) + (1-\theta) f(t_n),\\\frac{1}{\Delta t/2} \Bigl(u^{n+1}_\mathrm{high} - u^{n+1/2}_\mathrm{high}\Bigr) - \theta\Delta u^{n+1}_\mathrm{high} - (1-\theta)\Delta u^{n+1/2}_\mathrm{high} &= \theta f(t_{n+1}) + (1-\theta) f(t_{n+1/2}).\end{aligned}\end{align}

By Richardson extrapolation one can estimate the error of discretization (in time) by quantity:

(8)$\eta := \frac{\|u^{n+1}_\mathrm{high} - u^{n+1}_\mathrm{low}\|_{L^2(\Omega)}} {2^p - 1}$

where

(9)$p = \left\{\begin{array}{ll} 2 && \theta=\tfrac12, \newline 1 && \text{otherwise} \end{array}\right.$

is a theoretical order of accuracy of the $$\theta$$-scheme. Given a tolerance $$\mathrm{Tol}$$ set the new timestep to

(10)$\Delta t^* := \left( \frac{\rho\,\mathrm{Tol}}{\eta} \right)^\frac1p \Delta t.$

Here $$0<\rho\leq1$$ is a chosen safety factor. That asymptotically ensures that the error (or at least the estimator) committed with the new time step is $$\rho$$-multiple of the tolerance.

Now consider an algorithm:

1. compute $$u^{n+1}_\mathrm{low}$$ and $$u^{n+1}_\mathrm{high}$$

2. compute $$\eta$$

3. compute $$\Delta t^*$$

4. if $$\eta\leq\mathrm{Tol}$$:
$$u^{n+1}:=u^{n+1}_\mathrm{high}$$
$$n \mathrel{+}= 1$$
5. update timestep $$\Delta t:=\Delta t^*$$

Solve (1), (2)$$_{1,2,5}$$, (5) using the adaptive strategy described above.

Hint

You will need more than one Function and perform assignments between them. Having Functions f, g on the same space you can perform assignment $$f := g$$ by

f.vector()[:] = g.vector()


## Wave equation¶

Now consider problem

(11)\begin{align}\begin{aligned}u_{tt} - \Delta u &= f &&\quad\text{ in }\Omega\times(0, T),\\u &= 0 &&\quad\text{ on }\partial\Omega\times(0, T),\\u(\cdot,t) &= u_0 &&\quad\text{ in }\Omega,\\u_t(\cdot,t) &= v_0 &&\quad\text{ in }\Omega.\end{aligned}\end{align}

This problem can be discretized in time as

(12)$\frac{1}{(\Delta t)^2} \Bigl(u^{n+1} - 2u^n + u^{n-1}\Bigr) -\tfrac12 \Delta(u^{n+1} + u^{n-1}) = f(t^n).$

The iteration can be bootstrapped by

$u^1 := u^0 + \Delta t v^0.$

Implement solver for problem (11) by discretizing in time with (12). Solve the problem with data

\begin{align}\begin{aligned}f &= 0,\\u^0(x,y) &= \max(0, 1-4 r(x,y)) \qquad \text{where } r(x,y)=\operatorname{dist}((x,y),(\tfrac12,\tfrac{3}{10})),\\v^0(x,y) &= 0,\\T &= 5,\\\Omega &= (0,1)^2.\end{aligned}\end{align}

Visualize the result.

## Reference solution¶

Note

The reference solution follows the DRY principle. Hands-on participants are not expected to write such a structured code during the session.

Attention

For on-the-fly plotting, TkAgg Matplotlib backend has been tested. You can enforce its selection by a shell command

export MPLBACKEND=tkagg


Note that the the plotting is the bottleneck of the code. The code runs much faster without plots which can be ensured by

DOLFIN_NOPLOT=1 MPLBACKEND=template python3 heat.py


We leave as an exercise to add XDMF output for plotting in Paraview.

Show/Hide Code

Download Code

from dolfin import *
import matplotlib.pyplot as plt

def create_timestep_solver(get_data, dsN, theta, u_old, u_new):
"""Prepare timestep solver by theta-scheme for given
function get_data(t) returning data (f(t), g(t)), given
solution u_old at time t and unknown u_new at time t + dt.
Return a solve function taking (t, dt).
"""

# Initialize coefficients
f_n, g_n = get_data(0)
f_np1, g_np1 = get_data(0)
idt = Constant(0)

# Extract function space
V = u_new.function_space()

# Prepare weak formulation
u, v = TrialFunction(V), TestFunction(V)
theta = Constant(theta)
F = ( idt*(u - u_old)*v*dx
- (theta*f_np1 + (1-theta)*f_n)*v*dx
- (theta*g_np1 + (1-theta)*g_n)*v*dsN
)
a, L = lhs(F), rhs(F)

def solve_(t, dt):
"""Update problem data to interval (t, t+dt) and
run the solver"""

# Update coefficients to current t, dt
get_data(t, (f_n, g_n))
get_data(t+dt, (f_np1, g_np1))
idt.assign(1/dt)

# Push log level
old_level = get_log_level()
warning = LogLevel.WARNING if cpp.__version__ > '2017.2.0' else WARNING
set_log_level(warning)

# Run the solver
solve(a == L, u_new)

# Pop log level
set_log_level(old_level)

return solve_

def timestepping(V, dsN, theta, T, dt, u_0, get_data):
"""Perform timestepping using theta-scheme with
final time T, timestep dt, initial datum u_0 and
function get_data(t) returning (f(t), g(t))"""

# Initialize solution function
u = Function(V)

# Prepare solver for computing time step
solver = create_timestep_solver(get_data, dsN, theta, u, u)

# Set initial condition
u.interpolate(u_0)

# Open plot window
fig = init_plot()

print("{:10s} | {:10s} | {:10s}".format("t", "dt", "energy"))

# Perform timestepping
t = 0
while t < T:

# Report some numbers
energy = assemble(u*dx)
print("{:10.4f} | {:10.4f} | {:#10.4g}".format(t, dt, energy))

# Perform time step
solver(t, dt)
t += dt

# Update plot
update_plot(fig, u)

def timestepping_adaptive(V, dsN, theta, T, tol, u_0, get_data):
"""Perform adaptive timestepping using theta-scheme with
final time T, tolerance tol, initial datum u_0 and
function get_data(t) returning (f(t), g(t))"""

# Initialize needed functions
u_n = Function(V)
u_np1_low = Function(V)
u_np1_high = Function(V)

# Prepare solvers for computing tentative time steps
solver_low = create_timestep_solver(get_data, dsN, theta, u_n, u_np1_low)
solver_high_1 = create_timestep_solver(get_data, dsN, theta, u_n, u_np1_high)
solver_high_2 = create_timestep_solver(get_data, dsN, theta, u_np1_high, u_np1_high)

# Initial time step; the value does not really matter
dt = T/2

# Set initial conditions
u_n.interpolate(u_0)

# Open plot window
fig = init_plot()

print("{:10s} | {:10s} | {:10s}".format("t", "dt", "energy"))

# Perform timestepping
t = 0
while t < T:

# Report some numbers
energy = assemble(u_n*dx)
print("{:10.4f} | {:10.4f} | {:#10.4g}".format(t, dt, energy))

# Compute tentative time steps
solver_low(t, dt)
solver_high_1(t, dt/2)
solver_high_2(t+dt, dt/2)

# Compute error estimate and new timestep
est = compute_est(theta, u_np1_low, u_np1_high)
dt_new = compute_new_dt(theta, est, tol, dt)

if est > tol:
# Tolerance not met; repeat the step with new timestep
dt = dt_new
continue

# Move to next time step
u_n.vector()[:] = u_np1_high.vector()
t += dt
dt = dt_new

# Update plot
update_plot(fig, u_n)

def compute_est(theta, u_L, u_H):
"""Return error estimate by Richardson extrapolation"""
p = 2 if theta == 0.5 else 1
est = sqrt(assemble((u_L - u_H)**2*dx)) / (2**p - 1)
return est

def compute_new_dt(theta, est, tol, dt):
"""Return new time step"""
p = 2 if theta == 0.5 else 1
rho = 0.9
dt_new = dt * ( rho * tol / est )**(1/p)
return dt_new

def init_plot():
"""Open plot window and return its figure object"""
fig = plt.figure()
fig.show()
return fig

def update_plot(fig, u, zlims=(0, 2)):
"""Plot u in 3D warp mode with colorbar into figure fig;
use zlims as limits on z-axis"""
fig.clear()
p = plot(u, mode="warp")
if p is None:
return
fig.colorbar(p)
fig.gca().set_zlim(zlims)
fig.canvas.draw()

def create_function_space():
"""Return (arbitrary) H^1 conforming function space on
unit square domain"""
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "P", 1)
return V

def create_surface_measure_left(mesh):
"""Return surface measure on the left boundary of unit
square"""
class Left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0)
facets = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
Left().mark(facets, 1)
ds_left = Measure("ds", mesh, subdomain_data=facets, subdomain_id=1)
return ds_left

def get_data_2(t, result=None):
"""Create or update data for Task 2"""
f, g = result or (Constant(0), Constant(0))
f.assign(0)
g.assign(0)
return f, g

def get_data_3(t, result=None):
"""Create or update data for Task 3"""
f, g = result or (Constant(0), Constant(0))
f.assign(1)
g.assign(0)
return f, g

def get_data_4(t, result=None):
"""Create or update data for Task 4"""
f, g = result or (Constant(0), Constant(0))
f.assign(2-t)
g.assign(t)
return f, g

def get_data_5(t, result=None):
f, g = result or (Constant(0), Constant(0))
f.assign(0)
g.assign(max(0, 1-t)/2)
return f, g

if __name__ == '__main__':
# Common data
V = create_function_space()
ds_left = create_surface_measure_left(V.mesh())
T = 2
u_0 = Expression("x", degree=1)

# Run all problems
theta = 1
dt = 0.1
timestepping(V, ds_left, theta, T, dt, u_0, get_data_2)

theta = 1/2
dt = 0.1
timestepping(V, ds_left, theta, T, dt, u_0, get_data_2)

theta = 0
dt = 0.1
timestepping(V, ds_left, theta, T, dt, u_0, get_data_2)

theta = 1
dt = 0.1
timestepping(V, ds_left, theta, T, dt, u_0, get_data_3)

theta = 1
dt = 0.1
timestepping(V, ds_left, theta, T, dt, u_0, get_data_4)

theta = 1
dt = 0.1
timestepping(V, ds_left, theta, T, dt, u_0, get_data_5)

theta = 1