Intro to FEniCS - Part 3
Stefano Tronci
In the Intro to FEniCS - Part 2 episode we completed the first step in developing a FEniCS model. We characterised the domain. Then, we selected a governing PDE. We then provided suitable boundary conditions. Finally, we expressed the PDE in weak form. We chosen a generalised Helmholtz equation as our PDE. In this episode we will start using the resulting weak form to develop a model.
Project Files
All the files used for this project are available at the repositories below:
The Model
In this episode we will build the model of a rectangular room. We already built a series of Elmer models on this topic. This is why this system is well suited to explore a new solver: we have some understanding on it. This will help us figuring out whether the model is correct. In the Intro to FEniCS - Part 2 episode we developed a generalised formalism. We now restrict it to a rectangular room. If $l_{x_0}$, $l_{x_1}$ and $l_{x_2}$ are the dimension of the room along $x_0$, $x_1$ and $x_2$ respectively we have:
$$ l_{x_0} \doteq 4 \space \text{m} \\ l_{x_1} \doteq 5 \space \text{m} \\ l_{x_2} \doteq 3 \space \text{m} $$ $$ \Rightarrow \Omega \doteq \left( 0, l_{x_0} \right) \times \left( 0, l_{x_1} \right) \times \left( 0, l_{x_2} \right) $$
To make things easy, as this is our first FEniCS model, we will do some semplifications.
First, we remove all damping. So:
$$ \eta = 0 \space \frac{\text{rad}}{\text{s}} \text{m}^2 \\ \alpha = 0 \space \frac{\text{rad}}{\text{s}} $$
$$ \Rightarrow \gamma = 1, \space \beta = k^2; \space k = \frac{2 \pi \nu}{c} $$
Then, we only consider Neumann boundary conditions. The room wall identified by the equation $x_1=0$ will have an uniform velocity $w$ of $10$ $\frac{\text{m}}{\text{s}}$. All the other walls will be rigid. Hence, they will have a normal uniform velocity of $0$ $\frac{\text{m}}{\text{s}}$. Hence:
$$ N_{N} = 6 $$
If $\mathbf{x}$ is, again, the generic point in $\overline{\Omega}$:
$$ \mathbf{x} = \begin{bmatrix} x_{0} \\ x_{1} \\ x_{2} \end{bmatrix} \in \overline{\Omega}$$
$$ \begin{gather*} \\ \partial \Omega _{N} ^{0} \doteq \left\{ \mathbf{x} \in \overline{\Omega} : x_0 = 0 \right\} \\ \partial \Omega _{N} ^{1} \doteq \left\{ \mathbf{x} \in \overline{\Omega} : x_0 = l_{x_0} \right\} \\ \partial \Omega _{N} ^{2} \doteq \left\{ \mathbf{x} \in \overline{\Omega} : x_1 = 0 \right\} \\ \partial \Omega _{N} ^{3} \doteq \left\{ \mathbf{x} \in \overline{\Omega} : x_1 = l_{x_1} \right\} \\ \partial \Omega _{N} ^{4} \doteq \left\{ \mathbf{x} \in \overline{\Omega} : x_2 = 0 \right\} \\ \partial \Omega _{N} ^{5} \doteq \left\{ \mathbf{x} \in \overline{\Omega} : x_2 = l_{x_2} \right\} \end{gather*} $$
The normal boundary velocities for each $\partial\Omega_{N}^{l}$, for $l=0,1,2,3,4,5$, are:
$$ w_l \doteq \begin{cases} 10 \space \frac{\text{m}}{\text{s}} & \text{if } \space l = 2 \\ 0 \space \frac{\text{m}}{\text{s}} & \text{otherwise} \end{cases} $$
Finally, we consider the acoustic field as awoken by the Neumann boundary condition only. So, we put the source term $q$ to $0$:
$$ q = 0 $$
With all the choices above the weak form reduces to (after changing signs):
$$ \begin{aligned}
a\left(u,\phi\right) = \int_{\Omega} \nabla u \cdot \nabla \phi^\star d\mathbf{x} - k^2 \int_{\Omega} u \phi^{\star} d\mathbf{x} = \\
= \int_{\partial \Omega_{N}^2} g_{2} \phi^\star ds = L\left(\phi\right)
\end{aligned} $$
With $g$ being the acoustic flux for the boundary $\partial\Omega_{N}^{2}$:
$$ g_2 \doteq j \omega \rho w_2 = j \omega \rho 10 $$
As we noticed at the end of the Intro to FEniCS - Part 2 episode we need to expand the integrands in real and imaginary parts. This will make most sense when we start writing the code. For the time being, let’s just prepare the ingredients:
$$ \begin{gather*}
u = \Re \left( u \right) + j \Im \left( u \right) \\
\phi = \Re \left( \phi \right) + j \Im \left( \phi \right) \\
g_2 = \Re \left( g_2 \right) + j \Im \left( g_2 \right) \\
\Rightarrow \nabla u = \Re\left(\nabla u \right) + j \Im\left(\nabla u \right) = \nabla \Re \left(u\right) + j \nabla \Im \left(u\right) \\
\Rightarrow \nabla \phi = \Re\left(\nabla \phi \right) + j \Im\left(\nabla \phi \right) = \nabla \Re \left(\phi\right) + j \nabla \Im \left(\phi\right)
\end{gather*}$$
Hence, we are ready to rewrite the integrals. The first inegral for $a\left(u,\phi\right)$ reads:
$$ \begin{aligned}
\int_{\Omega} \nabla u \cdot \nabla \phi^\star d\mathbf{x} = \\
= \int_{\Omega} \nabla\Re\left( u \right) \cdot \nabla \Re\left( \phi \right) + \nabla\Im \left( u \right) \cdot \nabla\Im \left( \phi \right) d\mathbf{x} + \\
+ j \int_{\Omega} \nabla\Im \left( u \right) \cdot \nabla\Re\left( \phi \right) - \nabla\Re\left( u \right) \cdot \nabla\Im \left( \phi \right) d\mathbf{x} = \\
= \int_{\Omega} \nabla\Re\left( u \right) \cdot \nabla \Re\left( \phi \right) d\mathbf{x} + \int_{\Omega} \nabla\Im \left( u \right) \cdot \nabla\Im \left( \phi \right) d\mathbf{x} + \\
+ j\left( \int_{\Omega} \nabla\Im \left( u \right) \cdot \nabla\Re\left( \phi \right) d\mathbf{x} - \int_{\Omega} \nabla\Re\left( u \right) \cdot \nabla\Im \left( \phi \right) d\mathbf{x} \right)
\end{aligned} $$
The second integral for $a\left(u,\phi\right)$ reads:
$$ \begin{aligned}
\int_{\Omega} u \phi^{\star} d\mathbf{x} = \\
= \int_{\Omega} \Re\left(u\right) \Re\left(\phi\right) + \Im\left(u\right) \Im\left(\phi\right) d\mathbf{x} + \\
+ j \int_{\Omega} \Im\left(u\right) \Re\left(\phi\right) - \Re\left(u\right) \Im\left(\phi\right) d\mathbf{x} = \\
= \int_{\Omega} \Re\left(u\right) \Re\left(\phi\right) d\mathbf{x} + \int_{\Omega} \Im\left(u\right) \Im\left(\phi\right) d\mathbf{x} + \\
+ j \left( \int_{\Omega} \Im\left(u\right) \Re\left(\phi\right) d\mathbf{x} - \int_{\Omega} \Re\left(u\right) \Im\left(\phi\right) d\mathbf{x} \right)
\end{aligned} $$
While $L\left(\phi\right)$ reads:
$$ \begin{aligned}
\int_{\partial \Omega_{N}^2} g_{2} \phi^\star ds = \\
= \int_{\partial \Omega_{N}^2} \Re\left(g_{2}\right) \Re\left(\phi\right) + \Im\left(g_{2}\right) \Im\left(\phi\right) ds + \\
+ j \int_{\partial \Omega_{N}^2} \Im\left(g_{2}\right) \Re\left(\phi\right) - \Re\left(g_{2}\right) \Im\left(\phi\right) ds = \\
= \int_{\partial \Omega_{N}^2} \Re\left(g_{2}\right) \Re\left(\phi\right) ds + \int_{\partial \Omega_{N}^2} \Im\left(g_{2}\right) \Im\left(\phi\right) ds + \\
+ j \left( \int_{\partial \Omega_{N}^2} \Im\left(g_{2}\right) \Re\left(\phi\right) ds - \int_{\partial \Omega_{N}^2} \Re\left(g_{2}\right) \Im\left(\phi\right) ds \right)
\end{aligned} $$
Now that all integrals are split in real and imaginary parts we have no trouble putting together the real and imaginary parts for $a \left(u,\phi\right) $ and $L\left(\phi\right)$. We will do this last step right in the code. For now note how each real part and imaginary part is the addition of two terms.
Implementation
We will refer to the implementation as in the repository.
In the following we will assume you followed the installation instructions in the Intro to FEniCS - Part 1 episode
Preparation
First, let’s create a folder in your favorite location. Let’s cd into it and create a Python virtual environment:
mkdir project && cd project
python3 -m venv --system-site-packages venv
the --system-site-packages ensures our new virtual environment has FEniCS in it.
Then, you can open PyCharm. You want to open an existing project and select our project folder created above. PyCharm could launch in different ways. If it opens on the editor select File > Open. For more information, refer to the PyCharm Documentation.
However you end up opening the project folder with PyCharm you should have it loaded correctly with the Python interpreter from our new virtual environment. If not, check the settings under File > Settings > Project: project > Python Interpreter and File > Settings > Build, Execution, Deployment > Console > Python Console. As an example, see the pictures below.
Figure 1
Example PyCharm Interpreter Settings.
Figure 2
Example PyCharm Console Settings.
Note that it is entirely possible to create e new project all from the PyCharm GUI. Make sure you tick Inherit global site-packages when you do so. Creating the environment beforehand is typically easier when cloning projects, although also this can be done with the PyCharm GUI.
In the PyCharm window you will see, on the left, a browser. Right click project and select New > Python Package as shown below.
Figure 3
New Package Creation.
We can call this package simulation. You will see this crates a directory under project called simulation. In it you will find an empty __init__.py file. Right click the simulation folder in the browser and choose New > Python File. Let’s call the new file helmholtz. We will implement the solver in this file.
Note that it is not necessary to create a Python package to use FEniCS. We can simply code a script. However, it is useful to make a package as it makes things nice and modular and simple to reuse.
Code
The entire code is available here
First, let’s import the packages we need:
import numpy as np
import fenics
import pathlib
Then, we wrap all our simulation code in a function:
def simulate(
real_part_output_path: pathlib.Path,
imag_part_output_path: pathlib.Path
) -> None:
This function takes as arguments the paths to the real part and imaginary part files. That is, this code will produce two files, one for the real part of the field and one for the imaginary part. It will save them at the specified locations. These locations are provided as Path objects from pathlib. This is not necessary, but pathlib actually makes working with files very easy and convenient, and it is then recommended. We will use numpy for a few mathematical functions.
Now, let’s put all of our project parameters in a few variables:
nu = 57.17 # Frequency of the Simulation, Hz
c = 343 # Speed of sound in air, m/s
rho = 1.205 # Density of air, kg/m^3
w = 10 # Velocity normal to boundary, for inflow Neumann condition, m/s
l_x0 = 4 # Room size along x0, m
l_x1 = 5 # Room size along x1, m
l_x2 = 3 # Room size along x2, m
tol = 1e-10 # Tolerance for boundary condition definitions
The tol parameter will become clearer later. Note how these paramters could all be arguments of the function. That is, you can make a function to simulate any room. In this code we kept them in the body for simplicity.
We now compute the mesh size. It should be smaller than one tenth of the wavelenght. Once we have the size we can easily compute the number of element per side of our room.
# Computing some useful stuff
omega = 2 * np.pi * nu # Angular frequency, rad/s
k = omega / c # Wave number, rad/m
s = c / (10 * nu) # Element Size
n_x0 = np.intp(np.ceil(l_x0 / s)) # Number of elements for each direction
n_x1 = np.intp(np.ceil(l_x1 / s))
n_x2 = np.intp(np.ceil(l_x2 / s))
We use np.ceil to make sure that the number is the smallest integer bigger or equal to the exact size. np.intp is used to convert the type to integer (np.ceil will still return a floating point value).
Now it comes an important step. We are using FEniCS. But FEniCS does not support complex PDEs. But our PDE is complex. Turns out that a complex PDE is actually two coupled PDEs, one for the real part and one for the imaginary part. This is why we expanded the weak form in real and imaginary parts, so that we can easily see them. Solving coupled PDEs can be done in FEniCS by using Mixed Spaces. The Solving PDEs in Python book covers this in details in the chapter 3.5. For more information refer to said chapter. In a nutshell, we will be solving two coupled PDEs, one for the real part and one for the imaginary one. This means that $j$, the imaginary unit, will disappear.
First, we crate the mesh:
# Making a mesh, then an element. Then, Mixed Space so that we can solve for real and imaginary parts.
mesh = fenics.BoxMesh(
fenics.Point(0, 0, 0),
fenics.Point(l_x0, l_x1, l_x2),
n_x0,
n_x1,
n_x2
)
The BoxMesh function creates a box with one corner at the origin, the opposite corner at the chosen lengths. It also meshes it with the number of elements we previously computed. The mesh is tetrahedral.
Then, we define the type of element. We use a second order Lagrange element:
P = fenics.FiniteElement("Lagrange", mesh.ufl_cell(), 2)
This means that each cell in the mesh will be a second order element.
Finally, we can define our vector space V as a mixed space over the mesh:
V = fenics.FunctionSpace(mesh, P * P)
The mixed aspect is given by P * P, which signifies a mixed element for two fields, each of the P kind. Note that V is not the vector space $V$ we defined in the Intro to FEniCS - Part 2 episode. The space V is the space of functions represented as linear combinations of basis function (second order polynomials) over the mesh. This space is finite dimensional (finite number of basis functions). This is what allows to compute an approximate solution. That is: FEM computes a numerical solution by simplifying $V$. This step is automated by FEniCS, as you see. V largely works as $V$ in the sense that it is a linear space in which we search for a solution to the (now) simplified problem. So we still refer to it with the same letter. For more information refer to the FEM in a Nutshell episode.
In the V space the unknown function we are searching for is called u. This also goes by the name of Trial Function. The test function (which in $V$ was called $\phi$) will be called phi. Hence:
u_re, u_im = fenics.TrialFunctions(V)
phi_re, phi_im = fenics.TestFunctions(V)
Since we are using a mixed space these functions are returned already split in real and imaginary parts, as we did in the maths above. For convenience we also define $k^2$ as a Constant:
k_sq = fenics.Constant(k**2)
Now, we have all the ingredients ready for $a\left(u,\phi\right)$. Since we are using mixed space we simply sum the terms. That is, in the maths above we have split all terms that constitute $a\left(u,\phi\right)$ in real and imaginary parts. In mixed space formulation is like we had one $a\left(u,\phi\right)$ for the real part and one $a\left(u,\phi\right)$ for the imaginary part. These two are constituted of the real part and imaginary part of the original complex $a\left(u,\phi\right)$. The mixed space a is then simply the sum of them. We report it below:
a = \
fenics.inner(fenics.nabla_grad(u_re), fenics.nabla_grad(phi_re)) * fenics.dx - k_sq * phi_re * u_re * fenics.dx + \
fenics.inner(fenics.nabla_grad(u_im), fenics.nabla_grad(phi_im)) * fenics.dx - k_sq * phi_im * u_im * fenics.dx + \
fenics.inner(fenics.nabla_grad(u_im), fenics.nabla_grad(phi_re)) * fenics.dx - k_sq * u_im * phi_re * fenics.dx - \
fenics.inner(fenics.nabla_grad(u_re), fenics.nabla_grad(phi_im)) * fenics.dx + k_sq * u_re * phi_im * fenics.dx
You will remember how each real and imaginary part of the integrals for $a\left(u,\phi\right)$ is composed of two terms. The first line for a above puts together the first terms, the second one the second terms… and so on. a does not need to be constructed in this specific order, but this form is clearer. It shows that a mixed space a has a term with real parts only, one with an imaginary part only, and 2 mixed terms. The mixed terms are important, as they dictate the relationship the imaginary and real parts have: they are not independent. Note that a is defined without any symbol for integral. FEniCS already know it will need to take the integral of this thing, so we do not need it. We have to specify with respect what we have to integrate. Since these are volume integrals we use the default volume measure fenics.dx. One can define custom ones, for example if some integral are within a sub-volume. We will do something similar right now for the boundaries. Before leaving this section, note the use of inner and nabla_grad. inner is used to implement the scalar product between the gradients. As these are 1D tensors (vectors) then the result is the same as dot. Similarly, since u_re and u_im are scalar fields, using nabla or nabla_grad is not really different. However, for vector fields the choice between inner and dot, nabla and nabla_grad becomes important. These differences are explained at page 25 and page 58 of Solving PDEs in Python.
All is left to do is to specify the boundary conditions. One very useful resource (together with Solving PDEs in Python) is Douglas N Arnold page. What we want to do is to assign a different Neumann boundary condition to each wall (although only one wall has a non-zero normal velocity). So, we need to split the boundary. Since we are already dealing with the mesh we need to mark the nodes that belong to the each boundary. Then, on each boundary, we define a custom measure so that when we define $L\left(\phi\right)$ each integral is restricted to the correct boundary. In other words: since FEniCS does not have a symbol for the integral, but one for the measure (dx, ds…) we need the measure itself to be associated to the subset to which it is applied. Let’s go through.
First, we make a MeshFunction that returns integers. We initialise it to return 0 everywhere.
mf = fenics.MeshFunction("size_t", mesh, 2)
mf.set_all(0)
This function returns an unsigned integer of type size_t (the largest unsigned integer type on your machine). Note how size_t is a C\C++ type. This because all of these functions act on the C\C++ backend. Of course, this mesh function is defined over our mesh mesh which we defined above. Since we are using mf to mark boundaries, which are 2D entities, its dimensionality is 2. We will now edit it to return a different number for each boundary. This way, when our mesh function takes in a mesh node, it will return a number that tells FEniCS on which boundary that is (if any).
To do so, we need to define the boundaries. A boundary is a domain. So we use SubDomain. Defining a SubDomain is easy. We simply create a class which has SubDomain as a parent and override the inside method. This is the method that FEniCS will use to figure out if a mesh point is inside our subdomain. Below are all the subdomain for the various room walls:
class BX0(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[0], 0, tol)
class BXL(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[0], l_x0, tol)
class BY0(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[1], 0, tol)
class BYL(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[1], l_x1, tol)
class BZ0(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[2], 0, tol)
class BZL(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[2], l_x2, tol)
Note how the inside method takes the mesh point x, a 3D vector with the Cartesian coordinates of the point, and a boolean variable on_boundary as input. That is, in the method inside we will already know if x is on the boundary or not based on the value of on_boundary. We simply need to figure out which of our sub-boundaries that is. To do so, we use the function near. Since everything is in floating point, we cannot check for equality, so we need a tolerance. This is what our tol variable is for. For more information about this, check the excellent LEARN C++ tutorial about floating point comparison. As you can see, the SubDomain definitions directly mirror those for $\partial \Omega _{n}^l$ that we provided above.
Now that we have the sub-domains we just need to instruct our mesh function to return a different integer on each of them. We first instantiate one object for each SubDomain type we crated and use mf to mark the nodes with an integer of our choice.
bx0 = BX0()
bx0.mark(mf, 1)
bxl = BXL()
bxl.mark(mf, 2)
by0 = BY0()
by0.mark(mf, 3)
byl = BYL()
byl.mark(mf, 4)
bz0 = BZ0()
bz0.mark(mf, 5)
bzl = BZL()
bzl.mark(mf, 6)
Now, we define the fluxes $g_l$. Or better, we define their “digital” counterpart. We only have dead boundaries (zero velocity) and alive boundaries (non-zero velocity). So:
# The flux is 0 for all rigid walls
flux_rig = 0
g_rig_re = fenics.Constant(np.real(flux_rig))
g_rig_im = fenics.Constant(np.imag(flux_rig))
# The flux is this for the active walls
flux_in = 1j * omega * rho * w
g_in_re = fenics.Constant(np.real(flux_in))
g_in_im = fenics.Constant(np.imag(flux_in))
Finally, we define our custom measure ds over our mesh. We use our mesh function mf to ensure that ds has the various measures for each of the sub-domains:
ds = fenics.Measure('ds', domain=mesh, subdomain_data=mf)
Now, if we use ds(1) we will access the measure for boundary marked with 1 by mf, with ds(2) that marked with 2 and so forth. Hence, we can now put together L, the FEM analogue of $L\left(\phi\right)$:
L = \
(g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(1) + \
(g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(2) + \
(g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(6) + \
(g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(4) + \
(g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(5) + \
(g_in_re * v_re + g_in_im * v_im + g_in_im * v_re - g_in_re * v_im) * ds(3)
Since g_rig_re and g_rig_im are simply 0 we could simply omit the terms with them. In that case we would be left with the last term only, as our equations from above tells. However, it is instructive to show the fully expanded form for L, with all the different terms for each boundary.
All we have left to do is solve and save.
To solve, we define a new Function over the vector space V and invoke the solve method:
u = fenics.Function(V)
fenics.solve(a == L, u)
This will solve the equation and put the result in u. This is the easiest way to solve, but more advanced ways are also available. Control over the solver algorithm is provided by the more advanced methods. For now, let’s stick to the simple stuff.
Since u is a mixed space solution we need to split it to access real and imaginary parts. We do so and rename as well for easy access in the final VTK files. VTK files are simple to create with the File class. The code is below:
u_1, u_2 = u.split()
u_1.rename('re', 're')
u_2.rename('im', 'im')
vtk_file_1 = fenics.File(real_part_output_path.absolute().as_posix())
vtk_file_1 << u_1
vtk_file_2 = fenics.File(imag_part_output_path.absolute().as_posix())
vtk_file_2 << u_2
Running the code is simple. From your python environment simply run:
from simulation import helmholtz
import pathlib
path_re = pathlib.Path().joinpath('re.pvd')
path_im = pathlib.Path().joinpath('im.pvd')
helmholtz.run(path_re, path_im)
This will crate the files re.pvd and im.pvd in your current directory. Or, you can create a run.py file in your project folder with these contents. Then, simply run the file from your project folder:
venv/bin/python run.py
Results
We computed the results and solved the same problem with Elmer as well. We will not discuss the Elmer project setup as there are plenty of tutorials on this website about that. You can refer to the Elmer project files for more details. The results are comapred below.
Figure 4
Results (Elmer VS FEniCS).
Results are overall very similar. The real part is negligible, as expected from a damping free material with pure imaginary boundary flux. In fact, the weak form could be much simplified in this special case to yield a faster solution. However, we shown the entirety of it as it is more instructive. Interestingly the polarity of the real part is flipped in FEniCS. The imaginary parts agree, overall, fairly well. The agreement is not exact. However, since we validated Elmer before, we know that our solution with FEniCS is already a good starting point.
Conclusion
We solved a simplified equation for a damping free material in a rectangular room. We derived the equations to do so and written the FEniCS code to compute the solution. This solution was computed also with Elmer for benchmark. Results are very similar, a sign that our FEniCS solution was setup properly.
Getting to solve PDEs in FEniCS is more work with respect doing so with Elmer. However, we can setup to solve any PDE in a programmatic way. This means that we can write more complex simulation and analysis code all within Python. Additionally, getting intimate with the weak form provides more insight about the underlying formalism at the base of PDEs.
In the next episodes we will explore solving PDEs with FEniCS further.
License Information
This work is licensed under a Creative Commons Attribution 4.0 International License.