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 form-assembly code) and FIAT (finite element automatic tabulator). You can find the open-source 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.

fenics_ecossystem

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 MPI-GPU parallelism (DOLFIN only support MPI). This reliance on external solvers should not be overlooked: it allows the use of state-of-the-art 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 open-source 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 FEniCS-X 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 time-dependent heat equation problem. It shows how to:

  • convert the mesh in a DOLFIN’s-friendly 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, non-academic 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 dolfin-convert) 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.

dolfin_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 round-off 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 and boundary_xmax)
  • Constant specifies a constant value in all the nodes at the corresponding boundary. We could have use a float instead. Nevertheless, this way gives you already a hint on how to set up non-constant boundary conditions: use Expression 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 time-dependent 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 time-dependent 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.

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 medium-size 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 ten-fold. 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):

benchmark_coarse

For a fine mesh (~10M cells):

benchmark_refined

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 from dolfin.cpp which task is to print information (e.g. dolfin.cpp.la.list_linear_algebra_backends). This does not mean you should stop using jupyter 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 find set_log_level(DEBUG) and set_log_level(PROGRESS) as the go-to commands. Nevertheless, they do not work in version 2019.1.0, as DEBUG and PROGRESS cannot be imported. To access this variables you have to import LogLevel (from dolfin.cpp.log import LogLevel) and then access the desired variables (e.g. LogLevel.DEBUG). Alternatively , you can use an integer, being log_level=10 and log_level=16, for DEBUG and PROGRESS, respectively. Other options can be found here. I could tell you what each level prints, but why wouldn’t try it yourself? (not in a jupyter 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 script-based 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)

Like this post? Share on: TwitterFacebookEmail


Luís F. Pereira is an engineer that enjoys to develop science/engineering related software.
Antoine Dauptain is a research scientist focused on computer science and engineering topics for HPC.

Keep Reading


Published

Category

Tutorials

Tags

Stay in Touch