[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Freepooma-devel] ODE System of arbitrary length
From: |
Richard Guenther |
Subject: |
Re: [Freepooma-devel] ODE System of arbitrary length |
Date: |
Thu, 25 Aug 2005 20:55:13 +0200 |
User-agent: |
Debian Thunderbird 1.0.6 (X11/20050802) |
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