Automated Finite Elements: a comparison between Julia and C++

With Gridap, Julia has a Finite Element package that allows writing expressions that closely mimic the mathematical notation of the weak form of an equation and automate the assembly of the linear system from there. Rather than using macros, the equations are interpreted as regular Julia expressions. This approach is similar to what has been done in C++, e.g. in the Coolfluid 3 CFD code. In this talk, both methods will be compared, showing how Julia really is "C++ done right" for this use case.

This talk will be about comparing the implementation of the stabilized Navier-Stokes equations for incompressible flow, both using the Gridap package in Julia and using a Boost.Proto based C++ code. The concrete C++ implementation can be found here, while the equivalent Gridap code is here. Aside from the obvious differences due to the use of unicode and the fact that Gridap operates at a higher level of abstraction, there are also some striking similarities in both approaches. The main point is that both packages operate on expressions that are valid code in the programming language that is used (i.e. Julia for Gridap, C++ for Coolfluid 3). This is possible because both languages offer a lot of flexibility in terms of operator overloading and strong typing. In the case of C++, the Boost.Proto library helps with building a structured framework for the interpretation of the expressions, based on the idea of expression templates and thus avoiding runtime overhead of inheritance in C++. In Julia, this step is taken care of using the built-in type system and generated functions.

The whole objective of this type of machinery is to offer a simple interface to the user, but end up with a finite element assembly loop that is as fast as possible. To this end, information such as the size of element matrices and vector dimensions must be known to the compiler. We will show that both systems indeed achieve this, and result in good performance for the assembly loop.

Due to the similarity in approach, the experience visible to the end user is also similar: both systems exhibit long compilation times and long error messages in case of user errors such as mixing up incompatible matrix dimensions. This will be illustrated using examples.

Finally, more advanced numerical techniques, such as the stabilized methods used in fluid simulations, require the user to be able to define custom functions that are used during assembly, e.g. to calculate the value of stabilization coefficients. This is where Julia really shines, as it is possible to simply define a normal function, while in C++ some extensive boilerplate code is required, as will be shown.

The conclusion is that Gridap has reached a level of maturity that makes it very attractive to use Julia for this kind of work. Even if some performance optimization may still be needed, development of and experimenting with new numerical methods is much easier than in a complicated C++ code.