parallel_io.c 30.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
44
#include "gravity_io.h"
#include "hydro_io.h"
45
#include "hydro_properties.h"
46
#include "io_properties.h"
47
48
#include "kernel_hydro.h"
#include "part.h"
49
#include "stars_io.h"
50
#include "units.h"
51
52
53
54
55
56
57
58
59

/**
 * @brief Reads a data array from a given HDF5 group.
 *
 * @param grp The group from which to read.
 * @param name The name of the array to read.
 * @param type The #DATA_TYPE of the attribute.
 * @param N The number of particles.
 * @param dim The dimension of the data (1 for scalar, 3 for vector)
60
 * @param part_c A (char*) pointer on the first occurrence of the field of
61
62
63
 *interest in the parts array
 * @param importance If COMPULSORY, the data must be present in the IC file. If
 *OPTIONAL, the array will be zeroed when the data is not present.
64
 *
65
 * @todo A better version using HDF5 hyper-slabs to read the file directly into
66
 *the part array
67
 * will be written once the structures have been stabilized.
68
 *
69
70
 * Calls #error() if an error occurs.
 */
71
72
73
void readArray(hid_t grp, const struct io_props prop, size_t N,
               long long N_total, long long offset,
               const struct UnitSystem* internal_units,
74
75
76
77
78
               const struct UnitSystem* ic_units) {

  const size_t typeSize = sizeOfType(prop.type);
  const size_t copySize = typeSize * prop.dimension;
  const size_t num_elements = N * prop.dimension;
79
80

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

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

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

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

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

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

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

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

  /* Set collective reading properties */
137
  const hid_t h_plist_id = H5Pcreate(H5P_DATASET_XFER);
138
139
  H5Pset_dxpl_mpio(h_plist_id, H5FD_MPIO_COLLECTIVE);

140
141
142
  /* Read HDF5 dataspace in temporary buffer */
  /* Dirty version that happens to work for vectors but should be improved */
  /* Using HDF5 dataspaces would be better */
143
144
  const hid_t h_err = H5Dread(h_data, hdf5Type(prop.type), h_memspace,
                              h_filespace, h_plist_id, temp);
145
  if (h_err < 0) {
146
147
148
149
150
151
152
153
    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
154
    /* message("Converting ! factor=%e", factor); */
155
156
157
158
159
160
161
162

    if (isDoublePrecision(prop.type)) {
      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;
    }
163
  }
164
165

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

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

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

/**
 * @brief Writes a data array in given HDF5 group.
 *
186
 * @param e The #engine we are writing from.
187
188
189
190
191
192
 * @param grp The group in which to write.
 * @param fileName The name of the file in which the data is written
 * @param xmfFile The FILE used to write the XMF description
 * @param N The number of particles to write.
 * @param N_total Total number of particles across all cores
 * @param offset Offset in the array where this mpi task starts writing
193
194
 * @param internal_units The #UnitSystem used internally
 * @param snapshot_units The #UnitSystem used in the snapshots
195
 *
196
 * @todo A better version using HDF5 hyper-slabs to write the file directly from
197
 * the part array will be written once the structures have been stabilized.
198
199
 *
 */
200
void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
201
202
203
204
205
206
207
208
209
                char* partTypeGroupName, const struct io_props props, size_t N,
                long long N_total, int mpi_rank, long long offset,
                const struct UnitSystem* internal_units,
                const struct UnitSystem* snapshot_units) {

  const size_t typeSize = sizeOfType(props.type);
  const size_t copySize = typeSize * props.dimension;
  const size_t num_elements = N * props.dimension;

Matthieu Schaller's avatar
Matthieu Schaller committed
210
  /* message("Writing '%s' array...", props.name); */
211
212

  /* Allocate temporary buffer */
213
  void* temp = malloc(num_elements * sizeOfType(props.type));
214
  if (temp == NULL) error("Unable to allocate memory for temporary buffer");
215
216

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

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

    float* temp_f = temp;
    for (size_t i = 0; i < N; ++i)
234
      temp_f[i] = props.convert_gpart(e, &props.gparts[i]);
235
  }
236
237
238
239
240
241

  /* 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
242
    /* message("Converting ! factor=%e", factor); */
243
244
245
246
247
248
249
250
251

    if (isDoublePrecision(props.type)) {
      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;
    }
  }
252
253

  /* Create data space */
254
  const hid_t h_memspace = H5Screate(H5S_SIMPLE);
255
  if (h_memspace < 0) {
256
257
    error("Error while creating data space (memory) for field '%s'.",
          props.name);
258
  }
259

260
  hid_t h_filespace = H5Screate(H5S_SIMPLE);
261
  if (h_filespace < 0) {
262
    error("Error while creating data space (file) for field '%s'.", props.name);
263
264
  }

265
266
267
268
269
  int rank;
  hsize_t shape[2];
  hsize_t shape_total[2];
  hsize_t offsets[2];
  if (props.dimension > 1) {
270
271
    rank = 2;
    shape[0] = N;
272
    shape[1] = props.dimension;
273
    shape_total[0] = N_total;
274
    shape_total[1] = props.dimension;
275
276
277
278
279
280
281
282
283
284
285
286
    offsets[0] = offset;
    offsets[1] = 0;
  } else {
    rank = 1;
    shape[0] = N;
    shape[1] = 0;
    shape_total[0] = N_total;
    shape_total[1] = 0;
    offsets[0] = offset;
    offsets[1] = 0;
  }

287
  /* Change shape of memory data space */
288
  hid_t h_err = H5Sset_extent_simple(h_memspace, rank, shape, NULL);
289
290
  if (h_err < 0) {
    error("Error while changing data space (memory) shape for field '%s'.",
291
          props.name);
292
  }
293
294
295

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

301
  /* Create dataset */
302
303
304
  const hid_t h_data =
      H5Dcreate(grp, props.name, hdf5Type(props.type), h_filespace, H5P_DEFAULT,
                H5P_DEFAULT, H5P_DEFAULT);
305
  if (h_data < 0) {
306
    error("Error while creating dataset '%s'.", props.name);
307
308
  }

309
310
311
312
313
  H5Sclose(h_filespace);
  h_filespace = H5Dget_space(h_data);
  H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape, NULL);

  /* Create property list for collective dataset write.    */
314
  const hid_t h_plist_id = H5Pcreate(H5P_DATASET_XFER);
315
316
317
  H5Pset_dxpl_mpio(h_plist_id, H5FD_MPIO_COLLECTIVE);

  /* Write temporary buffer to HDF5 dataspace */
318
319
  h_err = H5Dwrite(h_data, hdf5Type(props.type), h_memspace, h_filespace,
                   h_plist_id, temp);
320
  if (h_err < 0) {
321
    error("Error while writing data array '%s'.", props.name);
322
  }
323
324

  /* Write XMF description for this data set */
Matthieu Schaller's avatar
Matthieu Schaller committed
325
  if (mpi_rank == 0)
326
327
    writeXMFline(xmfFile, fileName, partTypeGroupName, props.name, N_total,
                 props.dimension, props.type);
328

329
  /* Write unit conversion factors for this data set */
330
331
  char buffer[FIELD_BUFFER_SIZE];
  units_cgs_conversion_string(buffer, snapshot_units, props.units);
332
  writeAttribute_d(h_data, "CGS conversion factor",
333
334
335
336
337
                   units_cgs_conversion_factor(snapshot_units, props.units));
  writeAttribute_f(h_data, "h-scale exponent",
                   units_h_factor(snapshot_units, props.units));
  writeAttribute_f(h_data, "a-scale exponent",
                   units_a_factor(snapshot_units, props.units));
338
339
  writeAttribute_s(h_data, "Conversion factor", buffer);

340
341
342
343
344
  /* Free and close everything */
  free(temp);
  H5Dclose(h_data);
  H5Pclose(h_plist_id);
  H5Sclose(h_memspace);
345
  H5Sclose(h_filespace);
346
347
}

348
349
350
351
/**
 * @brief Reads an HDF5 initial condition file (GADGET-3 type) in parallel
 *
 * @param fileName The file to read.
352
 * @param internal_units The system units used internally
353
354
355
356
 * @param dim (output) The dimension of the volume read from the file.
 * @param parts (output) The array of #part read from the file.
 * @param N (output) The number of particles read from the file.
 * @param periodic (output) 1 if the volume is periodic, 0 if not.
357
 * @param flag_entropy (output) 1 if the ICs contained Entropy in the
358
359
360
361
362
 * InternalEnergy field
 * @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
363
 * @param dry_run If 1, don't read the particle. Only allocates the arrays.
364
365
366
367
368
369
370
371
372
373
374
 *
 * Opens the HDF5 file fileName and reads the particles contained
 * in the parts array. N is the returned number of particles found
 * in the file.
 *
 * @warning Can not read snapshot distributed over more than 1 file !!!
 * @todo Read snapshots distributed in more than one file.
 *
 * Calls #error() if an error occurs.
 *
 */
375
376
void read_ic_parallel(char* fileName, const struct UnitSystem* internal_units,
                      double dim[3], struct part** parts, struct gpart** gparts,
Matthieu Schaller's avatar
Matthieu Schaller committed
377
378
379
380
381
                      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) {
382

383
  hid_t h_file = 0, h_grp = 0;
384
  /* GADGET has only cubic boxes (in cosmological mode) */
Matthieu Schaller's avatar
Matthieu Schaller committed
385
  double boxSize[3] = {0.0, -1.0, -1.0};
386
387
  long long numParticles[NUM_PARTICLE_TYPES] = {0};
  long long numParticles_highWord[NUM_PARTICLE_TYPES] = {0};
388
389
390
  size_t N[NUM_PARTICLE_TYPES] = {0};
  long long N_total[NUM_PARTICLE_TYPES] = {0};
  long long offset[NUM_PARTICLE_TYPES] = {0};
391
  int dimension = 3; /* Assume 3D if nothing is specified */
392
  size_t Ndm = 0;
393
394
395
396
397
398
399
400
401
402
403
404

  /* 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..."); */
405
  h_grp = H5Gopen(h_file, "/RuntimePars", H5P_DEFAULT);
406
407
408
409
410
411
412
413
414
415
  if (h_grp < 0) error("Error while opening runtime parameters\n");

  /* Read the relevant information */
  readAttribute(h_grp, "PeriodicBoundariesOn", INT, periodic);

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

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

419
420
421
422
423
424
425
426
427
  /* 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");
  if (hid_dim > 0) readAttribute(h_grp, "Dimension", INT, &dimension);
  if (dimension != hydro_dimension)
    error("ICs dimensionality (%dD) does not match code dimensionality (%dD)",
          dimension, (int)hydro_dimension);

428
  /* Read the relevant information and print status */
429
430
431
  int flag_entropy_temp[6];
  readAttribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy_temp);
  *flag_entropy = flag_entropy_temp[0];
432
  readAttribute(h_grp, "BoxSize", DOUBLE, boxSize);
433
434
  readAttribute(h_grp, "NumPart_Total", LONGLONG, numParticles);
  readAttribute(h_grp, "NumPart_Total_HighWord", LONGLONG, numParticles_highWord);
435

436
  for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype)
437
438
    N_total[ptype] = (numParticles[ptype]) +
                     (numParticles_highWord[ptype] << 32);
439

440
441
442
443
  dim[0] = boxSize[0];
  dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
  dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];

444
445
446
  /* message("Found %lld particles in a %speriodic box of size [%f %f %f].", */
  /* 	  N_total[0], (periodic ? "": "non-"), dim[0], */
  /* 	  dim[1], dim[2]); */
447
448

  /* Divide the particles among the tasks. */
449
450
451
452
  for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) {
    offset[ptype] = mpi_rank * N_total[ptype] / mpi_size;
    N[ptype] = (mpi_rank + 1) * N_total[ptype] / mpi_size - offset[ptype];
  }
453
454
455
456

  /* Close header */
  H5Gclose(h_grp);

457
458
459
460
461
462
463
464
465
  /* Read the unit system used in the ICs */
  struct UnitSystem* ic_units = malloc(sizeof(struct UnitSystem));
  if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
  readUnitSystem(h_file, ic_units);

  /* 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
466
      message("IC and internal units match. No conversion needed.");
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492

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

493
  /* Convert the dimensions of the box */
Matthieu Schaller's avatar
Matthieu Schaller committed
494
495
496
  for (int j = 0; j < 3; j++)
    dim[j] *=
        units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
497

498
  /* Allocate memory to store SPH particles */
499
500
  if (with_hydro) {
    *Ngas = N[0];
Matthieu Schaller's avatar
Matthieu Schaller committed
501
502
    if (posix_memalign((void*)parts, part_align,
                       (*Ngas) * sizeof(struct part)) != 0)
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
      error("Error while allocating memory for particles");
    bzero(*parts, *Ngas * sizeof(struct part));
  }

  /* Allocate memory to store star particles */
  if (with_stars) {
    *Nstars = N[STAR];
    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];
    *Ngparts = (with_hydro ? N[GAS] : 0) + N[DM] + (with_stars ? N[STAR] : 0);
    if (posix_memalign((void*)gparts, gpart_align,
Matthieu Schaller's avatar
Matthieu Schaller committed
521
                       *Ngparts * sizeof(struct gpart)) != 0)
522
523
524
525
526
      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) /
527
528
   * (1024.*1024.)); */

529
  /* message("BoxSize = %lf", dim[0]); */
530
  /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */
531
532
533

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

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

538
539
540
    /* 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
541
             ptype);
542
543
544
545
    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
546

547
548
    int num_fields = 0;
    struct io_props list[100];
549
    size_t Nparticles = 0;
550

551
552
    /* Read particle fields into the particle structure */
    switch (ptype) {
Matthieu Schaller's avatar
Matthieu Schaller committed
553
554

      case GAS:
Matthieu Schaller's avatar
Matthieu Schaller committed
555
556
557
558
        if (with_hydro) {
          Nparticles = *Ngas;
          hydro_read_particles(*parts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
559
560
561
        break;

      case DM:
Matthieu Schaller's avatar
Matthieu Schaller committed
562
563
564
565
        if (with_gravity) {
          Nparticles = Ndm;
          darkmatter_read_particles(*gparts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
566
567
        break;

568
569
      case STAR:
        if (with_stars) {
Matthieu Schaller's avatar
Matthieu Schaller committed
570
571
572
573
          Nparticles = *Nstars;
          star_read_particles(*sparts, list, &num_fields);
        }
        break;
574

Matthieu Schaller's avatar
Matthieu Schaller committed
575
      default:
576
        message("Particle Type %d not yet supported. Particles ignored", ptype);
577
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
578

579
580
581
    /* Read everything */
    if (!dry_run)
      for (int i = 0; i < num_fields; ++i)
582
        readArray(h_grp, list[i], Nparticles, N_total[ptype], offset[ptype],
583
584
                  internal_units, ic_units);

585
586
587
    /* Close particle group */
    H5Gclose(h_grp);
  }
588

589
  /* Prepare the DM particles */
590
  if (!dry_run && with_gravity) prepare_dm_gparts(*gparts, Ndm);
591

592
593
594
595
596
597
598
  /* Duplicate the hydro particles into gparts */
  if (!dry_run && with_gravity && with_hydro)
    duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);

  /* Duplicate the star particles into gparts */
  if (!dry_run && with_gravity && with_stars)
    duplicate_star_gparts(*sparts, *gparts, *Nstars, Ndm + *Ngas);
599
600

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

602
603
604
  /* Clean up */
  free(ic_units);

605
606
607
608
609
610
611
  /* Close property handler */
  H5Pclose(h_plist_id);

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

612
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
613
614
 * @brief Writes an HDF5 output file (GADGET-3 type) with
 *its XMF descriptor
615
616
 *
 * @param e The engine containing all the system.
617
 * @param baseName The common part of the snapshot file name.
618
619
 * @param internal_units The #UnitSystem used internally
 * @param snapshot_units The #UnitSystem used in the snapshots
620
621
622
623
 * @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
624
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
625
 * Creates an HDF5 output file and writes the particles
626
627
 * contained in the engine. If such a file already exists, it is
 * erased and replaced by the new one.
628
629
630
631
632
 * The companion XMF file is also updated accordingly.
 *
 * Calls #error() if an error occurs.
 *
 */
633
void write_output_parallel(struct engine* e, const char* baseName,
634
635
636
637
638
                           const struct UnitSystem* internal_units,
                           const struct UnitSystem* snapshot_units,
                           int mpi_rank, int mpi_size, MPI_Comm comm,
                           MPI_Info info) {

639
  hid_t h_file = 0, h_grp = 0;
640
  const size_t Ngas = e->s->nr_parts;
641
  const size_t Nstars = e->s->nr_sparts;
642
  const size_t Ntot = e->s->nr_gparts;
643
644
645
  int periodic = e->s->periodic;
  int numFiles = 1;
  struct part* parts = e->s->parts;
646
647
  struct gpart* gparts = e->s->gparts;
  struct gpart* dmparts = NULL;
648
  struct spart* sparts = e->s->sparts;
649
  static int outputCount = 0;
650
651
  FILE* xmfFile = 0;

652
  /* Number of unassociated gparts */
653
  const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
654

655
  /* File name */
656
  char fileName[FILENAME_BUFFER_SIZE];
657
658
  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%03i.hdf5", baseName,
           outputCount);
659
660

  /* First time, we need to create the XMF file */
661
  if (outputCount == 0 && mpi_rank == 0) createXMFfile(baseName);
662

663
  /* Prepare the XMF file for the new entry */
664
  if (mpi_rank == 0) xmfFile = prepareXMFfile(baseName);
665

666
  /* Open HDF5 file */
667
668
669
  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);
670
671
672
  if (h_file < 0) {
    error("Error while opening file '%s'.", fileName);
  }
673

Matthieu Schaller's avatar
Matthieu Schaller committed
674
675
  /* Compute offset in the file and total number of
   * particles */
676
  size_t N[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 0, 0, Nstars, 0};
677
678
  long long N_total[NUM_PARTICLE_TYPES] = {0};
  long long offset[NUM_PARTICLE_TYPES] = {0};
679
  MPI_Exscan(&N, &offset, NUM_PARTICLE_TYPES, MPI_LONG_LONG_INT, MPI_SUM, comm);
680
681
682
  for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype)
    N_total[ptype] = offset[ptype] + N[ptype];

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

Matthieu Schaller's avatar
Matthieu Schaller committed
687
688
  /* Now everybody konws its offset and the total number of
   * particles of each
689
   * type */
690

Matthieu Schaller's avatar
Matthieu Schaller committed
691
692
  /* Write the part of the XMF file corresponding to this
   * specific output */
693
  if (mpi_rank == 0) writeXMFoutputheader(xmfFile, fileName, e->time);
694

695
696
  /* Open header to write simulation properties */
  /* message("Writing runtime parameters..."); */
Matthieu Schaller's avatar
Matthieu Schaller committed
697
698
  h_grp =
      H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
699
  if (h_grp < 0) error("Error while creating runtime parameters group\n");
700
701
702
703
704
705

  /* Write the relevant information */
  writeAttribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);

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

707
708
  /* Open header to write simulation properties */
  /* message("Writing file header..."); */
709
  h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
710
  if (h_grp < 0) error("Error while creating file header\n");
Matthieu Schaller's avatar
Matthieu Schaller committed
711

712
  /* Print the relevant information and print status */
Pedro Gonnet's avatar
Pedro Gonnet committed
713
  writeAttribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
714
715
  double dblTime = e->time;
  writeAttribute(h_grp, "Time", DOUBLE, &dblTime, 1);
716
717
  int dimension = (int)hydro_dimension;
  writeAttribute(h_grp, "Dimension", INT, &dimension, 1);
Matthieu Schaller's avatar
Matthieu Schaller committed
718

719
  /* GADGET-2 legacy values */
720
721
722
723
724
725
726
727
  /* Number of particles of each type */
  unsigned int numParticles[NUM_PARTICLE_TYPES] = {0};
  unsigned int numParticlesHighWord[NUM_PARTICLE_TYPES] = {0};
  for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) {
    numParticles[ptype] = (unsigned int)N_total[ptype];
    numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
  }
  writeAttribute(h_grp, "NumPart_ThisFile", LONGLONG, N_total,
Matthieu Schaller's avatar
Matthieu Schaller committed
728
                 NUM_PARTICLE_TYPES);
729
  writeAttribute(h_grp, "NumPart_Total", UINT, numParticles,
Matthieu Schaller's avatar
Matthieu Schaller committed
730
                 NUM_PARTICLE_TYPES);
731
  writeAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticlesHighWord,
Matthieu Schaller's avatar
Matthieu Schaller committed
732
                 NUM_PARTICLE_TYPES);
733
  double MassTable[6] = {0., 0., 0., 0., 0., 0.};
734
735
  writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, NUM_PARTICLE_TYPES);
  unsigned int flagEntropy[NUM_PARTICLE_TYPES] = {0};
736
  flagEntropy[0] = writeEntropyFlag();
737
  writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
Matthieu Schaller's avatar
Matthieu Schaller committed
738
                 NUM_PARTICLE_TYPES);
739
  writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
Matthieu Schaller's avatar
Matthieu Schaller committed
740

741
742
  /* Close header */
  H5Gclose(h_grp);
Matthieu Schaller's avatar
Matthieu Schaller committed
743

744
745
746
  /* Print the code version */
  writeCodeDescription(h_file);

747
  /* Print the SPH parameters */
748
749
  h_grp =
      H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
750
  if (h_grp < 0) error("Error while creating SPH group");
751
  hydro_props_print_snapshot(h_grp, e->hydro_properties);
752
753
754
755
756
757
758
759
760
  writeSPHflavour(h_grp);
  H5Gclose(h_grp);

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

762
763
764
765
  /* Print the system of Units used in the spashot */
  writeUnitSystem(h_file, snapshot_units, "Units");

  /* Print the system of Units used internally */
766
  writeUnitSystem(h_file, internal_units, "InternalCodeUnits");
767

768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
  /* 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);
    }
  }

801
802
  /* Loop over all particle types */
  for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
803

804
805
    /* Don't do anything if no particle of this kind */
    if (N_total[ptype] == 0) continue;
Matthieu Schaller's avatar
Matthieu Schaller committed
806
807
808

    /* Add the global information for that particle type to
     * the XMF meta-file */
809
    if (mpi_rank == 0)
810
811
      writeXMFgroupheader(xmfFile, fileName, N_total[ptype],
                          (enum PARTICLE_TYPE)ptype);
Matthieu Schaller's avatar
Matthieu Schaller committed
812

813
814
815
    /* 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
816
817
818
             ptype);
    h_grp = H5Gcreate(h_file, partTypeGroupName, H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
819
820
821
    if (h_grp < 0) {
      error("Error while opening particle group %s.", partTypeGroupName);
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
822

823
824
    int num_fields = 0;
    struct io_props list[100];
825
    size_t Nparticles = 0;
826
827

    /* Write particle fields from the particle structure */
828
    switch (ptype) {
Matthieu Schaller's avatar
Matthieu Schaller committed
829
830

      case GAS:
831
        Nparticles = Ngas;
832
        hydro_write_particles(parts, list, &num_fields);
Matthieu Schaller's avatar
Matthieu Schaller committed
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
        break;

      case DM:
        /* Allocate temporary array */
        if (posix_memalign((void*)&dmparts, gpart_align,
                           Ndm * sizeof(struct gpart)) != 0)
          error(
              "Error while allocating temporart memory for "
              "DM particles");
        bzero(dmparts, Ndm * sizeof(struct gpart));

        /* Collect the DM particles from gpart */
        collect_dm_gparts(gparts, Ntot, dmparts, Ndm);

        /* Write DM particles */
848
        Nparticles = Ndm;
849
        darkmatter_write_particles(dmparts, list, &num_fields);
Matthieu Schaller's avatar
Matthieu Schaller committed
850
851
        break;

852
      case STAR:
Matthieu Schaller's avatar
Matthieu Schaller committed
853
854
855
        Nparticles = Nstars;
        star_write_particles(sparts, list, &num_fields);
        break;
856

Matthieu Schaller's avatar
Matthieu Schaller committed
857
858
      default:
        error("Particle Type %d not yet supported. Aborting", ptype);
859
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
860

861
862
    /* Write everything */
    for (int i = 0; i < num_fields; ++i)
863
864
865
      writeArray(e, h_grp, fileName, xmfFile, partTypeGroupName, list[i],
                 Nparticles, N_total[ptype], mpi_rank, offset[ptype],
                 internal_units, snapshot_units);
Matthieu Schaller's avatar
Matthieu Schaller committed
866

867
    /* Free temporary array */
868
869
870
871
    if (dmparts) {
      free(dmparts);
      dmparts = 0;
    }
872

873
874
    /* Close particle group */
    H5Gclose(h_grp);
Matthieu Schaller's avatar
Matthieu Schaller committed
875

876
    /* Close this particle group in the XMF file as well */
877
    if (mpi_rank == 0) writeXMFgroupfooter(xmfFile, (enum PARTICLE_TYPE)ptype);
878
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
879

880
  /* Write LXMF file descriptor */
881
  if (mpi_rank == 0) writeXMFoutputfooter(xmfFile, outputCount, e->time);
Matthieu Schaller's avatar
Matthieu Schaller committed
882

883
  /* message("Done writing particles..."); */
Matthieu Schaller's avatar
Matthieu Schaller committed
884

885
886
  /* Close property descriptor */
  H5Pclose(plist_id);
Matthieu Schaller's avatar
Matthieu Schaller committed
887

888
889
  /* Close file */
  H5Fclose(h_file);
Matthieu Schaller's avatar
Matthieu Schaller committed
890

891
892
893
  ++outputCount;
}

894
#endif /* HAVE_HDF5 */