changa 3.5
Loading...
Searching...
No Matches
gravity.h
1#ifndef __GRAVITY_H__
2#define __GRAVITY_H__
3
4#include "TreeNode.h"
5#include "GenericTreeNode.h"
6#include "Space.h"
7#include "SSEdefs.h"
8
9extern cosmoType theta;
10extern cosmoType thetaMono;
11
12/*
13** see (A1) and (A2) of TREESPH: A UNIFICATION OF SPH WITH THE
14** HIERARCHICAL TREE METHOD by Lars Hernquist and Neal Katz.
15** APJ Supplemant Series 70:416-446, 1989
16**
17** Higher derivative terms c and d for use with quadrupole spline
18** softening (Joachim Stadel, Dec. 94).
19*/
20#if !CMK_SSE
21inline
22void SPLINEQ(cosmoType invr,cosmoType r2,cosmoType twoh,cosmoType& a,
23 cosmoType& b,cosmoType& c,cosmoType& d)
24{
25 cosmoType u,dih,dir=(invr);
26 if ((r2) < (twoh)*(twoh)) {
27 dih = COSMO_CONST(2.0)/(twoh);
28 u = dih/dir;
29 if (u < COSMO_CONST(1.0)) {
30 a = dih*(COSMO_CONST(7.0)/COSMO_CONST(5.0)
31 - COSMO_CONST(2.0)/COSMO_CONST(3.0)*u*u
32 + COSMO_CONST(3.0)/COSMO_CONST(10.0)*u*u*u*u
33 - COSMO_CONST(1.0)/COSMO_CONST(10.0)*u*u*u*u*u);
34 b = dih*dih*dih*(COSMO_CONST(4.0)/COSMO_CONST(3.0)
35 - COSMO_CONST(6.0)/COSMO_CONST(5.0)*u*u
36 + COSMO_CONST(1.0)/COSMO_CONST(2.0)*u*u*u);
37 c = dih*dih*dih*dih*dih*(COSMO_CONST(12.0)/COSMO_CONST(5.0)
38 - COSMO_CONST(3.0)/COSMO_CONST(2.0)*u);
39 d = COSMO_CONST(3.0)/COSMO_CONST(2.0)*dih*dih*dih*dih*dih*dih*dir;
40 }
41 else {
42 a = COSMO_CONST(-1.0)/COSMO_CONST(15.0)*dir
43 + dih*(COSMO_CONST(8.0)/COSMO_CONST(5.0)
44 - COSMO_CONST(4.0)/COSMO_CONST(3.0)*u*u + u*u*u
45 - COSMO_CONST(3.0)/COSMO_CONST(10.0)*u*u*u*u
46 + COSMO_CONST(1.0)/COSMO_CONST(30.0)*u*u*u*u*u);
47 b = COSMO_CONST(-1.0)/COSMO_CONST(15.0)*dir*dir*dir
48 + dih*dih*dih*(COSMO_CONST(8.0)/COSMO_CONST(3.0) - COSMO_CONST(3.0)*u
49 + COSMO_CONST(6.0)/COSMO_CONST(5.0)*u*u
50 - COSMO_CONST(1.0)/COSMO_CONST(6.0)*u*u*u);
51 c = COSMO_CONST(-1.0)/COSMO_CONST(5.0)*dir*dir*dir*dir*dir
52 + COSMO_CONST(3.0)*dih*dih*dih*dih*dir
53 + dih*dih*dih*dih*dih*(COSMO_CONST(-12.0)/COSMO_CONST(5.0)
54 + COSMO_CONST(1.0)/COSMO_CONST(2.0)*u);
55 d = -dir*dir*dir*dir*dir*dir*dir
56 + COSMO_CONST(3.0)*dih*dih*dih*dih*dir*dir*dir
57 - COSMO_CONST(1.0)/COSMO_CONST(2.0)*dih*dih*dih*dih*dih*dih*dir;
58 }
59 }
60 else {
61 a = dir;
62 b = a*a*a;
63 c = COSMO_CONST(3.0)*b*a*a;
64 d = COSMO_CONST(5.0)*c*a*a;
65 }
66}
67#else
68inline
69void SPLINEQ(SSEcosmoType invr, SSEcosmoType r2, SSEcosmoType twoh,
70 SSEcosmoType& a, SSEcosmoType& b,
71 SSEcosmoType& c, SSEcosmoType& d)
72{
73 SSEcosmoType dir;
74 dir = invr;
75 SSEcosmoType select0 = r2 < twoh * twoh;
76 int compare0 = movemask(select0);
77 // make common case fast, at the cost of some redundant code
78 // in the infrequent case
79 if (!compare0) {
80 a = dir;
81 b = a*a*a;
82 c = COSMO_CONST(3.0)*b*a*a;
83 d = COSMO_CONST(5.0)*c*a*a;
84 }
85 else {
86 SSEcosmoType u, dih;
87 SSEcosmoType select1;
88 SSEcosmoType a0,b0,c0,d0,a1,b1,c1,d1,a2,b2,c2,d2;
89 int compare1;
90
91 dih = COSMO_CONST(2.0)/twoh;
92 u = dih/dir;
93
94 if ((~compare0) & cosmoMask) {
95 a0 = dir;
96 b0 = a0*a0*a0;
97 c0 = COSMO_CONST(3.0)*b0*a0*a0;
98 d0 = COSMO_CONST(5.0)*c0*a0*a0;
99 }
100
101 select1 = u < COSMO_CONST(1.0);
102 compare1 = movemask(select1);
103 if (compare1) {
104 a1 = dih*(COSMO_CONST(7.0)/COSMO_CONST(5.0)
105 - COSMO_CONST(2.0)/COSMO_CONST(3.0)*u*u
106 + COSMO_CONST(3.0)/COSMO_CONST(10.0)*u*u*u*u
107 - COSMO_CONST(1.0)/COSMO_CONST(10.0)*u*u*u*u*u);
108 b1 = dih*dih*dih*(COSMO_CONST(4.0)/COSMO_CONST(3.0)
109 - COSMO_CONST(6.0)/COSMO_CONST(5.0)*u*u
110 + COSMO_CONST(1.0)/COSMO_CONST(2.0)*u*u*u);
111 c1 = dih*dih*dih*dih*dih*(COSMO_CONST(12.0)/COSMO_CONST(5.0)
112 - COSMO_CONST(3.0)/COSMO_CONST(2.0)*u);
113 d1 = COSMO_CONST(3.0)/COSMO_CONST(2.0)*dih*dih*dih*dih*dih*dih*dir;
114 }
115 if ((~compare1) & cosmoMask) {
116 a2 = COSMO_CONST(-1.0)/COSMO_CONST(15.0)*dir
117 + dih*(COSMO_CONST(8.0)/COSMO_CONST(5.0)
118 - COSMO_CONST(4.0)/COSMO_CONST(3.0)*u*u + u*u*u
119 - COSMO_CONST(3.0)/COSMO_CONST(10.0)*u*u*u*u
120 + COSMO_CONST(1.0)/COSMO_CONST(30.0)*u*u*u*u*u);
121 b2 = COSMO_CONST(-1.0)/COSMO_CONST(15.0)*dir*dir*dir
122 + dih*dih*dih*(COSMO_CONST(8.0)/COSMO_CONST(3.0) - COSMO_CONST(3.0)*u
123 + COSMO_CONST(6.0)/COSMO_CONST(5.0)*u*u
124 - COSMO_CONST(1.0)/COSMO_CONST(6.0)*u*u*u);
125 c2 = COSMO_CONST(-1.0)/COSMO_CONST(5.0)*dir*dir*dir*dir*dir
126 + COSMO_CONST(3.0)*dih*dih*dih*dih*dir
127 + dih*dih*dih*dih*dih*(COSMO_CONST(-12.0)/COSMO_CONST(5.0)
128 + COSMO_CONST(1.0)/COSMO_CONST(2.0)*u);
129 d2 = -dir*dir*dir*dir*dir*dir*dir
130 + COSMO_CONST(3.0)*dih*dih*dih*dih*dir*dir*dir
131 - COSMO_CONST(1.0)/COSMO_CONST(2.0)*dih*dih*dih*dih*dih*dih*dir;
132 }
133
134 a = andnot(select0, a0)
135 | (select0 & ((select1 & a1) | andnot(select1, a2)));
136 b = andnot(select0, b0)
137 | (select0 & ((select1 & b1) | andnot(select1, b2)));
138 c = andnot(select0, c0)
139 | (select0 & ((select1 & c1) | andnot(select1, c2)));
140 d = andnot(select0, d0)
141 | (select0 & ((select1 & d1) | andnot(select1, d2)));
142 }
143}
144#endif
145
146#if !CMK_SSE
147inline
148void SPLINE(cosmoType r2, cosmoType twoh, cosmoType &a, cosmoType &b)
149{
150 cosmoType r, u,dih,dir;
151 r = sqrt(r2);
152
153 if (r < (twoh)) {
154 dih = COSMO_CONST(2.0)/(twoh);
155 u = r*dih;
156 if (u < COSMO_CONST(1.0)) {
157 a = dih*(COSMO_CONST(7.0)/COSMO_CONST(5.0)
158 - COSMO_CONST(2.0)/COSMO_CONST(3.0)*u*u
159 + COSMO_CONST(3.0)/COSMO_CONST(10.0)*u*u*u*u
160 - COSMO_CONST(1.0)/COSMO_CONST(10.0)*u*u*u*u*u);
161 b = dih*dih*dih*(COSMO_CONST(4.0)/COSMO_CONST(3.0)
162 - COSMO_CONST(6.0)/COSMO_CONST(5.0)*u*u
163 + COSMO_CONST(1.0)/COSMO_CONST(2.0)*u*u*u);
164 }
165 else {
166 dir = COSMO_CONST(1.0)/r;
167 a = COSMO_CONST(-1.0)/COSMO_CONST(15.0)*dir
168 + dih*(COSMO_CONST(8.0)/COSMO_CONST(5.0)
169 - COSMO_CONST(4.0)/COSMO_CONST(3.0)*u*u + u*u*u
170 - COSMO_CONST(3.0)/COSMO_CONST(10.0)*u*u*u*u
171 + COSMO_CONST(1.0)/COSMO_CONST(30.0)*u*u*u*u*u);
172 b = COSMO_CONST(-1.0)/COSMO_CONST(15.0)*dir*dir*dir
173 + dih*dih*dih*(COSMO_CONST(8.0)/COSMO_CONST(3.0) - COSMO_CONST(3.0)*u
174 + COSMO_CONST(6.0)/COSMO_CONST(5.0)*u*u
175 - COSMO_CONST(1.0)/COSMO_CONST(6.0)*u*u*u);
176 }
177 }
178 else {
179 a = COSMO_CONST(1.0)/r;
180 b = a*a*a;
181 }
182}
183#else
184inline
185void SPLINE(SSEcosmoType r2, SSEcosmoType twoh,
186 SSEcosmoType &a, SSEcosmoType &b)
187{
188 SSEcosmoType r;
189 r = sqrt(r2);
190 SSEcosmoType select0 = r < twoh;
191 int compare0 = movemask(select0);
192 // make common case fast, at the cost of some redundant code
193 // in the infrequent case
194 if (!compare0) {
195 a = COSMO_CONST(1.0)/r;
196 b = a * a * a;
197 }
198 else {
199 SSEcosmoType u, dih, dir;
200 SSEcosmoType a0, b0, a1, b1, a2, b2;
201 SSEcosmoType select1;
202 int compare1;
203
204 dih = COSMO_CONST(2.0)/twoh;
205 u = r * dih;
206
207 if ((~compare0) & cosmoMask) {
208 a0 = COSMO_CONST(1.0)/r;
209 b0 = a0 * a0 * a0;
210 }
211
212 select1 = u < COSMO_CONST(1.0);
213 compare1 = movemask(select1);
214 if (compare1) {
215 a1 = dih*(COSMO_CONST(7.0)/COSMO_CONST(5.0)
216 - COSMO_CONST(2.0)/COSMO_CONST(3.0)*u*u
217 + COSMO_CONST(3.0)/COSMO_CONST(10.0)*u*u*u*u
218 - COSMO_CONST(1.0)/COSMO_CONST(10.0)*u*u*u*u*u);
219 b1 = dih*dih*dih*(COSMO_CONST(4.0)/COSMO_CONST(3.0)
220 - COSMO_CONST(6.0)/COSMO_CONST(5.0)*u*u
221 + COSMO_CONST(1.0)/COSMO_CONST(2.0)*u*u*u);
222 }
223 if ((~compare1) & cosmoMask) {
224 dir = COSMO_CONST(1.0)/r;
225 a2 = COSMO_CONST(-1.0)/COSMO_CONST(15.0)*dir
226 + dih*(COSMO_CONST(8.0)/COSMO_CONST(5.0)
227 - COSMO_CONST(4.0)/COSMO_CONST(3.0)*u*u + u*u*u
228 - COSMO_CONST(3.0)/COSMO_CONST(10.0)*u*u*u*u
229 + COSMO_CONST(1.0)/COSMO_CONST(30.0)*u*u*u*u*u);
230 b2 = COSMO_CONST(-1.0)/COSMO_CONST(15.0)*dir*dir*dir
231 + dih*dih*dih*(COSMO_CONST(8.0)/COSMO_CONST(3.0) - COSMO_CONST(3.0)*u
232 + COSMO_CONST(6.0)/COSMO_CONST(5.0)*u*u
233 - COSMO_CONST(1.0)/COSMO_CONST(6.0)*u*u*u);
234 }
235
236 a = andnot(select0, a0)
237 | (select0 & ((select1 & a1) | andnot(select1, a2)));
238 b = andnot(select0, b0)
239 | (select0 & ((select1 & b1) | andnot(select1, b2)));
240 }
241}
242#endif
243
244//
245// Return true if the soften nodes overlap, or if the source node's
246// softening overlaps the bounding box; i.e. the forces involve softening
247// @param node Node to test
248// @param myNode
249// @param offset Periodic offset applied to node
250//
251inline
252int openSoftening(Tree::GenericTreeNode *node, Tree::GenericTreeNode *myNode,
253 Vector3D<cosmoType> offset)
254{
255 Sphere<cosmoType> s(node->moments.cm + offset, 2.0*node->moments.soft);
256 Sphere<cosmoType> myS(myNode->moments.cm, 2.0*myNode->moments.soft);
257 if(Space::intersect(myS, s))
258 return true;
259 return Space::intersect(myNode->boundingBox, s);
260}
261
262#ifdef CMK_VERSION_BLUEGENE
263static int forProgress = 0;
264#endif
265
266#if !CMK_SSE
267inline int partBucketForce(ExternalGravityParticle *part,
268 Tree::GenericTreeNode *req,
269 GravityParticle *particles,
270 Vector3D<cosmoType> offset, int activeRung) {
271 int computed = 0;
272 Vector3D<cosmoType> r;
273 cosmoType rsq;
274 cosmoType twoh, a, b;
275
276 for(int j = req->firstParticle; j <= req->lastParticle; ++j) {
277 if (particles[j].rung >= activeRung) {
278#ifdef CMK_VERSION_BLUEGENE
279 if (++forProgress > 200) {
280 forProgress = 0;
281#ifdef COSMO_EVENTS
282 traceUserEvents(networkProgressUE);
283#endif
284 CmiNetworkProgress();
285 }
286#endif
287 computed++;
288 r = offset + part->position - particles[j].position;
289 rsq = r.lengthSquared();
290 twoh = part->soft + particles[j].soft;
291 if(rsq != 0) {
292 SPLINE(rsq, twoh, a, b);
293 cosmoType idt2 = (particles[j].mass + part->mass)*b; // (timescale)^-2
294 // of interaction
295 particles[j].treeAcceleration += r * (b * part->mass);
296 particles[j].potential -= part->mass * a;
297 if(idt2 > particles[j].dtGrav)
298 particles[j].dtGrav = idt2;
299 }
300 }
301 }
302 return computed;
303}
304#else
305inline int partBucketForce(ExternalGravityParticle *part,
306 Tree::GenericTreeNode *req,
307 GravityParticle **activeParticles,
308 Vector3D<cosmoType> offset, int nActiveParts) {
309 Vector3D<SSEcosmoType> r;
310 SSEcosmoType rsq;
311 SSEcosmoType twoh, a, b;
312
313 for (int i=0; i<nActiveParts; i+=SSE_VECTOR_WIDTH) {
314#ifdef CMK_VERSION_BLUEGENE
315 if (++forProgress > 200) {
316 forProgress = 0;
317#ifdef COSMO_EVENTS
318 traceUserEvents(networkProgressUE);
319#endif
320 CmiNetworkProgress();
321 }
322#endif
323 Vector3D<SSEcosmoType>
324 packedPos(SSELoad(SSEcosmoType, activeParticles, i, ->position.x),
325 SSELoad(SSEcosmoType, activeParticles, i, ->position.y),
326 SSELoad(SSEcosmoType, activeParticles, i, ->position.z));
327 SSELoad(SSEcosmoType packedSoft, activeParticles, i, ->soft);
328
329 r = -packedPos + offset + part->position;
330 rsq = r.lengthSquared();
331 twoh = part->soft + packedSoft;
332 SSEcosmoType select = rsq > COSMO_CONST(0.0);
333 int compare = movemask(select);
334 if(compare) {
335 SPLINE(rsq, twoh, a, b);
336 if ((~compare) & cosmoMask) {
337 a = select & a;
338 b = select & b;
339 }
340 SSEcosmoType SSELoad(packedMass, activeParticles, i, ->mass);
341 SSEcosmoType SSELoad(packedDtGrav, activeParticles, i, ->dtGrav);
342 Vector3D<SSEcosmoType>
343 packedAcc(SSELoad(SSEcosmoType,
344 activeParticles, i, ->treeAcceleration.x),
345 SSELoad(SSEcosmoType,
346 activeParticles, i, ->treeAcceleration.y),
347 SSELoad(SSEcosmoType,
348 activeParticles, i, ->treeAcceleration.z));
349 SSEcosmoType SSELoad(packedPotential, activeParticles, i, ->potential);
350 SSEcosmoType idt2 = (packedMass + part->mass) * b;
351 idt2 = max(idt2, packedDtGrav);
352 packedAcc += r * (b * part->mass);
353 packedPotential -= part->mass * a;
354 SSEStore(packedAcc.x, activeParticles, i, ->treeAcceleration.x);
355 SSEStore(packedAcc.y, activeParticles, i, ->treeAcceleration.y);
356 SSEStore(packedAcc.z, activeParticles, i, ->treeAcceleration.z);
357 SSEStore(packedPotential, activeParticles, i, ->potential);
358 SSEStore(idt2, activeParticles, i, ->dtGrav);
359 }
360 }
361 return nActiveParts;
362}
363
364inline int partBucketForce(ExternalGravityParticle *part,
365 Tree::GenericTreeNode *req,
366 GravityParticle *particles,
367 Vector3D<cosmoType> offset, int activeRung) {
368 int nActiveParts = 0;
369 GravityParticle dummyPart = particles[0];
370 dummyPart.soft = 0.0;
371
372 GravityParticle **activeParticles =
373 (GravityParticle**)alloca((req->lastParticle - req->firstParticle
374 + 1 + FORCE_INPUT_LIST_PAD)
375 * sizeof(GravityParticle*));
376 for (int j=req->firstParticle; j <= req->lastParticle; ++j) {
377 if (particles[j].rung >= activeRung)
378 activeParticles[nActiveParts++] = &particles[j];
379 }
380
381 activeParticles[nActiveParts] = &dummyPart;
382#if SSE_VECTOR_WIDTH == 4
383 activeParticles[nActiveParts+1] = &dummyPart;
384 activeParticles[nActiveParts+2] = &dummyPart;
385#endif
386
387 int ret = partBucketForce(part, req, activeParticles, offset, nActiveParts);
388 return ret;
389}
390
391#endif
392
393
394//
395// Calculated forces on active particles in a bucket due to the
396// multipole of a TreeNode. Return number of multipoles evaluated.
397//
398#if !CMK_SSE
399inline
400int nodeBucketForce(Tree::GenericTreeNode *node,
401 Tree::GenericTreeNode *req,
402 GravityParticle *particles,
403 Vector3D<cosmoType> offset,
404 int activeRung)
405
406{
407 int computed = 0;
408 Vector3D<cosmoType> r;
409 cosmoType rsq;
410 MultipoleMoments &m = node->moments;
411
412 Vector3D<cosmoType> cm(m.cm + offset);
413
414#ifdef HEXADECAPOLE
415 if(openSoftening(node, req, offset)) {
416 ExternalGravityParticle tmpPart;
417 tmpPart.mass = m.totalMass;
418 tmpPart.soft = m.soft;
419 tmpPart.position = m.cm;
420 return partBucketForce(&tmpPart, req, particles, offset, activeRung);
421 }
422#endif
423 for(int j = req->firstParticle; j <= req->lastParticle; ++j) {
424 if (particles[j].rung >= activeRung) {
425 particles[j].interMass += m.totalMass;
426#ifdef CMK_VERSION_BLUEGENE
427 if (++forProgress > 200) {
428 forProgress = 0;
429#ifdef COSMO_EVENTS
430 traceUserEvents(networkProgressUE);
431#endif
432 CmiNetworkProgress();
433 }
434#endif
435 computed++;
436 r = Vector3D<cosmoType>(particles[j].position) - cm;
437 rsq = r.lengthSquared();
438 cosmoType dir = COSMO_CONST(1.0)/sqrt(rsq);
439#ifdef HEXADECAPOLE
440 cosmoType magai;
441 momEvalFmomrcm(&m.mom, m.getRadius(), dir, r.x, r.y, r.z,
442 &particles[j].potential,
443 &particles[j].treeAcceleration.x,
444 &particles[j].treeAcceleration.y,
445 &particles[j].treeAcceleration.z, &magai);
446 cosmoType idt2 = (particles[j].mass + m.totalMass)*dir*dir*dir;
447#else
448 cosmoType twoh, a, b, c, d;
449 twoh = CONVERT_TO_COSMO_TYPE(m.soft + particles[j].soft);
450 SPLINEQ(dir, rsq, twoh, a, b, c, d);
451 cosmoType qirx = CONVERT_TO_COSMO_TYPE m.xx*r.x
452 + CONVERT_TO_COSMO_TYPE m.xy*r.y + CONVERT_TO_COSMO_TYPE m.xz*r.z;
453 cosmoType qiry = CONVERT_TO_COSMO_TYPE m.xy*r.x
454 + CONVERT_TO_COSMO_TYPE m.yy*r.y + CONVERT_TO_COSMO_TYPE m.yz*r.z;
455 cosmoType qirz = CONVERT_TO_COSMO_TYPE m.xz*r.x
456 + CONVERT_TO_COSMO_TYPE m.yz*r.y + CONVERT_TO_COSMO_TYPE m.zz*r.z;
457 cosmoType qir = COSMO_CONST(0.5)*(qirx*r.x + qiry*r.y + qirz*r.z);
458 cosmoType tr = COSMO_CONST(0.5)*(CONVERT_TO_COSMO_TYPE m.xx
459 + CONVERT_TO_COSMO_TYPE m.yy
460 + CONVERT_TO_COSMO_TYPE m.zz);
461 cosmoType qir3 = b*CONVERT_TO_COSMO_TYPE m.totalMass + d*qir - c*tr;
462 particles[j].potential -= CONVERT_TO_COSMO_TYPE m.totalMass * a
463 + c*qir - b*tr;
464 particles[j].treeAcceleration.x -= qir3*r.x - c*qirx;
465 particles[j].treeAcceleration.y -= qir3*r.y - c*qiry;
466 particles[j].treeAcceleration.z -= qir3*r.z - c*qirz;
467 cosmoType idt2 = (CONVERT_TO_COSMO_TYPE particles[j].mass
468 + CONVERT_TO_COSMO_TYPE m.totalMass)*b;
469#endif
470 if(idt2 > particles[j].dtGrav)
471 particles[j].dtGrav = idt2;
472 }
473 }
474 return computed;
475}
476
477#elif CMK_SSE
478inline
479int nodeBucketForce(Tree::GenericTreeNode *node,
480 Tree::GenericTreeNode *req,
481 GravityParticle *particles,
482 Vector3D<cosmoType> offset,
483 int activeRung)
484
485{
486 Vector3D<SSEcosmoType> r;
487 SSEcosmoType rsq;
488 SSEcosmoType twoh;
489 SSEcosmoType a,b,c,d;
490 MultipoleMoments m = node->moments;
491 Vector3D<cosmoType> cm(m.cm + offset);
492 int nActiveParts = 0;
493 GravityParticle dummyPart = particles[0];
494 dummyPart.soft = 0.0;
495
496 GravityParticle **activeParticles =
497 (GravityParticle**)alloca((req->lastParticle - req->firstParticle
498 + 1 + FORCE_INPUT_LIST_PAD)
499 * sizeof(GravityParticle*));
500 for (int j=req->firstParticle; j <= req->lastParticle; ++j) {
501 if (particles[j].rung >= activeRung)
502 activeParticles[nActiveParts++] = &particles[j];
503 }
504
505 activeParticles[nActiveParts] = &dummyPart;
506#if SSE_VECTOR_WIDTH == 4
507 activeParticles[nActiveParts+1] = &dummyPart;
508 activeParticles[nActiveParts+2] = &dummyPart;
509#endif
510
511#ifdef HEXADECAPOLE
512 if(openSoftening(node, req, offset)) {
513 ExternalGravityParticle tmpPart;
514 tmpPart.mass = m.totalMass;
515 tmpPart.soft = m.soft;
516 tmpPart.position = m.cm;
517 int ret = partBucketForce(&tmpPart, req, activeParticles,
518 offset, nActiveParts);
519 return ret;
520 }
521#endif
522 for (int i=0; i<nActiveParts; i+=SSE_VECTOR_WIDTH) {
523#ifdef CMK_VERSION_BLUEGENE
524 if (++forProgress > 200) {
525 forProgress = 0;
526#ifdef COSMO_EVENTS
527 traceUserEvents(networkProgressUE);
528#endif
529 CmiNetworkProgress();
530 }
531#endif
532 Vector3D<SSEcosmoType>
533 packedPos(SSELoad(SSEcosmoType, activeParticles, i, ->position.x),
534 SSELoad(SSEcosmoType, activeParticles, i, ->position.y),
535 SSELoad(SSEcosmoType, activeParticles, i, ->position.z));
536 r = packedPos - cm;
537 rsq = r.lengthSquared();
538 SSEcosmoType dir = COSMO_CONST(1.0)/sqrt(rsq);
539 Vector3D<SSEcosmoType>
540 packedAcc(SSELoad(SSEcosmoType,
541 activeParticles, i, ->treeAcceleration.x),
542 SSELoad(SSEcosmoType,
543 activeParticles, i, ->treeAcceleration.y),
544 SSELoad(SSEcosmoType,
545 activeParticles, i, ->treeAcceleration.z));
546 SSEcosmoType SSELoad(packedPotential, activeParticles, i, ->potential);
547 SSEcosmoType SSELoad(packedMass, activeParticles, i, ->mass);
548 SSEcosmoType SSELoad(packedDtGrav, activeParticles, i, ->dtGrav);
549#ifdef HEXADECAPOLE
550 SSEcosmoType magai;
551 momEvalFmomrcm(&m.mom, m.getRadius(), dir, r.x, r.y, r.z,
552 &packedPotential,
553 &packedAcc.x,
554 &packedAcc.y,
555 &packedAcc.z, &magai);
556 SSEcosmoType idt2 = (packedMass + m.totalMass)*dir*dir*dir;
557#else
558 SSELoad(SSEcosmoType packedSoft, activeParticles, i, ->soft);
559 twoh = CONVERT_TO_COSMO_TYPE m.soft + packedSoft;
560 SPLINEQ(dir, rsq, twoh, a, b, c, d);
561 SSEcosmoType qirx = CONVERT_TO_COSMO_TYPE m.xx*r.x
562 + CONVERT_TO_COSMO_TYPE m.xy*r.y + CONVERT_TO_COSMO_TYPE m.xz*r.z;
563 SSEcosmoType qiry = CONVERT_TO_COSMO_TYPE m.xy*r.x
564 + CONVERT_TO_COSMO_TYPE m.yy*r.y + CONVERT_TO_COSMO_TYPE m.yz*r.z;
565 SSEcosmoType qirz = CONVERT_TO_COSMO_TYPE m.xz*r.x
566 + CONVERT_TO_COSMO_TYPE m.yz*r.y + CONVERT_TO_COSMO_TYPE m.zz*r.z;
567 SSEcosmoType qir = COSMO_CONST(0.5)*(qirx*r.x + qiry*r.y + qirz*r.z);
568 SSEcosmoType tr = COSMO_CONST(0.5)*(CONVERT_TO_COSMO_TYPE m.xx
569 + CONVERT_TO_COSMO_TYPE m.yy
570 + CONVERT_TO_COSMO_TYPE m.zz);
571 SSEcosmoType qir3 = b*CONVERT_TO_COSMO_TYPE m.totalMass + d*qir - c*tr;
572 packedPotential -= CONVERT_TO_COSMO_TYPE m.totalMass *a
573 + c * qir - b * tr;
574 packedAcc.x -= qir3*r.x - c*qirx;
575 packedAcc.y -= qir3*r.y - c*qiry;
576 packedAcc.z -= qir3*r.z - c*qirz;
577 SSEcosmoType idt2 = (packedMass + CONVERT_TO_COSMO_TYPE m.totalMass)*b;
578#endif
579 SSEStore(packedAcc.x, activeParticles, i, ->treeAcceleration.x);
580 SSEStore(packedAcc.y, activeParticles, i, ->treeAcceleration.y);
581 SSEStore(packedAcc.z, activeParticles, i, ->treeAcceleration.z);
582 SSEStore(packedPotential, activeParticles, i, ->potential);
583 idt2 = max(idt2, packedDtGrav);
584 SSEStore(idt2, activeParticles, i, ->dtGrav);
585 }
586 return nActiveParts;
587}
588#endif
589
596
597inline bool
598openCriterionBucket(Tree::GenericTreeNode *node,
599 Tree::GenericTreeNode *bucketNode,
600 Vector3D<cosmoType> offset // Offset of node
601 ) {
602
603#if COSMO_STATS > 0
604 node->used = true;
605#endif
606 // Always open node if this many particles or fewer.
607 const int nMinParticleNode = 6;
608 if(node->particleCount <= nMinParticleNode) {
609 return true;
610 }
611
612 // Note that some of this could be pre-calculated into an "opening radius"
613 cosmoType radius = TreeStuff::opening_geometry_factor * node->moments.getRadius() / theta;
614 if(radius < node->moments.getRadius())
615 radius = node->moments.getRadius();
616
617 Sphere<cosmoType> s(node->moments.cm + offset, radius);
618
619#ifdef HEXADECAPOLE
620 if(!Space::intersect(bucketNode->boundingBox, s)) {
621 // Well separated, now check softening
622 if(!openSoftening(node, bucketNode, offset)) {
623 return false; // passed both tests: will be a Hex interaction
624 }
625 else { // Open as monopole?
626 radius = TreeStuff::opening_geometry_factor*node->moments.getRadius()/thetaMono;
627 Sphere<cosmoType> sM(node->moments.cm + offset, radius);
628 return Space::intersect(bucketNode->boundingBox, sM);
629 }
630 }
631 return true;
632#else
633 return Space::intersect(bucketNode->boundingBox, s);
634#endif
635}
636
651
652inline int openCriterionNode(Tree::GenericTreeNode *node,
653 Tree::GenericTreeNode *myNode,
654 Vector3D<cosmoType> offset
655 ) {
656
657#if COSMO_STATS > 0
658 node->used = true;
659#endif
660 // Always open node if this many particles or fewer.
661 const int nMinParticleNode = 6;
662 if(node->particleCount <= nMinParticleNode) {
663 return 1;
664 }
665
666 // Note that some of this could be pre-calculated into an "opening radius"
667 cosmoType radius = TreeStuff::opening_geometry_factor * node->moments.getRadius() / theta;
668 if(radius < node->moments.getRadius())
669 radius = node->moments.getRadius();
670
671 Sphere<cosmoType> s(node->moments.cm + offset, radius);
672
673 if(myNode->getType()==Tree::Bucket || myNode->getType()==Tree::CachedBucket || myNode->getType()==Tree::NonLocalBucket){
674 if(Space::intersect(myNode->boundingBox, s))
675 return 1;
676 else
677#ifdef HEXADECAPOLE
678 {
679 // Well separated, now check softening
680 if(!openSoftening(node, myNode, offset)) {
681 return 0; // passed both tests: will be a Hex interaction
682 }
683 else { // Open as monopole?
684 radius = TreeStuff::opening_geometry_factor*node->moments.getRadius()/thetaMono;
685 Sphere<cosmoType> sM(node->moments.cm + offset, radius);
686 if(Space::intersect(myNode->boundingBox, sM))
687 return 1;
688 else
689 return 0;
690 }
691 }
692#else
693 return 0;
694#endif
695 }
696 else{
697 if(Space::intersect(myNode->boundingBox, s)){
698 if(Space::contained(myNode->boundingBox,s))
699 return 1;
700 else
701 return -1;
702 }
703 else
704#ifdef HEXADECAPOLE
705 {
706 // Well separated, now check softening
707 if(!openSoftening(node, myNode, offset)) {
708 return 0; // passed both tests: will be a Hex interaction
709 }
710 else { // Open as monopole?
711 radius = TreeStuff::opening_geometry_factor*node->moments.getRadius()/thetaMono;
712 Sphere<cosmoType> sM(node->moments.cm + offset, radius);
713 if(Space::intersect(myNode->boundingBox, sM))
714 return 1;
715 else
716 return 0;
717 }
718 }
719#else
720 return 0;
721#endif
722 }
723}
724
725#endif
TreeStuff::opening_geometry_factor
const double opening_geometry_factor
Definition TreeNode.h:35
ExternalGravityParticle
Information needed to calculate gravity.
Definition GravityParticle.h:33
GravityParticle
Fundamental type for a particle.
Definition GravityParticle.h:364
GravityParticle::dtGrav
cosmoType dtGrav
Definition GravityParticle.h:383
MultipoleMoments
A representation of a multipole expansion.
Definition MultipoleMoments.h:163
MultipoleMoments::cm
Vector3D< cosmoType > cm
The center of mass (zeroth order multipole)
Definition MultipoleMoments.h:174
MultipoleMoments::totalMass
cosmoType totalMass
The total mass represented by this expansion.
Definition MultipoleMoments.h:172
Tree::GenericTreeNode
Base class for tree nodes.
Definition GenericTreeNode.h:59
Tree::GenericTreeNode::particleCount
unsigned int particleCount
Total number of particles contained (across all chares)
Definition GenericTreeNode.h:115
Tree::GenericTreeNode::getType
NodeType getType() const
return Tree::NodeType of node
Definition GenericTreeNode.h:166
Tree::GenericTreeNode::lastParticle
int lastParticle
An index to the last particle contained by this node, myNumParticles+1 means outside the node.
Definition GenericTreeNode.h:108
Tree::GenericTreeNode::boundingBox
OrientedBox< cosmoType > boundingBox
The axis-aligned bounding box of this node.
Definition GenericTreeNode.h:98
Tree::GenericTreeNode::firstParticle
int firstParticle
An index for the first particle contained by this node, 0 means outside the node.
Definition GenericTreeNode.h:106
Tree::GenericTreeNode::moments
MultipoleMoments moments
The moments for the gravity computation.
Definition GenericTreeNode.h:94
Generated by doxygen 1.13.2