GADGET-4
postproc_descendants.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 MERGERTREE
15
16#include <gsl/gsl_rng.h>
17#include <hdf5.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 <unistd.h>
26
27#include "../data/allvars.h"
28#include "../data/dtypes.h"
29#include "../data/mymalloc.h"
30#include "../fof/fof.h"
31#include "../fof/fof_io.h"
32#include "../io/hdf5_util.h"
33#include "../io/io.h"
34#include "../logs/timer.h"
35#include "../main/main.h"
36#include "../main/simulation.h"
37#include "../mergertree/io_readsnap.h"
38#include "../mergertree/mergertree.h"
39#include "../mpi_utils/mpi_utils.h"
40#include "../sort/parallel_sort.h"
41#include "../subfind/subfind.h"
42#include "../system/system.h"
43
44// this is used in FOF postprocessing to assign the previous subhalo length to particles
45void mergertree::get_previous_size_of_subhalo_for_each_particle(int num)
46{
47 if(num >= 0)
48 {
49 mpi_printf(
50 "SUBFIND / MERGERTREE: We are loading the previous group catalogue to set the size of previous subhalos for each "
51 "particle\n");
52
53 domain<simparticles> Domain{Communicator, Sp};
54 fof<simparticles> FoF{Communicator, Sp, &Domain};
55
56 /* load previous snapshot and group/subhalo catalogues */
57 fof_io<simparticles> FoF_IO{&FoF, Communicator, All.SnapFormat};
58 FoF_IO.fof_subfind_load_groups(num);
59
60 readsnap_io Snap_IO{this, this->Communicator, All.SnapFormat};
61 Snap_IO.mergertree_read_snap_ids(num);
62
63 /* make sure that group catalog and snapshot are in file order */
64 mycxxsort_parallel(MtrP, MtrP + MtrP_NumPart, compare_MtrP_FileOffset, Communicator);
65 mycxxsort_parallel(FoF.Group, FoF.Group + FoF.Ngroups, compare_Group_FileOffset, Communicator);
66 mycxxsort_parallel(FoF.Subhalo, FoF.Subhalo + FoF.Nsubhalos, compare_Subhalo_FileOffset, Communicator);
67
68 /* now assign (global) group and subhalo numbers to each particle belonging to the particular structure */
69 mergertree_assign_group_numbers(&FoF);
70
71 Mem.myfree_movable(FoF.Subhalo);
72 Mem.myfree_movable(FoF.Group);
73
74 /* get PrevSizeOfSubhalo from PrevMtrP by matching IDs */
75 mergertree_match_ids_of_current_snap();
76 }
77 else
78 {
79 for(int i = 0; i < Sp->NumPart; i++)
80 {
81 Sp->P[i].PrevSubhaloNr.set(HALONR_MAX);
82 Sp->P[i].PrevSizeOfSubhalo.set(0);
83 }
84 }
85}
86
87void mergertree::descendants_in_postprocessing(int num)
88{
89 if(num - 1 < 0)
90 Terminate("cannot execute for the given snapnum");
91
92 domain<simparticles> Domain{Communicator, Sp};
93 fof<simparticles> FoF{Communicator, Sp, &Domain};
94
95 /* load previous snapshot and group/subhalo catalogues */
96
97 fof_io<simparticles> FoF_IO{&FoF, Communicator, All.SnapFormat};
98 FoF_IO.fof_subfind_load_groups(num - 1);
99
100 readsnap_io Snap_IO{this, this->Communicator, All.SnapFormat};
101 Snap_IO.mergertree_read_snap_ids(num - 1);
102
103 /* make sure that group catalog and snapshot are in file order */
104 mycxxsort_parallel(MtrP, MtrP + MtrP_NumPart, compare_MtrP_FileOffset, Communicator);
105 mycxxsort_parallel(FoF.Group, FoF.Group + FoF.Ngroups, compare_Group_FileOffset, Communicator);
106 mycxxsort_parallel(FoF.Subhalo, FoF.Subhalo + FoF.Nsubhalos, compare_Subhalo_FileOffset, Communicator);
107
108 if(FoF_IO.LegacyFormat)
109 {
110 mpi_printf("\nFOF/SUBFIND: Legacy format from Arepo detected, trying to adjust for this.\n");
111 FoF.subfind_redetermine_groupnr();
112 }
113
114 /* now assign (global) group and subhalo numbers to each particle belonging to the particular structure */
115 mergertree_assign_group_numbers(&FoF);
116
117 Mem.myfree_movable(FoF.Subhalo);
118 Mem.myfree_movable(FoF.Group);
119
120 /* save previous data */
121 PrevTotNsubhalos = FoF.TotNsubhalos;
122 PrevNsubhalos = FoF.Nsubhalos;
123
124 PrevMtrP_NumPart = MtrP_NumPart;
125
126 PrevMtrP =
127 (mergertree_particle_data *)Mem.mymalloc_movable(&PrevMtrP, "PrevMtrP", PrevMtrP_NumPart * sizeof(mergertree_particle_data));
128 memcpy(PrevMtrP, MtrP, PrevMtrP_NumPart * sizeof(mergertree_particle_data));
129 Mem.myfree_movable(MtrP);
130
131 /* load new snapshot and group/subhalo catalogues */
132 FoF_IO.fof_subfind_load_groups(num);
133 CurrTotNsubhalos = FoF.TotNsubhalos;
134 CurrNsubhalos = FoF.Nsubhalos;
135
136 Snap_IO.mergertree_read_snap_ids(num);
137
138 /* make sure that group catalog and snapshot are in file order */
139 mycxxsort_parallel(MtrP, MtrP + MtrP_NumPart, compare_MtrP_FileOffset, Communicator);
140 mycxxsort_parallel(FoF.Group, FoF.Group + FoF.Ngroups, compare_Group_FileOffset, Communicator);
141 mycxxsort_parallel(FoF.Subhalo, FoF.Subhalo + FoF.Nsubhalos, compare_Subhalo_FileOffset, Communicator);
142
143 if(FoF_IO.LegacyFormat)
144 {
145 mpi_printf("\nFOF/SUBFIND: Legacy format from Arepo detected, trying to adjust for this.\n");
146 FoF.subfind_redetermine_groupnr();
147 }
148
149 /* now assign (global) group and subhalo numbers to each particle belonging to the particular structure */
150 mergertree_assign_group_numbers(&FoF);
151
152 Mem.myfree_movable(FoF.Subhalo);
153 Mem.myfree_movable(FoF.Group);
154
155 /* assign the determined new subhalonr */
156 for(int i = 0; i < MtrP_NumPart; i++)
157 {
158 MtrP[i].SubhaloNr = MtrP[i].PrevSubhaloNr;
159 MtrP[i].SubhaloLen = MtrP[i].PrevSubhaloLen;
160 MtrP[i].GroupNr = MtrP[i].PrevGroupNr;
161 MtrP[i].RankInSubhalo = MtrP[i].PrevRankInSubhalo;
162 MtrP[i].PrevSubhaloNr = HALONR_MAX;
163 MtrP[i].PrevGroupNr = HALONR_MAX;
164 MtrP[i].PrevSubhaloLen = 0;
165 MtrP[i].PrevRankInSubhalo = 0;
166 }
167
168 /* get PrevSubhaloNr from PrevMtrP by matching IDs */
169 mergertree_match_ids_of_previous_snap();
170
171 mergertree_determine_descendants_postproc(num);
172}
173
174void mergertree::mergertree_match_ids_of_previous_snap(void)
175{
176 int *Send_count = (int *)Mem.mymalloc("Send_count", sizeof(int) * NTask);
177 int *Send_offset = (int *)Mem.mymalloc("Send_offset", sizeof(int) * NTask);
178 int *Recv_count = (int *)Mem.mymalloc("Recv_count", sizeof(int) * NTask);
179 int *Recv_offset = (int *)Mem.mymalloc("Recv_offset", sizeof(int) * NTask);
180
181 mycxxsort_parallel(MtrP, MtrP + MtrP_NumPart, compare_MtrP_ID, Communicator);
182 mycxxsort_parallel(PrevMtrP, PrevMtrP + PrevMtrP_NumPart, compare_MtrP_ID, Communicator);
183
184 MyIDType *list_min_id = (MyIDType *)Mem.mymalloc("list_min_id", NTask * sizeof(MyIDType));
185 MyIDType *list_max_id = (MyIDType *)Mem.mymalloc("list_max_id", NTask * sizeof(MyIDType));
186 int *list_numpart = (int *)Mem.mymalloc("list_mumpart", NTask * sizeof(int));
187
188 MPI_Allgather(&MtrP[0].ID, sizeof(MyIDType), MPI_BYTE, list_min_id, sizeof(MyIDType), MPI_BYTE, Communicator);
189 MPI_Allgather(&MtrP[MtrP_NumPart > 0 ? MtrP_NumPart - 1 : 0].ID, sizeof(MyIDType), MPI_BYTE, list_max_id, sizeof(MyIDType), MPI_BYTE,
190 Communicator);
191 MPI_Allgather(&MtrP_NumPart, sizeof(int), MPI_BYTE, list_numpart, sizeof(int), MPI_BYTE, Communicator);
192
193 int nexport = 0, nimport = 0;
194 mergertree_particle_data *import_data = NULL, *export_data = NULL;
195
196 for(int mode = 0; mode < 2; mode++)
197 {
198 for(int i = 0; i < NTask; i++)
199 Send_count[i] = 0;
200
201 int target = 0;
202
203 for(int i = 0; i < PrevMtrP_NumPart; i++)
204 {
205 while(target < NTask - 1 && (list_numpart[target] == 0 || PrevMtrP[i].ID > list_max_id[target]))
206 target++;
207
208 if(list_numpart[target] != 0)
209 if(PrevMtrP[i].ID >= list_min_id[target] && PrevMtrP[i].ID <= list_max_id[target])
210 {
211 if(mode == 0)
212 Send_count[target]++;
213 else
214 {
215 int off = Send_offset[target] + Send_count[target]++;
216 export_data[off] = PrevMtrP[i];
217 }
218 }
219 }
220
221 if(mode == 0)
222 {
223 myMPI_Alltoall(Send_count, 1, MPI_INT, Recv_count, 1, MPI_INT, Communicator);
224 Recv_offset[0] = Send_offset[0] = 0;
225 for(int j = 0; j < NTask; j++)
226 {
227 nimport += Recv_count[j];
228 nexport += Send_count[j];
229 if(j > 0)
230 {
231 Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];
232 Recv_offset[j] = Recv_offset[j - 1] + Recv_count[j - 1];
233 }
234 }
235
236 export_data = (mergertree_particle_data *)Mem.mymalloc("export_data", nexport * sizeof(mergertree_particle_data));
237 import_data = (mergertree_particle_data *)Mem.mymalloc("import_data", nimport * sizeof(mergertree_particle_data));
238 }
239 }
240
241 for(int ngrp = 0; ngrp < (1 << PTask); ngrp++) /* note: here we also have a transfer from each task to itself (for ngrp=0) */
242 {
243 int recvTask = ThisTask ^ ngrp;
244 if(recvTask < NTask)
245 if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
246 myMPI_Sendrecv(&export_data[Send_offset[recvTask]], Send_count[recvTask] * sizeof(mergertree_particle_data), MPI_BYTE,
247 recvTask, TAG_DENS_B, &import_data[Recv_offset[recvTask]],
248 Recv_count[recvTask] * sizeof(mergertree_particle_data), MPI_BYTE, recvTask, TAG_DENS_B, Communicator,
249 MPI_STATUS_IGNORE);
250 }
251
252 /* incoming data should already be sorted, so now do the match */
253
254 for(int i = 0, j = 0; i < MtrP_NumPart && j < nimport;)
255 {
256 if(MtrP[i].ID < import_data[j].ID)
257 i++;
258 else if(MtrP[i].ID > import_data[j].ID)
259 j++;
260 else
261 {
262 MtrP[i].PrevSubhaloNr = import_data[j].PrevSubhaloNr;
263 MtrP[i].PrevSubhaloLen = import_data[j].PrevSubhaloLen;
264 MtrP[i].PrevGroupNr = import_data[j].PrevGroupNr;
265 MtrP[i].PrevRankInSubhalo = import_data[j].PrevRankInSubhalo;
266 i++;
267 j++;
268 }
269 }
270
271 Mem.myfree(import_data);
272 Mem.myfree(export_data);
273 Mem.myfree(list_numpart);
274 Mem.myfree(list_max_id);
275 Mem.myfree(list_min_id);
276
277 Mem.myfree(Recv_offset);
278 Mem.myfree(Recv_count);
279 Mem.myfree(Send_offset);
280 Mem.myfree(Send_count);
281}
282
283void mergertree::mergertree_assign_group_numbers(fof<simparticles> *FoF)
284{
285 int *Send_count = (int *)Mem.mymalloc("Send_count", sizeof(int) * NTask);
286 int *Send_offset = (int *)Mem.mymalloc("Send_offset", sizeof(int) * NTask);
287 int *Recv_count = (int *)Mem.mymalloc("Recv_count", sizeof(int) * NTask);
288 int *Recv_offset = (int *)Mem.mymalloc("Recv_offset", sizeof(int) * NTask);
289
290 /* Tell everybody, how many particles are stored by each processor */
291
292 FoF->fof_assign_group_offset();
293 FoF->subfind_assign_subhalo_offsettype();
294
295 int ntype_loc[NTYPES]; /* local particle number of each type */
296 long long ntype_tot[NTYPES]; /* total particle numbers of each type */
297 long long ntype_previous[NTYPES]; /* cumulative number of particles of each type on previous processors */
298
299 for(int i = 0; i < NTYPES; i++)
300 ntype_loc[i] = 0;
301
302 for(int i = 0; i < MtrP_NumPart; i++)
303 ntype_loc[MtrP[i].Type]++;
304
305 /* collect a table with the particle numbers of each type on each processors */
306 int *ntype_all = (int *)Mem.mymalloc("ntype_all", NTYPES * NTask * sizeof(int));
307 MPI_Allgather(ntype_loc, NTYPES, MPI_INT, ntype_all, NTYPES, MPI_INT, Communicator);
308
309 for(int i = 0; i < NTYPES; i++)
310 ntype_tot[i] = 0;
311
312 for(int i = 0; i < NTask; i++)
313 for(int j = 0; j < NTYPES; j++)
314 ntype_tot[j] += ntype_all[i * NTYPES + j];
315
316 for(int i = 0; i < NTYPES; i++)
317 ntype_previous[i] = 0;
318
319 for(int i = 0; i < ThisTask; i++)
320 for(int j = 0; j < NTYPES; j++)
321 ntype_previous[j] += ntype_all[i * NTYPES + j];
322
323 /* tell everybody how many groups each processor holds */
324
325 int *gcount = (int *)Mem.mymalloc("gcount", NTask * sizeof(int));
326 MPI_Allgather(&FoF->Ngroups, 1, MPI_INT, gcount, 1, MPI_INT, Communicator);
327
328 /* determine the number of the first group we hold */
329 long long first_groupnr = 0;
330 for(int i = 0; i < ThisTask; i++)
331 first_groupnr += gcount[i];
332
333 /* tell everybody how many subhalos each processor holds */
334
335 int *scount = (int *)Mem.mymalloc("scount", NTask * sizeof(int));
336 MPI_Allgather(&FoF->Nsubhalos, 1, MPI_INT, scount, 1, MPI_INT, Communicator);
337
338 /* determine the number of the first subhalo we hold */
339 long long first_subhalonr = 0;
340 for(int i = 0; i < ThisTask; i++)
341 first_subhalonr += scount[i];
342
343 /* let's now figure out which groups are needed by different processors to assign the group number information */
344
345 struct group_info
346 {
347 long long GroupNr;
348 long long First;
349 MyLenType Len;
350 };
351 group_info *export_group_data = NULL, *import_group_data = NULL;
352
353 for(int type = 0; type < NTYPES; type++)
354 {
355 int nexport = 0, nimport = 0;
356
357 for(int mode = 0; mode < 2; mode++)
358 {
359 for(int i = 0; i < NTask; i++)
360 Send_count[i] = 0;
361
362 int target = 0;
363 long long first_in_target = 0; /* this the first particle (of this type) on the current target processor */
364
365 for(int i = 0; i < FoF->Ngroups && target < NTask;)
366 {
367 int flag = 0;
368
369 /* check whether we have an overlap */
370 if(FoF->Group[i].OffsetType[type] + FoF->Group[i].LenType[type] > first_in_target &&
371 FoF->Group[i].OffsetType[type] < (first_in_target + ntype_all[target * NTYPES + type]))
372 {
373 flag = 1;
374
375 if(mode == 0)
376 Send_count[target]++;
377 else
378 {
379 int off = Send_offset[target] + Send_count[target]++;
380 export_group_data[off].GroupNr = first_groupnr + i;
381 export_group_data[off].First = FoF->Group[i].OffsetType[type];
382 export_group_data[off].Len = FoF->Group[i].LenType[type];
383 }
384 }
385
386 if(FoF->Group[i].OffsetType[type] + FoF->Group[i].LenType[type] > first_in_target + ntype_all[target * NTYPES + type])
387 {
388 first_in_target += ntype_all[target * NTYPES + type];
389 target++;
390 }
391 else
392 {
393 if(i < FoF->Ngroups && flag == 0 && FoF->Group[i].LenType[type] > 0)
394 {
395 Terminate(
396 "strange: type=%d mode=%d i=%d first_in_target=%lld target=%d FoF->Group[i].LenType[type]=%lld "
397 "Ngroups=%d",
398 type, mode, i, first_in_target, target, (long long)FoF->Group[i].LenType[type], FoF->Ngroups);
399 }
400
401 i++;
402 }
403 }
404
405 if(mode == 0)
406 {
407 myMPI_Alltoall(Send_count, 1, MPI_INT, Recv_count, 1, MPI_INT, Communicator);
408 Recv_offset[0] = Send_offset[0] = 0;
409 for(int j = 0; j < NTask; j++)
410 {
411 nimport += Recv_count[j];
412 nexport += Send_count[j];
413 if(j > 0)
414 {
415 Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];
416 Recv_offset[j] = Recv_offset[j - 1] + Recv_count[j - 1];
417 }
418 }
419
420 export_group_data = (group_info *)Mem.mymalloc("export_group_data", nexport * sizeof(group_info));
421 import_group_data = (group_info *)Mem.mymalloc("import_group_data", nimport * sizeof(group_info));
422 }
423 }
424
425 for(int ngrp = 0; ngrp < (1 << PTask); ngrp++) /* note: here we also have a transfer from each task to itself (for ngrp=0) */
426 {
427 int recvTask = ThisTask ^ ngrp;
428 if(recvTask < NTask)
429 if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
430 myMPI_Sendrecv(&export_group_data[Send_offset[recvTask]], Send_count[recvTask] * sizeof(group_info), MPI_BYTE, recvTask,
431 TAG_DENS_B, &import_group_data[Recv_offset[recvTask]], Recv_count[recvTask] * sizeof(group_info), MPI_BYTE,
432 recvTask, TAG_DENS_B, Communicator, MPI_STATUS_IGNORE);
433 }
434
435 /* now let's go through the local particles and assign the group numbers */
436
437 int p = ntype_previous[type];
438 int gr = 0;
439
440 for(int i = 0; i < MtrP_NumPart; i++)
441 {
442 if(MtrP[i].Type == type)
443 {
444 MtrP[i].PrevGroupNr = HALONR_MAX; /* default is not in a group */
445
446 while(gr < nimport && p > import_group_data[gr].First + import_group_data[gr].Len - 1)
447 gr++;
448
449 if(gr < nimport && p >= import_group_data[gr].First && p < import_group_data[gr].First + import_group_data[gr].Len)
450 MtrP[i].PrevGroupNr = import_group_data[gr].GroupNr;
451
452 p++;
453 }
454 }
455
456 Mem.myfree(import_group_data);
457 Mem.myfree(export_group_data);
458 }
459
460 /* let's now figure out which groups are needed by different processors to assign the group number information */
461
462 struct subhalo_info
463 {
464 long long SubhaloNr;
465 long long First;
466 MyLenType Len;
467 MyLenType SubhaloLen;
468 };
469 subhalo_info *export_subhalo_data = NULL, *import_subhalo_data = NULL;
470
471 for(int type = 0; type < NTYPES; type++)
472 {
473 int nexport = 0, nimport = 0;
474
475 for(int mode = 0; mode < 2; mode++)
476 {
477 for(int i = 0; i < NTask; i++)
478 Send_count[i] = 0;
479
480 int target = 0;
481 long long first_in_target = 0; /* this the first particle (of this type) on the current target processor */
482
483 for(int i = 0; i < FoF->Nsubhalos && target < NTask;)
484 {
485 /* check whether we have an overlap */
486 if(FoF->Subhalo[i].OffsetType[type] + FoF->Subhalo[i].LenType[type] > first_in_target &&
487 FoF->Subhalo[i].OffsetType[type] < (first_in_target + ntype_all[target * NTYPES + type]))
488 {
489 if(mode == 0)
490 Send_count[target]++;
491 else
492 {
493 int off = Send_offset[target] + Send_count[target]++;
494 export_subhalo_data[off].SubhaloNr = first_subhalonr + i;
495 export_subhalo_data[off].First = FoF->Subhalo[i].OffsetType[type];
496 export_subhalo_data[off].Len = FoF->Subhalo[i].LenType[type];
497 export_subhalo_data[off].SubhaloLen = FoF->Subhalo[i].Len;
498 }
499 }
500
501 if(FoF->Subhalo[i].OffsetType[type] + FoF->Subhalo[i].LenType[type] >
502 first_in_target + ntype_all[target * NTYPES + type])
503 {
504 first_in_target += ntype_all[target * NTYPES + type];
505 target++;
506 }
507 else
508 i++;
509 }
510
511 if(mode == 0)
512 {
513 myMPI_Alltoall(Send_count, 1, MPI_INT, Recv_count, 1, MPI_INT, Communicator);
514 Recv_offset[0] = Send_offset[0] = 0;
515 for(int j = 0; j < NTask; j++)
516 {
517 nimport += Recv_count[j];
518 nexport += Send_count[j];
519 if(j > 0)
520 {
521 Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];
522 Recv_offset[j] = Recv_offset[j - 1] + Recv_count[j - 1];
523 }
524 }
525
526 export_subhalo_data = (subhalo_info *)Mem.mymalloc("export_subhalo_data", nexport * sizeof(subhalo_info));
527 import_subhalo_data = (subhalo_info *)Mem.mymalloc("import_subhalo_data", nimport * sizeof(subhalo_info));
528 }
529 }
530
531 for(int ngrp = 0; ngrp < (1 << PTask); ngrp++) /* note: here we also have a transfer from each task to itself (for ngrp=0) */
532 {
533 int recvTask = ThisTask ^ ngrp;
534 if(recvTask < NTask)
535 if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
536 myMPI_Sendrecv(&export_subhalo_data[Send_offset[recvTask]], Send_count[recvTask] * sizeof(subhalo_info), MPI_BYTE,
537 recvTask, TAG_DENS_B, &import_subhalo_data[Recv_offset[recvTask]],
538 Recv_count[recvTask] * sizeof(subhalo_info), MPI_BYTE, recvTask, TAG_DENS_B, Communicator,
539 MPI_STATUS_IGNORE);
540 }
541
542 /* now let's go through the local particles and assign group numbers */
543
544 int p = ntype_previous[type];
545 int gr = 0;
546
547 for(int i = 0; i < MtrP_NumPart; i++)
548 {
549 if(MtrP[i].Type == type)
550 {
551 MtrP[i].PrevSubhaloNr = HALONR_MAX; /* default is not in a group */
552 MtrP[i].PrevSubhaloLen = 0;
553
554 while(gr < nimport && p > import_subhalo_data[gr].First + import_subhalo_data[gr].Len - 1)
555 gr++;
556
557 if(gr < nimport && p >= import_subhalo_data[gr].First && p < import_subhalo_data[gr].First + import_subhalo_data[gr].Len)
558 {
559 MtrP[i].PrevSubhaloNr = import_subhalo_data[gr].SubhaloNr;
560 MtrP[i].PrevSubhaloLen = import_subhalo_data[gr].SubhaloLen;
561
562 int rank = p - import_subhalo_data[gr].First; // Note: This is the rank within particles of the same type
563
564 if(rank > UCHAR_MAX) // restrict this to 1 byte (which is the storage we have reserved for this)
565 rank = UCHAR_MAX;
566
567 MtrP[i].PrevRankInSubhalo = static_cast<unsigned short>(rank);
568 }
569 p++;
570 }
571 }
572
573 Mem.myfree(import_subhalo_data);
574 Mem.myfree(export_subhalo_data);
575 }
576
577 Mem.myfree(scount);
578 Mem.myfree(gcount);
579 Mem.myfree(ntype_all);
580
581 Mem.myfree(Recv_offset);
582 Mem.myfree(Recv_count);
583 Mem.myfree(Send_offset);
584 Mem.myfree(Send_count);
585}
586
587void mergertree::mergertree_match_ids_of_current_snap(void)
588{
589 int *Send_count = (int *)Mem.mymalloc("Send_count", sizeof(int) * NTask);
590 int *Send_offset = (int *)Mem.mymalloc("Send_offset", sizeof(int) * NTask);
591 int *Recv_count = (int *)Mem.mymalloc("Recv_count", sizeof(int) * NTask);
592 int *Recv_offset = (int *)Mem.mymalloc("Recv_offset", sizeof(int) * NTask);
593
594 assign_particle_data *AsP = (assign_particle_data *)Mem.mymalloc("AsP", sizeof(assign_particle_data) * (Sp->NumPart + 1));
595
596 for(int i = 0; i < Sp->NumPart; i++)
597 {
598 AsP[i].OriginTask = ThisTask;
599 AsP[i].OriginIndex = i;
600 AsP[i].ID = Sp->P[i].ID.get();
601 }
602
603 mycxxsort_parallel(AsP, AsP + Sp->NumPart, compare_AssignP_ID, Communicator);
604 mycxxsort_parallel(MtrP, MtrP + MtrP_NumPart, compare_MtrP_ID, Communicator);
605
606 MyIDType *list_min_id = (MyIDType *)Mem.mymalloc("list_min_id", NTask * sizeof(MyIDType));
607 MyIDType *list_max_id = (MyIDType *)Mem.mymalloc("list_max_id", NTask * sizeof(MyIDType));
608 int *list_numpart = (int *)Mem.mymalloc("list_mumpart", NTask * sizeof(int));
609
610 MPI_Allgather(&AsP[0].ID, sizeof(MyIDType), MPI_BYTE, list_min_id, sizeof(MyIDType), MPI_BYTE, Communicator);
611 MPI_Allgather(&AsP[Sp->NumPart > 0 ? Sp->NumPart - 1 : 0].ID, sizeof(MyIDType), MPI_BYTE, list_max_id, sizeof(MyIDType), MPI_BYTE,
612 Communicator);
613 MPI_Allgather(&Sp->NumPart, sizeof(int), MPI_BYTE, list_numpart, sizeof(int), MPI_BYTE, Communicator);
614
615 int nexport = 0, nimport = 0;
616 mergertree_particle_data *import_data = NULL, *export_data = NULL;
617
618 for(int mode = 0; mode < 2; mode++)
619 {
620 for(int i = 0; i < NTask; i++)
621 Send_count[i] = 0;
622
623 int target = 0;
624
625 for(int i = 0; i < MtrP_NumPart; i++)
626 {
627 while(target < NTask - 1 && (list_numpart[target] == 0 || MtrP[i].ID > list_max_id[target]))
628 target++;
629
630 if(list_numpart[target] != 0)
631 if(MtrP[i].ID >= list_min_id[target] && MtrP[i].ID <= list_max_id[target])
632 {
633 if(mode == 0)
634 Send_count[target]++;
635 else
636 {
637 int off = Send_offset[target] + Send_count[target]++;
638 export_data[off] = MtrP[i];
639 }
640 }
641 }
642
643 if(mode == 0)
644 {
645 myMPI_Alltoall(Send_count, 1, MPI_INT, Recv_count, 1, MPI_INT, Communicator);
646 Recv_offset[0] = Send_offset[0] = 0;
647 for(int j = 0; j < NTask; j++)
648 {
649 nimport += Recv_count[j];
650 nexport += Send_count[j];
651 if(j > 0)
652 {
653 Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];
654 Recv_offset[j] = Recv_offset[j - 1] + Recv_count[j - 1];
655 }
656 }
657
658 export_data = (mergertree_particle_data *)Mem.mymalloc("export_data", nexport * sizeof(mergertree_particle_data));
659 import_data = (mergertree_particle_data *)Mem.mymalloc("import_data", nimport * sizeof(mergertree_particle_data));
660 }
661 }
662
663 for(int ngrp = 0; ngrp < (1 << PTask); ngrp++) /* note: here we also have a transfer from each task to itself (for ngrp=0) */
664 {
665 int recvTask = ThisTask ^ ngrp;
666 if(recvTask < NTask)
667 if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
668 myMPI_Sendrecv(&export_data[Send_offset[recvTask]], Send_count[recvTask] * sizeof(mergertree_particle_data), MPI_BYTE,
669 recvTask, TAG_DENS_B, &import_data[Recv_offset[recvTask]],
670 Recv_count[recvTask] * sizeof(mergertree_particle_data), MPI_BYTE, recvTask, TAG_DENS_B, Communicator,
671 MPI_STATUS_IGNORE);
672 }
673
674 /* incoming data should already be sorted, so now do the match */
675
676 for(int i = 0, j = 0; i < Sp->NumPart && j < nimport;)
677 {
678 if(AsP[i].ID < import_data[j].ID)
679 i++;
680 else if(AsP[i].ID > import_data[j].ID)
681 j++;
682 else
683 {
684 AsP[i].PrevSubhaloNr = import_data[j].PrevSubhaloNr;
685 AsP[i].PrevSubhaloLen = import_data[j].PrevSubhaloLen;
686 i++;
687 j++;
688 }
689 }
690
691 mycxxsort_parallel(AsP, AsP + Sp->NumPart, compare_AssignP_Origin, Communicator);
692
693 for(int i = 0; i < Sp->NumPart; i++)
694 {
695 Sp->P[i].PrevSizeOfSubhalo.set(AsP[i].PrevSubhaloLen);
696 Sp->P[i].PrevSubhaloNr.set(AsP[i].PrevSubhaloNr);
697 }
698
699 Mem.myfree(import_data);
700 Mem.myfree(export_data);
701 Mem.myfree(list_numpart);
702 Mem.myfree(list_max_id);
703 Mem.myfree(list_min_id);
704
705 Mem.myfree(AsP);
706
707 Mem.myfree(Recv_offset);
708 Mem.myfree(Recv_count);
709 Mem.myfree(Send_offset);
710 Mem.myfree(Send_count);
711}
712
713#endif
All
global_data_all_processes All
Definition: main.cc:40
domain
Definition: domain.h:31
fof_io
Definition: fof_io.h:20
NTYPES
#define NTYPES
Definition: constants.h:308
MyLenType
int MyLenType
Definition: dtypes.h:76
MyIDType
unsigned int MyIDType
Definition: dtypes.h:68
HALONR_MAX
#define HALONR_MAX
Definition: idstorage.h:20
Terminate
#define Terminate(...)
Definition: macros.h:15
myMPI_Alltoall
int myMPI_Alltoall(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, MPI_Comm comm)
Definition: myalltoall.cc:301
myMPI_Sendrecv
int myMPI_Sendrecv(void *sendbuf, size_t sendcount, MPI_Datatype sendtype, int dest, int sendtag, void *recvbuf, size_t recvcount, MPI_Datatype recvtype, int source, int recvtag, MPI_Comm comm, MPI_Status *status)
Definition: sizelimited_sendrecv.cc:25
TAG_DENS_B
#define TAG_DENS_B
Definition: mpi_utils.h:51
Mem
memory Mem
Definition: main.cc:44
mycxxsort_parallel
double mycxxsort_parallel(T *begin, T *end, Comp comp, MPI_Comm comm)
Definition: parallel_sort.h:162
global_data_all_processes::SnapFormat
int SnapFormat
Definition: allvars.h:46