bh.rar_visual c

/* * CODE.C: export version of hierarchical N-body code. * Copyright (c) 1991, Joshua E. Barnes, Honolulu, HI. * It's free because it's yours. */ #define global /* nada */ #include <cm/cmmd.h> #ifdef PLAIN #include <fcntl.h> #endif #include "future-cell.h" #include "mem-ref.h" #include "defs.h" #include "code.h" #ifdef PLAIN #define ALLOC(proc,sz) mymalloc(sz) #endif #ifdef VERIFY_AFFINITIES #include "affinity.h" CHECK8(cellptr,(Type(p)==BODY),subp[0],subp[1],subp[2],subp[3],subp[4],subp[5],subp[6],subp[7],tree) CHECK1_ACCUM(bodyptr,next,list) #endif int nbody; char *malloc(); double sqrt(), xrand(), my_rand(), exp, log; real pow(); extern icstruct intcoord(bodyptr p, treeptr t); extern int BhDebug; void computegrav(treeptr t, int nstep); nodeptr maketree(bodyptr btab, int nbody, treeptr t, int nsteps, int proc); void vp(bodyptr q, int nstep); typedef struct { vector cmr; vector cmv; bodyptr list; bodyptr tail; } datapoints; #ifdef OLD_FC typedef struct fc_datapoints { int state; struct fc_datapoints *next; word_t returnPC; word_t fp; datapoints value; } future_cell_datapoints; #else typedef struct fc_datapoints { future_cell_impl impl; datapoints value; } future_cell_datapoints; #endif bodyptr testdata(); datapoints uniform_testdata(int proc, int nbody, int seedfactor); void stepsystem(treeptr t, int nstep); treeptr old_main(); void my_free(nodeptr n); bodyptr ubody_alloc(int p); bodyptr movebodies(bodyptr list, int proc); void freetree(nodeptr n); void freetree1(nodeptr n); int arg1; /* Used to setup runtime system, get arguments-- see old_main */ main(int argc, char **argv) { treeptr t; /* Initialize the runtime system */ CMMD_sync_with_nodes(); #ifndef PLAIN SPMDInit(argc,argv); #else CMMD_fset_io_mode(stdout, CMMD_independent); filestuff(); dealwithargs(argc, argv); __MyNodeId = CMMD_self_address(); __InitRegs(__MyNodeId); if (CMMD_self_address() == 0) { #endif chatting("nbody = %d, numnodes = %d\n", nbody, __NumNodes); CMMD_node_timer_clear(0); CMMD_node_timer_start(0); t = old_main(); CMMD_node_timer_stop(0); chatting("Time = %f seconds\n", CMMD_node_timer_elapsed(0)); #ifdef VERIFY_AFFINITIES Print_Accumulated_list(); #endif #ifndef PLAIN __ShutDown(); #else } #endif } /* global! */ /* Main routine from original program */ treeptr old_main() { real tnow; real tout; int i, nsteps, nprocs; treeptr t; bodyptr bt0,p; long t1, t2; vector cmr, cmv; bodyptr prev=NULL; int tmp=0, range=((1<<NDIM) << NDIM) / __NumNodes; int bodiesper[MAX_NUM_NODES]; bodyptr ptrper[MAX_NUM_NODES]; future_cell_datapoints fc[32]; srand(123); /* set random generator */ /* Tree data structure is global, points to root, and bodylist, has size info */ t = (treeptr) ALLOC(0,sizeof(tree)); Root(t) = NULL; t->rmin[0] = -2.0; t->rmin[1] = -2.0; t->rmin[2] = -2.0; t->rsize = -2.0 * -2.0; /*t->rmin[0];*/ CLRV(cmr); CLRV(cmv); /* Creates a list of bodies */ #ifdef PLAIN for (i=0; i < 32; i++) { datapoints points; int processor= i/(32/__NumNodes); points=uniform_testdata(processor, nbody/32, i+1); t->bodytab[i]=points.list; if (prev) Next(prev)=points.list; prev = points.tail; ADDV(cmr,cmr,points.cmr); ADDV(cmv,cmv,points.cmv); } #else for (i=31; i >=0 ; i--) { int processor= i/(32/__NumNodes); FUTURE(processor, nbody/32, i+1, uniform_testdata, &(fc[i])); } for (i=31; i>=0; i--) { TOUCH(&fc[i]); t->bodytab[i]=fc[i].value.list; ADDV(cmr,cmr,fc[i].value.cmr); ADDV(cmv,cmv,fc[i].value.cmv); } for (i=0; i<32; i++) { if (prev) Next(prev)=fc[i].value.list; prev = fc[i].value.tail; } #endif chatting("bodies created \n"); DIVVS(cmr, cmr, (real) nbody); /* normalize cm coords */ DIVVS(cmv, cmv, (real) nbody); for (tmp=0; tmp<MAX_NUM_NODES; tmp++) { bodiesper[tmp]=0; ptrper[tmp]=NULL; } /* Recall particles are in lists by processor */ for (p = t->bodytab[0]; p != NULL; p = Next(p)) { /* loop over particles */ icstruct xqic; SUBV(Pos(p), Pos(p), cmr); /* offset by cm coords */ SUBV(Vel(p), Vel(p), cmv); xqic = intcoord(p,t); tmp = ((old_subindex(xqic, IMAX >> 1) << NDIM) + old_subindex(xqic, IMAX >> 2)); tmp = tmp / range; bodiesper[tmp]++; Proc_Next(p) = ptrper[tmp]; ptrper[tmp] = p; Proc(p) = tmp; } MLOCAL(t); for (tmp=0; tmp<__NumNodes; tmp++) { chatting("Bodies per %d = %d\n",tmp ,bodiesper[tmp]); t->bodiesperproc[tmp]=ptrper[tmp]; } #ifndef PLAIN for (tmp=__NumNodes-1; tmp >=0 ; tmp--) { bodyptr list; list = t->bodiesperproc[tmp]; list = movebodies(list,tmp); t->bodiesperproc[tmp]=list; } #endif #ifdef MCC_DEBUG { int i=0; bodyptr p = t->bodytab[0]; for (; i < nbody; i++, p=Next(p)) printf("%d -- %f %f %f\n", i, Pos(p)[0], Pos(p)[1], Pos(p)[2]); } #endif tmp = 0; /* t1 = mclock(); */ tnow = 0.0; i = 0; /* Go through sequence of iterations */ nsteps = NSTEPS; /*chatting("About to perform %d iters from %f to %f by %f\n",*/ /*nsteps,tnow,tstop,dtime);*/ while ((tnow < tstop + 0.1*dtime) && (i < nsteps)) { stepsystem(t, i); /* advance N-body system */ tnow = tnow + dtime; /*chatting("tnow = %f sp = 0x%x\n", tnow, __getsp());*/ i++; } /* t2 = mclock(); */ #ifdef MCC_DEBUG { int i=0; bodyptr p = t->bodytab[0]; for (; i < nbody; i++, p=Next(p)) printf("%d -- %f %f %f\n", i, Pos(p)[0], Pos(p)[1], Pos(p)[2]); } #endif /* Return the tree to calling routine */ return t; } /* Use 1 of testdata, uniform_testdata */ /* * TESTDATA: generate Plummer model initial conditions for test runs, * scaled to units such that M = -4E = G = 1 (Henon, Hegge, etc). * See Aarseth, SJ, Henon, M, & Wielen, R (1974) Astr & Ap, 37, 183. */ #define MFRAC 0.999 /* mass cut off at MFRAC of total */ /* don't use this unless it is fixed on random numbers, &c */ bodyptr testdata() { real rsc, vsc, r, v, x, y; vector cmr, cmv; bodyptr head, p, prev; register i; double temp, t1; double seed = 123.0; register int k; double rsq, rsc1; real rad; assert(0,99); rsc = 3 * PI / 16; /* set length scale factor */ vsc = sqrt(1.0 / rsc); /* and recip. speed scale */ CLRV(cmr); /* init cm pos, vel */ CLRV(cmv); head = (bodyptr) ubody_alloc(0); prev = head; for (i = 0; i < nbody; i++) { /* loop over particles */ p = ubody_alloc(0); /* alloc space for bodies */ if (p == NULL) /* check space is available */ error("testdata: not enough memory\n"); /* if not, cry */ Next(prev) = p; prev = p; Type(p) = BODY; /* tag as a body */ Mass(p) = 1.0 / nbody; /* set masses equal */ seed = my_rand(seed); t1 = xrand(0.0, MFRAC, seed); temp = pow(t1, /* pick r in struct units */ -2.0/3.0) - 1; r = 1 / sqrt(temp); /* pick scaled position */ rad = rsc * r; do { /* pick point in NDIM-space */ for (k = 0; k < NDIM; k++) { /* loop over dimensions */ seed = my_rand(seed); Pos(p)[k] = xrand(-1.0, 1.0, seed); /* pick from unit cube */ } DOTVP(rsq, Pos(p), Pos(p)); /* compute radius squared */ } while (rsq > 1.0); /* reject if outside sphere */ rsc1 = rad / sqrt(rsq); /* compute scaling factor */ MULVS(Pos(p), Pos(p), rsc1); /* rescale to radius given */ ADDV(cmr, cmr, Pos(p)); /* add to running sum */ do { /* select from fn g(x) */ seed = my_rand(seed); x = xrand(0.0, 1.0, seed); /* for x in range 0:1 */ seed = my_rand(seed); y = xrand(0.0, 0.1, seed); /* max of g(x) is 0