GADGET-4
subfind_properties.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
16#include <mpi.h>
17#include <algorithm>
18#include <cmath>
19#include <cstdio>
20#include <cstdlib>
21#include <cstring>
22
23#include "../data/allvars.h"
24#include "../data/dtypes.h"
25#include "../data/intposconvert.h"
26#include "../data/mymalloc.h"
27#include "../domain/domain.h"
28#include "../fof/fof.h"
29#include "../gravtree/gravtree.h"
30#include "../logs/timer.h"
31#include "../main/simulation.h"
32#include "../mpi_utils/mpi_utils.h"
33#include "../sort/cxxsort.h"
34#include "../sort/parallel_sort.h"
35#include "../subfind/subfind.h"
36#include "../system/system.h"
37#include "../time_integration/driftfac.h"
38
39template <>
40void fof<simparticles>::subfind_get_factors(double &fac_vel_to_phys, double &fac_hubbleflow, double &fac_comov_to_phys)
41{
43 {
44 fac_vel_to_phys = 1.0 / All.Time; // converts to physical velocity
45 fac_hubbleflow = Driftfac.hubble_function(All.Time);
46 fac_comov_to_phys = All.Time; // converts comoving distance to physical distance
47 }
48 else
49 { // non-comoving integration
50 fac_vel_to_phys = 1;
51 fac_hubbleflow = 0;
52 fac_comov_to_phys = 1;
53 }
54}
55
56#if defined(LIGHTCONE) && defined(LIGHTCONE_PARTICLES_GROUPS)
57template <>
58void fof<lcparticles>::subfind_get_factors(double &fac_vel_to_phys, double &fac_hubbleflow, double &fac_comov_to_phys)
59{
60 fac_vel_to_phys = 1.0; // the velocities on the lightcone are already peculiar velocities
61 fac_hubbleflow = Driftfac.hubble_function(Ascale);
62 fac_comov_to_phys = Ascale;
63}
64#endif
65
66template <typename partset>
67int fof<partset>::subfind_determine_sub_halo_properties(int *d, int num, subhalo_properties *subhalo, MPI_Comm Communicator)
68{
69 /* get the local communication context */
70 int commNTask, commThisTask;
71 MPI_Comm_size(Communicator, &commNTask);
72 MPI_Comm_rank(Communicator, &commThisTask);
73
74 double fac_vel_to_phys, fac_hubbleflow, fac_comov_to_phys;
75 subfind_get_factors(fac_vel_to_phys, fac_hubbleflow, fac_comov_to_phys);
76
77 typename partset::pdata *P = Tp->P;
78 subfind_data *PS = Tp->PS;
79
80 double mass_tab[NTYPES], halfmassradtype[NTYPES];
81 long long len_type[NTYPES];
82 for(int j = 0; j < NTYPES; j++)
83 {
84 mass_tab[j] = 0;
85 len_type[j] = 0;
86 halfmassradtype[j] = 0;
87 }
88
89#ifdef STARFORMATION
90 double sfr = 0;
91#endif
92
93 int minindex = -1;
94 double minpot = MAX_DOUBLE_NUMBER;
95 for(int i = 0; i < num; i++)
96 {
97 int p = d[i];
98 if(PS[p].u.s.u.DM_Potential < minpot || minindex == -1)
99 {
100 minpot = PS[p].u.s.u.DM_Potential;
101 minindex = p;
102 }
103 len_type[P[p].getType()]++;
104
105#ifdef STARFORMATION
106 if(P[p].getType() == 0)
107 sfr += Tp->SphP[p].Sfr;
108#endif
109 }
110
111 MPI_Allreduce(MPI_IN_PLACE, len_type, NTYPES, MPI_LONG_LONG, MPI_SUM, Communicator);
112
113 double *minpotlist = (double *)Mem.mymalloc("minpotlist", commNTask * sizeof(double));
114 MPI_Allgather(&minpot, 1, MPI_DOUBLE, minpotlist, 1, MPI_DOUBLE, Communicator);
115
116 int mincpu = -1;
117 minpot = MAX_DOUBLE_NUMBER;
118 for(int i = 0; i < commNTask; i++)
119 if(minpotlist[i] < minpot)
120 {
121 mincpu = i;
122 minpot = minpotlist[mincpu];
123 }
124
125 Mem.myfree(minpotlist);
126
127 if(mincpu < 0)
128 Terminate("mincpu < 0");
129
130 MyIntPosType intpos[3];
131
132 if(commThisTask == mincpu)
133 for(int j = 0; j < 3; j++)
134 intpos[j] = P[minindex].IntPos[j];
135
136 MPI_Bcast(intpos, 3 * sizeof(MyIntPosType), MPI_BYTE, mincpu, Communicator);
137
138#ifdef STARFORMATION
139 MPI_Allreduce(MPI_IN_PLACE, &sfr, 1, MPI_DOUBLE, MPI_SUM, Communicator);
140#endif
141
142 /* pos[] now holds the position of minimum potential */
143 /* we'll take it that as the center */
144
145 /* determine the particle ID with the smallest binding energy */
146 minindex = -1;
147 minpot = MAX_DOUBLE_NUMBER;
148 for(int i = 0; i < num; i++)
149 {
150 int p = d[i];
151 if(PS[p].v.DM_BindingEnergy < minpot || minindex == -1)
152 {
153 minpot = PS[p].v.DM_BindingEnergy;
154 minindex = p;
155 }
156 }
157
158 MyIDType mostboundid;
159
160 minpotlist = (double *)Mem.mymalloc("minpotlist", commNTask * sizeof(double));
161 MPI_Allgather(&minpot, 1, MPI_DOUBLE, minpotlist, 1, MPI_DOUBLE, Communicator);
162
163 mincpu = -1;
164 minpot = MAX_DOUBLE_NUMBER;
165 for(int i = 0; i < commNTask; i++)
166 if(minpotlist[i] < minpot)
167 {
168 mincpu = i;
169 minpot = minpotlist[mincpu];
170 }
171
172 Mem.myfree(minpotlist);
173
174 if(mincpu < 0)
175 Terminate("mincpu < 0");
176
177 int marked = 0;
178
179 if(commThisTask == mincpu)
180 {
181 mostboundid = P[minindex].ID.get();
182#if defined(SUBFIND_ORPHAN_TREATMENT)
183 if(!P[minindex].ID.is_previously_most_bound())
184 {
185 P[minindex].ID.mark_as_formerly_most_bound();
186 marked++;
187 }
188#endif
189 }
190
191 MPI_Bcast(&mostboundid, sizeof(mostboundid), MPI_BYTE, mincpu, Communicator);
192
193 /* let's get bulk velocity and the center-of-mass */
194 /* here we still take all particles */
195
196 double s[3] = {0, 0, 0}, v[3] = {0, 0, 0};
197 double mass = 0;
198
199#if defined(SUBFIND_ORPHAN_TREATMENT)
200 int totprevmostboundlen = 0;
201#endif
202
203 for(int i = 0; i < num; i++)
204 {
205 int p = d[i];
206
207 double dxyz[3];
208 Tp->nearest_image_intpos_to_pos(P[p].IntPos, intpos, dxyz);
209
210 s[0] += Tp->P[p].getMass() * dxyz[0];
211 s[1] += Tp->P[p].getMass() * dxyz[1];
212 s[2] += Tp->P[p].getMass() * dxyz[2];
213
214 for(int j = 0; j < 3; j++)
215 v[j] += Tp->P[p].getMass() * P[p].Vel[j];
216
217 mass += Tp->P[p].getMass();
218
219 int ptype = P[p].getType();
220 mass_tab[ptype] += Tp->P[p].getMass();
221#if defined(SUBFIND_ORPHAN_TREATMENT)
222 if(P[p].ID.is_previously_most_bound())
223 totprevmostboundlen++;
224#endif
225 }
226
227 double stot[3], vtot[3], masstot, mass_tabtot[NTYPES];
228
229 MPI_Allreduce(s, stot, 3, MPI_DOUBLE, MPI_SUM, Communicator);
230 MPI_Allreduce(&mass, &masstot, 1, MPI_DOUBLE, MPI_SUM, Communicator);
231 MPI_Allreduce(v, vtot, 3, MPI_DOUBLE, MPI_SUM, Communicator);
232 MPI_Allreduce(mass_tab, mass_tabtot, NTYPES, MPI_DOUBLE, MPI_SUM, Communicator);
233
234#if defined(SUBFIND_ORPHAN_TREATMENT)
235 MPI_Allreduce(MPI_IN_PLACE, &totprevmostboundlen, 1, MPI_INT, MPI_SUM, Communicator);
236#endif
237
238 mass = masstot;
239 for(int j = 0; j < 3; j++)
240 {
241 s[j] = stot[j];
242 v[j] = vtot[j];
243 }
244
245 for(int j = 0; j < NTYPES; j++)
246 mass_tab[j] = mass_tabtot[j];
247
248 double vel[3];
249 for(int j = 0; j < 3; j++)
250 {
251 s[j] /= mass; /* center of mass */
252 v[j] /= mass;
253 vel[j] = fac_vel_to_phys * v[j];
254 }
255
256 MySignedIntPosType off[3];
257 Tp->pos_to_signedintpos(s, off);
258
259 MyIntPosType int_cm[3];
260
261 int_cm[0] = off[0] + intpos[0];
262 int_cm[1] = off[1] + intpos[1];
263 int_cm[2] = off[2] + intpos[2];
264
265 double cm[3];
266 Tp->intpos_to_pos(int_cm, cm);
267
268 double disp = 0, lx = 0, ly = 0, lz = 0;
269
270 sort_r2list *rr_list = (sort_r2list *)Mem.mymalloc("rr_list", sizeof(sort_r2list) * (num + 1));
271
272 for(int i = 0; i < num; i++)
273 {
274 int p = d[i];
275
276 double dx[3], dv[3];
277
278 /* relative to center of mass */
279
280 Tp->nearest_image_intpos_to_pos(P[p].IntPos, int_cm, dx);
281
282 dx[0] *= fac_comov_to_phys;
283 dx[1] *= fac_comov_to_phys;
284 dx[2] *= fac_comov_to_phys;
285
286 for(int j = 0; j < 3; j++)
287 {
288 dv[j] = fac_vel_to_phys * (P[p].Vel[j] - v[j]);
289 dv[j] += fac_hubbleflow * dx[j];
290
291 disp += Tp->P[p].getMass() * dv[j] * dv[j];
292 }
293
294 lx += Tp->P[p].getMass() * (dx[1] * dv[2] - dx[2] * dv[1]);
295 ly += Tp->P[p].getMass() * (dx[2] * dv[0] - dx[0] * dv[2]);
296 lz += Tp->P[p].getMass() * (dx[0] * dv[1] - dx[1] * dv[0]);
297
298 /* for rotation curve computation, take minimum of potential as center */
299
300 double dxyz[3];
301 Tp->nearest_image_intpos_to_pos(P[p].IntPos, intpos, dxyz);
302
303 for(int j = 0; j < 3; j++)
304 dxyz[j] *= fac_comov_to_phys;
305
306 double r2 = dxyz[0] * dxyz[0] + dxyz[1] * dxyz[1] + dxyz[2] * dxyz[2];
307
308 double r = sqrt(r2);
309
310 rr_list[i].mass = Tp->P[p].getMass();
311 rr_list[i].r = r;
312 }
313
314 double spintot[3], disploc = disp;
315 double spinloc[3] = {lx, ly, lz};
316 MPI_Allreduce(spinloc, spintot, 3, MPI_DOUBLE, MPI_SUM, Communicator);
317 MPI_Allreduce(&disploc, &disp, 1, MPI_DOUBLE, MPI_SUM, Communicator);
318 lx = spintot[0];
319 ly = spintot[1];
320 lz = spintot[2];
321
322 double spin[3] = {lx / mass, ly / mass, lz / mass};
323
324 double veldisp = sqrt(disp / (3 * mass)); /* convert to 1d velocity dispersion */
325
326 mycxxsort_parallel(rr_list, rr_list + num, subfind_compare_dist_rotcurve, Communicator);
327
328 /* calculate cumulative mass */
329 for(int i = 1; i < num; i++)
330 rr_list[i].mass = rr_list[i - 1].mass + rr_list[i].mass;
331
332 double halfmassrad = 0;
333 double max = 0, maxrad = 0;
334
335 double mass_part = 0;
336 if(num)
337 mass_part = rr_list[num - 1].mass;
338 double *masslist = (double *)Mem.mymalloc("masslist", commNTask * sizeof(double));
339 MPI_Allgather(&mass_part, 1, MPI_DOUBLE, masslist, 1, MPI_DOUBLE, Communicator);
340
341 double massbefore = 0;
342 for(int i = 0; i < commThisTask; i++)
343 massbefore += masslist[i];
344
345 for(int i = 0; i < num; i++)
346 rr_list[i].mass += massbefore;
347
348 Mem.myfree(masslist);
349
350 /* now calculate rotation curve maximum and half mass radius */
351
352 double halfmassrad_loc = 0;
353 sort_r2list *rr_lowlist = (sort_r2list *)Mem.mymalloc("rr_lowlist", commNTask * sizeof(sort_r2list));
354 sort_r2list low_element;
355 if(num > 0)
356 low_element = rr_list[0];
357 else
358 {
359 low_element.mass = 0;
360 low_element.r = 0;
361 }
362 MPI_Allgather(&low_element, sizeof(sort_r2list), MPI_BYTE, rr_lowlist, sizeof(sort_r2list), MPI_BYTE, Communicator);
363
364 rr_list[num].mass = 0;
365 rr_list[num].r = 0;
366
367 for(int j = commThisTask + 1; j < commNTask; j++)
368 if(rr_lowlist[j].mass > 0)
369 {
370 rr_list[num] = rr_lowlist[j];
371 break;
372 }
373
374 Mem.myfree(rr_lowlist);
375
376 int *numlist = (int *)Mem.mymalloc("numlist", commNTask * sizeof(int));
377 MPI_Allgather(&num, 1, MPI_INT, numlist, 1, MPI_INT, Communicator);
378
379 int nbefore = 0;
380 for(int i = 0; i < commThisTask; i++)
381 nbefore += numlist[i];
382
383 for(int i = num - 1; i >= 0; i--)
384 {
385 if((i + nbefore) > 5 && rr_list[i].mass > max * rr_list[i].r)
386 {
387 max = rr_list[i].mass / rr_list[i].r;
388 maxrad = rr_list[i].r;
389 }
390
391 if(rr_list[i].mass < 0.5 * mass && rr_list[i + 1].mass >= 0.5 * mass)
392 halfmassrad_loc = 0.5 * (rr_list[i].r + rr_list[i + 1].r);
393 }
394
395 Mem.myfree(numlist);
396
397 MPI_Allreduce(&halfmassrad_loc, &halfmassrad, 1, MPI_DOUBLE, MPI_MAX, Communicator);
398 double *maxlist = (double *)Mem.mymalloc("maxlist", commNTask * sizeof(double));
399 double *maxradlist = (double *)Mem.mymalloc("maxradlist", commNTask * sizeof(double));
400 MPI_Allgather(&max, 1, MPI_DOUBLE, maxlist, 1, MPI_DOUBLE, Communicator);
401 MPI_Allgather(&maxrad, 1, MPI_DOUBLE, maxradlist, 1, MPI_DOUBLE, Communicator);
402
403 max = maxrad = 0;
404 for(int i = 0; i < commNTask; i++)
405 {
406 if(maxlist[i] > max)
407 {
408 max = maxlist[i];
409 maxrad = maxradlist[i];
410 }
411 }
412 Mem.myfree(maxradlist);
413 Mem.myfree(maxlist);
414
415 double vmax = sqrt(All.G * max);
416 double vmaxrad = maxrad;
417
418 Mem.myfree(rr_list);
419
420 /* half mass radii for different types */
421
422 int len_type_loc[NTYPES];
423 for(int j = 0; j < NTYPES; j++)
424 len_type_loc[j] = 0;
425
426 for(int i = 0; i < num; i++)
427 {
428 int p = d[i];
429 int ptype = P[p].getType();
430 len_type_loc[ptype]++;
431 }
432
433 for(int type = 0; type < NTYPES; type++)
434 {
435 sort_r2list *rr_list = (sort_r2list *)Mem.mymalloc("rr_list", sizeof(sort_r2list) * (len_type_loc[type] + 1));
436 int itmp = 0;
437 for(int i = 0; i < num; i++)
438 {
439 int p = d[i];
440
441 int ptype = P[p].getType();
442
443 if(ptype == type)
444 {
445 double dxyz[3];
446 Tp->nearest_image_intpos_to_pos(P[p].IntPos, intpos, dxyz);
447
448 for(int j = 0; j < 3; j++)
449 dxyz[j] *= fac_comov_to_phys;
450
451 double r2 = dxyz[0] * dxyz[0] + dxyz[1] * dxyz[1] + dxyz[2] * dxyz[2];
452 double r = sqrt(r2);
453
454 rr_list[itmp].mass = Tp->P[p].getMass();
455 rr_list[itmp].r = r;
456 itmp++;
457 }
458 }
459
460 if(itmp != len_type_loc[type])
461 Terminate("should not occur: %d %d", itmp, len_type_loc[type]);
462
463 mycxxsort_parallel(rr_list, rr_list + len_type_loc[type], subfind_compare_dist_rotcurve, Communicator);
464
465 /* calculate cumulative mass */
466 for(int i = 1; i < len_type_loc[type]; i++)
467 rr_list[i].mass = rr_list[i - 1].mass + rr_list[i].mass;
468
469 double mass_part = 0;
470 if(len_type_loc[type])
471 mass_part = rr_list[len_type_loc[type] - 1].mass;
472 double *masslist = (double *)Mem.mymalloc("masslist", commNTask * sizeof(double));
473 MPI_Allgather(&mass_part, 1, MPI_DOUBLE, masslist, 1, MPI_DOUBLE, Communicator);
474
475 double massbefore = 0;
476 for(int i = 0; i < commThisTask; i++)
477 massbefore += masslist[i];
478
479 for(int i = 0; i < len_type_loc[type]; i++)
480 rr_list[i].mass += massbefore;
481
482 Mem.myfree(masslist);
483
484 /* now calculate half mass radii */
485 double halfmassrad_loc = 0;
486 sort_r2list *rr_lowlist = (sort_r2list *)Mem.mymalloc("rr_lowlist", commNTask * sizeof(sort_r2list));
487 sort_r2list low_element;
488 if(len_type_loc[type] > 0)
489 low_element = rr_list[0];
490 else
491 {
492 low_element.mass = 0;
493 low_element.r = 0;
494 }
495
496 MPI_Allgather(&low_element, sizeof(sort_r2list), MPI_BYTE, rr_lowlist, sizeof(sort_r2list), MPI_BYTE, Communicator);
497
498 rr_list[len_type_loc[type]].mass = 0;
499 rr_list[len_type_loc[type]].r = 0;
500 for(int j = commThisTask + 1; j < commNTask; j++)
501 if(rr_lowlist[j].mass > 0)
502 {
503 rr_list[len_type_loc[type]] = rr_lowlist[j];
504 break;
505 }
506
507 Mem.myfree(rr_lowlist);
508
509 for(int i = len_type_loc[type] - 1; i >= 0; i--)
510 {
511 if(rr_list[i].mass < 0.5 * mass_tab[type] && rr_list[i + 1].mass >= 0.5 * mass_tab[type])
512 halfmassrad_loc = 0.5 * (rr_list[i].r + rr_list[i + 1].r);
513 }
514
515 MPI_Allreduce(&halfmassrad_loc, &halfmassradtype[type], 1, MPI_DOUBLE, MPI_MAX, Communicator);
516
517 Mem.myfree(rr_list);
518 }
519
520 /* properties of star forming gas */
521#ifdef STARFORMATION
522 double gasMassSfr = 0;
523 for(int i = 0; i < num; i++)
524 {
525 int p = d[i];
526
527 if(P[p].getType() == 0)
528 if(Tp->SphP[p].Sfr > 0)
529 gasMassSfr += P[p].getMass();
530 }
531#endif
532
533#ifdef STARFORMATION
534 double gasMassSfrtot;
535 MPI_Allreduce(&gasMassSfr, &gasMassSfrtot, 1, MPI_DOUBLE, MPI_SUM, Communicator);
536 gasMassSfr = gasMassSfrtot;
537#endif
538
539 long long totlen;
540 sumup_large_ints(1, &num, &totlen, Communicator);
541
542 /* now store the calculated properties in the subhalo structure */
543 if(commThisTask == 0)
544 {
545 subhalo->Len = totlen;
546 subhalo->Mass = mass;
547
548 for(int j = 0; j < NTYPES; j++)
549 {
550 subhalo->MassType[j] = mass_tab[j];
551 subhalo->LenType[j] = len_type[j];
552 subhalo->SubHalfMassRadType[j] = halfmassradtype[j];
553 }
554
555 double pos[3];
556 fof_get_halo_position(intpos, pos);
557
558 for(int j = 0; j < 3; j++)
559 {
560 subhalo->IntPos[j] = intpos[j];
561 subhalo->Pos[j] = pos[j];
562 subhalo->Vel[j] = vel[j];
563 subhalo->CM[j] = cm[j];
564 subhalo->Spin[j] = spin[j];
565 }
566
567 subhalo->SubMostBoundID = mostboundid;
568 subhalo->SubVelDisp = veldisp;
569 subhalo->SubVmax = vmax;
570 subhalo->SubVmaxRad = vmaxrad;
571 subhalo->SubHalfMassRad = halfmassrad;
572
573#if defined(SUBFIND_ORPHAN_TREATMENT)
574 subhalo->SubhaloLenPrevMostBnd = totprevmostboundlen;
575#endif
576#ifdef STARFORMATION
577 subhalo->Sfr = sfr;
578 subhalo->GasMassSfr = gasMassSfr;
579#endif
580 }
581
582 return marked;
583}
584
585/* now make sure that the following classes are really instantiated, otherwise we may get a linking problem */
586#include "../data/simparticles.h"
587template class fof<simparticles>;
588
589#if defined(LIGHTCONE) && defined(LIGHTCONE_PARTICLES_GROUPS)
590#include "../data/lcparticles.h"
591template class fof<lcparticles>;
592#endif
593
594#endif
global_data_all_processes All
Definition: main.cc:40
static double hubble_function(double a)
Definition: driftfac.h:39
#define NTYPES
Definition: constants.h:308
#define MAX_DOUBLE_NUMBER
Definition: constants.h:81
uint32_t MyIntPosType
Definition: dtypes.h:35
int32_t MySignedIntPosType
Definition: dtypes.h:36
unsigned int MyIDType
Definition: dtypes.h:68
#define Terminate(...)
Definition: macros.h:15
driftfac Driftfac
Definition: main.cc:41
void sumup_large_ints(int n, int *src, long long *res, MPI_Comm comm)
memory Mem
Definition: main.cc:44
expr sqrt(half arg)
Definition: half.hpp:2777
double mycxxsort_parallel(T *begin, T *end, Comp comp, MPI_Comm comm)
union subfind_data::@0 u