parallel_io.c 34.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 "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
55
56
/* The current limit of ROMIO (the underlying MPI-IO layer) is 2GB */
#define HDF5_PARALLEL_IO_MAX_BYTES 2000000000ll

57
58
59
60
/**
 * @brief Reads a data array from a given HDF5 group.
 *
 * @param grp The group from which to read.
61
62
63
64
65
66
 * @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 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.
67
 */
68
69
void readArray(hid_t grp, const struct io_props prop, size_t N,
               long long N_total, long long offset,
70
71
               const struct unit_system* internal_units,
               const struct unit_system* ic_units) {
72

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

  /* Check whether the dataspace exists or not */
78
  const htri_t exist = H5Lexists(grp, prop.name, 0);
79
  if (exist < 0) {
80
    error("Error while checking the existence of data set '%s'.", prop.name);
81
  } else if (exist == 0) {
82
83
    if (prop.importance == COMPULSORY) {
      error("Compulsory data set '%s' not present in the file.", prop.name);
84
    } else {
85
86
      for (size_t i = 0; i < N; ++i)
        memset(prop.field + i * prop.partSize, 0, copySize);
87
      return;
88
    }
89
  }
90

Matthieu Schaller's avatar
Matthieu Schaller committed
91
92
93
  /* message("Reading %s '%s' array...", */
  /*         prop.importance == COMPULSORY ? "compulsory" : "optional  ", */
  /*         prop.name); */
94

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

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

105
  /* Allocate temporary buffer */
106
  void* temp = malloc(num_elements * typeSize);
107
  if (temp == NULL) error("Unable to allocate memory for temporary buffer");
108

109
  /* Prepare information for hyper-slab */
110
111
112
  hsize_t shape[2], offsets[2];
  int rank;
  if (prop.dimension > 1) {
113
114
    rank = 2;
    shape[0] = N;
115
    shape[1] = prop.dimension;
116
117
118
    offsets[0] = offset;
    offsets[1] = 0;
  } else {
119
    rank = 2;
120
    shape[0] = N;
121
    shape[1] = 1;
122
123
124
    offsets[0] = offset;
    offsets[1] = 0;
  }
125
126

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

129
  /* Select hyper-slab in file */
130
  const hid_t h_filespace = H5Dget_space(h_data);
131
132
133
  H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape, NULL);

  /* Set collective reading properties */
134
  const hid_t h_plist_id = H5Pcreate(H5P_DATASET_XFER);
135
136
  H5Pset_dxpl_mpio(h_plist_id, H5FD_MPIO_COLLECTIVE);

137
138
139
  /* Read HDF5 dataspace in temporary buffer */
  /* Dirty version that happens to work for vectors but should be improved */
  /* Using HDF5 dataspaces would be better */
140
  const hid_t h_err = H5Dread(h_data, io_hdf5_type(prop.type), h_memspace,
141
                              h_filespace, h_plist_id, temp);
142
  if (h_err < 0) {
143
144
145
146
147
148
149
150
    error("Error while reading data array '%s'.", prop.name);
  }

  /* Unit conversion if necessary */
  const double factor =
      units_conversion_factor(ic_units, internal_units, prop.units);
  if (factor != 1. && exist != 0) {

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

153
    if (io_is_double_precision(prop.type)) {
154
155
156
157
158
159
      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;
    }
160
  }
161
162

  /* Copy temporary buffer to particle data */
163
164
165
  char* temp_c = temp;
  for (size_t i = 0; i < N; ++i)
    memcpy(prop.field + i * prop.partSize, &temp_c[i * copySize], copySize);
166

167
168
  /* Free and close everything */
  free(temp);
169
170
171
  H5Pclose(h_plist_id);
  H5Sclose(h_filespace);
  H5Sclose(h_memspace);
172
173
174
175
176
177
178
179
180
  H5Tclose(h_type);
  H5Dclose(h_data);
}

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

/**
181
 * @brief Writes a chunk of data in an open HDF5 dataset
182
 *
183
 * @param e The #engine we are writing from.
184
 * @param h_data The HDF5 dataset to write to.
185
 * @param h_plist_id the parallel HDF5 properties.
186
 * @param props The #io_props of the field to read.
187
 * @param N The number of particles to write.
188
189
190
 * @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.
191
 */
192
193
194
195
196
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) {

197
  const size_t typeSize = io_sizeof_type(props.type);
198
199
200
  const size_t copySize = typeSize * props.dimension;
  const size_t num_elements = N * props.dimension;

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

Matthieu Schaller's avatar
Matthieu Schaller committed
205
  /* message("Writing '%s' array...", props.name); */
206
207

  /* Allocate temporary buffer */
208
  void* temp = malloc(num_elements * typeSize);
209
  if (temp == NULL) error("Unable to allocate memory for temporary buffer");
210
211

  /* Copy particle data to temporary buffer */
212
213
214
215
216
217
218
219
220
221
222
  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)
223
      temp_f[i] = props.convert_part(e, &props.parts[i]);
224
225
226
227
228

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

    float* temp_f = temp;
    for (size_t i = 0; i < N; ++i)
229
      temp_f[i] = props.convert_gpart(e, &props.gparts[i]);
230
  }
231
232
233
234
235
236

  /* 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
237
    /* message("Converting ! factor=%e", factor); */
238

239
    if (io_is_double_precision(props.type)) {
240
241
242
243
244
245
246
      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;
    }
  }
247
248

  /* Create data space */
249
  const hid_t h_memspace = H5Screate(H5S_SIMPLE);
250
  if (h_memspace < 0) {
251
252
    error("Error while creating data space (memory) for field '%s'.",
          props.name);
253
  }
254

255
256
257
258
  int rank;
  hsize_t shape[2];
  hsize_t offsets[2];
  if (props.dimension > 1) {
259
260
    rank = 2;
    shape[0] = N;
261
    shape[1] = props.dimension;
262
263
264
265
266
267
268
269
270
271
    offsets[0] = offset;
    offsets[1] = 0;
  } else {
    rank = 1;
    shape[0] = N;
    shape[1] = 0;
    offsets[0] = offset;
    offsets[1] = 0;
  }

272
  /* Change shape of memory data space */
273
  hid_t h_err = H5Sset_extent_simple(h_memspace, rank, shape, NULL);
274
275
  if (h_err < 0) {
    error("Error while changing data space (memory) shape for field '%s'.",
276
          props.name);
277
  }
278

279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
  /* Select the hyper-salb corresponding to this rank */
  hid_t h_filespace = H5Dget_space(h_data);
  H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape, NULL);

  /* message("Writing %zd elements = %zd bytes (int=%d) at offset %zd", */
  /* 	  N * props.dimension, N * props.dimension * typeSize, */
  /* 	  (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,
                   h_plist_id, temp);
  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);
}
299
300
301
302
303
/**
 * @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.
304
305
 * @param fileName The name of the file in which the data is written.
 * @param xmfFile The FILE used to write the XMF description.
306
307
308
309
 * @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.
310
311
312
313
314
 * @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.
315
 */
316
317
318
319
320
321
322
323
324
325
326
void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
                char* partTypeGroupName, const struct io_props props, size_t N,
                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];
327
  hsize_t chunk_shape[2];
328
329
330
331
  if (props.dimension > 1) {
    rank = 2;
    shape_total[0] = N_total;
    shape_total[1] = props.dimension;
332
333
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = props.dimension;
334
335
336
337
  } else {
    rank = 1;
    shape_total[0] = N_total;
    shape_total[1] = 0;
338
339
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = 0;
340
341
  }

342
343
344
345
  /* Make sure the chunks are not larger than the dataset */
  if (chunk_shape[0] > N_total) chunk_shape[0] = N_total;

  /* Create the space in the file */
346
347
348
349
  hid_t h_filespace = H5Screate(H5S_SIMPLE);
  if (h_filespace < 0) {
    error("Error while creating data space (file) for field '%s'.", props.name);
  }
350

351
  /* Change shape of file data space */
352
  hid_t h_err = H5Sset_extent_simple(h_filespace, rank, shape_total, NULL);
353
  if (h_err < 0) {
354
355
    error("Error while changing data space (file) shape for field '%s'.",
          props.name);
356
357
  }

358
359
360
361
362
363
364
365
366
367
  /* Dataset properties */
  const hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);

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

368
  /* Create dataset */
369
370
  const hid_t h_data = H5Dcreate(grp, props.name, io_hdf5_type(props.type),
                                 h_filespace, H5P_DEFAULT, h_prop, H5P_DEFAULT);
371
  if (h_data < 0) {
372
    error("Error while creating dataset '%s'.", props.name);
373
374
  }

375
376
377
  H5Sclose(h_filespace);

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

381
382
383
  /* 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;
384
385
386
387
388
389
  while (redo) {

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

390
    /* Write the first chunk */
391
    const size_t this_chunk = (N > max_chunk_size) ? max_chunk_size : N;
392
    writeArray_chunk(e, h_data, h_plist_id, props, this_chunk, offset,
393
                     internal_units, snapshot_units);
394
395

    /* Compute how many items are left */
396
    if (N > max_chunk_size) {
397
398
399
      N -= max_chunk_size;
      offset += max_chunk_size;
      redo = 1;
400
    } else {
401
402
403
404
405
      N = 0;
      offset += 0;
      redo = 0;
    }

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

409
    /* Do we need to run again ? */
410
411
    MPI_Allreduce(MPI_IN_PLACE, &redo, 1, MPI_SIGNED_CHAR, MPI_MAX,
                  MPI_COMM_WORLD);
412
  }
413
414

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

419
  /* Write unit conversion factors for this data set */
420
421
  char buffer[FIELD_BUFFER_SIZE];
  units_cgs_conversion_string(buffer, snapshot_units, props.units);
422
423
424
425
426
427
428
429
  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);
430

431
432
  /* Close everything */
  H5Pclose(h_prop);
433
434
435
436
  H5Dclose(h_data);
  H5Pclose(h_plist_id);
}

437
438
439
440
/**
 * @brief Reads an HDF5 initial condition file (GADGET-3 type) in parallel
 *
 * @param fileName The file to read.
441
 * @param internal_units The system units used internally
442
443
 * @param dim (output) The dimension of the volume read from the file.
 * @param parts (output) The array of #part read from the file.
444
445
446
447
448
 * @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.
449
 * @param periodic (output) 1 if the volume is periodic, 0 if not.
450
 * @param flag_entropy (output) 1 if the ICs contained Entropy in the
451
 * InternalEnergy field
452
453
454
 * @param with_hydro Are we running with hydro ?
 * @param with_gravity Are we running with gravity ?
 * @param with_stars Are we running with stars ?
455
456
457
458
 * @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
459
 * @param dry_run If 1, don't read the particle. Only allocates the arrays.
460
461
 *
 */
462
void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
463
                      double dim[3], struct part** parts, struct gpart** gparts,
Matthieu Schaller's avatar
Matthieu Schaller committed
464
465
466
467
468
                      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) {
469

470
  hid_t h_file = 0, h_grp = 0;
471
  /* GADGET has only cubic boxes (in cosmological mode) */
Matthieu Schaller's avatar
Matthieu Schaller committed
472
  double boxSize[3] = {0.0, -1.0, -1.0};
473
474
475
476
477
  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};
478
  int dimension = 3; /* Assume 3D if nothing is specified */
479
  size_t Ndm = 0;
480
481
482
483
484
485
486
487
488
489
490
491

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

  /* Read the relevant information */
496
  io_read_attribute(h_grp, "PeriodicBoundariesOn", INT, periodic);
497
498
499
500
501
502

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

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

506
507
508
509
  /* 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");
510
  if (hid_dim > 0) io_read_attribute(h_grp, "Dimension", INT, &dimension);
511
512
513
514
  if (dimension != hydro_dimension)
    error("ICs dimensionality (%dD) does not match code dimensionality (%dD)",
          dimension, (int)hydro_dimension);

515
  /* Read the relevant information and print status */
516
  int flag_entropy_temp[6];
517
  io_read_attribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy_temp);
518
  *flag_entropy = flag_entropy_temp[0];
519
520
521
522
  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);
523

524
  for (int ptype = 0; ptype < swift_type_count; ++ptype)
Matthieu Schaller's avatar
Matthieu Schaller committed
525
526
    N_total[ptype] =
        (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);
527

528
  /* Get the box size if not cubic */
529
530
531
532
  dim[0] = boxSize[0];
  dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
  dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];

533
534
535
536
537
538
  /* 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];

539
540
541
  /* message("Found %lld particles in a %speriodic box of size [%f %f %f].", */
  /* 	  N_total[0], (periodic ? "": "non-"), dim[0], */
  /* 	  dim[1], dim[2]); */
542
543

  /* Divide the particles among the tasks. */
544
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
545
546
547
    offset[ptype] = mpi_rank * N_total[ptype] / mpi_size;
    N[ptype] = (mpi_rank + 1) * N_total[ptype] / mpi_size - offset[ptype];
  }
548
549
550
551

  /* Close header */
  H5Gclose(h_grp);

552
  /* Read the unit system used in the ICs */
553
  struct unit_system* ic_units = malloc(sizeof(struct unit_system));
554
  if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
555
  io_read_unit_system(h_file, ic_units);
556
557
558
559
560

  /* 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
561
      message("IC and internal units match. No conversion needed.");
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587

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

588
  /* Convert the dimensions of the box */
Matthieu Schaller's avatar
Matthieu Schaller committed
589
590
591
  for (int j = 0; j < 3; j++)
    dim[j] *=
        units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
592

593
  /* Allocate memory to store SPH particles */
594
595
  if (with_hydro) {
    *Ngas = N[0];
Matthieu Schaller's avatar
Matthieu Schaller committed
596
597
    if (posix_memalign((void*)parts, part_align,
                       (*Ngas) * sizeof(struct part)) != 0)
598
599
600
601
602
603
      error("Error while allocating memory for particles");
    bzero(*parts, *Ngas * sizeof(struct part));
  }

  /* Allocate memory to store star particles */
  if (with_stars) {
604
    *Nstars = N[swift_type_star];
605
606
607
608
609
610
611
612
613
    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];
614
615
616
    *Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
               N[swift_type_dark_matter] +
               (with_stars ? N[swift_type_star] : 0);
617
    if (posix_memalign((void*)gparts, gpart_align,
Matthieu Schaller's avatar
Matthieu Schaller committed
618
                       *Ngparts * sizeof(struct gpart)) != 0)
619
620
621
622
623
      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) /
624
625
   * (1024.*1024.)); */

626
  /* message("BoxSize = %lf", dim[0]); */
627
  /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */
628
629

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

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

635
636
637
    /* 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
638
             ptype);
639
640
641
642
    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
643

644
645
    int num_fields = 0;
    struct io_props list[100];
646
    size_t Nparticles = 0;
647

648
649
    /* Read particle fields into the particle structure */
    switch (ptype) {
Matthieu Schaller's avatar
Matthieu Schaller committed
650

651
      case swift_type_gas:
Matthieu Schaller's avatar
Matthieu Schaller committed
652
653
654
655
        if (with_hydro) {
          Nparticles = *Ngas;
          hydro_read_particles(*parts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
656
657
        break;

658
      case swift_type_dark_matter:
Matthieu Schaller's avatar
Matthieu Schaller committed
659
660
661
662
        if (with_gravity) {
          Nparticles = Ndm;
          darkmatter_read_particles(*gparts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
663
664
        break;

665
      case swift_type_star:
666
        if (with_stars) {
Matthieu Schaller's avatar
Matthieu Schaller committed
667
668
669
670
          Nparticles = *Nstars;
          star_read_particles(*sparts, list, &num_fields);
        }
        break;
671

Matthieu Schaller's avatar
Matthieu Schaller committed
672
      default:
673
        message("Particle Type %d not yet supported. Particles ignored", ptype);
674
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
675

676
677
678
    /* Read everything */
    if (!dry_run)
      for (int i = 0; i < num_fields; ++i)
679
        readArray(h_grp, list[i], Nparticles, N_total[ptype], offset[ptype],
680
681
                  internal_units, ic_units);

682
683
684
    /* Close particle group */
    H5Gclose(h_grp);
  }
685

686
  /* Prepare the DM particles */
687
  if (!dry_run && with_gravity) io_prepare_dm_gparts(*gparts, Ndm);
688

689
690
  /* Duplicate the hydro particles into gparts */
  if (!dry_run && with_gravity && with_hydro)
691
    io_duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
692
693
694

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

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

699
700
701
  /* Clean up */
  free(ic_units);

702
703
704
705
706
707
708
  /* Close property handler */
  H5Pclose(h_plist_id);

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

709
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
710
711
 * @brief Writes an HDF5 output file (GADGET-3 type) with
 *its XMF descriptor
712
713
 *
 * @param e The engine containing all the system.
714
 * @param baseName The common part of the snapshot file name.
715
716
 * @param internal_units The #unit_system used internally
 * @param snapshot_units The #unit_system used in the snapshots
717
718
719
720
 * @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
721
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
722
 * Creates an HDF5 output file and writes the particles
723
724
 * contained in the engine. If such a file already exists, it is
 * erased and replaced by the new one.
725
726
727
728
729
 * The companion XMF file is also updated accordingly.
 *
 * Calls #error() if an error occurs.
 *
 */
730
void write_output_parallel(struct engine* e, const char* baseName,
731
732
                           const struct unit_system* internal_units,
                           const struct unit_system* snapshot_units,
733
734
735
                           int mpi_rank, int mpi_size, MPI_Comm comm,
                           MPI_Info info) {

736
  hid_t h_file = 0, h_grp = 0;
737
  const size_t Ngas = e->s->nr_parts;
738
  const size_t Nstars = e->s->nr_sparts;
739
  const size_t Ntot = e->s->nr_gparts;
740
741
742
  int periodic = e->s->periodic;
  int numFiles = 1;
  struct part* parts = e->s->parts;
743
744
  struct gpart* gparts = e->s->gparts;
  struct gpart* dmparts = NULL;
745
  struct spart* sparts = e->s->sparts;
746
  static int outputCount = 0;
747
748
  FILE* xmfFile = 0;

749
  /* Number of unassociated gparts */
750
  const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
751

752
  /* File name */
753
  char fileName[FILENAME_BUFFER_SIZE];
754
  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
755
           outputCount);
756
757

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

760
  /* Prepare the XMF file for the new entry */
761
  if (mpi_rank == 0) xmfFile = xmf_prepare_file(baseName);
762

763
  /* Open HDF5 file */
764
765
766
  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);
767
768
769
  if (h_file < 0) {
    error("Error while opening file '%s'.", fileName);
  }
770

Matthieu Schaller's avatar
Matthieu Schaller committed
771
772
  /* Compute offset in the file and total number of
   * particles */
773
774
775
776
777
  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)
778
779
    N_total[ptype] = offset[ptype] + N[ptype];

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

Matthieu Schaller's avatar
Matthieu Schaller committed
784
785
  /* Now everybody konws its offset and the total number of
   * particles of each
786
   * type */
787

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

792
793
  /* Open header to write simulation properties */
  /* message("Writing runtime parameters..."); */
Matthieu Schaller's avatar
Matthieu Schaller committed
794
795
  h_grp =
      H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
796
  if (h_grp < 0) error("Error while creating runtime parameters group\n");
797
798

  /* Write the relevant information */
799
  io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
800
801
802

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

804
805
  /* Open header to write simulation properties */
  /* message("Writing file header..."); */
806
  h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
807
  if (h_grp < 0) error("Error while creating file header\n");
Matthieu Schaller's avatar
Matthieu Schaller committed
808

809
  /* Print the relevant information and print status */
810
  io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
811
  double dblTime = e->time;
812
  io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
813
  int dimension = (int)hydro_dimension;
814
  io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
Matthieu Schaller's avatar
Matthieu Schaller committed
815

816
  /* GADGET-2 legacy values */
817
  /* Number of particles of each type */
818
819
820
  unsigned int numParticles[swift_type_count] = {0};
  unsigned int numParticlesHighWord[swift_type_count] = {0};
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
821
822
823
    numParticles[ptype] = (unsigned int)N_total[ptype];
    numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
  }
824
825
826
827
828
829
  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);
830
  double MassTable[6] = {0., 0., 0., 0., 0., 0.};
831
832
  io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
  unsigned int flagEntropy[swift_type_count] = {0};
833
  flagEntropy[0] = writeEntropyFlag();
834
835
836
  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
837

838
839
  /* Close header */
  H5Gclose(h_grp);
Matthieu Schaller's avatar
Matthieu Schaller committed
840

841
  /* Print the code version */
842
  io_write_code_description(h_file);
843

844
  /* Print the SPH parameters */
845
846
847
  if (e->policy & engine_policy_hydro) {
    h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
848
849
850
851
852
    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);
  }
853

854
  /* Print the gravity parameters */
855
  if (e->policy & engine_policy_self_gravity) {
856
    h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
857
                      H5P_DEFAULT);
858
859
860
861
    if (h_grp < 0) error("Error while creating gravity group");
    gravity_props_print_snapshot(h_grp, e->gravity_properties);
    H5Gclose(h_grp);
  }
862
863
864
865
866
867
868

  /* 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);
869

870
  /* Print the system of Units used in the spashot */
871
  io_write_unit_system(h_file, snapshot_units, "Units");
872
873

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

876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
  /* Tell the user if a conversion will be needed */
  if (e->verbose && mpi_rank == 0) {
    if (units_are_equal(snapshot_units, internal_units)) {

      message("Snapshot and internal units match. No conversion needed.");

    } else {

      message("Conversion needed from:");
      message("(Snapshot) Unit system: U_M =      %e g.",
              snapshot_units->UnitMass_in_cgs);
      message("(Snapshot) Unit system: U_L =      %e cm.",
              snapshot_units->UnitLength_in_cgs);
      message("(Snapshot) Unit system: U_t =      %e s.",
              snapshot_units->UnitTime_in_cgs);
      message("(Snapshot) Unit system: U_I =      %e A.",
              snapshot_units->UnitCurrent_in_cgs);
      message("(Snapshot) Unit system: U_T =      %e K.",
              snapshot_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);
    }
  }

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

912
913
    /* Don't do anything if no particle of this kind */
    if (N_total[ptype] == 0) continue;
Matthieu Schaller's avatar
Matthieu Schaller committed
914
915
916

    /* Add the global information for that particle type to
     * the XMF meta-file */
917
    if (mpi_rank == 0)
918
919
      xmf_write_groupheader(xmfFile, fileName, N_total[ptype],
                            (enum part_type)ptype);
Matthieu Schaller's avatar
Matthieu Schaller committed
920

921
922
923
    /* 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
924
925
926
             ptype);
    h_grp = H5Gcreate(h_file, partTypeGroupName, H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
927
928
929
    if (h_grp < 0) {
      error("Error while opening particle group %s.", partTypeGroupName);
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
930

931
932
    int num_fields = 0;
    struct io_props list[100];
933
    size_t Nparticles = 0;