By: Christopher Rackauckas

Re-posted from: http://www.stochasticlifestyle.com/julia-ifem3/

This is the third part in the series for building a finite element method solver in Julia. Last time we used our mesh generation tools to assemble the stiffness matrix. The details for what will be outlined here can be found in this document. I do not want to dwell too much on the actual code details since they are quite nicely spelled out there, so instead I will focus on the porting of the code. The full code will be available soon on a github repository, and so since most of it follows a straight translation from the linked document, I’ll leave it out of the post for you to find on the github.

The Ups, Downs, and Remedies to Math in Julia

At this point I have been coding in Julia for over a week and have been loving it. I come into each new function knowing that if I just change the array dereferencing from () to [] and wrap vec() calls around vectors being used as indexes (and maybe int()), I am about 95% done with porting a function. Then I usually play with cosmetic details. There are a few little details which make code in Julia a lot prettier. For example, for the FEM solver, we need to specify a function. In MATLAB, specifying such the function for which we want to solve $-Delta u = f$ in-line would be done via anonymous functions like:

f = @(x)sin(2*pi.*x(:,1)).*cos(2*pi.*x(:,2));

I tend to have two problems with that code. For one, it’s not math, it’s programming, and so glancing at my equations and glancing at my code has a little bit of a translation step where errors can happen. Secondly, in many cases anonymous functions incur a huge performance decrease. This fact and the fact that MATLAB’s metaprogramming is restricted to string manipulation and eval destroyed a project I had tried a few years ago (general reaction-diffusion solver with a GUI, the GUI took in the reaction equations, but using anonymous functions killed the performance to where it was useless). However, in Julia functions can be defined inline and be first class, and have a nice appearance. For example, the same function in Julia is:

f(x) = sin(2π.*x[:,1]).*cos(2π.*x[:,2]);

Two major improvements here. For one, since the interpreter knows that variables cannot start with a number, it interprets 2π as the mathematical constant 2π. Secondly, yes, that’s a π! Julia uses allows unicode to be entered (In Juno you enter the latex pi and hit tab), and some of them are set to their appropriate mathematical constants. This is only cosmetic, but in the long-run I can see this as really beneficial for checking the implementation of equations.

However, not everything is rosy in Julia land. For one, many packages, including the FEM package I am working with, are in MATLAB. Luckily, the interfacing via MATLAB.jl tends to work really well. In the first post I showed how to do simple function calls. This did not work for what I needed to do since these function calls don’t know how to pass a function handle. However, digging into MATLAB.jl’s package I found out that I could do the following:

  put_variable(get_default_msession(),:node,node)
  put_variable(get_default_msession(),:elem,elem)
  put_variable(get_default_msession(),:u,u)
  eval_string(get_default_msession(),"sol = @(x)sin(2*pi.*x(:,1)).*cos(2*pi.*x(:,2))/(8*pi*pi);")
  eval_string(get_default_msession(),"Du = @(x)([cos(2*pi.*x(:,1)).*cos(2*pi.*x(:,2))./(4*pi) -sin(2*pi.*x(:,1)).*sin(2*pi.*x(:,2))./(4*pi)]);")
 
  eval_string(get_default_msession(),"h1 = getH1error(node,elem,Du,u);");
  eval_string(get_default_msession(),"l2 = getL2error(node,elem,sol,u);");
  h1 = jscalar(get_mvariable(get_default_msession(),:h1));
  l2 = jscalar(get_mvariable(get_default_msession(),:l2));

Here I send the variables node, elem, and u to MATLAB. Then I directly evaluate strings in MATLAB to make function handles. With all of the variables appropriately in MATLAB, I call the function getH1error and save its value (in MATLAB). I then use the get_mvariable to bring the result into MATLAB. That value is a value of MATLAB type, and so I use the MATLAB.jl’s conversion function jscalar to then get the scalar result h1. As you can see, using MATLAB.jl in this fashion is general enough to do any of your linking needs.

For very specialized packages, this is good. For testing the ported code for correctness, this is great. However, I hope not to do this in general. Sadly, every once in awhile I run into a missing function. In this example, I needed accumarray. It seems I am not the only MATLAB exile as once again Julia implementations are readily available. The lead Julia developer Stefan has a general answer:

function accumarray2(subs, val, fun=sum, fillval=0; sz=maximum(subs,1), issparse=false)
   counts = Dict()
   for i = 1:size(subs,1)
        counts[subs[i,:]]=[get(counts,subs[i,:],[]);val[i...]]
   end 
   A = fillval*ones(sz...) 
   for j = keys(counts)
        A[j...] = fun(counts[j])
   end
   issparse ? sparse(A) : A
end

  0.496260 seconds (2.94 M allocations: 123.156 MB, 8.01% gc time)
  0.536521 seconds (2.94 M allocations: 123.156 MB, 8.83% gc time)
  0.527007 seconds (2.94 M allocations: 123.156 MB, 9.41% gc time)
  0.544096 seconds (2.94 M allocations: 123.156 MB, 9.76% gc time)
  0.526110 seconds (2.94 M allocations: 123.156 MB, 12.22% gc time)

whereas Tim answer has less options but achieves better performance:

function accumarray(subs, val, sz=(maximum(subs),)) 
    A = zeros(eltype(val), sz...) 
    for i = 1:length(val) 
        @inbounds A[subs[i]] += val[i] 
    end 
    A 
end

Timings

  0.000355 seconds (10 allocations: 548.813 KB)
  0.000256 seconds (10 allocations: 548.813 KB)
  0.000556 seconds (10 allocations: 548.813 KB)
  0.000529 seconds (10 allocations: 548.813 KB)
  0.000536 seconds (10 allocations: 548.813 KB)
  0.000379 seconds (10 allocations: 548.813 KB)

Why is Julia “Missing” So Many Functions? And How Do We Fix It?

This is not the first MATLAB function I found myself missing. Off the top of my head I know I had to get/make versions of meshgrid, dot(var,dimension) just in the last week. To the dismay of many MATLAB exiles, many of the Julia developers are against “cluttering the base” with these types of functions. While it is easy to implement these routines yourself, many of these routines are simple and repeated by mathematical programmers around the world. By setting a standard name and implementation to the function, it helps code-reusability and interpretability.

However, the developers do make a good point that there is no reason for these functions to be in the Base. Julia’s Base is for functions that are required for general use and should be kept small in order to make it easier for the developers to focus on the core functionality and limit the resources required for a standard Julia install. This will increase the number of places where Julia could be used/adopted, and will help ensure the namespace isn’t too full (i.e. you’re not stepping on too many pre-made functions).

But mathematicians need these functions. That is why I will be starting a Julia Extended Mathematical Package. This package will be to hold the functions that are not essential language functions, but “essential math language” functions like you’d find in the base of MATLAB, R, numpy/scipy, etc., or even just really useful routines for mathematical programming. I plan on cleaning up the MATLAB implementations I have found/made ASAP and putting this package up on github for others to contribute to. My hope is to have a pretty strong package that contains the helper functions you’d expect to have in a mathematical programming language. Then just by typing using ExtendedMath, you will have access to all the special mathematical functionality you’re used to.

Conclusion

As of now I have a working FEM solver in Julia for Poisson’s equation with mixed Neumann and Dirichlet boundary conditions. This code has been tested for convergence and accuracy and is successful. However, this code has some calls to MATLAB. I hope to clean this up and after this portion is autonomous, I will open up the repository. The next stage will be to add more solvers: more equations, adaptive solvers, etc., as I go through the course. As, as mentioned before, I will be refactoring out the standard mathematical routines and putting that to a different library which I hope to get running ASAP for others to start contributing to. Stay tuned!

The post Julia iFEM3: Solving the Poisson Equation appeared first on Stochastic Lifestyle.