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

[Brian Merchant]

Hi Barry,

Points 1) and 2) are bang on -- you understand the situation very well.

1) with regards to being able to apply the RHS efficiently in C -- I would
be happy to do this, but I worry about how to cleanly get C arrays back
into NumPy arrays -- NumPy arrays are easy to manipulate for animation, and
general data analysis and visualization. Still, this is not so bad, because
all it requires is some learning from my part, since I know that NumPy does
have a C API.

2) b) Yes, the functions are discontinuous, in particular because some
parameters can jump in value from 0 to non-zero when certain conditions are
met ("is point outside of of boundary polygon? then turn on extra
inhibition").

a) I think if I ignored such discontinuities and just considered a simpler
problem where I had to write the derivative of r_n (one of the chemistry
variables at node n, whose inhibition is related to the magnitude of strain
in the edges connecting to the node in question) with respect to x_n (the
x-direction position variable of node n), things begin to get tricky and
complicated, manual differentiation wise.

Then, we also have to consider how r_n changes as x_(n+k) changes (how does
r at node n change, with respect to the change in position of another node,
n+k, which would affect the strain at this node), and things get even more
cumbersome and annoying -- at least in my head. Someone mentioned
"automatic differentiation", but I don't know enough about that to say
anything more -- is that the same as computing the Jacobian using finite
differences? Some quick Wikipedia reading seems to say "no".


> You could use explicit schemes with a small enough timestep to satisfy
any stability condition and forget about using implicit schemes. Then it
becomes crucial to have a very fast right hand side function (since you
will need many time steps). If you are lucky you can use something like a
4th order RK scheme (but again maybe with a non smooth right hand side may
you can't).

In the literature, this seems to be common approach for those who are only
modelling space variables at node n, and not any chemistry. I have a
feeling (although this is probably because I am a novice at numerical
analysis) that calculating some bound on the timesteps below which
stability is guaranteed would be very tough for the full problem including
space variables and chemistry variables. I suppose I could skip the math
and just try to reduce the timestep until things seem to work in comparison
with the results from an implicit scheme, and then decide "good enough!"?

How robust would that be though? If I have to repeat the comparison with
results from an implicit method quite often, then it may not be worth it?


Kind regards,
Brian








On Thu, Dec 10, 2015 at 11:51 AM, Barry Smith <bsmith at mcs.anl.gov> wrote:

>
>    Brian,
>
>      I see two distinct issues here
>
> 1) being apply to apply your right hand side efficiently and
>
> 2) what type of ODE integrators, if any, can work well for your problem
> with its funky, possibly discontinuous right hand side?
>
> 1)  Looking at the simplicity of your data structure and function
> evaluation I think you should just write your right hand side functions in
> C. The code would be no more complicated than it is now, from what you
> showed me. Even with "just in time" compilation likely you lose at least a
> factor of 10 by coding in python, maybe 100? I don't know. It looks to me
> like an easy translation of your current routines to C.
>
> 2) This is tricky. You a) cannot compute derivatives? and b) the function
> and derivatives may not be smooth?
>
>     If the function was well behaved you could use finite differences to
> compute the Jacobian (reasonably efficiently, the cost is just some number
> of function evaluations) but with a flaky right hand side function the
> finite differences can produce garbage.
>
>     You could use explicit schemes with a small enough timestep to satisfy
> any stability condition and forget about using implicit schemes. Then it
> becomes crucial to have a very fast right hand side function (since you
> will need many time steps). If you are lucky you can use something like a
> 4th order RK scheme (but again maybe with a non smooth right hand side may
> you can't).
>
>    I am no expert. Perhaps Emil who is far more knowledgable about these
> things has questions?
>
>    Barry
>
>
>
> > On Dec 10, 2015, at 12:52 PM, Brian Merchant <bhmerchant at gmail.com>
> wrote:
> >
> > Hi Barry,
> >
> > Here's some non-trivial example code:
> https://gist.github.com/bmer/2af429f88b0b696648a8
> >
> > I have still made some simplifications by removing some phase variables,
> expanding on variable names in general, and so on.
> >
> > The rhs function itself is defined on line 578. The functions referred
> to within it should be all defined above, so you can have a peek at them as
> necessary.
> >
> > Starting from line 634 I show how I use the rhs function. In particular,
> note the "disjointed" evaluation of the integral -- I don't just evaluate
> from 0 to t all at one go, but rather evaluate the integral in pieces
> (let's call the time spent between the end of one integral evaluation, and
> the start of the next integral evaluation a "pause"). This is so that if
> there were multiple amoebas, during the "pause", I can take into account
> changes in some of the parameters due to contact between one amoeba and
> another -- poor man's simplification.
> >
> > Please let me know if this is what you were looking for. I wouldn't be
> surprised if it wasn't, but instead would be happy to try to rework what
> I've got so it's more in line with what would be meaningful to you.
> >
> > Kind regards,
> > Brian
> >
> > On Wed, Dec 9, 2015 at 2:18 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> >   I prefer the actual code, not the mathematics or the explanation
> >
> > > On Dec 9, 2015, at 3:42 PM, Brian Merchant <bhmerchant at gmail.com>
> wrote:
> > >
> > > Hi Barry,
> > >
> > > > Could send an example of your "rhs" function; not a totally trivial
> example
> > >
> > > Sure thing! Although, did you check out the exam I tried to build up
> in this stackexchange question, along with a picture:
> http://scicomp.stackexchange.com/questions/21501/is-it-worth-switching-to-timesteppers-provided-by-petsc-if-i-cant-write-down-a
> > >
> > > I ask because that's probably the best I can do without using as
> little math as possible.
> > >
> > > Otherwise, what I'll do is take a couple of days to carefully look at
> my work, and write up a non-trivial example of a difficult-to-differentiate
> RHS, that still is a simplification of the whole mess -- expect a one or
> two page PDF?
> > >
> > > Kind regards,
> > > Brian
> > >
> > > On Mon, Dec 7, 2015 at 12:45 PM, Barry Smith <bsmith at mcs.anl.gov>
> wrote:
> > >
> > >    Brian,
> > >
> > >     Could send an example of your "rhs" function; not a totally
> trivial example
> > >
> > >    Barry
> > >
> > > > On Dec 7, 2015, at 11:21 AM, Brian Merchant <bhmerchant at gmail.com>
> wrote:
> > > >
> > > > 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/20151210/d4778374/attachment-0001.html>