On Thu, Dec 10, 2015 at 3:47 PM, Brian Merchant <[email protected]> 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. > > 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. > The memory is always in one place. The PETSc Vec and NumPy array are just two ways of looking into it. Matt > 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 >>>>>> >>>>>> >>>>> >>>>> >>>> >>>> >>> > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
