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