[SciPy-dev] the scipy mission, include finite element solver

[Ondrej Certik]

Hi,

the mission of scipy from scipy.org:

"
SciPy is open-source software for mathematics, science, and
engineering. [...] The SciPy library [...] provides many user-friendly
and efficient numerical routines such as routines for numerical
integration and optimization.
"

it has solvers for ODEs using several methods. Is having a finite
element (FEM) solvers for PDEs (and some ODEs too) among the scipy's
goals?

My group[0] would be very interested in that.

Of course, this is not something, that can happen overnight. But scipy
has interface to lot's of sparse solvers, so it might have some nice
pythonic interface to finite element solvers for PDEs too.

If you look at the topical software, PDE section:

http://scipy.org/Topical_Software#head-3df99e31c89f2e8ff60a2622805f6a304c50101f

there are currently 3 solvers listed, FiPy (finite volumes), SfePy
(regular finite elements, Robert Cimrman, some patches from me and
other people), Hermes (higher order adaptive finite elements,
developed at the University of Nevada/Reno[0], where I am now). There
are more good opensource FEM libraries out there, some list is here:

http://hpfem.org/femcodes.html

the most widely used is probably libmesh, though it's GPL and as far
as I know it doesn't yet have python wrappers (I wrote some swig
wrappers couple years back, but it was bad --- it's on my todo list to
write good wrappers in Cython for libmesh).

All FEM codes need a common functionality --- you need to load a mesh
(e.g. readers and writers for all the different mesh formats out
there, commercial or not), you need to define your PDE, then you need
to assemble the global matrices (this is the, where all the FEM codes
differ, e.g. you can have lots of different bases, etc. etc.), then
depending on your problem, in most cases you end up either with linear
system of equations:

A * x = b

or (eigenproblem):

A * x = lambda * B * x

where A, B are (sparse) matrices, "x" is the solution, "b" is the
right hand side (vector), lambda is the eigenvalue. So the goal of the
FE solver is to calculate A, B and "b". It should be general enough to
allow you to discretize any integral in the weak formulation of the
PDE, but in most cases, A, B or "b" is all you want.

Then I solve it using solvers that are in scipy, or pysparse, or
petsc4py, or slepc4py.

Next one has to take "x" and feed it to the finite element solver
again and then tell it to automatically adapt the mesh (there are
several approaches for this, you can use some not so good error
estimators, like some gradients of the solution, or you can use
sophisticated errror estimators, like a reference solution, calculated
on a uniformly refined mesh etc.). So you iteratively adapt, until you
are satisfied.

Most problems are nonlinear, so one also has to use some nonlinear
solvers, e.g. I wrote some Jacobian free nonlinear solvers in
scipy.optimize.nonlin, but this has to be polished and improved (and
also if your problem is simple enough to allow it, you can also
construct the Jacobian using FEM directly, that's the best thing).

Once you have the solution, you also need to save it to many different
file formats for visualization, and/or provide hooks to mayavi or
matplotlib (for 2D stuff).

And besides all of this, one also needs to have a common functionality
like doing arithmetics with the solution (e.g. sum/substract two
solutions, integrate it, etc.) , project it to different meshes, etc.
etc.



So do you think scipy could have some default, adaptive FEM solver,
that would do the job for any PDE, and then all the common
functionality above, so that one can easily hook in there any other
solver, just to construct the A, B, "b" matrices/vectors and then use
scipy for the rest (and optionally the solver again for the
adaptivity)?


This is a question where you want scipy to go in the longer term. But
I think a good mission statement is to provide a viable alternative to
Matlab. And for that, a good PDE solver and related tools are
necessary. And if scipy had all the common functionality above, so
that one can easily hook up any FEM solver in there, plus there would
be some default solver in scipy itself, imho that would rock.

In my group[0] (we also collaborate with Robert Cimrman) we are doing
all of the above and our goal is to have such a system anyway. But it
occurred to us that if this could be part of scipy, it could be a win
to everybody, since I could spend more time working on scipy and all
scipy users would have a good FEM solver. Plus I would then have
imminent interest to release soon, release often, which imho is not
bad for scipy either.

As to techinical side, the 2D solver is written in C++ and builds in
about 3s on 8 processors, and the python wrappers use Cython. The 3D
solver takes a bit longer to build I think 17s. So it would make the
scipy build a bit longer, but imho scipy should be able to solve PDEs.

Together with a project like SPD:

http://code.google.com/p/spdproject/

which easily builds scipy, numpy, atlas, blas, ..... on all platforms
from source easily, then imho we could have a viable alternative to
matlab.

Ondrej

P.S. Currently our FE solvers is GPL, but I think that could change to
BSD if there is enough interest in that (there is in me). SPD also has
to be GPL, because it's based on Sage, but if there is enough
interest, I guess it's not difficult to rewrite all the buildscripts
from scratch and use BSD (I myself don't mind using GPL for SPD and
save myself lots of time by taking advantage of Sage).

[0] http://hpfem.math.unr.edu/