[SciPy-User] Speeding up odeint for big problems

[Justin Bois]

Hello all,

I'm using odeint to solve a system of PDEs using the method of lines.  This works well for problems in one dimension where there are about ~10^2 coupled ODEs.  In 2D, I have ~10^4, and speed is a major issue.  Does anyone have any general tips on ways to optimally use odeint to solve large problems using method of lines?  Some useful info/questions about my problem:

--The Jacobian is dense, even if I use finite differencing for the spatial derivatives.
--Using something like PyDSTool might be difficult to implement because I need many of SciPy's handy functions like FFTs, etc.
--I'm worried I might be wasting time by passing a large object as an argument to the RHS function.  This is out of convenience because it contains useful information about the problem.  Does odeint do type conversions on all attributes/methods of the input argument, or only on what is needed?  Would paring down the input save time?

Any input is greatly appreciated.  I'd really like to avoid having to code this up in C.

Best,
Justin

[Rob Clewley]

> I'm using odeint to solve a system of PDEs using the method of lines. This
> works well for problems in one dimension where there are about ~10^2 coupled
> ODEs. In 2D, I have ~10^4, and speed is a major issue. Does anyone have
> any general tips on ways to optimally use odeint to solve large problems
> using method of lines? Some useful info/questions about my problem:

Justin,

If you could post some examples of the code that you'd like in the
right hand side of your ODEs then it will be easier to give you
advice.

-Rob
_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

[Sebastian Walter]

Well, my first guess would be that the Jacobian computation (+ its
factorizations) is the culprit since it is a 10000x10000 array for
10000 coupled ODEs.
Leaves the questions how many times the Jacobian is computed? At each
step of the integration?
Then it would be no surprise that the runtime is really high.

What kind of spatial discretization are you using to get a dense Jacobian?
Sebastian

[Justin Bois]

Thanks Ron and Sebastian for your help.  The code to calculate the RHS of the ODEs is really large so I can't really post it.  I am solving a reaction-diffusion-advection system in 2D.  The velocity (which contributes to the advection) is a complicated function of the concentrations of the species in the reaction-diffusion system.  At each time step, I have to solve a separate, time-independent ODE for the velocity, which is then incorporated into the RHS.  Because this ODE happens to be linear, I can solve it using a spectral method (for my boundary conditions of interest, regardless of my differencing scheme in the reaction-diffusion-advection system).  The presence of this velocity, manifest in the advection term, makes the Jacobian dense, even if I use low-order finite difference for my differencing.  Since I necessarily have a dense Jacobian, I typically use pseudospectral methods for computing the derivatives to get high accuracy with fewer grid points, though I have also implemented finite difference methods.

I have pasted below the first few lines of output of a profile I did using the magic function %prun with ipython -pylab.

      36447356 function calls in 1212.986 CPU seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1  903.674  903.674 1212.976 1212.976 {scipy.integrate._odepack.odeint}
   678300  157.372    0.000  239.023    0.000 basic.py:358(_raw_fftnd)
  1356600   24.225    0.000   35.504    0.000 basic.py:109(_fix_shape)

So, aside from odeint itself, most of the computational effort is spent computing derivatives and solving for the velocity, all of which is just using 2D FFTs.  Curiously, odeint spends a lot of time doing something that does not require evaluation of the RHS of the ODEs.  As I watch the integration proceed, it takes several time steps very rapidly, and then pauses at a particular time point, presumably to adjust the step size.  My guess is that Sebastian is right: this is slow because the Jacobian is large and dense.  Do you think just doing an explicit RK4 time stepping or something might help?  Or are there tricks we can play to get odeint to roll along more quickly?  I am passing an object with all sorts of attributes as an argument to the RHS function, though only some are used.  I am assuming the RHS function is all executed in Python and then the return values are all that undergo type conversion, right?  I'm guessing it's not so much type conversion that's the big slowdown, or do you think I'm mistaken?  odeint has been a nice little workhorse in my work so far!