1 #include <omp.h> 2 /******************************************************************************/ 3 /* */ 4 /* Copyright (c) 1996 University of California at Santa Barbara */ 5 /* */ 6 /* All rights reserved. */ 7 /* */ 8 /* Permission is given to use, copy, and modify this software for any */ 9 /* non-commercial purpose as long as this copyright notice is not removed */ 10 /* All other uses, including redistribution in whole or in part, are */ 11 /* forbidden without prior written permission. */ 12 /* */ 13 /* This software is provided with absolutely no warranty and no support. */ 14 /* */ 15 /* For further information contact: pedro@cs.ucsb.edu or martin@cs.ucsb.edu */ 16 /* */ 17 /******************************************************************************/ 18 19 #include <stdio.h> 20 #include <string.h> 21 #include <math.h> 22 #include <stdlib.h> 23 #include <iostream> 24 #include <fstream> 25 #include <limits.h> 26 27 #define CHECK 28 using namespace std; 29 30 // const int NULL = 0; 31 // const int FALSE = 0; 32 // const int TRUE = 1; 33 34 const double PI = 3.14159265358979323846; 35 // const double TWO_PI = 6.28318530717958647693; 36 // const double FOUR_PI = 12.56637061435917295385; 37 // const double HALF_PI = 1.57079632679489661923; 38 // const double FRTHRD_PI = 4.18879020478639098462; 39 40 const int MAX_NAME = 80; 41 42 extern "C" void reset_stat(); 43 extern "C" void print_stat(); 44 45 extern "C" void get_time(int *r); 46 extern "C" void get_ticks_time(int *r); 47 48 extern "C" double prand(); 49 extern "C" double xrand(double xl, double xh); 50 extern "C" void pranset(int seed); 51 52 extern "C" int scanbind(char **bvec, char *name); 53 extern "C" int matchname(char *bind, char *name); 54 extern "C" char* extrvalue(char *arg); 55 #define _OPENMP_ 56 #include "water.h" 57 //const int NDIM = 3; 58 inline void in_vector(ifstream& f, double *v){ 59 double x1, x2, x3; 60 f >> x1 >> x2 >> x3; 61 v[0] = x1; v[1] = x2; v[2] = x3; 62 } 63 64 inline void tvecClr(double *v) { int i; for(i=0; i < NDIM; i++) v[i] = 0.0; } 65 66 inline void tvecAdd(double *v1, double *v2) { int i; 67 for(i=0; i < NDIM; i++) v1[i] += v2[i]; } 68 69 inline void tvecSub(double *v1, double *v2) { int i; 70 for(i=0; i < NDIM; i++) v1[i] -= v2[i]; } 71 72 inline void tvecLoad(double *v1, double *v2) { int i; 73 for(i=0; i < NDIM; i++) v1[i] = v2[i]; } 74 75 inline void tvecStore(double *v1, double *v2) { int i; 76 for(i=0; i < NDIM; i++) v2[i] = v1[i]; } 77 78 inline void tvecDiv(double *v1, double *v2) { int i; 79 for(i=0; i < NDIM; i++) v1[i] /= v2[i]; } 80 81 inline void tvecProd(double *v1, double *v2) { int i; 82 for(i=0; i < NDIM; i++) v1[i] *= v2[i]; } 83 84 inline void tvecScale(double *v, double s) { int i; 85 for(i=0; i < NDIM; i++) v[i] *= s; } 86 87 inline double tvecDotProd(double *v1,double *v2) { int i; double d = 0.0; 88 for(i=0; i < NDIM; i++) d += (v1[i] * v2[i]); return d; } 89 90 inline void tvecCrossProd(double *v, double *u, double *w){ 91 v[0] = u[1]*w[2] - u[2]*w[1]; 92 v[1] = u[2]*w[0] - u[0]*w[2]; 93 v[2] = u[0]*w[1] - u[1]*w[0]; 94 } 95 96 inline void tvecMin(double *v1, double *v2){ int i; 97 for(i=0; i < NDIM; i++) if(v1[i] > v2[i]) v1[i] = v2[i]; } 98 99 inline double tvecNorm1(double *v){ int i; double d; 100 d = 0.0; for(i=0; i < NDIM; i++) d += v[i]; return d; } 101 102 103 // ----------- FILE: body.h ---------- 104 // BH-tree capacity constants 105 const int LEAFMAXBODIES = 10; 106 const int NSUB = 8; // (1 << NDIM) 107 const int IMAX = 1073741824; // (1 << MAXLEVEL) 108 109 const int CELL = 01; 110 const int LEAF = 02; 111 const int BODY = 04; 112 113 class node { // inherited by all CELL, LEAF and BODY classes 114 public: 115 int type; 116 double mass; /* total mass of node */ 117 vector pos; /* position of node */ 118 int level; 119 node *parent; /* ptr to parent of this node in tree */ 120 int child_num; /* Index that this node should be put at in parent cell */ 121 node() { type = 0; mass = 0.0; pos.vecClr(); level = 0; parent = NULL; child_num =0;} 122 ~node() { } 123 int Type() { return type; } 124 double Mass() { return mass; } 125 vector *Pos() { return &pos; } 126 int Level() { return level; } 127 node *Parent() { return parent; } 128 int ChildNum() { return child_num; } 129 void setChildNum(int n) { child_num = n; } 130 void setKind(int t) { type = t; } 131 void setMass(double m) { mass = m; } 132 void setLevel(int d) { level =d; } 133 void setParent(node *p) { parent = p; } 134 void PosvecStore(double *p) { pos.vecStore(p); } 135 void PosvecLoad(double *p) { pos.vecLoad(p); } 136 void PosvecAdd(double *p) { pos.vecAdd(p); } 137 void clear() { mass= 0.0; pos.vecClr(); level = 0; 138 parent = NULL; child_num =0; } 139 void adjLevel(node *p) { 140 if (p == NULL) level = (IMAX >> 1); 141 else level = (p->Level() >> 1); 142 } 143 int subindex(int *x, int lev); 144 void computecofm(); 145 }; 146 147 // ---------------------------------------------------------------------------- 148 // BODY: data structure used to represent particles. 149 // ---------------------------------------------------------------------------- 150 class body : public node { 151 vector vel; /* velocity of body */ 152 vector acc; /* acceleration of body */ 153 double phi; /* potential at body */ 154 vector tmp; /* Skratch vector */ 155 double subdiv(node *p, double *ai); 156 public: 157 body() { type = BODY; vel.vecClr(); acc.vecClr(); tmp.vecClr(); phi = 0.0; } 158 ~body() { } 159 vector *Vel() { return &vel; } 160 vector *Acc() { return &acc; } 161 double Phi() { return phi; } 162 vector *Tmp() { return &tmp; } 163 void BodyClear() { node::clear(); 164 vel.vecClr(); acc.vecClr(); tmp.vecClr(); 165 phi = 0.0; 166 } 167 void VelvecAdd(double *p) { vel.vecAdd(p); } 168 void VelvecStore(double *p) { vel.vecStore(p); } 169 void VelvecLoad(double *p) { vel.vecLoad(p); } 170 void AccvecStore(double *p) { acc.vecStore(p); } 171 172 double computeInter(node *p, double epsSq, double *res); 173 void openCell(node *n, double tolsq, double dsq, double epssq); 174 void openLeaf(node *n, double epssq, int iter); 175 void gravsub(node *p, double epssq); 176 void updatePhi(double inc); 177 void fastgravsub(node *p, double *ai, double drsq, double epssq); 178 void walksub(node *n, double tolsq, double dsq, double epssq); 179 void advance(double hts, double ts); 180 void swapAcc(); 181 void startVel(double f); 182 #ifdef CHECK 183 void dump(); 184 #endif 185 }; 186 187 // ---------------------------------------------------------------------------- 188 // LEAF: structure used to represent leaf nodes of tree. 189 // ---------------------------------------------------------------------------- 190 class cell; 191 class leaf : public node { 192 int num_bodies; 193 body *bodyp[LEAFMAXBODIES]; // bodies of leaf 194 public: 195 leaf() { type = LEAF; } 196 ~leaf() { } 197 body *Bodyp(int idx) { return bodyp[idx]; } 198 void setBodyp(int idx, body *p) { bodyp[idx] = p; } 199 int NumBody() { return num_bodies; } 200 void setNumBody(int n) { num_bodies = n; } 201 int full() { return (num_bodies == LEAFMAXBODIES); } 202 void LeafClear() { int i; 203 node::clear(); 204 num_bodies = 0; 205 for (i = 0; i < LEAFMAXBODIES; i++) 206 bodyp[i] = NULL; 207 } 208 209 void Subdivide(cell *ce, int idx); 210 void LeafAddBody(body *bp); 211 }; 212 213 // ---------------------------------------------------------------------------- 214 // CELL: structure used to represent internal nodes of tree. 215 // ---------------------------------------------------------------------------- 216 class cell : public node { 217 node *subp[NSUB]; // descendents of cell 218 public: 219 cell() { type = CELL; } 220 ~cell() { } 221 node *Subp(int idx) { return subp[idx]; } 222 void setSubp(int idx, node *p) { subp[idx] = p; } 223 void CellClear() { int i; 224 node::clear(); 225 for (i = 0; i < NSUB; i++) 226 subp[i] = NULL; 227 } 228 void CellAddBody(body *p, int *coords); 229 }; 230 231 // ---------------------------------------------------------------------------- 232 // CELLSET: set of cells. 233 // ---------------------------------------------------------------------------- 234 class cellset { 235 int numcell; /* no. cells used in celltab */ 236 int maxcell; /* max no. of cells allocated */ 237 cell *celltab; /* array of cells allocated */ 238 public: 239 cellset(int max) { int i; 240 celltab = new cell[max]; maxcell=max; numcell=0; 241 for(i=0; i < maxcell; i++) 242 celltab[i].clear(); 243 } 244 ~cellset() { delete[] celltab; maxcell=0; numcell=0; } 245 void CellSetclear() { int i; 246 for(i=0; i < maxcell; i++) 247 celltab[i].clear(); 248 } 249 void CellSetReset() { numcell = 0; } 250 cell *newcell(); 251 cell *makecell(node *parent); 252 cell *getcell(int idx) { return &celltab[idx]; } 253 }; 254 255 // ---------------------------------------------------------------------------- 256 // LEAFSET: set of leaves. 257 // ---------------------------------------------------------------------------- 258 class leafset { 259 int numleaf; /* no. leaves used in leafctab */ 260 int maxleaf; /* max no. of leaves allocated */ 261 leaf *leaftab; /* array of leaves allocated */ 262 public: 263 leafset(int max) { leaftab = new leaf[max]; maxleaf=max; numleaf=0; LeafSetClear(); } 264 ~leafset() { delete[] leaftab; maxleaf=0; numleaf=0;} 265 void LeafSetClear() { int i; for(i=0; i < maxleaf; i++) 266 leaftab[i].clear(); 267 } 268 void LeafSetReset() { numleaf = 0; } 269 leaf* newleaf(); 270 leaf* makeleaf(node *parent); 271 leaf *getleaf(int idx) { return &leaftab[idx]; } 272 // void add_body(body *bp); 273 }; 274 275 // ---------------------------------------------------------------------------- 276 // BODYSET: set of particle. 277 // ---------------------------------------------------------------------------- 278 class bodyset { 279 int numbody; /* no. bodies used in bodytab */ 280 int maxbody; /* max no. of bodies allocated */ 281 body *bodytab; /* array of bodies allocated */ 282 public: 283 bodyset(int m) { int i; bodytab = new body[m]; maxbody=m; numbody=0; for(i=0; i < maxbody; i++) bodytab[i].clear(); } 284 ~bodyset() { delete[] bodytab; maxbody=0; numbody=0;} 285 void BodySetClear() { int i; for(i=0; i < maxbody; i++) 286 bodytab[i].clear(); 287 } 288 void BodySetReset() { numbody = 0; } 289 body* getbody(int idx) { return &bodytab[idx]; } 290 body* newbody(); 291 int size() { return numbody; } 292 }; 293 294 class parms; 295 class NBodySystem { 296 int nbody; 297 cell* BH_tree; 298 cellset* cells; 299 leafset* leaves; 300 bodyset* bodies; 301 302 acc_double mtot; 303 vector *etot; 304 matrix *keten; 305 matrix *peten; 306 vector *cmphasePos; 307 vector *cmphaseVel; 308 vector *amvec; 309 310 vector *mincorner; 311 vector *maxcorner; 312 313 double size; 314 double tnow; 315 double tout; 316 317 int nstep; 318 int n2bcalc; 319 int nbccalc; 320 int computestart; 321 int computeend; 322 int tracktime; 323 int partitiontime; 324 int treebuildtime; 325 int forcecalctime; 326 int updatetime; 327 void pickshell(double *v, double rad); 328 public: 329 NBodySystem() { nbody = 0; cells = NULL; leaves = NULL; bodies = NULL; } 330 ~NBodySystem() { delete cells; delete leaves; delete bodies; nbody = 0; } 331 void init(); 332 void intcoord(int *xp, double *rp); 333 void updateInter(int t) { if(t == BODY) n2bcalc++; else nbccalc++; } 334 bodyset *getbodies() { return bodies; } 335 leafset *getleaves() { return leaves; } 336 cellset *getcells() { return cells; } 337 void loadData(char *name); 338 void loadTestData(); 339 void clearStats(); 340 void updateStats(); 341 void outputStats(); 342 void updateEnergy(); 343 void outputEnergy(); 344 void clearTiming(); 345 void outputTiming(); 346 void initOutput(parms *p); 347 void output(parms *p); 348 void startrun(parms *p); 349 void stepsystem(parms *p); 350 void setbound(); 351 void maketree(); 352 void SwapAccs(); 353 void Advance(double hts, double ts); 354 void ComputeAccels(double tolerance, double eps); 355 void StartVels(double halftimestep); 356 void computeForces(double tol, double hts, double eps); 357 #ifdef CHECK 358 void dump(); 359 #endif 360 }; 361 362 // ----------- FILE: simulation.h ---------- 363 class parms { 364 char headline[MAX_NAME];/* message describing calculation */ 365 char infile[MAX_NAME]; /* file name for snapshot input */ 366 char outfile[MAX_NAME]; /* file name for snapshot output */ 367 double dtime; /* timestep for leapfrog integrator */ 368 double dtout; /* time between data outputs */ 369 double tstop; /* time to stop calculation */ 370 double fcells; /* ratio of cells/leaves allocated */ 371 double fleaves; /* ratio of leaves/bodies allocated */ 372 double tol; /* accuracy parameter: 0.0 => exact */ 373 double tolsq; /* square of previous */ 374 double eps; /* potential softening parameter */ 375 double epssq; /* square of previous */ 376 double dthf; /* half time step */ 377 char *defaults[12]; /* pairs of name = value */ 378 379 public: 380 parms(); 381 ~parms(){ } 382 int getiparam(char *name); 383 double getdparam(char *name); 384 void getparam(char *name); 385 386 void setparms(char **fnames, double *args); 387 void loadInfile(); 388 void loadParms(); 389 char *getheadline() { return headline; } 390 char *getinfile() { return infile; } 391 char *getoutfile() { return outfile; } 392 double getdtime() { return dtime; } 393 double getdtout() { return dtout; } 394 double gettstop() { return tstop; } 395 double getfcells() { return fcells; } 396 double getfleaves() { return fleaves; } 397 double gettol() { return tol; } 398 double gettolsq() { return tolsq; } 399 double geteps() { return eps; } 400 double getepssq() { return epssq; } 401 double getdthf() { return dthf; } 402 }; 403 404 // --- Global variables --------------------------------------------------- 405 NBodySystem Nbody; 406 parms simparms; 407 char buf[128]; 408 409 // ### CODE BELOW ################################################### 410 // ### CODE BELOW ################################################### 411 // ### CODE BELOW ################################################### 412 413 414 #ifdef CHECK 415 // ---------------------------------------------------------------------------- 416 void NBodySystem::dump(){ 417 int i; 418 body *p; 419 printf ("## mass\t (position)\t (velocity)\t (acceleration)\n"); 420 for (i = 0; i < nbody ; i++){ 421 p = bodies->getbody(i); 422 p->dump(); 423 } 424 } 425 #endif 426 427 #ifdef CHECK 428 void body::dump(){ 429 printf("%8.4f ",mass); 430 pos.vecPrint(); 431 vel.vecPrint(); 432 acc.vecPrint(); 433 printf("\n"); 434 } 435 #endif 436 437 438 439 // --------------------------------------------------------------------------- 440 void NBodySystem::init() { 441 etot = new vector; 442 keten = new matrix; 443 peten = new matrix; 444 cmphasePos = new vector; 445 cmphaseVel = new vector; 446 amvec = new vector; 447 mincorner = new vector; 448 maxcorner = new vector; 449 } 450 451 // --------------------------------------------------------------------------- 452 void NBodySystem::clearStats(){ 453 nstep = 0; 454 n2bcalc = 0; 455 nbccalc = 0; 456 } 457 458 // --- NBodySystem --------------------------------------------------- 459 void NBodySystem::intcoord(int *xp, double *rp){ 460 int k; 461 double xsc, v[NDIM]; 462 463 mincorner->vecStore(v); 464 for(k = 0; k < NDIM; k++) { 465 xsc = (rp[k] - v[k]) / size; 466 xp[k] = (int)(floor(IMAX * xsc)); 467 } 468 } 469 470 // --------------------------------------------------------------------------- 471 void NBodySystem::updateStats(){ 472 int nttot, nbavg, ncavg; 473 nttot = n2bcalc + nbccalc; 474 nbavg = (int) ((double) n2bcalc / (double) nbody); 475 ncavg = (int) ((double) nbccalc / (double) nbody); 476 } 477 478 void NBodySystem::outputStats(){ 479 cout << "\n\nINTERACTIONS\tTOTAL\tITER\tAVG\n"; 480 cout << "body-body\t\t"<< n2bcalc << "\t" << n2bcalc/nstep<< "\t" << 481 (int)sqrt((double)(n2bcalc/nstep)) << "\n"; 482 cout << "body-cell\t\t"<< nbccalc << "\t" << nbccalc/nstep << "\t" << 483 (int)sqrt((double)(nbccalc/nstep)) << " \n";; 484 } 485 486 // ---------------------------------------------------------------------------- 487 void NBodySystem::updateEnergy(){ 488 int k; 489 double m, velsq; 490 body *p; 491 492 mtot.writeval(0.0); 493 etot->vecClr(); 494 keten->matClr(); 495 peten->matClr(); 496 cmphasePos->vecClr(); 497 cmphaseVel->vecClr(); 498 amvec->vecClr(); 499 500 #pragma omp parallel for private (velsq, m, p) schedule(static) 501 for(k = 0; k < nbody; k++) { 502 p = bodies->getbody(k); 503 matrix m1; 504 double v0[3], v1[3], v2[3], v3[3], v4[3]; 505 m = p->Mass(); 506 mtot.addvalRepl(m); 507 p->VelvecStore(v1); 508 velsq = tvecDotProd(v1,v1); 509 v4[1] = 0.5 * m * velsq; 510 v4[2] = 0.5 * m * p->Phi(); 511 v4[0] = v4[1] + v4[2]; 512 etot->vecAddRepl(v4); 513 p->PosvecStore(v1); 514 p->VelvecStore(v2); 515 p->AccvecStore(v3); 516 tvecLoad(v0,v2); 517 tvecScale(v0,0.5*m); 518 m1.matOutProd(v0,v2); 519 keten->matAddRepl(&m1); 520 tvecLoad(v0,v1); 521 tvecScale(v0,m); 522 m1.matOutProd(v0,v3); 523 peten->matAddRepl(&m1); 524 tvecLoad(v0,v1); 525 tvecScale(v0,m); 526 cmphasePos->vecAddRepl(v0); 527 tvecLoad(v0,v2); 528 tvecScale(v0,m); 529 cmphaseVel->vecAddRepl(v0); 530 tvecCrossProd(v0,v1,v2); 531 tvecScale(v0,m); 532 amvec->vecAddRepl(v0); 533 } 534 535 if (mtot.readval() != 0.0){ 536 cmphasePos->vecScale(1.0/mtot.readval()); 537 cmphaseVel->vecScale(1.0/mtot.readval()); 538 } 539 } 540 541 542 // ---------------------------------------------------------------------------- 543 void NBodySystem::outputEnergy(){ 544 cout << "\n\nENERGY TOTALS\n"; 545 cout << " MASS: "<< mtot.readval() << "\n"; 546 cout << " ENERGY: "; etot->vecPrint(); cout << "\n"; 547 cout << " KIN ENERGY: "; keten->matPrint(); cout << "\n"; 548 cout << " POT ENERGY: "; peten->matPrint(); cout << "\n"; 549 cout << " AGGREGATE CM POS: "; cmphasePos->vecPrint(); cout << "\n"; 550 cout << " AGGREGATE CM VEL: "; cmphaseVel->vecPrint(); cout << "\n"; 551 cout << " ANGULAR MOMENTUM: "; amvec->vecPrint(); cout << "\n\n"; 552 } 553 554 // ---------------------------------------------------------------------------- 555 void NBodySystem::output(parms *sim) { 556 if( (tout - 0.01 * sim->getdtime()) <= tnow) 557 tout += sim->getdtout(); 558 } 559 560 // ---------------------------------------------------------------------------- 561 void NBodySystem::clearTiming(){ 562 tracktime = 0; 563 updatetime = 0; 564 partitiontime = 0; 565 treebuildtime = 0; 566 forcecalctime = 0; 567 } 568 569 // ---------------------------------------------------------------------------- 570 void NBodySystem::initOutput(parms *p){ 571 cout << "\n\t\tOutput: " << p->getheadline() << "\n\n"; 572 cout << "nbody dtime eps tol dtout tstop fcells\n"; 573 cout << nbody << " " << p->getdtime() << " " << p->geteps() 574 << " " << p->gettol() << " " << p->getdtout() << " " 575 << p->gettstop() << " " << p->getfcells() << "\n\n"; 576 } 577 578 // ---------------------------------------------------------------------------- 579 void NBodySystem::outputTiming(){ 580 int ticks; 581 double dticks; 582 583 get_ticks_time(&ticks); 584 dticks = ticks; 585 cout << "COMPUTESTART = " << 586 (computestart / dticks) << "\n"; 587 cout << "COMPUTEEND = " << 588 (computeend / dticks) << "\n"; 589 cout << "COMPUTETIME = " << 590 ((computeend - computestart) / dticks) << "\n"; 591 cout << "TRACKTIME = " << 592 (tracktime / dticks) << "\n"; 593 cout << "PARTITIONTIME = " << 594 (partitiontime / dticks) << "\t" 595 << ((float)partitiontime)/tracktime << "\n"; 596 cout << "TREEBUILDTIME = " << 597 (treebuildtime / dticks) << "\t" 598 << ((float)treebuildtime)/tracktime << "\n"; 599 cout << "FORCECALCTIME = " << 600 (forcecalctime / dticks) << "\t" 601 << ((float)forcecalctime)/tracktime << "\n"; 602 cout << "UPDATEPOSTIME = " << 603 (updatetime / dticks) << "\t" 604 << ((float)updatetime)/tracktime << "\n"; 605 cout << "RESTTIME = " << 606 ((tracktime - partitiontime - 607 updatetime - treebuildtime - forcecalctime) 608 / dticks) << "\t" 609 << ((float)(tracktime-partitiontime- updatetime - treebuildtime- 610 forcecalctime))/tracktime<<"\n"; 611 cout.flush(); 612 } 613 614 // ---------------------------------------------------------------------------- 615 void NBodySystem::loadData(char *name){ 616 int i, ndim; 617 body *p; 618 double m, tmp[3]; 619 ifstream instr(name); 620 621 cerr << "Reading input file : " << name << '\n'; 622 cerr.flush(); 623 if (!instr) 624 cerr << "loadData: cannot find file " << name << '\n'; 625 instr >> nbody; 626 if (nbody < 1) 627 cerr << "loadData: nbody = " << nbody << " is absurd\n"; 628 instr >> ndim; 629 if(ndim != NDIM) 630 cerr << "inputdata: NDIM = " << NDIM << "ndim = " << ndim << " is absurd\n"; 631 instr >> tnow; 632 633 for(i = 0; i < nbody; i++){ 634 p = bodies->getbody(i); 635 instr >> m; 636 p->setMass(m); 637 } 638 for(i = 0; i < nbody; i++){ 639 p = bodies->getbody(i); 640 in_vector(instr,tmp); 641 p->PosvecLoad(tmp); 642 } 643 for(i = 0; i < nbody; i++){ 644 p = bodies->getbody(i); 645 in_vector(instr,tmp); 646 p->VelvecLoad(tmp); 647 } 648 instr.close(); 649 } 650 651 // ---------------------------------------------------------------------------- 652 void NBodySystem::loadTestData(){ 653 double rsc, vsc, r, v, x, y, offset, tmp; 654 double cmr[NDIM], cmv[NDIM], tmpv[NDIM], tmpv2[NDIM]; 655 body *p, *cp; 656 int halfnbody, i, k; 657 658 tnow = 0.0; 659 rsc = 9 * PI / 16; 660 vsc = sqrt(1.0 / rsc); 661 662 tvecClr(cmr); 663 tvecClr(cmv); 664 665 srand(123); 666 667 668 669 halfnbody = nbody / 2; 670 if (nbody % 2 != 0) halfnbody++; 671 672 for(k = 0; k < halfnbody; k++) { 673 p = bodies->newbody(); 674 //p->setMass((double)(1.0/nbody)); 675 //p->setMass(0.1 + 0.9*((double)rand()/INT_MAX)); 676 //p->setMass(0.5); 677 678 /* reject radii greater than 10 */ 679 do { 680 tmp = xrand(0.0, 0.999); 681 r = 1.0 / sqrt(pow(tmp, -2.0/3.0) - 1); 682 } while (r > 9.0); 683 684 pickshell(tmpv, rsc * r); 685 p->PosvecAdd(tmpv); 686 tvecAdd(cmr,tmpv); 687 do{ 688 x = xrand(0.0, 1.0); 689 y = xrand(0.0, 0.1); 690 } while (y > x*x * pow(1 - x*x, 3.5)); 691 v = sqrt(2.0) * x / pow(1 + r*r, 0.25); 692 pickshell(tmpv, vsc * v); 693 p->VelvecAdd(tmpv); 694 tvecAdd(cmv,tmpv); 695 } 696 697 offset = 4.0; 698 for(k = halfnbody; k < nbody; k++) { 699 p = bodies->newbody(); 700 //p->setMass((double)(1.0/nbody)); 701 //p->setMass(0.1 + 0.9*((double)rand()/INT_MAX)); 702 //p->setMass(0.5); 703 704 705 cp = bodies->getbody(k-halfnbody); 706 cp->PosvecStore(tmpv); 707 tvecClr(tmpv2); 708 for(i=0; i < NDIM; i++){ 709 tmpv2[i] = tmpv[i] + offset; 710 tvecAdd(cmr,tmpv2); 711 } 712 p->PosvecAdd(tmpv2); 713 cp->VelvecStore(tmpv); 714 tvecClr(tmpv2); 715 for(i=0; i < NDIM; i++){ 716 tmpv2[i] = tmpv[i]; 717 tvecAdd(cmv,tmpv2); 718 } 719 p->VelvecAdd(tmpv2); 720 } 721 722 tvecScale(cmr, 1.0/nbody); 723 tvecScale(cmv, 1.0/nbody); 724 725 for(k=0; k < nbody; k++) { 726 p = bodies->getbody(k); 727 //p->setMass(0.1 + 0.9*((double)rand()/INT_MAX)); 728 p->setMass(0.5); 729 730 p->PosvecStore(tmpv); 731 tvecSub(tmpv,cmr); 732 p->PosvecLoad(tmpv); 733 p->VelvecStore(tmpv); 734 tvecSub(tmpv,cmv); 735 p->VelvecLoad(tmpv); 736 } 737 } 738 739 // ---------------------------------------------------------------------------- 740 void NBodySystem::pickshell(double *vec, double rad){ 741 int k; 742 double rsq, rsc; 743 744 do{ 745 for(k = 0; k < NDIM; k++) 746 vec[k] = xrand(-1.0, 1.0); 747 rsq = tvecDotProd(vec,vec); 748 } while (rsq > 1.0); 749 rsc = rad / sqrt(rsq); 750 tvecScale(vec,rsc); 751 } 752 753 754 // ---------------------------------------------------------------------------- 755 void NBodySystem::startrun(parms *p){ 756 int seed; 757 char *name, headbuf[128]; 758 759 p->loadInfile(); 760 name = p->getinfile(); 761 if(strlen(name) > 1) { 762 sprintf(headbuf, "Hack code: input file %s\n", name); 763 strcpy(p->getheadline(),headbuf); 764 } else { 765 nbody = p->getiparam("nbody"); 766 if(nbody < 1) 767 cerr << "startrun: absurd nbody\n"; 768 seed = p->getiparam("seed"); 769 strcpy(p->getheadline(),"Hack code: Plummer model"); 770 } 771 772 p->loadParms(); 773 774 // allocation of data structures filing is done later 775 cells = new cellset((int)(p->getfcells()*p->getfleaves()*nbody)); 776 leaves = new leafset((int)(p->getfleaves()*nbody)); 777 bodies = new bodyset(nbody); 778 779 tout = tnow + p->getdtout(); 780 pranset(seed); 781 782 name = p->getinfile(); 783 if(strlen(name) > 1) { 784 loadData(name); 785 } else { 786 loadTestData(); 787 } 788 789 } 790 791 // ---------------------------------------------------------------------------- 792 void NBodySystem::stepsystem(parms *p){ 793 int i; 794 body *bptr; 795 double incr, hts, ts, stop, tol, eps; 796 797 int trackstart, trackend; 798 int treebuildstart, treebuildend; 799 int forcecalcstart, forcecalcend; 800 int updatestart, updateend; 801 int totalstep = 0; double now0 = tnow; 802 803 incr = p->getdtime(); 804 hts = p->getdthf(); 805 ts = p->getdtime(); 806 stop = p->gettstop(); 807 tol = p->gettolsq(); 808 eps = p->getepssq(); 809 810 for (int isteps = 0; isteps < (int) ((stop-now0)/incr + 0.1) ; isteps++) {//while (tnow < stop + 0.1 * incr){ 811 get_time(&trackstart); 812 leaves->LeafSetReset(); 813 cells->CellSetReset(); 814 BH_tree = cells->makecell(NULL); 815 816 get_time(&treebuildstart); 817 818 maketree(); 819 ((node*)BH_tree)->computecofm(); 820 get_time(&treebuildend); 821 treebuildtime += treebuildend - treebuildstart; 822 823 get_time(&forcecalcstart); 824 computeForces(tol,hts,eps); 825 get_time(&forcecalcend); 826 forcecalctime += forcecalcend - forcecalcstart; 827 828 /* advance bodies */ 829 get_time(&updatestart); 830 Advance(hts,ts); 831 /* 832 for(i = 0; i < nbody; i++) { 833 bptr = bodies->getbody(i); 834 bptr->advance(hts, ts); 835 } 836 */ 837 838 get_time(&updateend); 839 updatetime += updateend - updatestart; 840 get_time(&trackend); 841 tracktime += trackend - trackstart; 842 843 /* compute bounding box */ 844 setbound(); 845 846 nstep++; 847 tnow = tnow + incr; 848 totalstep++; 849 850 updateEnergy(); 851 outputEnergy(); 852 853 if (totalstep == 2) { 854 reset_stat(); 855 clearTiming(); 856 clearStats(); 857 get_time(&computestart); 858 } 859 } 860 get_time(&computeend); 861 print_stat(); 862 computeend = computeend - computestart; 863 computestart = 0; 864 865 } 866 867 868 // --- PARMS CODE ------------------------------------------------------------ 869 parms::parms() { 870 int i; 871 // SET DEFAULT PARAMETER VALUES 872 strcpy(headline,""); 873 strcpy(infile,""); 874 strcpy(outfile,""); 875 dtime = 0.0; 876 dtout = 0.0; 877 tstop = 0.0; 878 fcells = 0.0; 879 fleaves = 0.0; 880 tol = 0.0; 881 tolsq = 0.0; 882 eps = 0.0; 883 epssq = 0.0; 884 dthf = 0.0; 885 886 for(i = 0; i < 12; i++) 887 defaults[i] = (char*)malloc(MAX_NAME); 888 889 /* file names for input/output */ 890 strcpy(defaults[0],"in="); /* snapshot of initial conditions */ 891 strcpy(defaults[1],"out="); /* stream of output snapshots */ 892 /* params, used if no input specified, to make a Plummer Model */ 893 strcpy(defaults[2],"nbody=1024"); /* number of particles to generate */ 894 strcpy(defaults[3],"seed=123"); /* random number generator seed */ 895 /* params to control N-body integration */ 896 strcpy(defaults[4],"dtime=0.025"); /* integration time-step */ 897 strcpy(defaults[5],"eps=0.05"); /* usual potential softening */ 898 strcpy(defaults[6],"tol=1.0"); /* cell subdivision tolerence */ 899 strcpy(defaults[7],"fcells=2.0"); /* cell allocation parameter */ 900 strcpy(defaults[8],"fleaves=0.5"); /* leaf allocation parameter */ 901 strcpy(defaults[9],"tstop=0.075"); /* time to stop integration */ 902 strcpy(defaults[10],"dtout=0.25"); /* data-output interval */ 903 strcpy(defaults[11],""); 904 } 905 906 void parms::setparms(char **fnames, double *args){ 907 strcpy(headline,fnames[0]); 908 strcpy(infile,fnames[1]); 909 strcpy(outfile,fnames[2]); 910 dtime = args[0]; 911 dtout = args[1]; 912 tstop = args[2]; 913 fcells = args[3]; 914 fleaves = args[4]; 915 tol = args[5]; 916 tolsq = args[6]; 917 eps = args[7]; 918 epssq = args[8]; 919 dthf = args[9]; 920 } 921 922 923 int parms::getiparam(char *name){ 924 getparam(name); 925 return (atoi(buf)); 926 } 927 928 // --------------------------------------------------------------------------- 929 double parms::getdparam(char *name){ 930 getparam(name); 931 return (atof(buf)); 932 } 933 934 void parms::getparam(char *name) { 935 int i, leng; 936 char *temp; 937 938 buf[0] = '\0'; 939 if(defaults == NULL){ 940 cerr << "getparam: called before initparam\n"; 941 } 942 i = scanbind(defaults, name); 943 if(i < 0){ 944 cerr << "getparam: " << name << " unknown\n"; 945 } 946 temp = extrvalue(defaults[i]); 947 fgets(buf, 128, stdin); 948 leng = strlen(buf) + 1; 949 if(leng <= 1) 950 strcpy(buf,temp); 951 } 952 953 // --------------------------------------------------------------------------- 954 void parms::loadInfile(){ 955 getparam("in"); 956 strcpy(infile,buf); 957 } 958 959 // --------------------------------------------------------------------------- 960 void parms::loadParms(){ 961 getparam("out"); 962 strcpy(outfile,buf); 963 dtime = getdparam("dtime"); 964 dthf = 0.5 * dtime; 965 eps = getdparam("eps"); 966 epssq = eps*eps; 967 tol = getdparam("tol"); 968 tolsq = tol*tol; 969 fcells = getdparam("fcells"); 970 fleaves = getdparam("fleaves"); 971 tstop = getdparam("tstop"); 972 dtout = getdparam("dtout"); 973 } 974 975 976 977 // --- CELLSET CODE ------------------------------------------------------------ 978 cell *cellset::newcell(){ 979 cell *c; 980 if(numcell == maxcell){ 981 cerr << "newcell: More than " << maxcell << " cells; increase fcells\n"; 982 return NULL; 983 } 984 c = &celltab[numcell]; 985 numcell++; 986 c->CellClear(); 987 return c; 988 } 989 // --------------------------------------------------------------------------- 990 cell *cellset::makecell(node *parent) { 991 cell *c; 992 c = newcell(); 993 c->setParent(parent); 994 if(c == NULL) return NULL; 995 c->adjLevel(parent); 996 return c; 997 } 998 999 1000 // --- LEAFSET CODE ------------------------------------------------------------ 1001 leaf *leafset::newleaf() { 1002 leaf* p; 1003 if(numleaf == maxleaf){ 1004 cerr << "newleaf: More than " << maxleaf << " leaves; increase fleaves\n"; 1005 return NULL; 1006 } 1007 p = &leaftab[numleaf]; 1008 numleaf++; 1009 p->LeafClear(); 1010 return p; 1011 } 1012 1013 // --------------------------------------------------------------------------- 1014 leaf *leafset::makeleaf(node *parent){ 1015 leaf* l; 1016 l = newleaf(); 1017 if(l == NULL) return NULL; 1018 l->setParent(parent); 1019 l->adjLevel(parent); 1020 return (l); 1021 } 1022 1023 1024 // --- BODYSET CODE ------------------------------------------------------------ 1025 body* bodyset::newbody(){ 1026 body *p; 1027 if(numbody == maxbody){ 1028 cerr << "newbody: More than " << maxbody << "body; increase nbody\n"; 1029 return NULL; 1030 } 1031 p = &bodytab[numbody]; 1032 numbody++; 1033 return p; 1034 } 1035 1036 // --- LEAF CODE ------------------------------------------------------------ 1037 void leaf::LeafAddBody(body *bp){ // Assumes enough capacity 1038 bp->setParent(this); 1039 bp->setLevel(Level()); 1040 bp->setChildNum(num_bodies); 1041 bodyp[num_bodies] = bp; 1042 num_bodies++; 1043 } 1044 1045 void leaf::Subdivide(cell *parent, int idx){ 1046 cell *c; 1047 leaf *le; 1048 double tmpv[NDIM]; 1049 int i, index, xp[NDIM], num, Lev; 1050 body *loctab[LEAFMAXBODIES], *p; 1051 1052 num = num_bodies; 1053 for(i=0; i < num; i++) { 1054 loctab[i] = bodyp[i]; 1055 bodyp[i] = NULL; 1056 } 1057 num_bodies = 0; 1058 c = Nbody.getcells()->makecell(parent); 1059 c->setChildNum(child_num); 1060 parent->setSubp(idx,c); 1061 Lev = level; 1062 1063 p = loctab[0]; 1064 p->PosvecStore(tmpv); 1065 Nbody.intcoord(xp,tmpv); 1066 index = subindex(xp, Lev); 1067 c->setSubp(index,this); 1068 child_num = index; 1069 parent = c; 1070 level = (Lev >> 1); 1071 LeafAddBody(p); 1072 1073 for(i=1; i < num; i++) { 1074 p = loctab[i]; 1075 p->PosvecStore(tmpv); 1076 Nbody.intcoord(xp,tmpv); 1077 index = subindex(xp, Lev); 1078 le = (leaf*)(c->Subp(index)); 1079 if(le == NULL) { 1080 le = Nbody.getleaves()->makeleaf(c); 1081 le->setChildNum(index); 1082 c->setSubp(index,le); 1083 } 1084 le->LeafAddBody(p); 1085 } 1086 } 1087 1088 1089 // --- CELL CODE ------------------------------------------------------------ 1090 void cell::CellAddBody(body *p, int *coords) { 1091 int Lev, kidIndex; 1092 node *ptr; 1093 leaf *le; 1094 1095 kidIndex = subindex(coords,level); 1096 ptr = subp[kidIndex]; 1097 Lev = (level >> 1); 1098 if(Lev != 0){ 1099 if(ptr == NULL) { 1100 le = Nbody.getleaves()->makeleaf(this); 1101 le->setChildNum(kidIndex); 1102 subp[kidIndex] = le; 1103 le->LeafAddBody(p); 1104 } else { 1105 if(ptr->Type() == LEAF){ 1106 le = (leaf*)ptr; 1107 if(le->full()){ 1108 le->Subdivide(this,kidIndex); 1109 CellAddBody(p,coords); 1110 } else { 1111 le->LeafAddBody(p); 1112 } 1113 } else { 1114 ((cell*)ptr)->CellAddBody(p,coords); 1115 } 1116 } 1117 } else { 1118 cerr << " *** Error: Not enough levels in tree...(CellAddBody)\n"; 1119 } 1120 } 1121 1122 // ---------------------------------------------------------------------------- 1123 int node::subindex(int *x, int l){ 1124 int i, k, yes; 1125 1126 i = 0; 1127 yes = 0; 1128 for(k = 0; k < NDIM; k++) { 1129 if(((x[k] & l) && !yes) || (!(x[k] & l) && yes)) { 1130 i += NSUB >> (k + 1); 1131 yes = 1; 1132 } 1133 else yes = 0; 1134 } 1135 return (i); 1136 } 1137 1138 void NBodySystem::setbound(){ 1139 int i; 1140 double tmp1[NDIM], tmp2[NDIM], side; 1141 body *p; 1142 1143 mincorner->vecSetUnit(); 1144 mincorner->vecScale(1E99); 1145 maxcorner->vecSetUnit(); 1146 maxcorner->vecScale(-1E99); 1147 for(i=0; i < nbody; i++) { 1148 p = bodies->getbody(i); 1149 p->PosvecStore(tmp1); 1150 mincorner->vecMin(tmp1); 1151 maxcorner->vecMax(tmp1); 1152 } 1153 1154 side = 0.0; 1155 maxcorner->vecStore(tmp1); 1156 mincorner->vecStore(tmp2); 1157 tvecSub(tmp1,tmp2); 1158 for(i=0; i < NDIM; i++) 1159 if(side < tmp1[i]) 1160 side = tmp1[i]; 1161 1162 for(i=0; i < NDIM; i++) 1163 tmp1[i] = -side/100000; 1164 mincorner->vecAdd(tmp1); 1165 size = 1.00002*side; 1166 } 1167 1168 1169 // ---------------------------------------------------------------------------- 1170 // MAKETREE: initialize tree structure for force calculation. 1171 // ---------------------------------------------------------------------------- 1172 void NBodySystem::maketree(){ 1173 int i, xp[NDIM]; 1174 double v[NDIM]; 1175 body *p; 1176 1177 for(i=0; i < nbody; i++){ 1178 p = bodies->getbody(i); 1179 p->PosvecStore(v); 1180 intcoord(xp,v); 1181 BH_tree->CellAddBody(p,xp); 1182 } 1183 } 1184 1185 1186 void node::computecofm(){ // <--- Parallelizable... 1187 int i; 1188 node *nptr; 1189 body *p; 1190 leaf *le; 1191 cell *ce; 1192 double m; 1193 double tmpv[NDIM]; 1194 1195 if(type == LEAF){ 1196 le = (leaf*)this; 1197 for(i=0; i < le->NumBody(); i++) { 1198 p = le->Bodyp(i); 1199 m = p->Mass(); 1200 mass += m; 1201 p->PosvecStore(tmpv); 1202 tvecScale(tmpv,m); 1203 pos.vecAdd(tmpv); 1204 } 1205 pos.vecStore(tmpv); 1206 tvecScale(tmpv,1.0/mass); 1207 pos.vecLoad(tmpv); 1208 } else { 1209 ce = (cell*)this; 1210 for(i=0; i < NSUB; i++){ 1211 nptr = ce->Subp(i); 1212 if(nptr != NULL){ 1213 nptr->computecofm(); 1214 m = nptr->Mass(); 1215 mass += m; 1216 nptr->PosvecStore(tmpv); 1217 tvecScale(tmpv,m); 1218 pos.vecAdd(tmpv); 1219 } 1220 } 1221 pos.vecStore(tmpv); 1222 tvecScale(tmpv,1.0/mass); 1223 pos.vecLoad(tmpv); 1224 } 1225 } 1226 1227 double body::subdiv(node *p, double *res) { 1228 int i; 1229 double drsq, d; 1230 1231 drsq = 0.0; 1232 for (i = 0; i < NDIM; i++) { 1233 d = p->pos.vecField(i)-pos.vecField(i); 1234 drsq += d*d; 1235 res[i] = d; 1236 } 1237 return(drsq); 1238 } 1239 1240 inline double body::computeInter(node *p, double epsSq, double *res){ 1241 int i; 1242 double drabs, inc, mor3; 1243 double drsq, d; 1244 1245 drsq = epsSq; 1246 for (i = 0; i < NDIM; i++) { 1247 d = p->pos.vecField(i)-pos.vecField(i); 1248 drsq += d*d; 1249 res[i] = d; 1250 } 1251 drabs = sqrt(drsq); 1252 inc = p->Mass() / drabs; 1253 mor3 = inc / drsq; 1254 tvecScale(res,mor3); 1255 return inc; 1256 } 1257 1258 void body::updatePhi(double inc) { 1259 phi -= inc; 1260 } 1261 1262 void body::fastgravsub(node *p, double *ai, double pdrsq, double epssq) { 1263 double drsq; 1264 double drabs, inc, mor3; 1265 double tmpv[NDIM]; 1266 1267 drsq = pdrsq + epssq; 1268 drabs = sqrt(drsq); 1269 inc = p->Mass() / drabs; 1270 // phi -= inc; 1271 updatePhi(inc); 1272 mor3 = inc / drsq; 1273 tvecLoad(tmpv,ai); 1274 tvecScale(tmpv,mor3); 1275 acc.vecAdd(tmpv); 1276 } 1277 1278 void body::gravsub(node *p, double epsSq){ 1279 double phii; 1280 double ai[NDIM]; 1281 1282 phii = computeInter(p,epsSq,ai); 1283 // phi -= phii; 1284 updatePhi(phii); 1285 acc.vecAdd(ai); 1286 } 1287 1288 inline void body::openCell(node *n, double tolsq, double dsq, double epssq){ 1289 int i; 1290 node* nn; 1291 1292 for(i=0; i < NSUB; i++) { 1293 nn = ((cell*)n)->Subp(i); 1294 if(nn != NULL){ 1295 walksub(nn,tolsq,(dsq/4.0),epssq); 1296 } 1297 } 1298 } 1299 1300 inline void body::openLeaf(node *n,double epssq, int iter){ 1301 int i; 1302 body *bptr; 1303 1304 for(i=0; i < iter; i++){ 1305 bptr = ((leaf*)n)->Bodyp(i); 1306 if(bptr != this){ 1307 gravsub((node*)bptr, epssq); 1308 } 1309 } 1310 } 1311 1312 void body::walksub(node *n, double tolsq, double dsq, double epssq){ 1313 double drsq; 1314 double ai[NDIM]; 1315 drsq = subdiv(n,ai); 1316 if((tolsq * drsq) < dsq){ 1317 if((n->Type()) == CELL){ 1318 openCell(n,tolsq,dsq,epssq); 1319 } else { 1320 openLeaf(n,epssq,((leaf*)n)->NumBody()); 1321 } 1322 } else { 1323 fastgravsub(n, ai, drsq, epssq); 1324 } 1325 } 1326 1327 void body::advance(double hts, double ts){ 1328 double dvel[NDIM], vel1[NDIM], dpos[NDIM]; 1329 1330 acc.vecStore(dvel); 1331 vel.vecStore(vel1); 1332 tvecScale(dvel,hts); 1333 tvecAdd(vel1,dvel); 1334 tvecLoad(dpos,vel1); 1335 tvecScale(dpos,ts); 1336 tvecScale(dvel,2.0); 1337 1338 vel.vecAdd(dvel); 1339 pos.vecAdd(dpos); 1340 } 1341 1342 void NBodySystem::Advance(double hts, double ts){ 1343 int i; 1344 body *p; 1345 1346 for(i=0; i < nbody; i++) { 1347 p = bodies->getbody(i); 1348 p->advance(hts,ts); 1349 } 1350 } 1351 1352 void body::swapAcc() { 1353 double tmpv[NDIM]; 1354 1355 phi = 0.0; 1356 acc.vecStore(tmpv); 1357 tmp.vecLoad(tmpv); 1358 acc.vecClr(); 1359 } 1360 1361 void body::startVel(double f){ 1362 double tmpv1[NDIM], tmpv2[NDIM]; 1363 acc.vecStore(tmpv1); 1364 tmp.vecStore(tmpv2); 1365 tvecSub(tmpv1,tmpv2); 1366 tvecScale(tmpv1,f); 1367 vel.vecAdd(tmpv1); 1368 } 1369 1370 void NBodySystem::SwapAccs(){ 1371 int i; 1372 body *p; 1373 1374 for(i=0; i < nbody; i++) { 1375 p = bodies->getbody(i); 1376 p->swapAcc(); 1377 } 1378 } 1379 1380 void NBodySystem::ComputeAccels(double tol, double eps){ 1381 int i; 1382 body *p; 1383 int nbody_local = nbody; 1384 1385 #pragma omp parallel for private (p) schedule(static) 1386 for(i=0; i < nbody_local; i++) { 1387 p = bodies->getbody(i); 1388 p->walksub(BH_tree, tol, size*size, eps); 1389 } 1390 } 1391 1392 void NBodySystem::StartVels(double hts){ 1393 int i; 1394 body *p; 1395 1396 for(i=0; i < nbody; i++) { 1397 p = bodies->getbody(i); 1398 p->startVel(hts); 1399 } 1400 } 1401 1402 void NBodySystem::computeForces(double tol, double hts, double eps){ 1403 SwapAccs(); 1404 ComputeAccels(tol,eps); 1405 if(nstep > 0){ 1406 StartVels(hts); 1407 } 1408 } 1409 1410 // ---------------------------------------------------------------------------- 1411 int main(int argc, char **argv){ 1412 Nbody.init(); 1413 Nbody.startrun(&simparms); 1414 Nbody.initOutput(&simparms); 1415 Nbody.setbound(); 1416 Nbody.clearTiming(); 1417 Nbody.clearStats(); 1418 1419 Nbody.stepsystem(&simparms); 1420 1421 // bodyset * bset = Nbody.getbodies(); 1422 //for (int i = 0; i < bset->size(); i++) 1423 // bset-> 1424 1425 Nbody.outputTiming(); 1426 //Nbody.outputStats(); 1427 Nbody.output(&simparms); 1428 Nbody.outputEnergy(); 1429 //Nbody.dump(); 1430 } 1431