GADGET-4
gravtree.h
Go to the documentation of this file.
1/*******************************************************************************
2 * \copyright This file is part of the GADGET4 N-body/SPH code developed
3 * \copyright by Volker Springel. Copyright (C) 2014-2020 by Volker Springel
4 * \copyright (vspringel@mpa-garching.mpg.de) and all contributing authors.
5 *******************************************************************************/
6
12#ifndef FORCETREE_H
13#define FORCETREE_H
14
15#include "gadgetconfig.h"
16
17#ifndef NTAB
18#define NTAB 256 /* size of short-range look up */
19#endif
20
21#define TREE_DO_BASE_PM 1
22#define TREE_DO_HIGHRES_PM 2
23#define TREE_ACTIVE_CUTTOFF_BASE_PM 4
24#define TREE_ACTIVE_CUTTOFF_HIGHRES_PM 8
25
26#define TREE_MIN_WORKSTACK_SIZE 100000
27#define TREE_EXPECTED_CYCLES 80
28
29#define NODE_TYPE_LOCAL_PARTICLE 0
30#define NODE_TYPE_TREEPOINT_PARTICLE 1
31#define NODE_TYPE_FETCHED_PARTICLE 2
32#define NODE_TYPE_LOCAL_NODE 3
33#define NODE_TYPE_FETCHED_NODE 4
34
35#define NODE_USE 0
36#define NODE_OPEN 1
37#define NODE_DISCARD 2
38
39#define MAX_TREE_ALLOC_FACTOR 30.0
40
41#define TAKE_NSLOTS_IN_ONE_GO 32
42
43#include <string.h>
44
45#include "../data/simparticles.h"
46#include "../data/symtensors.h"
47#include "../domain/domain.h"
48#include "../mpi_utils/mpi_utils.h"
49#include "../tree/tree.h"
50
63struct gravnode : public basenode
64{
65 MyDouble mass;
66 vector<MyIntPosType> s;
67#if(MULTIPOLE_ORDER >= 3) || (MULTIPOLE_ORDER >= 2 && defined(EXTRAPOTTERM))
68 symtensor2<MyDouble> Q2Tensor;
69#endif
70#if(MULTIPOLE_ORDER >= 4) || (MULTIPOLE_ORDER >= 3 && defined(EXTRAPOTTERM))
71 symtensor3<MyDouble> Q3Tensor;
72#endif
73#if(MULTIPOLE_ORDER >= 5) || (MULTIPOLE_ORDER >= 4 && defined(EXTRAPOTTERM))
74 symtensor4<MyDouble> Q4Tensor;
75#endif
76#if(MULTIPOLE_ORDER >= 5 && defined(EXTRAPOTTERM))
77 symtensor5<MyDouble> Q5Tensor;
78#endif
79#ifdef FMM
80 float MinOldAcc;
81#endif
82#if(NSOFTCLASSES > 1)
83 unsigned char maxsofttype;
84 unsigned char minsofttype;
85#endif
86#if defined(PMGRID) && defined(PLACEHIGHRESREGION)
87 unsigned char overlap_flag : 2;
88#endif
89
90 inline int getSofteningClass(void)
91 {
92#if NSOFTCLASSES > 1
93 return maxsofttype;
94#else
95 return 0;
96#endif
97 }
98};
99
100struct gravpoint_data
101{
102 MyIntPosType IntPos[3];
103 MyDouble Mass;
104 float OldAcc;
105 int index;
106 int no;
107 unsigned char Type;
108#if NSOFTCLASSES > 1
109 unsigned char SofteningClass : 7;
110#endif
111#ifndef HIERARCHICAL_GRAVITY
112 unsigned char ActiveFlag : 1; /* we don't need this for hierarchical gravity as then the particles are always active */
113#endif
114#if defined(PMGRID) && defined(PLACEHIGHRESREGION)
115 unsigned char InsideOutsideFlag : 1;
116#endif
117
118 inline unsigned char getSofteningClass(void)
119 {
120#if NSOFTCLASSES > 1
121 return SofteningClass;
122#else
123 return 0;
124#endif
125 }
126};
127
128struct foreign_gravpoint_data
129{
130 MyIntPosType IntPos[3];
131 MyDouble Mass;
132 int Nextnode;
133 unsigned char Nextnode_shmrank;
134 unsigned char Type;
135 float OldAcc;
136#if NSOFTCLASSES > 1
137 unsigned char SofteningClass;
138#endif
139#if defined(PMGRID) && defined(PLACEHIGHRESREGION)
140 unsigned char InsideOutsideFlag : 1;
141#endif
142
143 inline unsigned char getSofteningClass(void)
144 {
145#if NSOFTCLASSES > 1
146 return SofteningClass;
147#else
148 return 0;
149#endif
150 }
151};
152
153#define HIGHEST_NEEDEDORDER_EWALD_DPHI 1
154
155#if(MULTIPOLE_ORDER >= 3) || (MULTIPOLE_ORDER >= 2 && defined(EXTRAPOTTERM) || (MULTIPOLE_ORDER >= 2 && defined(FMM)))
156#undef HIGHEST_NEEDEDORDER_EWALD_DPHI
157#define HIGHEST_NEEDEDORDER_EWALD_DPHI 2
158#endif
159
160#if(MULTIPOLE_ORDER >= 4) || (MULTIPOLE_ORDER >= 3 && defined(EXTRAPOTTERM) || (MULTIPOLE_ORDER >= 3 && defined(FMM)))
161#undef HIGHEST_NEEDEDORDER_EWALD_DPHI
162#define HIGHEST_NEEDEDORDER_EWALD_DPHI 3
163#endif
164
165#if(MULTIPOLE_ORDER >= 5) || (MULTIPOLE_ORDER >= 4 && defined(EXTRAPOTTERM) || (MULTIPOLE_ORDER >= 4 && defined(FMM)))
166#undef HIGHEST_NEEDEDORDER_EWALD_DPHI
167#define HIGHEST_NEEDEDORDER_EWALD_DPHI 4
168#endif
169
170#if((MULTIPOLE_ORDER >= 5 && defined(EXTRAPOTTERM) && defined(EVALPOTENTIAL)) || (MULTIPOLE_ORDER >= 5 && defined(FMM)) || \
171 defined(EWALD_TEST))
172#undef HIGHEST_NEEDEDORDER_EWALD_DPHI
173#define HIGHEST_NEEDEDORDER_EWALD_DPHI 5
174#endif
175
176#ifdef EXTRA_HIGH_EWALD_ACCURACY
177#define EWALD_TAYLOR_ORDER 3
178#else
179#define EWALD_TAYLOR_ORDER 2
180#endif
181
182/* variables for Ewald correction lookup table */
183struct ewald_data
184{
185 MyReal D0phi;
186 vector<MyReal> D1phi;
187 symtensor2<MyReal> D2phi;
188 symtensor3<MyReal> D3phi;
189#if(HIGHEST_NEEDEDORDER_EWALD_DPHI + EWALD_TAYLOR_ORDER) >= 4
190 symtensor4<MyReal> D4phi;
191#endif
192#if(HIGHEST_NEEDEDORDER_EWALD_DPHI + EWALD_TAYLOR_ORDER) >= 5
193 symtensor5<MyReal> D5phi;
194#endif
195#if(HIGHEST_NEEDEDORDER_EWALD_DPHI + EWALD_TAYLOR_ORDER) >= 6
196 symtensor6<MyReal> D6phi;
197#endif
198#if(HIGHEST_NEEDEDORDER_EWALD_DPHI + EWALD_TAYLOR_ORDER) >= 7
199 symtensor7<MyReal> D7phi;
200#endif
201};
202
203template <typename partset> /* partset will either be 'simparticles' or 'lightconeparticles' */
204class gravtree : public tree<gravnode, partset, gravpoint_data, foreign_gravpoint_data>
205{
206 public:
207 typedef tree<gravnode, partset, gravpoint_data, foreign_gravpoint_data> basetree;
208 using basetree::Buildtime;
209 using basetree::D; // this avoids that we have to use a "this->" when accessing these variables from the base class
210 using basetree::EndOfForeignNodes;
211 using basetree::EndOfTreePoints;
212 using basetree::Father;
213 using basetree::FirstNonTopLevelNode;
214 using basetree::Foreign_Nodes;
215 using basetree::get_nodep;
216 using basetree::ImportedNodeOffset;
217 using basetree::IndexList;
218 using basetree::MaxForeignNodes;
219 using basetree::MaxNodes;
220 using basetree::MaxPart;
221 using basetree::Nextnode;
222 using basetree::Ninsert;
223 using basetree::NodeIndex;
224 using basetree::Nodes;
225 using basetree::NumForeignNodes;
226 using basetree::NumForeignPoints;
227 using basetree::NumNodes;
228 using basetree::NumPartExported;
229 using basetree::NumPartImported;
230 using basetree::Points;
231 using basetree::Recv_count;
232 using basetree::Recv_offset;
233 using basetree::ResultIndexList;
234 using basetree::Send_count;
235 using basetree::Send_offset;
236 using basetree::TopNodes;
237 using basetree::Tp;
238 using basetree::TreeForeign_Nodes_offsets;
239 using basetree::TreeForeign_Points_offsets;
240 using basetree::TreeNextnode_offsets;
241 using basetree::TreeNodes_offsets;
242 using basetree::TreeP_offsets;
243 using basetree::TreePoints_offsets;
244 using basetree::TreePS_offsets;
245 using basetree::TreeSharedMem_ThisTask;
246 using basetree::TreeSharedMemBaseAddr;
247 using basetree::TreeSharedMemComm;
248 using basetree::TreeSphP_offsets;
249
250 struct mesh_factors
251 {
252 double rcut2;
253 double asmthfac;
254
255 double asmthinv1;
256 double asmthinv2;
257#if(MULTIPOLE_ORDER >= 2)
258 double asmthinv3;
259#endif
260#if(MULTIPOLE_ORDER >= 3)
261 double asmthinv4;
262#endif
263#if(MULTIPOLE_ORDER >= 4)
264 double asmthinv5;
265#endif
266#if(MULTIPOLE_ORDER >= 5)
267 double asmthinv6;
268#endif
269 MyIntPosType intrcut[3];
270 };
271 mesh_factors mf[2];
272
273 struct resultsactiveimported_data
274 {
275 vector<MyFloat> GravAccel;
276#ifdef EVALPOTENTIAL
277 MyFloat Potential;
278#endif
279 int GravCost;
280 int index;
281 };
282
283 struct index_data
284 {
285 int p;
286 int subnode;
287 };
288
291 char MeasureCostFlag;
292
293 resultsactiveimported_data *ResultsActiveImported;
294
295 ewald_data PotTaylor;
296 int num_layers = 0;
297
298 private:
307 private:
308 /* private member functions */
309
310 void update_node_recursive(int no, int sib, int mode) override;
311 void exchange_topleafdata(void) override;
312 void fill_in_export_points(gravpoint_data *exp_point, int i, int no) override;
313 void report_log_message(void) override;
314
315 public:
316 void set_softenings(void);
317
318#ifdef ALLOW_DIRECT_SUMMATION
319 void gravity_direct(partset *Pptr, domain<partset> *Dptr, int timebin);
320#endif
321
322 void gravity_exchange_forces(void);
323
325 public:
326#ifdef PMGRID
327 void short_range_init(void);
328 void set_mesh_factors(void);
329#endif
330
331 public:
332#ifdef PMGRID
338 struct
339 {
340 float fac0;
341 float fac1;
342#if(MULTIPOLE_ORDER >= 2)
343 float fac2;
344#endif
345#if(MULTIPOLE_ORDER >= 3)
346 float fac3;
347#endif
348#if(MULTIPOLE_ORDER >= 4)
349 float fac4;
350#endif
351#if(MULTIPOLE_ORDER >= 5)
352 float fac5;
353#endif
354 } shortrange_factors[NTAB + 1];
355#endif
356
357 public:
358#if defined(PMGRID)
359 char DoPM;
360#endif
361 char DoEwald;
362
363 struct gfactors
364 {
365 MyReal rinv2;
366#if(MULTIPOLE_ORDER >= 2)
367 MyReal rinv3;
368#endif
369
370 MyReal fac0; // g
371 MyReal fac1; // g'
372#if(MULTIPOLE_ORDER >= 2)
373 MyReal fac2; // (g'' - g'/r) = (g'' - g'/r)
374#endif
375#if(MULTIPOLE_ORDER >= 3)
376 MyReal fac3; // (g''' - 3*g''/r + 3 * g'/r^2) = (g''' - (3/r) * fac2)
377#endif
378#if(MULTIPOLE_ORDER >= 4)
379 MyReal fac4; // (g'''' - 6*g'''/r + 15*g''/r^2 - 15*g'/r^3) = (g'''' - (6/r) * fac3 - (3/r^2) * fac2)
380#endif
381#if(MULTIPOLE_ORDER >= 5)
382 MyReal fac5; // (g''''' - 10*g''''/r + 45*g'''/r^2 - 105*g''/r^3 + 105 g'/r^4) = (g''''' - (10/r) * fac4 - (15/r^2) * fac3)
383#endif
384
385 gfactors() /* constructor, initialize the factors to zero */
386 {
387 fac0 = 0;
388 fac1 = 0;
389#if(MULTIPOLE_ORDER >= 2)
390 fac2 = 0;
391#endif
392#if(MULTIPOLE_ORDER >= 3)
393 fac3 = 0;
394#endif
395#if(MULTIPOLE_ORDER >= 4)
396 fac4 = 0;
397#endif
398#if(MULTIPOLE_ORDER >= 5)
399 fac5 = 0;
400#endif
401 }
402 };
403
404 template <typename T>
405 inline void get_gfactors_multipole(gfactors &res, const T r, const T h_max, const T rinv)
406 {
407 res.rinv2 = rinv * rinv;
408
409#if(MULTIPOLE_ORDER >= 2)
410 res.rinv3 = res.rinv2 * rinv;
411#endif
412
413 if(r >= h_max)
414 {
415#ifdef EVALPOTENTIAL
416 res.fac0 += rinv;
417#endif
418 res.fac1 -= res.rinv2;
419#if(MULTIPOLE_ORDER >= 2)
420 res.fac2 += 3 * res.rinv3;
421#endif
422#if(MULTIPOLE_ORDER >= 3)
423 res.fac3 -= 15 * res.rinv3 * rinv;
424#endif
425#if(MULTIPOLE_ORDER >= 4)
426 res.fac4 += 105 * res.rinv3 * res.rinv2;
427#endif
428#if(MULTIPOLE_ORDER >= 5)
429 res.fac5 -= 945 * res.rinv3 * res.rinv3;
430#endif
431 }
432 else
433 {
434 T h_inv = 1 / h_max;
435 T h2_inv = h_inv * h_inv;
436 T u = r * h_inv;
437 T fac1;
438#if(MULTIPOLE_ORDER >= 2)
439 T h3_inv = h_inv * h2_inv;
440 T gpp;
441#endif
442#if(MULTIPOLE_ORDER >= 3)
443 T gppp;
444#endif
445#if(MULTIPOLE_ORDER >= 4)
446 T gpppp;
447#endif
448#if(MULTIPOLE_ORDER >= 5)
449 T gppppp;
450#endif
451 if(u < static_cast<T>(0.5))
452 {
453 T u2 = u * u;
454#ifdef EVALPOTENTIAL
455 res.fac0 -= h_inv * (static_cast<T>(SOFTFAC4) +
456 u2 * (static_cast<T>(SOFTFAC5) + u2 * (static_cast<T>(SOFTFAC6) * u + static_cast<T>(SOFTFAC7))));
457#endif
458 fac1 = -h2_inv * u * (static_cast<T>(SOFTFAC1) + u2 * (static_cast<T>(SOFTFAC2) * u + static_cast<T>(SOFTFAC3)));
459 res.fac1 += fac1;
460#if(MULTIPOLE_ORDER >= 2)
461 gpp = -h3_inv * (static_cast<T>(SOFTFAC30) + (static_cast<T>(SOFTFAC31) + static_cast<T>(SOFTFAC32) * u) * u2);
462#endif
463#if(MULTIPOLE_ORDER >= 3)
464 gppp = -h3_inv * h_inv * (static_cast<T>(SOFTFAC33) * u + static_cast<T>(SOFTFAC34) * u2);
465#endif
466#if(MULTIPOLE_ORDER >= 4)
467 gpppp = -h3_inv * h2_inv * (static_cast<T>(SOFTFAC33) + static_cast<T>(SOFTFAC35) * u);
468#endif
469#if(MULTIPOLE_ORDER >= 5)
470 gppppp = -h3_inv * h3_inv * static_cast<T>(SOFTFAC35);
471#endif
472 }
473 else
474 {
475 T u2 = u * u;
476 T u3 = u2 * u;
477 T u3inv = 1 / u3;
478#ifdef EVALPOTENTIAL
479 res.fac0 -=
480 h_inv * (static_cast<T>(SOFTFAC13) + static_cast<T>(SOFTFAC14) / u +
481 u2 * (static_cast<T>(SOFTFAC1) +
482 u * (static_cast<T>(SOFTFAC15) + u * (static_cast<T>(SOFTFAC16) + static_cast<T>(SOFTFAC17) * u))));
483#endif
484 fac1 = -h2_inv * u *
485 (static_cast<T>(SOFTFAC8) + static_cast<T>(SOFTFAC9) * u + static_cast<T>(SOFTFAC10) * u2 +
486 static_cast<T>(SOFTFAC11) * u3 + static_cast<T>(SOFTFAC12) * u3inv);
487 res.fac1 += fac1;
488#if(MULTIPOLE_ORDER >= 2)
489 gpp = -h3_inv * (static_cast<T>(SOFTFAC40) + static_cast<T>(SOFTFAC41) / (u * u2) + static_cast<T>(SOFTFAC42) * u +
490 static_cast<T>(SOFTFAC43) * u2 + static_cast<T>(SOFTFAC44) * u2 * u);
491#endif
492#if(MULTIPOLE_ORDER >= 3)
493 gppp = -h3_inv * h_inv *
494 (static_cast<T>(SOFTFAC45) + static_cast<T>(SOFTFAC46) / (u2 * u2) + static_cast<T>(SOFTFAC47) * u +
495 static_cast<T>(SOFTFAC48) * u2);
496#endif
497#if(MULTIPOLE_ORDER >= 4)
498 gpppp =
499 -h3_inv * h2_inv * (static_cast<T>(SOFTFAC49) / (u2 * u3) + static_cast<T>(SOFTFAC47) + static_cast<T>(SOFTFAC50) * u);
500#endif
501#if(MULTIPOLE_ORDER >= 5)
502 gppppp = -h3_inv * h3_inv * (static_cast<T>(SOFTFAC51) / (u3 * u3) + static_cast<T>(SOFTFAC50));
503#endif
504 }
505#if(MULTIPOLE_ORDER >= 2)
506 T fac2 = (gpp - rinv * fac1);
507 res.fac2 += fac2;
508#endif
509#if(MULTIPOLE_ORDER >= 3)
510 T fac3 = (gppp - 3 * rinv * fac2);
511 res.fac3 += fac3;
512#endif
513#if(MULTIPOLE_ORDER >= 4)
514 T fac4 = (gpppp - 6 * rinv * fac3 - 3 * rinv * rinv * fac2);
515 res.fac4 += fac4;
516#endif
517#if(MULTIPOLE_ORDER >= 5)
518 T fac5 = (gppppp - 10 * rinv * fac4 - 15 * rinv * rinv * fac3);
519 res.fac5 += fac5;
520#endif
521 }
522 }
523
524 template <typename T>
525 inline void get_gfactors_monopole(gfactors &res, const T r, const T h_max, const T rinv)
526 {
527 if(r > h_max)
528 {
529 res.fac1 -= rinv * rinv;
530#ifdef EVALPOTENTIAL
531 res.fac0 += rinv;
532#endif
533 }
534 else
535 {
536 T h_inv = 1 / h_max;
537 T h2_inv = h_inv * h_inv;
538 T u = r * h_inv;
539
540 if(u < 0.5f)
541 {
542 T u2 = u * u;
543 res.fac1 -= h2_inv * u * (static_cast<T>(SOFTFAC1) + u2 * (static_cast<T>(SOFTFAC2) * u + static_cast<T>(SOFTFAC3)));
544#ifdef EVALPOTENTIAL
545 res.fac0 -= h_inv * (static_cast<T>(SOFTFAC4) +
546 u2 * (static_cast<T>(SOFTFAC5) + u2 * (static_cast<T>(SOFTFAC6) * u + static_cast<T>(SOFTFAC7))));
547#endif
548 }
549 else
550 {
551 T u2 = u * u;
552 T u3 = u2 * u;
553 res.fac1 -= h2_inv * u *
554 (static_cast<T>(SOFTFAC8) + static_cast<T>(SOFTFAC9) * u + static_cast<T>(SOFTFAC10) * u2 +
555 static_cast<T>(SOFTFAC11) * u3 + static_cast<T>(SOFTFAC12) / u3);
556#ifdef EVALPOTENTIAL
557 res.fac0 -=
558 h_inv * (static_cast<T>(SOFTFAC13) + static_cast<T>(SOFTFAC14) / u +
559 u2 * (static_cast<T>(SOFTFAC1) +
560 u * (static_cast<T>(SOFTFAC15) + u * (static_cast<T>(SOFTFAC16) + static_cast<T>(SOFTFAC17) * u))));
561#endif
562 }
563 }
564 }
565
566 template <typename T>
567 inline void get_gfactors_potential(gfactors &res, const T r, const T hmax, const T rinv)
568 {
569 if(r >= hmax)
570 {
571 res.fac0 = rinv;
572#if(MULTIPOLE_ORDER >= 3)
573 T rinv3 = rinv * rinv * rinv;
574 res.fac1 = -rinv * rinv;
575 res.fac2 = 3 * rinv3;
576#endif
577 }
578 else
579 {
580 T h_inv = 1 / hmax;
581#if(MULTIPOLE_ORDER >= 3)
582 T h2_inv = h_inv * h_inv;
583 T h3_inv = h2_inv * h_inv;
584#endif
585 T u = r * h_inv;
586#if(MULTIPOLE_ORDER >= 3)
587 T fac1;
588#endif
589 if(u < 0.5)
590 {
591 T u2 = u * u;
592#if(MULTIPOLE_ORDER >= 3)
593 fac1 = -h2_inv * u * (static_cast<T>(SOFTFAC1) + u2 * (static_cast<T>(SOFTFAC2) * u + static_cast<T>(SOFTFAC3)));
594 res.fac1 += fac1;
595#endif
596 res.fac0 = -h_inv * (static_cast<T>(SOFTFAC4) +
597 u2 * (static_cast<T>(SOFTFAC5) + u2 * (static_cast<T>(SOFTFAC6) * u + static_cast<T>(SOFTFAC7))));
598 }
599 else
600 {
601 T u2 = u * u;
602#if(MULTIPOLE_ORDER >= 3)
603 T u3 = u2 * u;
604 fac1 = -h2_inv * u *
605 (static_cast<T>(SOFTFAC8) + static_cast<T>(SOFTFAC9) * u + static_cast<T>(SOFTFAC10) * u2 +
606 static_cast<T>(SOFTFAC11) * u3 + static_cast<T>(SOFTFAC12) / u3);
607 res.fac1 += fac1;
608#endif
609 res.fac0 =
610 -h_inv * (static_cast<T>(SOFTFAC13) + static_cast<T>(SOFTFAC14) / u +
611 u2 * (static_cast<T>(SOFTFAC1) +
612 u * (static_cast<T>(SOFTFAC15) + u * (static_cast<T>(SOFTFAC16) + static_cast<T>(SOFTFAC17) * u))));
613 }
614
615#if(MULTIPOLE_ORDER >= 3)
616 T gpp;
617 T u2 = u * u;
618
619 if(u < 0.5)
620 gpp = -h3_inv * (static_cast<T>(SOFTFAC30) + (static_cast<T>(SOFTFAC31) + static_cast<T>(SOFTFAC32) * u) * u2);
621 else
622 gpp = -h3_inv * (static_cast<T>(SOFTFAC40) + static_cast<T>(SOFTFAC41) / (u * u2) + static_cast<T>(SOFTFAC42) * u +
623 static_cast<T>(SOFTFAC43) * u2 + static_cast<T>(SOFTFAC44) * u2 * u);
624
625 res.fac2 = (gpp - rinv * fac1);
626#endif
627 }
628 }
629
630#ifdef PMGRID
631 inline bool modify_gfactors_pm_multipole(gfactors &res, const double r, const double rinv, const mesh_factors *mfp)
632 {
633 double tabentry = mfp->asmthfac * r;
634 int tabindex = (int)tabentry;
635
636 if(tabindex < NTAB)
637 {
638 double w1 = tabentry - tabindex;
639 double w0 = 1 - w1;
640
641#ifdef EVALPOTENTIAL
642 res.fac0 += mfp->asmthinv1 * (w0 * shortrange_factors[tabindex].fac0 + w1 * shortrange_factors[tabindex + 1].fac0);
643#endif
644 res.fac1 += mfp->asmthinv2 * (w0 * shortrange_factors[tabindex].fac1 + w1 * shortrange_factors[tabindex + 1].fac1);
645
646#if(MULTIPOLE_ORDER >= 2)
647 res.fac2 += mfp->asmthinv3 * (w0 * shortrange_factors[tabindex].fac2 + w1 * shortrange_factors[tabindex + 1].fac2);
648#endif
649#if(MULTIPOLE_ORDER >= 3)
650 res.fac3 += mfp->asmthinv4 * (w0 * shortrange_factors[tabindex].fac3 + w1 * shortrange_factors[tabindex + 1].fac3);
651#endif
652#if(MULTIPOLE_ORDER >= 4)
653 res.fac4 += mfp->asmthinv5 * (w0 * shortrange_factors[tabindex].fac4 + w1 * shortrange_factors[tabindex + 1].fac4);
654#endif
655#if(MULTIPOLE_ORDER >= 5)
656 res.fac5 += mfp->asmthinv6 * (w0 * shortrange_factors[tabindex].fac5 + w1 * shortrange_factors[tabindex + 1].fac5);
657#endif
658 return false;
659 }
660 else
661 return true;
662 }
663
664 inline bool modify_gfactors_pm_monopole(gfactors &res, const double r, const double rinv, const mesh_factors *mfp)
665 {
666 double tabentry = mfp->asmthfac * r;
667 int tabindex = (int)tabentry;
668
669 if(tabindex < NTAB)
670 {
671 double w1 = tabentry - tabindex;
672 double w0 = 1 - w1;
673
674#ifdef EVALPOTENTIAL
675 res.fac0 += mfp->asmthinv1 * (w0 * shortrange_factors[tabindex].fac0 + w1 * shortrange_factors[tabindex + 1].fac0);
676#endif
677 res.fac1 += mfp->asmthinv2 * (w0 * shortrange_factors[tabindex].fac1 + w1 * shortrange_factors[tabindex + 1].fac1);
678
679 return false;
680 }
681 else
682 return true;
683 }
684
685#endif
686};
687
688#endif
domain
Definition: domain.h:31
gravtree::gravity_exchange_forces
void gravity_exchange_forces(void)
Definition: gravtree.cc:134
gravtree::PotTaylor
ewald_data PotTaylor
Definition: gravtree.h:295
gravtree::basetree
tree< gravnode, partset, gravpoint_data, foreign_gravpoint_data > basetree
Definition: gravtree.h:207
gravtree::mf
mesh_factors mf[2]
Definition: gravtree.h:271
gravtree::num_layers
int num_layers
Definition: gravtree.h:296
gravtree::ResultsActiveImported
resultsactiveimported_data * ResultsActiveImported
Definition: gravtree.h:293
gravtree::DoEwald
char DoEwald
Definition: gravtree.h:361
gravtree::set_softenings
void set_softenings(void)
This function sets the (comoving) softening length of all particle types in the table All....
Definition: gravtree_build.cc:584
gravtree::get_gfactors_multipole
void get_gfactors_multipole(gfactors &res, const T r, const T h_max, const T rinv)
Definition: gravtree.h:405
gravtree::MeasureCostFlag
char MeasureCostFlag
Definition: gravtree.h:291
gravtree::get_gfactors_potential
void get_gfactors_potential(gfactors &res, const T r, const T hmax, const T rinv)
Definition: gravtree.h:567
gravtree::get_gfactors_monopole
void get_gfactors_monopole(gfactors &res, const T r, const T h_max, const T rinv)
Definition: gravtree.h:525
tree
Definition: tree.h:91
SOFTFAC3
#define SOFTFAC3
Definition: constants.h:391
SOFTFAC42
#define SOFTFAC42
Definition: constants.h:423
SOFTFAC13
#define SOFTFAC13
Definition: constants.h:401
SOFTFAC31
#define SOFTFAC31
Definition: constants.h:415
SOFTFAC48
#define SOFTFAC48
Definition: constants.h:429
SOFTFAC16
#define SOFTFAC16
Definition: constants.h:404
SOFTFAC35
#define SOFTFAC35
Definition: constants.h:419
SOFTFAC4
#define SOFTFAC4
Definition: constants.h:392
SOFTFAC14
#define SOFTFAC14
Definition: constants.h:402
SOFTFAC10
#define SOFTFAC10
Definition: constants.h:398
SOFTFAC11
#define SOFTFAC11
Definition: constants.h:399
SOFTFAC43
#define SOFTFAC43
Definition: constants.h:424
SOFTFAC46
#define SOFTFAC46
Definition: constants.h:427
SOFTFAC8
#define SOFTFAC8
Definition: constants.h:396
SOFTFAC2
#define SOFTFAC2
Definition: constants.h:390
SOFTFAC15
#define SOFTFAC15
Definition: constants.h:403
SOFTFAC32
#define SOFTFAC32
Definition: constants.h:416
SOFTFAC41
#define SOFTFAC41
Definition: constants.h:422
SOFTFAC1
#define SOFTFAC1
Definition: constants.h:389
SOFTFAC51
#define SOFTFAC51
Definition: constants.h:432
SOFTFAC40
#define SOFTFAC40
Definition: constants.h:421
SOFTFAC47
#define SOFTFAC47
Definition: constants.h:428
SOFTFAC33
#define SOFTFAC33
Definition: constants.h:417
SOFTFAC50
#define SOFTFAC50
Definition: constants.h:431
SOFTFAC45
#define SOFTFAC45
Definition: constants.h:426
SOFTFAC49
#define SOFTFAC49
Definition: constants.h:430
SOFTFAC44
#define SOFTFAC44
Definition: constants.h:425
SOFTFAC34
#define SOFTFAC34
Definition: constants.h:418
SOFTFAC30
#define SOFTFAC30
Definition: constants.h:414
SOFTFAC12
#define SOFTFAC12
Definition: constants.h:400
SOFTFAC9
#define SOFTFAC9
Definition: constants.h:397
SOFTFAC5
#define SOFTFAC5
Definition: constants.h:393
SOFTFAC6
#define SOFTFAC6
Definition: constants.h:394
SOFTFAC17
#define SOFTFAC17
Definition: constants.h:405
SOFTFAC7
#define SOFTFAC7
Definition: constants.h:395
MyDouble
float MyDouble
Definition: dtypes.h:87
MyFloat
float MyFloat
Definition: dtypes.h:86
MyIntPosType
uint32_t MyIntPosType
Definition: dtypes.h:35
MyReal
double MyReal
Definition: dtypes.h:82
NTAB
#define NTAB
Definition: gravtree.h:18
basenode
Definition: tree.h:48
ewald_data
Definition: gravtree.h:184
ewald_data::D1phi
vector< MyReal > D1phi
Definition: gravtree.h:186
ewald_data::D3phi
symtensor3< MyReal > D3phi
Definition: gravtree.h:188
ewald_data::D0phi
MyReal D0phi
Definition: gravtree.h:185
ewald_data::D2phi
symtensor2< MyReal > D2phi
Definition: gravtree.h:187
foreign_gravpoint_data
Definition: gravtree.h:129
foreign_gravpoint_data::getSofteningClass
unsigned char getSofteningClass(void)
Definition: gravtree.h:143
foreign_gravpoint_data::Mass
MyDouble Mass
Definition: gravtree.h:131
foreign_gravpoint_data::Type
unsigned char Type
Definition: gravtree.h:134
foreign_gravpoint_data::SofteningClass
unsigned char SofteningClass
Definition: gravtree.h:137
foreign_gravpoint_data::Nextnode_shmrank
unsigned char Nextnode_shmrank
Definition: gravtree.h:133
foreign_gravpoint_data::OldAcc
float OldAcc
Definition: gravtree.h:135
foreign_gravpoint_data::Nextnode
int Nextnode
Definition: gravtree.h:132
foreign_gravpoint_data::IntPos
MyIntPosType IntPos[3]
Definition: gravtree.h:130
gravnode::mass
MyDouble mass
Definition: gravtree.h:65
gravnode::s
vector< MyIntPosType > s
Definition: gravtree.h:66
gravnode::getSofteningClass
int getSofteningClass(void)
Definition: gravtree.h:90
gravnode::maxsofttype
unsigned char maxsofttype
Definition: gravtree.h:83
gravnode::minsofttype
unsigned char minsofttype
Definition: gravtree.h:84
gravpoint_data
Definition: gravtree.h:101
gravpoint_data::no
int no
Definition: gravtree.h:106
gravpoint_data::getSofteningClass
unsigned char getSofteningClass(void)
Definition: gravtree.h:118
gravpoint_data::index
int index
Definition: gravtree.h:105
gravpoint_data::ActiveFlag
unsigned char ActiveFlag
Definition: gravtree.h:112
gravpoint_data::Mass
MyDouble Mass
Definition: gravtree.h:103
gravpoint_data::Type
unsigned char Type
Definition: gravtree.h:107
gravpoint_data::SofteningClass
unsigned char SofteningClass
Definition: gravtree.h:109
gravpoint_data::OldAcc
float OldAcc
Definition: gravtree.h:104
gravpoint_data::IntPos
MyIntPosType IntPos[3]
Definition: gravtree.h:102
gravtree::index_data
Definition: gravtree.h:284
gravtree::index_data::subnode
int subnode
Definition: gravtree.h:286
gravtree::index_data::p
int p
Definition: gravtree.h:285
gravtree::mesh_factors
Definition: gravtree.h:251
gravtree::mesh_factors::rcut2
double rcut2
Definition: gravtree.h:252
gravtree::mesh_factors::asmthinv2
double asmthinv2
Definition: gravtree.h:256
gravtree::mesh_factors::asmthinv1
double asmthinv1
Definition: gravtree.h:255
gravtree::mesh_factors::intrcut
MyIntPosType intrcut[3]
Definition: gravtree.h:269
gravtree::mesh_factors::asmthfac
double asmthfac
Definition: gravtree.h:253
gravtree::resultsactiveimported_data
Definition: gravtree.h:274
gravtree::resultsactiveimported_data::index
int index
Definition: gravtree.h:280
gravtree::resultsactiveimported_data::GravCost
int GravCost
Definition: gravtree.h:279
gravtree::resultsactiveimported_data::GravAccel
vector< MyFloat > GravAccel
Definition: gravtree.h:275