/****************************************************************************/ /* TREEGRAV.C: routines to compute gravity. Public routines: gravcalc(). */ /* Copyright (c) 2001 by Joshua E. Barnes, Honolulu, Hawai`i. */ /****************************************************************************/ #include "stdinc.h" #include "mathfns.h" #include "vectmath.h" #include "treedefs.h" /* Local routines to perform force calculations. */ local void walktree(nodeptr *, nodeptr *, cellptr, cellptr, nodeptr, real, vector); local bool accept(nodeptr, real, vector); local void walksub(nodeptr *, nodeptr *, cellptr, cellptr, nodeptr, real, vector); local void gravsum(bodyptr, cellptr, cellptr); local void sumnode(cellptr, cellptr, vector, real *, vector); local void sumcell(cellptr, cellptr, vector, real *, vector); /* Lists of active nodes and interactions. */ #if !defined(FACTIVE) # define FACTIVE 0.75 /* active list fudge factor */ #endif local int actlen; /* length as allocated */ local nodeptr *active; /* list of nodes tested */ local cellptr interact; /* list of interactions */ /* * GRAVCALC: perform force calculation on all particles. */ void gravcalc(void) { double cpustart; vector rmid; actlen = FACTIVE * 216 * tdepth; /* estimate list length */ #if !defined(QUICKSCAN) actlen = actlen * rpow(theta, -2.5); /* allow for opening angle */ #endif active = (nodeptr *) allocate(actlen * sizeof(nodeptr)); interact = (cellptr) allocate(actlen * sizeof(cell)); cpustart = cputime(); /* record time, less alloc */ actmax = nbbcalc = nbccalc = 0; /* zero cumulative counters */ active[0] = (nodeptr) root; /* initialize active list */ CLRV(rmid); /* set center of root cell */ walktree(active, active + 1, interact, interact + actlen, (nodeptr) root, rsize, rmid); /* scan tree, update forces */ cpuforce = cputime() - cpustart; /* store CPU time w/o alloc */ free(active); free(interact); } /* * WALKTREE: do a complete walk of the tree, building the interaction * list level-by-level and computing the resulting force on each body. */ local void walktree(nodeptr *aptr, nodeptr *nptr, cellptr cptr, cellptr bptr, nodeptr p, real psize, vector pmid) { nodeptr *np, *ap, q; int actsafe; if (Update(p)) { /* are new forces needed? */ np = nptr; /* start new active list */ actsafe = actlen - NSUB; /* leave room for NSUB more */ for (ap = aptr; ap < nptr; ap++) /* loop over active nodes */ if (Type(*ap) == CELL) { /* is this node a cell? */ if (accept(*ap, psize, pmid)) { /* does it pass the test? */ Mass(cptr) = Mass(*ap); /* copy to interaction list */ SETV(Pos(cptr), Pos(*ap)); SETM(Quad(cptr), Quad(*ap)); cptr++; /* and bump cell array ptr */ } else { /* else it fails the test */ if (np - active >= actsafe) /* check list has room */ error("walktree: active list overflow\n"); for (q = More(*ap); q != Next(*ap); q = Next(q)) /* loop over all subcells */ *np++= q; /* put on new active list */ } } else /* else this node is a body */ if (*ap != p) { /* if not self-interaction */ --bptr; /* bump body array ptr */ Mass(bptr) = Mass(*ap); /* and copy data to array */ SETV(Pos(bptr), Pos(*ap)); } actmax = MAX(actmax, np - active); /* keep track of max active */ if (np != nptr) /* if new actives listed */ walksub(nptr, np, cptr, bptr, p, psize, pmid); /* then visit next level */ else { /* else no actives left, so */ if (Type(p) != BODY) /* must have found a body */ error("walktree: recursion terminated with cell\n"); gravsum((bodyptr) p, cptr, bptr); /* sum force on the body */ } } } #if defined(QUICKSCAN) /* * ACCEPT: quick criterion accepts any cell not touching cell p. */ local bool accept(nodeptr c, real psize, vector pmid) { real p15, dk; p15 = ((real) 1.5) * psize; /* premultiply cell size */ dk = Pos(c)[0] - pmid[0]; /* find distance to midpnt */ if (ABS(dk) > p15) /* if c does not touch p */ return (TRUE); /* then accept interaction */ dk = Pos(c)[1] - pmid[1]; /* find distance to midpnt */ if (ABS(dk) > p15) /* if c does not touch p */ return (TRUE); /* then accept interaction */ dk = Pos(c)[2] - pmid[2]; /* find distance to midpnt */ if (ABS(dk) > p15) /* if c does not touch p */ return (TRUE); /* then accept interaction */ return (FALSE); /* else do not accept it */ } #else /* * ACCEPT: standard criterion accepts cell if its critical radius * does not intersect cell p, and also imposes above condition. */ local bool accept(nodeptr c, real psize, vector pmid) { real dmax, dsq, dk; int k; dmax = psize; /* init maximum distance */ dsq = 0.0; /* and squared min distance */ for (k = 0; k < NDIM; k++) { /* loop over space dims */ dk = Pos(c)[k] - pmid[k]; /* form distance to midpnt */ if (dk < 0) /* and get absolute value */ dk = - dk; if (dk > dmax) /* keep track of max value */ dmax = dk; dk -= ((real) 0.5) * psize; /* allow for size of cell */ if (dk > 0) dsq += dk * dk; /* sum min dist to cell ^2 */ } return (dsq > Rcrit2(c) && /* test angular criterion */ dmax > ((real) 1.5) * psize); /* and adjacency criterion */ } #endif /* * WALKSUB: test next level's active list against subnodes of p. */ local void walksub(nodeptr *nptr, nodeptr *np, cellptr cptr, cellptr bptr, nodeptr p, real psize, vector pmid) { real poff; nodeptr q; int k; vector nmid; poff = psize / 4; /* precompute mid. offset */ if (Type(p) == CELL) { /* fanout over descendents */ for (q = More(p); q != Next(p); q = Next(q)) { /* loop over all subcells */ for (k = 0; k < NDIM; k++) /* locate each's midpoint */ nmid[k] = pmid[k] + (Pos(q)[k] < pmid[k] ? - poff : poff); walktree(nptr, np, cptr, bptr, q, psize / 2, nmid); /* recurse on subcell */ } } else { /* extend virtual tree */ for (k = 0; k < NDIM; k++) /* locate next midpoint */ nmid[k] = pmid[k] + (Pos(p)[k] < pmid[k] ? - poff : poff); walktree(nptr, np, cptr, bptr, p, psize / 2, nmid); /* and search next level */ } } /* * GRAVSUM: compute gravitational field at body p0. */ local void gravsum(bodyptr p0, cellptr cptr, cellptr bptr) { vector pos0, acc0; real phi0; SETV(pos0, Pos(p0)); /* copy position of body */ phi0 = 0.0; /* init total potential */ CLRV(acc0); /* and total acceleration */ if (usequad) /* if using quad moments */ sumcell(interact, cptr, pos0, &phi0, acc0); /* sum cell forces w quads */ else /* not using quad moments */ sumnode(interact, cptr, pos0, &phi0, acc0); /* sum cell forces wo quads */ sumnode(bptr, interact + actlen, pos0, &phi0, acc0); /* sum forces from bodies */ Phi(p0) = phi0; /* store total potential */ SETV(Acc(p0), acc0); /* and total acceleration */ nbbcalc += interact + actlen - bptr; /* count body-body forces */ nbccalc += cptr - interact; /* count body-cell forces */ } /* * SUMNODE: add up body-node interactions. */ local void sumnode(cellptr start, cellptr finish, vector pos0, real *phi0, vector acc0) { cellptr p; real eps2, dr2, drab, phi_p, mr3i; vector dr; eps2 = eps * eps; /* avoid extra multiplys */ for (p = start; p < finish; p++) { /* loop over node list */ DOTPSUBV(dr2, dr, Pos(p), pos0); /* compute separation */ /* and distance squared */ dr2 += eps2; /* add standard softening */ drab = rsqrt(dr2); /* form scalar "distance" */ phi_p = Mass(p) / drab; /* get partial potential */ *phi0 -= phi_p; /* decrement tot potential */ mr3i = phi_p / dr2; /* form scale factor for dr */ ADDMULVS(acc0, dr, mr3i); /* sum partial acceleration */ } } /* * SUMCELL: add up body-cell interactions. */ local void sumcell(cellptr start, cellptr finish, vector pos0, real *phi0, vector acc0) { cellptr p; real eps2, dr2, drab, phi_p, mr3i, drqdr, dr5i, phi_q; vector dr, qdr; eps2 = eps * eps; for (p = start; p < finish; p++) { /* loop over node list */ DOTPSUBV(dr2, dr, Pos(p), pos0); /* do mono part of force */ dr2 += eps2; drab = rsqrt(dr2); phi_p = Mass(p) / drab; mr3i = phi_p / dr2; DOTPMULMV(drqdr, qdr, Quad(p), dr); /* do quad part of force */ dr5i = ((real) 1.0) / (dr2 * dr2 * drab); phi_q = ((real) 0.5) * dr5i * drqdr; *phi0 -= phi_p + phi_q; /* add mono and quad pot */ mr3i += ((real) 5.0) * phi_q / dr2; ADDMULVS2(acc0, dr, mr3i, qdr, -dr5i); /* add mono and quad acc */ } }