GADGET-4
io.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#include <errno.h>
15#include <hdf5.h>
16#include <math.h>
17#include <mpi.h>
18#include <stdio.h>
19#include <stdlib.h>
20#include <string.h>
21#include <sys/stat.h>
22#include <algorithm>
23
24#include "../cooling_sfr/cooling.h"
25#include "../data/allvars.h"
26#include "../data/dtypes.h"
27#include "../data/mymalloc.h"
28#include "../fof/fof.h"
29#include "../io/hdf5_util.h"
30#include "../io/io.h"
31#include "../io/parameters.h"
32#include "../lightcone/lightcone.h"
33#include "../logs/timer.h"
34#include "../main/main.h"
35#include "../main/simulation.h"
36#include "../mergertree/mergertree.h"
37#include "../mpi_utils/mpi_utils.h"
38#include "../subfind/subfind.h"
39#include "../system/system.h"
40
41#define HALF_ROUND_STYLE 1
42#include "../half/half.hpp"
43using half_float::half;
44
45/* local functions */
46
47IO_Def::~IO_Def()
48{
49 if(N_IO_Fields > 0)
50 Mem.myfree_movable(IO_Fields);
51}
52
66void IO_Def::init_field(const char *label, const char *datasetname, enum types_in_memory type_in_memory,
67 enum types_in_file type_in_file_output, enum read_flags read_flag, int values_per_block, enum arrays array,
68 void *pointer_to_field, void (*io_func)(IO_Def *, int, int, void *, int), int typelist_bitmask, int flagunit,
69 double a, double h, double L, double M, double V, double c, bool compression_on)
70{
71 const int alloc_step = 5;
72
73 if(values_per_block < 1) // if we have no values, we don't register the field
74 return;
75
76 if(Max_IO_Fields == 0)
77 {
78 IO_Fields = (IO_Field *)Mem.mymalloc_movable(&IO_Fields, "IO_Fields", alloc_step * sizeof(IO_Field));
79 Max_IO_Fields = alloc_step;
80 }
81 else if(Max_IO_Fields == N_IO_Fields)
82 {
83 Max_IO_Fields = ((Max_IO_Fields / alloc_step) + 1) * alloc_step;
84 IO_Fields = (IO_Field *)Mem.myrealloc_movable(IO_Fields, Max_IO_Fields * sizeof(IO_Field));
85 }
86
87 int n = N_IO_Fields++;
88 IO_Field *field = &IO_Fields[n];
89
90 strncpy(field->label, label, LABEL_LEN);
91 field->label[LABEL_LEN] = 0;
92
93 strncpy(field->datasetname, datasetname, DATASETNAME_LEN);
94 field->datasetname[DATASETNAME_LEN] = 0;
95
96 field->type_in_memory = type_in_memory;
97 field->type_in_file_output = type_in_file_output;
98 field->read_flag = read_flag;
99 field->values_per_block = values_per_block;
100 field->typelist = typelist_bitmask;
101#ifdef ALLOW_HDF5_COMPRESSION
102 field->compression_on = compression_on;
103#else
104 field->compression_on = false;
105#endif
106
107 field->array = array;
108 field->io_func = io_func;
109
110 if(array == A_NONE)
111 {
112 field->offset = 0;
113 }
114 else
115 {
116 field->offset = (size_t)pointer_to_field - (size_t)get_base_address_of_structure(array, 0);
117 }
118
119 field->hasunit = flagunit;
120 field->a = a;
121 field->h = h;
122 field->L = L;
123 field->M = M;
124 field->V = V;
125 field->c = c;
126}
127
132int IO_Def::find_files(const char *fname, const char *fname_multiple)
133{
134 FILE *fd;
135 char buf[200], buf1[200];
136 int dummy, files_found = 0;
137
138 if(file_format == FILEFORMAT_HDF5)
139 {
140 sprintf(buf, "%s.%d.hdf5", fname_multiple, 0);
141 sprintf(buf1, "%s.hdf5", fname);
142 }
143 else
144 {
145 sprintf(buf, "%s.%d", fname_multiple, 0);
146 sprintf(buf1, "%s", fname);
147 }
148
149 memset(header_buf, 0, header_size);
150
151 if(ThisTask == 0)
152 {
153 if((fd = fopen(buf, "r")))
154 {
155 if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
156 {
157 if(file_format == FILEFORMAT_LEGACY2)
158 {
159 my_fread(&dummy, sizeof(dummy), 1, fd);
160 my_fread(&dummy, sizeof(dummy), 1, fd);
161 my_fread(&dummy, sizeof(dummy), 1, fd);
162 my_fread(&dummy, sizeof(dummy), 1, fd);
163 }
164
165 my_fread(&dummy, sizeof(dummy), 1, fd);
166 my_fread(header_buf, header_size, 1, fd);
167 my_fread(&dummy, sizeof(dummy), 1, fd);
168 }
169 fclose(fd);
170
171 if(file_format == FILEFORMAT_HDF5)
172 read_header_fields(buf);
173
174 files_found = 1;
175 }
176 }
177
178 MPI_Bcast(header_buf, header_size, MPI_BYTE, 0, Communicator);
179 MPI_Bcast(&files_found, 1, MPI_INT, 0, Communicator);
180
181 if(get_filenr_from_header() > 0)
182 return get_filenr_from_header();
183
184 if(ThisTask == 0)
185 {
186 if((fd = fopen(buf1, "r")))
187 {
188 if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
189 {
190 if(file_format == FILEFORMAT_LEGACY2)
191 {
192 my_fread(&dummy, sizeof(dummy), 1, fd);
193 my_fread(&dummy, sizeof(dummy), 1, fd);
194 my_fread(&dummy, sizeof(dummy), 1, fd);
195 my_fread(&dummy, sizeof(dummy), 1, fd);
196 }
197
198 my_fread(&dummy, sizeof(dummy), 1, fd);
199 my_fread(header_buf, header_size, 1, fd);
200 my_fread(&dummy, sizeof(dummy), 1, fd);
201 }
202 fclose(fd);
203
204 if(file_format == FILEFORMAT_HDF5)
205 read_header_fields(buf1);
206
207 set_filenr_in_header(1);
208
209 files_found = 1;
210 }
211 }
212
213 MPI_Bcast(header_buf, header_size, MPI_BYTE, 0, Communicator);
214 MPI_Bcast(&files_found, 1, MPI_INT, 0, Communicator);
215
216 if(get_filenr_from_header() > 0)
217 return get_filenr_from_header();
218
219 if(files_found != 0)
220 Terminate("Have found IC files, but number of files in header seems to be zero\n");
221
222 Terminate("\nCan't find files, neither as '%s'\nnor as '%s'\n", buf, buf1);
223
224 return 0;
225}
226
227void IO_Def::read_files_driver(const char *fname, int rep, int num_files)
228{
229 if(rep == 0)
230 {
231 ntype_in_files =
232 (long long *)Mem.mymalloc_movable(&ntype_in_files, "ntype_in_files", num_files * N_DataGroups * sizeof(long long));
233 memset(ntype_in_files, 0, num_files * N_DataGroups * sizeof(long long));
234 }
235
236 void *CommBuffer = Mem.mymalloc("CommBuffer", COMMBUFFERSIZE);
237
238 int rest_files = num_files;
239
240 while(rest_files > NTask)
241 {
242 char buf[MAXLEN_PATH];
243
244 sprintf(buf, "%s.%d", fname, ThisTask + (rest_files - NTask));
245 if(file_format == FILEFORMAT_HDF5)
246 sprintf(buf, "%s.%d.hdf5", fname, ThisTask + (rest_files - NTask));
247
248 int ngroups = NTask / All.MaxFilesWithConcurrentIO;
249 if((NTask % All.MaxFilesWithConcurrentIO))
250 ngroups++;
251 int groupMaster = (ThisTask / ngroups) * ngroups;
252
253 for(int gr = 0; gr < ngroups; gr++)
254 {
255 if(ThisTask == (groupMaster + gr)) /* ok, it's this processor's turn */
256 {
257 if(rep == 0)
258 share_particle_number_in_file(buf, ThisTask + (rest_files - NTask), ThisTask, ThisTask);
259 else
260 read_file(buf, ThisTask + (rest_files - NTask), ThisTask, ThisTask, CommBuffer);
261 }
262 MPI_Barrier(Communicator);
263 }
264
265 rest_files -= NTask;
266 }
267
268 if(rest_files > 0)
269 {
270 int masterTask, filenr, lastTask;
271
272 distribute_file(rest_files, &filenr, &masterTask, &lastTask);
273
274 char buf[MAXLEN_PATH];
275
276 if(num_files > 1)
277 {
278 sprintf(buf, "%s.%d", fname, filenr);
279 if(file_format == FILEFORMAT_HDF5)
280 sprintf(buf, "%s.%d.hdf5", fname, filenr);
281 }
282 else
283 {
284 sprintf(buf, "%s", fname);
285 if(file_format == FILEFORMAT_HDF5)
286 sprintf(buf, "%s.hdf5", fname);
287 }
288
289 int ngroups = rest_files / All.MaxFilesWithConcurrentIO;
290 if((rest_files % All.MaxFilesWithConcurrentIO))
291 ngroups++;
292
293 for(int gr = 0; gr < ngroups; gr++)
294 {
295 if((filenr / All.MaxFilesWithConcurrentIO) == gr) /* ok, it's this processor's turn */
296 {
297 if(rep == 0)
298 share_particle_number_in_file(buf, filenr, masterTask, lastTask);
299 else
300 read_file(buf, filenr, masterTask, lastTask, CommBuffer);
301 }
302 MPI_Barrier(Communicator);
303 }
304 }
305
306 Mem.myfree(CommBuffer);
307
308 if(rep == 0)
309 MPI_Allreduce(MPI_IN_PLACE, ntype_in_files, num_files * N_DataGroups, MPI_LONG_LONG, MPI_MAX, Communicator);
310 else
311 {
312 /* we are done */
313 Mem.myfree_movable(ntype_in_files);
314 ntype_in_files = NULL;
315 }
316}
317
330void IO_Def::share_particle_number_in_file(const char *fname, int filenr, int readTask, int lastTask)
331{
332 long long n_type[N_DataGroups], npart[N_DataGroups];
333 unsigned int blksize1, blksize2;
334
335 if(ThisTask == readTask)
336 {
337 if(file_format == FILEFORMAT_HDF5)
338 {
339 read_header_fields(fname);
340 }
341 else if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
342 {
343 FILE *fd = 0;
344
345 if(!(fd = fopen(fname, "r")))
346 Terminate("can't open file `%s' for reading initial conditions.\n", fname);
347
348 if(file_format == FILEFORMAT_LEGACY2)
349 {
350 char label[LABEL_LEN];
351 int nextblock;
352 my_fread(&blksize1, sizeof(int), 1, fd);
353 my_fread(&label, sizeof(char), LABEL_LEN, fd);
354 my_fread(&nextblock, sizeof(int), 1, fd);
355 printf("%s Reading header => '%c%c%c%c' (%d byte)\n", info, label[0], label[1], label[2], label[3], nextblock);
356 my_fread(&blksize2, sizeof(int), 1, fd);
357 }
358
359 my_fread(&blksize1, sizeof(int), 1, fd);
360 my_fread(header_buf, header_size, 1, fd);
361 my_fread(&blksize2, sizeof(int), 1, fd);
362
363#ifdef GADGET2_HEADER
364 if(blksize1 != 256 || blksize2 != 256)
365 Terminate("incorrect GADGET2 header format, blksize1=%d blksize2=%d header_size=%d\n", blksize1, blksize2,
366 (int)header_size);
367#else
368 if(blksize1 != blksize2)
369 Terminate("incorrect header format, blksize1=%d blksize2=%d header_size=%d \n%s \n", blksize1, blksize2, (int)header_size,
370 blksize1 == 256 ? "You may need to set GADGET2_HEADER" : "");
371#endif
372 fclose(fd);
373 }
374 else
375 Terminate("illegal format");
376
377 for(int task = readTask + 1; task <= lastTask; task++)
378 MPI_Ssend(header_buf, header_size, MPI_BYTE, task, TAG_HEADER, Communicator);
379 }
380 else
381 MPI_Recv(header_buf, header_size, MPI_BYTE, readTask, TAG_HEADER, Communicator, MPI_STATUS_IGNORE);
382
383 read_file_header(fname, filenr, readTask, lastTask, n_type, npart, NULL);
384
385 if(ThisTask == readTask)
386 {
387 mpi_printf("READIC: Reading file `%s' on task=%d and distribute it to %d to %d.\n", fname, ThisTask, readTask, lastTask);
388 myflush(stdout);
389 }
390
391 for(int type = 0; type < N_DataGroups; type++)
392 {
393 ntype_in_files[filenr * N_DataGroups + type] = npart[type];
394
395 long long n_in_file = npart[type];
396 int ntask = lastTask - readTask + 1;
397 int n_for_this_task = n_in_file / ntask;
398 if((ThisTask - readTask) < (n_in_file % ntask))
399 n_for_this_task++;
400
401 read_increase_numbers(type, n_for_this_task);
402 }
403}
404
414void IO_Def::fill_write_buffer(int blocknr, int *startindex, int pc, int type, void *CommBuffer)
415{
416 if(blocknr < 0 || blocknr >= N_IO_Fields)
417 Terminate("something is wrong here: blocknr=%d N_IO_Fields=%d", blocknr, N_IO_Fields);
418
419 IO_Field *field = &IO_Fields[blocknr];
420
421 int *intp = (int *)CommBuffer;
422 long long *longp = (long long *)CommBuffer;
423 float *floatp = (float *)CommBuffer;
424 double *doublep = (double *)CommBuffer;
425
426 MyIDType *ip = (MyIDType *)CommBuffer;
427 MyFloat *fp = (MyFloat *)CommBuffer;
428 MyDouble *dp = (MyDouble *)CommBuffer;
429 MyIntPosType *intposp = (MyIntPosType *)CommBuffer;
430
431 int pindex = *startindex;
432
433 for(int n = 0; n < pc; pindex++)
434 {
435 if(type_of_file == FILE_IS_SNAPSHOT && type < NTYPES && get_type_of_element(pindex) != type)
436 continue;
437#ifdef LIGHTCONE
438 if(type_of_file == FILE_IS_LIGHTCONE && type < NTYPES && get_type_of_element(pindex) != type)
439 continue;
440#endif
441 if(field->io_func)
442 {
443 switch(field->type_in_memory)
444 {
445 case MEM_INT:
446 field->io_func(this, pindex, field->values_per_block, intp, 0);
447 intp += field->values_per_block;
448 break;
449 case MEM_INT64:
450 field->io_func(this, pindex, field->values_per_block, longp, 0);
451 longp += field->values_per_block;
452 break;
453 case MEM_MY_ID_TYPE:
454 field->io_func(this, pindex, field->values_per_block, ip, 0);
455 ip += field->values_per_block;
456 break;
457 case MEM_MY_INTPOS_TYPE:
458 field->io_func(this, pindex, field->values_per_block, intposp, 0);
459 intposp += field->values_per_block;
460 break;
461 case MEM_FLOAT:
462 field->io_func(this, pindex, field->values_per_block, floatp, 0);
463 floatp += field->values_per_block;
464 break;
465 case MEM_DOUBLE:
466 field->io_func(this, pindex, field->values_per_block, doublep, 0);
467 doublep += field->values_per_block;
468 break;
469 case MEM_MY_FLOAT:
470 field->io_func(this, pindex, field->values_per_block, fp, 0);
471 fp += field->values_per_block;
472 break;
473 case MEM_MY_DOUBLE:
474 field->io_func(this, pindex, field->values_per_block, dp, 0);
475 dp += field->values_per_block;
476 break;
477 default:
478 Terminate("ERROR in fill_write_buffer: Type not found!\n");
479
480 break;
481 }
482 }
483 else
484 {
485 void *array_pos;
486
487 switch(field->array)
488 {
489 case A_NONE:
490 array_pos = NULL;
491 break;
492
493 default:
494 array_pos = get_base_address_of_structure(field->array, pindex);
495 break;
496 }
497
498 for(int k = 0; k < field->values_per_block; k++)
499 {
500 switch(field->type_in_memory)
501 {
502 case MEM_INT:
503 *intp++ = *((int *)((size_t)array_pos + field->offset + k * sizeof(int)));
504 break;
505
506 case MEM_INT64:
507 *longp++ = *((long long *)((size_t)array_pos + field->offset + k * sizeof(long long)));
508 break;
509
510 case MEM_MY_ID_TYPE:
511 *ip++ = *((MyIDType *)((size_t)array_pos + field->offset + k * sizeof(MyIDType)));
512 break;
513
514 case MEM_MY_INTPOS_TYPE:
515 *intposp++ = *((MyIntPosType *)((size_t)array_pos + field->offset + k * sizeof(MyIntPosType)));
516 break;
517
518 case MEM_FLOAT:
519 *floatp++ = *((float *)((size_t)array_pos + field->offset + k * sizeof(float)));
520 break;
521
522 case MEM_DOUBLE:
523 *doublep++ = *((double *)((size_t)array_pos + field->offset + k * sizeof(double)));
524 break;
525
526 case MEM_MY_FLOAT:
527 *fp++ = *((MyFloat *)((size_t)array_pos + field->offset + k * sizeof(MyFloat)));
528 break;
529
530 case MEM_MY_DOUBLE:
531 *dp++ = *((MyDouble *)((size_t)array_pos + field->offset + k * sizeof(MyDouble)));
532 break;
533
534 default:
535 Terminate("ERROR in fill_write_buffer: Type not found!\n");
536 break;
537 }
538 }
539 }
540
541 n++;
542 }
543
544 *startindex = pindex;
545}
546
556void IO_Def::empty_read_buffer(int blocknr, int offset, int pc, int type, long long nprevious, void *CommBuffer)
557{
558 IO_Field *field = &IO_Fields[blocknr];
559
560 int *intp = (int *)CommBuffer;
561 long long *longp = (long long *)CommBuffer;
562 float *floatp = (float *)CommBuffer;
563 double *doublep = (double *)CommBuffer;
564
565 MyIDType *ip = (MyIDType *)CommBuffer;
566 MyFloat *fp = (MyFloat *)CommBuffer;
567 MyDouble *dp = (MyDouble *)CommBuffer;
568 MyIntPosType *intposp = (MyIntPosType *)CommBuffer;
569
570 if(field->read_flag != SKIP_ON_READ || field->type_in_memory == MEM_MY_FILEOFFSET)
571 {
572 for(int n = 0; n < pc; n++)
573 {
574 if(field->io_func)
575 {
576 switch(field->type_in_memory)
577 {
578 case MEM_INT:
579 field->io_func(this, offset + n, field->values_per_block, intp, 1);
580 intp += field->values_per_block;
581 break;
582 case MEM_INT64:
583 field->io_func(this, offset + n, field->values_per_block, longp, 1);
584 longp += field->values_per_block;
585 break;
586 case MEM_MY_ID_TYPE:
587 field->io_func(this, offset + n, field->values_per_block, ip, 1);
588 ip += field->values_per_block;
589 break;
590 case MEM_MY_INTPOS_TYPE:
591 field->io_func(this, offset + n, field->values_per_block, intposp, 1);
592 intposp += field->values_per_block;
593 break;
594 case MEM_FLOAT:
595 field->io_func(this, offset + n, field->values_per_block, floatp, 1);
596 floatp += field->values_per_block;
597 break;
598 case MEM_DOUBLE:
599 field->io_func(this, offset + n, field->values_per_block, doublep, 1);
600 doublep += field->values_per_block;
601 break;
602 case MEM_MY_FLOAT:
603 field->io_func(this, offset + n, field->values_per_block, fp, 1);
604 fp += field->values_per_block;
605 break;
606 case MEM_MY_DOUBLE:
607 field->io_func(this, offset + n, field->values_per_block, dp, 1);
608 dp += field->values_per_block;
609 break;
610 case MEM_MY_FILEOFFSET:
611 Terminate("undefined");
612 break;
613 }
614 }
615 else
616 {
617 void *array_pos;
618 switch(field->array)
619 {
620 case A_NONE:
621 array_pos = 0;
622 break;
623
624 default:
625 array_pos = get_base_address_of_structure(field->array, offset + n);
626 break;
627 }
628
629 for(int k = 0; k < field->values_per_block; k++)
630 {
631 switch(field->type_in_memory)
632 {
633 case MEM_INT:
634 *((int *)((size_t)array_pos + field->offset + k * sizeof(int))) = *intp++;
635 break;
636 case MEM_INT64:
637 *((long long *)((size_t)array_pos + field->offset + k * sizeof(long long))) = *longp++;
638 break;
639 case MEM_MY_ID_TYPE:
640 *((MyIDType *)((size_t)array_pos + field->offset + k * sizeof(MyIDType))) = *ip++;
641 break;
642 case MEM_MY_INTPOS_TYPE:
643 *((MyIntPosType *)((size_t)array_pos + field->offset + k * sizeof(MyIntPosType))) = *intposp++;
644 break;
645 case MEM_FLOAT:
646 *((float *)((size_t)array_pos + field->offset + k * sizeof(float))) = *floatp++;
647 break;
648 case MEM_DOUBLE:
649 *((double *)((size_t)array_pos + field->offset + k * sizeof(double))) = *doublep++;
650 break;
651 case MEM_MY_FLOAT:
652 *((MyFloat *)((size_t)array_pos + field->offset + k * sizeof(MyFloat))) = *fp++;
653 break;
654 case MEM_MY_DOUBLE:
655 *((MyDouble *)((size_t)array_pos + field->offset + k * sizeof(MyDouble))) = *dp++;
656 break;
657 case MEM_MY_FILEOFFSET:
658 *((long long *)((size_t)array_pos + field->offset + k * sizeof(long long))) = nprevious++;
659 break;
660 default:
661 Terminate("ERROR: Type not found!\n");
662 break;
663 }
664 }
665 }
666 }
667 }
668}
669
670void IO_Def::polling(int numfilesperdump)
671{
672 if(ThisTask == 0)
673 if(files_completed < numfilesperdump)
674 {
675 MPI_Status status;
676 int flag;
677
678 /* now check for a completion message */
679 MPI_Iprobe(MPI_ANY_SOURCE, TAG_KEY, Communicator, &flag, &status);
680
681 if(flag)
682 {
683 int source = status.MPI_SOURCE;
684
685 int dummy;
686 MPI_Recv(&dummy, 1, MPI_INT, source, TAG_KEY, Communicator, MPI_STATUS_IGNORE);
687 files_completed++;
688
689 if(files_started < numfilesperdump)
690 {
691 /* send start signal */
692 MPI_Ssend(&ThisTask, 1, MPI_INT, seq[files_started++].thistask, TAG_N, Communicator);
693 }
694 }
695 }
696}
697
698/* driver routine for outputting multiple files, scheduled in optimal order under the constraint not to write more than a certain
699 * number of files simultanuously */
700void IO_Def::write_multiple_files(char *fname, int numfilesperdump, int append_flag, int chunksize)
701{
702 if(ThisTask == 0)
703 if(!(seq = (seq_data *)Mem.mymalloc("seq", NTask * sizeof(seq_data))))
704 Terminate("can't allocate seq_data");
705
706 void *CommBuffer = Mem.mymalloc("CommBuffer", COMMBUFFERSIZE);
707
708 /* assign processors to output files */
709 int filenr, masterTask, lastTask;
710 distribute_file(numfilesperdump, &filenr, &masterTask, &lastTask);
711
712 char buf[MAXLEN_PATH];
713 if(numfilesperdump > 1)
714 sprintf(buf, "%s.%d", fname, filenr);
715 else
716 sprintf(buf, "%s", fname);
717
718 seq_data seq_loc;
719 seq_loc.thistask = ThisTask;
720 seq_loc.rankinnode = RankInThisNode;
721 seq_loc.thisnode = ThisNode;
722
723 if(masterTask != ThisTask)
724 seq_loc.thistask = -1;
725
726 MPI_Gather(&seq_loc, sizeof(seq_data), MPI_BYTE, seq, sizeof(seq_data), MPI_BYTE, 0, Communicator);
727
728 if(ThisTask == 0)
729 {
730 int count = NTask;
731 for(int i = 0; i < count; i++)
732 {
733 if(seq[i].thistask < 0)
734 {
735 count--;
736 seq[i] = seq[count];
737 i--;
738 }
739 }
740 if(count != numfilesperdump)
741 Terminate("count=%d != numfilesperdump=%d", count, numfilesperdump);
742
743 std::sort(seq, seq + numfilesperdump);
744
745 files_started = 0;
746 files_completed = 0;
747
748 for(int i = 1; i < std::min<int>(All.MaxFilesWithConcurrentIO, numfilesperdump); i++)
749 {
750 files_started++;
751 MPI_Ssend(&ThisTask, 1, MPI_INT, seq[i].thistask, TAG_N, Communicator);
752 }
753
754 files_started++;
755 if(append_flag)
756 append_file(buf, masterTask, lastTask, CommBuffer, numfilesperdump, chunksize);
757 else
758 write_file(buf, masterTask, lastTask, CommBuffer, numfilesperdump, chunksize);
759 files_completed++;
760
761 if(files_started < numfilesperdump)
762 {
763 /* send start signal */
764 MPI_Ssend(&ThisTask, 1, MPI_INT, seq[files_started++].thistask, TAG_N, Communicator);
765 }
766
767 while(files_completed < numfilesperdump)
768 polling(numfilesperdump);
769 }
770 else if(masterTask == ThisTask)
771 {
772 /* wait for start signal */
773 int dummy;
774 MPI_Recv(&dummy, 1, MPI_INT, 0, TAG_N, Communicator, MPI_STATUS_IGNORE); /* wait until we are told to start */
775
776 if(append_flag)
777 append_file(buf, masterTask, lastTask, CommBuffer, numfilesperdump, chunksize);
778 else
779 write_file(buf, masterTask, lastTask, CommBuffer, numfilesperdump, chunksize);
780
781 /* send back completion notice */
782 MPI_Ssend(&ThisTask, 1, MPI_INT, 0, TAG_KEY, Communicator);
783 }
784 else
785 {
786 if(append_flag)
787 append_file(buf, masterTask, lastTask, CommBuffer, numfilesperdump, chunksize);
788 else
789 write_file(buf, masterTask, lastTask, CommBuffer, numfilesperdump, chunksize);
790 }
791
792 Mem.myfree(CommBuffer);
793
794 if(ThisTask == 0)
795 Mem.myfree(seq);
796}
797
817void IO_Def::write_file(char *fname, int writeTask, int lastTask, void *CommBuffer, int numfilesperdump, int chunksize)
818{
819 int typelist[N_DataGroups];
820 long long n_type[N_DataGroups], npart[N_DataGroups], pcsum = 0;
821 char label[LABEL_LEN + 1];
822 unsigned int blksize, bytes_per_blockelement_in_file = 0;
823 FILE *fd = 0;
824 hid_t hdf5_file = 0, hdf5_grp[N_DataGroups], hdf5_headergrp = 0, hdf5_dataspace_memory;
825 hid_t hdf5_dataspace_in_file = 0, hdf5_dataset = 0, hdf5_prop = 0;
826 hsize_t dims[2], count[2], start[2];
827 int rank = 0;
828 hid_t hdf5_paramsgrp = 0;
829 hid_t hdf5_configgrp = 0;
830
831#define SKIP \
832 { \
833 my_fwrite(&blksize, sizeof(int), 1, fd); \
834 }
835
836 fill_file_header(writeTask, lastTask, n_type, npart);
837
838 /* open file and write header */
839 if(ThisTask == writeTask)
840 {
841 if(file_format == FILEFORMAT_HDF5)
842 {
843 char buf[MAXLEN_PATH];
844 sprintf(buf, "%s.hdf5", fname);
845 mpi_printf("%s file: '%s' (file 1 of %d)\n", info, fname, numfilesperdump);
846
847 rename_file_to_bak_if_it_exists(buf);
848
849 hdf5_file = my_H5Fcreate(buf, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
850
851 hdf5_headergrp = my_H5Gcreate(hdf5_file, "/Header", 0);
852 write_header_fields(hdf5_headergrp);
853
854 hdf5_paramsgrp = my_H5Gcreate(hdf5_file, "/Parameters", 0);
855 write_parameters_attributes_in_hdf5(hdf5_paramsgrp);
856
857 hdf5_configgrp = my_H5Gcreate(hdf5_file, "/Config", 0);
858 write_compile_time_options_in_hdf5(hdf5_configgrp);
859
860 for(int type = 0; type < N_DataGroups; type++)
861 {
862 if(npart[type] > 0)
863 {
864 get_datagroup_name(type, buf);
865 hdf5_grp[type] = my_H5Gcreate(hdf5_file, buf, 0);
866 }
867 }
868 }
869 else
870 {
871 rename_file_to_bak_if_it_exists(fname);
872
873 if(!(fd = fopen(fname, "w")))
874 Terminate("can't open file `%s' for writing.\n", fname);
875
876 mpi_printf("%s file: '%s' (file 1 of %d)\n", info, fname, numfilesperdump);
877
878 if(file_format == FILEFORMAT_LEGACY2)
879 {
880 blksize = sizeof(int) + 4 * sizeof(char);
881 SKIP;
882 my_fwrite((const void *)"HEAD", sizeof(char), 4, fd);
883 int nextblock = header_size + 2 * sizeof(int);
884 my_fwrite(&nextblock, sizeof(int), 1, fd);
885 SKIP;
886 }
887
888 blksize = header_size;
889 SKIP;
890 my_fwrite(header_buf, header_size, 1, fd);
891 SKIP;
892 }
893 }
894
895 for(int blocknr = 0; blocknr < N_IO_Fields; blocknr++)
896 {
897 polling(numfilesperdump);
898
899 if(IO_Fields[blocknr].type_in_file_output != FILE_NONE)
900 {
901 unsigned int bytes_per_blockelement = get_bytes_per_memory_blockelement(blocknr, 0);
902 int blockmaxlen = (int)(COMMBUFFERSIZE / bytes_per_blockelement);
903 long long npart_in_block = get_particles_in_block(blocknr, npart, typelist);
904 hid_t hdf5_memory_datatype = get_hdf5_memorytype_of_block(blocknr);
905 char dname[MAXLEN_PATH];
906 get_dataset_name(blocknr, dname);
907
908 if(npart_in_block > 0)
909 {
910 mpi_printf("%s block %d (%s)...\n", info, blocknr, dname);
911
912 if(ThisTask == writeTask)
913 {
914 if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
915 {
916 bytes_per_blockelement_in_file =
917 IO_Fields[blocknr].values_per_block * H5Tget_size(get_hdf5_outputtype_of_block(blocknr));
918
919 if(file_format == FILEFORMAT_LEGACY2)
920 {
921 blksize = sizeof(int) + LABEL_LEN * sizeof(char);
922 SKIP;
923 get_Tab_IO_Label(blocknr, label);
924 my_fwrite(label, sizeof(char), LABEL_LEN, fd);
925 int nextblock = npart_in_block * bytes_per_blockelement_in_file + 2 * sizeof(int);
926 my_fwrite(&nextblock, sizeof(int), 1, fd);
927 SKIP;
928 }
929
930 blksize = npart_in_block * bytes_per_blockelement_in_file;
931 SKIP;
932 }
933 }
934
935 for(int type = 0; type < N_DataGroups; type++)
936 {
937 if(typelist[type])
938 {
939 if(ThisTask == writeTask && file_format == FILEFORMAT_HDF5 && npart[type] > 0)
940 {
941 hid_t hdf5_file_datatype = get_hdf5_outputtype_of_block(blocknr);
942
943 dims[0] = npart[type];
944 dims[1] = get_values_per_blockelement(blocknr);
945 if(dims[1] == 1)
946 rank = 1;
947 else
948 rank = 2;
949
950 if(chunksize > 0)
951 {
952 hsize_t maxdims[2] = {H5S_UNLIMITED, dims[1]};
953 hdf5_dataspace_in_file = my_H5Screate_simple(rank, dims, maxdims);
954
955 /* Modify dataset creation properties, i.e. enable chunking */
956 hdf5_prop = H5Pcreate(H5P_DATASET_CREATE);
957 hsize_t chunk_dims[2] = {0, dims[1]};
958 chunk_dims[0] = chunksize;
959 H5Pset_chunk(hdf5_prop, rank, chunk_dims);
960
961 hdf5_dataset =
962 my_H5Dcreate(hdf5_grp[type], dname, hdf5_file_datatype, hdf5_dataspace_in_file, hdf5_prop);
963 }
964 else if(IO_Fields[blocknr].compression_on)
965 {
966 hdf5_dataspace_in_file = my_H5Screate_simple(rank, dims, NULL);
967
968 /* Modify dataset creation properties, i.e. enable compression */
969 hdf5_prop = H5Pcreate(H5P_DATASET_CREATE);
970 hsize_t chunk_dims[2] = {COMPRESSION_CHUNKSIZE, dims[1]};
971 if(chunk_dims[0] > dims[0])
972 chunk_dims[0] = dims[0];
973 H5Pset_chunk(hdf5_prop, rank, chunk_dims); /* set chunk size */
974 H5Pset_shuffle(hdf5_prop); /* reshuffle bytes to get better compression ratio */
975 H5Pset_deflate(hdf5_prop, 9); /* gzip compression level 9 */
976 if(H5Pall_filters_avail(hdf5_prop))
977 hdf5_dataset =
978 my_H5Dcreate(hdf5_grp[type], dname, hdf5_file_datatype, hdf5_dataspace_in_file, hdf5_prop);
979 else
980 Terminate("HDF5: Compression not available!\n");
981 }
982 else
983 {
984 hdf5_dataspace_in_file = my_H5Screate_simple(rank, dims, NULL);
985 hdf5_dataset =
986 my_H5Dcreate(hdf5_grp[type], dname, hdf5_file_datatype, hdf5_dataspace_in_file, H5P_DEFAULT);
987 }
988
989 write_dataset_attributes(hdf5_dataset, blocknr);
990
991 byte_count += dims[0] * dims[1] * my_H5Tget_size(hdf5_file_datatype); /* for I/O performance measurement */
992
993 pcsum = 0;
994 }
995
996 for(int task = writeTask, offset = 0; task <= lastTask; task++)
997 {
998 long long n_for_this_task;
999
1000 if(task == ThisTask)
1001 {
1002 n_for_this_task = n_type[type];
1003
1004 for(int p = writeTask; p <= lastTask; p++)
1005 if(p != ThisTask)
1006 MPI_Send(&n_for_this_task, sizeof(n_for_this_task), MPI_BYTE, p, TAG_NFORTHISTASK, Communicator);
1007 }
1008 else
1009 MPI_Recv(&n_for_this_task, sizeof(n_for_this_task), MPI_BYTE, task, TAG_NFORTHISTASK, Communicator,
1010 MPI_STATUS_IGNORE);
1011
1012 while(n_for_this_task > 0)
1013 {
1014 long long pc = n_for_this_task;
1015
1016 if(pc > blockmaxlen)
1017 pc = blockmaxlen;
1018
1019 if(ThisTask == task)
1020 fill_write_buffer(blocknr, &offset, pc, type, CommBuffer);
1021
1022 if(ThisTask == writeTask && task != writeTask)
1023 MPI_Recv(CommBuffer, bytes_per_blockelement * pc, MPI_BYTE, task, TAG_PDATA, Communicator,
1024 MPI_STATUS_IGNORE);
1025
1026 if(ThisTask != writeTask && task == ThisTask)
1027 MPI_Ssend(CommBuffer, bytes_per_blockelement * pc, MPI_BYTE, writeTask, TAG_PDATA, Communicator);
1028
1029 if(ThisTask == writeTask)
1030 {
1031 if(file_format == FILEFORMAT_HDF5)
1032 {
1033 start[0] = pcsum;
1034 start[1] = 0;
1035
1036 count[0] = pc;
1037 count[1] = get_values_per_blockelement(blocknr);
1038 pcsum += pc;
1039
1040 my_H5Sselect_hyperslab(hdf5_dataspace_in_file, H5S_SELECT_SET, start, NULL, count, NULL);
1041
1042 dims[0] = pc;
1043 dims[1] = get_values_per_blockelement(blocknr);
1044 hdf5_dataspace_memory = my_H5Screate_simple(rank, dims, NULL);
1045
1046 my_H5Dwrite(hdf5_dataset, hdf5_memory_datatype, hdf5_dataspace_memory, hdf5_dataspace_in_file,
1047 H5P_DEFAULT, CommBuffer, dname);
1048
1049 my_H5Sclose(hdf5_dataspace_memory, H5S_SIMPLE);
1050 }
1051 else
1052 {
1053 if(bytes_per_blockelement_in_file != bytes_per_blockelement)
1054 {
1055 char *CommAuxBuffer =
1056 (char *)Mem.mymalloc("CommAuxBuffer", bytes_per_blockelement_in_file * pc);
1057 type_cast_data((char *)CommBuffer, bytes_per_blockelement, (char *)CommAuxBuffer,
1058 bytes_per_blockelement_in_file, pc, blocknr);
1059 my_fwrite(CommAuxBuffer, bytes_per_blockelement_in_file, pc, fd);
1060 Mem.myfree(CommAuxBuffer);
1061 }
1062 else
1063 my_fwrite(CommBuffer, bytes_per_blockelement, pc, fd);
1064 }
1065 }
1066
1067 n_for_this_task -= pc;
1068 }
1069 }
1070
1071 if(ThisTask == writeTask && file_format == FILEFORMAT_HDF5 && npart[type] > 0)
1072 {
1073 if(file_format == FILEFORMAT_HDF5)
1074 {
1075 my_H5Dclose(hdf5_dataset, dname);
1076 if(chunksize > 0 || IO_Fields[blocknr].compression_on)
1077 my_H5Pclose(hdf5_prop);
1078 my_H5Sclose(hdf5_dataspace_in_file, H5S_SIMPLE);
1079 }
1080 }
1081 }
1082 }
1083
1084 if(ThisTask == writeTask)
1085 {
1086 if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
1087 SKIP;
1088 }
1089 }
1090 }
1091 }
1092
1093 if(ThisTask == writeTask)
1094 {
1095 if(file_format == FILEFORMAT_HDF5)
1096 {
1097 char buf[MAXLEN_PATH];
1098
1099 for(int type = N_DataGroups - 1; type >= 0; type--)
1100 if(npart[type] > 0)
1101 {
1102 get_datagroup_name(type, buf);
1103 my_H5Gclose(hdf5_grp[type], buf);
1104 }
1105
1106 my_H5Gclose(hdf5_headergrp, "/Header");
1107 my_H5Gclose(hdf5_paramsgrp, "/Parameters");
1108 my_H5Gclose(hdf5_configgrp, "/Config");
1109
1110 sprintf(buf, "%s.hdf5", fname);
1111 my_H5Fclose(hdf5_file, buf);
1112 }
1113 else
1114 fclose(fd);
1115 }
1116}
1117
1118void IO_Def::append_file(char *fname, int writeTask, int lastTask, void *CommBuffer, int numfilesperdump, int chunksize)
1119{
1120 int typelist[N_DataGroups];
1121 long long n_type[N_DataGroups], npart[N_DataGroups], n_previous[N_DataGroups], pcsum = 0;
1122 hid_t hdf5_file = 0, hdf5_grp[N_DataGroups], hdf5_headergrp = 0, hdf5_dataspace_memory;
1123 hid_t hdf5_dataspace_in_file = 0, hdf5_dataset = 0;
1124 hsize_t dims[2], count[2], start[2];
1125 int rank = 0;
1126
1127 if(file_format != FILEFORMAT_HDF5)
1128 Terminate("appending to files only works with HDF5 format\n");
1129
1130 read_file_header(NULL, 0, writeTask, lastTask, n_type, npart, NULL);
1131
1132 for(int n = 0; n < N_DataGroups; n++)
1133 n_previous[n] = n_type[n];
1134
1135 fill_file_header(writeTask, lastTask, n_type, npart);
1136
1137 /* open file and write header */
1138 if(ThisTask == writeTask)
1139 {
1140 char buf[MAXLEN_PATH];
1141 sprintf(buf, "%s.hdf5", fname);
1142
1143 hdf5_file = my_H5Fopen(buf, H5F_ACC_RDWR, H5P_DEFAULT);
1144
1145 hdf5_headergrp = my_H5Gopen(hdf5_file, "/Header");
1146 write_header_fields(hdf5_headergrp);
1147
1148 for(int type = 0; type < N_DataGroups; type++)
1149 {
1150 if(npart[type] > 0)
1151 {
1152 get_datagroup_name(type, buf);
1153 if(n_previous[type] == 0)
1154 hdf5_grp[type] = my_H5Gcreate(hdf5_file, buf, 0);
1155 else
1156 hdf5_grp[type] = my_H5Gopen(hdf5_file, buf);
1157 }
1158 }
1159 }
1160
1161 for(int blocknr = 0; blocknr < N_IO_Fields; blocknr++)
1162 {
1163 polling(numfilesperdump);
1164
1165 if(IO_Fields[blocknr].type_in_file_output != FILE_NONE)
1166 {
1167 unsigned int bytes_per_blockelement = get_bytes_per_memory_blockelement(blocknr, 0);
1168 int blockmaxlen = (int)(COMMBUFFERSIZE / bytes_per_blockelement);
1169 long long npart_in_block = get_particles_in_block(blocknr, npart, typelist);
1170 hid_t hdf5_memory_datatype = get_hdf5_memorytype_of_block(blocknr);
1171 char dname[MAXLEN_PATH];
1172 get_dataset_name(blocknr, dname);
1173
1174 if(npart_in_block > 0)
1175 {
1176 mpi_printf("%s block %d (%s)...\n", info, blocknr, dname);
1177
1178 for(int type = 0; type < N_DataGroups; type++)
1179 {
1180 if(typelist[type])
1181 {
1182 if(ThisTask == writeTask && file_format == FILEFORMAT_HDF5 && npart[type] > 0)
1183 {
1184 hid_t hdf5_file_datatype = get_hdf5_outputtype_of_block(blocknr);
1185
1186 dims[0] = npart[type] + n_previous[type];
1187 dims[1] = get_values_per_blockelement(blocknr);
1188 if(dims[1] == 1)
1189 rank = 1;
1190 else
1191 rank = 2;
1192
1193 hdf5_dataspace_in_file = my_H5Screate_simple(rank, dims, NULL);
1194
1195 if(n_previous[type] == 0)
1196 {
1197 if(chunksize > 0)
1198 {
1199 hsize_t maxdims[2] = {H5S_UNLIMITED, dims[1]};
1200 hdf5_dataspace_in_file = my_H5Screate_simple(rank, dims, maxdims);
1201
1202 /* Modify dataset creation properties, i.e. enable chunking */
1203 hid_t prop = H5Pcreate(H5P_DATASET_CREATE);
1204 hsize_t chunk_dims[2] = {0, dims[1]};
1205 chunk_dims[0] = chunksize;
1206 H5Pset_chunk(prop, rank, chunk_dims);
1207
1208 hdf5_dataset = my_H5Dcreate(hdf5_grp[type], dname, hdf5_file_datatype, hdf5_dataspace_in_file, prop);
1209 }
1210 else
1211 {
1212 hdf5_dataset =
1213 my_H5Dcreate(hdf5_grp[type], dname, hdf5_file_datatype, hdf5_dataspace_in_file, H5P_DEFAULT);
1214 write_dataset_attributes(hdf5_dataset, blocknr);
1215 }
1216 }
1217 else
1218 {
1219 hdf5_dataset = my_H5Dopen_if_existing(hdf5_grp[type], dname);
1220 my_H5Dset_extent(hdf5_dataset, dims);
1221 }
1222
1223 byte_count += dims[0] * dims[1] * my_H5Tget_size(hdf5_file_datatype); /* for I/O performance measurement */
1224
1225 pcsum = 0;
1226 }
1227
1228 for(int task = writeTask, offset = 0; task <= lastTask; task++)
1229 {
1230 long long n_for_this_task;
1231
1232 if(task == ThisTask)
1233 {
1234 n_for_this_task = n_type[type];
1235
1236 for(int p = writeTask; p <= lastTask; p++)
1237 if(p != ThisTask)
1238 MPI_Send(&n_for_this_task, sizeof(n_for_this_task), MPI_BYTE, p, TAG_NFORTHISTASK, Communicator);
1239 }
1240 else
1241 MPI_Recv(&n_for_this_task, sizeof(n_for_this_task), MPI_BYTE, task, TAG_NFORTHISTASK, Communicator,
1242 MPI_STATUS_IGNORE);
1243
1244 while(n_for_this_task > 0)
1245 {
1246 long long pc = n_for_this_task;
1247
1248 if(pc > blockmaxlen)
1249 pc = blockmaxlen;
1250
1251 if(ThisTask == task)
1252 fill_write_buffer(blocknr, &offset, pc, type, CommBuffer);
1253
1254 if(ThisTask == writeTask && task != writeTask)
1255 MPI_Recv(CommBuffer, bytes_per_blockelement * pc, MPI_BYTE, task, TAG_PDATA, Communicator,
1256 MPI_STATUS_IGNORE);
1257
1258 if(ThisTask != writeTask && task == ThisTask)
1259 MPI_Ssend(CommBuffer, bytes_per_blockelement * pc, MPI_BYTE, writeTask, TAG_PDATA, Communicator);
1260
1261 if(ThisTask == writeTask)
1262 {
1263 start[0] = pcsum + n_previous[type];
1264 start[1] = 0;
1265
1266 count[0] = pc;
1267 count[1] = get_values_per_blockelement(blocknr);
1268 pcsum += pc;
1269
1270 my_H5Sselect_hyperslab(hdf5_dataspace_in_file, H5S_SELECT_SET, start, NULL, count, NULL);
1271
1272 dims[0] = pc;
1273 dims[1] = get_values_per_blockelement(blocknr);
1274 hdf5_dataspace_memory = my_H5Screate_simple(rank, dims, NULL);
1275
1276 my_H5Dwrite(hdf5_dataset, hdf5_memory_datatype, hdf5_dataspace_memory, hdf5_dataspace_in_file,
1277 H5P_DEFAULT, CommBuffer, dname);
1278
1279 my_H5Sclose(hdf5_dataspace_memory, H5S_SIMPLE);
1280 }
1281
1282 n_for_this_task -= pc;
1283 }
1284 }
1285
1286 if(ThisTask == writeTask && npart[type] > 0)
1287 {
1288 my_H5Dclose(hdf5_dataset, dname);
1289 my_H5Sclose(hdf5_dataspace_in_file, H5S_SIMPLE);
1290 }
1291 }
1292 }
1293 }
1294 }
1295 }
1296
1297 if(ThisTask == writeTask)
1298 {
1299 char buf[MAXLEN_PATH];
1300
1301 for(int type = N_DataGroups - 1; type >= 0; type--)
1302 if(npart[type] > 0)
1303 {
1304 get_datagroup_name(type, buf);
1305 my_H5Gclose(hdf5_grp[type], buf);
1306 }
1307
1308 my_H5Gclose(hdf5_headergrp, "/Header");
1309
1310 sprintf(buf, "%s.hdf5", fname);
1311 my_H5Fclose(hdf5_file, buf);
1312 }
1313}
1314
1315void IO_Def::type_cast_data(char *src, int src_bytes_per_element, char *target, int target_bytes_per_element, int len, int blocknr)
1316{
1317 switch(IO_Fields[blocknr].type_in_memory)
1318 {
1319 case MEM_INT:
1320 case MEM_INT64:
1321 case MEM_MY_ID_TYPE:
1322 case MEM_MY_INTPOS_TYPE:
1323 case MEM_MY_FILEOFFSET:
1324 if(target_bytes_per_element > src_bytes_per_element)
1325 {
1326 if(target_bytes_per_element != 2 * src_bytes_per_element)
1327 Terminate("something is odd here: target_bytes_per_element=%d != 2 * src_bytes_per_element=%d", target_bytes_per_element,
1328 src_bytes_per_element);
1329
1330 int fac = src_bytes_per_element / sizeof(int);
1331 int *sp = (int *)src;
1332 long long *tp = (long long *)target;
1333 for(int i = 0; i < len * fac; i++)
1334 *tp++ = *sp++;
1335 }
1336 else
1337 {
1338 if(src_bytes_per_element != 2 * target_bytes_per_element)
1339 Terminate("something is odd here: src_bytes_per_element=%d != 2 * target_bytes_per_element=%d", src_bytes_per_element,
1340 target_bytes_per_element);
1341 int fac = src_bytes_per_element / sizeof(long long);
1342 long long *sp = (long long *)src;
1343 int *tp = (int *)target;
1344 for(int i = 0; i < len * fac; i++)
1345 *tp++ = *sp++;
1346 }
1347 break;
1348
1349 case MEM_FLOAT:
1350 case MEM_DOUBLE:
1351 case MEM_MY_FLOAT:
1352 case MEM_MY_DOUBLE:
1353 if((target_bytes_per_element % 8) == 0 &&
1354 (src_bytes_per_element % 4) == 0) // target_bytes_per_element multiply of 8, src_bytes_per_element mulitple of 4
1355 {
1356 int fac = src_bytes_per_element / sizeof(float);
1357 float *sp = (float *)src;
1358 double *tp = (double *)target;
1359 for(int i = 0; i < len * fac; i++)
1360 *tp++ = *sp++;
1361 }
1362 else if((target_bytes_per_element % 4) == 0 && (src_bytes_per_element % 8) == 0)
1363 {
1364 int fac = src_bytes_per_element / sizeof(double);
1365 double *sp = (double *)src;
1366 float *tp = (float *)target;
1367 for(int i = 0; i < len * fac; i++)
1368 *tp++ = *sp++;
1369 }
1370 else if((target_bytes_per_element & 8) == 0 && (src_bytes_per_element % 2) == 0)
1371 {
1372 int fac = src_bytes_per_element / sizeof(half);
1373 half *sp = (half *)src;
1374 double *tp = (double *)target;
1375 for(int i = 0; i < len * fac; i++)
1376 *tp++ = *sp++;
1377 }
1378 else if((target_bytes_per_element % 2) == 0 && (src_bytes_per_element % 8) == 0)
1379 {
1380 int fac = src_bytes_per_element / sizeof(double);
1381 double *sp = (double *)src;
1382 half *tp = (half *)target;
1383 for(int i = 0; i < len * fac; i++)
1384 *tp++ = *sp++;
1385 }
1386 else
1387 {
1388 Terminate("Strange conversion requested: target_bytes_per_element=%d src_bytes_per_element=%d\n",
1389 target_bytes_per_element, src_bytes_per_element);
1390 }
1391 break;
1392 }
1393}
1394
1404void IO_Def::read_file(const char *fname, int filenr, int readTask, int lastTask, void *CommBuffer)
1405{
1406 long long n_type[N_DataGroups], npart[N_DataGroups];
1407 int typelist[N_DataGroups];
1408 unsigned int blksize1, blksize2, bytes_per_blockelement_in_file = 0;
1409 hid_t hdf5_file = 0, hdf5_grp[N_DataGroups], hdf5_dataspace_in_file;
1410 hid_t hdf5_dataspace_in_memory, hdf5_dataset;
1411 FILE *fd = 0;
1412
1413 /* open file and read header */
1414
1415 if(ThisTask == readTask)
1416 {
1417 if(file_format == FILEFORMAT_HDF5)
1418 {
1419 read_header_fields(fname);
1420 }
1421 else if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
1422 {
1423 if(!(fd = fopen(fname, "r")))
1424 Terminate("can't open file `%s' for reading initial conditions.\n", fname);
1425
1426 if(file_format == FILEFORMAT_LEGACY2)
1427 {
1428 int nextblock;
1429 char label[LABEL_LEN];
1430 my_fread(&blksize1, sizeof(int), 1, fd);
1431 my_fread(&label, sizeof(char), LABEL_LEN, fd);
1432 my_fread(&nextblock, sizeof(int), 1, fd);
1433 my_fread(&blksize2, sizeof(int), 1, fd);
1434 }
1435
1436 my_fread(&blksize1, sizeof(int), 1, fd);
1437 my_fread(header_buf, header_size, 1, fd);
1438 my_fread(&blksize2, sizeof(int), 1, fd);
1439 }
1440 else
1441 Terminate("unknown ICFormat");
1442
1443 for(int task = readTask + 1; task <= lastTask; task++)
1444 MPI_Ssend(header_buf, header_size, MPI_BYTE, task, TAG_HEADER, Communicator);
1445 }
1446 else
1447 MPI_Recv(header_buf, header_size, MPI_BYTE, readTask, TAG_HEADER, Communicator, MPI_STATUS_IGNORE);
1448
1449 int nstart;
1450 read_file_header(fname, filenr, readTask, lastTask, n_type, npart, &nstart);
1451
1452 if(ThisTask == readTask)
1453 {
1454 if(file_format == FILEFORMAT_HDF5)
1455 {
1456 hdf5_file = my_H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
1457
1458 for(int type = 0; type < N_DataGroups; type++)
1459 {
1460 if(npart[type] > 0)
1461 {
1462 char buf[MAXLEN_PATH];
1463 get_datagroup_name(type, buf);
1464 hdf5_grp[type] = my_H5Gopen(hdf5_file, buf);
1465 }
1466 }
1467 }
1468 }
1469
1470 for(int blocknr = 0; blocknr < N_IO_Fields; blocknr++)
1471 {
1472 if((IO_Fields[blocknr].read_flag != SKIP_ON_READ &&
1473 !(file_format == FILEFORMAT_LEGACY1 && All.RestartFlag == RST_BEGIN && type_of_file == FILE_IS_SNAPSHOT &&
1474 blocknr > 4) /* this second conditions allows short legacy ICs to be read in */
1475 ) ||
1476 IO_Fields[blocknr].type_in_memory == MEM_MY_FILEOFFSET)
1477 {
1478 unsigned int bytes_per_blockelement = get_bytes_per_memory_blockelement(blocknr, 1);
1479 int blockmaxlen = (int)(COMMBUFFERSIZE / bytes_per_blockelement);
1480 long long npart_in_block = get_particles_in_block(blocknr, npart, &typelist[0]);
1481 hid_t hdf5_memory_datatype = get_hdf5_memorytype_of_block(blocknr);
1482 char dname[MAXLEN_PATH];
1483 get_dataset_name(blocknr, dname);
1484
1485 if(npart_in_block > 0)
1486 {
1487 if(filenr == 0)
1488 mpi_printf("READIC: reading block %d (%s)...\n", blocknr, dname);
1489
1490 if(ThisTask == readTask && IO_Fields[blocknr].type_in_memory != MEM_MY_FILEOFFSET)
1491 {
1492 if(file_format == FILEFORMAT_LEGACY2)
1493 {
1494 char expected_label[LABEL_LEN + 1];
1495 get_Tab_IO_Label(blocknr, expected_label);
1496 find_block(expected_label, fd);
1497 }
1498
1499 if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
1500 {
1501 my_fread(&blksize1, sizeof(int), 1, fd);
1502 bytes_per_blockelement_in_file = blksize1 / npart_in_block;
1503 }
1504 }
1505
1506 int offset = 0;
1507
1508 for(int type = 0; type < N_DataGroups; type++)
1509 {
1510 if(type_of_file != FILE_IS_SNAPSHOT)
1511 offset = 0;
1512
1513 int n_in_file = npart[type];
1514 int pcsum = 0;
1515
1516 long long nprevious = 0;
1517 for(int t = 0; t < type; t++)
1518 for(int f = 0; f < get_filenr_from_header(); f++)
1519 nprevious += ntype_in_files[f * N_DataGroups + t];
1520
1521 for(int nr = 0; nr < filenr; nr++)
1522 nprevious += ntype_in_files[nr * N_DataGroups + type];
1523
1524 if(typelist[type] == 0)
1525 {
1526 int ntask = lastTask - readTask + 1;
1527 int n_for_this_task = n_in_file / ntask;
1528 if((ThisTask - readTask) < (n_in_file % ntask))
1529 n_for_this_task++;
1530
1531 offset += n_for_this_task;
1532 }
1533 else
1534 {
1535 for(int task = readTask; task <= lastTask; task++)
1536 {
1537 int ntask = lastTask - readTask + 1;
1538 int n_for_this_task = n_in_file / ntask;
1539 if((task - readTask) < (n_in_file % ntask))
1540 n_for_this_task++;
1541
1542 do
1543 {
1544 int pc = n_for_this_task;
1545
1546 if(pc > blockmaxlen)
1547 pc = blockmaxlen;
1548
1549 if(ThisTask == readTask && IO_Fields[blocknr].type_in_memory != MEM_MY_FILEOFFSET)
1550 {
1551 if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
1552 {
1553 if(bytes_per_blockelement_in_file != bytes_per_blockelement)
1554 {
1555 char *CommAuxBuffer =
1556 (char *)Mem.mymalloc("CommAuxBuffer", bytes_per_blockelement_in_file * pc);
1557 my_fread(CommAuxBuffer, bytes_per_blockelement_in_file, pc, fd);
1558 type_cast_data((char *)CommAuxBuffer, bytes_per_blockelement_in_file, (char *)CommBuffer,
1559 bytes_per_blockelement, pc, blocknr);
1560 Mem.myfree(CommAuxBuffer);
1561 }
1562 else
1563 my_fread(CommBuffer, bytes_per_blockelement, pc, fd);
1564 }
1565
1566 if(file_format == FILEFORMAT_HDF5 && pc > 0)
1567 {
1568 hdf5_dataset = my_H5Dopen_if_existing(hdf5_grp[type], dname);
1569 int rank;
1570 hsize_t dims[2], count[2], start[2];
1571
1572 dims[0] = npart[type];
1573 dims[1] = get_values_per_blockelement(blocknr);
1574 if(dims[1] == 1)
1575 rank = 1;
1576 else
1577 rank = 2;
1578
1579 hdf5_dataspace_in_file = my_H5Screate_simple(rank, dims, NULL);
1580
1581 dims[0] = pc;
1582 hdf5_dataspace_in_memory = my_H5Screate_simple(rank, dims, NULL);
1583
1584 start[0] = pcsum;
1585 start[1] = 0;
1586
1587 count[0] = pc;
1588 count[1] = get_values_per_blockelement(blocknr);
1589 pcsum += pc;
1590
1591 my_H5Sselect_hyperslab(hdf5_dataspace_in_file, H5S_SELECT_SET, start, NULL, count, NULL);
1592
1593 // Test if dataset was present
1594 if(hdf5_dataset < 0)
1595 {
1596 // no, pad with zeros
1597 if((ThisTask == readTask) && (task == ThisTask))
1598 printf("\tDataset %s not present for particle type %d, using zero.\n", dname, type);
1599 memset(CommBuffer, 0, dims[0] * dims[1] * my_H5Tget_size(hdf5_memory_datatype));
1600 }
1601 else
1602 {
1603 hid_t hdf5_file_datatype = H5Dget_type(hdf5_dataset);
1604 byte_count += dims[0] * dims[1] *
1605 my_H5Tget_size(hdf5_file_datatype); /* for I/O performance measurement */
1606 H5Tclose(hdf5_file_datatype);
1607
1608 my_H5Dread(hdf5_dataset, hdf5_memory_datatype, hdf5_dataspace_in_memory,
1609 hdf5_dataspace_in_file, H5P_DEFAULT, CommBuffer, dname);
1610 my_H5Dclose(hdf5_dataset, dname);
1611 }
1612 my_H5Sclose(hdf5_dataspace_in_memory, H5S_SIMPLE);
1613 my_H5Sclose(hdf5_dataspace_in_file, H5S_SIMPLE);
1614 }
1615 }
1616
1617 if(ThisTask == readTask && task != readTask && pc > 0)
1618 MPI_Ssend(CommBuffer, bytes_per_blockelement * pc, MPI_BYTE, task, TAG_PDATA, Communicator);
1619
1620 if(ThisTask != readTask && task == ThisTask && pc > 0)
1621 MPI_Recv(CommBuffer, bytes_per_blockelement * pc, MPI_BYTE, readTask, TAG_PDATA, Communicator,
1622 MPI_STATUS_IGNORE);
1623
1624 if(ThisTask == task)
1625 {
1626 if(blocknr == 0 && IO_Fields[blocknr].array == A_P)
1627 {
1628 for(int n = 0; n < pc; n++)
1629 set_type_of_element(nstart + offset + n, type); /* initialize type */
1630 }
1631#ifdef MERGERTREE
1632 if(blocknr == 0 && IO_Fields[blocknr].array == A_MTRP)
1633 {
1634 for(int n = 0; n < pc; n++)
1635 {
1636 set_type_of_element(nstart + offset + n, type); /* initialize type */
1637 // MtrP[nstart + offset + n].Type = type; /* initialize type */
1638 }
1639 }
1640#endif
1641 empty_read_buffer(blocknr, nstart + offset, pc, type, nprevious, CommBuffer);
1642
1643 offset += pc;
1644 }
1645
1646 n_for_this_task -= pc;
1647 nprevious += pc;
1648 }
1649 while(n_for_this_task > 0);
1650 }
1651 }
1652 }
1653 if(ThisTask == readTask && IO_Fields[blocknr].type_in_memory != MEM_MY_FILEOFFSET)
1654 {
1655 if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
1656 {
1657 my_fread(&blksize2, sizeof(int), 1, fd);
1658 ;
1659 if(blksize1 != blksize2)
1660 {
1661 char buf[MAXLEN_PATH];
1662 sprintf(buf, "incorrect block-sizes detected!\n Task=%d blocknr=%d blksize1=%d blksize2=%d\n", ThisTask,
1663 blocknr, blksize1, blksize2);
1664 if(blocknr == 2) /* block number 2 is always IDs */
1665 strcat(buf, "Possible mismatch of 32bit and 64bit ID's in IC file and GADGET compilation !\n");
1666 Terminate(buf);
1667 }
1668 }
1669 }
1670 }
1671 }
1672 }
1673
1674 for(int type = 0; type < N_DataGroups; type++)
1675 {
1676 long long n_in_file = npart[type];
1677 int ntask = lastTask - readTask + 1;
1678 int n_for_this_task = n_in_file / ntask;
1679 if((ThisTask - readTask) < (n_in_file % ntask))
1680 n_for_this_task++;
1681
1682 read_increase_numbers(type, n_for_this_task);
1683 }
1684
1685 if(ThisTask == readTask)
1686 {
1687 if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
1688 fclose(fd);
1689
1690 if(file_format == FILEFORMAT_HDF5)
1691 {
1692 for(int type = N_DataGroups - 1; type >= 0; type--)
1693 if(npart[type] > 0)
1694 {
1695 char buf[MAXLEN_PATH];
1696 get_datagroup_name(type, buf);
1697 my_H5Gclose(hdf5_grp[type], buf);
1698 }
1699 my_H5Fclose(hdf5_file, fname);
1700 }
1701 }
1702}
1703
1715void IO_Def::distribute_file(int nfiles, int *filenr, int *master, int *last)
1716{
1717 int tasks_per_file = NTask / nfiles;
1718 int tasks_left = NTask % nfiles;
1719
1720 if(tasks_left == 0)
1721 {
1722 int group = ThisTask / tasks_per_file;
1723 *master = group * tasks_per_file;
1724 *last = (group + 1) * tasks_per_file - 1;
1725 *filenr = group;
1726 return;
1727 }
1728
1729 double tpf = ((double)NTask) / nfiles;
1730
1731 *last = -1;
1732 for(int i = 0; i < nfiles; i++)
1733 {
1734 *master = *last + 1;
1735 *last = (i + 1) * tpf;
1736 if(*last >= NTask)
1737 *last = *last - 1;
1738 if(*last < *master)
1739 Terminate("last < master");
1740 *filenr = i;
1741
1742 if(i == nfiles - 1)
1743 *last = NTask - 1;
1744
1745 if(ThisTask >= *master && ThisTask <= *last)
1746 return;
1747 }
1748}
1749
1759int IO_Def::get_bytes_per_memory_blockelement(int blocknr, int mode)
1760{
1761 if(blocknr < 0 || blocknr >= N_IO_Fields)
1762 Terminate("something is wrong here: blocknr=%d N_IO_Fields=%d", blocknr, N_IO_Fields);
1763
1764 int bytes_per_blockelement = 0;
1765
1766 IO_Field *field = &IO_Fields[blocknr];
1767
1768 switch(field->type_in_memory)
1769 {
1770 case MEM_INT:
1771 bytes_per_blockelement = field->values_per_block * sizeof(int);
1772 break;
1773 case MEM_INT64:
1774 bytes_per_blockelement = field->values_per_block * sizeof(long long);
1775 break;
1776 case MEM_MY_ID_TYPE:
1777 bytes_per_blockelement = field->values_per_block * sizeof(MyIDType);
1778 break;
1779 case MEM_MY_INTPOS_TYPE:
1780 bytes_per_blockelement = field->values_per_block * sizeof(MyIntPosType);
1781 break;
1782 case MEM_FLOAT:
1783 bytes_per_blockelement = field->values_per_block * sizeof(float);
1784 break;
1785 case MEM_DOUBLE:
1786 bytes_per_blockelement = field->values_per_block * sizeof(double);
1787 break;
1788 case MEM_MY_FLOAT:
1789 bytes_per_blockelement = field->values_per_block * sizeof(MyFloat);
1790 break;
1791 case MEM_MY_DOUBLE:
1792 bytes_per_blockelement = field->values_per_block * sizeof(MyDouble);
1793 break;
1794 case MEM_MY_FILEOFFSET:
1795 bytes_per_blockelement = field->values_per_block * sizeof(long long);
1796 break;
1797 }
1798
1799 return bytes_per_blockelement;
1800}
1801
1802hid_t IO_Def::get_hdf5_outputtype_of_block(int blocknr)
1803{
1804 hid_t hdf5_datatype = 0;
1805
1806 switch(IO_Fields[blocknr].type_in_file_output)
1807 {
1808 case FILE_INT:
1809 hdf5_datatype = H5T_NATIVE_INT;
1810 break;
1811 case FILE_INT64:
1812 hdf5_datatype = H5T_NATIVE_INT64;
1813 break;
1814 case FILE_MY_IO_FLOAT:
1815#ifdef OUTPUT_IN_DOUBLEPRECISION
1816 hdf5_datatype = H5T_NATIVE_DOUBLE;
1817#else
1818 hdf5_datatype = H5T_NATIVE_FLOAT;
1819#endif
1820 break;
1821 case FILE_MY_ID_TYPE:
1822#if defined(IDS_32BIT)
1823 hdf5_datatype = H5T_NATIVE_UINT32;
1824#elif defined(IDS_48BIT)
1825 hdf5_datatype = Int48_memtype;
1826#else
1827 hdf5_datatype = H5T_NATIVE_UINT64;
1828#endif
1829 break;
1830 case FILE_MY_INTPOS_TYPE:
1831#if defined(POSITIONS_IN_32BIT)
1832 hdf5_datatype = H5T_NATIVE_UINT32;
1833#elif defined(POSITIONS_IN_64BIT)
1834 hdf5_datatype = H5T_NATIVE_UINT64;
1835#else
1836 hdf5_datatype = Int128_memtype;
1837#endif
1838 break;
1839 case FILE_DOUBLE:
1840 hdf5_datatype = H5T_NATIVE_DOUBLE;
1841 break;
1842 case FILE_FLOAT:
1843 hdf5_datatype = H5T_NATIVE_FLOAT;
1844 break;
1845 case FILE_HALF:
1846 hdf5_datatype = Halfprec_memtype;
1847 break;
1848 case FILE_NONE:
1849 Terminate("undefined type");
1850 break;
1851 }
1852
1853 return hdf5_datatype;
1854}
1855
1856hid_t IO_Def::get_hdf5_memorytype_of_block(int blocknr)
1857{
1858 hid_t hdf5_datatype = 0;
1859
1860 switch(IO_Fields[blocknr].type_in_memory)
1861 {
1862 case MEM_INT:
1863 hdf5_datatype = H5T_NATIVE_INT;
1864 break;
1865 case MEM_INT64:
1866 hdf5_datatype = H5T_NATIVE_INT64;
1867 break;
1868 case MEM_MY_ID_TYPE:
1869#ifdef IDS_32BIT
1870 hdf5_datatype = H5T_NATIVE_UINT32;
1871#else
1872 hdf5_datatype = H5T_NATIVE_UINT64;
1873#endif
1874 break;
1875 case MEM_MY_INTPOS_TYPE:
1876#ifdef POSITIONS_IN_32BIT
1877 hdf5_datatype = H5T_NATIVE_UINT32;
1878#elif defined(POSITIONS_IN_64BIT)
1879 hdf5_datatype = H5T_NATIVE_UINT64;
1880#else
1881 hdf5_datatype = Int128_memtype;
1882#endif
1883 break;
1884 case MEM_FLOAT:
1885 hdf5_datatype = H5T_NATIVE_FLOAT;
1886 break;
1887 case MEM_DOUBLE:
1888 hdf5_datatype = H5T_NATIVE_DOUBLE;
1889 break;
1890 case MEM_MY_FLOAT:
1891 hdf5_datatype = H5T_NATIVE_MYFLOAT;
1892 break;
1893 case MEM_MY_DOUBLE:
1894 hdf5_datatype = H5T_NATIVE_MYDOUBLE;
1895 break;
1896 case MEM_MY_FILEOFFSET:
1897 hdf5_datatype = H5T_NATIVE_INT64;
1898 break;
1899 }
1900
1901 return hdf5_datatype;
1902}
1903
1912int IO_Def::get_values_per_blockelement(int blocknr)
1913{
1914 if(blocknr < 0 || blocknr >= N_IO_Fields)
1915 Terminate("something is wrong here: blocknr=%d N_IO_Fields=%d", blocknr, N_IO_Fields);
1916
1917 return IO_Fields[blocknr].values_per_block;
1918}
1919
1930long long IO_Def::get_particles_in_block(int blocknr, long long *npart_file, int *typelist)
1931{
1932 long long npart = 0;
1933
1934 for(int i = 0; i < N_DataGroups; i++)
1935 typelist[i] = 0;
1936
1937 switch(IO_Fields[blocknr].typelist)
1938 {
1939 case MASS_BLOCK:
1940 for(int i = 0; i < NTYPES; i++)
1941 {
1942 typelist[i] = (All.MassTable[i] == 0 && npart_file[i] > 0);
1943 npart += npart_file[i] * typelist[i];
1944 }
1945 return npart;
1946 break;
1947
1948 case AGE_BLOCK:
1949 for(int i = 0; i < NTYPES; i++)
1950 {
1951 typelist[i] = (npart_file[i] > 0);
1952
1953 if((file_format == FILEFORMAT_HDF5 && (i == 0 || i == 1 || i == 5)) || (file_format != FILEFORMAT_HDF5 && i != 4))
1954 typelist[i] = 0;
1955
1956 npart += npart_file[i] * typelist[i];
1957 }
1958 return npart;
1959 break;
1960
1961 case Z_BLOCK:
1962 for(int i = 0; i < NTYPES; i++)
1963 {
1964 typelist[i] = (npart_file[i] > 0);
1965 if((file_format == FILEFORMAT_HDF5 && (i == 1 || i == 5)) || (file_format != FILEFORMAT_HDF5 && (i != 0 && i != 4)))
1966 typelist[i] = 0;
1967
1968 npart += npart_file[i] * typelist[i];
1969 }
1970 return npart;
1971 break;
1972
1973 case GROUPS:
1974 npart = npart_file[0];
1975 typelist[0] = 1;
1976 return npart;
1977 break;
1978
1979 case SUBGROUPS:
1980 npart = npart_file[1];
1981 typelist[1] = 1;
1982 return npart;
1983 break;
1984
1985 case ID_BLOCK:
1986 npart = npart_file[2];
1987 typelist[2] = 1;
1988 return npart;
1989 break;
1990
1991 case CURRSUBS:
1992 case PREVSUBS:
1993 case TREELINK:
1994 case TREELENGTH:
1995 case MASSMAPS:
1996 case GALSNAPS:
1997 npart = npart_file[0];
1998 typelist[0] = 1;
1999 return npart;
2000 break;
2001
2002 case TREEHALOS:
2003 npart = npart_file[1];
2004 typelist[1] = 1;
2005 return npart;
2006 break;
2007
2008 case TREETIMES:
2009 npart = npart_file[2];
2010 typelist[2] = 1;
2011 return npart;
2012 break;
2013
2014 case TREETABLE:
2015 npart = npart_file[NTYPES];
2016 typelist[NTYPES] = 1;
2017 return npart;
2018 break;
2019
2020 case HEALPIXTAB:
2021 npart = npart_file[NTYPES + 1];
2022 typelist[NTYPES + 1] = 1;
2023 return npart;
2024 break;
2025
2026 default:
2027 for(int i = 0; i < N_DataGroups; i++)
2028 {
2029 if((IO_Fields[blocknr].typelist & (1 << i)) && npart_file[i] > 0)
2030 {
2031 typelist[i] = 1;
2032 npart += npart_file[i];
2033 }
2034 else
2035 typelist[i] = 0;
2036 }
2037 return npart;
2038 break;
2039 }
2040
2041 return 0;
2042}
2043
2051void IO_Def::get_Tab_IO_Label(int blocknr, char *label)
2052{
2053 if(blocknr < 0 || blocknr >= N_IO_Fields)
2054 Terminate("something is wrong here: blocknr=%d N_IO_Fields=%d", blocknr, N_IO_Fields);
2055
2056 strcpy(label, IO_Fields[blocknr].label);
2057}
2058
2066void IO_Def::get_dataset_name(int blocknr, char *buf)
2067{
2068 if(blocknr < 0 || blocknr >= N_IO_Fields)
2069 Terminate("something is wrong here: blocknr=%d N_IO_Fields=%d", blocknr, N_IO_Fields);
2070
2071 strcpy(buf, IO_Fields[blocknr].datasetname);
2072}
2073
2074void IO_Def::write_dataset_attributes(hid_t hdf5_dataset, int blocknr)
2075{
2076 if(blocknr < 0 || blocknr >= N_IO_Fields)
2077 Terminate("something is wrong here: blocknr=%d N_IO_Fields=%d", blocknr, N_IO_Fields);
2078
2079 if(IO_Fields[blocknr].hasunit == 0)
2080 return;
2081
2082 if(All.ComovingIntegrationOn)
2083 {
2084 write_scalar_attribute(hdf5_dataset, "a_scaling", &IO_Fields[blocknr].a, H5T_NATIVE_DOUBLE);
2085 write_scalar_attribute(hdf5_dataset, "h_scaling", &IO_Fields[blocknr].h, H5T_NATIVE_DOUBLE);
2086 }
2087 else
2088 {
2089 double zero = 0;
2090 write_scalar_attribute(hdf5_dataset, "a_scaling", &zero, H5T_NATIVE_DOUBLE);
2091 write_scalar_attribute(hdf5_dataset, "h_scaling", &zero, H5T_NATIVE_DOUBLE);
2092 }
2093
2094 write_scalar_attribute(hdf5_dataset, "length_scaling", &IO_Fields[blocknr].L, H5T_NATIVE_DOUBLE);
2095 write_scalar_attribute(hdf5_dataset, "mass_scaling", &IO_Fields[blocknr].M, H5T_NATIVE_DOUBLE);
2096 write_scalar_attribute(hdf5_dataset, "velocity_scaling", &IO_Fields[blocknr].V, H5T_NATIVE_DOUBLE);
2097
2098 write_scalar_attribute(hdf5_dataset, "to_cgs", &IO_Fields[blocknr].c, H5T_NATIVE_DOUBLE);
2099}
2100
2108void IO_Def::write_parameters_attributes_in_hdf5(hid_t handle)
2109{
2110 for(int i = 0; i < All.NParameters; i++)
2111 {
2112 switch(All.ParametersType[i])
2113 {
2114 case PARAM_DOUBLE:
2115 write_scalar_attribute(handle, All.ParametersTag[i], All.ParametersValue[i], H5T_NATIVE_DOUBLE);
2116 break;
2117 case PARAM_STRING:
2118 write_string_attribute(handle, All.ParametersTag[i], (const char *)All.ParametersValue[i]);
2119 break;
2120 case PARAM_INT:
2121 write_scalar_attribute(handle, All.ParametersTag[i], All.ParametersValue[i], H5T_NATIVE_INT);
2122 break;
2123 }
2124 }
2125}
2126
2127/*---------------------- Routine find a block in a snapfile -------------------*/
2128void IO_Def::find_block(char *label, FILE *fd)
2129{
2130 unsigned int blocksize = 0, blksize;
2131 char blocklabel[5] = {" "};
2132
2133#define FBSKIP \
2134 { \
2135 my_fread(&blksize, sizeof(int), 1, fd); \
2136 }
2137
2138 rewind(fd);
2139
2140 while(!feof(fd) && blocksize == 0)
2141 {
2142 FBSKIP;
2143
2144 if(blksize != 8)
2145 {
2146 Terminate("Incorrect Format (blksize=%u)!\n", blksize);
2147 }
2148 else
2149 {
2150 my_fread(blocklabel, LABEL_LEN * sizeof(char), 1, fd);
2151 my_fread(&blocksize, sizeof(int), 1, fd);
2152
2153 FBSKIP;
2154
2155 if(strncmp(label, blocklabel, LABEL_LEN) != 0)
2156 {
2157 fseek(fd, blocksize, 1);
2158 blocksize = 0;
2159 }
2160 }
2161 }
2162 if(feof(fd))
2163 Terminate("Block '%c%c%c%c' not found !\n", label[0], label[1], label[2], label[3]);
2164}
2165
2166void IO_Def::read_segment(const char *fname, int type, long long offset, long long count, int num_files)
2167{
2168 long long nleft = count;
2169 long long offleft = offset;
2170 long long nskip = 0;
2171
2172 for(int filenr = 0; filenr < num_files && nleft > 0; filenr++)
2173 {
2174 long long nloc = ntype_in_files[filenr * N_DataGroups + type];
2175
2176 if(nloc > offleft) // we may have something in this file
2177 {
2178 long long nread;
2179
2180 if(nloc - offleft > nleft) // there are more particles in the file then we need
2181 nread = nleft;
2182 else
2183 nread = nloc - offleft;
2184
2185 /* now read partial list */
2186 read_single_file_segment(fname, filenr, type, offleft, nread, nskip, num_files);
2187
2188 nleft -= nread;
2189 nskip += nread;
2190 offleft += nread;
2191 }
2192
2193 offleft -= nloc;
2194 }
2195
2196 if(nleft > 0)
2197 {
2198 for(int filenr = 0; filenr < num_files; filenr++)
2199 {
2200 long long nloc = ntype_in_files[filenr * N_DataGroups + type];
2201 printf("filenr=%d: nloc=%lld\n", filenr, nloc);
2202 }
2203 Terminate("Not all desired entries read: nleft=%lld type=%d\n", nleft, type);
2204 }
2205}
2206
2207void IO_Def::read_single_file_segment(const char *basename, int filenr, int type, long long offset, unsigned long long count,
2208 long long storage_offset, int num_files)
2209{
2210 int bytes_per_blockelement_in_file = 0;
2211 hid_t hdf5_file = 0, hdf5_grp = 0, hdf5_dataspace_in_file;
2212 hid_t hdf5_dataspace_in_memory, hdf5_dataset;
2213 FILE *fd = 0;
2214 char fname[MAXLEN_PATH];
2215
2216 if(num_files > 1)
2217 {
2218 if(file_format == FILEFORMAT_HDF5)
2219 sprintf(fname, "%s.%d.hdf5", basename, filenr);
2220 else
2221 sprintf(fname, "%s.%d", basename, filenr);
2222 }
2223 else
2224 {
2225 if(file_format == FILEFORMAT_HDF5)
2226 sprintf(fname, "%s.hdf5", basename);
2227 else
2228 sprintf(fname, "%s", basename);
2229 }
2230
2231 /* open file */
2232 if(file_format == FILEFORMAT_HDF5)
2233 {
2234 hdf5_file = my_H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
2235 char buf[MAXLEN_PATH];
2236 get_datagroup_name(type, buf);
2237 hdf5_grp = my_H5Gopen(hdf5_file, buf);
2238 }
2239 else if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
2240 {
2241 if(!(fd = fopen(fname, "r")))
2242 Terminate("can't open file `%s' for reading initial conditions.\n", fname);
2243
2244 unsigned int blksize1, blksize2;
2245
2246 if(file_format == FILEFORMAT_LEGACY2)
2247 {
2248 int nextblock;
2249 char label[LABEL_LEN];
2250 my_fread(&blksize1, sizeof(int), 1, fd);
2251 my_fread(&label, sizeof(char), LABEL_LEN, fd);
2252 my_fread(&nextblock, sizeof(int), 1, fd);
2253 my_fread(&blksize2, sizeof(int), 1, fd);
2254 }
2255
2256 my_fread(&blksize1, sizeof(int), 1, fd);
2257 my_fread(header_buf, header_size, 1, fd);
2258 my_fread(&blksize2, sizeof(int), 1, fd);
2259 }
2260 else
2261 Terminate("unknown ICFormat");
2262
2263 long long npart[N_DataGroups];
2264 for(int i = 0; i < N_DataGroups; i++)
2265 npart[i] = ntype_in_files[filenr * N_DataGroups + i];
2266
2267 for(int blocknr = 0; blocknr < N_IO_Fields; blocknr++)
2268 {
2269 if(IO_Fields[blocknr].type_in_memory != MEM_MY_FILEOFFSET)
2270 {
2271 unsigned int blksize1, blksize2;
2272 int typelist[N_DataGroups];
2273 int bytes_per_blockelement = get_bytes_per_memory_blockelement(blocknr, 1);
2274 long long npart_in_block = get_particles_in_block(blocknr, npart, &typelist[0]);
2275 hid_t hdf5_memory_datatype = get_hdf5_memorytype_of_block(blocknr);
2276 char dname[MAXLEN_PATH];
2277 get_dataset_name(blocknr, dname);
2278
2279 if(npart_in_block > 0 && typelist[type] > 0 && IO_Fields[blocknr].read_flag != SKIP_ON_READ)
2280 {
2281 if(file_format == FILEFORMAT_LEGACY2)
2282 {
2283 char expected_label[LABEL_LEN + 1];
2284 get_Tab_IO_Label(blocknr, expected_label);
2285 find_block(expected_label, fd);
2286 }
2287
2288 if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
2289 {
2290 my_fread(&blksize1, sizeof(int), 1, fd);
2291 bytes_per_blockelement_in_file = blksize1 / npart_in_block;
2292 }
2293
2294 void *CommBuffer = (char *)Mem.mymalloc("CommBuffer", bytes_per_blockelement * count);
2295
2296 if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
2297 {
2298 fseek(fd, bytes_per_blockelement_in_file * offset, SEEK_CUR);
2299
2300 if(bytes_per_blockelement_in_file != bytes_per_blockelement)
2301 {
2302 char *CommAuxBuffer = (char *)Mem.mymalloc("CommAuxBuffer", bytes_per_blockelement_in_file * count);
2303 my_fread(CommAuxBuffer, bytes_per_blockelement_in_file, count, fd);
2304 type_cast_data((char *)CommAuxBuffer, bytes_per_blockelement_in_file, (char *)CommBuffer, bytes_per_blockelement,
2305 count, blocknr);
2306 Mem.myfree(CommAuxBuffer);
2307 }
2308 else
2309 my_fread(CommBuffer, bytes_per_blockelement, count, fd);
2310
2311 fseek(fd, bytes_per_blockelement_in_file * (npart_in_block - offset - count), SEEK_CUR);
2312
2313 my_fread(&blksize2, sizeof(int), 1, fd);
2314 if(blksize1 != blksize2)
2315 {
2316 char buf[MAXLEN_PATH];
2317 sprintf(buf, "incorrect block-sizes detected!\n Task=%d blocknr=%d blksize1=%d blksize2=%d\n", ThisTask,
2318 blocknr, blksize1, blksize2);
2319 if(blocknr == 2) /* block number 2 is always IDs */
2320 strcat(buf, "Possible mismatch of 32bit and 64bit ID's in IC file and GADGET compilation !\n");
2321 Terminate(buf);
2322 }
2323 }
2324
2325 if(file_format == FILEFORMAT_HDF5)
2326 {
2327 hdf5_dataset = my_H5Dopen_if_existing(hdf5_grp, dname);
2328 int rank;
2329 hsize_t dims[2], nelem[2], start[2];
2330
2331 dims[0] = npart[type];
2332 dims[1] = get_values_per_blockelement(blocknr);
2333 if(dims[1] == 1)
2334 rank = 1;
2335 else
2336 rank = 2;
2337
2338 hdf5_dataspace_in_file = my_H5Screate_simple(rank, dims, NULL);
2339
2340 dims[0] = count;
2341 hdf5_dataspace_in_memory = my_H5Screate_simple(rank, dims, NULL);
2342
2343 start[0] = offset;
2344 start[1] = 0;
2345
2346 nelem[0] = count;
2347 nelem[1] = get_values_per_blockelement(blocknr);
2348
2349 my_H5Sselect_hyperslab(hdf5_dataspace_in_file, H5S_SELECT_SET, start, NULL, nelem, NULL);
2350
2351 // Test if dataset was present
2352 if(hdf5_dataset < 0)
2353 {
2354 // no, pad with zeros
2355 printf("\tDataset %s not present for particle type %d, using zero.\n", dname, type);
2356 memset(CommBuffer, 0, dims[0] * dims[1] * my_H5Tget_size(hdf5_memory_datatype));
2357 }
2358 else
2359 {
2360 hid_t hdf5_file_datatype = H5Dget_type(hdf5_dataset);
2361 byte_count += dims[0] * dims[1] * my_H5Tget_size(hdf5_file_datatype); /* for I/O performance measurement */
2362 H5Tclose(hdf5_file_datatype);
2363
2364 my_H5Dread(hdf5_dataset, hdf5_memory_datatype, hdf5_dataspace_in_memory, hdf5_dataspace_in_file, H5P_DEFAULT,
2365 CommBuffer, dname);
2366 my_H5Dclose(hdf5_dataset, dname);
2367 }
2368 my_H5Sclose(hdf5_dataspace_in_memory, H5S_SIMPLE);
2369 my_H5Sclose(hdf5_dataspace_in_file, H5S_SIMPLE);
2370 }
2371
2372 empty_read_buffer(blocknr, storage_offset, count, type, 0, CommBuffer);
2373
2374 Mem.myfree(CommBuffer);
2375 }
2376 else
2377 {
2378 if(file_format == FILEFORMAT_LEGACY1 && npart_in_block > 0)
2379 {
2380 my_fread(&blksize1, sizeof(int), 1, fd);
2381 bytes_per_blockelement_in_file = blksize1 / npart_in_block;
2382 fseek(fd, bytes_per_blockelement_in_file * npart_in_block, SEEK_CUR);
2383 my_fread(&blksize2, sizeof(int), 1, fd);
2384 if(blksize1 != blksize2)
2385 {
2386 char buf[MAXLEN_PATH];
2387 sprintf(buf, "incorrect block-sizes detected!\n Task=%d blocknr=%d blksize1=%d blksize2=%d\n", ThisTask,
2388 blocknr, blksize1, blksize2);
2389 if(blocknr == 2) /* block number 2 is always IDs */
2390 strcat(buf, "Possible mismatch of 32bit and 64bit ID's in IC file and GADGET compilation !\n");
2391 Terminate(buf);
2392 }
2393 }
2394 }
2395 }
2396 }
2397
2398 if(file_format == FILEFORMAT_LEGACY1 || file_format == FILEFORMAT_LEGACY2)
2399 fclose(fd);
2400
2401 if(file_format == FILEFORMAT_HDF5)
2402 {
2403 char buf[MAXLEN_PATH];
2404 get_datagroup_name(type, buf);
2405 my_H5Gclose(hdf5_grp, buf);
2406 my_H5Fclose(hdf5_file, fname);
2407 }
2408
2409 read_increase_numbers(type, count);
2410}
2411
2412void IO_Def::rename_file_to_bak_if_it_exists(char *fname)
2413{
2414 char fin[MAXLEN_PATH], buf[2 * MAXLEN_PATH];
2415
2416 strcpy(fin, fname);
2417
2418 char *p = strrchr(fin, '/');
2419 if(p)
2420 {
2421 *p = 0;
2422 sprintf(buf, "%s/bak-%s", fin, p + 1);
2423 }
2424 else
2425 sprintf(buf, "bak-%s", fname);
2426
2427 if(FILE *fcheck = fopen(fname, "r")) // check if file already exists, if yes, try to rename the existing file
2428 {
2429 fclose(fcheck);
2430
2431 if(!fopen(buf, "r")) // only do the rename of the old file if the back-up file doesn't exist yet
2432 {
2433 mpi_printf("%s rename '%s' to '%s'\n", info, fname, buf);
2434 rename(fname, buf);
2435 }
2436 }
2437}
2438
2439void IO_Def::alloc_and_read_ntype_in_files(const char *fname, int num_files)
2440{
2441 ntype_in_files = (long long *)Mem.mymalloc_movable(&ntype_in_files, "ntype_in_files", num_files * N_DataGroups * sizeof(long long));
2442
2443 for(int filenr = 0; filenr < num_files; filenr++)
2444 {
2445 char buf[3 * MAXLEN_PATH];
2446
2447 if(num_files > 1)
2448 {
2449 sprintf(buf, "%s.%d", fname, filenr);
2450 if(file_format == 3)
2451 sprintf(buf, "%s.%d.hdf5", fname, filenr);
2452 }
2453 else
2454 {
2455 sprintf(buf, "%s", fname);
2456 if(file_format == 3)
2457 sprintf(buf, "%s.hdf5", fname);
2458 }
2459
2460 if(file_format == 3)
2461 {
2462 read_header_fields(buf);
2463 }
2464 else if(file_format == 1 || file_format == 2)
2465 {
2466 FILE *fd = 0;
2467
2468 if(!(fd = fopen(buf, "r")))
2469 Terminate("can't open file `%s' for reading initial conditions.\n", fname);
2470
2471 int blksize1, blksize2;
2472
2473 if(file_format == 2)
2474 {
2475 char label[4];
2476 int nextblock;
2477 my_fread(&blksize1, sizeof(int), 1, fd);
2478 my_fread(&label, sizeof(char), 4, fd);
2479 my_fread(&nextblock, sizeof(int), 1, fd);
2480 printf("READIC: Reading header => '%c%c%c%c' (%d byte)\n", label[0], label[1], label[2], label[3], nextblock);
2481 my_fread(&blksize2, sizeof(int), 1, fd);
2482 }
2483
2484 my_fread(&blksize1, sizeof(int), 1, fd);
2485 my_fread(header_buf, header_size, 1, fd);
2486 my_fread(&blksize2, sizeof(int), 1, fd);
2487
2488 if(blksize1 != blksize2)
2489 Terminate("incorrect header format, blksize1=%d blksize2=%d header_size=%d\n", blksize1, blksize2, (int)header_size);
2490
2491 fclose(fd);
2492 }
2493
2494 long long n_type[N_DataGroups], npart[N_DataGroups];
2495
2496 read_file_header(fname, filenr, 0, 0, n_type, npart, NULL);
2497
2498 for(int type = 0; type < N_DataGroups; type++)
2499 ntype_in_files[filenr * N_DataGroups + type] = npart[type];
2500 }
2501}
All
global_data_all_processes All
Definition: main.cc:40
IO_Def
Definition: io.h:129
IO_Def::get_filenr_from_header
virtual int get_filenr_from_header(void)=0
IO_Def::N_DataGroups
int N_DataGroups
Definition: io.h:139
IO_Def::alloc_and_read_ntype_in_files
void alloc_and_read_ntype_in_files(const char *fname, int num_files)
Definition: io.cc:2439
IO_Def::write_multiple_files
void write_multiple_files(char *fname, int numfilesperdump, int append_flag=0, int chunk_size=0)
Definition: io.cc:700
IO_Def::ntype_in_files
long long * ntype_in_files
Definition: io.h:176
IO_Def::Max_IO_Fields
int Max_IO_Fields
Definition: io.h:141
IO_Def::~IO_Def
virtual ~IO_Def()
Definition: io.cc:47
IO_Def::header_buf
void * header_buf
Definition: io.h:174
IO_Def::get_base_address_of_structure
virtual void * get_base_address_of_structure(enum arrays array, int index)=0
IO_Def::read_segment
void read_segment(const char *fname, int type, long long offset, long long count, int numfiles)
Definition: io.cc:2166
IO_Def::set_type_of_element
virtual void set_type_of_element(int index, int type)=0
IO_Def::read_increase_numbers
virtual void read_increase_numbers(int type, int n_for_this_task)=0
IO_Def::fill_file_header
virtual void fill_file_header(int writeTask, int lastTask, long long *nloc_part, long long *npart)=0
IO_Def::header_size
size_t header_size
Definition: io.h:173
IO_Def::init_field
void init_field(const char *label, const char *datasetname, enum types_in_memory type_in_memory, enum types_in_file type_in_file_output, enum read_flags read_flag, int values_per_block, enum arrays array, void *pointer_to_field, void(*io_func)(IO_Def *, int, int, void *, int), int typelist_bitmask, int hasunits, double a, double h, double L, double M, double V, double c, bool compression_on=false)
Definition: io.cc:66
IO_Def::N_IO_Fields
int N_IO_Fields
Definition: io.h:140
IO_Def::read_header_fields
virtual void read_header_fields(const char *fname)=0
IO_Def::read_file_header
virtual void read_file_header(const char *fname, int filenr, int readTask, int lastTask, long long *nloc_part, long long *npart, int *nstart)=0
IO_Def::info
char info[100]
Definition: io.h:177
IO_Def::set_filenr_in_header
virtual void set_filenr_in_header(int)=0
IO_Def::type_of_file
enum file_contents type_of_file
Definition: io.h:187
IO_Def::write_header_fields
virtual void write_header_fields(hid_t)=0
IO_Def::read_files_driver
void read_files_driver(const char *fname, int rep, int numfiles)
Definition: io.cc:227
IO_Def::get_datagroup_name
virtual void get_datagroup_name(int grnr, char *gname)=0
IO_Def::read_single_file_segment
void read_single_file_segment(const char *fname, int filenr, int type, long long offset, unsigned long long count, long long storage_offset, int numfiles)
Definition: io.cc:2207
IO_Def::get_type_of_element
virtual int get_type_of_element(int index)=0
IO_Def::find_files
int find_files(const char *fname, const char *fname_multiple)
This function determines on how many files a given snapshot or group/desc catalogue is distributed.
Definition: io.cc:132
IO_Def::write_compile_time_options_in_hdf5
void write_compile_time_options_in_hdf5(hid_t handle)
half_float::half
Definition: half.hpp:1336
io_streamcount::byte_count
long long byte_count
Definition: io_streamcount.h:20
io_streamcount::my_fread
size_t my_fread(void *ptr, size_t size, size_t nmemb, FILE *stream)
A wrapper for the fread() function.
Definition: io_streamcount.h:66
io_streamcount::my_fwrite
size_t my_fwrite(const void *ptr, size_t size, size_t nmemb, FILE *stream)
A wrapper for the fwrite() function.
Definition: io_streamcount.h:33
parameters::ParametersType
char ParametersType[MAX_PARAMETERS]
Definition: parameters.h:46
parameters::ParametersTag
char ParametersTag[MAX_PARAMETERS][MAXLEN_PARAM_TAG]
Definition: parameters.h:44
parameters::NParameters
int NParameters
Definition: parameters.h:42
parameters::ParametersValue
void * ParametersValue[MAX_PARAMETERS]
Definition: parameters.h:45
setcomm::ThisNode
int ThisNode
Definition: setcomm.h:36
setcomm::mpi_printf
void mpi_printf(const char *fmt,...)
Definition: setcomm.h:55
setcomm::ThisTask
int ThisTask
Definition: setcomm.h:33
setcomm::NTask
int NTask
Definition: setcomm.h:32
setcomm::RankInThisNode
int RankInThisNode
Definition: setcomm.h:39
setcomm::Communicator
MPI_Comm Communicator
Definition: setcomm.h:31
FILEFORMAT_LEGACY2
#define FILEFORMAT_LEGACY2
Definition: constants.h:18
FILEFORMAT_LEGACY1
#define FILEFORMAT_LEGACY1
Definition: constants.h:17
COMMBUFFERSIZE
#define COMMBUFFERSIZE
Definition: constants.h:47
NTYPES
#define NTYPES
Definition: constants.h:308
FILEFORMAT_HDF5
#define FILEFORMAT_HDF5
Definition: constants.h:19
MAXLEN_PATH
#define MAXLEN_PATH
Definition: constants.h:300
H5T_NATIVE_MYDOUBLE
#define H5T_NATIVE_MYDOUBLE
Definition: dtypes.h:92
MyDouble
float MyDouble
Definition: dtypes.h:87
MyFloat
float MyFloat
Definition: dtypes.h:86
MyIntPosType
uint32_t MyIntPosType
Definition: dtypes.h:35
H5T_NATIVE_MYFLOAT
#define H5T_NATIVE_MYFLOAT
Definition: dtypes.h:91
RST_BEGIN
@ RST_BEGIN
Definition: dtypes.h:313
MyIDType
unsigned int MyIDType
Definition: dtypes.h:68
my_H5Gopen
hid_t my_H5Gopen(hid_t loc_id, const char *groupname)
Definition: hdf5_util.cc:311
my_H5Sselect_hyperslab
herr_t my_H5Sselect_hyperslab(hid_t space_id, H5S_seloper_t op, const hsize_t *start, const hsize_t *stride, const hsize_t *count, const hsize_t *block)
Definition: hdf5_util.cc:549
my_H5Dread
herr_t my_H5Dread(hid_t dataset_id, hid_t mem_type_id, hid_t mem_space_id, hid_t file_space_id, hid_t xfer_plist_id, void *buf, const char *datasetname)
Definition: hdf5_util.cc:385
my_H5Fopen
hid_t my_H5Fopen(const char *fname, unsigned int flags, hid_t fapl_id)
Definition: hdf5_util.cc:297
my_H5Gclose
herr_t my_H5Gclose(hid_t group_id, const char *groupname)
Definition: hdf5_util.cc:457
Halfprec_memtype
hid_t Halfprec_memtype
Definition: hdf5_util.cc:29
my_H5Tget_size
size_t my_H5Tget_size(hid_t datatype_id)
Definition: hdf5_util.cc:564
my_H5Dset_extent
herr_t my_H5Dset_extent(hid_t dset_id, const hsize_t size[])
Definition: hdf5_util.cc:336
Int48_memtype
hid_t Int48_memtype
Definition: hdf5_util.cc:30
my_H5Gcreate
hid_t my_H5Gcreate(hid_t loc_id, const char *groupname, size_t size_hint)
Definition: hdf5_util.cc:174
my_H5Pclose
herr_t my_H5Pclose(hid_t plist)
Definition: hdf5_util.cc:468
my_H5Screate_simple
hid_t my_H5Screate_simple(int rank, const hsize_t *current_dims, const hsize_t *maximum_dims)
Definition: hdf5_util.cc:283
my_H5Fclose
herr_t my_H5Fclose(hid_t file_id, const char *fname)
Definition: hdf5_util.cc:482
my_H5Dclose
herr_t my_H5Dclose(hid_t dataset_id, const char *datasetname)
Definition: hdf5_util.cc:443
my_H5Fcreate
hid_t my_H5Fcreate(const char *fname, unsigned int flags, hid_t fcpl_id, hid_t fapl_id)
Definition: hdf5_util.cc:158
my_H5Dopen_if_existing
hid_t my_H5Dopen_if_existing(hid_t file_id, const char *datasetname)
Definition: hdf5_util.cc:352
my_H5Dcreate
hid_t my_H5Dcreate(hid_t loc_id, const char *datasetname, hid_t type_id, hid_t space_id, hid_t dcpl_id)
Definition: hdf5_util.cc:190
my_H5Sclose
herr_t my_H5Sclose(hid_t dataspace_id, H5S_class_t type)
Definition: hdf5_util.cc:496
Int128_memtype
hid_t Int128_memtype
Definition: hdf5_util.cc:31
write_scalar_attribute
void write_scalar_attribute(hid_t handle, const char *attr_name, const void *buf, hid_t mem_type_id)
Definition: hdf5_util.cc:621
my_H5Dwrite
herr_t my_H5Dwrite(hid_t dataset_id, hid_t mem_type_id, hid_t mem_space_id, hid_t file_space_id, hid_t xfer_plist_id, const void *buf, const char *datasetname)
Definition: hdf5_util.cc:206
write_string_attribute
void write_string_attribute(hid_t handle, const char *attr_name, const char *buf)
Definition: hdf5_util.cc:654
COMPRESSION_CHUNKSIZE
#define COMPRESSION_CHUNKSIZE
Definition: hdf5_util.h:17
SKIP
#define SKIP
FBSKIP
#define FBSKIP
read_flags
read_flags
Definition: io.h:123
SKIP_ON_READ
@ SKIP_ON_READ
Definition: io.h:125
GALSNAPS
@ GALSNAPS
Definition: io.h:118
MASSMAPS
@ MASSMAPS
Definition: io.h:115
TREELENGTH
@ TREELENGTH
Definition: io.h:112
PREVSUBS
@ PREVSUBS
Definition: io.h:110
AGE_BLOCK
@ AGE_BLOCK
Definition: io.h:105
TREETIMES
@ TREETIMES
Definition: io.h:117
TREELINK
@ TREELINK
Definition: io.h:114
Z_BLOCK
@ Z_BLOCK
Definition: io.h:106
GROUPS
@ GROUPS
Definition: io.h:107
MASS_BLOCK
@ MASS_BLOCK
Definition: io.h:104
CURRSUBS
@ CURRSUBS
Definition: io.h:111
SUBGROUPS
@ SUBGROUPS
Definition: io.h:108
ID_BLOCK
@ ID_BLOCK
Definition: io.h:109
TREEHALOS
@ TREEHALOS
Definition: io.h:113
TREETABLE
@ TREETABLE
Definition: io.h:116
HEALPIXTAB
@ HEALPIXTAB
Definition: io.h:119
arrays
arrays
Definition: io.h:30
A_MTRP
@ A_MTRP
Definition: io.h:40
A_NONE
@ A_NONE
Definition: io.h:31
A_P
@ A_P
Definition: io.h:33
types_in_file
types_in_file
Definition: io.h:65
FILE_FLOAT
@ FILE_FLOAT
Definition: io.h:74
FILE_MY_INTPOS_TYPE
@ FILE_MY_INTPOS_TYPE
Definition: io.h:72
FILE_DOUBLE
@ FILE_DOUBLE
Definition: io.h:73
FILE_HALF
@ FILE_HALF
Definition: io.h:70
FILE_NONE
@ FILE_NONE
Definition: io.h:66
FILE_MY_IO_FLOAT
@ FILE_MY_IO_FLOAT
Definition: io.h:69
FILE_INT64
@ FILE_INT64
Definition: io.h:68
FILE_MY_ID_TYPE
@ FILE_MY_ID_TYPE
Definition: io.h:71
FILE_INT
@ FILE_INT
Definition: io.h:67
LABEL_LEN
#define LABEL_LEN
Definition: io.h:26
DATASETNAME_LEN
#define DATASETNAME_LEN
Definition: io.h:27
types_in_memory
types_in_memory
Definition: io.h:78
MEM_INT
@ MEM_INT
Definition: io.h:79
MEM_DOUBLE
@ MEM_DOUBLE
Definition: io.h:84
MEM_MY_FLOAT
@ MEM_MY_FLOAT
Definition: io.h:85
MEM_MY_FILEOFFSET
@ MEM_MY_FILEOFFSET
Definition: io.h:87
MEM_MY_ID_TYPE
@ MEM_MY_ID_TYPE
Definition: io.h:81
MEM_MY_INTPOS_TYPE
@ MEM_MY_INTPOS_TYPE
Definition: io.h:82
MEM_FLOAT
@ MEM_FLOAT
Definition: io.h:83
MEM_MY_DOUBLE
@ MEM_MY_DOUBLE
Definition: io.h:86
MEM_INT64
@ MEM_INT64
Definition: io.h:80
FILE_IS_SNAPSHOT
@ FILE_IS_SNAPSHOT
Definition: io.h:53
FILE_IS_LIGHTCONE
@ FILE_IS_LIGHTCONE
Definition: io.h:59
Terminate
#define Terminate(...)
Definition: macros.h:15
TAG_HEADER
#define TAG_HEADER
Definition: mpi_utils.h:26
TAG_NFORTHISTASK
#define TAG_NFORTHISTASK
Definition: mpi_utils.h:39
TAG_PDATA
#define TAG_PDATA
Definition: mpi_utils.h:27
TAG_N
#define TAG_N
Definition: mpi_utils.h:25
TAG_KEY
#define TAG_KEY
Definition: mpi_utils.h:29
Mem
memory Mem
Definition: main.cc:44
PARAM_INT
#define PARAM_INT
Definition: parameters.h:20
PARAM_STRING
#define PARAM_STRING
Definition: parameters.h:19
PARAM_DOUBLE
#define PARAM_DOUBLE
Definition: parameters.h:18
global_data_all_processes::MaxFilesWithConcurrentIO
int MaxFilesWithConcurrentIO
Definition: allvars.h:49
global_data_all_processes::RestartFlag
enum restart_options RestartFlag
Definition: allvars.h:68
global_data_all_processes::MassTable
double MassTable[NTYPES]
Definition: allvars.h:269
global_data_all_processes::ComovingIntegrationOn
int ComovingIntegrationOn
Definition: allvars.h:148
myflush
void myflush(FILE *fstream)
Definition: system.cc:125