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 for(k = 0; k < nbody; k++) { 501 p = bodies->getbody(k); 502 matrix m1; 503 double v0[3], v1[3], v2[3], v3[3], v4[3]; 504 m = p->Mass(); 505 mtot.addval(m); 506 p->VelvecStore(v1); 507 velsq = tvecDotProd(v1,v1); 508 v4[1] = 0.5 * m * velsq; 509 v4[2] = 0.5 * m * p->Phi(); 510 v4[0] = v4[1] + v4[2]; 511 etot->vecAdd(v4); 512 p->PosvecStore(v1); 513 p->VelvecStore(v2); 514 p->AccvecStore(v3); 515 tvecLoad(v0,v2); 516 tvecScale(v0,0.5*m); 517 m1.matOutProd(v0,v2); 518 keten->matAdd(&m1); 519 tvecLoad(v0,v1); 520 tvecScale(v0,m); 521 m1.matOutProd(v0,v3); 522 peten->matAdd(&m1); 523 tvecLoad(v0,v1); 524 tvecScale(v0,m); 525 cmphasePos->vecAdd(v0); 526 tvecLoad(v0,v2); 527 tvecScale(v0,m); 528 cmphaseVel->vecAdd(v0); 529 tvecCrossProd(v0,v1,v2); 530 tvecScale(v0,m); 531 amvec->vecAdd(v0); 532 } 533 534 if (mtot.readval() != 0.0){ 535 cmphasePos->vecScale(1.0/mtot.readval()); 536 cmphaseVel->vecScale(1.0/mtot.readval()); 537 } 538 } 539 540 541 // ---------------------------------------------------------------------------- 542 void NBodySystem::outputEnergy(){ 543 cout << "\n\nENERGY TOTALS\n"; 544 cout << " MASS: "<< mtot.readval() << "\n"; 545 cout << " ENERGY: "; etot->vecPrint(); cout << "\n"; 546 cout << " KIN ENERGY: "; keten->matPrint(); cout << "\n"; 547 cout << " POT ENERGY: "; peten->matPrint(); cout << "\n"; 548 cout << " AGGREGATE CM POS: "; cmphasePos->vecPrint(); cout << "\n"; 549 cout << " AGGREGATE CM VEL: "; cmphaseVel->vecPrint(); cout << "\n"; 550 cout << " ANGULAR MOMENTUM: "; amvec->vecPrint(); cout << "\n\n"; 551 } 552 553 // ---------------------------------------------------------------------------- 554 void NBodySystem::output(parms *sim) { 555 if( (tout - 0.01 * sim->getdtime()) <= tnow) 556 tout += sim->getdtout(); 557 } 558 559 // ---------------------------------------------------------------------------- 560 void NBodySystem::clearTiming(){ 561 tracktime = 0; 562 updatetime = 0; 563 partitiontime = 0; 564 treebuildtime = 0; 565 forcecalctime = 0; 566 } 567 568 // ---------------------------------------------------------------------------- 569 void NBodySystem::initOutput(parms *p){ 570 cout << "\n\t\tOutput: " << p->getheadline() << "\n\n"; 571 cout << "nbody dtime eps tol dtout tstop fcells\n"; 572 cout << nbody << " " << p->getdtime() << " " << p->geteps() 573 << " " << p->gettol() << " " << p->getdtout() << " " 574 << p->gettstop() << " " << p->getfcells() << "\n\n"; 575 } 576 577 // ---------------------------------------------------------------------------- 578 void NBodySystem::outputTiming(){ 579 int ticks; 580 double dticks; 581 582 get_ticks_time(&ticks); 583 dticks = ticks; 584 cout << "COMPUTESTART = " << 585 (computestart / dticks) << "\n"; 586 cout << "COMPUTEEND = " << 587 (computeend / dticks) << "\n"; 588 cout << "COMPUTETIME = " << 589 ((computeend - computestart) / dticks) << "\n"; 590 cout << "TRACKTIME = " << 591 (tracktime / dticks) << "\n"; 592 cout << "PARTITIONTIME = " << 593 (partitiontime / dticks) << "\t" 594 << ((float)partitiontime)/tracktime << "\n"; 595 cout << "TREEBUILDTIME = " << 596 (treebuildtime / dticks) << "\t" 597 << ((float)treebuildtime)/tracktime << "\n"; 598 cout << "FORCECALCTIME = " << 599 (forcecalctime / dticks) << "\t" 600 << ((float)forcecalctime)/tracktime << "\n"; 601 cout << "UPDATEPOSTIME = " << 602 (updatetime / dticks) << "\t" 603 << ((float)updatetime)/tracktime << "\n"; 604 cout << "RESTTIME = " << 605 ((tracktime - partitiontime - 606 updatetime - treebuildtime - forcecalctime) 607 / dticks) << "\t" 608 << ((float)(tracktime-partitiontime- updatetime - treebuildtime- 609 forcecalctime))/tracktime<<"\n"; 610 cout.flush(); 611 } 612 613 // ---------------------------------------------------------------------------- 614 void NBodySystem::loadData(char *name){ 615 int i, ndim; 616 body *p; 617 double m, tmp[3]; 618 ifstream instr(name); 619 620 cerr << "Reading input file : " << name << '\n'; 621 cerr.flush(); 622 if (!instr) 623 cerr << "loadData: cannot find file " << name << '\n'; 624 instr >> nbody; 625 if (nbody < 1) 626 cerr << "loadData: nbody = " << nbody << " is absurd\n"; 627 instr >> ndim; 628 if(ndim != NDIM) 629 cerr << "inputdata: NDIM = " << NDIM << "ndim = " << ndim << " is absurd\n"; 630 instr >> tnow; 631 632 for(i = 0; i < nbody; i++){ 633 p = bodies->getbody(i); 634 instr >> m; 635 p->setMass(m); 636 } 637 for(i = 0; i < nbody; i++){ 638 p = bodies->getbody(i); 639 in_vector(instr,tmp); 640 p->PosvecLoad(tmp); 641 } 642 for(i = 0; i < nbody; i++){ 643 p = bodies->getbody(i); 644 in_vector(instr,tmp); 645 p->VelvecLoad(tmp); 646 } 647 instr.close(); 648 } 649 650 // ---------------------------------------------------------------------------- 651 void NBodySystem::loadTestData(){ 652 double rsc, vsc, r, v, x, y, offset, tmp; 653 double cmr[NDIM], cmv[NDIM], tmpv[NDIM], tmpv2[NDIM]; 654 body *p, *cp; 655 int halfnbody, i, k; 656 657 tnow = 0.0; 658 rsc = 9 * PI / 16; 659 vsc = sqrt(1.0 / rsc); 660 661 tvecClr(cmr); 662 tvecClr(cmv); 663 664 srand(123); 665 666 667 668 halfnbody = nbody / 2; 669 if (nbody % 2 != 0) halfnbody++; 670 671 for(k = 0; k < halfnbody; k++) { 672 p = bodies->newbody(); 673 //p->setMass((double)(1.0/nbody)); 674 //p->setMass(0.1 + 0.9*((double)rand()/INT_MAX)); 675 //p->setMass(0.5); 676 677 /* reject radii greater than 10 */ 678 do { 679 tmp = xrand(0.0, 0.999); 680 r = 1.0 / sqrt(pow(tmp, -2.0/3.0) - 1); 681 } while (r > 9.0); 682 683 pickshell(tmpv, rsc * r); 684 p->PosvecAdd(tmpv); 685 tvecAdd(cmr,tmpv); 686 do{ 687 x = xrand(0.0, 1.0); 688 y = xrand(0.0, 0.1); 689 } while (y > x*x * pow(1 - x*x, 3.5)); 690 v = sqrt(2.0) * x / pow(1 + r*r, 0.25); 691 pickshell(tmpv, vsc * v); 692 p->VelvecAdd(tmpv); 693 tvecAdd(cmv,tmpv); 694 } 695 696 offset = 4.0; 697 for(k = halfnbody; k < nbody; k++) { 698 p = bodies->newbody(); 699 //p->setMass((double)(1.0/nbody)); 700 //p->setMass(0.1 + 0.9*((double)rand()/INT_MAX)); 701 //p->setMass(0.5); 702 703 704 cp = bodies->getbody(k-halfnbody); 705 cp->PosvecStore(tmpv); 706 tvecClr(tmpv2); 707 for(i=0; i < NDIM; i++){ 708 tmpv2[i] = tmpv[i] + offset; 709 tvecAdd(cmr,tmpv2); 710 } 711 p->PosvecAdd(tmpv2); 712 cp->VelvecStore(tmpv); 713 tvecClr(tmpv2); 714 for(i=0; i < NDIM; i++){ 715 tmpv2[i] = tmpv[i]; 716 tvecAdd(cmv,tmpv2); 717 } 718 p->VelvecAdd(tmpv2); 719 } 720 721 tvecScale(cmr, 1.0/nbody); 722 tvecScale(cmv, 1.0/nbody); 723 724 for(k=0; k < nbody; k++) { 725 p = bodies->getbody(k); 726 //p->setMass(0.1 + 0.9*((double)rand()/INT_MAX)); 727 p->setMass(0.5); 728 729 p->PosvecStore(tmpv); 730 tvecSub(tmpv,cmr); 731 p->PosvecLoad(tmpv); 732 p->VelvecStore(tmpv); 733 tvecSub(tmpv,cmv); 734 p->VelvecLoad(tmpv); 735 } 736 } 737 738 // ---------------------------------------------------------------------------- 739 void NBodySystem::pickshell(double *vec, double rad){ 740 int k; 741 double rsq, rsc; 742 743 do{ 744 for(k = 0; k < NDIM; k++) 745 vec[k] = xrand(-1.0, 1.0); 746 rsq = tvecDotProd(vec,vec); 747 } while (rsq > 1.0); 748 rsc = rad / sqrt(rsq); 749 tvecScale(vec,rsc); 750 } 751 752 753 // ---------------------------------------------------------------------------- 754 void NBodySystem::startrun(parms *p){ 755 int seed; 756 char *name, headbuf[128]; 757 758 p->loadInfile(); 759 name = p->getinfile(); 760 if(strlen(name) > 1) { 761 sprintf(headbuf, "Hack code: input file %s\n", name); 762 strcpy(p->getheadline(),headbuf); 763 } else { 764 nbody = p->getiparam("nbody"); 765 if(nbody < 1) 766 cerr << "startrun: absurd nbody\n"; 767 seed = p->getiparam("seed"); 768 strcpy(p->getheadline(),"Hack code: Plummer model"); 769 } 770 771 p->loadParms(); 772 773 // allocation of data structures filing is done later 774 cells = new cellset((int)(p->getfcells()*p->getfleaves()*nbody)); 775 leaves = new leafset((int)(p->getfleaves()*nbody)); 776 bodies = new bodyset(nbody); 777 778 tout = tnow + p->getdtout(); 779 pranset(seed); 780 781 name = p->getinfile(); 782 if(strlen(name) > 1) { 783 loadData(name); 784 } else { 785 loadTestData(); 786 } 787 788 } 789 790 // ---------------------------------------------------------------------------- 791 void NBodySystem::stepsystem(parms *p){ 792 int i; 793 body *bptr; 794 double incr, hts, ts, stop, tol, eps; 795 796 int trackstart, trackend; 797 int treebuildstart, treebuildend; 798 int forcecalcstart, forcecalcend; 799 int updatestart, updateend; 800 int totalstep = 0; double now0 = tnow; 801 802 incr = p->getdtime(); 803 hts = p->getdthf(); 804 ts = p->getdtime(); 805 stop = p->gettstop(); 806 tol = p->gettolsq(); 807 eps = p->getepssq(); 808 809 for (int isteps = 0; isteps < (int) ((stop-now0)/incr + 0.1) ; isteps++) {//while (tnow < stop + 0.1 * incr){ 810 get_time(&trackstart); 811 leaves->LeafSetReset(); 812 cells->CellSetReset(); 813 BH_tree = cells->makecell(NULL); 814 815 get_time(&treebuildstart); 816 817 maketree(); 818 ((node*)BH_tree)->computecofm(); 819 get_time(&treebuildend); 820 treebuildtime += treebuildend - treebuildstart; 821 822 get_time(&forcecalcstart); 823 computeForces(tol,hts,eps); 824 get_time(&forcecalcend); 825 forcecalctime += forcecalcend - forcecalcstart; 826 827 /* advance bodies */ 828 get_time(&updatestart); 829 Advance(hts,ts); 830 /* 831 for(i = 0; i < nbody; i++) { 832 bptr = bodies->getbody(i); 833 bptr->advance(hts, ts); 834 } 835 */ 836 837 get_time(&updateend); 838 updatetime += updateend - updatestart; 839 get_time(&trackend); 840 tracktime += trackend - trackstart; 841 842 /* compute bounding box */ 843 setbound(); 844 845 nstep++; 846 tnow = tnow + incr; 847 totalstep++; 848 849 updateEnergy(); 850 outputEnergy(); 851 852 if (totalstep == 2) { 853 reset_stat(); 854 clearTiming(); 855 clearStats(); 856 get_time(&computestart); 857 } 858 } 859 get_time(&computeend); 860 print_stat(); 861 computeend = computeend - computestart; 862 computestart = 0; 863 864 } 865 866 867 // --- PARMS CODE ------------------------------------------------------------ 868 parms::parms() { 869 int i; 870 // SET DEFAULT PARAMETER VALUES 871 strcpy(headline,""); 872 strcpy(infile,""); 873 strcpy(outfile,""); 874 dtime = 0.0; 875 dtout = 0.0; 876 tstop = 0.0; 877 fcells = 0.0; 878 fleaves = 0.0; 879 tol = 0.0; 880 tolsq = 0.0; 881 eps = 0.0; 882 epssq = 0.0; 883 dthf = 0.0; 884 885 for(i = 0; i < 12; i++) 886 defaults[i] = (char*)malloc(MAX_NAME); 887 888 /* file names for input/output */ 889 strcpy(defaults[0],"in="); /* snapshot of initial conditions */ 890 strcpy(defaults[1],"out="); /* stream of output snapshots */ 891 /* params, used if no input specified, to make a Plummer Model */ 892 strcpy(defaults[2],"nbody=1024"); /* number of particles to generate */ 893 strcpy(defaults[3],"seed=123"); /* random number generator seed */ 894 /* params to control N-body integration */ 895 strcpy(defaults[4],"dtime=0.025"); /* integration time-step */ 896 strcpy(defaults[5],"eps=0.05"); /* usual potential softening */ 897 strcpy(defaults[6],"tol=1.0"); /* cell subdivision tolerence */ 898 strcpy(defaults[7],"fcells=2.0"); /* cell allocation parameter */ 899 strcpy(defaults[8],"fleaves=0.5"); /* leaf allocation parameter */ 900 strcpy(defaults[9],"tstop=0.075"); /* time to stop integration */ 901 strcpy(defaults[10],"dtout=0.25"); /* data-output interval */ 902 strcpy(defaults[11],""); 903 } 904 905 void parms::setparms(char **fnames, double *args){ 906 strcpy(headline,fnames[0]); 907 strcpy(infile,fnames[1]); 908 strcpy(outfile,fnames[2]); 909 dtime = args[0]; 910 dtout = args[1]; 911 tstop = args[2]; 912 fcells = args[3]; 913 fleaves = args[4]; 914 tol = args[5]; 915 tolsq = args[6]; 916 eps = args[7]; 917 epssq = args[8]; 918 dthf = args[9]; 919 } 920 921 922 int parms::getiparam(char *name){ 923 getparam(name); 924 return (atoi(buf)); 925 } 926 927 // --------------------------------------------------------------------------- 928 double parms::getdparam(char *name){ 929 getparam(name); 930 return (atof(buf)); 931 } 932 933 void parms::getparam(char *name) { 934 int i, leng; 935 char *temp; 936 937 buf[0] = '\0'; 938 if(defaults == NULL){ 939 cerr << "getparam: called before initparam\n"; 940 } 941 i = scanbind(defaults, name); 942 if(i < 0){ 943 cerr << "getparam: " << name << " unknown\n"; 944 } 945 temp = extrvalue(defaults[i]); 946 fgets(buf, 128, stdin); 947 leng = strlen(buf) + 1; 948 if(leng <= 1) 949 strcpy(buf,temp); 950 } 951 952 // --------------------------------------------------------------------------- 953 void parms::loadInfile(){ 954 getparam("in"); 955 strcpy(infile,buf); 956 } 957 958 // --------------------------------------------------------------------------- 959 void parms::loadParms(){ 960 getparam("out"); 961 strcpy(outfile,buf); 962 dtime = getdparam("dtime"); 963 dthf = 0.5 * dtime; 964 eps = getdparam("eps"); 965 epssq = eps*eps; 966 tol = getdparam("tol"); 967 tolsq = tol*tol; 968 fcells = getdparam("fcells"); 969 fleaves = getdparam("fleaves"); 970 tstop = getdparam("tstop"); 971 dtout = getdparam("dtout"); 972 } 973 974 975 976 // --- CELLSET CODE ------------------------------------------------------------ 977 cell *cellset::newcell(){ 978 cell *c; 979 if(numcell == maxcell){ 980 cerr << "newcell: More than " << maxcell << " cells; increase fcells\n"; 981 return NULL; 982 } 983 c = &celltab[numcell]; 984 numcell++; 985 c->CellClear(); 986 return c; 987 } 988 // --------------------------------------------------------------------------- 989 cell *cellset::makecell(node *parent) { 990 cell *c; 991 c = newcell(); 992 c->setParent(parent); 993 if(c == NULL) return NULL; 994 c->adjLevel(parent); 995 return c; 996 } 997 998 999 // --- LEAFSET CODE ------------------------------------------------------------ 1000 leaf *leafset::newleaf() { 1001 leaf* p; 1002 if(numleaf == maxleaf){ 1003 cerr << "newleaf: More than " << maxleaf << " leaves; increase fleaves\n"; 1004 return NULL; 1005 } 1006 p = &leaftab[numleaf]; 1007 numleaf++; 1008 p->LeafClear(); 1009 return p; 1010 } 1011 1012 // --------------------------------------------------------------------------- 1013 leaf *leafset::makeleaf(node *parent){ 1014 leaf* l; 1015 l = newleaf(); 1016 if(l == NULL) return NULL; 1017 l->setParent(parent); 1018 l->adjLevel(parent); 1019 return (l); 1020 } 1021 1022 1023 // --- BODYSET CODE ------------------------------------------------------------ 1024 body* bodyset::newbody(){ 1025 body *p; 1026 if(numbody == maxbody){ 1027 cerr << "newbody: More than " << maxbody << "body; increase nbody\n"; 1028 return NULL; 1029 } 1030 p = &bodytab[numbody]; 1031 numbody++; 1032 return p; 1033 } 1034 1035 // --- LEAF CODE ------------------------------------------------------------ 1036 void leaf::LeafAddBody(body *bp){ // Assumes enough capacity 1037 bp->setParent(this); 1038 bp->setLevel(Level()); 1039 bp->setChildNum(num_bodies); 1040 bodyp[num_bodies] = bp; 1041 num_bodies++; 1042 } 1043 1044 void leaf::Subdivide(cell *parent, int idx){ 1045 cell *c; 1046 leaf *le; 1047 double tmpv[NDIM]; 1048 int i, index, xp[NDIM], num, Lev; 1049 body *loctab[LEAFMAXBODIES], *p; 1050 1051 num = num_bodies; 1052 for(i=0; i < num; i++) { 1053 loctab[i] = bodyp[i]; 1054 bodyp[i] = NULL; 1055 } 1056 num_bodies = 0; 1057 c = Nbody.getcells()->makecell(parent); 1058 c->setChildNum(child_num); 1059 parent->setSubp(idx,c); 1060 Lev = level; 1061 1062 p = loctab[0]; 1063 p->PosvecStore(tmpv); 1064 Nbody.intcoord(xp,tmpv); 1065 index = subindex(xp, Lev); 1066 c->setSubp(index,this); 1067 child_num = index; 1068 parent = c; 1069 level = (Lev >> 1); 1070 LeafAddBody(p); 1071 1072 for(i=1; i < num; i++) { 1073 p = loctab[i]; 1074 p->PosvecStore(tmpv); 1075 Nbody.intcoord(xp,tmpv); 1076 index = subindex(xp, Lev); 1077 le = (leaf*)(c->Subp(index)); 1078 if(le == NULL) { 1079 le = Nbody.getleaves()->makeleaf(c); 1080 le->setChildNum(index); 1081 c->setSubp(index,le); 1082 } 1083 le->LeafAddBody(p); 1084 } 1085 } 1086 1087 1088 // --- CELL CODE ------------------------------------------------------------ 1089 void cell::CellAddBody(body *p, int *coords) { 1090 int Lev, kidIndex; 1091 node *ptr; 1092 leaf *le; 1093 1094 kidIndex = subindex(coords,level); 1095 ptr = subp[kidIndex]; 1096 Lev = (level >> 1); 1097 if(Lev != 0){ 1098 if(ptr == NULL) { 1099 le = Nbody.getleaves()->makeleaf(this); 1100 le->setChildNum(kidIndex); 1101 subp[kidIndex] = le; 1102 le->LeafAddBody(p); 1103 } else { 1104 if(ptr->Type() == LEAF){ 1105 le = (leaf*)ptr; 1106 if(le->full()){ 1107 le->Subdivide(this,kidIndex); 1108 CellAddBody(p,coords); 1109 } else { 1110 le->LeafAddBody(p); 1111 } 1112 } else { 1113 ((cell*)ptr)->CellAddBody(p,coords); 1114 } 1115 } 1116 } else { 1117 cerr << " *** Error: Not enough levels in tree...(CellAddBody)\n"; 1118 } 1119 } 1120 1121 // ---------------------------------------------------------------------------- 1122 int node::subindex(int *x, int l){ 1123 int i, k, yes; 1124 1125 i = 0; 1126 yes = 0; 1127 for(k = 0; k < NDIM; k++) { 1128 if(((x[k] & l) && !yes) || (!(x[k] & l) && yes)) { 1129 i += NSUB >> (k + 1); 1130 yes = 1; 1131 } 1132 else yes = 0; 1133 } 1134 return (i); 1135 } 1136 1137 void NBodySystem::setbound(){ 1138 int i; 1139 double tmp1[NDIM], tmp2[NDIM], side; 1140 body *p; 1141 1142 mincorner->vecSetUnit(); 1143 mincorner->vecScale(1E99); 1144 maxcorner->vecSetUnit(); 1145 maxcorner->vecScale(-1E99); 1146 for(i=0; i < nbody; i++) { 1147 p = bodies->getbody(i); 1148 p->PosvecStore(tmp1); 1149 mincorner->vecMin(tmp1); 1150 maxcorner->vecMax(tmp1); 1151 } 1152 1153 side = 0.0; 1154 maxcorner->vecStore(tmp1); 1155 mincorner->vecStore(tmp2); 1156 tvecSub(tmp1,tmp2); 1157 for(i=0; i < NDIM; i++) 1158 if(side < tmp1[i]) 1159 side = tmp1[i]; 1160 1161 for(i=0; i < NDIM; i++) 1162 tmp1[i] = -side/100000; 1163 mincorner->vecAdd(tmp1); 1164 size = 1.00002*side; 1165 } 1166 1167 1168 // ---------------------------------------------------------------------------- 1169 // MAKETREE: initialize tree structure for force calculation. 1170 // ---------------------------------------------------------------------------- 1171 void NBodySystem::maketree(){ 1172 int i, xp[NDIM]; 1173 double v[NDIM]; 1174 body *p; 1175 1176 for(i=0; i < nbody; i++){ 1177 p = bodies->getbody(i); 1178 p->PosvecStore(v); 1179 intcoord(xp,v); 1180 BH_tree->CellAddBody(p,xp); 1181 } 1182 } 1183 1184 1185 void node::computecofm(){ // <--- Parallelizable... 1186 int i; 1187 node *nptr; 1188 body *p; 1189 leaf *le; 1190 cell *ce; 1191 double m; 1192 double tmpv[NDIM]; 1193 1194 if(type == LEAF){ 1195 le = (leaf*)this; 1196 for(i=0; i < le->NumBody(); i++) { 1197 p = le->Bodyp(i); 1198 m = p->Mass(); 1199 mass += m; 1200 p->PosvecStore(tmpv); 1201 tvecScale(tmpv,m); 1202 pos.vecAdd(tmpv); 1203 } 1204 pos.vecStore(tmpv); 1205 tvecScale(tmpv,1.0/mass); 1206 pos.vecLoad(tmpv); 1207 } else { 1208 ce = (cell*)this; 1209 for(i=0; i < NSUB; i++){ 1210 nptr = ce->Subp(i); 1211 if(nptr != NULL){ 1212 nptr->computecofm(); 1213 m = nptr->Mass(); 1214 mass += m; 1215 nptr->PosvecStore(tmpv); 1216 tvecScale(tmpv,m); 1217 pos.vecAdd(tmpv); 1218 } 1219 } 1220 pos.vecStore(tmpv); 1221 tvecScale(tmpv,1.0/mass); 1222 pos.vecLoad(tmpv); 1223 } 1224 } 1225 1226 double body::subdiv(node *p, double *res) { 1227 int i; 1228 double drsq, d; 1229 1230 drsq = 0.0; 1231 for (i = 0; i < NDIM; i++) { 1232 d = p->pos.vecField(i)-pos.vecField(i); 1233 drsq += d*d; 1234 res[i] = d; 1235 } 1236 return(drsq); 1237 } 1238 1239 inline double body::computeInter(node *p, double epsSq, double *res){ 1240 int i; 1241 double drabs, inc, mor3; 1242 double drsq, d; 1243 1244 drsq = epsSq; 1245 for (i = 0; i < NDIM; i++) { 1246 d = p->pos.vecField(i)-pos.vecField(i); 1247 drsq += d*d; 1248 res[i] = d; 1249 } 1250 drabs = sqrt(drsq); 1251 inc = p->Mass() / drabs; 1252 mor3 = inc / drsq; 1253 tvecScale(res,mor3); 1254 return inc; 1255 } 1256 1257 void body::updatePhi(double inc) { 1258 phi -= inc; 1259 } 1260 1261 void body::fastgravsub(node *p, double *ai, double pdrsq, double epssq) { 1262 double drsq; 1263 double drabs, inc, mor3; 1264 double tmpv[NDIM]; 1265 1266 drsq = pdrsq + epssq; 1267 drabs = sqrt(drsq); 1268 inc = p->Mass() / drabs; 1269 // phi -= inc; 1270 updatePhi(inc); 1271 mor3 = inc / drsq; 1272 tvecLoad(tmpv,ai); 1273 tvecScale(tmpv,mor3); 1274 acc.vecAdd(tmpv); 1275 } 1276 1277 void body::gravsub(node *p, double epsSq){ 1278 double phii; 1279 double ai[NDIM]; 1280 1281 phii = computeInter(p,epsSq,ai); 1282 // phi -= phii; 1283 updatePhi(phii); 1284 acc.vecAdd(ai); 1285 } 1286 1287 inline void body::openCell(node *n, double tolsq, double dsq, double epssq){ 1288 int i; 1289 node* nn; 1290 1291 for(i=0; i < NSUB; i++) { 1292 nn = ((cell*)n)->Subp(i); 1293 if(nn != NULL){ 1294 walksub(nn,tolsq,(dsq/4.0),epssq); 1295 } 1296 } 1297 } 1298 1299 inline void body::openLeaf(node *n,double epssq, int iter){ 1300 int i; 1301 body *bptr; 1302 1303 for(i=0; i < iter; i++){ 1304 bptr = ((leaf*)n)->Bodyp(i); 1305 if(bptr != this){ 1306 gravsub((node*)bptr, epssq); 1307 } 1308 } 1309 } 1310 1311 void body::walksub(node *n, double tolsq, double dsq, double epssq){ 1312 double drsq; 1313 double ai[NDIM]; 1314 drsq = subdiv(n,ai); 1315 if((tolsq * drsq) < dsq){ 1316 if((n->Type()) == CELL){ 1317 openCell(n,tolsq,dsq,epssq); 1318 } else { 1319 openLeaf(n,epssq,((leaf*)n)->NumBody()); 1320 } 1321 } else { 1322 fastgravsub(n, ai, drsq, epssq); 1323 } 1324 } 1325 1326 void body::advance(double hts, double ts){ 1327 double dvel[NDIM], vel1[NDIM], dpos[NDIM]; 1328 1329 acc.vecStore(dvel); 1330 vel.vecStore(vel1); 1331 tvecScale(dvel,hts); 1332 tvecAdd(vel1,dvel); 1333 tvecLoad(dpos,vel1); 1334 tvecScale(dpos,ts); 1335 tvecScale(dvel,2.0); 1336 1337 vel.vecAdd(dvel); 1338 pos.vecAdd(dpos); 1339 } 1340 1341 void NBodySystem::Advance(double hts, double ts){ 1342 int i; 1343 body *p; 1344 1345 for(i=0; i < nbody; i++) { 1346 p = bodies->getbody(i); 1347 p->advance(hts,ts); 1348 } 1349 } 1350 1351 void body::swapAcc() { 1352 double tmpv[NDIM]; 1353 1354 phi = 0.0; 1355 acc.vecStore(tmpv); 1356 tmp.vecLoad(tmpv); 1357 acc.vecClr(); 1358 } 1359 1360 void body::startVel(double f){ 1361 double tmpv1[NDIM], tmpv2[NDIM]; 1362 acc.vecStore(tmpv1); 1363 tmp.vecStore(tmpv2); 1364 tvecSub(tmpv1,tmpv2); 1365 tvecScale(tmpv1,f); 1366 vel.vecAdd(tmpv1); 1367 } 1368 1369 void NBodySystem::SwapAccs(){ 1370 int i; 1371 body *p; 1372 1373 for(i=0; i < nbody; i++) { 1374 p = bodies->getbody(i); 1375 p->swapAcc(); 1376 } 1377 } 1378 1379 void NBodySystem::ComputeAccels(double tol, double eps){ 1380 int i; 1381 body *p; 1382 int nbody_local = nbody; 1383 1384 #pragma omp parallel for private (p) schedule(dynamic) 1385 for(i=0; i < nbody_local; i++) { 1386 p = bodies->getbody(i); 1387 p->walksub(BH_tree, tol, size*size, eps); 1388 } 1389 } 1390 1391 void NBodySystem::StartVels(double hts){ 1392 int i; 1393 body *p; 1394 1395 for(i=0; i < nbody; i++) { 1396 p = bodies->getbody(i); 1397 p->startVel(hts); 1398 } 1399 } 1400 1401 void NBodySystem::computeForces(double tol, double hts, double eps){ 1402 SwapAccs(); 1403 ComputeAccels(tol,eps); 1404 if(nstep > 0){ 1405 StartVels(hts); 1406 } 1407 } 1408 1409 // ---------------------------------------------------------------------------- 1410 int main(int argc, char **argv){ 1411 Nbody.init(); 1412 Nbody.startrun(&simparms); 1413 Nbody.initOutput(&simparms); 1414 Nbody.setbound(); 1415 Nbody.clearTiming(); 1416 Nbody.clearStats(); 1417 1418 Nbody.stepsystem(&simparms); 1419 1420 // bodyset * bset = Nbody.getbodies(); 1421 //for (int i = 0; i < bset->size(); i++) 1422 // bset-> 1423 1424 Nbody.outputTiming(); 1425 //Nbody.outputStats(); 1426 Nbody.output(&simparms); 1427 Nbody.outputEnergy(); 1428 //Nbody.dump(); 1429 } 1430