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

[Matthew Knepley]

On Thu, Dec 10, 2015 at 2:56 PM, Brian Merchant <bhmerchant at gmail.com>
wrote:

> 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.
>

v.getArray() gives back a NumPy array.

  Thanks,

    Matt


> 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
>> > > >
>> > >
>> > >
>> >
>> >
>>
>>
>


-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20151210/ee78d8ba/attachment.html>