GADGET-4
gravtree_build.cc
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#include "gadgetconfig.h"
13
14#include <math.h>
15#include <mpi.h>
16#include <stdio.h>
17#include <stdlib.h>
18#include <string.h>
19
20#include "../data/allvars.h"
21#include "../data/dtypes.h"
22#include "../data/intposconvert.h"
23#include "../data/mymalloc.h"
24#include "../domain/domain.h"
25#include "../gravtree/gravtree.h"
26#include "../io/io.h"
27#include "../logs/logs.h"
28#include "../logs/timer.h"
29#include "../main/simulation.h"
30#include "../mpi_utils/mpi_utils.h"
31#include "../ngbtree/ngbtree.h"
32#include "../pm/pm.h"
33#include "../sort/cxxsort.h"
34#include "../sort/peano.h"
35#include "../system/system.h"
36#include "../time_integration/timestep.h"
37
52template <typename partset>
53void gravtree<partset>::report_log_message(void)
54{
55 TIMER_START(CPU_LOGS);
56
57 int max_imported;
58 long long tot_imported;
59 sumup_large_ints(1, &NumPartImported, &tot_imported, D->Communicator);
60 MPI_Reduce(&NumPartImported, &max_imported, 1, MPI_INT, MPI_MAX, 0, D->Communicator);
61
62 MyReal numnodes = NumNodes, tot_numnodes;
63 MPI_Reduce(&numnodes, &tot_numnodes, 1, MPI_DOUBLE, MPI_SUM, 0, D->Communicator);
64
65 if(Ninsert == Tp->NumPart)
66 D->mpi_printf(
67 "GRAVTREE: Tree construction done. took %g sec <numnodes>=%g NTopnodes=%d NTopleaves=%d tree-build-scalability=%g\n",
68 Buildtime, tot_numnodes / D->NTask, D->NTopnodes, D->NTopleaves,
69 ((double)((tot_numnodes - D->NTask * ((double)D->NTopnodes)) + D->NTopnodes)) / (tot_numnodes > 0 ? tot_numnodes : 1));
70
71 TIMER_STOP(CPU_LOGS);
72}
73
74template <>
75void gravtree<simparticles>::fill_in_export_points(gravpoint_data *exp_point, int i, int no)
76{
77 /* this point has to go to another task */
78 for(int j = 0; j < 3; j++)
79 exp_point->IntPos[j] = Tp->P[i].IntPos[j];
80
81 exp_point->Mass = Tp->P[i].getMass();
82 exp_point->OldAcc = Tp->P[i].OldAcc;
83 exp_point->index = i;
84 exp_point->Type = Tp->P[i].getType();
85#if NSOFTCLASSES > 1
86 exp_point->SofteningClass = Tp->P[i].getSofteningClass();
87#endif
88 exp_point->no = no;
89#ifndef HIERARCHICAL_GRAVITY
90 if(Tp->TimeBinSynchronized[Tp->P[i].TimeBinGrav])
91 exp_point->ActiveFlag = 1;
92 else
93 exp_point->ActiveFlag = 0;
94#endif
95#if defined(PMGRID) && defined(PLACEHIGHRESREGION)
96 exp_point->InsideOutsideFlag = Tp->P[i].InsideOutsideFlag;
97#endif
98}
99
100#if defined(LIGHTCONE) && (defined(LIGHTCONE_PARTICLES_GROUPS) || defined(LIGHTCONE_IMAGE_COMP_HSML_VELDISP))
101
102template <>
103void gravtree<lcparticles>::fill_in_export_points(gravpoint_data *exp_point, int i, int no)
104{
105 /* this point has to go to another task */
106 for(int j = 0; j < 3; j++)
107 exp_point->IntPos[j] = Tp->P[i].IntPos[j];
108
109 exp_point->Mass = Tp->P[i].getMass();
110 exp_point->OldAcc = 0;
111 exp_point->index = i;
112 exp_point->Type = Tp->P[i].getType();
113#if NSOFTCLASSES > 1
114 exp_point->SofteningClass = Tp->P[i].getSofteningClass();
115#endif
116 exp_point->no = no;
117}
118
119#endif
120
125template <typename partset>
126void gravtree<partset>::exchange_topleafdata(void)
127{
128 struct leafnode_data
129 {
130 vector<MyIntPosType> s;
131 MyDouble mass;
132#if(MULTIPOLE_ORDER >= 3) || (MULTIPOLE_ORDER >= 2 && defined(EXTRAPOTTERM))
133 symtensor2<MyDouble> Q2Tensor;
134#endif
135#if(MULTIPOLE_ORDER >= 4) || (MULTIPOLE_ORDER >= 3 && defined(EXTRAPOTTERM))
136 symtensor3<MyDouble> Q3Tensor;
137#endif
138#if(MULTIPOLE_ORDER >= 5) || (MULTIPOLE_ORDER >= 4 && defined(EXTRAPOTTERM))
139 symtensor4<MyDouble> Q4Tensor;
140#endif
141#if(MULTIPOLE_ORDER >= 5 && defined(EXTRAPOTTERM))
142 symtensor5<MyDouble> Q5Tensor;
143#endif
144#ifdef FMM
145 float MinOldAcc;
146#endif
147 unsigned char not_empty;
148#if NSOFTCLASSES > 1
149 unsigned char maxsofttype;
150 unsigned char minsofttype;
151#endif
152 unsigned char level;
153#if defined(PMGRID) && defined(PLACEHIGHRESREGION)
154 unsigned char overlap_flag : 2;
155#endif
156 };
157 leafnode_data *glob_leaf_node_data, *loc_leaf_node_data;
158
159 glob_leaf_node_data = (leafnode_data *)Mem.mymalloc("glob_leaf_node_data", D->NTopleaves * sizeof(leafnode_data));
160
161 /* share the pseudo-particle data accross CPUs */
162 int *recvcounts = (int *)Mem.mymalloc("recvcounts", sizeof(int) * D->NTask);
163 int *recvoffset = (int *)Mem.mymalloc("recvoffset", sizeof(int) * D->NTask);
164 int *bytecounts = (int *)Mem.mymalloc("bytecounts", sizeof(int) * D->NTask);
165 int *byteoffset = (int *)Mem.mymalloc("byteoffset", sizeof(int) * D->NTask);
166
167 for(int task = 0; task < D->NTask; task++)
168 recvcounts[task] = 0;
169
170 for(int n = 0; n < D->NTopleaves; n++)
171 recvcounts[D->TaskOfLeaf[n]]++;
172
173 for(int task = 0; task < D->NTask; task++)
174 bytecounts[task] = recvcounts[task] * sizeof(leafnode_data);
175
176 recvoffset[0] = 0;
177 byteoffset[0] = 0;
178
179 for(int task = 1; task < D->NTask; task++)
180 {
181 recvoffset[task] = recvoffset[task - 1] + recvcounts[task - 1];
182 byteoffset[task] = byteoffset[task - 1] + bytecounts[task - 1];
183 }
184
185 loc_leaf_node_data = (leafnode_data *)Mem.mymalloc("loc_leaf_node_data", recvcounts[D->ThisTask] * sizeof(leafnode_data));
186
187 int idx = 0;
188
189 for(int n = 0; n < D->NTopleaves; n++)
190 {
191 if(D->TaskOfLeaf[n] == D->ThisTask)
192 {
193 int no = NodeIndex[n];
194 gravnode *nop = &TopNodes[no];
195
196 leafnode_data *locp = &loc_leaf_node_data[idx];
197
198 /* read out the multipole moments from the local base cells */
199 locp->s = nop->s;
200 locp->mass = nop->mass;
201#if(MULTIPOLE_ORDER >= 3) || (MULTIPOLE_ORDER >= 2 && defined(EXTRAPOTTERM))
202 locp->Q2Tensor = nop->Q2Tensor;
203#endif
204#if(MULTIPOLE_ORDER >= 4) || (MULTIPOLE_ORDER >= 3 && defined(EXTRAPOTTERM))
205 locp->Q3Tensor = nop->Q3Tensor;
206#endif
207#if(MULTIPOLE_ORDER >= 5) || (MULTIPOLE_ORDER >= 4 && defined(EXTRAPOTTERM))
208 locp->Q4Tensor = nop->Q4Tensor;
209#endif
210#if(MULTIPOLE_ORDER >= 5 && defined(EXTRAPOTTERM))
211 locp->Q5Tensor = nop->Q5Tensor;
212#endif
213#if NSOFTCLASSES > 1
214 locp->maxsofttype = nop->maxsofttype;
215 locp->minsofttype = nop->minsofttype;
216#endif
217#ifdef FMM
218 locp->MinOldAcc = nop->MinOldAcc;
219#endif
220 locp->not_empty = nop->not_empty;
221 locp->level = nop->level;
222#if defined(PMGRID) && defined(PLACEHIGHRESREGION)
223 locp->overlap_flag = nop->overlap_flag;
224#endif
225 idx++;
226 }
227 }
228
229 // optimise this step - only need to update this once per shared memory node
230
231 myMPI_Allgatherv(loc_leaf_node_data, bytecounts[D->ThisTask], MPI_BYTE, glob_leaf_node_data, bytecounts, byteoffset, MPI_BYTE,
232 D->Communicator);
233
234 for(int task = 0; task < D->NTask; task++)
235 recvcounts[task] = 0;
236
237 for(int n = 0; n < D->NTopleaves; n++)
238 {
239 int task = D->TaskOfLeaf[n];
240 if(task != D->ThisTask)
241 {
242 int no = NodeIndex[n];
243 gravnode *nop = &TopNodes[no];
244
245 idx = recvoffset[task] + recvcounts[task]++;
246 leafnode_data *globp = &glob_leaf_node_data[idx];
247
248 nop->s = globp->s;
249 nop->mass = globp->mass;
250#if(MULTIPOLE_ORDER >= 3) || (MULTIPOLE_ORDER >= 2 && defined(EXTRAPOTTERM))
251 nop->Q2Tensor = globp->Q2Tensor;
252#endif
253#if(MULTIPOLE_ORDER >= 4) || (MULTIPOLE_ORDER >= 3 && defined(EXTRAPOTTERM))
254 nop->Q3Tensor = globp->Q3Tensor;
255#endif
256#if(MULTIPOLE_ORDER >= 5) || (MULTIPOLE_ORDER >= 4 && defined(EXTRAPOTTERM))
257 nop->Q4Tensor = globp->Q4Tensor;
258#endif
259#if(MULTIPOLE_ORDER >= 5 && defined(EXTRAPOTTERM))
260 nop->Q5Tensor = globp->Q5Tensor;
261#endif
262#if NSOFTCLASSES > 1
263 nop->maxsofttype = globp->maxsofttype;
264 nop->minsofttype = globp->minsofttype;
265#endif
266 nop->not_empty = globp->not_empty;
267#ifdef FMM
268 nop->MinOldAcc = globp->MinOldAcc;
269#endif
270 nop->level = globp->level;
271#if defined(PMGRID) && defined(PLACEHIGHRESREGION)
272 nop->overlap_flag = globp->overlap_flag;
273#endif
274 }
275 }
276
277 Mem.myfree(loc_leaf_node_data);
278 Mem.myfree(byteoffset);
279 Mem.myfree(bytecounts);
280 Mem.myfree(recvoffset);
281 Mem.myfree(recvcounts);
282 Mem.myfree(glob_leaf_node_data);
283}
284
291template <typename partset>
292void gravtree<partset>::update_node_recursive(int no, int sib, int mode)
293{
294 if(!(no >= MaxPart && no < MaxPart + MaxNodes)) /* are we an internal node? */
295 Terminate("no internal node\n");
296
297 gravnode *nop = get_nodep(no);
298
299 if(mode == TREE_MODE_TOPLEVEL)
300 {
301 int p = nop->nextnode;
302
303 /* if the next node is not a top-level node, we have reached a leaf node, and we need to do nothing */
304 if(p < MaxPart || p >= FirstNonTopLevelNode)
305 return;
306 }
307
308 MyReal mass = 0;
309 vector<MyReal> s(0.0);
310
311#if(MULTIPOLE_ORDER >= 3) || (MULTIPOLE_ORDER >= 2 && defined(EXTRAPOTTERM))
312 symtensor2<MyReal> Q2Tensor(0.0);
313#endif
314#if(MULTIPOLE_ORDER >= 4) || (MULTIPOLE_ORDER >= 3 && defined(EXTRAPOTTERM))
315 symtensor3<MyReal> Q3Tensor(0.0);
316#endif
317#if(MULTIPOLE_ORDER >= 5) || (MULTIPOLE_ORDER >= 4 && defined(EXTRAPOTTERM))
318 symtensor4<MyReal> Q4Tensor(0.0);
319#endif
320#if(MULTIPOLE_ORDER >= 5 && defined(EXTRAPOTTERM))
321 symtensor5<MyReal> Q5Tensor(0.0);
322#endif
323#if NSOFTCLASSES > 1
324 unsigned char maxsofttype = NSOFTCLASSES + NSOFTCLASSES_HYDRO;
325 unsigned char minsofttype = NSOFTCLASSES + NSOFTCLASSES_HYDRO + 1;
326#endif
327 unsigned char not_empty = 0;
328#ifdef FMM
329 float minOldAcc = MAX_FLOAT_NUMBER;
330#endif
331
332 int p = nop->nextnode;
333
334 while(p != nop->sibling)
335 {
336 if(p >= 0)
337 {
338 if(p >= MaxPart && p < MaxPart + MaxNodes) /* we have an internal node */
339 {
340 int nextsib = get_nodep(p)->sibling;
341
342 update_node_recursive(p, nextsib, mode);
343 }
344
345 if(p < MaxPart) /* a particle */
346 {
347 vector<MyReal> dxyz;
348 Tp->nearest_image_intpos_to_pos(Tp->P[p].IntPos, nop->center.da, dxyz.da);
349
350 MyReal m = Tp->P[p].getMass();
351
352 mass += m;
353 s += m * dxyz;
354
355#if(MULTIPOLE_ORDER >= 3) || (MULTIPOLE_ORDER >= 2 && defined(EXTRAPOTTERM))
356 vector<MyReal> mdxyz = m * dxyz;
357 symtensor2<MyReal> mdxyz2(mdxyz, dxyz);
358 Q2Tensor += mdxyz2;
359#endif
360#if(MULTIPOLE_ORDER >= 4) || (MULTIPOLE_ORDER >= 3 && defined(EXTRAPOTTERM))
361 symtensor3<MyReal> mdxyz3(dxyz, mdxyz2);
362 Q3Tensor += mdxyz3;
363#endif
364#if(MULTIPOLE_ORDER >= 5) || (MULTIPOLE_ORDER >= 4 && defined(EXTRAPOTTERM))
365 symtensor4<MyReal> mdxyz4(dxyz, mdxyz3);
366 Q4Tensor += mdxyz4;
367#endif
368#if(MULTIPOLE_ORDER >= 5 && defined(EXTRAPOTTERM))
369 symtensor5<MyReal> mdxyz5(dxyz, mdxyz4);
370 Q5Tensor += mdxyz5;
371#endif
372 not_empty = 1;
373#ifndef HIERARCHICAL_GRAVITY
374 if(Tp->getTimeBinSynchronized(Tp->P[p].getTimeBinGrav()))
375#endif
376 {
377#ifdef FMM
378 if(minOldAcc > Tp->P[p].getOldAcc())
379 minOldAcc = Tp->P[p].getOldAcc();
380#endif
381 }
382
383#if NSOFTCLASSES > 1
384 if(All.ForceSoftening[maxsofttype] < All.ForceSoftening[Tp->P[p].getSofteningClass()])
385 maxsofttype = Tp->P[p].getSofteningClass();
386 if(All.ForceSoftening[minsofttype] > All.ForceSoftening[Tp->P[p].getSofteningClass()])
387 minsofttype = Tp->P[p].getSofteningClass();
388#endif
389 p = Nextnode[p];
390 }
391 else if(p < MaxPart + MaxNodes) /* an internal node */
392 {
393 gravnode *noptr = get_nodep(p);
394
395 vector<MyReal> dxyz;
396 Tp->nearest_image_intpos_to_pos(noptr->s.da, nop->center.da, dxyz.da);
397
398 MyReal m = noptr->mass;
399
400 mass += m;
401 s += m * dxyz;
402
403#if(MULTIPOLE_ORDER >= 3) || (MULTIPOLE_ORDER >= 2 && defined(EXTRAPOTTERM))
404 Q2Tensor += noptr->Q2Tensor;
405
406 vector<MyReal> mdxyz = m * dxyz;
407 symtensor2<MyReal> mdxyz2(mdxyz, dxyz);
408 Q2Tensor += mdxyz2;
409#endif
410#if(MULTIPOLE_ORDER >= 4) || (MULTIPOLE_ORDER >= 3 && defined(EXTRAPOTTERM))
411 Q3Tensor += noptr->Q3Tensor;
412
413 symtensor3<MyReal> mdxyz3(dxyz, mdxyz2);
414 Q3Tensor += mdxyz3;
415
416 Q3Tensor += outer_prod_sum(noptr->Q2Tensor, dxyz);
417#endif
418#if(MULTIPOLE_ORDER >= 5) || (MULTIPOLE_ORDER >= 4 && defined(EXTRAPOTTERM))
419 Q4Tensor += noptr->Q4Tensor;
420
421 symtensor4<MyReal> mdxyz4(dxyz, mdxyz3);
422 Q4Tensor += mdxyz4;
423
424 Q4Tensor += outer_prod_sum(noptr->Q3Tensor, dxyz);
425 symtensor2<MyReal> dxyz2(dxyz, dxyz);
426 Q4Tensor += outer_prod_sum(noptr->Q2Tensor, dxyz2);
427#endif
428#if(MULTIPOLE_ORDER >= 5 && defined(EXTRAPOTTERM))
429 Q5Tensor += noptr->Q5Tensor;
430
431 symtensor5<MyReal> mdxyz5(dxyz, mdxyz4);
432 Q5Tensor += mdxyz5;
433
434 Q5Tensor += outer_prod_sum(noptr->Q4Tensor, dxyz);
435 Q5Tensor += outer_prod_sum(noptr->Q3Tensor, dxyz2);
436 symtensor3<MyReal> dxyz3(dxyz, dxyz2);
437 Q5Tensor += outer_prod_sum(dxyz3, noptr->Q2Tensor);
438#endif
439
440#if NSOFTCLASSES > 1
441 if(All.ForceSoftening[maxsofttype] < All.ForceSoftening[noptr->maxsofttype])
442 maxsofttype = noptr->maxsofttype;
443 if(All.ForceSoftening[minsofttype] > All.ForceSoftening[noptr->minsofttype])
444 minsofttype = noptr->minsofttype;
445#endif
446 not_empty |= noptr->not_empty;
447#ifdef FMM
448 if(minOldAcc > noptr->MinOldAcc)
449 minOldAcc = noptr->MinOldAcc;
450#endif
451
452 p = noptr->sibling;
453 }
454 else if(p < MaxPart + MaxNodes + D->NTopleaves) /* a pseudo particle */
455 {
456 /* we are processing a local leaf-node which does not have any particles.
457 * can continue to the next element, which should end the work.
458 */
459 p = Nextnode[p - MaxNodes];
460 }
461 else
462 {
463 /* an imported point */
464 int n = p - ImportedNodeOffset;
465
466 if(n >= NumPartImported)
467 Terminate("n=%d >= NumPartImported=%d MaxPart=%d MaxNodes=%d D->NTopleaves=%d", n, NumPartImported, MaxPart,
468 MaxNodes, D->NTopleaves);
469
470 vector<MyReal> dxyz;
471 Tp->nearest_image_intpos_to_pos(Points[n].IntPos, nop->center.da, dxyz.da);
472
473 MyReal m = Points[n].Mass;
474
475 mass += m;
476 s += m * dxyz;
477#if(MULTIPOLE_ORDER >= 3) || (MULTIPOLE_ORDER >= 2 && defined(EXTRAPOTTERM))
478 vector<MyReal> mdxyz = m * dxyz;
479 symtensor2<MyReal> mdxyz2(mdxyz, dxyz);
480 Q2Tensor += mdxyz2;
481#endif
482#if(MULTIPOLE_ORDER >= 4) || (MULTIPOLE_ORDER >= 3 && defined(EXTRAPOTTERM))
483 symtensor3<MyReal> mdxyz3(dxyz, mdxyz2);
484 Q3Tensor += mdxyz3;
485#endif
486#if(MULTIPOLE_ORDER >= 5) || (MULTIPOLE_ORDER >= 4 && defined(EXTRAPOTTERM))
487 symtensor4<MyReal> mdxyz4(dxyz, mdxyz3);
488 Q4Tensor += mdxyz4;
489#endif
490#if(MULTIPOLE_ORDER >= 5 && defined(EXTRAPOTTERM))
491 symtensor5<MyReal> mdxyz5(dxyz, mdxyz4);
492 Q5Tensor += mdxyz5;
493#endif
494 not_empty = 1;
495#ifndef HIERARCHICAL_GRAVITY
496 if(Points[n].ActiveFlag)
497#endif
498 {
499#ifdef FMM
500 if(minOldAcc > Points[n].OldAcc)
501 minOldAcc = Points[n].OldAcc;
502#endif
503 }
504#if NSOFTCLASSES > 1
505 if(All.ForceSoftening[maxsofttype] < All.ForceSoftening[Points[n].SofteningClass])
506 maxsofttype = Points[n].SofteningClass;
507 if(All.ForceSoftening[minsofttype] > All.ForceSoftening[Points[n].SofteningClass])
508 minsofttype = Points[n].SofteningClass;
509#endif
510 p = Nextnode[p - MaxNodes];
511 }
512 }
513 }
514
515 if(mass)
516 {
517 s *= (1 / mass);
518 }
519
520 nop->mass = mass;
521
522 vector<MySignedIntPosType> off;
523 Tp->pos_to_signedintpos(s.da, off.da);
524
525 nop->s[0] = off[0] + nop->center[0];
526 nop->s[1] = off[1] + nop->center[1];
527 nop->s[2] = off[2] + nop->center[2];
528
529#if(MULTIPOLE_ORDER >= 3) || (MULTIPOLE_ORDER >= 2 && defined(EXTRAPOTTERM))
530 vector<MyReal> ms = mass * s;
531 symtensor2<MyReal> ms2(ms, s);
532 Q2Tensor -= ms2;
533 nop->Q2Tensor = Q2Tensor;
534#endif
535
536#if(MULTIPOLE_ORDER >= 4) || (MULTIPOLE_ORDER >= 3 && defined(EXTRAPOTTERM))
537 symtensor3<MyReal> ms3(s, ms2);
538 Q3Tensor -= ms3;
539 Q3Tensor -= outer_prod_sum(Q2Tensor, s);
540 nop->Q3Tensor = Q3Tensor;
541#endif
542
543#if(MULTIPOLE_ORDER >= 5) || (MULTIPOLE_ORDER >= 4 && defined(EXTRAPOTTERM))
544 symtensor4<MyReal> ms4(s, ms3);
545 Q4Tensor -= ms4;
546 Q4Tensor -= outer_prod_sum(Q3Tensor, s);
547 symtensor2<MyReal> s2(s, s);
548 Q4Tensor -= outer_prod_sum(Q2Tensor, s2);
549 nop->Q4Tensor = Q4Tensor;
550#endif
551
552#if(MULTIPOLE_ORDER >= 5 && defined(EXTRAPOTTERM))
553 symtensor5<MyReal> ms5(s, ms4);
554 Q5Tensor -= ms5;
555 Q5Tensor -= outer_prod_sum(Q4Tensor, s);
556 Q5Tensor -= outer_prod_sum(Q3Tensor, s2);
557 symtensor3<MyReal> s3(s, s2);
558 Q5Tensor -= outer_prod_sum(s3, Q2Tensor);
559 nop->Q5Tensor = Q5Tensor;
560#endif
561
562#if NSOFTCLASSES > 1
563 nop->maxsofttype = maxsofttype;
564 nop->minsofttype = minsofttype;
565#endif
566 nop->cannot_be_opened_locally = 0;
567 nop->not_empty = not_empty;
568#ifdef FMM
569 nop->MinOldAcc = minOldAcc;
570#endif
571#if defined(PMGRID) && defined(PLACEHIGHRESREGION)
572 MyIntPosType halflen = ((MyIntPosType)1) << ((BITS_FOR_POSITIONS - 1) - nop->level);
573 nop->overlap_flag = Tp->check_high_res_overlap(nop->center.da, halflen);
574#endif
575}
576
583template <typename partset>
584void gravtree<partset>::set_softenings(void)
585{
586 if(All.ComovingIntegrationOn)
587 {
588 for(int i = 0; i < NSOFTCLASSES; i++)
589 if(All.SofteningComoving[i] * All.Time > All.SofteningMaxPhys[i])
590 All.SofteningTable[i] = All.SofteningMaxPhys[i] / All.Time;
591 else
592 All.SofteningTable[i] = All.SofteningComoving[i];
593 }
594 else
595 {
596 for(int i = 0; i < NSOFTCLASSES; i++)
597 All.SofteningTable[i] = All.SofteningComoving[i];
598 }
599
600#ifdef ADAPTIVE_HYDRO_SOFTENING
601 for(int i = 0; i < NSOFTCLASSES_HYDRO; i++)
602 All.SofteningTable[i + NSOFTCLASSES] = All.MinimumComovingHydroSoftening * pow(All.AdaptiveHydroSofteningSpacing, i);
603
604 if(All.AdaptiveHydroSofteningSpacing < 1)
605 Terminate("All.AdaptiveHydroSofteningSpacing < 1");
606
607 /* we check that type=0 has its own slot 0 in the softening types, so that only gas masses are stored there */
608 if(All.SofteningClassOfPartType[0] != 0)
609 Terminate("All.SofteningClassOfPartType[0] != 0");
610
611 for(int i = 1; i < NTYPES; i++)
612 if(All.SofteningClassOfPartType[i] == All.SofteningClassOfPartType[0])
613 Terminate("i=%d: All.SofteningClassOfPartType[i] == All.SofteningClassOfPartType[0]", i);
614#endif
615
616 for(int i = 0; i < NSOFTCLASSES + NSOFTCLASSES_HYDRO; i++)
617 All.ForceSoftening[i] = 2.8 * All.SofteningTable[i];
618
619 All.ForceSoftening[NSOFTCLASSES + NSOFTCLASSES_HYDRO] =
620 0; /* important - this entry is actually used in the tree construction for the search of the maximum softening in a node */
621 All.ForceSoftening[NSOFTCLASSES + NSOFTCLASSES_HYDRO + 1] =
622 MAX_FLOAT_NUMBER; /* important - this entry is actually used in the tree construction for the search of the maximum softening in
623 a node */
624}
625
626/* make sure that we instantiate the template */
627#include "../data/simparticles.h"
628template class gravtree<simparticles>;
629
630/* make sure that we instantiate the template */
631#if defined(LIGHTCONE) && (defined(LIGHTCONE_PARTICLES_GROUPS) || defined(LIGHTCONE_IMAGE_COMP_HSML_VELDISP))
632#include "../data/lcparticles.h"
633template class gravtree<lcparticles>;
634#endif
All
global_data_all_processes All
Definition: main.cc:40
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
vector::da
T da[3]
Definition: symtensors.h:106
NSOFTCLASSES_HYDRO
#define NSOFTCLASSES_HYDRO
Definition: constants.h:321
NSOFTCLASSES
#define NSOFTCLASSES
Definition: constants.h:312
NTYPES
#define NTYPES
Definition: constants.h:308
MAX_FLOAT_NUMBER
#define MAX_FLOAT_NUMBER
Definition: constants.h:79
MyDouble
float MyDouble
Definition: dtypes.h:87
MyIntPosType
uint32_t MyIntPosType
Definition: dtypes.h:35
MyReal
double MyReal
Definition: dtypes.h:82
BITS_FOR_POSITIONS
#define BITS_FOR_POSITIONS
Definition: dtypes.h:37
myMPI_Allgatherv
int myMPI_Allgatherv(void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int *recvcount, int *displs, MPI_Datatype recvtype, MPI_Comm comm)
Definition: hypercube_allgatherv.cc:26
TIMER_START
#define TIMER_START(counter)
Starts the timer counter.
Definition: logs.h:197
TIMER_STOP
#define TIMER_STOP(counter)
Stops the timer counter.
Definition: logs.h:220
Terminate
#define Terminate(...)
Definition: macros.h:15
sumup_large_ints
void sumup_large_ints(int n, int *src, long long *res, MPI_Comm comm)
Definition: sums_and_minmax.cc:33
Mem
memory Mem
Definition: main.cc:44
half_float::detail::pow
expr pow(half base, half exp)
Definition: half.hpp:2803
basenode::level
unsigned char level
Definition: tree.h:64
basenode::cannot_be_opened_locally
std::atomic< unsigned char > cannot_be_opened_locally
Definition: tree.h:70
basenode::not_empty
unsigned char not_empty
Definition: tree.h:73
basenode::nextnode
int nextnode
Definition: tree.h:57
basenode::center
vector< MyIntPosType > center
Definition: tree.h:51
basenode::sibling
int sibling
Definition: tree.h:53
global_data_all_processes::SofteningClassOfPartType
int SofteningClassOfPartType[NTYPES]
Definition: allvars.h:250
global_data_all_processes::SofteningMaxPhys
double SofteningMaxPhys[NSOFTCLASSES]
Definition: allvars.h:253
global_data_all_processes::ComovingIntegrationOn
int ComovingIntegrationOn
Definition: allvars.h:148
global_data_all_processes::SofteningComoving
double SofteningComoving[NSOFTCLASSES]
Definition: allvars.h:252
global_data_all_processes::ForceSoftening
double ForceSoftening[NSOFTCLASSES+NSOFTCLASSES_HYDRO+2]
Definition: allvars.h:258
global_data_all_processes::Time
double Time
Definition: allvars.h:179
global_data_all_processes::SofteningTable
double SofteningTable[NSOFTCLASSES+NSOFTCLASSES_HYDRO]
Definition: allvars.h:256
gravnode::mass
MyDouble mass
Definition: gravtree.h:65
gravnode::s
vector< MyIntPosType > s
Definition: gravtree.h:66
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::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
outer_prod_sum
symtensor3< typename which_return< T1, T2 >::type > outer_prod_sum(const symtensor2< T1 > &D, const vector< T2 > &v)
Definition: symtensors.h:1621
TREE_MODE_TOPLEVEL
#define TREE_MODE_TOPLEVEL
Definition: tree.h:20