parallel_io.c 38 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 "common_io.h"
40
#include "dimension.h"
41
#include "engine.h"
42
#include "error.h"
43
#include "gravity_io.h"
44
#include "gravity_properties.h"
45
#include "hydro_io.h"
46
#include "hydro_properties.h"
47
#include "io_properties.h"
48
49
#include "kernel_hydro.h"
#include "part.h"
50
#include "stars_io.h"
51
#include "units.h"
52
#include "xmf.h"
53

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

57
/**
58
 * @brief Reads a chunk of data from an open HDF5 dataset
59
 *
60
61
62
63
64
 * @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.
65
 * @param internal_units The #unit_system used internally.
66
 * @param ic_units The #unit_system used in the snapshots.
67
 */
68
void readArray_chunk(hid_t h_data, hid_t h_plist_id,
69
70
71
                     const struct io_props props, size_t N, long long offset,
                     const struct unit_system* internal_units,
                     const struct unit_system* ic_units) {
72

73
74
75
  const size_t typeSize = io_sizeof_type(props.type);
  const size_t copySize = typeSize * props.dimension;
  const size_t num_elements = N * props.dimension;
76

77
78
79
  /* 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!");
80

81
  /* Allocate temporary buffer */
82
  void* temp = malloc(num_elements * typeSize);
83
  if (temp == NULL) error("Unable to allocate memory for temporary buffer");
84

85
  /* Prepare information for hyper-slab */
86
87
  hsize_t shape[2], offsets[2];
  int rank;
88
  if (props.dimension > 1) {
89
90
    rank = 2;
    shape[0] = N;
91
    shape[1] = props.dimension;
92
93
94
    offsets[0] = offset;
    offsets[1] = 0;
  } else {
95
    rank = 2;
96
    shape[0] = N;
97
    shape[1] = 1;
98
99
100
    offsets[0] = offset;
    offsets[1] = 0;
  }
101
102

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

105
  /* Select hyper-slab in file */
106
  const hid_t h_filespace = H5Dget_space(h_data);
107
108
  H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape, NULL);

109
110
111
  /* Read HDF5 dataspace in temporary buffer */
  /* Dirty version that happens to work for vectors but should be improved */
  /* Using HDF5 dataspaces would be better */
112
  const hid_t h_err = H5Dread(h_data, io_hdf5_type(props.type), h_memspace,
113
                              h_filespace, h_plist_id, temp);
114
  if (h_err < 0) {
115
    error("Error while reading data array '%s'.", props.name);
116
117
118
119
  }

  /* Unit conversion if necessary */
  const double factor =
120
121
      units_conversion_factor(ic_units, internal_units, props.units);
  if (factor != 1.) {
122

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

125
    if (io_is_double_precision(props.type)) {
126
127
128
129
130
131
      double* temp_d = temp;
      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= factor;
    } else {
      float* temp_f = temp;
      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor;
    }
132
  }
133
134

  /* Copy temporary buffer to particle data */
135
136
  char* temp_c = temp;
  for (size_t i = 0; i < N; ++i)
137
    memcpy(props.field + i * props.partSize, &temp_c[i * copySize], copySize);
138

139
140
  /* Free and close everything */
  free(temp);
141
142
  H5Sclose(h_filespace);
  H5Sclose(h_memspace);
143
144
145
146
147
148
}

/**
 * @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
149
 * @param props The #io_props of the field to read.
150
151
152
153
154
155
156
 * @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.
 */
157
158
void readArray(hid_t grp, struct io_props props, size_t N, long long N_total,
               int mpi_rank, long long offset,
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
               const struct unit_system* internal_units,
               const struct unit_system* ic_units) {

  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 =
200
        HDF5_PARALLEL_IO_MAX_BYTES / (props.dimension * typeSize);
201
202
203
204

    /* 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,
205
                    internal_units, ic_units);
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228

    /* Compute how many items are left */
    if (N > max_chunk_size) {
      N -= max_chunk_size;
      props.field += max_chunk_size * props.partSize;
      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);
229
230
231
232
233
234
235
236
237
  H5Tclose(h_type);
  H5Dclose(h_data);
}

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

/**
238
 * @brief Writes a chunk of data in an open HDF5 dataset
239
 *
240
 * @param e The #engine we are writing from.
241
 * @param h_data The HDF5 dataset to write to.
242
 * @param h_plist_id the parallel HDF5 properties.
243
 * @param props The #io_props of the field to read.
244
 * @param N The number of particles to write.
245
246
247
 * @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.
248
 */
249
250
251
252
253
void writeArray_chunk(struct engine* e, hid_t h_data, hid_t h_plist_id,
                      const struct io_props props, size_t N, long long offset,
                      const struct unit_system* internal_units,
                      const struct unit_system* snapshot_units) {

254
  const size_t typeSize = io_sizeof_type(props.type);
255
256
  const size_t num_elements = N * props.dimension;

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

Matthieu Schaller's avatar
Matthieu Schaller committed
261
  /* message("Writing '%s' array...", props.name); */
262
263

  /* Allocate temporary buffer */
264
  void* temp = malloc(num_elements * typeSize);
265
  if (posix_memalign((void**)&temp, IO_BUFFER_ALIGNMENT,
266
267
                     num_elements * typeSize) != 0)
    error("Unable to allocate temporary i/o buffer");
268

269
270
  /* MPI_Barrier(MPI_COMM_WORLD); */
  /* ticks tic = getticks(); */
271

272
273
  /* Copy the particle data to the temporary buffer */
  io_copy_temp_buffer(temp, e, props, N, internal_units, snapshot_units);
274

275
276
277
278
  /* MPI_Barrier(MPI_COMM_WORLD); */
  /* if(engine_rank == 0) */
  /*   message( "Copying for '%s' took %.3f %s." , props.name, */
  /* 	     clocks_from_ticks(getticks() - tic), clocks_getunit()); */
279

280
  /* Create data space */
281
  const hid_t h_memspace = H5Screate(H5S_SIMPLE);
282
  if (h_memspace < 0) {
283
284
    error("Error while creating data space (memory) for field '%s'.",
          props.name);
285
  }
286

287
288
289
290
  int rank;
  hsize_t shape[2];
  hsize_t offsets[2];
  if (props.dimension > 1) {
291
292
    rank = 2;
    shape[0] = N;
293
    shape[1] = props.dimension;
294
295
296
297
298
299
300
301
302
303
    offsets[0] = offset;
    offsets[1] = 0;
  } else {
    rank = 1;
    shape[0] = N;
    shape[1] = 0;
    offsets[0] = offset;
    offsets[1] = 0;
  }

304
  /* Change shape of memory data space */
305
  hid_t h_err = H5Sset_extent_simple(h_memspace, rank, shape, NULL);
306
307
  if (h_err < 0) {
    error("Error while changing data space (memory) shape for field '%s'.",
308
          props.name);
309
  }
310

311
312
  /* Select the hyper-salb corresponding to this rank */
  hid_t h_filespace = H5Dget_space(h_data);
Matthieu Schaller's avatar
Matthieu Schaller committed
313
314
315
  if (N > 0) {
    H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape,
                        NULL);
316
317
318
  } else {
    H5Sselect_none(h_filespace);
  }
319

Matthieu Schaller's avatar
Matthieu Schaller committed
320
321
  /* message("Writing %lld '%s', %zd elements = %zd bytes (int=%d) at offset
   * %zd", */
322
  /* 	  N, props.name, N * props.dimension, N * props.dimension * typeSize, */
323
324
  /* 	  (int)(N * props.dimension * typeSize), offset); */

325
326
  /* MPI_Barrier(MPI_COMM_WORLD); */
  /* tic = getticks(); */
327

328
329
  /* Write temporary buffer to HDF5 dataspace */
  h_err = H5Dwrite(h_data, io_hdf5_type(props.type), h_memspace, h_filespace,
Matthieu Schaller's avatar
Matthieu Schaller committed
330
                   h_plist_id, temp);
331
332
333
334
  if (h_err < 0) {
    error("Error while writing data array '%s'.", props.name);
  }

335
336
337
338
339
340
341
342
343
344
345
  /* 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.));

346
347
348
349
350
  /* Free and close everything */
  free(temp);
  H5Sclose(h_memspace);
  H5Sclose(h_filespace);
}
351

352
353
354
355
356
/**
 * @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.
357
358
 * @param fileName The name of the file in which the data is written.
 * @param xmfFile The FILE used to write the XMF description.
359
360
361
362
 * @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.
363
364
365
366
367
 * @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.
368
 */
369
void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
370
                char* partTypeGroupName, struct io_props props, size_t N,
371
372
373
374
375
                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);
376
  /* const ticks tic = getticks(); */
377
378
379
380

  /* Work out properties of the array in the file */
  int rank;
  hsize_t shape_total[2];
381
  hsize_t chunk_shape[2];
382
383
384
385
  if (props.dimension > 1) {
    rank = 2;
    shape_total[0] = N_total;
    shape_total[1] = props.dimension;
386
387
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = props.dimension;
388
389
390
391
  } else {
    rank = 1;
    shape_total[0] = N_total;
    shape_total[1] = 0;
392
393
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = 0;
394
395
  }

396
  /* Make sure the chunks are not larger than the dataset */
lhausamm's avatar
lhausamm committed
397
  if (chunk_shape[0] > (hsize_t)N_total) chunk_shape[0] = N_total;
398
399

  /* Create the space in the file */
400
401
402
403
  hid_t h_filespace = H5Screate(H5S_SIMPLE);
  if (h_filespace < 0) {
    error("Error while creating data space (file) for field '%s'.", props.name);
  }
404

405
  /* Change shape of file data space */
406
  hid_t h_err = H5Sset_extent_simple(h_filespace, rank, shape_total, NULL);
407
  if (h_err < 0) {
408
409
    error("Error while changing data space (file) shape for field '%s'.",
          props.name);
410
411
  }

412
413
414
415
  /* Dataset properties */
  const hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);

  /* Set chunk size */
416
417
418
419
420
  /* 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); */
  /* } */
421

422
  /* Create dataset */
423
424
  const hid_t h_data = H5Dcreate(grp, props.name, io_hdf5_type(props.type),
                                 h_filespace, H5P_DEFAULT, h_prop, H5P_DEFAULT);
425
  if (h_data < 0) {
426
    error("Error while creating dataset '%s'.", props.name);
427
428
  }

429
430
431
  H5Sclose(h_filespace);

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

435
436
437
  /* 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;
438
439
440
441
442
443
  while (redo) {

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

444
    /* Write the first chunk */
445
    const size_t this_chunk = (N > max_chunk_size) ? max_chunk_size : N;
446
    writeArray_chunk(e, h_data, h_plist_id, props, this_chunk, offset,
447
                     internal_units, snapshot_units);
448
449

    /* Compute how many items are left */
450
    if (N > max_chunk_size) {
451
      N -= max_chunk_size;
452
      props.field += max_chunk_size * props.partSize;
453
454
      offset += max_chunk_size;
      redo = 1;
455
    } else {
456
457
458
459
460
461
      N = 0;
      offset += 0;
      redo = 0;
    }

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

465
    if (redo /* && e->verbose*/ && mpi_rank == 0)
466
      message("Need to redo one iteration for array '%s'", props.name);
467
  }
468
469

  /* Write XMF description for this data set */
Matthieu Schaller's avatar
Matthieu Schaller committed
470
  if (mpi_rank == 0)
471
472
    xmf_write_line(xmfFile, fileName, partTypeGroupName, props.name, N_total,
                   props.dimension, props.type);
473

474
  /* Write unit conversion factors for this data set */
475
476
  char buffer[FIELD_BUFFER_SIZE];
  units_cgs_conversion_string(buffer, snapshot_units, props.units);
477
478
479
480
481
482
483
484
  io_write_attribute_d(
      h_data, "CGS conversion factor",
      units_cgs_conversion_factor(snapshot_units, props.units));
  io_write_attribute_f(h_data, "h-scale exponent",
                       units_h_factor(snapshot_units, props.units));
  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);
485

486
487
  /* Close everything */
  H5Pclose(h_prop);
488
489
  H5Dclose(h_data);
  H5Pclose(h_plist_id);
490

491
492
493
494
  /* MPI_Barrier(MPI_COMM_WORLD); */
  /* if(engine_rank == 0) */
  /*   message( "'%s' took %.3f %s." , props.name, */
  /* 	     clocks_from_ticks(getticks() - tic), clocks_getunit()); */
495
496
}

497
498
499
500
/**
 * @brief Reads an HDF5 initial condition file (GADGET-3 type) in parallel
 *
 * @param fileName The file to read.
501
 * @param internal_units The system units used internally
502
503
 * @param dim (output) The dimension of the volume read from the file.
 * @param parts (output) The array of #part read from the file.
504
505
506
507
508
 * @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.
509
 * @param periodic (output) 1 if the volume is periodic, 0 if not.
510
 * @param flag_entropy (output) 1 if the ICs contained Entropy in the
511
 * InternalEnergy field
512
513
514
 * @param with_hydro Are we running with hydro ?
 * @param with_gravity Are we running with gravity ?
 * @param with_stars Are we running with stars ?
515
516
517
518
 * @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
519
 * @param dry_run If 1, don't read the particle. Only allocates the arrays.
520
521
 *
 */
522
void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
523
                      double dim[3], struct part** parts, struct gpart** gparts,
Matthieu Schaller's avatar
Matthieu Schaller committed
524
525
526
527
528
                      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,
                      int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
                      int dry_run) {
529

530
  hid_t h_file = 0, h_grp = 0;
531
  /* GADGET has only cubic boxes (in cosmological mode) */
Matthieu Schaller's avatar
Matthieu Schaller committed
532
  double boxSize[3] = {0.0, -1.0, -1.0};
533
534
535
536
537
  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};
538
  int dimension = 3; /* Assume 3D if nothing is specified */
539
  size_t Ndm = 0;
540
541
542
543
544
545
546
547
548
549
550
551

  /* 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..."); */
552
  h_grp = H5Gopen(h_file, "/RuntimePars", H5P_DEFAULT);
553
554
555
  if (h_grp < 0) error("Error while opening runtime parameters\n");

  /* Read the relevant information */
556
  io_read_attribute(h_grp, "PeriodicBoundariesOn", INT, periodic);
557
558
559
560
561
562

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

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

566
567
568
569
  /* 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");
570
  if (hid_dim > 0) io_read_attribute(h_grp, "Dimension", INT, &dimension);
571
572
573
574
  if (dimension != hydro_dimension)
    error("ICs dimensionality (%dD) does not match code dimensionality (%dD)",
          dimension, (int)hydro_dimension);

575
  /* Read the relevant information and print status */
576
  int flag_entropy_temp[6];
577
  io_read_attribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy_temp);
578
  *flag_entropy = flag_entropy_temp[0];
579
580
581
582
  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);
583

584
  for (int ptype = 0; ptype < swift_type_count; ++ptype)
Matthieu Schaller's avatar
Matthieu Schaller committed
585
586
    N_total[ptype] =
        (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);
587

588
  /* Get the box size if not cubic */
589
590
591
592
  dim[0] = boxSize[0];
  dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
  dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];

593
594
595
596
597
598
  /* 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];

599
600
601
  /* message("Found %lld particles in a %speriodic box of size [%f %f %f].", */
  /* 	  N_total[0], (periodic ? "": "non-"), dim[0], */
  /* 	  dim[1], dim[2]); */
602
603

  /* Divide the particles among the tasks. */
604
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
605
606
607
    offset[ptype] = mpi_rank * N_total[ptype] / mpi_size;
    N[ptype] = (mpi_rank + 1) * N_total[ptype] / mpi_size - offset[ptype];
  }
608
609
610
611

  /* Close header */
  H5Gclose(h_grp);

612
  /* Read the unit system used in the ICs */
613
  struct unit_system* ic_units = malloc(sizeof(struct unit_system));
614
  if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
615
  io_read_unit_system(h_file, ic_units, mpi_rank);
616
617
618
619
620

  /* 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
621
      message("IC and internal units match. No conversion needed.");
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647

    } 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);
    }
  }

648
  /* Convert the dimensions of the box */
Matthieu Schaller's avatar
Matthieu Schaller committed
649
650
651
  for (int j = 0; j < 3; j++)
    dim[j] *=
        units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
652

653
  /* Allocate memory to store SPH particles */
654
655
  if (with_hydro) {
    *Ngas = N[0];
Matthieu Schaller's avatar
Matthieu Schaller committed
656
657
    if (posix_memalign((void*)parts, part_align,
                       (*Ngas) * sizeof(struct part)) != 0)
658
659
660
661
662
663
      error("Error while allocating memory for particles");
    bzero(*parts, *Ngas * sizeof(struct part));
  }

  /* Allocate memory to store star particles */
  if (with_stars) {
664
    *Nstars = N[swift_type_star];
665
666
667
668
669
670
671
672
673
    if (posix_memalign((void*)sparts, spart_align,
                       *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];
674
675
676
    *Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
               N[swift_type_dark_matter] +
               (with_stars ? N[swift_type_star] : 0);
677
    if (posix_memalign((void*)gparts, gpart_align,
Matthieu Schaller's avatar
Matthieu Schaller committed
678
                       *Ngparts * sizeof(struct gpart)) != 0)
679
680
681
682
683
      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) /
684
685
   * (1024.*1024.)); */

686
  /* message("BoxSize = %lf", dim[0]); */
687
  /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */
688
689

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

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

695
696
697
    /* 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
698
             ptype);
699
700
701
702
    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
703

704
705
    int num_fields = 0;
    struct io_props list[100];
706
    size_t Nparticles = 0;
707

708
709
    /* Read particle fields into the particle structure */
    switch (ptype) {
Matthieu Schaller's avatar
Matthieu Schaller committed
710

711
      case swift_type_gas:
Matthieu Schaller's avatar
Matthieu Schaller committed
712
713
714
715
        if (with_hydro) {
          Nparticles = *Ngas;
          hydro_read_particles(*parts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
716
717
        break;

718
      case swift_type_dark_matter:
Matthieu Schaller's avatar
Matthieu Schaller committed
719
720
721
722
        if (with_gravity) {
          Nparticles = Ndm;
          darkmatter_read_particles(*gparts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
723
724
        break;

725
      case swift_type_star:
726
        if (with_stars) {
Matthieu Schaller's avatar
Matthieu Schaller committed
727
728
729
730
          Nparticles = *Nstars;
          star_read_particles(*sparts, list, &num_fields);
        }
        break;
731

Matthieu Schaller's avatar
Matthieu Schaller committed
732
      default:
733
734
735
        if (mpi_rank == 0)
          message("Particle Type %d not yet supported. Particles ignored",
                  ptype);
736
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
737

738
739
740
    /* Read everything */
    if (!dry_run)
      for (int i = 0; i < num_fields; ++i)
741
742
        readArray(h_grp, list[i], Nparticles, N_total[ptype], mpi_rank,
                  offset[ptype], internal_units, ic_units);
743

744
745
746
    /* Close particle group */
    H5Gclose(h_grp);
  }
747

748
  /* Prepare the DM particles */
749
  if (!dry_run && with_gravity) io_prepare_dm_gparts(*gparts, Ndm);
750

751
752
  /* Duplicate the hydro particles into gparts */
  if (!dry_run && with_gravity && with_hydro)
753
    io_duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
754
755
756

  /* Duplicate the star particles into gparts */
  if (!dry_run && with_gravity && with_stars)
757
    io_duplicate_star_gparts(*sparts, *gparts, *Nstars, Ndm + *Ngas);
758
759

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

761
762
763
  /* Clean up */
  free(ic_units);

764
765
766
767
768
769
770
  /* Close property handler */
  H5Pclose(h_plist_id);

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

771
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
772
773
 * @brief Writes an HDF5 output file (GADGET-3 type) with
 *its XMF descriptor
774
775
 *
 * @param e The engine containing all the system.
776
 * @param baseName The common part of the snapshot file name.
777
778
 * @param internal_units The #unit_system used internally
 * @param snapshot_units The #unit_system used in the snapshots
779
780
781
782
 * @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
783
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
784
 * Creates an HDF5 output file and writes the particles
785
786
 * contained in the engine. If such a file already exists, it is
 * erased and replaced by the new one.
787
788
789
790
791
 * The companion XMF file is also updated accordingly.
 *
 * Calls #error() if an error occurs.
 *
 */
792
void write_output_parallel(struct engine* e, const char* baseName,
793
794
                           const struct unit_system* internal_units,
                           const struct unit_system* snapshot_units,
795
796
797
                           int mpi_rank, int mpi_size, MPI_Comm comm,
                           MPI_Info info) {

798
  hid_t h_file = 0, h_grp = 0;
799
  const size_t Ngas = e->s->nr_parts;
800
  const size_t Nstars = e->s->nr_sparts;
801
  const size_t Ntot = e->s->nr_gparts;
802
803
804
  int periodic = e->s->periodic;
  int numFiles = 1;
  struct part* parts = e->s->parts;
805
806
  struct gpart* gparts = e->s->gparts;
  struct gpart* dmparts = NULL;
807
  struct spart* sparts = e->s->sparts;
808
  static int outputCount = 0;
809
810
  FILE* xmfFile = 0;

811
  /* Number of unassociated gparts */
812
  const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
813

814
  /* File name */
815
  char fileName[FILENAME_BUFFER_SIZE];
816
  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
817
           outputCount);
818
819

  /* First time, we need to create the XMF file */
820
  if (outputCount == 0 && mpi_rank == 0) xmf_create_file(baseName);
821

822
  /* Prepare the XMF file for the new entry */
823
  if (mpi_rank == 0) xmfFile = xmf_prepare_file(baseName);
824

825
  /* Prepare some file-access properties */
826
  hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
827
828
829

  /* Set some MPI-IO parameters */
  // MPI_Info_set(info, "IBM_largeblock_io", "true");
830
831
  MPI_Info_set(info, "romio_cb_write", "enable");
  MPI_Info_set(info, "romio_ds_write", "disable");
832
833
834
835
836
837
838
839

  /* Activate parallel i/o */
  hid_t h_err = H5Pset_fapl_mpio(plist_id, comm, info);
  if (h_err < 0) error("Error setting parallel i/o");
  
  /* Align on 4k pages. */
  h_err = H5Pset_alignment(plist_id, 1024, 4096);
  if (h_err < 0) error("Error setting Hdf5 alignment");
840
841
842
843
844

  /* Disable meta-data cache eviction */
  H5AC_cache_config_t mdc_config;
  mdc_config.version = H5AC__CURR_CACHE_CONFIG_VERSION;
  h_err = H5Pget_mdc_config(plist_id, &mdc_config);
845
  if (h_err < 0) error("Error getting the MDC config");
846
847
848
849

  mdc_config.evictions_enabled = 0; /* false */
  mdc_config.incr_mode = H5C_incr__off;
  mdc_config.decr_mode = H5C_decr__off;
850
  mdc_config.flash_incr_mode = H5C_flash_incr__off;
851
  h_err = H5Pset_mdc_config(plist_id, &mdc_config);
852
  if (h_err < 0) error("Error setting the MDC config");
853

854
  /* Open HDF5 file with the chosen parameters */
855
  h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
856
857
858
  if (h_file < 0) {
    error("Error while opening file '%s'.", fileName);
  }
859

860
  /* Compute offset in the file and total number of particles */
861
862
863
864
865
  size_t N[swift_type_count] = {Ngas, Ndm, 0, 0, Nstars, 0};
  long long N_total[swift_type_count] = {0};
  long long offset[swift_type_count] = {0};
  MPI_Exscan(&N, &offset, swift_type_count, MPI_LONG_LONG_INT, MPI_SUM, comm);
  for (int ptype = 0; ptype < swift_type_count; ++ptype)
866
867
    N_total[ptype] = offset[ptype] + N[ptype];

Matthieu Schaller's avatar
Matthieu Schaller committed
868
869
  /* The last rank now has the correct N_total. Let's
   * broadcast from there */
870
  MPI_Bcast(&N_total, 6, MPI_LONG_LONG_INT, mpi_size - 1, comm);
871

Matthieu Schaller's avatar
Matthieu Schaller committed
872
  /* Now everybody konws its offset and the total number of
873
   * particles of each type */
874

Matthieu Schaller's avatar
Matthieu Schaller committed
875
876
  /* Write the part of the XMF file corresponding to this
   * specific output */
877
  if (mpi_rank == 0) xmf_write_outputheader(xmfFile, fileName, e->time);
878

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

  /* Write the relevant information */
886
  io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
887
888
889

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

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

896
  /* Print the relevant information and print status */
897
  io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
898
  double dblTime = e->time;
899
  io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
900
  int dimension = (int)hydro_dimension;
901
  io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
Matthieu Schaller's avatar
Matthieu Schaller committed
902

903
  /* GADGET-2 legacy values */
904
  /* Number of particles of each type */
905
906
907
  unsigned int numParticles[swift_type_count] = {0};
  unsigned int numParticlesHighWord[swift_type_count] = {0};
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
908
909
910
    numParticles[ptype] = (unsigned int)N_total[ptype];
    numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
  }
911
912
913
914
915
916
  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);
917
  double MassTable[6] = {0., 0., 0., 0., 0., 0.};
918
919
  io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
  unsigned int flagEntropy[swift_type_count] = {0};
920
  flagEntropy[0] = writeEntropyFlag();
921
922
923
  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
924

925
926
  /* Close header */
  H5Gclose(h_grp);
Matthieu Schaller's avatar
Matthieu Schaller committed
927

928
  /* Print the code version */
929
  io_write_code_description(h_file);
930

931
  /* Print the SPH parameters */