parallel_io.c 45.7 KB
Newer Older
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
4
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk).
5
 *
6
7
8
9
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
10
 *
11
12
13
14
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
15
 *
16
17
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
18
 *
19
20
21
22
23
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

24
#if defined(HAVE_HDF5) && defined(WITH_MPI) && defined(HAVE_PARALLEL_HDF5)
25
26
27
28
29

/* Some standard headers. */
#include <hdf5.h>
#include <math.h>
#include <mpi.h>
30
31
32
33
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
34

35
36
37
38
/* This object's header. */
#include "parallel_io.h"

/* Local includes. */
39
#include "chemistry_io.h"
40
#include "common_io.h"
41
#include "cooling.h"
42
#include "dimension.h"
43
#include "engine.h"
44
#include "error.h"
45
#include "gravity_io.h"
46
#include "gravity_properties.h"
47
#include "hydro_io.h"
48
#include "hydro_properties.h"
49
#include "io_properties.h"
50
51
#include "kernel_hydro.h"
#include "part.h"
52
#include "stars_io.h"
53
#include "units.h"
54
#include "xmf.h"
55

56
/* The current limit of ROMIO (the underlying MPI-IO layer) is 2GB */
57
#define HDF5_PARALLEL_IO_MAX_BYTES 2000000000LL
58

59
/* Are we timing the i/o? */
60
//#define IO_SPEED_MEASUREMENT
61

62
/**
63
 * @brief Reads a chunk of data from an open HDF5 dataset
64
 *
65
66
67
68
69
 * @param h_data The HDF5 dataset to write to.
 * @param h_plist_id the parallel HDF5 properties.
 * @param props The #io_props of the field to read.
 * @param N The number of particles to write.
 * @param offset Offset in the array where this mpi task starts writing.
70
 * @param internal_units The #unit_system used internally.
71
 * @param ic_units The #unit_system used in the snapshots.
72
73
 * @param cleanup_h Are we removing h-factors from the ICs?
 * @param h The value of the reduced Hubble constant to use for cleaning.
74
 */
75
void readArray_chunk(hid_t h_data, hid_t h_plist_id,
76
77
                     const struct io_props props, size_t N, long long offset,
                     const struct unit_system* internal_units,
78
79
                     const struct unit_system* ic_units, int cleanup_h,
                     double h) {
80

81
82
83
  const size_t typeSize = io_sizeof_type(props.type);
  const size_t copySize = typeSize * props.dimension;
  const size_t num_elements = N * props.dimension;
84

85
86
87
  /* Can't handle writes of more than 2GB */
  if (N * props.dimension * typeSize > HDF5_PARALLEL_IO_MAX_BYTES)
    error("Dataset too large to be written in one pass!");
88

89
  /* Allocate temporary buffer */
90
  void* temp = malloc(num_elements * typeSize);
91
  if (temp == NULL) error("Unable to allocate memory for temporary buffer");
92

93
  /* Prepare information for hyper-slab */
94
95
  hsize_t shape[2], offsets[2];
  int rank;
96
  if (props.dimension > 1) {
97
98
    rank = 2;
    shape[0] = N;
99
    shape[1] = props.dimension;
100
101
102
    offsets[0] = offset;
    offsets[1] = 0;
  } else {
103
    rank = 2;
104
    shape[0] = N;
105
    shape[1] = 1;
106
107
108
    offsets[0] = offset;
    offsets[1] = 0;
  }
109
110

  /* Create data space in memory */
111
  const hid_t h_memspace = H5Screate_simple(rank, shape, NULL);
112

113
  /* Select hyper-slab in file */
114
  const hid_t h_filespace = H5Dget_space(h_data);
115
116
  H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape, NULL);

117
118
119
  /* Read HDF5 dataspace in temporary buffer */
  /* Dirty version that happens to work for vectors but should be improved */
  /* Using HDF5 dataspaces would be better */
120
  const hid_t h_err = H5Dread(h_data, io_hdf5_type(props.type), h_memspace,
121
                              h_filespace, h_plist_id, temp);
122
  if (h_err < 0) {
123
    error("Error while reading data array '%s'.", props.name);
124
125
126
127
  }

  /* Unit conversion if necessary */
  const double factor =
128
129
      units_conversion_factor(ic_units, internal_units, props.units);
  if (factor != 1.) {
130

Matthieu Schaller's avatar
Matthieu Schaller committed
131
    /* message("Converting ! factor=%e", factor); */
132

133
    if (io_is_double_precision(props.type)) {
Matthieu Schaller's avatar
Matthieu Schaller committed
134
      double* temp_d = (double*)temp;
135
136
      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= factor;
    } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
137
      float* temp_f = (float*)temp;
138
139
      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor;
    }
140
  }
141

142
143
  /* Clean-up h if necessary */
  const float h_factor_exp = units_h_factor(internal_units, props.units);
144
  if (cleanup_h && h_factor_exp != 0.f && exist != 0) {
145
146
147
148
149
150
151
152
153
154
155
156
157
158
    const double h_factor = pow(h, h_factor_exp);

    /* message("Multipltying '%s' by h^%f=%f", props.name, h_factor_exp,
     * h_factor); */

    if (io_is_double_precision(props.type)) {
      double* temp_d = (double*)temp;
      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= h_factor;
    } else {
      float* temp_f = (float*)temp;
      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= h_factor;
    }
  }

159
  /* Copy temporary buffer to particle data */
Matthieu Schaller's avatar
Matthieu Schaller committed
160
  char* temp_c = (char*)temp;
161
  for (size_t i = 0; i < N; ++i)
162
    memcpy(props.field + i * props.partSize, &temp_c[i * copySize], copySize);
163

164
165
  /* Free and close everything */
  free(temp);
166
167
  H5Sclose(h_filespace);
  H5Sclose(h_memspace);
168
169
170
171
172
173
}

/**
 * @brief Reads a data array from a given HDF5 group.
 *
 * @param grp The group from which to read.
Matthieu Schaller's avatar
Matthieu Schaller committed
174
 * @param props The #io_props of the field to read.
175
176
177
178
179
180
 * @param N The number of particles on that rank.
 * @param N_total The total number of particles.
 * @param mpi_rank The MPI rank of this node.
 * @param offset The offset in the array on disk for this rank.
 * @param internal_units The #unit_system used internally.
 * @param ic_units The #unit_system used in the ICs.
181
182
 * @param cleanup_h Are we removing h-factors from the ICs?
 * @param h The value of the reduced Hubble constant to use for cleaning.
183
 */
184
185
void readArray(hid_t grp, struct io_props props, size_t N, long long N_total,
               int mpi_rank, long long offset,
186
               const struct unit_system* internal_units,
187
               const struct unit_system* ic_units, int cleanup_h, double h) {
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226

  const size_t typeSize = io_sizeof_type(props.type);
  const size_t copySize = typeSize * props.dimension;

  /* Check whether the dataspace exists or not */
  const htri_t exist = H5Lexists(grp, props.name, 0);
  if (exist < 0) {
    error("Error while checking the existence of data set '%s'.", props.name);
  } else if (exist == 0) {
    if (props.importance == COMPULSORY) {
      error("Compulsory data set '%s' not present in the file.", props.name);
    } else {
      for (size_t i = 0; i < N; ++i)
        memset(props.field + i * props.partSize, 0, copySize);
      return;
    }
  }

  /* Open data space in file */
  const hid_t h_data = H5Dopen2(grp, props.name, H5P_DEFAULT);
  if (h_data < 0) error("Error while opening data space '%s'.", props.name);

  /* Check data type */
  const hid_t h_type = H5Dget_type(h_data);
  if (h_type < 0) error("Unable to retrieve data type from the file");
  /* if (!H5Tequal(h_type, hdf5_type(type))) */
  /*   error("Non-matching types between the code and the file"); */

  /* Create property list for collective dataset read. */
  const hid_t h_plist_id = H5Pcreate(H5P_DATASET_XFER);
  H5Pset_dxpl_mpio(h_plist_id, H5FD_MPIO_COLLECTIVE);

  /* Given the limitations of ROM-IO we will need to read the data in chunk of
     HDF5_PARALLEL_IO_MAX_BYTES bytes per node until all the nodes are done. */
  char redo = 1;
  while (redo) {

    /* Maximal number of elements */
    const size_t max_chunk_size =
227
        HDF5_PARALLEL_IO_MAX_BYTES / (props.dimension * typeSize);
228
229
230
231

    /* Write the first chunk */
    const size_t this_chunk = (N > max_chunk_size) ? max_chunk_size : N;
    readArray_chunk(h_data, h_plist_id, props, this_chunk, offset,
232
                    internal_units, ic_units);
233
234
235
236

    /* Compute how many items are left */
    if (N > max_chunk_size) {
      N -= max_chunk_size;
237
238
      props.field += max_chunk_size * props.partSize; /* char* on the field */
      props.parts += max_chunk_size;                  /* part* on the part */
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
      offset += max_chunk_size;
      redo = 1;
    } else {
      N = 0;
      offset += 0;
      redo = 0;
    }

    /* Do we need to run again ? */
    MPI_Allreduce(MPI_IN_PLACE, &redo, 1, MPI_SIGNED_CHAR, MPI_MAX,
                  MPI_COMM_WORLD);

    if (redo && mpi_rank == 0)
      message("Need to redo one iteration for array '%s'", props.name);
  }

  /* Close everything */
  H5Pclose(h_plist_id);
257
258
259
260
261
262
263
264
265
  H5Tclose(h_type);
  H5Dclose(h_data);
}

/*-----------------------------------------------------------------------------
 * Routines writing an output file
 *-----------------------------------------------------------------------------*/

/**
266
 * @brief Prepares an array in the snapshot.
267
 *
268
 * @param e The #engine we are writing from.
269
270
271
272
273
274
275
 * @param grp The HDF5 grp to write to.
 * @param fileName The name of the file we are writing to.
 * @param xmfFile The (opened) XMF file we are appending to.
 * @param partTypeGroupName The name of the group we are writing to.
 * @param props The #io_props of the field to write.
 * @param N_total The total number of particles to write in this array.
 * @param snapshot_units The units used for the data in this snapshot.
276
 */
277
void prepareArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
Matthieu Schaller's avatar
Matthieu Schaller committed
278
279
                  char* partTypeGroupName, struct io_props props,
                  long long N_total, const struct unit_system* snapshot_units) {
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303

  /* Create data space */
  const hid_t h_space = H5Screate(H5S_SIMPLE);
  if (h_space < 0)
    error("Error while creating data space for field '%s'.", props.name);

  int rank = 0;
  hsize_t shape[2];
  hsize_t chunk_shape[2];
  if (props.dimension > 1) {
    rank = 2;
    shape[0] = N_total;
    shape[1] = props.dimension;
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = props.dimension;
  } else {
    rank = 1;
    shape[0] = N_total;
    shape[1] = 0;
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = 0;
  }

  /* Make sure the chunks are not larger than the dataset */
Matthieu Schaller's avatar
Matthieu Schaller committed
304
  if ((long long)chunk_shape[0] > N_total) chunk_shape[0] = N_total;
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322

  /* Change shape of data space */
  hid_t h_err = H5Sset_extent_simple(h_space, rank, shape, NULL);
  if (h_err < 0)
    error("Error while changing data space shape for field '%s'.", props.name);

  /* Create property list for collective dataset write.    */
  const hid_t h_plist_id = H5Pcreate(H5P_DATASET_XFER);
  H5Pset_dxpl_mpio(h_plist_id, H5FD_MPIO_COLLECTIVE);

  /* Set chunk size */
  /* h_err = H5Pset_chunk(h_prop, rank, chunk_shape); */
  /* if (h_err < 0) { */
  /*   error("Error while setting chunk size (%llu, %llu) for field '%s'.", */
  /*         chunk_shape[0], chunk_shape[1], props.name); */
  /* } */

  /* Create dataset */
Matthieu Schaller's avatar
Matthieu Schaller committed
323
324
325
326
  const hid_t h_data =
      H5Dcreate(grp, props.name, io_hdf5_type(props.type), h_space, H5P_DEFAULT,
                H5P_DEFAULT, H5P_DEFAULT);
  if (h_data < 0) error("Error while creating dataspace '%s'.", props.name);
327
328
329
330
331
332
333

  /* Write unit conversion factors for this data set */
  char buffer[FIELD_BUFFER_SIZE];
  units_cgs_conversion_string(buffer, snapshot_units, props.units);
  io_write_attribute_d(
      h_data, "CGS conversion factor",
      units_cgs_conversion_factor(snapshot_units, props.units));
334
  io_write_attribute_f(h_data, "h-scale exponent", 0);
335
336
337
338
339
340
  io_write_attribute_f(h_data, "a-scale exponent",
                       units_a_factor(snapshot_units, props.units));
  io_write_attribute_s(h_data, "Conversion factor", buffer);

  /* Add a line to the XMF */
  xmf_write_line(xmfFile, fileName, partTypeGroupName, props.name, N_total,
Matthieu Schaller's avatar
Matthieu Schaller committed
341
                 props.dimension, props.type);
342
343
344
345
346
347
348

  /* Close everything */
  H5Pclose(h_plist_id);
  H5Dclose(h_data);
  H5Sclose(h_space);
}

349
350
351
352
353
354
355
356
357
358
359
/**
 * @brief Writes a chunk of data in an open HDF5 dataset
 *
 * @param e The #engine we are writing from.
 * @param h_data The HDF5 dataset to write to.
 * @param props The #io_props of the field to write.
 * @param N The number of particles to write.
 * @param offset Offset in the array where this mpi task starts writing.
 * @param internal_units The #unit_system used internally.
 * @param snapshot_units The #unit_system used in the snapshots.
 */
360
void writeArray_chunk(struct engine* e, hid_t h_data,
361
362
363
364
                      const struct io_props props, size_t N, long long offset,
                      const struct unit_system* internal_units,
                      const struct unit_system* snapshot_units) {

365
  const size_t typeSize = io_sizeof_type(props.type);
366
367
  const size_t num_elements = N * props.dimension;

368
  /* Can't handle writes of more than 2GB */
369
  if (N * props.dimension * typeSize > HDF5_PARALLEL_IO_MAX_BYTES)
370
    error("Dataset too large to be written in one pass!");
371

Matthieu Schaller's avatar
Matthieu Schaller committed
372
  /* message("Writing '%s' array...", props.name); */
373
374

  /* Allocate temporary buffer */
375
  void* temp = NULL;
376
  if (posix_memalign((void**)&temp, IO_BUFFER_ALIGNMENT,
377
378
                     num_elements * typeSize) != 0)
    error("Unable to allocate temporary i/o buffer");
379

380
381
382
383
#ifdef IO_SPEED_MEASUREMENT
  MPI_Barrier(MPI_COMM_WORLD);
  ticks tic = getticks();
#endif
384

385
386
  /* Copy the particle data to the temporary buffer */
  io_copy_temp_buffer(temp, e, props, N, internal_units, snapshot_units);
387

388
389
#ifdef IO_SPEED_MEASUREMENT
  MPI_Barrier(MPI_COMM_WORLD);
390
391
392
  if (engine_rank == 0)
    message("Copying for '%s' took %.3f %s.", props.name,
            clocks_from_ticks(getticks() - tic), clocks_getunit());
393
#endif
394

395
  /* Create data space */
396
  const hid_t h_memspace = H5Screate(H5S_SIMPLE);
397
  if (h_memspace < 0) {
398
399
    error("Error while creating data space (memory) for field '%s'.",
          props.name);
400
  }
401

402
403
404
405
  int rank;
  hsize_t shape[2];
  hsize_t offsets[2];
  if (props.dimension > 1) {
406
407
    rank = 2;
    shape[0] = N;
408
    shape[1] = props.dimension;
409
410
411
412
413
414
415
416
417
418
    offsets[0] = offset;
    offsets[1] = 0;
  } else {
    rank = 1;
    shape[0] = N;
    shape[1] = 0;
    offsets[0] = offset;
    offsets[1] = 0;
  }

419
  /* Change shape of memory data space */
420
  hid_t h_err = H5Sset_extent_simple(h_memspace, rank, shape, NULL);
421
422
  if (h_err < 0) {
    error("Error while changing data space (memory) shape for field '%s'.",
423
          props.name);
424
  }
425

426
427
  /* Select the hyper-salb corresponding to this rank */
  hid_t h_filespace = H5Dget_space(h_data);
Matthieu Schaller's avatar
Matthieu Schaller committed
428
429
430
  if (N > 0) {
    H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape,
                        NULL);
431
432
433
  } else {
    H5Sselect_none(h_filespace);
  }
434

435
436
437
/* message("Writing %lld '%s', %zd elements = %zd bytes (int=%d) at offset
 * %zd", N, props.name, N * props.dimension, N * props.dimension * typeSize, */
/* 	  (int)(N * props.dimension * typeSize), offset); */
438

439
440
441
442
#ifdef IO_SPEED_MEASUREMENT
  MPI_Barrier(MPI_COMM_WORLD);
  tic = getticks();
#endif
443

444
445
  /* Write temporary buffer to HDF5 dataspace */
  h_err = H5Dwrite(h_data, io_hdf5_type(props.type), h_memspace, h_filespace,
446
                   H5P_DEFAULT, temp);
447
448
449
450
  if (h_err < 0) {
    error("Error while writing data array '%s'.", props.name);
  }

451
452
453
454
455
456
457
458
459
460
461
#ifdef IO_SPEED_MEASUREMENT
  MPI_Barrier(MPI_COMM_WORLD);
  ticks toc = getticks();
  float ms = clocks_from_ticks(toc - tic);
  int megaBytes = N * props.dimension * typeSize / (1024 * 1024);
  int total = 0;
  MPI_Reduce(&megaBytes, &total, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
  if (engine_rank == 0)
    message("H5Dwrite for '%s' (%d MB) took %.3f %s (speed = %f MB/s).",
            props.name, total, ms, clocks_getunit(), total / (ms / 1000.));
#endif
462

463
464
465
466
467
  /* Free and close everything */
  free(temp);
  H5Sclose(h_memspace);
  H5Sclose(h_filespace);
}
468

469
470
471
472
473
/**
 * @brief Writes a data array in given HDF5 group.
 *
 * @param e The #engine we are writing from.
 * @param grp The group in which to write.
474
 * @param fileName The name of the file in which the data is written.
475
476
477
478
 * @param partTypeGroupName The name of the group containing the particles in
 * the HDF5 file.
 * @param props The #io_props of the field to read
 * @param N The number of particles to write.
479
480
481
482
483
 * @param N_total Total number of particles across all cores.
 * @param mpi_rank The rank of this node.
 * @param offset Offset in the array where this mpi task starts writing.
 * @param internal_units The #unit_system used internally.
 * @param snapshot_units The #unit_system used in the snapshots.
484
 */
485
void writeArray(struct engine* e, hid_t grp, char* fileName,
486
                char* partTypeGroupName, struct io_props props, size_t N,
487
488
489
490
491
                long long N_total, int mpi_rank, long long offset,
                const struct unit_system* internal_units,
                const struct unit_system* snapshot_units) {

  const size_t typeSize = io_sizeof_type(props.type);
492
493
494
495

#ifdef IO_SPEED_MEASUREMENT
  const ticks tic = getticks();
#endif
496

497
498
  /* Open dataset */
  const hid_t h_data = H5Dopen(grp, props.name, H5P_DEFAULT);
Matthieu Schaller's avatar
Matthieu Schaller committed
499
  if (h_data < 0) error("Error while opening dataset '%s'.", props.name);
500

501
502
503
  /* Given the limitations of ROM-IO we will need to write the data in chunk of
     HDF5_PARALLEL_IO_MAX_BYTES bytes per node until all the nodes are done. */
  char redo = 1;
504
505
506
507
508
509
  while (redo) {

    /* Maximal number of elements */
    const size_t max_chunk_size =
        HDF5_PARALLEL_IO_MAX_BYTES / (props.dimension * typeSize);

510
    /* Write the first chunk */
511
    const size_t this_chunk = (N > max_chunk_size) ? max_chunk_size : N;
512
513
    writeArray_chunk(e, h_data, props, this_chunk, offset, internal_units,
                     snapshot_units);
514
515

    /* Compute how many items are left */
516
    if (N > max_chunk_size) {
517
      N -= max_chunk_size;
518
519
      props.field += max_chunk_size * props.partSize; /* char* on the field */
      props.parts += max_chunk_size;                  /* part* on the part */
520
521
      offset += max_chunk_size;
      redo = 1;
522
    } else {
523
524
525
526
527
528
      N = 0;
      offset += 0;
      redo = 0;
    }

    /* Do we need to run again ? */
529
530
    MPI_Allreduce(MPI_IN_PLACE, &redo, 1, MPI_SIGNED_CHAR, MPI_MAX,
                  MPI_COMM_WORLD);
531

532
    if (redo && e->verbose && mpi_rank == 0)
533
      message("Need to redo one iteration for array '%s'", props.name);
534
  }
535

536
  /* Close everything */
537
  H5Dclose(h_data);
538

539
540
#ifdef IO_SPEED_MEASUREMENT
  MPI_Barrier(MPI_COMM_WORLD);
541
542
543
  if (engine_rank == 0)
    message("'%s' took %.3f %s.", props.name,
            clocks_from_ticks(getticks() - tic), clocks_getunit());
544
#endif
545
546
}

547
548
549
550
/**
 * @brief Reads an HDF5 initial condition file (GADGET-3 type) in parallel
 *
 * @param fileName The file to read.
551
 * @param internal_units The system units used internally
552
553
 * @param dim (output) The dimension of the volume read from the file.
 * @param parts (output) The array of #part read from the file.
554
555
556
557
558
 * @param gparts (output) The array of #gpart read from the file.
 * @param sparts (output) The array of #spart read from the file.
 * @param Ngas (output) The number of particles read from the file.
 * @param Ngparts (output) The number of particles read from the file.
 * @param Nstars (output) The number of particles read from the file.
559
 * @param periodic (output) 1 if the volume is periodic, 0 if not.
560
 * @param flag_entropy (output) 1 if the ICs contained Entropy in the
561
 * InternalEnergy field
562
563
564
 * @param with_hydro Are we running with hydro ?
 * @param with_gravity Are we running with gravity ?
 * @param with_stars Are we running with stars ?
565
566
 * @param cleanup_h Are we cleaning-up h-factors from the quantities we read?
 * @param h The value of the reduced Hubble constant to use for correction.
567
568
569
570
 * @param mpi_rank The MPI rank of this node
 * @param mpi_size The number of MPI ranks
 * @param comm The MPI communicator
 * @param info The MPI information object
571
 * @param n_threads The number of threads to use for local operations.
572
 * @param dry_run If 1, don't read the particle. Only allocates the arrays.
573
574
 *
 */
575
void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
576
                      double dim[3], struct part** parts, struct gpart** gparts,
Matthieu Schaller's avatar
Matthieu Schaller committed
577
578
579
                      struct spart** sparts, size_t* Ngas, size_t* Ngparts,
                      size_t* Nstars, int* periodic, int* flag_entropy,
                      int with_hydro, int with_gravity, int with_stars,
580
581
582
                      int cleanup_h, double h, int mpi_rank, int mpi_size,
                      MPI_Comm comm, MPI_Info info, int n_threads,
                      int dry_run) {
583

584
  hid_t h_file = 0, h_grp = 0;
585
  /* GADGET has only cubic boxes (in cosmological mode) */
Matthieu Schaller's avatar
Matthieu Schaller committed
586
  double boxSize[3] = {0.0, -1.0, -1.0};
587
588
589
590
591
  long long numParticles[swift_type_count] = {0};
  long long numParticles_highWord[swift_type_count] = {0};
  size_t N[swift_type_count] = {0};
  long long N_total[swift_type_count] = {0};
  long long offset[swift_type_count] = {0};
592
  int dimension = 3; /* Assume 3D if nothing is specified */
593
  size_t Ndm = 0;
594
595
596
597
598
599
600
601
602
603
604
605

  /* Open file */
  /* message("Opening file '%s' as IC.", fileName); */
  hid_t h_plist_id = H5Pcreate(H5P_FILE_ACCESS);
  H5Pset_fapl_mpio(h_plist_id, comm, info);
  h_file = H5Fopen(fileName, H5F_ACC_RDONLY, h_plist_id);
  if (h_file < 0) {
    error("Error while opening file '%s'.", fileName);
  }

  /* Open header to read simulation properties */
  /* message("Reading runtime parameters..."); */
606
  h_grp = H5Gopen(h_file, "/RuntimePars", H5P_DEFAULT);
607
608
609
  if (h_grp < 0) error("Error while opening runtime parameters\n");

  /* Read the relevant information */
610
  io_read_attribute(h_grp, "PeriodicBoundariesOn", INT, periodic);
611
612
613
614
615
616

  /* Close runtime parameters */
  H5Gclose(h_grp);

  /* Open header to read simulation properties */
  /* message("Reading file header..."); */
617
  h_grp = H5Gopen(h_file, "/Header", H5P_DEFAULT);
618
619
  if (h_grp < 0) error("Error while opening file header\n");

620
621
622
623
  /* Check the dimensionality of the ICs (if the info exists) */
  const hid_t hid_dim = H5Aexists(h_grp, "Dimension");
  if (hid_dim < 0)
    error("Error while testing existance of 'Dimension' attribute");
624
  if (hid_dim > 0) io_read_attribute(h_grp, "Dimension", INT, &dimension);
625
626
627
628
  if (dimension != hydro_dimension)
    error("ICs dimensionality (%dD) does not match code dimensionality (%dD)",
          dimension, (int)hydro_dimension);

629
  /* Read the relevant information and print status */
630
  int flag_entropy_temp[6];
631
  io_read_attribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy_temp);
632
  *flag_entropy = flag_entropy_temp[0];
633
634
635
636
  io_read_attribute(h_grp, "BoxSize", DOUBLE, boxSize);
  io_read_attribute(h_grp, "NumPart_Total", LONGLONG, numParticles);
  io_read_attribute(h_grp, "NumPart_Total_HighWord", LONGLONG,
                    numParticles_highWord);
637

638
  for (int ptype = 0; ptype < swift_type_count; ++ptype)
Matthieu Schaller's avatar
Matthieu Schaller committed
639
640
    N_total[ptype] =
        (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);
641

642
  /* Get the box size if not cubic */
643
644
645
646
  dim[0] = boxSize[0];
  dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
  dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];

647
648
649
650
651
652
  /* Change box size in the 1D and 2D case */
  if (hydro_dimension == 2)
    dim[2] = min(dim[0], dim[1]);
  else if (hydro_dimension == 1)
    dim[2] = dim[1] = dim[0];

653
654
655
656
657
658
659
  /* Convert the box size if we want to clean-up h-factors */
  if (cleanup_h) {
    dim[0] /= h;
    dim[1] /= h;
    dim[2] /= h;
  }

660
661
662
  /* message("Found %lld particles in a %speriodic box of size [%f %f %f].", */
  /* 	  N_total[0], (periodic ? "": "non-"), dim[0], */
  /* 	  dim[1], dim[2]); */
663
664

  /* Divide the particles among the tasks. */
665
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
666
667
668
    offset[ptype] = mpi_rank * N_total[ptype] / mpi_size;
    N[ptype] = (mpi_rank + 1) * N_total[ptype] / mpi_size - offset[ptype];
  }
669
670
671
672

  /* Close header */
  H5Gclose(h_grp);

673
  /* Read the unit system used in the ICs */
Matthieu Schaller's avatar
Matthieu Schaller committed
674
675
  struct unit_system* ic_units =
      (struct unit_system*)malloc(sizeof(struct unit_system));
676
  if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
677
  io_read_unit_system(h_file, ic_units, mpi_rank);
678
679
680
681
682

  /* Tell the user if a conversion will be needed */
  if (mpi_rank == 0) {
    if (units_are_equal(ic_units, internal_units)) {

Matthieu Schaller's avatar
Matthieu Schaller committed
683
      message("IC and internal units match. No conversion needed.");
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709

    } else {

      message("Conversion needed from:");
      message("(ICs) Unit system: U_M =      %e g.", ic_units->UnitMass_in_cgs);
      message("(ICs) Unit system: U_L =      %e cm.",
              ic_units->UnitLength_in_cgs);
      message("(ICs) Unit system: U_t =      %e s.", ic_units->UnitTime_in_cgs);
      message("(ICs) Unit system: U_I =      %e A.",
              ic_units->UnitCurrent_in_cgs);
      message("(ICs) Unit system: U_T =      %e K.",
              ic_units->UnitTemperature_in_cgs);
      message("to:");
      message("(internal) Unit system: U_M = %e g.",
              internal_units->UnitMass_in_cgs);
      message("(internal) Unit system: U_L = %e cm.",
              internal_units->UnitLength_in_cgs);
      message("(internal) Unit system: U_t = %e s.",
              internal_units->UnitTime_in_cgs);
      message("(internal) Unit system: U_I = %e A.",
              internal_units->UnitCurrent_in_cgs);
      message("(internal) Unit system: U_T = %e K.",
              internal_units->UnitTemperature_in_cgs);
    }
  }

710
  /* Convert the dimensions of the box */
Matthieu Schaller's avatar
Matthieu Schaller committed
711
712
713
  for (int j = 0; j < 3; j++)
    dim[j] *=
        units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
714

715
  /* Allocate memory to store SPH particles */
716
717
  if (with_hydro) {
    *Ngas = N[0];
718
    if (posix_memalign((void**)parts, part_align,
Matthieu Schaller's avatar
Matthieu Schaller committed
719
                       (*Ngas) * sizeof(struct part)) != 0)
720
721
722
723
724
725
      error("Error while allocating memory for particles");
    bzero(*parts, *Ngas * sizeof(struct part));
  }

  /* Allocate memory to store star particles */
  if (with_stars) {
726
    *Nstars = N[swift_type_star];
727
    if (posix_memalign((void**)sparts, spart_align,
728
729
730
731
732
733
734
735
                       *Nstars * sizeof(struct spart)) != 0)
      error("Error while allocating memory for star particles");
    bzero(*sparts, *Nstars * sizeof(struct spart));
  }

  /* Allocate memory to store gravity particles */
  if (with_gravity) {
    Ndm = N[1];
736
737
738
    *Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
               N[swift_type_dark_matter] +
               (with_stars ? N[swift_type_star] : 0);
739
    if (posix_memalign((void**)gparts, gpart_align,
Matthieu Schaller's avatar
Matthieu Schaller committed
740
                       *Ngparts * sizeof(struct gpart)) != 0)
741
742
743
744
745
      error("Error while allocating memory for gravity particles");
    bzero(*gparts, *Ngparts * sizeof(struct gpart));
  }

  /* message("Allocated %8.2f MB for particles.", *N * sizeof(struct part) /
746
747
   * (1024.*1024.)); */

748
  /* message("BoxSize = %lf", dim[0]); */
749
  /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */
750
751

  /* Loop over all particle types */
752
  for (int ptype = 0; ptype < swift_type_count; ptype++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
753

754
755
    /* Don't do anything if no particle of this kind */
    if (N_total[ptype] == 0) continue;
Matthieu Schaller's avatar
Matthieu Schaller committed
756

757
758
759
    /* Open the particle group in the file */
    char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE];
    snprintf(partTypeGroupName, PARTICLE_GROUP_BUFFER_SIZE, "/PartType%d",
Matthieu Schaller's avatar
Matthieu Schaller committed
760
             ptype);
761
762
763
764
    h_grp = H5Gopen(h_file, partTypeGroupName, H5P_DEFAULT);
    if (h_grp < 0) {
      error("Error while opening particle group %s.", partTypeGroupName);
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
765

766
767
    int num_fields = 0;
    struct io_props list[100];
768
    size_t Nparticles = 0;
769

770
771
    /* Read particle fields into the particle structure */
    switch (ptype) {
Matthieu Schaller's avatar
Matthieu Schaller committed
772

773
      case swift_type_gas:
Matthieu Schaller's avatar
Matthieu Schaller committed
774
775
776
        if (with_hydro) {
          Nparticles = *Ngas;
          hydro_read_particles(*parts, list, &num_fields);
777
          num_fields += chemistry_read_particles(*parts, list + num_fields);
Matthieu Schaller's avatar
Matthieu Schaller committed
778
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
779
780
        break;

781
      case swift_type_dark_matter:
Matthieu Schaller's avatar
Matthieu Schaller committed
782
783
784
785
        if (with_gravity) {
          Nparticles = Ndm;
          darkmatter_read_particles(*gparts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
786
787
        break;

788
      case swift_type_star:
789
        if (with_stars) {
Matthieu Schaller's avatar
Matthieu Schaller committed
790
791
792
793
          Nparticles = *Nstars;
          star_read_particles(*sparts, list, &num_fields);
        }
        break;
794

Matthieu Schaller's avatar
Matthieu Schaller committed
795
      default:
796
797
798
        if (mpi_rank == 0)
          message("Particle Type %d not yet supported. Particles ignored",
                  ptype);
799
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
800

801
802
803
    /* Read everything */
    if (!dry_run)
      for (int i = 0; i < num_fields; ++i)
804
        readArray(h_grp, list[i], Nparticles, N_total[ptype], mpi_rank,
805
                  offset[ptype], internal_units, ic_units, cleanup_h, h);
806

807
808
809
    /* Close particle group */
    H5Gclose(h_grp);
  }
810

811
  if (!dry_run && with_gravity) {
812

813
814
    /* Let's initialise a bit of thread parallelism here */
    struct threadpool tp;
815
    threadpool_init(&tp, n_threads);
816

817
    /* Prepare the DM particles */
818
    io_prepare_dm_gparts(&tp, *gparts, Ndm);
819
820

    /* Duplicate the hydro particles into gparts */
821
    if (with_hydro) io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas, Ndm);
822
823
824

    /* Duplicate the star particles into gparts */
    if (with_stars)
825
      io_duplicate_star_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
826
827
828

    threadpool_clean(&tp);
  }
829
830

  /* message("Done Reading particles..."); */
831

832
833
834
  /* Clean up */
  free(ic_units);

835
836
837
838
839
840
841
  /* Close property handler */
  H5Pclose(h_plist_id);

  /* Close file */
  H5Fclose(h_file);
}

842
843
844
845
846
847
848
849
850
/**
 * @brief Prepares a file for a parallel write.
 *
 * @param e The #engine.
 * @param baseName The base name of the snapshots.
 * @param N_total The total number of particles of each type to write.
 * @param internal_units The #unit_system used internally.
 * @param snapshot_units The #unit_system used in the snapshots.
 */
851
void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
Matthieu Schaller's avatar
Matthieu Schaller committed
852
853
                  const struct unit_system* internal_units,
                  const struct unit_system* snapshot_units) {
854

855
856
857
858
  const struct part* parts = e->s->parts;
  const struct xpart* xparts = e->s->xparts;
  const struct gpart* gparts = e->s->gparts;
  const struct spart* sparts = e->s->sparts;
859
  FILE* xmfFile = 0;
860
861
  int periodic = e->s->periodic;
  int numFiles = 1;
862
863

  /* First time, we need to create the XMF file */
864
  if (e->snapshotOutputCount == 0) xmf_create_file(baseName);
865

866
  /* Prepare the XMF file for the new entry */
867
  xmfFile = xmf_prepare_file(baseName);
868

869
870
871
  /* HDF5 File name */
  char fileName[FILENAME_BUFFER_SIZE];
  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
872
           e->snapshotOutputCount);
873

874
  /* Open HDF5 file with the chosen parameters */
875
  hid_t h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
Matthieu Schaller's avatar
Matthieu Schaller committed
876
  if (h_file < 0) error("Error while opening file '%s'.", fileName);
877

Matthieu Schaller's avatar
Matthieu Schaller committed
878
879
  /* Write the part of the XMF file corresponding to this
   * specific output */
880
  xmf_write_outputheader(xmfFile, fileName, e->time);
881

882
883
  /* Open header to write simulation properties */
  /* message("Writing runtime parameters..."); */
884
  hid_t h_grp =
Matthieu Schaller's avatar
Matthieu Schaller committed
885
      H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
886
  if (h_grp < 0) error("Error while creating runtime parameters group\n");
887
888

  /* Write the relevant information */
889
  io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
890
891
892

  /* Close runtime parameters */
  H5Gclose(h_grp);
893

894
895
  /* Open header to write simulation properties */
  /* message("Writing file header..."); */
896
  h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
897
  if (h_grp < 0) error("Error while creating file header\n");
Matthieu Schaller's avatar
Matthieu Schaller committed
898

899
  /* Print the relevant information and print status */
900
  io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
901
  double dblTime = e->time;
902
  io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
903
  int dimension = (int)hydro_dimension;
904
  io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
905
906
  io_write_attribute(h_grp, "Redshift", DOUBLE, &e->cosmology->z, 1);
  io_write_attribute(h_grp, "Scale-factor", DOUBLE, &e->cosmology->a, 1);
Matthieu Schaller's avatar
Matthieu Schaller committed
907

908
  /* GADGET-2 legacy values */
909
  /* Number of particles of each type */
910
911
912
  unsigned int numParticles[swift_type_count] = {0};
  unsigned int numParticlesHighWord[swift_type_count] = {0};
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
913
914
915
    numParticles[ptype] = (unsigned int)N_total[ptype];
    numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
  }
916
917
918
919
920
921
  io_write_attribute(h_grp, "NumPart_ThisFile", LONGLONG, N_total,
                     swift_type_count);
  io_write_attribute(h_grp, "NumPart_Total", UINT, numParticles,
                     swift_type_count);
  io_write_attribute(h_grp, "NumPart_Total_HighWord", UINT,
                     numParticlesHighWord, swift_type_count);
922
  double MassTable[6] = {0., 0., 0., 0., 0., 0.};
923
924
  io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
  unsigned int flagEntropy[swift_type_count] = {0};
925
  flagEntropy[0] = writeEntropyFlag();
926
927
928
  io_write_attribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
                     swift_type_count);
  io_write_attribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
Matthieu Schaller's avatar
Matthieu Schaller committed
929

930
931
  /* Close header */
  H5Gclose(h_grp);
Matthieu Schaller's avatar
Matthieu Schaller committed
932

933
  /* Print the code version */
934
  io_write_code_description(h_file);
935

936
937
938
  /* Print the run's policy */
  io_write_engine_policy(h_file, e);

939
  /* Print the SPH parameters */
940
941
942
  if (e->policy & engine_policy_hydro) {
    h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
943
944
    if (h_grp < 0) error("Error while creating SPH group");
    hydro_props_print_snapshot(h_grp, e->hydro_properties);
945
    hydro_write_flavour(h_grp);
946
947
    H5Gclose(h_grp);
  }
948

949
950
951
952
  /* Print the subgrid parameters */
  h_grp = H5Gcreate(h_file, "/SubgridScheme", H5P_DEFAULT, H5P_DEFAULT,
                    H5P_DEFAULT);
  if (h_grp < 0) error("Error while creating subgrid group");
953
954
  cooling_write_flavour(h_grp);
  chemistry_write_flavour(h_grp);
955
956
  H5Gclose(h_grp);

957
  /* Print the gravity parameters */
958
  if (e->policy & engine_policy_self_gravity) {
959
    h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
960
                      H5P_DEFAULT);
961
962
963
964
    if (h_grp < 0) error("Error while creating gravity group");
    gravity_props_print_snapshot(h_grp, e->gravity_properties);
    H5Gclose(h_grp);
  }
965

966
967
  /* Print the gravity parameters */
  if (e->policy & engine_policy_cosmology) {
Matthieu Schaller's avatar