# Heat equation¶

We will be interested in solving heat equation:

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:

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. |

Task 1

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

## First steps¶

Consider data

Task 2

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.

Define all relevant data from(2). Use`Constant`

or`Expression`

classes to define \(f\), \(g\), \(u_0\).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)

Define bilinear and linear forms describing Galerkin descretization of the weak formulation derived inTask 1on 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.

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)

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(...)`

.

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.

Task 3

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

Task 3

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

- subclass
`SubDomain`

, - define
`MeshFunction`

, - mark the mesh function using
`SubDomain.mark`

method, - define integration
`Measure`

.

- subclass

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], 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)
```

## Time-dependent BC¶

Consider time-dependent data

Task 4

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

Task 5

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

## Adaptive time-stepping¶

Consider solution of *low* precision generated by timestep
\(\Delta t\):

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

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

where

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

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:

compute \(u^{n+1}_\mathrm{low}\) and \(u^{n+1}_\mathrm{high}\)

compute \(\eta\)

compute \(\Delta t^*\)

if \(\eta\leq\mathrm{Tol}\):\(u^{n+1}:=u^{n+1}_\mathrm{high}\)\(n \mathrel{+}= 1\)update timestep \(\Delta t:=\Delta t^*\)

## Wave equation¶

Now consider problem

This problem can be discretized in time as

The iteration can be bootstrapped by

Task 7

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

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**

```
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
+ inner(grad(theta*u + (1-theta)*u_old), grad(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 table header
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 table header
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], 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):
"""Create or update data for Task 5 and Task 6"""
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[0]", degree=1)
# Run all problems
theta = 1
print("Task 2, theta =", theta)
dt = 0.1
timestepping(V, ds_left, theta, T, dt, u_0, get_data_2)
theta = 1/2
print("Task 2, theta =", theta)
dt = 0.1
timestepping(V, ds_left, theta, T, dt, u_0, get_data_2)
theta = 0
print("Task 2, theta =", theta)
dt = 0.1
timestepping(V, ds_left, theta, T, dt, u_0, get_data_2)
theta = 1
print("Task 3, theta =", theta)
dt = 0.1
timestepping(V, ds_left, theta, T, dt, u_0, get_data_3)
theta = 1
print("Task 4, theta =", theta)
dt = 0.1
timestepping(V, ds_left, theta, T, dt, u_0, get_data_4)
theta = 1
print("Task 5, theta =", theta)
dt = 0.1
timestepping(V, ds_left, theta, T, dt, u_0, get_data_5)
theta = 1
print("Task 6, theta =", theta)
tol = 1e-3
timestepping_adaptive(V, ds_left, theta, T, tol, u_0, get_data_5)
theta = 1/2
print("Task 6, theta =", theta)
tol = 1e-3
timestepping_adaptive(V, ds_left, theta, T, tol, u_0, get_data_5)
# Hold plots before quitting
plt.show()
```