[petsc-users] Is it still worth switching to PETSc if I can't write a Jacobian for my problem?

[Brian Merchant]

Hi all,

I am considering using petsc4py instead of scipy.integrate.odeint (which is
a wrapper for Fortran solvers) for a problem involving the solution of a
system of ODEs. The problem has the potential to be stiff. Writing down its
Jacobian is very hard.

So far, I have been able to produce reasonable speed gains by writing the
RHS functions in "something like C" (using either numba or Cython). I'd
like to get even more performance out, hence my consideration of PETSc.

Due to the large number of equations involved, it is already tedious to
think about writing down a Jacobian. Even worse though, is that some of the
functions governing a particular interaction do not have neat analytical
forms (let alone whether or not their derivatives have neat analytical
forms), so we might have a mess of piecewise functions needed to
approximate them if we were to go about still trying to produce a
Jacobian...

All the toy examples I see of PETSc time stepping problems have Jacobians
defined, so I wonder if I would even get a speed gain going from switching
to it, if perhaps one of the reasons why I have a high computational cost
is due to not being able to provide a Jacobian function?

I described the sort of problem I am working with in more detail in this
scicomp.stackexchange question, which is where most of this question is
copied from, except it also comes with a toy version of the problem I am
dealing with:
http://scicomp.stackexchange.com/questions/21501/is-it-worth-switching-to-timesteppers-provided-by-petsc-if-i-cant-write-down-a

All your advice would be most helpful :)

Kind regards,Brian
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20151207/25b8816d/attachment.html>