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>]
21st Workshop on Numerical Methods for Problems with Layer Phenomena
Niall Madden and Sean Tobin, April 2025
You can access this notebook at:
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:
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:
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).
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:
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:
import numpy as np
import pandas as pd
# %matplotlib notebook # may not work on colab
import matplotlib.pyplot as plt
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.
I like to use lambda functions for defining functions
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)
We'll use a uniform mesh with $N=16$ intervals.
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):
print(f"mesh1D has {mesh1D.num_vertices()} vertices and {mesh1D.num_cells()} edges.")
mesh1D has 17 vertices and 16 edges.
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.
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
x, which has only one member: x[0]u_true = Function(V)
x = SpatialCoordinate(V)
u_true.interpolate( u_e(x[0]) )
plot(u_true)
[<firedrake.pyplot.mpl._FiredrakeFunctionPath at 0x7fb212cf3c50>]
Define trial and test functions:
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). $$
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.a = (epsilon**2)*u.dx(0)*v.dx(0)*dx + r*u*v*dx
l = f(x[0])*v*dx
It is easy to impose any reasonable boundary conditions you can think of. Here we'll use set homogeneous Dirichlet:
BCs = DirichletBC(V, 0.0, "on_boundary")
And we are ready to solve! Define a function to store the solution, and call the solver:
uN = Function(V)
solve(a == l, uN, bcs=BCs)
plot(uN, label="$u^N$")
plt.legend()
<matplotlib.legend.Legend at 0x7fb212bab850>
We'll define a function $d(x)=u(x)-u^N(x)$ and compute the $L_2$, energy, and $\max$ norms of that.
dx operator to form an integral, and then assemble that to evaluate.d at those, using the at method. (There are other approaches).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
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:
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()
<matplotlib.legend.Legend at 0x7fb2125c9550>
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.
This function also returns both the FEM and true solution for the provided epsilon, on the same function space.
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.
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()
<matplotlib.legend.Legend at 0x7fb21242c310>
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.
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
Smesh1D = UnitIntervalMesh(N)
Smesh1D.coordinates.dat.data[:]=Smesh_points
Let's test:
uN,uE = Solve_1DRD(epsilon, Smesh1D, p)
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()
<matplotlib.legend.Legend at 0x7fb2124f8410>
Of course, we'd like to check the errors. We'll define a function to do that:
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)
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
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.
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
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.
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:
table.T.plot(loglog=True, style='--o')
<Axes: >
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}$.
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.
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.
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).
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)
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:
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:
[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)
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)
Text(0.5, 0.5, 'y')
matplotlib: use of Paraview recommended).The End