Galerkin

From Sven Kreiss

Created on Feb 15, 2010

Being far away from home, I thought I post some pretty pictures of what I do here.

A homework for my E&M course requires to write a finite element analysis using Galerkin's method. The results are quite colorful so I decided to show them here. Also, I want to show it because I put in some effort to write C++ code that is then wrapped into a Python module through C bindings. The Fortran library LAPACK is also used from the C++ code. So it is a nice mix of languages all working together and each doing what it is doing best. The original reason to use Python was to easily plot the result, but as it turns out, it is almost trivial to also add video output. Python is really great for this. In science, a Python module should be the new executable.

Contents 1 Video: Dependence on Grid Size 2 For Physicists 3 For Programmers 3.1 LAPACK: Calling Fortran from C++ 4 More Videos: Scaling

[edit] Video: Dependence on Grid Size

Solving this hugh system of linear equations is a big task. Still, the code works reasonably fast on my laptop. It calculated all images and then stitched together the following video in 192 seconds.

UPDATE: Turns out LAPACK is much much faster on 64bit operating systems. I used a 64bit live cd of Ubuntu on the same laptop to calculate a 100x100 image in the same time as a 60x60 image on 32bit Ubuntu.

[edit] For Physicists

The physical problem is to solve for the electrostatic potential in a 2D box with uniform unit charge distribution and zero potential on all four sides.

The bottom images on the right extend this idea to charge distributions given in matrix form (see Python code below).

[edit] For Programmers

This exercise was the first time I did not use SWIG for C bindings into Python. It is actually not that hard to write them by hand. Going through C to access C++ functions is a bit of a pain especially as in the ideal case the python module would have the same interface as the public part of the C++ class. Unfortunately, I don't know of an easy way to write the Python to C++ interface directly.

Did I mention that Python is awesome? You get so much for free once the data is back in Python.

With my code, I can do something like this in Python:

import galerkin, pylab
nx = 30

#points = galerkin.run(nx) # implies uniform charge distribution
#points = galerkin.run(nx, [[1,0],[0,1]])
#points = galerkin.run(nx, [[1,0,1], [0,0,0], [1,0,1]])
points = galerkin.run(nx, [
 [1,0,0,0,1],
 [0,0,0,0,0],
 [0,0.3,0,0.3,0],
 [0,0,0,0,0],
 [1,0,0,0,1]
])
pylab.imshow(points, interpolation='nearest')
pylab.show()

This is how I created the images on the right.

[edit] LAPACK: Calling Fortran from C++

How to use LAPACK caused quite some discussions, so here is how I did it. I did not use CLAPACK. Just the plain LAPACK library which is also installed on our machines. As there is no header file, I declared the following functions in my C++ code:

extern "C" void dgetri_(int &, double *, int &, int *, double *, int &, int &);
extern "C" void dgesv_(int &, int &, double *, int &, int *, double *, int &, int &);
extern "C" void dgetrf_(int &, int &, double *, int &, int *, int &);

For plain C, those declarations should be the same without the "extern C" and all references replaced by pointers (Fortran only knows pass-by-reference, ie C pointers). I ended up using dgesv_(...) only. The function call looks like:

dgesv_(n, nrhs, matrix, n, ipiv, b, n, info);

The only other thing to do is to append -llapack when linking the program or making the Python module.

[edit] More Videos: Scaling

Scaling problem:

Scaling done right:

Update Feb 21, 2010: Processed on 64bit machine with matrices up to 10,000x10,000 entries. Charge distribution was given in a non-square matrix as:

[
 [0,0,1,0],
 [1,0,0,0]
]