parallel_io.c 36.9 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
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
257
  const size_t copySize = typeSize * props.dimension;
  const size_t num_elements = N * props.dimension;

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

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

  /* Allocate temporary buffer */
265
  void* temp = malloc(num_elements * typeSize);
266
  if (temp == NULL) error("Unable to allocate memory for temporary buffer");
267
268

  /* Copy particle data to temporary buffer */
269
270
271
272
273
274
275
276
277
278
279
  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)
280
      temp_f[i] = props.convert_part(e, &props.parts[i]);
281
282
283
284
285

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

    float* temp_f = temp;
    for (size_t i = 0; i < N; ++i)
286
      temp_f[i] = props.convert_gpart(e, &props.gparts[i]);
287
  }
288
289
290
291
292
293

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

296
    if (io_is_double_precision(props.type)) {
297
298
299
300
301
302
303
      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;
    }
  }
304
305

  /* Create data space */
306
  const hid_t h_memspace = H5Screate(H5S_SIMPLE);
307
  if (h_memspace < 0) {
308
309
    error("Error while creating data space (memory) for field '%s'.",
          props.name);
310
  }
311

312
313
314
315
  int rank;
  hsize_t shape[2];
  hsize_t offsets[2];
  if (props.dimension > 1) {
316
317
    rank = 2;
    shape[0] = N;
318
    shape[1] = props.dimension;
319
320
321
322
323
324
325
326
327
328
    offsets[0] = offset;
    offsets[1] = 0;
  } else {
    rank = 1;
    shape[0] = N;
    shape[1] = 0;
    offsets[0] = offset;
    offsets[1] = 0;
  }

329
  /* Change shape of memory data space */
330
  hid_t h_err = H5Sset_extent_simple(h_memspace, rank, shape, NULL);
331
332
  if (h_err < 0) {
    error("Error while changing data space (memory) shape for field '%s'.",
333
          props.name);
334
  }
335

336
337
  /* Select the hyper-salb corresponding to this rank */
  hid_t h_filespace = H5Dget_space(h_data);
Matthieu Schaller's avatar
Matthieu Schaller committed
338
339
340
  if (N > 0) {
    H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape,
                        NULL);
341
342
343
  } else {
    H5Sselect_none(h_filespace);
  }
344

Matthieu Schaller's avatar
Matthieu Schaller committed
345
346
  /* message("Writing %lld '%s', %zd elements = %zd bytes (int=%d) at offset
   * %zd", */
347
  /* 	  N, props.name, N * props.dimension, N * props.dimension * typeSize, */
348
349
350
351
  /* 	  (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
352
                   h_plist_id, temp);
353
354
355
356
357
358
359
360
361
  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);
}
362

363
364
365
366
367
/**
 * @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.
368
369
 * @param fileName The name of the file in which the data is written.
 * @param xmfFile The FILE used to write the XMF description.
370
371
372
373
 * @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.
374
375
376
377
378
 * @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.
379
 */
380
void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
381
                char* partTypeGroupName, struct io_props props, size_t N,
382
383
384
385
386
387
388
389
390
                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];
391
  hsize_t chunk_shape[2];
392
393
394
395
  if (props.dimension > 1) {
    rank = 2;
    shape_total[0] = N_total;
    shape_total[1] = props.dimension;
396
397
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = props.dimension;
398
399
400
401
  } else {
    rank = 1;
    shape_total[0] = N_total;
    shape_total[1] = 0;
402
403
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = 0;
404
405
  }

406
  /* Make sure the chunks are not larger than the dataset */
lhausamm's avatar
lhausamm committed
407
  if (chunk_shape[0] > (hsize_t)N_total) chunk_shape[0] = N_total;
408
409

  /* Create the space in the file */
410
411
412
413
  hid_t h_filespace = H5Screate(H5S_SIMPLE);
  if (h_filespace < 0) {
    error("Error while creating data space (file) for field '%s'.", props.name);
  }
414

415
  /* Change shape of file data space */
416
  hid_t h_err = H5Sset_extent_simple(h_filespace, rank, shape_total, NULL);
417
  if (h_err < 0) {
418
419
    error("Error while changing data space (file) shape for field '%s'.",
          props.name);
420
421
  }

422
423
424
425
  /* Dataset properties */
  const hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);

  /* Set chunk size */
426
427
428
429
430
  /* 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); */
  /* } */
431

432
  /* Create dataset */
433
434
  const hid_t h_data = H5Dcreate(grp, props.name, io_hdf5_type(props.type),
                                 h_filespace, H5P_DEFAULT, h_prop, H5P_DEFAULT);
435
  if (h_data < 0) {
436
    error("Error while creating dataset '%s'.", props.name);
437
438
  }

439
440
441
  H5Sclose(h_filespace);

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

445
446
447
  /* 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;
448
449
450
451
452
453
  while (redo) {

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

454
    /* Write the first chunk */
455
    const size_t this_chunk = (N > max_chunk_size) ? max_chunk_size : N;
456
    writeArray_chunk(e, h_data, h_plist_id, props, this_chunk, offset,
457
                     internal_units, snapshot_units);
458
459

    /* Compute how many items are left */
460
    if (N > max_chunk_size) {
461
      N -= max_chunk_size;
462
      props.field += max_chunk_size * props.partSize;
463
464
      offset += max_chunk_size;
      redo = 1;
465
    } else {
466
467
468
469
470
471
      N = 0;
      offset += 0;
      redo = 0;
    }

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

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

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

484
  /* Write unit conversion factors for this data set */
485
486
  char buffer[FIELD_BUFFER_SIZE];
  units_cgs_conversion_string(buffer, snapshot_units, props.units);
487
488
489
490
491
492
493
494
  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);
495

496
497
  /* Close everything */
  H5Pclose(h_prop);
498
499
500
501
  H5Dclose(h_data);
  H5Pclose(h_plist_id);
}

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

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

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

  /* Read the relevant information */
561
  io_read_attribute(h_grp, "PeriodicBoundariesOn", INT, periodic);
562
563
564
565
566
567

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

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

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

580
  /* Read the relevant information and print status */
581
  int flag_entropy_temp[6];
582
  io_read_attribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy_temp);
583
  *flag_entropy = flag_entropy_temp[0];
584
585
586
587
  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);
588

589
  for (int ptype = 0; ptype < swift_type_count; ++ptype)
Matthieu Schaller's avatar
Matthieu Schaller committed
590
591
    N_total[ptype] =
        (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);
592

593
  /* Get the box size if not cubic */
594
595
596
597
  dim[0] = boxSize[0];
  dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
  dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];

598
599
600
601
602
603
  /* 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];

604
605
606
  /* message("Found %lld particles in a %speriodic box of size [%f %f %f].", */
  /* 	  N_total[0], (periodic ? "": "non-"), dim[0], */
  /* 	  dim[1], dim[2]); */
607
608

  /* Divide the particles among the tasks. */
609
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
610
611
612
    offset[ptype] = mpi_rank * N_total[ptype] / mpi_size;
    N[ptype] = (mpi_rank + 1) * N_total[ptype] / mpi_size - offset[ptype];
  }
613
614
615
616

  /* Close header */
  H5Gclose(h_grp);

617
  /* Read the unit system used in the ICs */
618
  struct unit_system* ic_units = malloc(sizeof(struct unit_system));
619
  if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
620
  io_read_unit_system(h_file, ic_units);
621
622
623
624
625

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

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

653
  /* Convert the dimensions of the box */
Matthieu Schaller's avatar
Matthieu Schaller committed
654
655
656
  for (int j = 0; j < 3; j++)
    dim[j] *=
        units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
657

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

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

691
  /* message("BoxSize = %lf", dim[0]); */
692
  /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */
693
694

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

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

700
701
702
    /* 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
703
             ptype);
704
705
706
707
    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
708

709
710
    int num_fields = 0;
    struct io_props list[100];
711
    size_t Nparticles = 0;
712

713
714
    /* Read particle fields into the particle structure */
    switch (ptype) {
Matthieu Schaller's avatar
Matthieu Schaller committed
715

716
      case swift_type_gas:
Matthieu Schaller's avatar
Matthieu Schaller committed
717
718
719
720
        if (with_hydro) {
          Nparticles = *Ngas;
          hydro_read_particles(*parts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
721
722
        break;

723
      case swift_type_dark_matter:
Matthieu Schaller's avatar
Matthieu Schaller committed
724
725
726
727
        if (with_gravity) {
          Nparticles = Ndm;
          darkmatter_read_particles(*gparts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
728
729
        break;

730
      case swift_type_star:
731
        if (with_stars) {
Matthieu Schaller's avatar
Matthieu Schaller committed
732
733
734
735
          Nparticles = *Nstars;
          star_read_particles(*sparts, list, &num_fields);
        }
        break;
736

Matthieu Schaller's avatar
Matthieu Schaller committed
737
      default:
738
        message("Particle Type %d not yet supported. Particles ignored", ptype);
739
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
740

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

747
748
749
    /* Close particle group */
    H5Gclose(h_grp);
  }
750

751
  /* Prepare the DM particles */
752
  if (!dry_run && with_gravity) io_prepare_dm_gparts(*gparts, Ndm);
753

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

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

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

764
765
766
  /* Clean up */
  free(ic_units);

767
768
769
770
771
772
773
  /* Close property handler */
  H5Pclose(h_plist_id);

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

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

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

814
  /* Number of unassociated gparts */
815
  const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
816

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

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

825
  /* Prepare the XMF file for the new entry */
826
  if (mpi_rank == 0) xmfFile = xmf_prepare_file(baseName);
827

828
  /* Open HDF5 file */
829
830
831
  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);
832
833
834
  if (h_file < 0) {
    error("Error while opening file '%s'.", fileName);
  }
835

Matthieu Schaller's avatar
Matthieu Schaller committed
836
837
  /* Compute offset in the file and total number of
   * particles */
838
839
840
841
842
  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)
843
844
    N_total[ptype] = offset[ptype] + N[ptype];

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

Matthieu Schaller's avatar
Matthieu Schaller committed
849
850
  /* Now everybody konws its offset and the total number of
   * particles of each
851
   * type */
852

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

857
858
  /* Open header to write simulation properties */
  /* message("Writing runtime parameters..."); */
Matthieu Schaller's avatar
Matthieu Schaller committed
859
860
  h_grp =
      H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
861
  if (h_grp < 0) error("Error while creating runtime parameters group\n");
862
863

  /* Write the relevant information */
864
  io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
865
866
867

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

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

874
  /* Print the relevant information and print status */
875
  io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
876
  double dblTime = e->time;
877
  io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
878
  int dimension = (int)hydro_dimension;
879
  io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
Matthieu Schaller's avatar
Matthieu Schaller committed
880

881
  /* GADGET-2 legacy values */
882
  /* Number of particles of each type */
883
884
885
  unsigned int numParticles[swift_type_count] = {0};
  unsigned int numParticlesHighWord[swift_type_count] = {0};
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
886
887
888
    numParticles[ptype] = (unsigned int)N_total[ptype];
    numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
  }
889
890
891
892
893
894
  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);
895
  double MassTable[6] = {0., 0., 0., 0., 0., 0.};
896
897
  io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
  unsigned int flagEntropy[swift_type_count] = {0};
898
  flagEntropy[0] = writeEntropyFlag();
899
900
901
  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
902

903
904
  /* Close header */
  H5Gclose(h_grp);
Matthieu Schaller's avatar
Matthieu Schaller committed
905

906
  /* Print the code version */
907
  io_write_code_description(h_file);
908

909
  /* Print the SPH parameters */
910
911
912
  if (e->policy & engine_policy_hydro) {
    h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
913
914
915
916
917
    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);
  }
918

919
  /* Print the gravity parameters */
920
  if (e->policy & engine_policy_self_gravity) {
921
    h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
922
                      H5P_DEFAULT);
923
924
925
926
    if (h_grp < 0) error("Error while creating gravity group");
    gravity_props_print_snapshot(h_grp, e->gravity_properties);
    H5Gclose(h_grp);
  }
927
928
929
930
931
932
933

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

935
  /* Print the system of Units used in the spashot */
936
  io_write_unit_system(h_file, snapshot_units, "Units");
937
938

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

941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
  /* Tell the user if a conversion will be needed */
  if (e->verbose && mpi_rank == 0) {
    if (units_are_equal(snapshot_units, inter