GADGET-4
io_descendant.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 "../io/hdf5_util.h"
32#include "../io/io.h"
33#include "../logs/timer.h"
34#include "../main/simulation.h"
35#include "../mergertree/io_descendant.h"
36#include "../mergertree/mergertree.h"
37#include "../mpi_utils/mpi_utils.h"
38#include "../sort/parallel_sort.h"
39#include "../subfind/subfind.h"
40#include "../system/system.h"
41
42descendant_io::descendant_io(mergertree *MergerTree_ptr, MPI_Comm comm, int format) : IO_Def(comm, format)
43{
44 MergerTree = MergerTree_ptr;
45
46 this->N_IO_Fields = 0;
47 this->N_DataGroups = 1;
48 this->header_size = sizeof(header);
49 this->header_buf = &header;
50 this->type_of_file = FILE_IS_DESCCAT;
51 sprintf(this->info, "MERGERTREE: writing descendant information");
52
53 init_field("DSNR", "DescSubhaloNr", MEM_INT64, FILE_INT64, READ_IF_PRESENT, 1, A_DESC, NULL, io_func_descsubhalonr, PREVSUBS, 0, 0,
54 0, 0, 0, 0, 0, true);
55
56 init_field("FDNR", "FirstDescSubhaloNr", MEM_INT64, FILE_INT64, READ_IF_PRESENT, 1, A_DESC, NULL, io_func_firstdescsubhalonr,
57 PREVSUBS, 0, 0, 0, 0, 0, 0, 0, true);
58
59 init_field("NSNR", "NextProgSubhaloNr", MEM_INT64, FILE_INT64, READ_IF_PRESENT, 1, A_DESC, NULL, io_func_nextsubhalonr, PREVSUBS, 0,
60 0, 0, 0, 0, 0, 0, true);
61
62 init_field("SHNR", "SubhaloNr", MEM_INT64, FILE_INT64, READ_IF_PRESENT, 1, A_DESC, &MergerTree->Descendants[0].PrevSubhaloNr, NULL,
63 PREVSUBS, 0, 0, 0, 0, 0, 0, 0, true); // this field can in principle be deleted
64
65 init_field("DFLO", "DescFileOffset", MEM_MY_FILEOFFSET, FILE_NONE, SKIP_ON_READ, 1, A_DESC, &MergerTree->Descendants[0].FileOffset,
66 NULL, PREVSUBS, 0, 0, 0, 0, 0, 0, 0, true);
67}
68
69void descendant_io::mergertree_save_descendants(int num)
70{
71 char buf[MAXLEN_PATH_EXTRA];
72
73 /* write Descendants and Nextsubhalos */
74
76 {
77 if(ThisTask == 0)
78 {
79 sprintf(buf, "%s/groups_%03d", All.OutputDir, num - 1);
80 mkdir(buf, 02755);
81 }
82 MPI_Barrier(Communicator);
83 }
84
86 sprintf(buf, "%s/groups_%03d/%s_%03d", All.OutputDir, num - 1, "subhalo_desc", num - 1);
87 else
88 sprintf(buf, "%s%s_%03d", All.OutputDir, "subhalo_desc", num - 1);
89
90 write_multiple_files(buf, All.NumFilesPerSnapshot);
91}
92
93void descendant_io::mergertree_read_descendants(int num)
94{
95 char fname[MAXLEN_PATH_EXTRA], fname_multiple[MAXLEN_PATH_EXTRA];
96
97 sprintf(fname_multiple, "%s/groups_%03d/%s_%03d", All.OutputDir, num, "subhalo_desc", num);
98 sprintf(fname, "%s%s_%03d", All.OutputDir, "subhalo_desc", num);
99
100 TotNsubhalos = 0;
101
102 int num_files = find_files(fname, fname_multiple);
103
104 if(num_files > 1)
105 strcpy(fname, fname_multiple);
106
107 /* we repeat reading the headers of the files two times. In the first iteration, only the
108 * particle numbers ending up on each processor are assembled, followed by memory allocation.
109 * In the second iteration, the data is actually read in.
110 */
111 for(int rep = 0; rep < 2; rep++)
112 {
113 Nsubhalos = 0;
114
115 read_files_driver(fname, rep, num_files);
116
117 /* now do the memory allocation */
118 if(rep == 0)
119 {
120 MergerTree->Descendants = (mergertree::desc_list *)Mem.mymalloc_movable(&MergerTree->Descendants, "Descendants",
121 Nsubhalos * sizeof(mergertree::desc_list));
122 }
123 }
124
125 MPI_Barrier(Communicator);
126}
127
128void descendant_io::fill_file_header(int writeTask, int lastTask, long long *n_type, long long *ntot_type)
129{
130 /* determine group/id numbers of each type in file */
131
132 n_type[0] = MergerTree->PrevNsubhalos;
133
134 if(ThisTask == writeTask)
135 {
136 for(int n = 0; n < 1; n++)
137 ntot_type[n] = n_type[n];
138
139 for(int task = writeTask + 1; task <= lastTask; task++)
140 {
141 long long nn[3];
142 MPI_Recv(&nn[0], 1, MPI_LONG_LONG, task, TAG_LOCALN, Communicator, MPI_STATUS_IGNORE);
143 for(int n = 0; n < 1; n++)
144 ntot_type[n] += nn[n];
145 }
146
147 for(int task = writeTask + 1; task <= lastTask; task++)
148 MPI_Send(&ntot_type[0], 1, MPI_LONG_LONG, task, TAG_N, Communicator);
149 }
150 else
151 {
152 MPI_Send(&n_type[0], 1, MPI_LONG_LONG, writeTask, TAG_LOCALN, Communicator);
153 MPI_Recv(&ntot_type[0], 1, MPI_LONG_LONG, writeTask, TAG_N, Communicator, MPI_STATUS_IGNORE);
154 }
155
156 /* fill file header */
157
158 header.Nsubhalos = ntot_type[0];
159 header.TotNsubhalos = MergerTree->PrevTotNsubhalos;
160 header.num_files = All.NumFilesPerSnapshot;
161}
162
163void descendant_io::read_file_header(const char *fname, int filenr, int readTask, int lastTask, long long *n_type,
164 long long *ntot_type, int *nstart)
165{
166 if(ThisTask == readTask)
167 {
168 if(filenr == 0 && nstart == NULL)
169 {
170 mpi_printf("\nREAD-DESCENDANTS: filenr=%d, '%s' contains (subhalos): %8d\n", filenr, fname, header.Nsubhalos);
171 }
172 }
173
174 if(TotNsubhalos == 0)
175 TotNsubhalos = header.TotNsubhalos;
176
177 for(int k = 0; k < 1; k++)
178 n_type[k] = ntot_type[k] = 0;
179
180 /* to collect the gas particles all at the beginning (in case several
181 snapshot files are read on the current CPU) we move the collisionless
182 particles such that a gap of the right size is created */
183
184 ntot_type[0] = header.Nsubhalos;
185
186 long long n_in_file = header.Nsubhalos;
187 int ntask = lastTask - readTask + 1;
188 int n_for_this_task = n_in_file / ntask;
189 if((ThisTask - readTask) < (n_in_file % ntask))
190 n_for_this_task++;
191
192 n_type[0] = n_for_this_task;
193
194 if(nstart)
195 {
196 memmove(&MergerTree->Descendants[n_for_this_task], &MergerTree->Descendants[0], Nsubhalos * sizeof(mergertree::desc_list));
197 *nstart = 0;
198 }
199}
200
201void descendant_io::write_header_fields(hid_t handle)
202{
203 write_scalar_attribute(handle, "Nsubhalos_ThisFile", &header.Nsubhalos, H5T_NATIVE_UINT64);
204
205 write_scalar_attribute(handle, "Nsubhalos_Total", &header.TotNsubhalos, H5T_NATIVE_UINT64);
206
207 write_scalar_attribute(handle, "NumFiles", &header.num_files, H5T_NATIVE_INT);
208}
209
214void descendant_io::read_header_fields(const char *fname)
215{
216 memset(&header, 0, sizeof(io_header));
217
218 hid_t hdf5_file = my_H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
219 hid_t handle = my_H5Gopen(hdf5_file, "/Header");
220
221 /* now read the header fields */
222 read_scalar_attribute(handle, "Nsubhalos_ThisFile", "Nsubgroups_ThisFile", &header.Nsubhalos, H5T_NATIVE_UINT64);
223
224 read_scalar_attribute(handle, "Nsubhalos_Total", "Nsubgroups_Total", &header.TotNsubhalos, H5T_NATIVE_UINT64);
225
226 read_scalar_attribute(handle, "NumFiles", &header.num_files, H5T_NATIVE_INT);
227
228 my_H5Gclose(handle, "/Header");
229 my_H5Fclose(hdf5_file, fname);
230}
231
232int descendant_io::get_filenr_from_header(void) { return header.num_files; }
233
234void descendant_io::set_filenr_in_header(int numfiles) { header.num_files = numfiles; }
235
236void descendant_io::read_increase_numbers(int type, int n_for_this_task)
237{
238 switch(type)
239 {
240 case 0:
241 Nsubhalos += n_for_this_task;
242 break;
243 default:
244 Terminate("wrong group");
245 break;
246 }
247}
248
249void descendant_io::get_datagroup_name(int type, char *buf)
250{
251 switch(type)
252 {
253 case 0:
254 sprintf(buf, "/Subhalo");
255 break;
256 default:
257 Terminate("wrong group");
258 break;
259 }
260}
261
262int descendant_io::get_type_of_element(int index)
263{
264 /* empty */
265 return 0;
266}
267
268void descendant_io::set_type_of_element(int index, int type)
269{ /* empty */
270}
271
272void *descendant_io::get_base_address_of_structure(enum arrays array, int index)
273{
274 switch(array)
275 {
276 case A_DESC:
277 return (void *)(MergerTree->Descendants + index);
278 default:
279 Terminate("strange, we don't expect to get here");
280 }
281}
282
283#endif
global_data_all_processes All
Definition: main.cc:40
Definition: io.h:129
#define MAXLEN_PATH_EXTRA
Definition: constants.h:301
hid_t my_H5Gopen(hid_t loc_id, const char *groupname)
Definition: hdf5_util.cc:311
hid_t my_H5Fopen(const char *fname, unsigned int flags, hid_t fapl_id)
Definition: hdf5_util.cc:297
herr_t my_H5Gclose(hid_t group_id, const char *groupname)
Definition: hdf5_util.cc:457
herr_t my_H5Fclose(hid_t file_id, const char *fname)
Definition: hdf5_util.cc:482
void write_scalar_attribute(hid_t handle, const char *attr_name, const void *buf, hid_t mem_type_id)
Definition: hdf5_util.cc:621
void read_scalar_attribute(hid_t handle, const char *attr_name, void *buf, hid_t mem_type_id)
Definition: hdf5_util.cc:589
@ READ_IF_PRESENT
Definition: io.h:124
@ SKIP_ON_READ
Definition: io.h:125
@ PREVSUBS
Definition: io.h:110
arrays
Definition: io.h:30
@ A_DESC
Definition: io.h:38
@ FILE_NONE
Definition: io.h:66
@ FILE_INT64
Definition: io.h:68
@ MEM_MY_FILEOFFSET
Definition: io.h:87
@ MEM_INT64
Definition: io.h:80
@ FILE_IS_DESCCAT
Definition: io.h:55
#define Terminate(...)
Definition: macros.h:15
#define TAG_LOCALN
Definition: mpi_utils.h:52
#define TAG_N
Definition: mpi_utils.h:25
memory Mem
Definition: main.cc:44
char OutputDir[MAXLEN_PATH]
Definition: allvars.h:272