GADGET-4
foftree_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/mymalloc.h"
23#include "../domain/domain.h"
24#include "../fof/foftree.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 "../sort/peano.h"
31#include "../system/system.h"
32#include "../time_integration/driftfac.h"
33#include "../time_integration/timestep.h"
34
35template <typename partset>
37{
38 double numnodes = NumNodes, tot_numnodes;
39 MPI_Reduce(&numnodes, &tot_numnodes, 1, MPI_DOUBLE, MPI_SUM, 0, D->Communicator);
40
41 D->mpi_printf("FOFTREE: Ngb-tree construction done. took %g sec <numnodes>=%g NTopnodes=%d NTopleaves=%d\n", Buildtime,
42 tot_numnodes / D->NTask, D->NTopnodes, D->NTopleaves);
43}
44
45template <typename partset>
47{
48 Terminate("we don't expect to get here");
49}
50
51template <typename partset>
53{
54 struct leafnode_data
55 {
56 MyIntPosType range_min[3];
57 MyIntPosType range_max[3];
58 };
59 leafnode_data *glob_leaf_node_data, *loc_leaf_node_data;
60
61 glob_leaf_node_data = (leafnode_data *)Mem.mymalloc("glob_leaf_node_data", D->NTopleaves * sizeof(leafnode_data));
62
63 /* share the pseudo-particle data accross CPUs */
64 int *recvcounts = (int *)Mem.mymalloc("recvcounts", sizeof(int) * D->NTask);
65 int *recvoffset = (int *)Mem.mymalloc("recvoffset", sizeof(int) * D->NTask);
66 int *bytecounts = (int *)Mem.mymalloc("bytecounts", sizeof(int) * D->NTask);
67 int *byteoffset = (int *)Mem.mymalloc("byteoffset", sizeof(int) * D->NTask);
68
69 for(int task = 0; task < D->NTask; task++)
70 recvcounts[task] = 0;
71
72 for(int n = 0; n < D->NTopleaves; n++)
73 recvcounts[D->TaskOfLeaf[n]]++;
74
75 for(int task = 0; task < D->NTask; task++)
76 bytecounts[task] = recvcounts[task] * sizeof(leafnode_data);
77
78 recvoffset[0] = 0;
79 byteoffset[0] = 0;
80
81 for(int task = 1; task < D->NTask; task++)
82 {
83 recvoffset[task] = recvoffset[task - 1] + recvcounts[task - 1];
84 byteoffset[task] = byteoffset[task - 1] + bytecounts[task - 1];
85 }
86
87 loc_leaf_node_data = (leafnode_data *)Mem.mymalloc("loc_leaf_node_data", recvcounts[D->ThisTask] * sizeof(leafnode_data));
88
89 int idx = 0;
90
91 for(int n = 0; n < D->NTopleaves; n++)
92 {
93 if(D->TaskOfLeaf[n] == D->ThisTask)
94 {
95 int no = NodeIndex[n];
96 fofnode *nop = &TopNodes[no];
97
98 leafnode_data *locp = &loc_leaf_node_data[idx];
99
100 for(int k = 0; k < 3; k++)
101 {
102 locp->range_min[k] = nop->range_min[k];
103 locp->range_max[k] = nop->range_max[k];
104 }
105
106 idx++;
107 }
108 }
109
110 myMPI_Allgatherv(loc_leaf_node_data, bytecounts[D->ThisTask], MPI_BYTE, glob_leaf_node_data, bytecounts, byteoffset, MPI_BYTE,
111 D->Communicator);
112
113 for(int task = 0; task < D->NTask; task++)
114 recvcounts[task] = 0;
115
116 for(int n = 0; n < D->NTopleaves; n++)
117 {
118 int task = D->TaskOfLeaf[n];
119 if(task != D->ThisTask)
120 {
121 int no = NodeIndex[n];
122 fofnode *nop = &TopNodes[no];
123
124 int idx = recvoffset[task] + recvcounts[task]++;
125 leafnode_data *globp = &glob_leaf_node_data[idx];
126
127 for(int k = 0; k < 3; k++)
128 {
129 nop->range_min[k] = globp->range_min[k];
130 nop->range_max[k] = globp->range_max[k];
131 }
132 }
133 }
134
135 Mem.myfree(loc_leaf_node_data);
136 Mem.myfree(byteoffset);
137 Mem.myfree(bytecounts);
138 Mem.myfree(recvoffset);
139 Mem.myfree(recvcounts);
140 Mem.myfree(glob_leaf_node_data);
141}
142
148template <typename partset>
149void foftree<partset>::update_node_recursive(int no, int sib, int mode)
150{
151 MyIntPosType range_min[3];
152 MyIntPosType range_max[3];
153
154 if(!(no >= MaxPart && no < MaxPart + MaxNodes)) /* are we an internal node? */
155 Terminate("no internal node\n");
156
157 fofnode *nop = get_nodep(no);
158
159 if(mode == TREE_MODE_TOPLEVEL)
160 {
161 int p = nop->nextnode;
162
163 /* if the next node is not a top-level node, we have reached a leaf node, and we need to do nothing */
164 if(p < MaxPart || p >= FirstNonTopLevelNode)
165 return;
166 }
167
168 for(int k = 0; k < 3; k++)
169 {
170 range_min[k] = ~((MyIntPosType)0);
171 range_max[k] = 0;
172 }
173
174 int p = nop->nextnode;
175
176 while(p != nop->sibling)
177 {
178 if(p >= 0)
179 {
180 if(p >= MaxPart && p < MaxPart + MaxNodes) /* we have an internal node */
181 {
182 int nextsib = get_nodep(p)->sibling;
183
184 update_node_recursive(p, nextsib, mode);
185 }
186
187 if(p < MaxPart) /* a particle */
188 {
189 for(int k = 0; k < 3; k++)
190 {
191 if(range_min[k] > Tp->P[p].IntPos[k])
192 range_min[k] = Tp->P[p].IntPos[k];
193
194 if(range_max[k] < Tp->P[p].IntPos[k])
195 range_max[k] = Tp->P[p].IntPos[k];
196 }
197
198 p = Nextnode[p];
199 }
200 else if(p < MaxPart + MaxNodes) /* an internal node */
201 {
202 fofnode *nop = get_nodep(p);
203
204 for(int k = 0; k < 3; k++)
205 {
206 if(range_min[k] > nop->range_min[k])
207 range_min[k] = nop->range_min[k];
208
209 if(range_max[k] < nop->range_max[k])
210 range_max[k] = nop->range_max[k];
211 }
212
213 p = nop->sibling;
214 }
215 else if(p < MaxPart + MaxNodes + D->NTopleaves) /* a pseudo particle */
216 {
217 /* we are processing a local leaf-node which does not have any particles.
218 * can continue to the next element, which should end the work.
219 */
220 p = Nextnode[p - MaxNodes];
221 }
222 else
223 {
224 /* an imported point */
225
226 Terminate("Ups!");
227 p = Nextnode[p - MaxNodes];
228 }
229 }
230 }
231
232 for(int k = 0; k < 3; k++)
233 {
234 nop->range_min[k] = range_min[k];
235 nop->range_max[k] = range_max[k];
236 }
237}
238
239/* make sure that we instantiate the template */
240#include "../data/simparticles.h"
241template class foftree<simparticles>;
242
243#if defined(LIGHTCONE) && (defined(LIGHTCONE_PARTICLES_GROUPS) || defined(LIGHTCONE_IMAGE_COMP_HSML_VELDISP))
244/* make sure that we instantiate the template */
245#include "../data/lcparticles.h"
246template class foftree<lcparticles>;
247#endif
void exchange_topleafdata(void) override
void fill_in_export_points(fofpoint_data *exp_point, int i, int no) override
void update_node_recursive(int no, int sib, int mode) override
void report_log_message(void) override
uint32_t MyIntPosType
Definition: dtypes.h:35
int myMPI_Allgatherv(void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int *recvcount, int *displs, MPI_Datatype recvtype, MPI_Comm comm)
#define Terminate(...)
Definition: macros.h:15
memory Mem
Definition: main.cc:44
int nextnode
Definition: tree.h:57
int sibling
Definition: tree.h:53
MyIntPosType range_max[3]
Definition: foftree.h:33
MyIntPosType range_min[3]
Definition: foftree.h:32
#define TREE_MODE_TOPLEVEL
Definition: tree.h:20