GADGET-4
subfind_density.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#ifdef SUBFIND
15#ifndef SUBFIND_HBT
16
17#include <gsl/gsl_math.h>
18#include <math.h>
19#include <mpi.h>
20#include <stdio.h>
21#include <stdlib.h>
22#include <string.h>
23#include <sys/stat.h>
24#include <sys/types.h>
25#include <algorithm>
26
27#include "../data/allvars.h"
28#include "../data/dtypes.h"
29#include "../data/intposconvert.h"
30#include "../data/mymalloc.h"
31#include "../domain/domain.h"
32#include "../fof/fof.h"
33#include "../gravtree/gravtree.h"
34#include "../logs/timer.h"
35#include "../main/simulation.h"
36#include "../mpi_utils/generic_comm.h"
37#include "../mpi_utils/mpi_utils.h"
38#include "../subfind/subfind.h"
39#include "../system/system.h"
40
41// Structure for communication during the density computation. Holds data that is sent to other processors.
42struct subdens_in : data_in_generic
43{
44 MyIntPosType IntPos[3];
45 MyFloat Hsml;
46};
47
48/* local data structure that holds results acquired on remote processors */
49struct subdens_out
50{
51 int Ngb;
52 MyFloat Rho;
53#ifdef SUBFIND_STORE_LOCAL_DENSITY
54 MyFloat VelDisp, Vx, Vy, Vz;
55#endif
56};
57
58static int *DM_NumNgb;
59#ifdef SUBFIND_STORE_LOCAL_DENSITY
60static MyFloat *Vx, *Vy, *Vz;
61#endif
62
63template <typename T_tree, typename T_domain, typename T_partset>
64class subdens_comm : public generic_comm<subdens_in, subdens_out, T_tree, T_domain, T_partset>
65{
66 public:
68 using gcomm::D;
69 using gcomm::Thread;
70 using gcomm::Tp; // This makes sure that we can access Tp from the base class without having to use "this->Tp"
71 using gcomm::Tree;
72
73 /* need to call the base class constructor explicitly */
74 subdens_comm(T_domain *dptr, T_tree *tptr, T_partset *pptr) : gcomm(dptr, tptr, pptr) {}
75
76 /* routine that fills the relevant particle/cell data into the input structure defined above */
77 void particle2in(subdens_in *in, int i) override
78 {
79 in->IntPos[0] = Tp->P[i].IntPos[0];
80 in->IntPos[1] = Tp->P[i].IntPos[1];
81 in->IntPos[2] = Tp->P[i].IntPos[2];
82 in->Hsml = Tp->PS[i].v.DM_Hsml;
83 }
84
85 /* routine to store or combine result data */
86 void out2particle(subdens_out *out, int i, int mode) override
87 {
88 if(mode == MODE_LOCAL_PARTICLES) /* initial store */
89 {
90 DM_NumNgb[i] = out->Ngb;
91 Tp->PS[i].u.s.u.DM_Density = out->Rho;
92#ifdef SUBFIND_STORE_LOCAL_DENSITY
93 Vx[i] = out->Vx;
94 Vy[i] = out->Vy;
95 Vz[i] = out->Vz;
96 Tp->PS[i].SubfindVelDisp = out->VelDisp;
97#endif
98 }
99 else /* combine */
100 {
101 DM_NumNgb[i] += out->Ngb;
102 Tp->PS[i].u.s.u.DM_Density += out->Rho;
103#ifdef SUBFIND_STORE_LOCAL_DENSITY
104 Vx[i] += out->Vx;
105 Vy[i] += out->Vy;
106 Vz[i] += out->Vz;
107 Tp->PS[i].SubfindVelDisp += out->VelDisp;
108#endif
109 }
110 }
111
116 int evaluate(int target, int mode, int thread_id, int action, subdens_in *in, int numnodes, node_info *firstnode,
117 subdens_out &out) override
118 {
119 MyIntPosType *intpos = in->IntPos;
120 double hsml = in->Hsml;
121
122 double h2 = hsml * hsml;
123 double hinv = 1.0 / hsml;
124 double hinv3 = hinv * hinv * hinv;
125
126 int numngb = 0;
127 double rhosum = 0;
128#ifdef SUBFIND_STORE_LOCAL_DENSITY
129 double vxsum = 0;
130 double vysum = 0;
131 double vzsum = 0;
132 double v2sum = 0;
133#endif
134
135 for(int k = 0; k < numnodes; k++)
136 {
137 int no;
138 if(mode == MODE_LOCAL_PARTICLES)
139 {
140 no = Tree->MaxPart; /* root node */
141 }
142 else
143 {
144 no = firstnode[k].Node;
145 no = Tree->get_nodep(no)->nextnode; /* open it */
146 }
147
148 unsigned int shmrank = Tree->TreeSharedMem_ThisTask;
149
150 while(no >= 0)
151 {
152 int p, type;
153 double mass, r2;
154 typename T_partset::pdata *P;
155
156 if(no < Tree->MaxPart) /* single particle */
157 {
158 p = no;
159 P = Tree->get_Pp(no, shmrank);
160
161 no = Tree->get_nextnodep(shmrank)[no]; /* note: here shmrank cannot change */
162
163 double dxyz[3];
164 Tp->nearest_image_intpos_to_pos(P->IntPos, intpos, dxyz); /* converts the integer distance to floating point */
165
166 r2 = dxyz[0] * dxyz[0];
167 if(r2 > h2)
168 continue;
169
170 r2 += dxyz[1] * dxyz[1];
171 if(r2 > h2)
172 continue;
173
174 r2 += dxyz[2] * dxyz[2];
175 if(r2 > h2)
176 continue;
177
178 mass = P->getMass();
179 type = P->getType();
180 }
181 else if(no < Tree->MaxPart + Tree->MaxNodes) /* internal node */
182 {
183 if(mode == MODE_IMPORTED_PARTICLES)
184 {
185 if(no < Tree->FirstNonTopLevelNode) /* we reached a top-level node again, which means that we are done with the
186 branch */
187 break;
188 }
189
190 gravnode *current = Tree->get_nodep(no, shmrank);
191
192 no = current->sibling; /* in case the node can be discarded */
193 shmrank = current->sibling_shmrank;
194
195 double dxyz[3];
196 Tp->nearest_image_intpos_to_pos(current->center.da, intpos,
197 dxyz); /* converts the integer distance to floating point */
198
199 double lenhalf = (((MyIntPosType)1) << (BITS_FOR_POSITIONS - 1 - current->level)) * Tp->FacIntToCoord;
200
201 double dist = hsml + lenhalf;
202
203 if(fabs(dxyz[0]) > dist)
204 continue;
205 if(fabs(dxyz[1]) > dist)
206 continue;
207 if(fabs(dxyz[2]) > dist)
208 continue;
209
210 /* now test against the minimal sphere enclosing everything */
211 dist += FACT1 * 2.0 * lenhalf;
212 if(dxyz[0] * dxyz[0] + dxyz[1] * dxyz[1] + dxyz[2] * dxyz[2] > dist * dist)
213 continue;
214
215 no = current->nextnode; /* ok, we need to open the node */
216 shmrank = current->nextnode_shmrank;
217
218 continue;
219 }
220 else if(no >= Tree->ImportedNodeOffset) /* point from imported nodelist */
221 {
222 int n = no - Tree->ImportedNodeOffset;
223 no = Tree->Nextnode[no - Tree->MaxNodes];
224 /* note: here shmrank cannot change */
225
226 double dxyz[3];
227 Tp->nearest_image_intpos_to_pos(Tree->Points[n].IntPos, intpos,
228 dxyz); /* converts the integer distance to floating point */
229
230 r2 = dxyz[0] * dxyz[0];
231 if(r2 > h2)
232 continue;
233
234 r2 += dxyz[1] * dxyz[1];
235 if(r2 > h2)
236 continue;
237
238 r2 += dxyz[2] * dxyz[2];
239 if(r2 > h2)
240 continue;
241
242 mass = Tree->Points[n].Mass;
243 type = Tree->Points[n].Type;
244
245 p = -1;
246 }
247 else /* pseudo particle */
248 {
249 if(mode == MODE_LOCAL_PARTICLES)
250 if(target >= 0) /* if no target is given, export will not occur */
251 Tree->tree_export_node_threads(no, target, &Thread);
252
253 no = Tree->Nextnode[no - Tree->MaxNodes];
254 continue;
255 }
256
257 if(r2 < h2)
258 {
259 if(is_type_primary_link_type(type))
260 {
261 numngb++;
262
263 if(p < 0)
264 Terminate("this should not occur");
265
266#ifdef SUBFIND_STORE_LOCAL_DENSITY
267 vxsum += P->Vel[0];
268 vysum += P->Vel[1];
269 vzsum += P->Vel[2];
270 v2sum += P->Vel[0] * P->Vel[0] + P->Vel[1] * P->Vel[1] + P->Vel[2] * P->Vel[2];
271#endif
272 }
273
274 if(is_type_primary_link_type(type) || is_type_secondary_link_type(type))
275 {
276 double r = sqrt(r2);
277
278 double u = r * hinv, wk;
279
280 if(u < 0.5)
281 wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
282 else
283 wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);
284
285 rhosum += mass * wk;
286 }
287 }
288 }
289 }
290
291 out.Ngb = numngb;
292 out.Rho = rhosum;
293#ifdef SUBFIND_STORE_LOCAL_DENSITY
294 out.Vx = vxsum;
295 out.Vy = vysum;
296 out.Vz = vzsum;
297 out.VelDisp = v2sum;
298#endif
299 return 0;
300 }
301};
302
303template <typename partset>
304double fof<partset>::subfind_density(void)
305{
306 mpi_printf("SUBFIND: finding total densities around all particles\n");
307
308 double tstart = Logs.second();
309
310 DM_NumNgb = (int *)Mem.mymalloc_movable(&DM_NumNgb, "DM_NumNgb", sizeof(int) * Tp->NumPart);
311 MyFloat *Left = (MyFloat *)Mem.mymalloc_movable(&Left, "Left", sizeof(MyFloat) * Tp->NumPart);
312 MyFloat *Right = (MyFloat *)Mem.mymalloc_movable(&Right, "Right", sizeof(MyFloat) * Tp->NumPart);
313
314 int *targetlist = (int *)Mem.mymalloc_movable(&targetlist, "targetlist", sizeof(int) * Tp->NumPart);
315
316#ifdef SUBFIND_STORE_LOCAL_DENSITY
317 Vx = (MyFloat *)Mem.mymalloc("Vx", sizeof(MyFloat) * Tp->NumPart);
318 Vy = (MyFloat *)Mem.mymalloc("Vy", sizeof(MyFloat) * Tp->NumPart);
319 Vz = (MyFloat *)Mem.mymalloc("Vz", sizeof(MyFloat) * Tp->NumPart);
320#endif
321
322 int ntodo = 0;
323
324 for(int i = 0; i < Tp->NumPart; i++)
325 {
326 Left[i] = Right[i] = 0;
327 DM_NumNgb[i] = 0;
328
329 Tp->PS[i].u.s.u.DM_Density = 0;
330
331#ifdef SUBFIND_STORE_LOCAL_DENSITY
332 Tp->PS[i].SubfindHsml = 0;
333 Tp->PS[i].SubfindDensity = 0;
334 Tp->PS[i].SubfindVelDisp = 0;
335 if(is_type_primary_link_type(Tp->P[i].getType()) || is_type_secondary_link_type(Tp->P[i].getType()))
336 targetlist[ntodo++] = i;
337#else
338 if(Tp->PS[i].GroupNr.get() != HALONR_MAX) /* do it only for particles is in a group */
339 targetlist[ntodo++] = i;
340#endif
341 }
342
343 subdens_comm<gravtree<partset>, domain<partset>, partset> commpattern{FoFDomain, &FoFGravTree, Tp};
344
345 /* we will repeat the whole thing for those particles where we didn't find enough neighbours */
346 long long ntot;
347 int iter = 0;
348 do
349 {
350 double t0 = Logs.second();
351
352 commpattern.execute(ntodo, targetlist, MODE_DEFAULT);
353
354 /* do final operations on results */
355 int npleft = 0;
356 for(int n = 0; n < ntodo; n++)
357 {
358 int i = targetlist[n];
359
360 /* now check whether we had enough neighbours */
361 if(abs(DM_NumNgb[i] - All.DesNumNgb) > All.MaxNumNgbDeviation &&
362 ((Right[i] - Left[i]) > 1.0e-4 * Left[i] || Left[i] == 0 || Right[i] == 0))
363 {
364 /* need to redo this particle */
365 targetlist[npleft++] = i;
366
367 if(DM_NumNgb[i] < All.DesNumNgb)
368 Left[i] = (MyFloat)std::max<double>(Tp->PS[i].v.DM_Hsml, Left[i]);
369 else
370 {
371 if(Right[i] != 0)
372 {
373 if(Tp->PS[i].v.DM_Hsml < Right[i])
374 Right[i] = Tp->PS[i].v.DM_Hsml;
375 }
376 else
377 Right[i] = Tp->PS[i].v.DM_Hsml;
378 }
379
380 if(iter >= MAXITER - 10)
381 {
382 double pos[3];
383 Tp->intpos_to_pos(Tp->P[i].IntPos, pos);
384
385 printf("SUBFIND: i=%d task=%d ID=%lld Hsml=%g Left=%g Right=%g Ngbs=%g Right-Left=%g\n pos=(%g|%g|%g)\n", i,
386 ThisTask, (long long)Tp->P[i].ID.get(), Tp->PS[i].v.DM_Hsml, Left[i], Right[i], (double)DM_NumNgb[i],
387 Right[i] - Left[i], pos[0], pos[1], pos[2]);
388 myflush(stdout);
389 }
390
391 if(Right[i] > 0 && Left[i] > 0)
392 Tp->PS[i].v.DM_Hsml = (MyFloat)pow(0.5 * (pow(Left[i], 3) + pow(Right[i], 3)), 1.0 / 3.0);
393 else
394 {
395 if(Right[i] == 0 && Left[i] == 0)
396 Terminate("can't occur");
397
398 if(Right[i] == 0 && Left[i] > 0)
399 Tp->PS[i].v.DM_Hsml *= 1.26;
400
401 if(Right[i] > 0 && Left[i] == 0)
402 Tp->PS[i].v.DM_Hsml /= 1.26;
403 }
404 }
405 }
406
407 ntodo = npleft;
408
409 sumup_large_ints(1, &npleft, &ntot, Communicator);
410
411 double t1 = Logs.second();
412
413 if(ntot > 0)
414 {
415 iter++;
416
417 if(iter > 0)
418 mpi_printf("SUBFIND: ngb iteration %2d: need to repeat for %15lld particles. (took %g sec)\n", iter, ntot,
419 Logs.timediff(t0, t1));
420
421 if(iter > MAXITER)
422 Terminate("failed to converge in neighbour iteration in subfind_density()\n");
423 }
424 }
425 while(ntot > 0);
426
427#ifdef SUBFIND_STORE_LOCAL_DENSITY
428
429 double vel_to_phys = subfind_vel_to_phys_factor();
430
431 for(int i = 0; i < Tp->NumPart; i++)
432 {
433 if(DM_NumNgb[i] > 0)
434 {
435 Vx[i] /= DM_NumNgb[i];
436 Vy[i] /= DM_NumNgb[i];
437 Vz[i] /= DM_NumNgb[i];
438 Tp->PS[i].SubfindVelDisp /= DM_NumNgb[i];
439 Tp->PS[i].SubfindVelDisp = vel_to_phys * sqrt(Tp->PS[i].SubfindVelDisp - Vx[i] * Vx[i] - Vy[i] * Vy[i] - Vz[i] * Vz[i]);
440 }
441 else
442 Tp->PS[i].SubfindVelDisp = 0;
443 }
444
445 Mem.myfree_movable(Vz);
446 Mem.myfree_movable(Vy);
447 Mem.myfree_movable(Vx);
448#endif
449
450 Mem.myfree_movable(targetlist);
451 Mem.myfree_movable(Right);
452 Mem.myfree_movable(Left);
453 Mem.myfree_movable(DM_NumNgb);
454
455#ifdef SUBFIND_STORE_LOCAL_DENSITY
456 for(int i = 0; i < Tp->NumPart; i++)
457 {
458 Tp->PS[i].SubfindHsml = Tp->PS[i].v.DM_Hsml;
459 Tp->PS[i].SubfindDensity = Tp->PS[i].u.s.u.DM_Density;
460 }
461#endif
462
463 double tend = Logs.second();
464 return Logs.timediff(tstart, tend);
465}
466
467template <>
468double fof<simparticles>::subfind_vel_to_phys_factor(void)
469{
471 return 1.0 / All.Time;
472 else
473 return 1.0;
474}
475
476#if defined(LIGHTCONE) && defined(LIGHTCONE_PARTICLES_GROUPS)
477template <>
478double fof<lcparticles>::subfind_vel_to_phys_factor(void)
479{
480 return 1.0;
481}
482#endif
483
484template <typename partset>
485void fof<partset>::subfind_density_hsml_guess(void) /* set the initial guess for the smoothing length */
486{
487 double hsml_prev = 0;
488
489 for(int i = 0; i < Tp->NumPart; i++)
490 {
491 if(is_type_primary_link_type(Tp->P[i].getType()))
492 {
493 int no = FoFGravTree.Father[i];
494
495 while(8.0 * All.DesNumNgb * Tp->P[i].getMass() > FoFGravTree.get_nodep(no)->mass)
496 {
497 int p = FoFGravTree.get_nodep(no)->father;
498
499 if(p < 0)
500 break;
501
502 no = p;
503 }
504
505 double len = (((MyIntPosType)1) << (BITS_FOR_POSITIONS - FoFGravTree.get_nodep(no)->level)) * Tp->FacIntToCoord;
506
507 Tp->PS[i].v.DM_Hsml = hsml_prev =
508 (pow(3.0 / (4.0 * M_PI) * All.DesNumNgb * Tp->P[i].getMass() / FoFGravTree.get_nodep(no)->mass, 1.0 / 3.0) * len);
509
510 if(Tp->PS[i].v.DM_Hsml == 0)
511 {
512 double pos[3];
513 Tp->intpos_to_pos(Tp->P[i].IntPos, pos);
514
515 Terminate(
516 "zero hsml guess: Hsml=0 task=%d i=%d no=%d Nodes[no].len=%g Nodes[no].mass=%g P[i].Mass=%g type=%d ID=%llu "
517 "pos=(%g|%g|%g)\n",
518 ThisTask, i, no, len, FoFGravTree.get_nodep(no)->mass, Tp->P[i].getMass(), Tp->P[i].getType(),
519 (long long)Tp->P[i].ID.get(), pos[0], pos[1], pos[2]);
520 }
521 }
522 else
523 {
524 if(hsml_prev)
525 Tp->PS[i].v.DM_Hsml = hsml_prev;
526 else
527 Tp->PS[i].v.DM_Hsml = All.SofteningTable[All.SofteningClassOfPartType[Tp->P[i].getType()]];
528 }
529 }
530}
531
532/* now make sure that the following classes are really instantiated, otherwise we may get a linking problem */
533#include "../data/simparticles.h"
534template class fof<simparticles>;
535
536#if defined(LIGHTCONE) && defined(LIGHTCONE_PARTICLES_GROUPS)
537#include "../data/lcparticles.h"
538template class fof<lcparticles>;
539#endif
540
541#endif
542#endif
global_data_all_processes All
Definition: main.cc:40
Definition: domain.h:31
virtual void out2particle(T_out *out, int i, int mode)=0
virtual void particle2in(T_in *in, int i)=0
virtual int evaluate(int target, int mode, int thread_id, int action, T_in *in, int numnodes, node_info *firstnode, T_out &out)=0
double timediff(double t0, double t1)
Definition: logs.cc:488
double second(void)
Definition: logs.cc:471
#define KERNEL_COEFF_2
Definition: constants.h:371
#define KERNEL_COEFF_5
Definition: constants.h:374
#define MAXITER
Definition: constants.h:305
#define MODE_DEFAULT
Definition: constants.h:23
#define KERNEL_COEFF_1
Definition: constants.h:370
#define MODE_IMPORTED_PARTICLES
Definition: constants.h:22
#define FACT1
Definition: constants.h:435
#define MODE_LOCAL_PARTICLES
Definition: constants.h:21
#define M_PI
Definition: constants.h:56
float MyFloat
Definition: dtypes.h:86
uint32_t MyIntPosType
Definition: dtypes.h:35
#define BITS_FOR_POSITIONS
Definition: dtypes.h:37
#define HALONR_MAX
Definition: idstorage.h:20
logs Logs
Definition: main.cc:43
#define Terminate(...)
Definition: macros.h:15
void sumup_large_ints(int n, int *src, long long *res, MPI_Comm comm)
memory Mem
Definition: main.cc:44
half abs(half arg)
Definition: half.hpp:2620
half fabs(half arg)
Definition: half.hpp:2627
expr pow(half base, half exp)
Definition: half.hpp:2803
expr sqrt(half arg)
Definition: half.hpp:2777
unsigned char level
Definition: tree.h:64
unsigned char nextnode_shmrank
Definition: tree.h:66
int nextnode
Definition: tree.h:57
vector< MyIntPosType > center
Definition: tree.h:51
unsigned char sibling_shmrank
Definition: tree.h:65
int sibling
Definition: tree.h:53
int SofteningClassOfPartType[NTYPES]
Definition: allvars.h:250
double SofteningTable[NSOFTCLASSES+NSOFTCLASSES_HYDRO]
Definition: allvars.h:256
Definition: tree.h:77
int Node
Definition: tree.h:78
void myflush(FILE *fstream)
Definition: system.cc:125