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. 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 <[email protected]> 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 <[email protected]> >>> 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 <[email protected]> wrote: >>> >>> I prefer the actual code, not the mathematics or the explanation >>> >>> On Dec 9, 2015, at 3:42 PM, Brian Merchant <[email protected]> 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 <[email protected]> >>>> 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 <[email protected]> >>>>> 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 >>>>> >>>>> >>>> >>>> >>> >>> >>
