Table of Contents

  • 1  What is Firedrake?
    • 1.1  Background
    • 1.2  Getting started on Google's colab
  • 2  A one-dimensional problem
    • 2.1  Problem data
    • 2.2  Define the mesh
    • 2.3  The FEM space
    • 2.4  The Weak Form
    • 2.5  Boundary conditions
    • 2.6  Solve
    • 2.7  Compute errors
    • 2.8  Plot solutions
  • 3  A fitted mesh in 1D
    • 3.1  Solver as a function
    • 3.2  A 1D fitted mesh
  • 4  Convergence tables
  • 5  A two-dimensional problem
    • 5.1  The mesh
    • 5.2  The function space and weak form.
    • 5.3  Plotting
  • 6  Other things Firedrake can and cannot do

A tutorial on solving singularly perturbed problems in Firedrake¶

21st Workshop on Numerical Methods for Problems with Layer Phenomena

Niall Madden and Sean Tobin, April 2025

You can access this notebook at:

  • Jupyter notebook: https://www.niallmadden.ie/SPPs_Firedrake.ipynb
  • HTML version: https://www.niallmadden.ie/SPPs_Firedrake.html

In this session, we introduce how to solve one-dimensional and two-dimensional singularly perturbed differential equations (SPDEs) on $(0,1)^d$, $d=1,2$ in Firedrake.

For a very basic introduction to Firedrake, see https://niallmadden.ie/IntroToFiredrake. For more advanced topics, see https://www.firedrakeproject.org/documentation.html

While this session will include the very fundamentals of using Firedrake for solving PDEs, there will be an emphasis on topics particular to SPDEs:

  • constructing layer-adapted meshes
  • error estimation in various norms
  • convergence tables
  • extension to higher dimensions

What is Firedrake?¶

According to the documentation: Firedrake is an automated system for the solution of partial differential equations using the finite element method (FEM). Firedrake uses sophisticated code generation to provide mathematicians, scientists, and engineers with a very high productivity way to create sophisticated high performance simulations.

Crucial:

  • Available as a Python module;
  • You can use one of the many built-in mesh generators, or define a mesh (by hand, or import from, say, gmsh);
  • Define a elements, and a function space on that mesh;
  • Express the weak form of a PDE
  • Call a solver
  • Compute errors and/or visualise solutions.

Background¶

Firedrake began in 2014 is as a "reimagining" of FEniCS.

  • It uses the FEniCS Unified Form Language (UFL) to express PDEs.

  • Boasts sophisticated, programmable solvers through "seamless" coupling with PETSc.

  • Triangular, quadrilateral, and tetrahedral unstructured meshes.

  • Vast array of finite element spaces.

Etc.

(To the user, Firedrake may seem a little closer the FEniCS than FEniCSx).

Getting started on Google's colab¶

For serious work, you should install Firedrake on your own system. That is not too hard to do: as of last month, there is a pip install approach.

But to get started, we'll try Google colab: https://colab.google/

Note the incantation needed to install Firedrake. Can take several minutes to run:

In [1]:
try:
    from firedrake import *
except ImportError:
    !wget "https://fem-on-colab.github.io/releases/firedrake-install-release-real.sh" -O "/tmp/firedrake-install.sh" && bash "/tmp/firedrake-install.sh"
    from firedrake import *
firedrake:WARNING OMP_NUM_THREADS is not set or is set to a value greater than 1, we suggest setting OMP_NUM_THREADS=1 to improve performance

We'll also use the numpy, pandas, and matplotlib modules:

In [2]:
import numpy as np
import pandas as pd
# %matplotlib notebook # may not work  on colab
import matplotlib.pyplot as plt

A one-dimensional problem¶

The first problem we'll solve is this ODE: find $u(x)$ such that $$-\varepsilon^2 u'' + r(x)u = f(x) \text{ on } (0,1), \qquad\text{ with } u(0)=u(1)=0.$$

The weak form is $$\underbrace{\varepsilon^2 (u',v') + (ru,v)}_{a(u,v)} = \underbrace{(f,v)}_{l(v)},$$ where here $(u,v):=\int_0^1 u(x)v(x) dx$, is the usual $L^2$ inner product.

The specific example we'll consider will have $r(x) \equiv 1$ and $f(x)=e^x$.

We'll solve this problem using a standard $P_1$ (continuous) Galerkin method.

Problem data¶

I like to use lambda functions for defining functions

In [3]:
epsilon = Constant(0.1)
r = Constant("1.0")
f = lambda x : exp(x)
u_e = lambda x : exp(x)/(1-epsilon**2) + (exp((x-1)/epsilon)*(exp(1.0)-exp(-1.0/epsilon)) + \
      exp(-x/epsilon)*(1 - exp((epsilon-1)/epsilon)))/(1-epsilon**2)*(exp(-2/epsilon) -1) 

Define the mesh¶

We'll use a uniform mesh with $N=16$ intervals.

In [4]:
N = 16 
mesh1D = UnitIntervalMesh(N)

It can be a good idea to check the mesh properties (especially in higher dimensions, or with hierarchically constructed meshes):

In [5]:
print(f"mesh1D has {mesh1D.num_vertices()} vertices and {mesh1D.num_cells()} edges.")
mesh1D has 17 vertices and 16 edges.

The FEM space¶

We define a function space on our mesh. The basic syntax is
V = FunctionSpace( <name of mesh> , <element family>, <element degree>)

In the following example we'll use Continuous Galerkin ("CG" = "Lagrange") elements of degree 1.

In [6]:
V = FunctionSpace(mesh1D, "CG", 1)

While such spaces are used, primarily, for defining trial and test functions, and functions for storing solution, they can be used to represent any function you like (via interpolation), and preform operations on them (arithmetic, differential, integral, etc).

Here we'll define functions on the space V and use that to represent the true solution. To do that we'll

  • extract the spatial coordinates, x, which has only one member: x[0]
  • express the solution in terms of that.
In [7]:
u_true = Function(V)
x = SpatialCoordinate(V)
u_true.interpolate( u_e(x[0]) )
plot(u_true)
Out[7]:
[<firedrake.pyplot.mpl._FiredrakeFunctionPath at 0x7fb212cf3c50>]

The Weak Form¶

Define trial and test functions:

In [8]:
u = TrialFunction(V)
v = TestFunction(V)

We now define the bilinear and linear forms $$ a(u,v) := \varepsilon^2(u',v') + (ru,v) \qquad and \quad l(v) := (f,v). $$

  • one can represent $u'(x)$ and $v'(x)$ as u.dx(0) and v.dx(0), where here the index says we are differentiating with respect to the first (and only) dimension.
  • dx is also used to represents the (symbolic) integral operator.
In [9]:
a = (epsilon**2)*u.dx(0)*v.dx(0)*dx + r*u*v*dx
l = f(x[0])*v*dx

Boundary conditions¶

It is easy to impose any reasonable boundary conditions you can think of. Here we'll use set homogeneous Dirichlet:

In [10]:
BCs = DirichletBC(V, 0.0, "on_boundary")

Solve¶

And we are ready to solve! Define a function to store the solution, and call the solver:

In [11]:
uN = Function(V)
solve(a == l, uN, bcs=BCs)
plot(uN, label="$u^N$")
plt.legend()
Out[11]:
<matplotlib.legend.Legend at 0x7fb212bab850>

Compute errors¶

We'll define a function $d(x)=u(x)-u^N(x)$ and compute the $L_2$, energy, and $\max$ norms of that.

  • For $L_2$-type norms, apply the dx operator to form an integral, and then assemble that to evaluate.
  • For maximum pointwise error, choose a set of points (e.g., the mesh points) and evaluate d at those, using the at method. (There are other approaches).
In [12]:
d = Function(V).interpolate(uN-u_true)

L2_error = sqrt(assemble( d**2*dx)) # same as L2=
H1_error = sqrt(assemble( d.dx(0)**2*dx))
a_error  = sqrt( float(epsilon)*H1_error**2.0 + L2_error**2.0)

xp = np.linspace(0,1,N+1) # points at which we compute errors
max_pw_error = max(np.abs(d.at(xp)))

print(f"Errors: L2={L2_error : 8.3e}, H1={H1_error : 8.3e}, Energy={a_error : 8.3e}, Max={max_pw_error: 8.3e}")
Errors: L2= 7.218e-03, H1= 7.047e-02, Energy= 2.342e-02, Max= 1.604e-02

Plot solutions¶

While we can plot solutions using the plot() function, which is an interface to matplotlib. But if you are familiar with matplotlib it might be better to use it directly:

In [13]:
x_vals = mesh1D.coordinates.dat.data[:]

plot(u_true, label='$u(x)$') # true
plt.plot(x_vals, uN.at(x_vals),'r--o', label='$u^N$')  # use for small N
plt.xlabel("x")
plt.legend()
Out[13]:
<matplotlib.legend.Legend at 0x7fb2125c9550>

A fitted mesh in 1D¶

We now what to compute solutions (and corresponding errors) on a fitted mesh. At this stage, it makes sense to wrap our solver in a function, so we can call it for different $\varepsilon$ and meshes.

Solver as a function¶

This function also returns both the FEM and true solution for the provided epsilon, on the same function space.

In [14]:
def Solve_1DRD(epsilon, my_mesh, my_degree):
    V = FunctionSpace(my_mesh, "CG", my_degree)
    x = SpatialCoordinate(V)
    u_true = Function(V)
    u_true.interpolate( u_e(x[0]) )
    u = TrialFunction(V)
    v = TestFunction(V)
    a = (epsilon**2)*u.dx(0)*v.dx(0)*dx + r*u*v*dx
    l = f(x[0])*v*dx
    uN = Function(V)
    BCs = DirichletBC(V, 0.0, "on_boundary")
    solve(a == l, uN, bcs=BCs)
    return(uN, u_true)

Let's check if that works. We'll solve for $\varepsilon^2 = 10^{-4}$ and $N=16$; expect oscillations.

In [15]:
N = np.int32(16)
epsilon = 1.0e-2
mesh1D = UnitIntervalMesh(N)
uN,uE = Solve_1DRD(epsilon, mesh1D, 1)
x_vals = mesh1D.coordinates.dat.data[:]
plot(uE, label='$u(x)$') # true
plt.plot(x_vals, uN.at(x_vals),'r--o', label='$u^N$')  # use for small N
plt.legend()
Out[15]:
<matplotlib.legend.Legend at 0x7fb21242c310>

A 1D fitted mesh¶

In this example we'll construct a simple mesh piecewise uniform (Shishkin) mesh, and solve on that. There are various ways to do this. In this example, we construct a uniform mesh, and then change the mesh values.

First define the mesh points as a numpy array.

In [16]:
N=16
p=2 # degree of the element
tau = min( (p+1)*epsilon*np.log(float(N)), 0.25) # mesh transition point
mesh1 = np.linspace(0.0, tau, int(N/4)+1) # left mesh; assuming N divisible by 4
mesh2 = np.linspace(tau, 1-tau, int(N/2)+1) # interior
mesh3 = np.linspace(1-tau, 1.0, int(N/4)+1) # right mesh

Smesh_points = np.concatenate( (mesh1, mesh2[1:-1], mesh3) ) 

Now make a (Firedrake) mesh, and overwrite the mesh points

In [17]:
Smesh1D = UnitIntervalMesh(N)
Smesh1D.coordinates.dat.data[:]=Smesh_points

Let's test:

In [18]:
uN,uE = Solve_1DRD(epsilon, Smesh1D, p)
In [19]:
plot(uE, label='$u(x)$') # true
plt.plot(Smesh_points, uN.at(Smesh_points),'r--o', label='$u^N$')  # use for small N
plt.legend()
Out[19]:
<matplotlib.legend.Legend at 0x7fb2124f8410>

Of course, we'd like to check the errors. We'll define a function to do that:

In [20]:
def ErrorNorms(u1, u2):
    d = Function(V).interpolate(u1-u2)
    L2 = sqrt(assemble( d**2*dx)) # same as L2=
    H1 = sqrt(assemble( d.dx(0)**2*dx))
    Energy  = sqrt( float(epsilon)*H1**2.0 + L2**2.0)
    xp = np.linspace(0,1,N+1) # points at which we compute errors
    Max = max(np.abs(d.at(xp)))
    return (L2, H1, Energy, Max)
In [21]:
L2, H1, Energy, Max = ErrorNorms(uN, uE)
print(f"epsilon={epsilon:0.1e}, N={N:4d}, Errors: L2={L2:.3e}, H1={H1:.3e}, Energy={Energy:.3e}, max={Max:.3e}")
epsilon=1.0e-02, N=  16, Errors: L2=1.120e-04, H1=3.417e-03, Energy=3.596e-04, max=5.256e-04

Convergence tables¶

Many studies of (robust) methods for singularly perturbed problems feature convergence tables showing how the errors depend on both $N$ and $\varepsilon$.

Here is how we can do that (though it is not Firedrake-specific).

First: a function for defining a mesh.

In [22]:
def ShishkinMesh(epsilon, N, p):
    tau = min( (p+1)*epsilon*np.log(float(N)), 0.25) # mesh transition point
    mesh1 = np.linspace(0.0, tau, int(N/4)+1) # left mesh; assuming N divisible by 4
    mesh2 = np.linspace(tau, 1-tau, int(N/2)+1) # interior
    mesh3 = np.linspace(1-tau, 1.0, int(N/4)+1) # right mesh
    Smesh_points = np.concatenate( (mesh1, mesh2[1:-1], mesh3) ) 
    mesh1d = UnitIntervalMesh(N)
    mesh1d.coordinates.dat.data[:]=Smesh_points
    return(mesh1d)

Now we'll solve the problem, and compute errors for

  • $\varepsilon = 10^{-1}, 10^{-2}, \dots, 10^{-6}$
  • $N = 32, 64, \dots, 1024$
In [23]:
p = 1
Epsilons = 10.0**np.arange(-1,-7,-1)
Ns = 2**np.arange(5,11)

L2 = np.zeros( (len(Epsilons), len(Ns)) )
Energy = np.zeros_like(L2)
Max = np.zeros_like(L2)

e=-1
for epsilon in [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6]:
    e+=1
    n=-1
    for N in Ns:
        n+=1
        Smesh1 = ShishkinMesh(epsilon, N, p)
        uN, uE = Solve_1DRD(epsilon, Smesh1, p)
        L2[e,n], H1, Energy[e,n], Max[e,n] = ErrorNorms(uN, uE)

We'll display the table using a pandas data frame. You could also easily make a log-log plot.

In [24]:
pd.set_option('display.float_format', lambda x: '%10.2e' % x)
table = pd.DataFrame(Max, index=Epsilons, columns=Ns)
print('Maximum pointwise errors')
print(table)
Maximum pointwise errors
               32         64         128        256        512        1024
1.00e-01   3.91e-03   9.70e-04   2.42e-04   6.05e-05   1.51e-05   3.78e-06
1.00e-02   1.15e-03   3.61e-04   1.24e-04   4.04e-05   1.28e-05   3.95e-06
1.00e-03   7.67e-04   2.00e-04   4.96e-05   1.24e-05   3.09e-06   7.67e-07
1.00e-04   8.01e-04   2.08e-04   5.18e-05   1.29e-05   3.23e-06   8.07e-07
1.00e-05   8.31e-04   2.08e-04   5.19e-05   1.30e-05   3.25e-06   8.11e-07
1.00e-06   8.34e-04   2.08e-04   5.19e-05   1.30e-05   3.25e-06   8.12e-07

The mandatory log-log plot:

In [25]:
table.T.plot(loglog=True, style='--o')
Out[25]:
<Axes: >

A two-dimensional problem¶

We'll finish with the two-dimensional problem: find $u=u(x,y)$ such that $$-\varepsilon^2 \Delta u + r(x,y)u = f(x,y) \text{ on } (0,1)^2,$$ with homogeneous Dirichlet boundary conditions.

The specific example we'll consider will have $r(x,y) = \sqrt{1+x+y}$ and $f(x,y)=e^{x-y}$.

The mesh¶

We'll solve with quadrilateral elements, on an $N \times N$ tensor product (Shishkin) mesh. Handily, the TensorRectangleMesh method can take two numpy arrays as inputs, and generate the corresponding tensor product mesh.

In [26]:
N = 32
p = 1 # degree of the element

tau = min( (p+1)*epsilon*np.log(float(N)), 0.25) # mesh transition point
mesh1 = np.linspace(0.0, tau, int(N/4)+1) # left mesh; assuming N divisible by 4
mesh2 = np.linspace(tau, 1-tau, int(N/2)+1) # interior
mesh3 = np.linspace(1-tau, 1.0, int(N/4)+1) # right mesh
xp = np.concatenate( (mesh1, mesh2[1:-1], mesh3) )
yp = xp

mesh2D = TensorRectangleMesh(xp, yp, quadrilateral=True)
print(f"mesh2D has {mesh2D.num_vertices()} vertices and {mesh2D.num_cells()} edges.")
mesh2D has 1089 vertices and 1024 edges.

The function space and weak form.¶

Since we are now in 2D, we can use the grad operator and inner() method. (We could have used these in 1D too, if we wished).

In [27]:
V = FunctionSpace(mesh2D, "CG", p)
x,y = SpatialCoordinate(V)
u = TrialFunction(V)
v = TestFunction(V)

a = (epsilon**2)*inner(grad(u), grad(v))*dx + inner(sqrt(1+x+y)*u,v)*dx
l = exp(x+y)*v*dx
uN = Function(V)
BCs = DirichletBC(V, 0.0, "on_boundary")
solve(a == l, uN, bcs=BCs)

Plotting¶

Finally we plot the solution. Plotting is a little limited (I think) and it is better to export for use in some dedicated tool, such as Paraview. But we'll make do:

In [28]:
fig = plt.figure()
axes = fig.add_subplot(projection='3d')
trisurf(uN, axes=axes);

To make an interactive plot, first make a 2D np array of solution values:

In [29]:
[X,Y] = np.meshgrid(xp,yp)
uN_vals = np.zeros_like(X)
for i in range(len(xp)):
    for j in range(len(yp)):
        uN_vals[i,j] = uN.at(xp[i],yp[j])

The use plot_surface from matplotlib: (note: might not work on colab)

In [30]:
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X,Y, uN_vals, cmap='bwr')
fig.colorbar(surf)
ax.set_xlabel('x', labelpad=20)
ax.set_ylabel('y', labelpad=20)
Out[30]:
Text(0.5, 0.5, 'y')

Other things Firedrake can and cannot do¶

  1. Many other elements are implemented; see https://www.firedrakeproject.org/variational-problems.html#supported-finite-elements. These include elemenst on triangles, tets, quads, and hex elements.
  2. dG elements are implemented, but only on conforming meshes.
  3. The 'irksome' module adds support for RK-methods in time (and solver support).
  4. Geometric multigrid
  5. Visualisation is limited (just an interface to matplotlib: use of Paraview recommended).
  6. Nonlinear problems are represented in a natural way, but care is needed with solver options.

The End