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

[Brian Merchant]

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 <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>
>>> 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/25e0ca49/attachment-0001.html>