Heat equation

Goals

Learn how to deal with time-dependent problems. Solve heat equation by \(\theta\)-scheme. Solve wave equation with central differences. Plot some nice figures.

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.

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

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

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.

  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.

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

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

Task 3

  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], 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} \]

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

(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} \]

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\):

(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^*\)

Task 6

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

Task 7

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
        + 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()