[SciPy-Dev] GSoC' 15 Ideas: define general integrators; parallelize integration; higher-order ODEs

Tue Mar 17 05:26:45 EDT 2015

2015-03-17 0:42 GMT+01:00 Geordie McBain <gdmcbain at freeshell.org>:

> 2015-03-17 3:37 GMT+11:00 Max Mertens <max.mail at dameweb.de>:
> > Hello everyone,
> >
> > I am a bit late to write to this mailing list, therefore sorry for that.
> >
> > My name is Max Mertens; I am a Electrical Engineering student at the
> > University of Ulm, Germany.
> > I have done a few projects in C(++) on physical simulations of multibody
> > dynamics with collision detection, and molecular dynamics with
> > interatomic potentials.
> >
> > I got interested in SciPy as it facilitates physical simulations in an
> > easy way, but have not worked with SciPy so far. I would like to
> > contribute to SciPy as part of the GSoC, and am willing to contribute by
> > providing an "easy-fix" to one of the open issues.
> > Is there already a possibility to perform multibody/molecular dynamics
> > in SciPy or are there contributions needed in this scope?
> >
> > The scipy.integrate module with ode and odeint provides several methods
> > to integrate differential equations. For this I have the following ideas:
> >
> >   * add a interface/class/method to define general integration methods:
> >     For now you can specify various integrators like "vode" and "dopri5"
> >     as a string. The new code would allow to enter a Butcher tableau to
> >     define implicit/explicit (embedded) Runge-Kutta methods, which would
> >     cover "dopri5" and "dopri853" (see [0]); and possibly other general
> >     integrator descriptions as needed.
>
> Excellent idea.  I would like this very much.  I had to write
> something like this with Butcher tableaux for general implicit
> Runge-Kutta methods (for differential-algebraic equations rather than
> ordinary differential equations), but it's a bit in-house and
> application-specific for publishing.  Having this as part of SciPy
> would be great.
>
> >   * add distributed integration: Linear multistep integrators like
> >     Runge-Kutta with multiple differential equations can be parallelized
> >     to speed up calculations. This would distribute integrations to
> >     multiple threads; and/or, if needed for complicated equations like
> >     in multibody dynamics, distribute to multiple physical machines (I
> >     have developed a few ideas on how to accomplish this).
>
> Another excellent idea and something I've had on my own TODO list for
> a while.  I'm not sure what the SciPy policy on parallelization is
> though; SciPy doesn't seem to deal with it much.
>    I would be very interested to read your ideas on this.  Ultimately
> I suppose I would like my integrator to use a sparse parallel linear
> solver.
>

This for stiff problems and implicit methods! The python code would not
care as long as normal arrays can be passed. Most specific implementations
have their own datastructures however.

Code like fipy casts PDE problems to equations which can be solved via
pytrilinos, http://trilinos.org/packages/pytrilinos/ in parallel via MPI.
You get out numpy arrays, internally something else is used.

The cvode solver of sundials for stiff problems has a fully parallel
implementation, and sundials has some RK methods in ARKode. I don't think
any of the python bindings expose that at the moment. I suppose the mapping
to parallel array is not straightforward. Having a look at the parallel
implementation there might also give ideas though for a general interface.
Pytrilinos interfaces parts of sundials via Rythmos (
http://trilinos.org/docs/dev/packages/rythmos/doc/html/classRythmos_1_1ImplicitBDFStepper.html),
but it's unclear to me if that is parallel or not. In any case, the data
structures used are not ndarray, but must be mapped to Epetra.Vector. All
far too high level abstraction to be usable in scipy. Scipy should focus on
simple interfaces for ode/dae, but the existing examples of parallel
implementation seem to indicate simple is no longer possible then.

Whatever the approach, scipy should not redo work present in high level
packages like pytrilinos, but instead offer the basis to start from, so
people can evolve to those packages if needed.

Benny

> >   * provide methods to integrate higher-order differential equations: is
> >     this needed, or are users of the library expected to express these
> >     as multiple first-order differential equations? Could this step be
> >     automated?
>
> I think this could be worthwhile too.  At the very least for
> second-order equations which arise so often in physical applications
> and for which there are popular special methods like
> https://en.wikipedia.org/wiki/Newmark-beta_method and
> https://en.wikipedia.org/wiki/Verlet_integration.
>
> Another reason I had to write my own implicit Runge-Kutta methods was
> because the differential equations that I was dealing with weren't
> 'ordinary' in the sense that they could be solved for the first
> derivative so were 'implicit' or 'differential-algebraic'.  This
> wasn't a particularly difficult generalization but is not currently
> covered by anything in scipy.optimize; e.g. scipy.optimize.ode insists
> on a form y' = f (t, y) but my equation was more like f (t, y, y') =
> 0, as in
>
> Hairer, E., C. Lubich, & M. Roche (1989). The numerical solution of
> differential-algebraic systems by Runge-Kutta methods, Volume 1409 of
> Lecture Notes in Mathematics. Berlin: Springer
>
> Brenan, K. E., S. L. Campbell, & L. R. Petzold (1996). Numerical
> solution of initial-value problems in differential-algebraic
> equations, Volume 14 of Classics in Applied Mathematics. Philadelphia:
> Society for Industrial and Applied Mathematics
>
> DAEs can be much trickier than ODEs though (as described in those two
> books and Petzold's other papers), so it is harder to write robust
> general-purpose programs for them for inclusion in something as
> high-level as SciPy; e.g. this is a large part of why I describe my
> solutions as too application-specific.
>
> Good luck.  There is much worthwhile work to be done here.
> _______________________________________________
> 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/20150317/4769b5fb/attachment.html>