12#include "gadgetconfig.h"
16#include <gsl/gsl_rng.h>
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"
33#include "../logs/timer.h"
34#include "../main/simulation.h"
35#include "../mergertree/io_readtrees_mbound.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"
51readtrees_mbound_io::readtrees_mbound_io(mergertree *MergerTree_ptr, MPI_Comm comm,
int format) :
IO_Def(comm, format)
53 MergerTree = MergerTree_ptr;
55 this->N_IO_Fields = 0;
56 this->N_DataGroups = 2;
57 this->header_size =
sizeof(header);
58 this->header_buf = &header;
60 sprintf(this->info,
"MERGERTREE: reading/writing mergertrees");
62 init_field(
"MTRL",
"Length",
MEM_INT,
FILE_INT,
READ_IF_PRESENT, 1,
A_TT, &MergerTree->TreeTable[0].HaloCount, NULL,
TREELENGTH, 0,
66 init_field(
"MTRI",
"TreeID",
MEM_INT64,
FILE_INT64,
READ_IF_PRESENT, 1,
A_TT, &MergerTree->TreeTable[0].TreeID, NULL,
TREELENGTH, 0,
71 init_field(
"TRID",
"TreeID",
MEM_INT64,
FILE_INT64,
READ_IF_PRESENT, 1,
A_TID, &MergerTree->HaloIDdata[0].TreeID, NULL,
TREEHALOS, 0,
75 &MergerTree->HaloIDdata[0].SubMostBoundID, NULL,
TREEHALOS, 0, 0, 0, 0, 0, 0, 0);
78void readtrees_mbound_io::read_trees_mostbound(
void)
82 MergerTree->TotNtrees = 0;
83 MergerTree->TotNhalos = 0;
88 sprintf(fname,
"%s/treedata/%s",
All.
OutputDir,
"trees");
92 int num_files = find_files(fname, fname);
94 reset_io_byte_count();
96 for(
int rep = 0; rep < 2; rep++)
98 MergerTree->Ntrees = 0;
99 MergerTree->Nhalos = 0;
101 read_files_driver(fname, rep, num_files);
106 MergerTree->TreeTable = (
halotrees_table *)
Mem.mymalloc_movable(&MergerTree->TreeTable,
"TreeTable",
108 MergerTree->HaloIDdata = (treehalo_ids_type *)
Mem.mymalloc_movable(&MergerTree->HaloIDdata,
"HaloIDdata",
109 (MergerTree->Nhalos + 1) *
sizeof(treehalo_ids_type));
113 MPI_Barrier(Communicator);
115 long long byte_count = get_io_byte_count(), byte_count_all;
116 sumup_longs(1, &byte_count, &byte_count_all, Communicator);
120 mpi_printf(
"MERGERTREE-READ: reading done. Took %g sec, total size %g MB, corresponds to effective I/O rate of %g MB/sec\n",
121 Logs.
timediff(t0, t1), byte_count_all / (1024.0 * 1024.0), byte_count_all / (1024.0 * 1024.0) /
Logs.
timediff(t0, t1));
123 mpi_printf(
"\nMERGERTREE-READ: Total number of trees=%ldd, total number of halos=%lld\n\n", MergerTree->TotNtrees,
124 MergerTree->TotNhalos);
126 MPI_Barrier(Communicator);
128 for(
int i = 0; i < MergerTree->Ntrees; i++)
130 if(MergerTree->HaloIDdata[i].TreeID > MergerTree->TotNtrees)
131 Terminate(
"i=%d MergerTree->Ntrees=%d MergerTree->HaloIDdata[i].TreeID=%lld MergerTree->TotNtrees=%lld ", i,
132 MergerTree->Ntrees, MergerTree->HaloIDdata[i].TreeID, MergerTree->TotNtrees);
136void readtrees_mbound_io::fill_file_header(
int writeTask,
int lastTask,
long long *n_type,
long long *ntot_type) {}
138void readtrees_mbound_io::write_header_fields(hid_t handle) {}
140int readtrees_mbound_io::get_filenr_from_header(
void) {
return header.num_files; }
142void readtrees_mbound_io::set_filenr_in_header(
int numfiles) { header.num_files = numfiles; }
144void readtrees_mbound_io::read_increase_numbers(
int type,
int n_for_this_task)
149 MergerTree->Ntrees += n_for_this_task;
152 MergerTree->Nhalos += n_for_this_task;
161void readtrees_mbound_io::read_file_header(
const char *fname,
int filenr,
int readTask,
int lastTask,
long long *n_type,
162 long long *ntot_type,
int *nstart)
164 if(ThisTask == readTask)
166 if(filenr == 0 && nstart == NULL)
169 "\nREAD-TREES: filenr=%d, '%s' contains %lld trees out of a total of %lld, and %lld halos out of a total of %lld\n",
170 filenr, fname, header.Ntrees, header.TotNtrees, header.Nhalos, header.TotNhalos);
174 if(MergerTree->TotNtrees == 0)
176 MergerTree->TotNtrees = header.TotNtrees;
177 MergerTree->TotNhalos = header.TotNhalos;
180 for(
int k = 0; k < 2; k++)
181 n_type[k] = ntot_type[k] = 0;
184 ntot_type[0] = header.Ntrees;
185 long long n_in_file = header.Ntrees;
186 int ntask = lastTask - readTask + 1;
187 int n_for_this_task = n_in_file / ntask;
188 if((ThisTask - readTask) < (n_in_file % ntask))
190 n_type[0] = n_for_this_task;
193 memmove(&MergerTree->TreeTable[n_for_this_task], &MergerTree->TreeTable[0], MergerTree->Ntrees *
sizeof(
halotrees_table));
197 ntot_type[1] = header.Nhalos;
198 long long n_in_file = header.Nhalos;
199 int ntask = lastTask - readTask + 1;
200 int n_for_this_task = n_in_file / ntask;
201 if((ThisTask - readTask) < (n_in_file % ntask))
203 n_type[1] = n_for_this_task;
206 memmove(&MergerTree->HaloIDdata[n_for_this_task], &MergerTree->HaloIDdata[0], MergerTree->Nhalos *
sizeof(treehalo_ids_type));
213void readtrees_mbound_io::read_header_fields(
const char *fname)
215 memset(&header, 0,
sizeof(header));
217 hid_t hdf5_file =
my_H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
218 hid_t handle =
my_H5Gopen(hdf5_file,
"/Header");
234void readtrees_mbound_io::get_datagroup_name(
int type,
char *buf)
239 sprintf(buf,
"/TreeTable");
242 sprintf(buf,
"/TreeHalos");
250int readtrees_mbound_io::get_type_of_element(
int index)
256void readtrees_mbound_io::set_type_of_element(
int index,
int type)
260void *readtrees_mbound_io::get_base_address_of_structure(
enum arrays array,
int index)
265 return (
void *)(MergerTree->TreeTable + index);
267 return (
void *)(MergerTree->HaloIDdata + index);
269 Terminate(
"strange, we don't expect to get here");
global_data_all_processes All
double timediff(double t0, double t1)
#define MAXLEN_PATH_EXTRA
hid_t my_H5Gopen(hid_t loc_id, const char *groupname)
hid_t my_H5Fopen(const char *fname, unsigned int flags, hid_t fapl_id)
herr_t my_H5Gclose(hid_t group_id, const char *groupname)
herr_t my_H5Fclose(hid_t file_id, const char *fname)
void read_scalar_attribute(hid_t handle, const char *attr_name, void *buf, hid_t mem_type_id)
void sumup_longs(int n, long long *src, long long *res, MPI_Comm comm)
char OutputDir[MAXLEN_PATH]