About FEniCS
FEniCS is the project name for a set of computational components, which include libraries such as DOLFIN, UFL (unified form language, which is a special language for implementing the variational formulation of a problem), FFC (FEniCS form compiler), UFC (unified formassembly code) and FIAT (finite element automatic tabulator). You can find the opensource repositories here.
In a nutshell, FFC transforms UFL language into the corresponding C++ interface code (UFC). FIAT contains finite element basis functions and the corresponding quadrature rules. DOLFIN contains the data structures, assembly routines and linear solvers.
The next diagram, although slightly outdated (e.g. FErari is deprecated), gives an excellent overview of the ecosystem.
The normal user interacts directly only with DOLFIN and UFL. DOLFIN is written is C++, but wrapped by Python. UFL is pure Python code.
DOLFIN interfaces with several linear algebra packages, known in FEniCS lingo as linear algebra backends. By default, PETSc is used. PETSc supports MPI, GPUs through CUDA or OpenCL and hybrid MPIGPU parallelism (DOLFIN only support MPI). This reliance on external solvers should not be overlooked: it allows the use of stateoftheart solvers with minimal effort, while being extremely flexible. If you are brave enough and really want to use your own linear algebra library, you can probably change the source code to interact with it (that’s the beauty of opensource software).
To check which are the backends available in your installation, just do:
from dolfin.cpp.la import list_linear_algebra_backends
list_linear_algebra_backends()
You can also get more information about the solvers and preconditioners by doing:
from dolfin.cpp.la import list_lu_solver_methods
from dolfin.cpp.la import list_krylov_solver_methods
from dolfin.cpp.la import list_krylov_solver_preconditioners
list_lu_solver_methods()
list_krylov_solver_methods()
list_krylov_solver_preconditioners()
This page shows how to install FEniCS in your machine. The easiest way for a Python regular user is to rely on Anaconda (Linux and MacOS only):
conda install fenics
Recently, FEniCS core libraries started undergoing a major redevelopment. That’s why you may have found information about FEniCSX in your searches. These new libraries should be used in the future, as the old project will “see much less development and less maintenance”. Nevertheless, since the new project was officially released only very recently, this tutorial still uses the old libraries. You can find out everything about it here. The new libraries repositories can be found here.
Our project
In this tutorial, the goal is to show how you can take an existing Cerfacs mesh (created with e.g. hip
) and solve the timedependent heat equation problem. It shows how to:
 convert the mesh in a DOLFIN’sfriendly format
 apply boundary and initial conditions
 translate the variational formulation into code
 solve the problem relying on different solvers
 parallelize your code
The heat equation tutorial found in Langtangen’s 2017 book (web version) is closely followed. You should go there if you struggle with any of the steps.
Formulating the problem
TODO: in time also TODO: add formulation
After we have the variational problem formulated, implementing it in FEniCS is straightforward.
A first working simulation
Convert and load the mesh
We can create simple meshes inside FEniCS, such as unit squares (dolfin.cpp.generation.UnitSquareMesh
) and 3D boxes (dolfin.cpp.generation.BoxMesh
). Nevertheless, nonacademic applications normally require more complex geometries. For such cases, you are advised to use your own favorite CAD tool and mesh generator and rely on mesh conversion tools (such as pyhip
, meshio
and/or dolfinconvert
) to convert it to the desired format (Dolfin’s .xml
).
In this example, we start with a typical Cerfacs .hdf5
file. Then, we convert it to gmsh
format using pyhip
and finally meshio
helps us converting it to a DOLFIN’s friendly format.
We start by converting from .hdf5
to .msh
:
import os
from pyhip.commands.readers import read_hdf5_mesh
from pyhip.commands.writers import write_gmsh
from pyhip.commands.operations import hip_exit
mesh_dir = os.path.join('mesh')
mesh_filename = os.path.join(mesh_dir, 'trappedvtx.mesh.h5')
gmsh_filename = os.path.join(mesh_dir, 'trappedvtx')
read_hdf5_mesh(mesh_filename)
write_gmsh(gmsh_filename)
hip_exit()
Then, we convert from .msh
to .xml
(meshio
infers the new file format from the file name):
import os
import meshio
mesh_dir = os.path.join('mesh')
mesh_filename = os.path.join(mesh_dir, 'trappedvtx.mesh.msh')
meshio_mesh = meshio.Mesh.read(mesh_filename)
new_mesh_filename = os.path.join(mesh_dir, 'trappedvtx.xml')
meshio_mesh.write(new_mesh_filename)
And finally, we load the mesh in DOLFIN:
import os
from dolfin.cpp.mesh import Mesh
mesh_filename = os.path.join('mesh', 'trappedvtx.xml')
mesh = Mesh(mesh_filename)
If your using a notebook, we can be visualize the imported mesh by simply typing mesh
.
This mesh has 34684 nodes.
Define the function space
TODO: add this info to the problem formulation
As seen in the problem formulation, a variational problem consists in finding the solution that lies in the trial space. To define a function space (which will be then used to build test and a trial spaces), you just need to specify the element family (look here for more options). For a simple Lagrange space of degree 1:
from dolfin.function.functionspace import FunctionSpace
V = FunctionSpace(mesh, 'Lagrange', 1)
Define the boundary conditions
In this case, we specify Dirichlet boundary conditions at the inlet and the outlet. All the other boundaries assume null flux.
from dolfin.function.constant import Constant
from dolfin.cpp.math import near
from dolfin.fem.dirichletbc import DirichletBC
x_min, x_max = 0., 0.2
u_xmin, u_xmax = Constant(300.), Constant(300.)
def boundary_xmin(x, on_boundary):
return on_boundary and near(x[0], x_min)
def boundary_xmax(x, on_boundary):
return on_boundary and near(x[0], x_max)
bcs = [DirichletBC(V, u_xmin, boundary_xmin),
DirichletBC(V, u_xmax, boundary_xmax)]
Here some aspects worth noticing:
 The boolean
on_boundary
checks if a node belongs to a boundary (it is optional).  Use
near
to avoid roundoff errors. dolfin
passes the node coordinates as first argument to the functions created to verify if a node belongs to a given boundary (boundary_xmin
andboundary_xmax
)Constant
specifies a constant value in all the nodes at the corresponding boundary. We could have use afloat
instead. Nevertheless, this way gives you already a hint on how to set up nonconstant boundary conditions: useExpression
instead. All the boundary conditions are stored in an array that will then be passed to the solver.
Define the initial condition
Since this is a timedependent problem, we have to specify an initial condition.
from dolfin.function.expression import Expression
from dolfin.function.function import Function
c = (x_max  x_min) / 2
param = u_xmin.values()[0] / c**2
u_initial = Expression('param * (x[0]  c) * (x[0]  c)',
degree=2, param=param, c=c, name='T')
u_n = Function(V, name='T')
u_n.interpolate(u_initial)
A small caveat: string expressions must have valid C++ syntax.
u_n
represents our initial condition (a parabolic temperature field). We are getting it via interpolation. Alternatively, this could have been solved by projection (by turning the initial value problem into a variational problem), via projection
method.
Translate the problem formulation into code
The last setup step is to “translate” the variational mathematical formulation into code. The unified form language (UFL) library is of great help here.
from dolfin.function.argument import TestFunction
from dolfin.function.argument import TrialFunction
from dolfin.function.function import Function
from dolfin.function.constant import Constant
from ufl import grad
from ufl import dx
from ufl import dot
u_h = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.)
a = u_h*v*dx + dt*dot(grad(u_h), grad(v))*dx
L = (u_n + dt*f)*v*dx
Don’t forget to define dt
first.
T = 0.05 # final time
num_steps = 20 # number of time steps
dt = T / num_steps
TODO: comment about choice of u_h as variable name
Solve the problem
And finally we arrived to the last destination. Nevertheless, let’s not go there so fast… let’s start simple. If our problem did not depend on time, this is how we would solve it:
from dolfin.function.function import Function
from dolfin.fem.solving import solve
u = Function(V, name='T')
solve(a == L, u, bcs)
Simple, isn’t it? You can also store the results in a format amenable to paraview
:
import os
from dolfin.cpp.io import File
vtkfile = File(os.path.join('outputs', 'trappedvtx_heat_eq.pvd'))
vtkfile << u
Nevertheless, our problem does depend on time. This means we will have to perform a loop and update the timedependent variables each cycle:
from dolfin.cpp.io import File
from dolfin.function.function import Function
solver_parameters = {}
vtkfile = File(os.path.join('outputs', 'trappedvtx_heat_eq.pvd'))
# store initial solution
t = 0
vtkfile << (u_n, t)
u = Function(V, name='T')
for n in range(num_steps):
# update current time
t += dt
# compute solution
solve(a == L, u, bcs, solver_parameters=solver_parameters)
# save to vtk
vtkfile << (u, t)
# update previous solutions
u_n.assign(u)
And this is all. You can now go to paraview
and play with your results.
TODO: better version of the gif (it stays to much time without changing)
What about a fastest solver?
But wait… have you followed all the steps and now you are very unhappy with the performance of the solver? Are you wondering why such a slow code would be worth too look at? Well, the slowness is your fault… you are solving a mediumsize problem with LU decomposition, which is the default solver in FEniCS and which is suggested for problems with “a few thousand of unknowns”. Just do the following, and your doubts will start vanishing.
solver_parameters = {'linear_solver': 'gmres',
'preconditioner': 'ilu'}
We have just changed the solver to a Krylov subspace method and our performance increased about tenfold. In the implementation, only the line containing the solver changes, since now we want to have control over its parameters.
Let’s go parallel
DOLFIN supports parallel computing through MPI and is supposed to perform well in massively parallel supercomputers. It is designed such that the parallel simulations use the same code as for serial computations. The difference is the way you run your script:
mpirun n <n_cpus> python <your_script>
Don’t forget to use a linear algebra backend that supports parallel operations (e.g. PETSc
).
Multithreading was previously supported, but it was removed around 2014.
TODO: speak about mesh partitioning
Benchmarking against AVTP
TODO: complete case description: outlet Dirichlet 1300K, everything else 300K, fixed timestep
TODO: show how to run FEniCS within kraken
TODO: speak about conda problems + fenics problems when too little dof per cpu
For a coarse mesh (~150k cells):
For a fine mesh (~10M cells):
Useful tips
Here, a list of useful tips that, hopefully, will make your life easier.

Some functions don’t have the expected behavior within a
jupyter notebook
. This applies mainly to functions imported fromdolfin.cpp
which task is to print information (e.g.dolfin.cpp.la.list_linear_algebra_backends
). This does not mean you should stop usingjupyter notebook
(they’re very useful!), it simply means you have to go for a script to get this kind of information. 
You can change the verbosity of your code:
set_log_level(log_level)
. In old documentation, you can findset_log_level(DEBUG)
andset_log_level(PROGRESS)
as the goto commands. Nevertheless, they do not work in version2019.1.0
, asDEBUG
andPROGRESS
cannot be imported. To access this variables you have to importLogLevel
(from dolfin.cpp.log import LogLevel
) and then access the desired variables (e.g.LogLevel.DEBUG
). Alternatively , you can use an integer, beinglog_level=10
andlog_level=16
, forDEBUG
andPROGRESS
, respectively. Other options can be found here. I could tell you what each level prints, but why wouldn’t try it yourself? (not in ajupyter notebook
!)
TODO: pydoc
TODO: linear variational problem and solver objects; speak about the source code
Useful references
TODO: add fenics (e.g. forum) web books papers include Barbacourses
Closing remarks
There are several ways to implement a problem in FEniCS. The example presented here follows a scriptbased way, which is easier to read in a blog post. Nevertheless, if you start using FEniCS often, you will crave for creating your own objects that encapsulate the tasks you perform recurrently (and which are not particular to your very specific problem). If you ever felt such a crave, you are in the right path.
As this process seemed long to you? Just look how short is the script with all the pieces we have used above.
TODO: add full script to make it easier to replicate
TODO: ensure your strategy can be parallelized (show list)