parallel_io.c 37.2 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
149
150
151
152
153
154
155
156
}

/**
 * @brief Reads a data array from a given HDF5 group.
 *
 * @param grp The group from which to read.
 * @param prop The #io_props of the field to read.
 * @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
               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;
    }
  }

  /* Work out properties of the array in the file */
  int rank;
  hsize_t shape_total[2];
  if (props.dimension > 1) {
    rank = 2;
    shape_total[0] = N_total;
    shape_total[1] = props.dimension;
  } else {
    rank = 1;
    shape_total[0] = N_total;
    shape_total[1] = 0;
  }
191

192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
  /* 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 =
213
        HDF5_PARALLEL_IO_MAX_BYTES / (props.dimension * typeSize);
214
215
216
217

    /* 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,
218
                    internal_units, ic_units);
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241

    /* 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);
242
243
244
245
246
247
248
249
250
  H5Tclose(h_type);
  H5Dclose(h_data);
}

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

/**
251
 * @brief Writes a chunk of data in an open HDF5 dataset
252
 *
253
 * @param e The #engine we are writing from.
254
 * @param h_data The HDF5 dataset to write to.
255
 * @param h_plist_id the parallel HDF5 properties.
256
 * @param props The #io_props of the field to read.
257
 * @param N The number of particles to write.
258
259
260
 * @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.
261
 */
262
263
264
265
266
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) {

267
  const size_t typeSize = io_sizeof_type(props.type);
268
269
270
  const size_t copySize = typeSize * props.dimension;
  const size_t num_elements = N * props.dimension;

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

Matthieu Schaller's avatar
Matthieu Schaller committed
275
  /* message("Writing '%s' array...", props.name); */
276
277

  /* Allocate temporary buffer */
278
  void* temp = malloc(num_elements * typeSize);
279
  if (temp == NULL) error("Unable to allocate memory for temporary buffer");
280
281

  /* Copy particle data to temporary buffer */
282
283
284
285
286
287
288
289
290
291
292
  if (props.convert_part == NULL &&
      props.convert_gpart == NULL) { /* No conversion */

    char* temp_c = temp;
    for (size_t i = 0; i < N; ++i)
      memcpy(&temp_c[i * copySize], props.field + i * props.partSize, copySize);

  } else if (props.convert_part != NULL) { /* conversion (for parts)*/

    float* temp_f = temp;
    for (size_t i = 0; i < N; ++i)
293
      temp_f[i] = props.convert_part(e, &props.parts[i]);
294
295
296
297
298

  } else if (props.convert_gpart != NULL) { /* conversion (for gparts)*/

    float* temp_f = temp;
    for (size_t i = 0; i < N; ++i)
299
      temp_f[i] = props.convert_gpart(e, &props.gparts[i]);
300
  }
301
302
303
304
305
306

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

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

309
    if (io_is_double_precision(props.type)) {
310
311
312
313
314
315
316
      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;
    }
  }
317
318

  /* Create data space */
319
  const hid_t h_memspace = H5Screate(H5S_SIMPLE);
320
  if (h_memspace < 0) {
321
322
    error("Error while creating data space (memory) for field '%s'.",
          props.name);
323
  }
324

325
326
327
328
  int rank;
  hsize_t shape[2];
  hsize_t offsets[2];
  if (props.dimension > 1) {
329
330
    rank = 2;
    shape[0] = N;
331
    shape[1] = props.dimension;
332
333
334
335
336
337
338
339
340
341
    offsets[0] = offset;
    offsets[1] = 0;
  } else {
    rank = 1;
    shape[0] = N;
    shape[1] = 0;
    offsets[0] = offset;
    offsets[1] = 0;
  }

342
  /* Change shape of memory data space */
343
  hid_t h_err = H5Sset_extent_simple(h_memspace, rank, shape, NULL);
344
345
  if (h_err < 0) {
    error("Error while changing data space (memory) shape for field '%s'.",
346
          props.name);
347
  }
348

349
350
  /* Select the hyper-salb corresponding to this rank */
  hid_t h_filespace = H5Dget_space(h_data);
Matthieu Schaller's avatar
Matthieu Schaller committed
351
352
353
  if (N > 0) {
    H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape,
                        NULL);
354
355
356
  } else {
    H5Sselect_none(h_filespace);
  }
357

Matthieu Schaller's avatar
Matthieu Schaller committed
358
359
  /* message("Writing %lld '%s', %zd elements = %zd bytes (int=%d) at offset
   * %zd", */
360
  /* 	  N, props.name, N * props.dimension, N * props.dimension * typeSize, */
361
362
363
364
  /* 	  (int)(N * props.dimension * typeSize), offset); */

  /* 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
365
                   h_plist_id, temp);
366
367
368
369
370
371
372
373
374
  if (h_err < 0) {
    error("Error while writing data array '%s'.", props.name);
  }

  /* Free and close everything */
  free(temp);
  H5Sclose(h_memspace);
  H5Sclose(h_filespace);
}
375

376
377
378
379
380
/**
 * @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.
381
382
 * @param fileName The name of the file in which the data is written.
 * @param xmfFile The FILE used to write the XMF description.
383
384
385
386
 * @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.
387
388
389
390
391
 * @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.
392
 */
393
void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
394
                char* partTypeGroupName, struct io_props props, size_t N,
395
396
397
398
399
400
401
402
403
                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);

  /* Work out properties of the array in the file */
  int rank;
  hsize_t shape_total[2];
404
  hsize_t chunk_shape[2];
405
406
407
408
  if (props.dimension > 1) {
    rank = 2;
    shape_total[0] = N_total;
    shape_total[1] = props.dimension;
409
410
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = props.dimension;
411
412
413
414
  } else {
    rank = 1;
    shape_total[0] = N_total;
    shape_total[1] = 0;
415
416
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = 0;
417
418
  }

419
  /* Make sure the chunks are not larger than the dataset */
lhausamm's avatar
lhausamm committed
420
  if (chunk_shape[0] > (hsize_t)N_total) chunk_shape[0] = N_total;
421
422

  /* Create the space in the file */
423
424
425
426
  hid_t h_filespace = H5Screate(H5S_SIMPLE);
  if (h_filespace < 0) {
    error("Error while creating data space (file) for field '%s'.", props.name);
  }
427

428
  /* Change shape of file data space */
429
  hid_t h_err = H5Sset_extent_simple(h_filespace, rank, shape_total, NULL);
430
  if (h_err < 0) {
431
432
    error("Error while changing data space (file) shape for field '%s'.",
          props.name);
433
434
  }

435
436
437
438
  /* Dataset properties */
  const hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);

  /* Set chunk size */
439
440
441
442
443
  /* 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); */
  /* } */
444

445
  /* Create dataset */
446
447
  const hid_t h_data = H5Dcreate(grp, props.name, io_hdf5_type(props.type),
                                 h_filespace, H5P_DEFAULT, h_prop, H5P_DEFAULT);
448
  if (h_data < 0) {
449
    error("Error while creating dataset '%s'.", props.name);
450
451
  }

452
453
454
  H5Sclose(h_filespace);

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

458
459
460
  /* 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;
461
462
463
464
465
466
  while (redo) {

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

467
    /* Write the first chunk */
468
    const size_t this_chunk = (N > max_chunk_size) ? max_chunk_size : N;
469
    writeArray_chunk(e, h_data, h_plist_id, props, this_chunk, offset,
470
                     internal_units, snapshot_units);
471
472

    /* Compute how many items are left */
473
    if (N > max_chunk_size) {
474
      N -= max_chunk_size;
475
      props.field += max_chunk_size * props.partSize;
476
477
      offset += max_chunk_size;
      redo = 1;
478
    } else {
479
480
481
482
483
484
      N = 0;
      offset += 0;
      redo = 0;
    }

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

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

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

497
  /* Write unit conversion factors for this data set */
498
499
  char buffer[FIELD_BUFFER_SIZE];
  units_cgs_conversion_string(buffer, snapshot_units, props.units);
500
501
502
503
504
505
506
507
  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);
508

509
510
  /* Close everything */
  H5Pclose(h_prop);
511
512
513
514
  H5Dclose(h_data);
  H5Pclose(h_plist_id);
}

515
516
517
518
/**
 * @brief Reads an HDF5 initial condition file (GADGET-3 type) in parallel
 *
 * @param fileName The file to read.
519
 * @param internal_units The system units used internally
520
521
 * @param dim (output) The dimension of the volume read from the file.
 * @param parts (output) The array of #part read from the file.
522
523
524
525
526
 * @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.
527
 * @param periodic (output) 1 if the volume is periodic, 0 if not.
528
 * @param flag_entropy (output) 1 if the ICs contained Entropy in the
529
 * InternalEnergy field
530
531
532
 * @param with_hydro Are we running with hydro ?
 * @param with_gravity Are we running with gravity ?
 * @param with_stars Are we running with stars ?
533
534
535
536
 * @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
537
 * @param dry_run If 1, don't read the particle. Only allocates the arrays.
538
539
 *
 */
540
void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
541
                      double dim[3], struct part** parts, struct gpart** gparts,
Matthieu Schaller's avatar
Matthieu Schaller committed
542
543
544
545
546
                      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) {
547

548
  hid_t h_file = 0, h_grp = 0;
549
  /* GADGET has only cubic boxes (in cosmological mode) */
Matthieu Schaller's avatar
Matthieu Schaller committed
550
  double boxSize[3] = {0.0, -1.0, -1.0};
551
552
553
554
555
  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};
556
  int dimension = 3; /* Assume 3D if nothing is specified */
557
  size_t Ndm = 0;
558
559
560
561
562
563
564
565
566
567
568
569

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

  /* Read the relevant information */
574
  io_read_attribute(h_grp, "PeriodicBoundariesOn", INT, periodic);
575
576
577
578
579
580

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

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

584
585
586
587
  /* 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");
588
  if (hid_dim > 0) io_read_attribute(h_grp, "Dimension", INT, &dimension);
589
590
591
592
  if (dimension != hydro_dimension)
    error("ICs dimensionality (%dD) does not match code dimensionality (%dD)",
          dimension, (int)hydro_dimension);

593
  /* Read the relevant information and print status */
594
  int flag_entropy_temp[6];
595
  io_read_attribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy_temp);
596
  *flag_entropy = flag_entropy_temp[0];
597
598
599
600
  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);
601

602
  for (int ptype = 0; ptype < swift_type_count; ++ptype)
Matthieu Schaller's avatar
Matthieu Schaller committed
603
604
    N_total[ptype] =
        (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);
605

606
  /* Get the box size if not cubic */
607
608
609
610
  dim[0] = boxSize[0];
  dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
  dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];

611
612
613
614
615
616
  /* 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];

617
618
619
  /* message("Found %lld particles in a %speriodic box of size [%f %f %f].", */
  /* 	  N_total[0], (periodic ? "": "non-"), dim[0], */
  /* 	  dim[1], dim[2]); */
620
621

  /* Divide the particles among the tasks. */
622
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
623
624
625
    offset[ptype] = mpi_rank * N_total[ptype] / mpi_size;
    N[ptype] = (mpi_rank + 1) * N_total[ptype] / mpi_size - offset[ptype];
  }
626
627
628
629

  /* Close header */
  H5Gclose(h_grp);

630
  /* Read the unit system used in the ICs */
631
  struct unit_system* ic_units = malloc(sizeof(struct unit_system));
632
  if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
633
  io_read_unit_system(h_file, ic_units);
634
635
636
637
638

  /* 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
639
      message("IC and internal units match. No conversion needed.");
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665

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

666
  /* Convert the dimensions of the box */
Matthieu Schaller's avatar
Matthieu Schaller committed
667
668
669
  for (int j = 0; j < 3; j++)
    dim[j] *=
        units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
670

671
  /* Allocate memory to store SPH particles */
672
673
  if (with_hydro) {
    *Ngas = N[0];
Matthieu Schaller's avatar
Matthieu Schaller committed
674
675
    if (posix_memalign((void*)parts, part_align,
                       (*Ngas) * sizeof(struct part)) != 0)
676
677
678
679
680
681
      error("Error while allocating memory for particles");
    bzero(*parts, *Ngas * sizeof(struct part));
  }

  /* Allocate memory to store star particles */
  if (with_stars) {
682
    *Nstars = N[swift_type_star];
683
684
685
686
687
688
689
690
691
    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];
692
693
694
    *Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
               N[swift_type_dark_matter] +
               (with_stars ? N[swift_type_star] : 0);
695
    if (posix_memalign((void*)gparts, gpart_align,
Matthieu Schaller's avatar
Matthieu Schaller committed
696
                       *Ngparts * sizeof(struct gpart)) != 0)
697
698
699
700
701
      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) /
702
703
   * (1024.*1024.)); */

704
  /* message("BoxSize = %lf", dim[0]); */
705
  /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */
706
707

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

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

713
714
715
    /* 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
716
             ptype);
717
718
719
720
    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
721

722
723
    int num_fields = 0;
    struct io_props list[100];
724
    size_t Nparticles = 0;
725

726
727
    /* Read particle fields into the particle structure */
    switch (ptype) {
Matthieu Schaller's avatar
Matthieu Schaller committed
728

729
      case swift_type_gas:
Matthieu Schaller's avatar
Matthieu Schaller committed
730
731
732
733
        if (with_hydro) {
          Nparticles = *Ngas;
          hydro_read_particles(*parts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
734
735
        break;

736
      case swift_type_dark_matter:
Matthieu Schaller's avatar
Matthieu Schaller committed
737
738
739
740
        if (with_gravity) {
          Nparticles = Ndm;
          darkmatter_read_particles(*gparts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
741
742
        break;

743
      case swift_type_star:
744
        if (with_stars) {
Matthieu Schaller's avatar
Matthieu Schaller committed
745
746
747
748
          Nparticles = *Nstars;
          star_read_particles(*sparts, list, &num_fields);
        }
        break;
749

Matthieu Schaller's avatar
Matthieu Schaller committed
750
      default:
751
        message("Particle Type %d not yet supported. Particles ignored", ptype);
752
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
753

754
755
756
    /* Read everything */
    if (!dry_run)
      for (int i = 0; i < num_fields; ++i)
757
758
        readArray(h_grp, list[i], Nparticles, N_total[ptype], mpi_rank,
                  offset[ptype], internal_units, ic_units);
759

760
761
762
    /* Close particle group */
    H5Gclose(h_grp);
  }
763

764
  /* Prepare the DM particles */
765
  if (!dry_run && with_gravity) io_prepare_dm_gparts(*gparts, Ndm);
766

767
768
  /* Duplicate the hydro particles into gparts */
  if (!dry_run && with_gravity && with_hydro)
769
    io_duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
770
771
772

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

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

777
778
779
  /* Clean up */
  free(ic_units);

780
781
782
783
784
785
786
  /* Close property handler */
  H5Pclose(h_plist_id);

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

787
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
788
789
 * @brief Writes an HDF5 output file (GADGET-3 type) with
 *its XMF descriptor
790
791
 *
 * @param e The engine containing all the system.
792
 * @param baseName The common part of the snapshot file name.
793
794
 * @param internal_units The #unit_system used internally
 * @param snapshot_units The #unit_system used in the snapshots
795
796
797
798
 * @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
799
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
800
 * Creates an HDF5 output file and writes the particles
801
802
 * contained in the engine. If such a file already exists, it is
 * erased and replaced by the new one.
803
804
805
806
807
 * The companion XMF file is also updated accordingly.
 *
 * Calls #error() if an error occurs.
 *
 */
808
void write_output_parallel(struct engine* e, const char* baseName,
809
810
                           const struct unit_system* internal_units,
                           const struct unit_system* snapshot_units,
811
812
813
                           int mpi_rank, int mpi_size, MPI_Comm comm,
                           MPI_Info info) {

814
  hid_t h_file = 0, h_grp = 0;
815
  const size_t Ngas = e->s->nr_parts;
816
  const size_t Nstars = e->s->nr_sparts;
817
  const size_t Ntot = e->s->nr_gparts;
818
819
820
  int periodic = e->s->periodic;
  int numFiles = 1;
  struct part* parts = e->s->parts;
821
822
  struct gpart* gparts = e->s->gparts;
  struct gpart* dmparts = NULL;
823
  struct spart* sparts = e->s->sparts;
824
  static int outputCount = 0;
825
826
  FILE* xmfFile = 0;

827
  /* Number of unassociated gparts */
828
  const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
829

830
  /* File name */
831
  char fileName[FILENAME_BUFFER_SIZE];
832
  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
833
           outputCount);
834
835

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

838
  /* Prepare the XMF file for the new entry */
839
  if (mpi_rank == 0) xmfFile = xmf_prepare_file(baseName);
840

841
  /* Open HDF5 file */
842
843
844
  hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
  H5Pset_fapl_mpio(plist_id, comm, info);
  h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
845
846
847
  if (h_file < 0) {
    error("Error while opening file '%s'.", fileName);
  }
848

Matthieu Schaller's avatar
Matthieu Schaller committed
849
850
  /* Compute offset in the file and total number of
   * particles */
851
852
853
854
855
  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)
856
857
    N_total[ptype] = offset[ptype] + N[ptype];

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

Matthieu Schaller's avatar
Matthieu Schaller committed
862
863
  /* Now everybody konws its offset and the total number of
   * particles of each
864
   * type */
865

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

870
871
  /* Open header to write simulation properties */
  /* message("Writing runtime parameters..."); */
Matthieu Schaller's avatar
Matthieu Schaller committed
872
873
  h_grp =
      H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
874
  if (h_grp < 0) error("Error while creating runtime parameters group\n");
875
876

  /* Write the relevant information */
877
  io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
878
879
880

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

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

887
  /* Print the relevant information and print status */
888
  io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
889
  double dblTime = e->time;
890
  io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
891
  int dimension = (int)hydro_dimension;
892
  io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
Matthieu Schaller's avatar
Matthieu Schaller committed
893

894
  /* GADGET-2 legacy values */
895
  /* Number of particles of each type */
896
897
898
  unsigned int numParticles[swift_type_count] = {0};
  unsigned int numParticlesHighWord[swift_type_count] = {0};
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
899
900
901
    numParticles[ptype] = (unsigned int)N_total[ptype];
    numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
  }
902
903
904
905
906
907
  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);
908
  double MassTable[6] = {0., 0., 0., 0., 0., 0.};
909
910
  io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
  unsigned int flagEntropy[swift_type_count] = {0};
911
  flagEntropy[0] = writeEntropyFlag();
912
913
914
  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
915

916
917
  /* Close header */
  H5Gclose(h_grp);
Matthieu Schaller's avatar
Matthieu Schaller committed
918

919
  /* Print the code version */
920
  io_write_code_description(h_file);
921

922
  /* Print the SPH parameters */
923
924
925
  if (e->policy & engine_policy_hydro) {
    h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
926
927
928
929
930
    if (h_grp < 0) error("Error while creating SPH group");
    hydro_props_print_snapshot(h_grp, e->hydro_properties);
    writeSPHflavour(h_grp);
    H5Gclose(h_grp);
  }
931

932
  /* Print the gravity parameters */
933
  if (e->policy & engine_policy_self_gravity) {
934
    h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
935
                      H5P_DEFAULT);
936
937
938
939
    if (h_grp < 0) error("Error while creating gravity group");
    gravity_props_print_snapshot(h_grp, e->gravity_properties);
    H5Gclose(h_grp);
  }
940
941
942
943
944
945
946

  /* Print the runtime parameters */
  h_grp =
      H5Gcreate(h_file, "/Parameters", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
  if (h_grp < 0) error("Error while creating parameters group");
  parser_write_params_to_hdf5(e->parameter_file, h_grp);
  H5Gclose(h_grp);
947

948
  /* Print the system of Units used in the spashot */
949
  io_write_unit_system(h_file, snapshot_units, "Units");
950
951

  /* Print the system of Units used internally */
952
  io_write_unit_system(h_file, internal_units, "InternalCodeUnits");