DAPPLE Example: nbody
Source code:
// @TITLE "nbody.cc -- parallel nbody solution"
//
// nbody -- a parallel program to solve the n-bodies problem (in 2 dimensions)
// written using the Data-Parallel Programming Library for Education (DAPPLE)
//
// This program has lots of parallelism. It is an n^2 algorithm, in that
// each body must compute the force on itself caused by all the other n bodies.
// So, we loop over the n nodies, and in each iteration we compute the
// forces caused by the other bodies, in parallel. Then, we update the
// position and velocity of all n bodies in parallel, using the forces
// from the previous step. Over all this we iterate through timesteps.
//
// The interesting parallel parts are in ComputeForces and MoveBodies.
//
// David Kotz 1994
// $Id: nbody.cc,v 1.3 94/09/29 12:18:10 dfk Exp Locker: dfk $
#include
#include
#include
#include "dapple.h"
void ComputeForces(const int N, const doubleVector& Mass,
const doubleVector& Xpos, const doubleVector& Ypos,
const doubleVector& Xvel, const doubleVector& Yvel,
doubleVector& Xforce, doubleVector& Yforce);
void MoveBodies(const int N, const doubleVector& Mass,
doubleVector& Xpos, doubleVector& Ypos,
doubleVector& Xvel, doubleVector& Yvel,
const doubleVector& Xforce, const doubleVector& Yforce,
double& time, const double delta_t);
void PrintBodies(const int N, const doubleVector& Mass,
const doubleVector& Xpos, const doubleVector& Ypos,
const doubleVector& Xvel, const doubleVector& Yvel,
const double time);
boolean ReadBodies(FILE *fp, const int N, const doubleVector& Mass,
const doubleVector& Xpos, const doubleVector& Ypos,
const doubleVector& Xvel, const doubleVector& Yvel);
/////////////////////////////////////////////////////////////
// @SUBTITLE "Main program"
int
main(int argc, char **argv)
{
char *filename; // name of input file
double delta_t; // timestep
int steps; // number of timesteps
int N; // number of bodies
if (argc != 4) {
fprintf(stderr, "usage: nbody filename delta_t #steps\n");
exit(1);
}
filename = argv[1];
delta_t = atof(argv[2]);
if (delta_t <= 0) {
fprintf(stderr, "nbody: timestep (%f) should be > 0\n", delta_t);
exit(-1);
}
steps = atoi(argv[3]);
if (steps <= 0) {
fprintf(stderr, "nbody: number of timesteps (%d) should be > 0\n", steps);
exit(-1);
}
// open the file
FILE *fp;
if ((fp = fopen(filename, "r")) == (FILE *)NULL) {
fprintf(stderr, "nbody: cannot open input file '%s'\n", filename);
exit(1);
}
// read the number of bodies
if (fscanf(fp, "%d\n", &N) != 1) {
fprintf(stderr, "nbody: file '%s': missing value of N on first line\n",
filename);
fclose(fp);
exit(1);
}
// allocate vectors to represent the state of N bodies
doubleVector Mass(N); // mass in kilograms
doubleVector Xpos(N), Ypos(N); // x,y position in meters
doubleVector Xvel(N), Yvel(N); // x,y velocity in meters/second
doubleVector Xforce(N), Yforce(N); // x,y force components in Newtons
double time = 0; // the current time in seconds
// read in the N bodies.
if (! ReadBodies(fp, N, Mass, Xpos, Ypos, Xvel, Yvel)) {
fclose(fp);
exit(1);
} else
fclose(fp);
////////////////////////////////////////////////////
// OK, time to do the simulation
for (int i = 0; i < steps; i++) {
ComputeForces(N, Mass, Xpos, Ypos, Xvel, Yvel, Xforce, Yforce);
MoveBodies(N, Mass, Xpos, Ypos, Xvel, Yvel, Xforce, Yforce, time, delta_t);
PrintBodies(N, Mass, Xpos, Ypos, Xvel, Yvel, time);
}
return(0);
}
/////////////////////////////////////////////////////////////
// @SUBTITLE "ReadBodies: read an input file"
// Data file format:
// N
// body 0
// body 1
// body 2
// ...
// body N-1
//
// where each body is 5 floats on a line, separated by spaces:
// mass, x, y, x velocity, y velocity
// units: kg m m m/s m/s
//
boolean // true if success
ReadBodies(FILE *fp, const int N, const doubleVector& Mass,
const doubleVector& Xpos, const doubleVector& Ypos,
const doubleVector& Xvel, const doubleVector& Yvel)
{
for (int i = 0; i < N; i++) {
// 5 floats for each body: mass, x, y, x velocity, y velocity
if (fscanf(fp, "%lf %lf %lf %lf %lf\n",
&Mass[i], &Xpos[i], &Ypos[i], &Xvel[i], &Yvel[i]) != 5) {
fprintf(stderr, "nbody: bad or missing data for body %d\n", i);
return(false);
}
if (Mass[i] <= 0) {
fprintf(stderr, "nbody: zero or negative mass for body %d\n", i);
return(false);
}
}
return(true);
}
/////////////////////////////////////////////////////////////
// @SUBTITLE "ComputeForces: compute the total force on each body"
// The total force on each body is the sum of the forces
// on it by all other bodies.
// We compute the force for each body sequentially, but we can
// compute, in parallel, the force on it by all other bodies.
void ComputeForces(const int N, const doubleVector& Mass,
const doubleVector& Xpos, const doubleVector& Ypos,
const doubleVector& Xvel, const doubleVector& Yvel,
doubleVector& Xforce, doubleVector& Yforce)
{
// distance between one body and each other body
doubleVector dx(N), dy(N), d(N); // x and y components, and total distance
doubleVector f(N); // scalar component of force
const intVector VP(N, Identity); // 0, 1, 2, ..., N-1
// gravitational constant G: negative since it is attractive
// Units: Nm^2/kg^2; thus force is in Newtons
const double GRAVITY = -6.670e-11;
for (int body = 0; body < N; body++)
ifp (VP != body) {
// find the distance to all other bodies
dx = Xpos[body] - Xpos;
dy = Ypos[body] - Ypos;
d = apply(sqrt, dx*dx + dy*dy);
// find the scalar component of the force caused by all other bodies
f = GRAVITY * Mass[body] * Mass / (d * d * d);
// sum all the forces in both the x and y directions
Xforce[body] = sum(dx * f);
Yforce[body] = sum(dy * f);
}
}
/////////////////////////////////////////////////////////////
// @SUBTITLE "MoveBodies: update the position/velocity from accumulated force"
void MoveBodies(const int N, const doubleVector& Mass,
doubleVector& Xpos, doubleVector& Ypos,
doubleVector& Xvel, doubleVector& Yvel,
const doubleVector& Xforce, const doubleVector& Yforce,
double& time, const double delta_t)
{
doubleVector ax(N), ay(N); // the acceleration due to the force
// a = F/m
ax = Xforce / Mass;
ay = Yforce / Mass;
// x' = x + vt + 1/2 at^2
Xpos = Xpos + Xvel * delta_t + 0.5 * ax * delta_t * delta_t;
Ypos = Ypos + Yvel * delta_t + 0.5 * ay * delta_t * delta_t;
// v' = v + at
Xvel = Xvel + ax * delta_t;
Yvel = Yvel + ay * delta_t;
// t' = t + dt
time += delta_t;
}
/////////////////////////////////////////////////////////////
// @SUBTITLE "PrintBodies: print the state of all bodies"
void PrintBodies(const int N, const doubleVector& Mass,
const doubleVector& Xpos, const doubleVector& Ypos,
const doubleVector& Xvel, const doubleVector& Yvel,
const double time)
{
printf("State of the %d bodies at time %g seconds:\n", N, time);
printf("%4s %12s %12s %12s %12s %12s\n",
"Body", "Mass", "X", "Y", "Xvel", "Yvel");
for (int body = 0; body < N; body++)
printf("%4d %12g %12g %12g %12g %12g\n",
body, Mass[body], Xpos[body], Ypos[body], Xvel[body], Yvel[body]);
printf("\n");
fflush(stdout);
}
Demonstration:
cat nbody.data/10
10
417053 18 27 0 0
100785 46 83 0 0
85307 19 81 0 0
214703 51 26 0 0
401032 25 2 0 0
97961 65 51 0 0
60540 75 77 0 0
401469 57 64 0 0
466678 51 19 0 0
416242 63 15 0 0
nbody nbody.data/10 1000 10
State of the 10 bodies at time 1000 seconds:
Body Mass X Y Xvel Yvel
0 417053 18.0363 26.9816 7.25803e-05 -3.68819e-05
1 100785 46.0126 82.9595 2.52204e-05 -8.10575e-05
2 85307 19.0168 80.9834 3.36482e-05 -3.31017e-05
3 214703 51.0212 25.6547 4.23873e-05 -0.00069054
4 401032 25.0223 2.03956 4.45635e-05 7.91256e-05
5 97961 64.9523 51.0192 -9.54805e-05 3.83001e-05
6 60540 74.966 76.9689 -6.80892e-05 -6.21225e-05
7 401469 56.9999 63.9735 -2.71262e-07 -5.30206e-05
8 466678 51.0608 19.1246 0.000121562 0.000249299
9 416242 62.8729 15.0575 -0.000254179 0.000114973
State of the 10 bodies at time 2000 seconds:
Body Mass X Y Xvel Yvel
0 417053 18.1453 26.9261 0.000145359 -7.39599e-05
1 100785 46.0505 82.8378 5.04948e-05 -0.000162245
2 85307 19.0673 80.9338 6.73167e-05 -6.62347e-05
3 214703 51.089 24.5702 9.32955e-05 -0.00147844
4 401032 25.0892 2.15838 8.92047e-05 0.000158515
5 97961 64.8089 51.077 -0.000191184 7.74785e-05
6 60540 74.8637 76.8757 -0.000136326 -0.000124416
7 401469 56.9995 63.8939 -4.33298e-07 -0.000106241
8 466678 51.244 19.5189 0.000244783 0.000539207
9 416242 62.4884 15.2323 -0.000514937 0.000234627
State of the 10 bodies at time 3000 seconds:
Body Mass X Y Xvel Yvel
0 417053 18.3273 26.8333 0.000218717 -0.000111635
1 100785 46.1137 82.6348 7.59331e-05 -0.000243828
2 85307 19.1515 80.8509 0.000101048 -9.94633e-05
3 214703 51.2319 22.4494 0.000192491 -0.00276315
4 401032 25.2008 2.35699 0.00013411 0.000238702
5 97961 64.5695 51.1755 -0.000287628 0.000119352
6 60540 74.6931 76.7198 -0.000205011 -0.000187228
7 401469 56.9992 63.7607 -2.67495e-07 -0.000160103
8 466678 51.5484 20.3103 0.000364148 0.00104366
9 416242 61.8323 15.5346 -0.00079724 0.000370004
State of the 10 bodies at time 4000 seconds:
Body Mass X Y Xvel Yvel
0 417053 18.5831 26.7024 0.000292962 -0.000150329
1 100785 46.2025 82.3498 0.000101652 -0.000326098
2 85307 19.2695 80.7348 0.000134892 -0.000132855
3 214703 51.9641 16.353 0.00127192 -0.00942969
4 401032 25.3577 2.63645 0.000179575 0.000320222
5 97961 64.2329 51.3181 -0.000385565 0.000165913
6 60540 74.4533 76.5008 -0.000274469 -0.000250933
7 401469 56.9993 63.573 4.48759e-07 -0.00021518
8 466678 51.7654 22.8298 6.98913e-05 0.00399521
9 416242 60.8721 15.9886 -0.00112301 0.000537924
State of the 10 bodies at time 5000 seconds:
Body Mass X Y Xvel Yvel
0 417053 18.9136 26.5321 0.00036789 -0.000190129
1 100785 46.3172 81.982 0.000127793 -0.000409478
2 85307 19.4214 80.5851 0.000168956 -0.000166517
3 214703 53.3763 7.29308 0.00155245 -0.00869019
4 401032 25.5606 2.99804 0.000226279 0.000402946
5 97961 63.7967 51.5105 -0.000486845 0.000218875
6 60540 74.1435 76.2173 -0.000345124 -0.000316011
7 401469 57.0004 63.3291 1.8416e-06 -0.000272659
8 466678 51.9063 26.5974 0.000211854 0.00353998
9 416242 59.5469 16.6108 -0.00152749 0.000706578
State of the 10 bodies at time 6000 seconds:
Body Mass X Y Xvel Yvel
0 417053 19.3187 26.3223 0.000442458 -0.000229511
1 100785 46.4584 81.5299 0.000154537 -0.000494824
2 85307 19.6076 80.4015 0.000203531 -0.000200706
3 214703 54.963 -1.25412 0.00162094 -0.00840422
4 401032 25.8105 3.44131 0.000273469 0.000483594
5 97961 63.2551 51.759 -0.000596344 0.000278242
6 60540 73.7621 75.8677 -0.000417726 -0.000383225
7 401469 57.0031 63.0251 3.6117e-06 -0.000335478
8 466678 52.1552 30.0572 0.000285942 0.00337977
9 416242 57.9106 17.3564 -0.001745 0.000784519
State of the 10 bodies at time 7000 seconds:
Body Mass X Y Xvel Yvel
0 417053 19.7984 26.0735 0.000516958 -0.000268035
1 100785 46.6267 80.991 0.00018207 -0.000582954
2 85307 19.8288 80.1833 0.000238859 -0.000235658
3 214703 56.5683 -9.59193 0.00158956 -0.00827139
4 401032 26.1057 3.96416 0.000316985 0.000562106
5 97961 62.5976 52.0709 -0.000718719 0.000345537
6 60540 73.3068 75.4494 -0.000493007 -0.000453391
7 401469 57.0078 62.6539 5.65991e-06 -0.00040685
8 466678 52.4559 33.376 0.000315501 0.00325774
9 416242 56.1108 18.2024 -0.00185466 0.000907486
State of the 10 bodies at time 8000 seconds:
Body Mass X Y Xvel Yvel
0 417053 20.353 25.7863 0.000592129 -0.000306366
1 100785 46.823 80.3622 0.000210595 -0.000674728
2 85307 20.0858 79.9297 0.000275156 -0.000271586
3 214703 58.1419 -17.8241 0.00155777 -0.00819296
4 401032 26.4416 4.56641 0.000354791 0.000642397
5 97961 61.8084 52.4558 -0.000859759 0.000424238
6 60540 72.7744 74.959 -0.000571744 -0.000527432
7 401469 57.0146 62.2049 7.99049e-06 -0.000491129
8 466678 52.7728 36.5898 0.000318267 0.0031699
9 416242 54.221 19.171 -0.00192505 0.00102968
State of the 10 bodies at time 9000 seconds:
Body Mass X Y Xvel Yvel
0 417053 20.9835 25.4603 0.000668988 -0.000345774
1 100785 47.0485 79.6392 0.000240344 -0.000771138
2 85307 20.3797 79.6395 0.000312638 -0.000308693
3 214703 59.6881 -25.9901 0.00153447 -0.00813897
4 401032 26.8129 5.25183 0.000387821 0.000728441
5 97961 60.8645 52.9286 -0.00102789 0.000521273
6 60540 72.1611 74.392 -0.000654813 -0.000606445
7 401469 57.0239 61.6617 1.0595e-05 -0.00059522
8 466678 53.0866 39.7338 0.000309342 0.0031181
9 416242 52.2696 20.254 -0.00197761 0.00113649
State of the 10 bodies at time 10000 seconds:
Body Mass X Y Xvel Yvel
0 417053 21.6925 25.0934 0.000748994 -0.000387896
1 100785 47.3044 78.8169 0.000271581 -0.000873415
2 85307 20.7118 79.3116 0.000351536 -0.000347191
3 214703 61.2143 -34.1091 0.00151803 -0.00809913
4 401032 27.2156 6.0277 0.000417599 0.0008233
5 97961 59.7324 53.5149 -0.00123646 0.000651398
6 60540 71.4621 73.7429 -0.000743259 -0.000691784
7 401469 57.0358 60.9983 1.31568e-05 -0.000731659
8 466678 53.3915 42.8461 0.000300341 0.00310648
9 416242 50.2689 21.4352 -0.00202394 0.00122582