Mid term accomplishments

I will try to give a comprehensive feel of what I achieved in this first part of the Google Summer of Code, since it is time for the mid term evaluation. Let's start with an example: as usual, it is the Poisson equation, but today, as a twist, we consider a fully Neumann problem. In order for such a problem to be well posed there is the need of an additional constraint, otherwise the solution would not be unique, so in the Octave code there is the Lagrange multiplier c. Here you can find more details and the C++ and Python code, I will just write down the differential problem for convenience:

- Δu = f in Ω

∇u ⋅ n = g on ∂Ω

Here is the Octave code that solves the above mentioned problem:

pkg load fem-fenics msh

ufl start NeumannPoisson
ufl CG = FiniteElement '("CG", triangle, 1)'
ufl R = FiniteElement '("R", triangle, 0)'
ufl W = CG * R
ufl
ufl "(u, c)" = TrialFunctions (W)
ufl "(v, d)" = TestFunctions (W)
ufl f = Coefficient (CG)
ufl g = Coefficient (CG)
ufl
ufl a = "(inner (grad (u), grad (v)) + c*v + u*d)*dx"
ufl L = f*v*dx + g*v*ds
ufl end

# Create mesh and function space
x = y = linspace (0, 1, 33);
mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4));

W = FunctionSpace ("NeumannPoisson", mesh);

# Define variational problem
f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02));
g = Expression ('g', @(x,y) - sin (5.0 * x));

a = BilinearForm ("NeumannPoisson", W, W);
L = LinearForm ("NeumannPoisson", W, f, g);

# Compute solution
[A, b] = assemble_system (a, L);
sol = A \ b;
solution = Function ('solution', W, sol);

u = Function ('u', solution, 1);

# Plot solution
[X, Y] = meshgrid (x, y);
U = u (X, Y);
surf (X, Y, U);

At the very beginning you can see a block with every line starting with ufl. That is what you would have to put in a separate UFL file before this summer. In a sense it is not plain UFL, but there are extra quotes and apices. They are needed because, using the current version of Octave, those brackets with commas inside would otherwise be interpreted as function calls. After this blocks closes with the ufl end line, the resulting UFL file is compiled to obtain a FunctionSpace, a BilinearForm and a LinearForm. These are oct-files that fem-fenics will use later on to define the corresponding variables in Octave. A robust implementation of ufl.m, the function that provides this binding to the UFL language, is one of the results of the first term.

In the end of the snippet you can see that the solution u is evaluated in its domain exactly as you expect to do with a regular function taking two arguments and returning one value. This is due to the new subsref method of the function class, which is used to represent the elements of a function space. Aside from surface plots, this feature can be of interest to generalise methods that rely on analytical solutions to differential problems, or to apply basically any algorithm to such functions. Here is the plot you will obtain with this script:

I wrote in an earlier post of the interpolate function: with this you can get the representation of a Function or Expression on a given FunctionSpace. It is useful, for instance, to compare your numerical solution with an exact one you happen to know. Or, in the example above, you might want to view what is the forcing term like:

f_cg = interpolate ("f_cg", f, u);
F = f_cg (X, Y);
surf (X, Y, F);

There is one last achievement to highlight for the mid term evaluation: currently both the initial compilation of the package and all the ones performed just-in-time when importing UFL instructions proceed smoothly without user intervention. To this end, now the build system relies on pkg-config to get at once all the flags needed for proper compilation and linking, since some dependencies of dolfin, the FEniCS interface, are not to be found in standard directories. In order to exploit the extracted information also for the subsequent run time builds, the autoconf substitution is performed also in the get_vars.m auxiliary function, which in turn provides it to generate_makefile.m. An implementation detail that proved quite tricky is how to pass all the preprocessor flags to mkoctfile: only a subset of the options of g++ are hard-coded in it, so I needed to resort to a workaround. Indeed, CPPFLAGS are always passed as environment variables and not as command line flags, so that mkoctfile will just copy and deliver them to the real compiler.

To further enhance the build system, I implemented other internal functions that hash the UFL file that was compiled and, later, check it to understand if it changed between the previous and the freshly requested build. In the example above, you will find in your working directory four new files after a run: the three already mentioned oct-files and a text file storing the md5 sum of the UFL that has been imported. Until one of these files gets somehow deleted or the problem in the ufl block changes, you will not need to take on a time consuming compilation any more.

Function evaluation

As said in the previous post, my latest result is the possibility to evaluate a fem-fenics function on a point of its domain. This way it is possible to generalise methods that, otherwise, would rely on analytical solutions. In [1] we have an example of such a method.

The paper deals with the deposition of nanoparticles in tissues, for the treatment of cancer. The phenomenon is described with a Monte Carlo simulation of these particles' trajectories, assuming that the velocity field of the carrying fluid is known. In this study, some simplifying hypotheses about the geometry of the cells and the fluid layer nearby allow for an analytical solution of the Stokes equation. Unfortunately, these assumptions do not hold generally in human tissues: for instance, in the liver cells have cubic shape, contrasting to the spherical one used in this paper. Now, if this method is implemented in Octave, we can solve numerically the Stokes equation on a realistic domain and obtain right away a more general approach to this significant application.

The evaluation of a fem-fenics function was already possible via the feval method, but it had some glitches. One aspect is that the solution of a differential problem could not be used as if it was a regular Octave function, then a user should have adapted his/her algorithms to take advantage of it. One more critical issue is that the previous implementation did not handle the exception raised by the underlying FEniCS method when it gets as argument the coordinates of a point outside of the domain, thus leading to a crash of Octave.

In order to address these problems, I added the subsref method to the function class and implemented the proper exception handling in feval. To avoid code duplication, the former relies on the latter for the real computation, so it basically just forwards the parameters after checking that the right type of indexing was used. As a result, it is now possible to solve the equations:

- ν Δu + ∇p = 0

∇ ⋅ u = 0

with relevant border conditions, on a proper mesh and finite element space, and then evaluate the solution with the Octave expression  values = u (points), where points is a matrix holding the coordinates of every point where to do so, one per column. Moreover, a careless evaluation will not result in your Octave session crashing any more.

Even if this feature of the package underwent some improvement, there is still room for more. Two issues I have not addressed yet are the somehow weird interface and the possibility to create a function handle to perform evaluations with. Regarding the former, we might observe that the above mentioned expression remains exactly the same no matter what the geometrical dimension of the domain is. I should modify the implementation so that a vectorial function on a 3D space is evaluated with [ux, uy, uz] = velocity (x, y, z). Moving to the latter, in my understanding the class design should be modified to allow the exploitation of the Octave internals managing functions, so this would require a careful reflection on all the possible collateral effects of such a change.

[1] http://link.springer.com/article/10.1007%2Fs11517-010-0615-0

Goals for future development

The mid-term review is approaching, so it is time to highlight what is done, what is underway and what are the future goals. In this post I will try to do so as clearly as possible.

Mid-term review

The main effort during the first part of the Google Summer of Code was the implementation of the bindings for the UFL language in Octave. Now UFL code can be written directly in m-files, without the need of a separate file to define the problem. To this end, ufl has been implemented for opening a file, writing to it and importing the variational problem when it is complete.

Further, I implemented interpolate, which allows the interpolation of a Function or an Expression on a given FunctionSpace. This can be of interest to test the validity of a discretisation method, for instance if an analytical solution is available in closed form, so that it is possible to compare it with the numerically obtained one.

Lately, I focused on the build system, both for the package compilation and for the just-in-time ones needed to import variational problems in Octave. The former is now backed by pkg-config, so that all the proper compiling and linking options required by the dependencies are obtained at once. Thinking about the latter, this information is used to accordingly configure the get_vars function, which provides it to the one that generates the Makefiles used to compile oct-files just-in-time. In the end, currently these oct-files are compiled again only when necessity arises, for example if one of them has been deleted or if the UFL file has been changed.

In the upcoming week I will add another feature: it will be possible to get a function handle for the evaluation of a Function. This way the solution of a variational problem can be used exactly as any other function in Octave, for instance allowing the generalisation of algorithms relying on exact solutions of differential problems, which are thus limited to simple cases. I will provide some details on an application in my post about this feature.

Final review

In the second part of the project I will be mainly committed to the parallelisation of the package execution via MPI. As noted in an earlier post, the parallelisation through the OpenMP paradigm has been quickly abandoned because it does not provide a significant performance gain, while opening the way to bugs and errors. Parallelism is, anyway, an interesting feature for the package's use cases, so it will be the main goal of the final hand in.

Just-in-time compilation

One of the known issues of the fem-fenics package was related to the errors during the just-in-time compilation, due to missing include directories. Among my preliminary contributions there is a changeset addressing the problem in the initial build of the package, when installing it into Octave. Now I went on and solved it also during the usage of fem-fenics.

At the moment of the first build, autoconf is in charge of finding out the relevant compiler and linker flags through pkg-config. They are, then, substituted in the Makefile, which compiles the package making use of them. This piece of information is needed also when an UFL file is imported and transformed into the oct-files used to transform the weak formulation at hand into an algebraic system, but until now the user had to supply it by the means of an environment variable.

Currently, I added a new utility function that provides those flags. In the configuration process they are substituted also in get_vars.m, which is called by generate_makefile.m when a differential problem is imported. The latter replaces two placeholders and writes the ad hoc Makefile with all the necessary compile and link options. This way users will not need to provide compilation flags anymore, instead the package will manage this aspect on its own.

As noted in a previous post, however, this just-in-time build is relatively time consuming, taking around half a minute each time. Nonetheless, a common usage pattern could entail the resolution of the same weak formulation on varying meshes or with different boundary conditions, forcing terms or physical parameters. Every mentioned situation does not need the recompilation of the problem's oct-files, since they carry information only about the function space and the formal expressions of the bilinear form and the linear operator. It is useful, then, to take on the build process only when needed.

To add this feature, I created three function to perform appropriate actions. After every successful just-in-time compilation, save_hash.m takes care of hashing the imported UFL file and writing the result to <ufl-filename>.md5sum. On the other hand, at the beginning of every import_ufl_*.m function, a check like this is performed:

         if (check_hash (var_prob) ||

             ! check_oct_files (var_prob, "Problem"))

You can see in it the remaining two functions implemented lately. The first one, check_hash.m, receives as argument the name of the variational problem, reconstructs the UFL file name, looks for a saved hash sum and compares it with the current file's. It returns true if the proper .md5sum file is not found or if the new and old hashes are not the same. Clearly, the oct-files should be rebuilt if one of them is missing: check_oct_files.m looks for the relevant files, with its second option stating which import is underway (thus, which files are expected as output), and returns true if they are all available.

OpenMP

These days I have worked on the implementation of the interface to the OpenMP-powered assembly offered by FEniCS. Despite being potentially a one-line intervention, it proved quite tricky: indeed, with the needed addition, the fem-fenics function for system assembly broke with a huge number of run time errors, probably due to a change in the underlying data structure that is transparent to the library users, but does not go unnoticed if you need to access it directly, as fem-fenics does. This led me to leave this functionality behind.

My choice is backed by some computational experiments. They show that the approach enacted by the FEniCS library is quite effective, with times required for assembly reduced by half using four threads instead of just one. However, they are negligible compared to the linear system solve phase, even when the OpenMP parallelisation is disabled. I used a great number of mesh nodes in order to have meaningful timings: even if linear systems took some minutes for resolution, the assembly phase lasted as much as a couple of hundredth of a second in serial code. If we add to these findings that the fem-fenics package requires a just-in-time compilation lasting around half a minute, we understand that there is no point in devoting effort for the implementation of this feature.

Towards a parallel fem-fenics

These days I have been studying the documentation of the FEniCS project, mainly the FEniCS book, in order to understand the features related to parallel execution that it boasts. This preliminary study is aimed at adding them to the fem-fenics package. First of all I will summarise my findings, then I will comment the problems I need to address to implement this functionality.

Parallelism in FEniCS

FEniCS implements parallelism in such a way to be transparent to the user of the library. Moreover, it scales on different architectures, ranging from multi-core personal computers to distributed clusters. To this end, FEniCS makes use of two paradigms, which can be exploited both separately and together.

The first approach is tailored for shared memory architectures, such as the vast majority of the PCs nowadays, but also in many cases each node of a computational cluster. The implementation is based on OpenMP, and adding a simple instruction one can enable parallelisation to speed up the matrix assembly phase. It should be noted that this paradigm has little support in the underlying linear algebra libraries, so the resolution phase can take advantage of multi-threading only with the PaStiX solver. Since in a shared memory model parallel programs might suffer race conditions, the mesh is coloured to identify subsets, so that no two neighbouring elements belong to the same set. Obviously, the notion of proximity depends on the particular function space, then this is considered in the colouring algorithm. The assembly proceeds iterating over colours and splitting their nodes among threads: with this technique race conditions are avoided and the user can enjoy the benefits of parallelisation without incurring in unpredictable behaviour.

Contrasting to the first approach, the second paradigm is based on MPI and addresses the needs of distributed memory architectures. Unfortunately, the latter is less immediate than the former, requiring a DOLFIN program to be launched with the MPI execution utility, but in this case the code need not be modified. In this implementation, the mesh is split so that each process gets its part of it, with an algorithm striving to minimise inter-process communication. With scalability in mind, no single process holds the full matrix and, moreover, everything happens behind the scenes: this way the user has no need of taking care of low level issues. The distributed memory paradigm is diffusely supported in the algebraic back-ends, so it allows the usage of several solvers, both indirect and direct. As already noted, this and the previous approach can be combined, for instance distributing the computation on a cluster and further speeding up the assembly process enabling multi-threading within each node, provided they are multi-core machines. 

The implementation in fem-fenics

The shared memory paradigm should be quite straightforward to implement in fem-fenics. I expect to operate on a couple of functions: the private generate_makefile and the two assemble and assemble_system. The former should have the proper compilation flag (-fopenmp) added. The latter should have a new line reading like:

dolfin::parameters["num_threads"] = femfenicsthreads;

The number of threads could be passed to those functions as an argument, but this would ruin the interface compatibility with FEniCS, so this is a poor approach. Another way of addressing the issue is to define a global Octave variable in PKG_ADD and store in it the desired number of concurrent threads to use for the assembly.

The implementation of the distributed memory paradigm, instead, seems quite tricky. Basically, Octave does not use MPI, at least not Octave core. Nonetheless, there are two Forge packages with this goal, mpi and parallel. I will go through the documentation of these packages to understand if and, in case, how they address the problem of launching the oct-file with mpirun or mpiexec. Even leaving this aspect aside, I still do not know how easily the distributed objects storing matrices and vectors can be accessed to obtain the whole data.

In conclusion, I will initially work to add shared memory parallelism, at the same time looking deeper into the issues related to the distributed memory paradigm, which I suspect of being more than the ones highlighted.