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

[Richard Guenther]

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.

#include "tramp.h"
#include "nbody.h"


namespace NBody {

  NBVArray_t nb_pos(Interval<1>(0));
  NBVArray_t nb_vel(Interval<1>(0));
  NBSArray_t nb_mass(Interval<1>(0));
  std::vector<int> nb_int;
  int nb_n = 0;

  void nb_add(const Vector<Dim, double> &pos,
              const Vector<Dim, double> &vel,
              double mass, bool integrate)
  {
    nb_pos.create(1);
    nb_vel.create(1);
    nb_mass.create(1);
    NBVArray_t::Base_t x = nb_pos.array();
    NBVArray_t::Base_t v = nb_vel.array();
    NBSArray_t::Base_t m = nb_mass.array();
  
    x(nb_n) = pos;
    v(nb_n) = vel;
    m(nb_n) = mass;
    nb_int.push_back(integrate);
    nb_n++;
  }

  template <class X>
  static void nb_force(const NBVArray_t::Base_t &f,
                       const X &x,
                       const NBSArray_t::Base_t &m)
  {
    Interval<1> I(f.domain());
    f = 0.0;
    for (int i = 0; i<I.size(); ++i) {
      if (!nb_int[i])
        continue;
      for (int j = 0; j<I.size(); ++j)
        if (i != j)
          f(i) += (x.read(j)-x.read(i))/(norm(x.read(j)-x.read(i)))
            * m.read(j)/norm2(x.read(j)-x.read(i));
    }
  }

  /// RK4
  void nb_integrate(double dt)
  {
    NBVArray_t::Base_t x = nb_pos.array();
    NBVArray_t::Base_t v = nb_vel.array();
    NBSArray_t::Base_t m = nb_mass.array();

    Interval<1> I(x.domain());
    NBVArray_t::Base_t f(I); // temporary for forces
    NBVArray_t::Base_t k1x(I), k2x(I), k3x(I), k4x(I);
    NBVArray_t::Base_t k1v(I), k2v(I), k3v(I), k4v(I);

    nb_force(f, x, m);
    k1x = v*dt;
    k1v = f*dt;

    nb_force(f, x+k1x/2, m);
    k2x = (v+k1v/2)*dt;
    k2v = f*dt;

    nb_force(f, x+k2x/2, m);
    k3x = (v+k2v/2)*dt;
    k3v = f*dt;

    nb_force(f, x+k3x, m);
    //k4x = (v+k3v)*dt;
    //k4v = f*dt;

    x += k1x/6 + k2x/3 + k3x/3 + (v+k3v)*dt/6;
    v += k1v/6 + k2v/3 + k3v/3 + f*dt/6;
  }

}

// -*- C++ -*-

#ifndef TRAMP_NBODY_H
#define TRAMP_NBODY_H

#include "Pooma/DynamicArrays.h"


namespace NBody {

  typedef DynamicArray<Vector<Dim, double> > NBVArray_t;
  typedef DynamicArray<double> NBSArray_t;

  extern NBVArray_t nb_pos;
  extern NBVArray_t nb_vel;
  extern NBSArray_t nb_mass;
  extern std::vector<int> nb_int;
  extern int nb_n;

  void nb_add(const Vector<Dim, double> &pos,
              const Vector<Dim, double> &vel,
              double mass, bool integrate = true);

  void nb_integrate(double dt);

}

#endif