parallel_io.c 36 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
  /* Copy the particle data to the temporary buffer */
  io_copy_temp_buffer(temp, e, props, N, internal_units, snapshot_units);
271
272

  /* Create data space */
273
  const hid_t h_memspace = H5Screate(H5S_SIMPLE);
274
  if (h_memspace < 0) {
275
276
    error("Error while creating data space (memory) for field '%s'.",
          props.name);
277
  }
278

279
280
281
282
  int rank;
  hsize_t shape[2];
  hsize_t offsets[2];
  if (props.dimension > 1) {
283
284
    rank = 2;
    shape[0] = N;
285
    shape[1] = props.dimension;
286
287
288
289
290
291
292
293
294
295
    offsets[0] = offset;
    offsets[1] = 0;
  } else {
    rank = 1;
    shape[0] = N;
    shape[1] = 0;
    offsets[0] = offset;
    offsets[1] = 0;
  }

296
  /* Change shape of memory data space */
297
  hid_t h_err = H5Sset_extent_simple(h_memspace, rank, shape, NULL);
298
299
  if (h_err < 0) {
    error("Error while changing data space (memory) shape for field '%s'.",
300
          props.name);
301
  }
302

303
304
  /* Select the hyper-salb corresponding to this rank */
  hid_t h_filespace = H5Dget_space(h_data);
Matthieu Schaller's avatar
Matthieu Schaller committed
305
306
307
  if (N > 0) {
    H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape,
                        NULL);
308
309
310
  } else {
    H5Sselect_none(h_filespace);
  }
311

Matthieu Schaller's avatar
Matthieu Schaller committed
312
313
  /* message("Writing %lld '%s', %zd elements = %zd bytes (int=%d) at offset
   * %zd", */
314
  /* 	  N, props.name, N * props.dimension, N * props.dimension * typeSize, */
315
316
317
318
  /* 	  (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
319
                   h_plist_id, temp);
320
321
322
323
324
325
326
327
328
  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);
}
329

330
331
332
333
334
/**
 * @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.
335
336
 * @param fileName The name of the file in which the data is written.
 * @param xmfFile The FILE used to write the XMF description.
337
338
339
340
 * @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.
341
342
343
344
345
 * @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.
346
 */
347
void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
348
                char* partTypeGroupName, struct io_props props, size_t N,
349
350
351
352
353
354
355
356
357
                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];
358
  hsize_t chunk_shape[2];
359
360
361
362
  if (props.dimension > 1) {
    rank = 2;
    shape_total[0] = N_total;
    shape_total[1] = props.dimension;
363
364
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = props.dimension;
365
366
367
368
  } else {
    rank = 1;
    shape_total[0] = N_total;
    shape_total[1] = 0;
369
370
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = 0;
371
372
  }

373
  /* Make sure the chunks are not larger than the dataset */
lhausamm's avatar
lhausamm committed
374
  if (chunk_shape[0] > (hsize_t)N_total) chunk_shape[0] = N_total;
375
376

  /* Create the space in the file */
377
378
379
380
  hid_t h_filespace = H5Screate(H5S_SIMPLE);
  if (h_filespace < 0) {
    error("Error while creating data space (file) for field '%s'.", props.name);
  }
381

382
  /* Change shape of file data space */
383
  hid_t h_err = H5Sset_extent_simple(h_filespace, rank, shape_total, NULL);
384
  if (h_err < 0) {
385
386
    error("Error while changing data space (file) shape for field '%s'.",
          props.name);
387
388
  }

389
390
391
392
  /* Dataset properties */
  const hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);

  /* Set chunk size */
393
394
395
396
397
  /* 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); */
  /* } */
398

399
  /* Create dataset */
400
401
  const hid_t h_data = H5Dcreate(grp, props.name, io_hdf5_type(props.type),
                                 h_filespace, H5P_DEFAULT, h_prop, H5P_DEFAULT);
402
  if (h_data < 0) {
403
    error("Error while creating dataset '%s'.", props.name);
404
405
  }

406
407
408
  H5Sclose(h_filespace);

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

412
413
414
  /* 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;
415
416
417
418
419
420
  while (redo) {

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

421
    /* Write the first chunk */
422
    const size_t this_chunk = (N > max_chunk_size) ? max_chunk_size : N;
423
    writeArray_chunk(e, h_data, h_plist_id, props, this_chunk, offset,
424
                     internal_units, snapshot_units);
425
426

    /* Compute how many items are left */
427
    if (N > max_chunk_size) {
428
      N -= max_chunk_size;
429
      props.field += max_chunk_size * props.partSize;
430
431
      offset += max_chunk_size;
      redo = 1;
432
    } else {
433
434
435
436
437
438
      N = 0;
      offset += 0;
      redo = 0;
    }

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

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

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

451
  /* Write unit conversion factors for this data set */
452
453
  char buffer[FIELD_BUFFER_SIZE];
  units_cgs_conversion_string(buffer, snapshot_units, props.units);
454
455
456
457
458
459
460
461
  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);
462

463
464
  /* Close everything */
  H5Pclose(h_prop);
465
466
467
468
  H5Dclose(h_data);
  H5Pclose(h_plist_id);
}

469
470
471
472
/**
 * @brief Reads an HDF5 initial condition file (GADGET-3 type) in parallel
 *
 * @param fileName The file to read.
473
 * @param internal_units The system units used internally
474
475
 * @param dim (output) The dimension of the volume read from the file.
 * @param parts (output) The array of #part read from the file.
476
477
478
479
480
 * @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.
481
 * @param periodic (output) 1 if the volume is periodic, 0 if not.
482
 * @param flag_entropy (output) 1 if the ICs contained Entropy in the
483
 * InternalEnergy field
484
485
486
 * @param with_hydro Are we running with hydro ?
 * @param with_gravity Are we running with gravity ?
 * @param with_stars Are we running with stars ?
487
488
489
490
 * @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
491
 * @param dry_run If 1, don't read the particle. Only allocates the arrays.
492
493
 *
 */
494
void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
495
                      double dim[3], struct part** parts, struct gpart** gparts,
Matthieu Schaller's avatar
Matthieu Schaller committed
496
497
498
499
500
                      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) {
501

502
  hid_t h_file = 0, h_grp = 0;
503
  /* GADGET has only cubic boxes (in cosmological mode) */
Matthieu Schaller's avatar
Matthieu Schaller committed
504
  double boxSize[3] = {0.0, -1.0, -1.0};
505
506
507
508
509
  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};
510
  int dimension = 3; /* Assume 3D if nothing is specified */
511
  size_t Ndm = 0;
512
513
514
515
516
517
518
519
520
521
522
523

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

  /* Read the relevant information */
528
  io_read_attribute(h_grp, "PeriodicBoundariesOn", INT, periodic);
529
530
531
532
533
534

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

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

538
539
540
541
  /* 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");
542
  if (hid_dim > 0) io_read_attribute(h_grp, "Dimension", INT, &dimension);
543
544
545
546
  if (dimension != hydro_dimension)
    error("ICs dimensionality (%dD) does not match code dimensionality (%dD)",
          dimension, (int)hydro_dimension);

547
  /* Read the relevant information and print status */
548
  int flag_entropy_temp[6];
549
  io_read_attribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy_temp);
550
  *flag_entropy = flag_entropy_temp[0];
551
552
553
554
  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);
555

556
  for (int ptype = 0; ptype < swift_type_count; ++ptype)
Matthieu Schaller's avatar
Matthieu Schaller committed
557
558
    N_total[ptype] =
        (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);
559

560
  /* Get the box size if not cubic */
561
562
563
564
  dim[0] = boxSize[0];
  dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
  dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];

565
566
567
568
569
570
  /* 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];

571
572
573
  /* message("Found %lld particles in a %speriodic box of size [%f %f %f].", */
  /* 	  N_total[0], (periodic ? "": "non-"), dim[0], */
  /* 	  dim[1], dim[2]); */
574
575

  /* Divide the particles among the tasks. */
576
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
577
578
579
    offset[ptype] = mpi_rank * N_total[ptype] / mpi_size;
    N[ptype] = (mpi_rank + 1) * N_total[ptype] / mpi_size - offset[ptype];
  }
580
581
582
583

  /* Close header */
  H5Gclose(h_grp);

584
  /* Read the unit system used in the ICs */
585
  struct unit_system* ic_units = malloc(sizeof(struct unit_system));
586
  if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
587
  io_read_unit_system(h_file, ic_units, mpi_rank);
588
589
590
591
592

  /* 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
593
      message("IC and internal units match. No conversion needed.");
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619

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

620
  /* Convert the dimensions of the box */
Matthieu Schaller's avatar
Matthieu Schaller committed
621
622
623
  for (int j = 0; j < 3; j++)
    dim[j] *=
        units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
624

625
  /* Allocate memory to store SPH particles */
626
627
  if (with_hydro) {
    *Ngas = N[0];
Matthieu Schaller's avatar
Matthieu Schaller committed
628
629
    if (posix_memalign((void*)parts, part_align,
                       (*Ngas) * sizeof(struct part)) != 0)
630
631
632
633
634
635
      error("Error while allocating memory for particles");
    bzero(*parts, *Ngas * sizeof(struct part));
  }

  /* Allocate memory to store star particles */
  if (with_stars) {
636
    *Nstars = N[swift_type_star];
637
638
639
640
641
642
643
644
645
    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];
646
647
648
    *Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
               N[swift_type_dark_matter] +
               (with_stars ? N[swift_type_star] : 0);
649
    if (posix_memalign((void*)gparts, gpart_align,
Matthieu Schaller's avatar
Matthieu Schaller committed
650
                       *Ngparts * sizeof(struct gpart)) != 0)
651
652
653
654
655
      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) /
656
657
   * (1024.*1024.)); */

658
  /* message("BoxSize = %lf", dim[0]); */
659
  /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */
660
661

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

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

667
668
669
    /* 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
670
             ptype);
671
672
673
674
    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
675

676
677
    int num_fields = 0;
    struct io_props list[100];
678
    size_t Nparticles = 0;
679

680
681
    /* Read particle fields into the particle structure */
    switch (ptype) {
Matthieu Schaller's avatar
Matthieu Schaller committed
682

683
      case swift_type_gas:
Matthieu Schaller's avatar
Matthieu Schaller committed
684
685
686
687
        if (with_hydro) {
          Nparticles = *Ngas;
          hydro_read_particles(*parts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
688
689
        break;

690
      case swift_type_dark_matter:
Matthieu Schaller's avatar
Matthieu Schaller committed
691
692
693
694
        if (with_gravity) {
          Nparticles = Ndm;
          darkmatter_read_particles(*gparts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
695
696
        break;

697
      case swift_type_star:
698
        if (with_stars) {
Matthieu Schaller's avatar
Matthieu Schaller committed
699
700
701
702
          Nparticles = *Nstars;
          star_read_particles(*sparts, list, &num_fields);
        }
        break;
703

Matthieu Schaller's avatar
Matthieu Schaller committed
704
      default:
705
706
707
        if (mpi_rank == 0)
          message("Particle Type %d not yet supported. Particles ignored",
                  ptype);
708
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
709

710
711
712
    /* Read everything */
    if (!dry_run)
      for (int i = 0; i < num_fields; ++i)
713
714
        readArray(h_grp, list[i], Nparticles, N_total[ptype], mpi_rank,
                  offset[ptype], internal_units, ic_units);
715

716
717
718
    /* Close particle group */
    H5Gclose(h_grp);
  }
719

720
  /* Prepare the DM particles */
721
  if (!dry_run && with_gravity) io_prepare_dm_gparts(*gparts, Ndm);
722

723
724
  /* Duplicate the hydro particles into gparts */
  if (!dry_run && with_gravity && with_hydro)
725
    io_duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
726
727
728

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

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

733
734
735
  /* Clean up */
  free(ic_units);

736
737
738
739
740
741
742
  /* Close property handler */
  H5Pclose(h_plist_id);

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

743
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
744
745
 * @brief Writes an HDF5 output file (GADGET-3 type) with
 *its XMF descriptor
746
747
 *
 * @param e The engine containing all the system.
748
 * @param baseName The common part of the snapshot file name.
749
750
 * @param internal_units The #unit_system used internally
 * @param snapshot_units The #unit_system used in the snapshots
751
752
753
754
 * @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
755
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
756
 * Creates an HDF5 output file and writes the particles
757
758
 * contained in the engine. If such a file already exists, it is
 * erased and replaced by the new one.
759
760
761
762
763
 * The companion XMF file is also updated accordingly.
 *
 * Calls #error() if an error occurs.
 *
 */
764
void write_output_parallel(struct engine* e, const char* baseName,
765
766
                           const struct unit_system* internal_units,
                           const struct unit_system* snapshot_units,
767
768
769
                           int mpi_rank, int mpi_size, MPI_Comm comm,
                           MPI_Info info) {

770
  hid_t h_file = 0, h_grp = 0;
771
  const size_t Ngas = e->s->nr_parts;
772
  const size_t Nstars = e->s->nr_sparts;
773
  const size_t Ntot = e->s->nr_gparts;
774
775
776
  int periodic = e->s->periodic;
  int numFiles = 1;
  struct part* parts = e->s->parts;
777
778
  struct gpart* gparts = e->s->gparts;
  struct gpart* dmparts = NULL;
779
  struct spart* sparts = e->s->sparts;
780
  static int outputCount = 0;
781
782
  FILE* xmfFile = 0;

783
  /* Number of unassociated gparts */
784
  const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
785

786
  /* File name */
787
  char fileName[FILENAME_BUFFER_SIZE];
788
  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
789
           outputCount);
790
791

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

794
  /* Prepare the XMF file for the new entry */
795
  if (mpi_rank == 0) xmfFile = xmf_prepare_file(baseName);
796

797
  /* Open HDF5 file */
798
799
800
  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);
801
802
803
  if (h_file < 0) {
    error("Error while opening file '%s'.", fileName);
  }
804

Matthieu Schaller's avatar
Matthieu Schaller committed
805
806
  /* Compute offset in the file and total number of
   * particles */
807
808
809
810
811
  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)
812
813
    N_total[ptype] = offset[ptype] + N[ptype];

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

Matthieu Schaller's avatar
Matthieu Schaller committed
818
819
  /* Now everybody konws its offset and the total number of
   * particles of each
820
   * type */
821

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

826
827
  /* Open header to write simulation properties */
  /* message("Writing runtime parameters..."); */
Matthieu Schaller's avatar
Matthieu Schaller committed
828
829
  h_grp =
      H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
830
  if (h_grp < 0) error("Error while creating runtime parameters group\n");
831
832

  /* Write the relevant information */
833
  io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
834
835
836

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

838
839
  /* Open header to write simulation properties */
  /* message("Writing file header..."); */
840
  h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
841
  if (h_grp < 0) error("Error while creating file header\n");
Matthieu Schaller's avatar
Matthieu Schaller committed
842

843
  /* Print the relevant information and print status */
844
  io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
845
  double dblTime = e->time;
846
  io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
847
  int dimension = (int)hydro_dimension;
848
  io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
Matthieu Schaller's avatar
Matthieu Schaller committed
849

850
  /* GADGET-2 legacy values */
851
  /* Number of particles of each type */
852
853
854
  unsigned int numParticles[swift_type_count] = {0};
  unsigned int numParticlesHighWord[swift_type_count] = {0};
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
855
856
857
    numParticles[ptype] = (unsigned int)N_total[ptype];
    numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
  }
858
859
860
861
862
863
  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);
864
  double MassTable[6] = {0., 0., 0., 0., 0., 0.};
865
866
  io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
  unsigned int flagEntropy[swift_type_count] = {0};
867
  flagEntropy[0] = writeEntropyFlag();
868
869
870
  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
871

872
873
  /* Close header */
  H5Gclose(h_grp);
Matthieu Schaller's avatar
Matthieu Schaller committed
874

875
  /* Print the code version */
876
  io_write_code_description(h_file);
877

878
  /* Print the SPH parameters */
879
880
881
  if (e->policy & engine_policy_hydro) {
    h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
882
883
884
885
886
    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);
  }
887

888
  /* Print the gravity parameters */
889
  if (e->policy & engine_policy_self_gravity) {
890
    h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
891
                      H5P_DEFAULT);
892
893
894
895
    if (h_grp < 0) error("Error while creating gravity group");
    gravity_props_print_snapshot(h_grp, e->gravity_properties);
    H5Gclose(h_grp);
  }
896
897
898
899
900
901
902

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

904
  /* Print the system of Units used in the spashot */
905
  io_write_unit_system(h_file, snapshot_units, "Units");
906
907

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

910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
  /* 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);
    }
  }

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

946
947
    /* Don't do anything if no particle of this kind */
    if (N_total[ptype] == 0) continue;
Matthieu Schaller's avatar
Matthieu Schaller committed
948
949
950

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