Article on finite element programming in Julia

[Amuthan Arunkumar Ramabathiran]

Hello all,

I recently wrote a small article on the usefulness of Julia for developing finite element codes: http://www.codeproject.com/Articles/579983/Finite-Element-programming-in-Julia

I'm finding Julia more and more useful for my research. A big thanks to the Julia's developers!

Best regards,

Amuthan.

[Benjamin Silbaugh]

Interesting.

As a side note, I thought Julia did support OOP; i.e. multiple dispatch is a form of OOP. I believe Julia's treatment of OOP is similar to CLOS?

[Amuthan]

maybe, but the absence of member functions in Julia's type makes it different in usage compared to how objects are used in C++ or Python. I don't have a detailed knowledge the technical definitions related to OOP and so I could well be wrong here. Can anyone clarify? 

[Benjamin Silbaugh]

A few questions:

1. What algorithms does FEniCS and FreeFem++ use? That is, does either methods 1, 2, or 3 of the Julia code correspond to the exact same numerical methods used by FEniCS and FreeFem++? Might it be that FEniCS or FreeFem++ are employing more robust numerical algorithms, which may carry a higher performance penalty?

2. Why the disparity in performance between FEniCS and FreeFem++? They seem to be scaling very differently with respect to N. (This actually suggests that we're not comparing the same algorithms.)

3. Do you happen to know how Julia compares to some of the Common Lisp FEM tools such as FEMLISP? Common Lisp is another language that I think is frequently overlooked by the scientific community. (OCaml and Haskell might also be interesting to compare to...)

Also, I didn't see any verification of the Julia code. I realize that this is a blog, not a formal research publication. However, it would be helpful to see some comparison between the numerical accuracy of your Julia implementations, FEniCS, and FreeFem++.

Thanks for putting this together. I think it's a good start.

[Benjamin Silbaugh]

One minor nit-picky comment: a log-log plot of the time versus N would help to identify if each codes are using algorithms of the same complexity.

On Thursday, April 18, 2013 10:14:49 AM UTC-4, Amuthan Arunkumar Ramabathiran wrote:

[Amuthan]

Hi Benjamin,

Thanks for your comments. Yes, this article is _not_ a scientific publication (which should anyway be obvious from the freestyle presentation), but just a quick assessment of how Julia fares when it comes to handling sparse matrix computations.

I used the code given in the FEniCS tutorial (which selects the best solver?) and the FreeFem++ code uses GMRES. It is possible to assign a variety of solvers in both FEniCS and FreeFem++ and the performance results could be very different depending on the choice of the solver. The performance is also dependent on how the stiffness matrix is calculated. For instance, in this case the stiffness matrix need not be calculated for each element since its anaytical form is known. Using this could also improve performance significantly. I'm not sure how FEniCS and FreeFem++ convert the higher level weak form into efficient low level code. Anyway my intention was to see how a Julia code which exlicitly computes the stiffness matrix compares with a black box implementation of the same problem with other open source fem packages. I agree that this is not the best comparison, but just a start.

Regarding Lisp, I had a look at femlisp (http://www.femlisp.org/), but did not have the time to explore it in detail. If you have tried that earlier, could you throw some light on its performance? 

And the verification part in the Julia code is straightforward. For the case when f(x,y) = -6 and u0 = 1 + x^2 + 2y^2, the exact solution is u = 1 + x^2 + 2y^2. So one quick check would be to compare the nodal error obtained from the three codes. Just for completeness, here is what you need to add to the julia code to get this:

err = 0.0

for i = 1:mesh.n_nd

  du = u.node_val[i] - dbc(mesh.x_node[i,1], mesh.x_node[i,2])

  err += du*du

end

println(sqrt(err))

Of course, this code is quite elementary as I pointed out in the article, but based on the results from this simple comparison I think it would be worthwhile investing time to develop a full fledged finite element package in Julia. 

[Stefan Karpinski]

This is a lovely article. Thanks for writing it up! It's really exciting to see that your Julia implementation basically perfectly matches the FEniCS solver for this application. That's a pretty solid real-world validation, IMO. It will definitely be great to build a nice infrastructure for finite element problems – I know there's a lot of people here with interest in experience in this area.

[Stefan Karpinski]

On Thu, Apr 18, 2013 at 4:38 PM, Benjamin Silbaugh <ben.si...@gmail.com> wrote:

Interesting.

As a side note, I thought Julia did support OOP; i.e. multiple dispatch is a form of OOP. I believe Julia's treatment of OOP is similar to CLOS?

 

On Thu, Apr 18, 2013 at 10:45 AM, Amuthan <apa...@gmail.com> wrote:

maybe, but the absence of member functions in Julia's type makes it different in usage compared to how objects are used in C++ or Python. I don't have a detailed knowledge the technical definitions related to OOP and so I could well be wrong here. Can anyone clarify?

It all depends on what you mean by object-oriented. If you want a named bag of functions stuffed into data types like C++, Java, Ruby and Python, then no, Julia isn't object-oriented like that. Nor does it fit Alan Kay's conception of OOP as programs consisting of objects encapsulating state and interacting by sending messages to each other (which is really more like the Actor model of concurrency than what modern OO has evolved into). However, some of the concepts commonly cited as object-oriented – namely dynamic dispatch, abstraction, and polymorphism – are extremely Julian. In fact, Julia takes all three of these ideas further than most OO languages do: it does dynamic dispatch on all arguments, rather than just the first argument (or receiver); the fact that abstract types are "traits" in OO lingo, means that you can write large amounts of generic code to abstractions; and of course, Julia supports all kinds of polymorphism.

So yeah, there's no simple answer to whether Julia is object-oriented. Kind of? I actually think that dynamic multiple-dispatch bridges functional and object-oriented programming in a sense. Generic functions can be passed around as first-class objects, allowing functional style programming. However, unlike classic functional languages, these functions can be extended dynamically and made to work on new types that didn't even exist when generic code was written using them. Both CLOS and Dylan have multiple dispatch and generic functions as well, but in those languages generic functions are the exception rather than the rule, whereas in Julia generic functions are the norm and are completely everyday things upon which the entire language is based.

[Amuthan]

Thanks for your feedback! 

There are many interesting experiments one could do with Julia in the context of finite element programming considering the flexibility Julia provides. In any case, it would be quite cool to have a finite element package built entirely using Julia, or atleast develop good interfaces to widely mesh generation software like Triangle, gmsh, netgen etc. and visualization software like Paraview, to start with. I'm quite keen on taking this further ahead and if you know anyone else who is serious about experimenting with finite elements using Julia, I would be very interested in hearing from them.

[Amuthan]

Thanks for the clarification!

[Benjamin Silbaugh]

Unfortunately, I don't have any first had experience with femlisp, and my knowledge of Common Lisp (CL) is purely academic at this point. (I've been spending my coffee breaks researching language alternatives to Python/C++ for doing HPC work. So far, Julia, Common Lisp, Haskell, and OCaml have earned my attention.)

The only reason I mention femlisp in this context, is that it looks like Julia derives a number of concepts from lisp (who hasn't?), and I find myself wondering if the decision by the Julia developers to use Matlab-like syntax is helping or hurting their ability to generate optimized machine code. My understanding is that most Common Lisp implementations (e.g. SBCL) can usually match the performance of C++, and in some cases even beat the performance of C. For example there was some work done by Didier Verna back in 2003 (?) which showed that a person could match the performance of C using Common Lisp. I think the paper was titled, "Beating C in Scientific Computing". (I can look up the details if you're interested.) Of course, Lisp reads less like written English than Matlab. I suspect there many who would find the Lisp syntax disconcerting. (I only raise these questions because I'm trying to come to an objective understanding of what Julia has to offer. It is not my intention to discourage anyone from using Julia. For all I know, Julia may be the Next Big Thing.)

[Stefan Karpinski]

I don't really understand how choice of syntax would have any impact on machine code generation. High-quality Lisp implementations are fast because they're high-quality, not because they use Lisp syntax. The advantages that Julia brings to the table over Lisp for scientific computing probably does include the friendlier syntax, but there are other more substantive ones – like the expressive type system and dynamic multiple dispatch. I also think Lisp started the whole game of trying to isolate the programmer from the machine – an experiment which hasn't panned out very well for scientific computing. Julia takes a different approach and tries to make it easy to work with what the machine provides without hiding it. This includes fast fixed-size integers with wrap-around on overflow, fast IEEE 754 arithmetic, and the ability to control memory usage significantly better than other high-level languages, while not forcing the programmer to do manual memory management. It's sort of a compromise between C and Lisp – with Matlabish syntax and a dynamic type system.

[Benjamin Silbaugh]

Stefan,

Thanks for clarifying this.

I guess the reason for my confusion regarding syntax and efficiency stems from the fact that MathWorks (or anyone else for that matter) has yet to produce a Matlab JIT that yields highly performant code. Likewise, projects like PyPy and Unladen swallow seem to have had difficulties developing a Python compiler that comes close to a native C compiler. My inference was that certain code optimizations could not be performed due to ambiguities introduced by the language syntax. Perhaps I'm confusing "syntax" with some other aspect of language design.

[Stefan Karpinski]

I see. The syntax is purely the *form* of the language and has very little impact on performance. [Although certain things can make a difference – for example, the fact that in Matlab array indexing and function calls have identical syntax ends up forcing the language to determine certain things at run-time that could otherwise be decided statically when parsing code (e.g. which object an `end` used as an index corresponds to – it's the one that happens to be bound to an array object at run time, as in `f(g(h(end)))`).]

Semantics, on the other hand, is about what things *mean* and has huge impact on performance. Much of Julia's performance comes from making choices about the semantics of the language that make it relatively easy to implement in an efficient way – e.g. using machine-supported arithmetic. We're already avoiding many of the problems that languages like Python and Matlab have in that regard. Honestly, Lisp implementations are not fast because of their semantics, but despite them – Lisp is really old and a huge amount of time an effort has gone into inventing brilliant tricks that can let it be executed efficiently while maintaining its high-level behavior.

JMW's blog post from a while back and some of the back-and-forth in the comments address some of this.

[Benjamin Silbaugh]

Interesting. Thanks!

[Amuthan]

Hi Benjamin,

Sorry to interrupt the very interesting conversation between you and Stefan, but I have a quick answer to the question you raised earlier regarding the huge disparity between the FEniCS and FreeFem++ codes. Turns out that it was a bad idea to specify the solver explicitly as GMRES in the FreeFem++ code. I removed this and let FreeFem++ choose the best solver and now the FreeFem++ code is the fastest of the lot for mesh size ranging from 10^2 to 10^6 nodes! I have updated the results (http://www.codeproject.com/Articles/579983/Finite-Element-programming-in-Julia). The C++ code is the fastest now, as one would expect, but the conclusions about Julia's performance and its usefulness for developing large finite element codes remain true nevertheless.

Thanks for pointing out the discrepancy in the results between FEniCS and FreeFem++.   

[Benjamin Silbaugh]

Glad I could help.

[Viral Shah]

Hi Amuthan,

This is a really interesting read. A couple of orthogonal comments:

The sparse matrix implementation in julia is fairly young, but it is nice to see that it is already usable. Much of the sparse matrix implementation (except the solvers) is written in julia, and I expect it to get faster as the compiler improves.

In julia, you are using \ to solve the system. This currently defaults to using LU factorization. However, if your system is SPD, there are gains from using a cholesky solver. Matlab does this automatically, and Julia does this automatically in the dense case, but not yet in the sparse case. I suspect that your system is not SPD from your choice of GMRES. 

We currently do not have the preconditioned iterative solvers. These should be easy to add, and I suspect this will happen as soon as there is demand.

By default, on 64-bit systems, you get sparse matrices with 64-bit integers. I suspect that there could be a modest performance gain (10-15%) by using sparse matrices with 32-bit integers. Some of this needs to be documented. 

-viral

[Amuthan]

Hi Viral,

Thanks for your comments. 

 

The sparse matrix implementation in julia is fairly young, but it is nice to see that it is already usable. Much of the sparse matrix implementation (except the solvers) is written in julia, and I expect it to get faster as the compiler improves.

This is fantastic. Looking forward to it!

 

In julia, you are using \ to solve the system. This currently defaults to using LU factorization. However, if your system is SPD, there are gains from using a cholesky solver. Matlab does this automatically, and Julia does this automatically in the dense case, but not yet in the sparse case. I suspect that your system is not SPD from your choice of GMRES. 

Oh ok, I thought Julia chooses the best solver in the sparse case too. I'll give a quick try with a cholesky solver along the lines of how LU decomposition is implemented in Julia and see how it improves the performance.

And the matrix coming out of this finite element implementation is indeed SPD. I had taken a default code from the FreeFem++ example list and failed to notice that the solver type was set to GMRES initially. My bad!

 

We currently do not have the preconditioned iterative solvers. These should be easy to add, and I suspect this will happen as soon as there is demand.

This would be quite a welcome addition. Finite element matrices are almost always sparse and having good preconditioned iterative solvers becomes crucial for high performance codes. Is anyone currently working on this? I would be very happy to contribute whatever little I can.

 

By default, on 64-bit systems, you get sparse matrices with 64-bit integers. I suspect that there could be a modest performance gain (10-15%) by using sparse matrices with 32-bit integers. Some of this needs to be documented. 

Interesting. I'll try that out too.

Thanks once again for your detailed comments!

Amuthan.

[Viral Shah]

Great. I'll be happy to help in any way I can. Here's a Cholesky example - which also allows reusing of the cholesky factor to solve multiple RHS.

julia> a = speye(100);
julia> a = speye(100); b = rand(100);
julia> r = cholfact(a);
julia> x = r \ b

CHOLMOD dense: : 100-by-1,
leading dimension 100, nzmax 100, real, double
col 0:
0: 0.73927
1: 0.38795
2: 0.21423
3: 0.058064
...
96: 0.81107
97: 0.6066
98: 0.35818
99: 0.96451
OK

julia> norm (a*x.mat-b) # Note: This will be fixed so that r \ b gives a regular dense vector.
0.0

-viral

[Amuthan]

Thanks for the tip Viral. Using cholfact results in a significant performance improvement for large N - for a mesh with 1 million nodes the computational time came down from about 64s to 40s. This is really awesome!

It would be great to see some fast iterative solvers (for both dense and sparse matrices) and hopefully a good interface to high performance packages like PetSc (http://www.mcs.anl.gov/petsc/) 

[Viral Shah]

There is a discussion thread on julia + PETSc.
https://github.com/JuliaLang/julia/issues/2645

There is also some preliminary work:
https://github.com/qsnake/petsc/blob/master/bin/julia/PETSc.jl

-viral

[Viral Shah]

I am curious if using SparseMatrixCSC{Float64,Int32} will lead to further speedups. This could lower memory usage, and indexing overhead inside CHOLMOD. On the other hand, the supernodes inside CHOLMOD may already be removing much of this indexing overhead, and it may not matter.

-viral



On 22-Apr-2013, at 2:21 PM, Amuthan <apa...@gmail.com> wrote:

[Amuthan]

Thanks for the links to the earlier discussions. 

Using SparseMatrixCSC{Float64,Int32} seems to slow down the computations (to roughly the same same speed as that obtained using an LU decomposition). I'm not sure if I'm implementing it in the best possible way, but the essential changes I made (in the fe_matrices function in finite_element.jl) were:

try 

   iu[sparse_count] = convert(Int32,ig)

   ju[sparse_count] = convert(Int32,jg)

   vu[sparse_count] = ke[i,j]

catch BoundsError()

   iu = vcat(iu,convert(Int32,ig))

   ju = vcat(ju,convert(Int32,jg))

   vu = vcat(vu,ke[i,j])

end

This doesn't seem to help. Changing all the ints to Int32 is difficult to work with since every integer operation needs to be converted to Int32. Or maybe there is a simpler way to do that? 

[Viral Shah]

This is a good enough way to do it - and it answers the question. 

Thanks,

-viral

[Petr Krysl]

Hello everybody,

I found Amuthan 's blog a while back, but only about two weeks ago I found the time to look seriously at Julia. What I found was very encouraging.

For a variety of teaching and research purposes I maintain a Matlab FEA toolkit called FinEALE. It is about 80,000 lines of code with all the examples and tutorials. In the past week I rewrote the bits and pieces that allow me to run a comparison with Amuthan 's code. Here are the results:

For 1000 x 1000 grid (2 million triangles):

Amuthan's code: 29 seconds

J FinEALE: 86 seconds

FinEALE: 810 seconds

Mind you, we are not doing the same thing in these codes. FinEALE and J FinEALE implement code to solve the heat conduction problem with arbitrarily anisotropic materials. The calculation of the FE space is also not vectorized as in Amuthan's code. The code is written to be legible and general: the same code that calculates the matrices and vectors for a triangle mesh would also work for quadrilaterals, linear and quadratic, both in the pure 2-D and the axially symmetric set up, and tetrahedral and hexahedral elements in 3-D. There is obviously a price to pay for all this generality.

Concerning Amuthan 's comparison with the two compiled FEA codes: it really depends how the problem is set up for those codes. I believe that Fenics has a form compiler which can spit out an optimized code that in this case would be entirely equivalent to the simplified calculation (isotropic material with conductivity equal to 1.0), and linear triangles. I'm not sure about freefem++, but since it has a domain-specific language, it can also presumably optimize its operations. So in my opinion it is rather impressive that Amuthan 's code in Julia can do so well.

Petr

 

[Rob J. Goedman]

Petr,

Are you referring to http://www.codeproject.com/Articles/579983/Finite-Element-programming-in-Julia ?

Or is this from another blog?

Rob J. Goedman
goe...@mac.com

[Petr Krysl]

Correct. That is the blog cited by the OP in this thread.

[Petr Krysl]

Sorry: minor correction.  I mistakenly timed also output to a VTK file for paraview postprocessing.  This being an ASCII file, it takes a few seconds.

With J FinEALE the timing is now 78 seconds.

Petr

[Amuthan]

Hi Petr: thanks for the note. It is definitely true that making the code more general to accommodate a larger class of elements and/or PDEs would add some computational overhead. I wrote that code primarily to assess the ease with which one could implement a simple FE solver in Julia; the code is admittedly restrictive in its scope and not very extensible either since a lot of the details are hard-coded, but that was a deliberate choice. I had some plans of developing a full fledged FE solver in Julia last year but dropped the idea since I shifted to atomistic simulation and no longer work with finite elements. I would still be interested in hearing about how a general FE solver in Julia performs in comparison to an equivalent FE software. If you have some pointers on that please drop in a note. Thanks!
Best regards,
Amuthan

[Mauro]

There was one talk about FEM at JuliaCon:
https://www.youtube.com/watch?v=8wx6apk3xQs

On Mon, 2014-12-08 at 04:48, Petr Krysl <krysl...@gmail.com> wrote:
> Bit more optimization:

>
> Amuthan's code: 29 seconds

> J FinEALE: 54 seconds
> Matlab FinEALE: 810 seconds

[Petr Krysl]

Thanks!