Re: Re: [Freepooma-devel] ODE System of arbitrary length

[Richard Guenther]

On Fri, 26 Aug 2005, Frank Gollas wrote:

> Richard Guenther wrote:
> >Harris wrote:
> >>Hi,
> >>
> >>I just stumbled on POOMA and I think it's pretty
> >>cool and better suited for my needs than blitz++.
> >>
> >>What I want to do is a numerical 4 step Runge-Kutta
> >>integration of an ODE system
> >>
> >>dx1/dt = f1(x1,x2, x3 ..., Y )
> >>dx2/dt = f2(x1,x2, x3 ..., Y )
> >>dx3/dt = f3(x1,x2, x3 ..., Y )
> >>...
> >>
> >>this should be rather straight forward.
> >>But: the number of state equations is only
> >>known at _runtime_, the class that does
> >>the integration holds a dynamical list
> >>of arrays. Is there an efficient
> >>way of doing this using stencils?
> >
> >I don't exactly understand what you mean by "dynamical list
> >of arrays", but I can show you the code-snipped I'm using to
> >integrate a gravitating n-Body system in time using a 4th order
> >Runge-Kutta (attached).
> >
> >Basically, you need to use DynamicArray to store the state (and
> >temporaries) for a variable number of bodies.  The integration
> >part is data-parallel, while the force calculation (which is
> >brute force n^2 algorithm) is using manual loops.
> >
> >If I understand you correctly, you want to use a stencil for
> >calculating the integration in time.  This is possible, if
> >your "force" can be calculated locally.  Maybe the code
> >snipped helps answering more questions, otherwise just ask.
> >
> >Richard.
>
>
> thank you very much for your fast answer and the code snipsets that were
> pretty helpful. I'm also attaching some of my testing code to make clear my
> intention. x is a container of Array's, each representing a state variable
> field of the ODE system. f1 to fn are represented by stencils. Unfortunatly
> I didn't manage to express f1(x1,x2,x3 ... ), f2(x1,x2,x3 ... ) in POOMA. In
> my code only one equation f1(x[i]) with one array as parameter is applied.
>
> The goal is to have as many Stencils as state equations in the ODE System,
> each taking all arrays representing the state variables (and some constant
> input ..) as argument, such that interaction between different state
> variables can be implemented.

(Sorry for the late response)

The only way to do this is to use the ScalarCode type of "stencils"
where you use the MultiArg<> object to encapsulate multiple arguments
by function overloading.  You need to provide overloads for different
number of arguments in your stencils operator(), too, of course.

But it doesn't play nicely with std::vector (aka runtime variable
number of arguments).  It's all to be set statically at compile time.
The only way around it is to use DynamicArray, which I don't know if
you can easily use with Stencil or ScalarCode.  Just try doing
DynamicArray< Array<2> > instead of using std::vector.


> Is it possible to do this with layers? Does my code look ok up to this
> point? BTW the Laplace Stencil example did compute nonsense until I replaced
> "world(interior_2) = stencil(world);" with "temp(interior_2) =
> stencil(world); world = temp;" what might be the reason for that?

yeah, world(interior_2) = stencil(world) is broken.  FreePOOMA doesn't
allow aliasing of lhs/rhs.

Richard.

--
Richard Guenther <richard dot guenther at uni-tuebingen dot de>
WWW: http://www.tat.physik.uni-tuebingen.de/~rguenth/