GADGET-4
subfind_processing.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 <gsl/gsl_math.h>
17#include <mpi.h>
18#include <algorithm>
19#include <cmath>
20#include <cstdio>
21#include <cstdlib>
22#include <cstring>
23
24#include "../data/allvars.h"
25#include "../data/dtypes.h"
26#include "../data/intposconvert.h"
27#include "../data/mymalloc.h"
28#include "../domain/domain.h"
29#include "../fof/fof.h"
30#include "../gravtree/gravtree.h"
31#include "../logs/timer.h"
32#include "../main/simulation.h"
33#include "../mpi_utils/mpi_utils.h"
34#include "../sort/cxxsort.h"
35#include "../sort/parallel_sort.h"
36#include "../sort/peano.h"
37#include "../subfind/subfind.h"
38#include "../system/system.h"
39
40template <typename partset>
41void fof<partset>::subfind_processing(domain<partset> *SubDomain, domain_options mode)
42{
43 double t0 = Logs.second();
44
45 if(mode == COLL_SUBFIND)
46 {
47 /* make a sanity check: We should have exactly 1 group, stored on the root of the processor subset if we collectively do a group
48 */
49 if(SubThisTask == 0 && Ngroups != 1)
50 Terminate("Ngroups=%d != 1 SubNTask=%d SubThisTask=%d", Ngroups, SubNTask, SubThisTask);
51
52 if(SubThisTask != 0 && Ngroups != 0)
53 Terminate("Ngroups=%d != 0 SubNTask=%d SubThisTask=%d", Ngroups, SubNTask, SubThisTask);
54
55 subfind_collective_printf("SUBFIND: root-task=%d: Collectively doing halo %lld of length %lld on %d processors.\n", ThisTask,
56 Group[0].GroupNr, (long long)Group[0].Len, SubNTask);
57
58 if(SubThisTask == 0)
59 {
60 GroupNr = Group[0].GroupNr;
61 Ascale = Group[0].Ascale;
62 }
63
64 for(int i = 0; i < Tp->NumPart; i++)
65 if(Tp->PS[i].GroupNr.get() == GroupNr)
66 Tp->PS[i].DomainFlag = 1;
67 else
68 Tp->PS[i].DomainFlag = 0;
69
70 /* tell everybody in the set the group number */
71 MPI_Bcast(&GroupNr, 1, MPI_LONG_LONG, 0, SubComm);
72
73 /* tell everybody in the set the group's scale factor */
74 MPI_Bcast(&Ascale, 1, MPI_DOUBLE, 0, SubComm);
75 }
76 else
77 {
78 for(int i = 0; i < Tp->NumPart; i++)
79 Tp->PS[i].DomainFlag = 1;
80
81 if(SubNTask != 1)
82 Terminate("Strange: SubNTask=%d Ngroups=%d SubThisTask=%d (expect to be a single processor here)", SubNTask, Ngroups,
83 SubThisTask);
84 }
85
86 /* Create a domain decomposition for the sub-communicator and the particles in it.
87 * For the serial algorithm, this will be trivial, for collectively treated groups, only particles in the group get a gravity weight.
88 * The outline of the toplevel tree nodes resulting from the domain decomposition can be used to create a gravity tree.
89 */
90
91 if(SubDomain->NumNodes != 0)
92 Terminate("SubDomain.NumNodes=%d\n", SubDomain->NumNodes);
93
94 double ta = Logs.second();
95 SubDomain->domain_decomposition(mode);
96 double tb = Logs.second();
97
98 mpi_printf("SUBFIND: subdomain decomposition took %g sec\n", Logs.timediff(ta, tb));
99
100 if(mode == COLL_SUBFIND)
101 SubDomain->particle_exchange_based_on_PS(SubComm);
102
103 for(int i = 0; i < Tp->NumPart; i++)
104 Tp->PS[i].SubRankInGr = INT_MAX; /* set a default that is larger than any reasonable group number */
105
106 /* now let us sort according to GroupNr and Density. This step will temporarily break the association with SphP[] and other arrays!
107 */
108 submp = (submp_data *)Mem.mymalloc_movable(&submp, "submp", sizeof(submp_data) * Tp->NumPart);
109 for(int i = 0; i < Tp->NumPart; i++)
110 {
111 submp[i].index = i;
112 submp[i].GroupNr = Tp->PS[i].GroupNr.get();
113#ifndef SUBFIND_HBT
114 submp[i].DM_Density = Tp->PS[i].u.s.u.DM_Density;
115#endif
116 }
117 mycxxsort(submp, submp + Tp->NumPart, subfind_compare_submp_GroupNr_DM_Density);
118
119 /* In this index list, we store the indices of the local particles in the group that we currently need to process (it is for sure
120 * shorter than NumPart) */
121 IndexList = (int *)Mem.mymalloc_movable(&IndexList, "IndexList", Tp->NumPart * sizeof(int));
122
123 MPI_Comm SingleComm;
124 int thistask;
125 MPI_Comm_rank(SubDomain->Communicator, &thistask);
126 MPI_Comm_split(SubDomain->Communicator, thistask, thistask, &SingleComm); // create a communicator for single ranks
127
128 /* prepare a domain decomposition for unbinding locally independent ones */
129 domain<partset> SingleDomain{SingleComm, Tp};
130
131 if(SingleDomain.NumNodes != 0)
132 Terminate("SubDomain.NumNodes=%d\n", SingleDomain.NumNodes);
133
134 double taa = Logs.second();
135 SingleDomain.domain_decomposition(SERIAL_SUBFIND);
136 double tbb = Logs.second();
137
138 mpi_printf("SUBFIND: serial subfind subdomain decomposition took %g sec\n", Logs.timediff(taa, tbb));
139
140 double ta0 = Logs.second();
141 if(mode == COLL_SUBFIND)
142 {
143 /* determine the number of local group particles, and fill the list of the indices */
144 NumPartGroup = 0;
145 for(int i = 0; i < Tp->NumPart; i++)
146 if(Tp->PS[i].GroupNr.get() == GroupNr)
147 IndexList[NumPartGroup++] = i;
148
149 /* call the processing of the group */
150#ifdef SUBFIND_HBT
151 subfind_hbt_single_group(SubDomain, &SingleDomain, mode, 0);
152#else
153 subfind_process_single_group(SubDomain, &SingleDomain, mode, 0);
154#endif
155 }
156 else
157 {
158 int i = 0;
159 for(int gr = 0; gr < Ngroups; gr++) /* process all local groups */
160 {
161 GroupNr = Group[gr].GroupNr;
162 Ascale = Group[gr].Ascale;
163
164 /* determine the number of local group particles, and set up the list of the indices */
165 NumPartGroup = 0;
166 for(; i < Tp->NumPart; i++)
167 if(Tp->PS[submp[i].index].GroupNr.get() == GroupNr)
168 IndexList[NumPartGroup++] = submp[i].index;
169 else
170 break;
171
172 /* do local group with Group[] index 'gr' */
173#ifdef SUBFIND_HBT
174 subfind_hbt_single_group(SubDomain, &SingleDomain, mode, gr);
175#else
176 subfind_process_single_group(SubDomain, &SingleDomain, mode, gr);
177#endif
178 }
179 }
180 double tb0 = Logs.second();
181 mpi_printf("SUBFIND: subfind_hbt_single_group() processing for Ngroups=%d took %g sec\n", Ngroups, Logs.timediff(ta0, tb0));
182
183 SingleDomain.domain_free();
184 MPI_Comm_free(&SingleComm);
185
186 Mem.myfree(IndexList);
187 Mem.myfree(submp);
188
189 SubDomain->domain_free();
190
191 double t1 = Logs.second();
192
193 subfind_collective_printf("SUBFIND: root-task=%d: Collective processing of halo %d took %g\n", ThisTask, Group[0].GroupNr,
194 Logs.timediff(t0, t1));
195
196 if(!(mode == COLL_SUBFIND) && ThisTask == 0)
197 mpi_printf("SUBFIND: root-task=%d: Serial processing of halo %d took %g\n", ThisTask, Group[0].GroupNr, Logs.timediff(t0, t1));
198}
199
200/* now make sure that the following classes are really instantiated, otherwise we may get a linking problem */
201#include "../data/simparticles.h"
202template class fof<simparticles>;
203
204#if defined(LIGHTCONE) && defined(LIGHTCONE_PARTICLES_GROUPS)
205#include "../data/lcparticles.h"
206template class fof<lcparticles>;
207#endif
208
209#endif
Definition: domain.h:31
void domain_free(void)
Definition: domain.cc:179
void particle_exchange_based_on_PS(MPI_Comm Communicator)
void domain_decomposition(domain_options mode)
Definition: domain.cc:51
double timediff(double t0, double t1)
Definition: logs.cc:488
double second(void)
Definition: logs.cc:471
int NumNodes
Definition: setcomm.h:37
MPI_Comm Communicator
Definition: setcomm.h:31
double mycxxsort(T *begin, T *end, Tcomp comp)
Definition: cxxsort.h:39
domain_options
Definition: domain.h:22
@ COLL_SUBFIND
Definition: domain.h:24
@ SERIAL_SUBFIND
Definition: domain.h:25
logs Logs
Definition: main.cc:43
#define Terminate(...)
Definition: macros.h:15
memory Mem
Definition: main.cc:44