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

Mon Apr 6 05:16:36 EDT 2009

As a user, I would very much like to see scipy include some FEM base.
As you say, a generic framework would allow others to just write a different
core, support other mesh formats, ..., eg they could be done as scikits.
The fragmentation now is large.

I am using freefem++ (http://www.freefem.org/ ) recently for quick reference
solutions in 2D. Would you plan such a high level parsing of weak forms?

Anyway, the question is always to use C/C++ directly or python to do the
job. Personally I am leaning to using DUNE http://www.dune-project.org for
my next project as it is very versatile. I know they where playing with a
python interface too, but I find nothing on their site now about this.

Benny

2009/4/5 Ondrej Certik <ondrej at certik.cz>

> 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/
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20090406/2208a781/attachment.html>