[petsc-users] Is it still worth switching to PETSc if I can't write a Jacobian for my problem?
Emil Constantinescu
emconsta at mcs.anl.gov
Thu Dec 10 16:22:24 CST 2015
On 12/10/15 3:47 PM, Brian Merchant wrote:
> Hi Emil:
>
> I had a look at odeint's output -- it seems that while Adams is used
> initially, bdf is used most often overall (especially in later
> timesteps), some 80% of the timesteps in total.
>
> This is likely because I ran the test on the full problem, which
> includes interactions between multiple amoebas -- contact between
> amoebas turns certain parameters on temporarily.
Hmm, so it may be that the problem is not stiff, but stepping over
discontinuities makes it appear stiff and hence the switch to BDF or it
may be that the problem has periods when it's really stiff. Like Hong
said, by using PETSc as the time stepper you may be able to diagnose
better what's going on b/c you have a choice of using many integrators
and even events should you try to address those discontinuities directly.
Emil
> Matt: is this the getArray method of petsc4py? It converts a PETSc array
> to a NumPy array? How did the PETSc array (let's say from a C program)
> get loaded by petsc4py in the first place -- is there some sort of a
> "load" function? I am sure there are straightforward answers to this
> question, but these are the sorts of things I don't know right now.
>
>
> Kind regards,
> Brian
>
> On Thu, Dec 10, 2015 at 12:55 PM, Emil Constantinescu
> <emconsta at mcs.anl.gov <mailto:emconsta at mcs.anl.gov>> wrote:
>
> Brian,
>
> Can you take a look at what odeint returns? Specifically, at:
>
> ‘mused’ a vector of method indicators for each successful
> time step: 1: adams (nonstiff), 2: bdf (stiff)
>
> I suspect it's using Adams all the way, which means it's doesn't
> need a Jacobian.
>
> Emil
>
>
>
> On 12/10/15 1:51 PM, Barry Smith 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 <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto: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
>
>
>
>
>
>
>
More information about the petsc-users
mailing list