| 1 | /****************************************************************************/ |
|---|
| 2 | /* TREEGRAV.C: routines to compute gravity. Public routines: gravcalc(). */ |
|---|
| 3 | /* Copyright (c) 2001 by Joshua E. Barnes, Honolulu, Hawai`i. */ |
|---|
| 4 | /****************************************************************************/ |
|---|
| 5 | |
|---|
| 6 | #include "stdinc.h" |
|---|
| 7 | #include "mathfns.h" |
|---|
| 8 | #include "vectmath.h" |
|---|
| 9 | #include "treedefs.h" |
|---|
| 10 | |
|---|
| 11 | /* Local routines to perform force calculations. */ |
|---|
| 12 | |
|---|
| 13 | local void walktree(nodeptr *, nodeptr *, cellptr, cellptr, |
|---|
| 14 | nodeptr, real, vector); |
|---|
| 15 | local bool accept(nodeptr, real, vector); |
|---|
| 16 | local void walksub(nodeptr *, nodeptr *, cellptr, cellptr, |
|---|
| 17 | nodeptr, real, vector); |
|---|
| 18 | local void gravsum(bodyptr, cellptr, cellptr); |
|---|
| 19 | local void sumnode(cellptr, cellptr, vector, real *, vector); |
|---|
| 20 | local void sumcell(cellptr, cellptr, vector, real *, vector); |
|---|
| 21 | |
|---|
| 22 | /* Lists of active nodes and interactions. */ |
|---|
| 23 | |
|---|
| 24 | #if !defined(FACTIVE) |
|---|
| 25 | # define FACTIVE 0.75 /* active list fudge factor */ |
|---|
| 26 | #endif |
|---|
| 27 | |
|---|
| 28 | local int actlen; /* length as allocated */ |
|---|
| 29 | |
|---|
| 30 | local nodeptr *active; /* list of nodes tested */ |
|---|
| 31 | |
|---|
| 32 | local cellptr interact; /* list of interactions */ |
|---|
| 33 | |
|---|
| 34 | /* |
|---|
| 35 | * GRAVCALC: perform force calculation on all particles. |
|---|
| 36 | */ |
|---|
| 37 | |
|---|
| 38 | void gravcalc(void) |
|---|
| 39 | { |
|---|
| 40 | double cpustart; |
|---|
| 41 | vector rmid; |
|---|
| 42 | |
|---|
| 43 | actlen = FACTIVE * 216 * tdepth; /* estimate list length */ |
|---|
| 44 | #if !defined(QUICKSCAN) |
|---|
| 45 | actlen = actlen * rpow(theta, -2.5); /* allow for opening angle */ |
|---|
| 46 | #endif |
|---|
| 47 | active = (nodeptr *) allocate(actlen * sizeof(nodeptr)); |
|---|
| 48 | interact = (cellptr) allocate(actlen * sizeof(cell)); |
|---|
| 49 | cpustart = cputime(); /* record time, less alloc */ |
|---|
| 50 | actmax = nbbcalc = nbccalc = 0; /* zero cumulative counters */ |
|---|
| 51 | active[0] = (nodeptr) root; /* initialize active list */ |
|---|
| 52 | CLRV(rmid); /* set center of root cell */ |
|---|
| 53 | walktree(active, active + 1, interact, interact + actlen, |
|---|
| 54 | (nodeptr) root, rsize, rmid); /* scan tree, update forces */ |
|---|
| 55 | cpuforce = cputime() - cpustart; /* store CPU time w/o alloc */ |
|---|
| 56 | free(active); |
|---|
| 57 | free(interact); |
|---|
| 58 | } |
|---|
| 59 | |
|---|
| 60 | /* |
|---|
| 61 | * WALKTREE: do a complete walk of the tree, building the interaction |
|---|
| 62 | * list level-by-level and computing the resulting force on each body. |
|---|
| 63 | */ |
|---|
| 64 | |
|---|
| 65 | local void walktree(nodeptr *aptr, nodeptr *nptr, cellptr cptr, cellptr bptr, |
|---|
| 66 | nodeptr p, real psize, vector pmid) |
|---|
| 67 | { |
|---|
| 68 | nodeptr *np, *ap, q; |
|---|
| 69 | int actsafe; |
|---|
| 70 | |
|---|
| 71 | if (Update(p)) { /* are new forces needed? */ |
|---|
| 72 | np = nptr; /* start new active list */ |
|---|
| 73 | actsafe = actlen - NSUB; /* leave room for NSUB more */ |
|---|
| 74 | for (ap = aptr; ap < nptr; ap++) /* loop over active nodes */ |
|---|
| 75 | if (Type(*ap) == CELL) { /* is this node a cell? */ |
|---|
| 76 | if (accept(*ap, psize, pmid)) { /* does it pass the test? */ |
|---|
| 77 | Mass(cptr) = Mass(*ap); /* copy to interaction list */ |
|---|
| 78 | SETV(Pos(cptr), Pos(*ap)); |
|---|
| 79 | SETM(Quad(cptr), Quad(*ap)); |
|---|
| 80 | cptr++; /* and bump cell array ptr */ |
|---|
| 81 | } else { /* else it fails the test */ |
|---|
| 82 | if (np - active >= actsafe) /* check list has room */ |
|---|
| 83 | error("walktree: active list overflow\n"); |
|---|
| 84 | for (q = More(*ap); q != Next(*ap); q = Next(q)) |
|---|
| 85 | /* loop over all subcells */ |
|---|
| 86 | *np++= q; /* put on new active list */ |
|---|
| 87 | } |
|---|
| 88 | } else /* else this node is a body */ |
|---|
| 89 | if (*ap != p) { /* if not self-interaction */ |
|---|
| 90 | --bptr; /* bump body array ptr */ |
|---|
| 91 | Mass(bptr) = Mass(*ap); /* and copy data to array */ |
|---|
| 92 | SETV(Pos(bptr), Pos(*ap)); |
|---|
| 93 | } |
|---|
| 94 | actmax = MAX(actmax, np - active); /* keep track of max active */ |
|---|
| 95 | if (np != nptr) /* if new actives listed */ |
|---|
| 96 | walksub(nptr, np, cptr, bptr, p, psize, pmid); |
|---|
| 97 | /* then visit next level */ |
|---|
| 98 | else { /* else no actives left, so */ |
|---|
| 99 | if (Type(p) != BODY) /* must have found a body */ |
|---|
| 100 | error("walktree: recursion terminated with cell\n"); |
|---|
| 101 | gravsum((bodyptr) p, cptr, bptr); /* sum force on the body */ |
|---|
| 102 | } |
|---|
| 103 | } |
|---|
| 104 | } |
|---|
| 105 | |
|---|
| 106 | #if defined(QUICKSCAN) |
|---|
| 107 | |
|---|
| 108 | /* |
|---|
| 109 | * ACCEPT: quick criterion accepts any cell not touching cell p. |
|---|
| 110 | */ |
|---|
| 111 | |
|---|
| 112 | local bool accept(nodeptr c, real psize, vector pmid) |
|---|
| 113 | { |
|---|
| 114 | real p15, dk; |
|---|
| 115 | |
|---|
| 116 | p15 = ((real) 1.5) * psize; /* premultiply cell size */ |
|---|
| 117 | dk = Pos(c)[0] - pmid[0]; /* find distance to midpnt */ |
|---|
| 118 | if (ABS(dk) > p15) /* if c does not touch p */ |
|---|
| 119 | return (TRUE); /* then accept interaction */ |
|---|
| 120 | dk = Pos(c)[1] - pmid[1]; /* find distance to midpnt */ |
|---|
| 121 | if (ABS(dk) > p15) /* if c does not touch p */ |
|---|
| 122 | return (TRUE); /* then accept interaction */ |
|---|
| 123 | dk = Pos(c)[2] - pmid[2]; /* find distance to midpnt */ |
|---|
| 124 | if (ABS(dk) > p15) /* if c does not touch p */ |
|---|
| 125 | return (TRUE); /* then accept interaction */ |
|---|
| 126 | return (FALSE); /* else do not accept it */ |
|---|
| 127 | } |
|---|
| 128 | |
|---|
| 129 | #else |
|---|
| 130 | |
|---|
| 131 | /* |
|---|
| 132 | * ACCEPT: standard criterion accepts cell if its critical radius |
|---|
| 133 | * does not intersect cell p, and also imposes above condition. |
|---|
| 134 | */ |
|---|
| 135 | |
|---|
| 136 | local bool accept(nodeptr c, real psize, vector pmid) |
|---|
| 137 | { |
|---|
| 138 | real dmax, dsq, dk; |
|---|
| 139 | int k; |
|---|
| 140 | |
|---|
| 141 | dmax = psize; /* init maximum distance */ |
|---|
| 142 | dsq = 0.0; /* and squared min distance */ |
|---|
| 143 | for (k = 0; k < NDIM; k++) { /* loop over space dims */ |
|---|
| 144 | dk = Pos(c)[k] - pmid[k]; /* form distance to midpnt */ |
|---|
| 145 | if (dk < 0) /* and get absolute value */ |
|---|
| 146 | dk = - dk; |
|---|
| 147 | if (dk > dmax) /* keep track of max value */ |
|---|
| 148 | dmax = dk; |
|---|
| 149 | dk -= ((real) 0.5) * psize; /* allow for size of cell */ |
|---|
| 150 | if (dk > 0) |
|---|
| 151 | dsq += dk * dk; /* sum min dist to cell ^2 */ |
|---|
| 152 | } |
|---|
| 153 | return (dsq > Rcrit2(c) && /* test angular criterion */ |
|---|
| 154 | dmax > ((real) 1.5) * psize); /* and adjacency criterion */ |
|---|
| 155 | } |
|---|
| 156 | |
|---|
| 157 | #endif |
|---|
| 158 | |
|---|
| 159 | /* |
|---|
| 160 | * WALKSUB: test next level's active list against subnodes of p. |
|---|
| 161 | */ |
|---|
| 162 | |
|---|
| 163 | local void walksub(nodeptr *nptr, nodeptr *np, cellptr cptr, cellptr bptr, |
|---|
| 164 | nodeptr p, real psize, vector pmid) |
|---|
| 165 | { |
|---|
| 166 | real poff; |
|---|
| 167 | nodeptr q; |
|---|
| 168 | int k; |
|---|
| 169 | vector nmid; |
|---|
| 170 | |
|---|
| 171 | poff = psize / 4; /* precompute mid. offset */ |
|---|
| 172 | if (Type(p) == CELL) { /* fanout over descendents */ |
|---|
| 173 | for (q = More(p); q != Next(p); q = Next(q)) { |
|---|
| 174 | /* loop over all subcells */ |
|---|
| 175 | for (k = 0; k < NDIM; k++) /* locate each's midpoint */ |
|---|
| 176 | nmid[k] = pmid[k] + (Pos(q)[k] < pmid[k] ? - poff : poff); |
|---|
| 177 | walktree(nptr, np, cptr, bptr, q, psize / 2, nmid); |
|---|
| 178 | /* recurse on subcell */ |
|---|
| 179 | } |
|---|
| 180 | } else { /* extend virtual tree */ |
|---|
| 181 | for (k = 0; k < NDIM; k++) /* locate next midpoint */ |
|---|
| 182 | nmid[k] = pmid[k] + (Pos(p)[k] < pmid[k] ? - poff : poff); |
|---|
| 183 | walktree(nptr, np, cptr, bptr, p, psize / 2, nmid); |
|---|
| 184 | /* and search next level */ |
|---|
| 185 | } |
|---|
| 186 | } |
|---|
| 187 | |
|---|
| 188 | /* |
|---|
| 189 | * GRAVSUM: compute gravitational field at body p0. |
|---|
| 190 | */ |
|---|
| 191 | |
|---|
| 192 | local void gravsum(bodyptr p0, cellptr cptr, cellptr bptr) |
|---|
| 193 | { |
|---|
| 194 | vector pos0, acc0; |
|---|
| 195 | real phi0; |
|---|
| 196 | |
|---|
| 197 | SETV(pos0, Pos(p0)); /* copy position of body */ |
|---|
| 198 | phi0 = 0.0; /* init total potential */ |
|---|
| 199 | CLRV(acc0); /* and total acceleration */ |
|---|
| 200 | if (usequad) /* if using quad moments */ |
|---|
| 201 | sumcell(interact, cptr, pos0, &phi0, acc0); |
|---|
| 202 | /* sum cell forces w quads */ |
|---|
| 203 | else /* not using quad moments */ |
|---|
| 204 | sumnode(interact, cptr, pos0, &phi0, acc0); |
|---|
| 205 | /* sum cell forces wo quads */ |
|---|
| 206 | sumnode(bptr, interact + actlen, pos0, &phi0, acc0); |
|---|
| 207 | /* sum forces from bodies */ |
|---|
| 208 | Phi(p0) = phi0; /* store total potential */ |
|---|
| 209 | SETV(Acc(p0), acc0); /* and total acceleration */ |
|---|
| 210 | nbbcalc += interact + actlen - bptr; /* count body-body forces */ |
|---|
| 211 | nbccalc += cptr - interact; /* count body-cell forces */ |
|---|
| 212 | } |
|---|
| 213 | |
|---|
| 214 | /* |
|---|
| 215 | * SUMNODE: add up body-node interactions. |
|---|
| 216 | */ |
|---|
| 217 | |
|---|
| 218 | local void sumnode(cellptr start, cellptr finish, |
|---|
| 219 | vector pos0, real *phi0, vector acc0) |
|---|
| 220 | { |
|---|
| 221 | cellptr p; |
|---|
| 222 | real eps2, dr2, drab, phi_p, mr3i; |
|---|
| 223 | vector dr; |
|---|
| 224 | |
|---|
| 225 | eps2 = eps * eps; /* avoid extra multiplys */ |
|---|
| 226 | for (p = start; p < finish; p++) { /* loop over node list */ |
|---|
| 227 | DOTPSUBV(dr2, dr, Pos(p), pos0); /* compute separation */ |
|---|
| 228 | /* and distance squared */ |
|---|
| 229 | dr2 += eps2; /* add standard softening */ |
|---|
| 230 | drab = rsqrt(dr2); /* form scalar "distance" */ |
|---|
| 231 | phi_p = Mass(p) / drab; /* get partial potential */ |
|---|
| 232 | *phi0 -= phi_p; /* decrement tot potential */ |
|---|
| 233 | mr3i = phi_p / dr2; /* form scale factor for dr */ |
|---|
| 234 | ADDMULVS(acc0, dr, mr3i); /* sum partial acceleration */ |
|---|
| 235 | } |
|---|
| 236 | } |
|---|
| 237 | |
|---|
| 238 | /* |
|---|
| 239 | * SUMCELL: add up body-cell interactions. |
|---|
| 240 | */ |
|---|
| 241 | |
|---|
| 242 | local void sumcell(cellptr start, cellptr finish, |
|---|
| 243 | vector pos0, real *phi0, vector acc0) |
|---|
| 244 | { |
|---|
| 245 | cellptr p; |
|---|
| 246 | real eps2, dr2, drab, phi_p, mr3i, drqdr, dr5i, phi_q; |
|---|
| 247 | vector dr, qdr; |
|---|
| 248 | |
|---|
| 249 | eps2 = eps * eps; |
|---|
| 250 | for (p = start; p < finish; p++) { /* loop over node list */ |
|---|
| 251 | DOTPSUBV(dr2, dr, Pos(p), pos0); /* do mono part of force */ |
|---|
| 252 | dr2 += eps2; |
|---|
| 253 | drab = rsqrt(dr2); |
|---|
| 254 | phi_p = Mass(p) / drab; |
|---|
| 255 | mr3i = phi_p / dr2; |
|---|
| 256 | DOTPMULMV(drqdr, qdr, Quad(p), dr); /* do quad part of force */ |
|---|
| 257 | dr5i = ((real) 1.0) / (dr2 * dr2 * drab); |
|---|
| 258 | phi_q = ((real) 0.5) * dr5i * drqdr; |
|---|
| 259 | *phi0 -= phi_p + phi_q; /* add mono and quad pot */ |
|---|
| 260 | mr3i += ((real) 5.0) * phi_q / dr2; |
|---|
| 261 | ADDMULVS2(acc0, dr, mr3i, qdr, -dr5i); /* add mono and quad acc */ |
|---|
| 262 | } |
|---|
| 263 | } |
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/PP2009/chrome/site/your_project_logo.png){kind=logo}
