Introduction
There is a good deal of discussion in the scientific computation circles regarding the "ideal" language for developing efficient finite element codes. A cursory view of open source finite element software reveals that the dominant languages are C/C++, Fortran and Python. Each of these languages have their pros and cons: C/C++ is fast and powerful but very difficult to work with (for the typical scientific researcher who does not have an extensive background in programming), Fortran is quite outdated inspite of its amazing speed, Python is incredibly flexible and powerful for testing new ideas but does not fare well when it comes to speed. A typical compromise is to develop and test prototypes in Python and either partially or fully re-implement the computational hotspots using high performance languages like Fortran or C/C++. However mixed language programming brings its own share of difficulties. Moreover, developing efficient finite element code is crucially dependent on fast and robust linear algebra libraries and it is a non-trivial task to identify and integrate such libraries with a finite element code. From the point of view of a researcher, what is required is a language that is as easy to use as Python and as fast as C/C++ so that the major focus would be on testing new scientific ideas as against learning difficult programming paradigms and code optimization techniques. While no single language meets this requirement, a very interesting contender that comes close to the golden ratio is Julia. This article explores the utility of Julia for developing robust and efficient finite element codes.
A simple finite element code in Julia for solving the 2D Poisson equation over the unit square is presented in this article. The code uses only Julia's standard library and thus does not have any external dependencies. It solves the partial differential equation u,xx + u,yy = -6 over the unit square with the Dirichlet boundary condition u = 1 + x^2 + 2y^2 on the sides of the unit square. The code can however be easily modified to study other instances of the Poisson equation. The unit square is discretized with linear triangular finite elements and the stiffness matrix and load vectors are calculated using a one-point quadrature rule.
Background
Details of the finite element formulation for the Poisson equation and the type of mesh used can be found here. The performance of the code is compared with that of two good open source finite element software: FEniCS and FreeFem++. FEniCS is a collection of software for high level finite element code development written in Python and C++. FreeFem++ is a partial differential equations solver written in C++ and uses its own DSL (Domain Specific Language) with a C++ like syntax. There are many other excellent software which are not considered here.
About the code
Thanks to Julia's elegant syntax the code is largely self-explanatory. The finite element algorithm is presented in three files: poisson_fem.jl, finite_element.jl, mesh.jl. The implementation of a simple triangular mesh over the unit square, the mesh data types and other helper functions are provided in mesh.jl. Details of the finite element implementation for a given weak form, external load and boundary conditions are provided in finite_element.jl. The top level code that specifies the details of the problem to be solved is presented in poisson_fem.jl.
The Dirichlet boundary, the boundary condition and the external load are speficied in poisson_fem.jl using the functions on_boundary, dbc and f_ext respectively. The kernel of the "stiffness" term for the Poisson equation is obtained from the variational formulation as grad(u).grad(v). This is specified in the function poisson_stiffness, as shown below, in terms of the shape function vector N and the matrix containing the partial derivatives of the shape functions dN where dN[i,j] = dN_i/dx_j.
# Stiffness matrix
function poisson_stiffness(N, dN, i, j)
dN[i,1]*dN[j,1] + dN[i,2]*dN[j,2]
endA triangular mesh over the unit square is generated using
mesh = UnitSquare(100)The argument of UnitSquare denotes the number of uniform divisions into which each edge of the unit square is divided. The linear triangular finite element space over the mesh is initialized as
u = fe_space_C0(mesh, dbc)Finally, the solution of the finite element problem is obtained as:
solve_fe(mesh, u, poisson_stiffness, f_ext)Functions dealing with the finite element implementation are gathered in finite_element.jl. Linear shape functions over a triangle with a one-point quadrature rule are used for the formulation of the finite element matrices. The crucial part of the code that determines the efficiency of the entire analysis is the specific way in which the stiffness matrix is handled. Three strategies, referred henceforth as methods 1, 2 and 3, are explored in the function solve_fe:
# Method 1
#K = zeros(n_size,n_size)
#fe_matrices(mesh,u,stiffness,fext,K,F)
# Method 2
#K = spzeros(n_size,n_size)
#fe_matrices(mesh,u,stiffness,fext,K,F)
# Method 3
iu,ju,vu = fe_matrices(mesh,u,stiffness,fext,F)
K = sparse(iu,ju,vu,n_size,n_size)Method 1 uses a full matrix representation of the stiffness matrix. This is almost never used in practise since finite element stiffness matrices are sparse and algorithms that exploit this are significantly faster. A nice feature of Julia is its support for sparse matrices as part of its standard library. For cases when the number of non-zero elements of the sparse matrix is unknown, an empty sparse matrix can be initialized as shown in Method 2. The syntax for handling sparse matrices is identical to that of a regular dense matrix. However, inserting elements into a sparse matrix can be costly as the number of such operations increases. A further optimization can be effected by using an approximate estimate of the number of non-zero entries in a sparse matrix K to initialized three vectors I,J and V such that K[I[n],J[n]] = V[n], as illustrated in method 3. Further details regarding the implementation of method 3 can be found in the function fe_matrices.
Performance
The code execution time (in seconds), as a function of the total number of nodes (N), for the three strategies implemented in the Julia code are compared with Poisson solvers implemented in FEniCS and FreeFem++ in the table below. The FEniCS and FreeFem++ codes can be found in the references provided earlier. The values given below are averages over a few runs and would be different depending on the machines on which the codes are run.
| N | Julia | FEniCS | FreeFem++ | ||
| Method 1 | Method 2 | Method 3 | |||
| 121 | 0.78 | 0.64 | 0.99 | 0.67 | 0.01 |
| 2601 | 1.77 | 0.74 | 1.07 | 0.76 | 0.05 |
| 10201 | 44.98 | 1.43 | 1.37 | 1.00 | 0.23 |
| 40401 | - | 10.68 | 2.63 | 2.09 | 1.05 |
| 123201 | - | 87.86 | 6.29 | 5.88 | 4.03 |
| 251001 | - | - | 12.28 | 12.16 | 9.09 |
A few observations regarding the three strategies used in the Julia code are due here. The performance results highlight the reason why it is almost always a bad idea to use a full matrix representation of the stiffness matrix (method 1). Initializing the sparse stiffness matrix as in method 2 works fine for small N, but the cost of inserting new elements into the stiffness matrix is significant as the number of nodes increases. The most efficient strategy for large N is method 3 which initializes the sparse matrix with a (guessed) number of non-zero elements. The performance of method 3 in comparison with that of the FEniCS and FreeFem++ codes for much larger values of N is shown below:
It can be seen that the performance of the Julia code is as good as FEniCS' and close to that of FreeFem++ for N ranging from 10^2 to 10^6. This is quite impressive considering the ease with which the finite element solver can be implemented in Julia.
Final thoughts
The code presented here is quite elementary in comparsion to what is actually used in research. But what is really interesting about Julia is the relative ease with which various strategies can be implemented and tested without leading to code swell, while at the same time resulting in high performance code. As a preliminary assessment based on the results given here, Julia appears to be a very good choice for developing research oriented finite element software that is both fast and easy to develop. Julia does not support object-orientation, but its type system, use of multiple dispatch etc. provide a very clean and easier alternative that suits all practical needs for finite element computations. Its native support of dense and sparse matrices along with fast linear algebra solvers as part of the standard library are very attractive for finite element programming. The primary disadvantage of Julia is the absence of a large ecosystem of external libraries that other languages like Python or C++ have. This is likely to be a short term phenomenon since the language is quite new.
Trivia: the surface plot at the beginning of the article is the solution of the Poisson equation over the unit square with zero Dirichlet boundary condition on all its edges and an external load f(x,y) = 2sin(2pi*x)sin(2pi*y). Gnuplot was used to produce all the graphs in this article.
