parallel_io.c 34.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
59
60
/**
 * @brief Reads a data array from a given HDF5 group.
 *
 * @param grp The group from which to read.
61
62
63
64
65
66
 * @param prop The #io_props of the field to read.
 * @param N The number of particles on that rank.
 * @param N_total The total number of particles.
 * @param offset The offset in the array on disk for this rank.
 * @param internal_units The #unit_system used internally.
 * @param ic_units The #unit_system used in the ICs.
67
 */
68
69
void readArray(hid_t grp, const struct io_props prop, size_t N,
               long long N_total, long long offset,
70
71
               const struct unit_system* internal_units,
               const struct unit_system* ic_units) {
72

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

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

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

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

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

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

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

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

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

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

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

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

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

153
    if (io_is_double_precision(prop.type)) {
154
155
156
157
158
159
      double* temp_d = temp;
      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= factor;
    } else {
      float* temp_f = temp;
      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor;
    }
160
  }
161
162

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

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

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

/**
181
 * @brief Writes a chunk of data in an open HDF5 dataset
182
 *
183
 * @param e The #engine we are writing from.
184
 * @param h_data The HDF5 dataset to write to.
185
 * @param h_plist_id the parallel HDF5 properties.
186
 * @param props The #io_props of the field to read.
187
 * @param N The number of particles to write.
188
189
190
 * @param offset Offset in the array where this mpi task starts writing.
 * @param internal_units The #unit_system used internally.
 * @param snapshot_units The #unit_system used in the snapshots.
191
 */
192
193
194
195
196
void writeArray_chunk(struct engine* e, hid_t h_data, hid_t h_plist_id,
                      const struct io_props props, size_t N, long long offset,
                      const struct unit_system* internal_units,
                      const struct unit_system* snapshot_units) {

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

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

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

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

  /* Copy particle data to temporary buffer */
212
213
214
215
216
217
218
219
220
221
222
  if (props.convert_part == NULL &&
      props.convert_gpart == NULL) { /* No conversion */

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

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

    float* temp_f = temp;
    for (size_t i = 0; i < N; ++i)
223
      temp_f[i] = props.convert_part(e, &props.parts[i]);
224
225
226
227
228

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

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

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

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

239
    if (io_is_double_precision(props.type)) {
240
241
242
243
244
245
246
      double* temp_d = temp;
      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= factor;
    } else {
      float* temp_f = temp;
      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor;
    }
  }
247
248

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

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

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

279
280
  /* Select the hyper-salb corresponding to this rank */
  hid_t h_filespace = H5Dget_space(h_data);
Matthieu Schaller's avatar
Matthieu Schaller committed
281
282
283
  if (N > 0) {
    H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape,
                        NULL);
284
285
286
  } else {
    H5Sselect_none(h_filespace);
  }
287

Matthieu Schaller's avatar
Matthieu Schaller committed
288
289
  /* message("Writing %lld '%s', %zd elements = %zd bytes (int=%d) at offset
   * %zd", */
290
  /* 	  N, props.name, N * props.dimension, N * props.dimension * typeSize, */
291
292
293
294
  /* 	  (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
295
                   h_plist_id, temp);
296
297
298
299
300
301
302
303
304
  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);
}
305

306
307
308
309
310
/**
 * @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.
311
312
 * @param fileName The name of the file in which the data is written.
 * @param xmfFile The FILE used to write the XMF description.
313
314
315
316
 * @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.
317
318
319
320
321
 * @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.
322
 */
323
void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
324
                char* partTypeGroupName, struct io_props props, size_t N,
325
326
327
328
329
330
331
332
333
                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];
334
  hsize_t chunk_shape[2];
335
336
337
338
  if (props.dimension > 1) {
    rank = 2;
    shape_total[0] = N_total;
    shape_total[1] = props.dimension;
339
340
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = props.dimension;
341
342
343
344
  } else {
    rank = 1;
    shape_total[0] = N_total;
    shape_total[1] = 0;
345
346
    chunk_shape[0] = 1 << 16; /* Just a guess...*/
    chunk_shape[1] = 0;
347
348
  }

349
350
351
352
  /* Make sure the chunks are not larger than the dataset */
  if (chunk_shape[0] > N_total) chunk_shape[0] = N_total;

  /* Create the space in the file */
353
354
355
356
  hid_t h_filespace = H5Screate(H5S_SIMPLE);
  if (h_filespace < 0) {
    error("Error while creating data space (file) for field '%s'.", props.name);
  }
357

358
  /* Change shape of file data space */
359
  hid_t h_err = H5Sset_extent_simple(h_filespace, rank, shape_total, NULL);
360
  if (h_err < 0) {
361
362
    error("Error while changing data space (file) shape for field '%s'.",
          props.name);
363
364
  }

365
366
367
368
  /* Dataset properties */
  const hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);

  /* Set chunk size */
369
370
371
372
373
  /* 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); */
  /* } */
374

375
  /* Create dataset */
376
377
  const hid_t h_data = H5Dcreate(grp, props.name, io_hdf5_type(props.type),
                                 h_filespace, H5P_DEFAULT, h_prop, H5P_DEFAULT);
378
  if (h_data < 0) {
379
    error("Error while creating dataset '%s'.", props.name);
380
381
  }

382
383
384
  H5Sclose(h_filespace);

  /* Create property list for collective dataset write.    */
385
  const hid_t h_plist_id = H5Pcreate(H5P_DATASET_XFER);
386
387
  H5Pset_dxpl_mpio(h_plist_id, H5FD_MPIO_COLLECTIVE);

388
389
390
  /* 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;
391
392
393
394
395
396
  while (redo) {

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

397
    /* Write the first chunk */
398
    const size_t this_chunk = (N > max_chunk_size) ? max_chunk_size : N;
399
    writeArray_chunk(e, h_data, h_plist_id, props, this_chunk, offset,
400
                     internal_units, snapshot_units);
401
402

    /* Compute how many items are left */
403
    if (N > max_chunk_size) {
404
      N -= max_chunk_size;
405
      props.field += max_chunk_size * props.partSize;
406
407
      offset += max_chunk_size;
      redo = 1;
408
    } else {
409
410
411
412
413
414
      N = 0;
      offset += 0;
      redo = 0;
    }

    /* Do we need to run again ? */
415
416
    MPI_Allreduce(MPI_IN_PLACE, &redo, 1, MPI_SIGNED_CHAR, MPI_MAX,
                  MPI_COMM_WORLD);
417
418
419

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

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

427
  /* Write unit conversion factors for this data set */
428
429
  char buffer[FIELD_BUFFER_SIZE];
  units_cgs_conversion_string(buffer, snapshot_units, props.units);
430
431
432
433
434
435
436
437
  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);
438

439
440
  /* Close everything */
  H5Pclose(h_prop);
441
442
443
444
  H5Dclose(h_data);
  H5Pclose(h_plist_id);
}

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

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

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

  /* Read the relevant information */
504
  io_read_attribute(h_grp, "PeriodicBoundariesOn", INT, periodic);
505
506
507
508
509
510

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

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

514
515
516
517
  /* 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");
518
  if (hid_dim > 0) io_read_attribute(h_grp, "Dimension", INT, &dimension);
519
520
521
522
  if (dimension != hydro_dimension)
    error("ICs dimensionality (%dD) does not match code dimensionality (%dD)",
          dimension, (int)hydro_dimension);

523
  /* Read the relevant information and print status */
524
  int flag_entropy_temp[6];
525
  io_read_attribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy_temp);
526
  *flag_entropy = flag_entropy_temp[0];
527
528
529
530
  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);
531

532
  for (int ptype = 0; ptype < swift_type_count; ++ptype)
Matthieu Schaller's avatar
Matthieu Schaller committed
533
534
    N_total[ptype] =
        (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);
535

536
  /* Get the box size if not cubic */
537
538
539
540
  dim[0] = boxSize[0];
  dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
  dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];

541
542
543
544
545
546
  /* 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];

547
548
549
  /* message("Found %lld particles in a %speriodic box of size [%f %f %f].", */
  /* 	  N_total[0], (periodic ? "": "non-"), dim[0], */
  /* 	  dim[1], dim[2]); */
550
551

  /* Divide the particles among the tasks. */
552
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
553
554
555
    offset[ptype] = mpi_rank * N_total[ptype] / mpi_size;
    N[ptype] = (mpi_rank + 1) * N_total[ptype] / mpi_size - offset[ptype];
  }
556
557
558
559

  /* Close header */
  H5Gclose(h_grp);

560
  /* Read the unit system used in the ICs */
561
  struct unit_system* ic_units = malloc(sizeof(struct unit_system));
562
  if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
563
  io_read_unit_system(h_file, ic_units);
564
565
566
567
568

  /* 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
569
      message("IC and internal units match. No conversion needed.");
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595

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

596
  /* Convert the dimensions of the box */
Matthieu Schaller's avatar
Matthieu Schaller committed
597
598
599
  for (int j = 0; j < 3; j++)
    dim[j] *=
        units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
600

601
  /* Allocate memory to store SPH particles */
602
603
  if (with_hydro) {
    *Ngas = N[0];
Matthieu Schaller's avatar
Matthieu Schaller committed
604
605
    if (posix_memalign((void*)parts, part_align,
                       (*Ngas) * sizeof(struct part)) != 0)
606
607
608
609
610
611
      error("Error while allocating memory for particles");
    bzero(*parts, *Ngas * sizeof(struct part));
  }

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

634
  /* message("BoxSize = %lf", dim[0]); */
635
  /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */
636
637

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

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

643
644
645
    /* 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
646
             ptype);
647
648
649
650
    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
651

652
653
    int num_fields = 0;
    struct io_props list[100];
654
    size_t Nparticles = 0;
655

656
657
    /* Read particle fields into the particle structure */
    switch (ptype) {
Matthieu Schaller's avatar
Matthieu Schaller committed
658

659
      case swift_type_gas:
Matthieu Schaller's avatar
Matthieu Schaller committed
660
661
662
663
        if (with_hydro) {
          Nparticles = *Ngas;
          hydro_read_particles(*parts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
664
665
        break;

666
      case swift_type_dark_matter:
Matthieu Schaller's avatar
Matthieu Schaller committed
667
668
669
670
        if (with_gravity) {
          Nparticles = Ndm;
          darkmatter_read_particles(*gparts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
671
672
        break;

673
      case swift_type_star:
674
        if (with_stars) {
Matthieu Schaller's avatar
Matthieu Schaller committed
675
676
677
678
          Nparticles = *Nstars;
          star_read_particles(*sparts, list, &num_fields);
        }
        break;
679

Matthieu Schaller's avatar
Matthieu Schaller committed
680
      default:
681
        message("Particle Type %d not yet supported. Particles ignored", ptype);
682
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
683

684
685
686
    /* Read everything */
    if (!dry_run)
      for (int i = 0; i < num_fields; ++i)
687
        readArray(h_grp, list[i], Nparticles, N_total[ptype], offset[ptype],
688
689
                  internal_units, ic_units);

690
691
692
    /* Close particle group */
    H5Gclose(h_grp);
  }
693

694
  /* Prepare the DM particles */
695
  if (!dry_run && with_gravity) io_prepare_dm_gparts(*gparts, Ndm);
696

697
698
  /* Duplicate the hydro particles into gparts */
  if (!dry_run && with_gravity && with_hydro)
699
    io_duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
700
701
702

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

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

707
708
709
  /* Clean up */
  free(ic_units);

710
711
712
713
714
715
716
  /* Close property handler */
  H5Pclose(h_plist_id);

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

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

744
  hid_t h_file = 0, h_grp = 0;
745
  const size_t Ngas = e->s->nr_parts;
746
  const size_t Nstars = e->s->nr_sparts;
747
  const size_t Ntot = e->s->nr_gparts;
748
749
750
  int periodic = e->s->periodic;
  int numFiles = 1;
  struct part* parts = e->s->parts;
751
752
  struct gpart* gparts = e->s->gparts;
  struct gpart* dmparts = NULL;
753
  struct spart* sparts = e->s->sparts;
754
  static int outputCount = 0;
755
756
  FILE* xmfFile = 0;

757
  /* Number of unassociated gparts */
758
  const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
759

760
  /* File name */
761
  char fileName[FILENAME_BUFFER_SIZE];
762
  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
763
           outputCount);
764
765

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

768
  /* Prepare the XMF file for the new entry */
769
  if (mpi_rank == 0) xmfFile = xmf_prepare_file(baseName);
770

771
  /* Open HDF5 file */
772
773
774
  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);
775
776
777
  if (h_file < 0) {
    error("Error while opening file '%s'.", fileName);
  }
778

Matthieu Schaller's avatar
Matthieu Schaller committed
779
780
  /* Compute offset in the file and total number of
   * particles */
781
782
783
784
785
  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)
786
787
    N_total[ptype] = offset[ptype] + N[ptype];

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

Matthieu Schaller's avatar
Matthieu Schaller committed
792
793
  /* Now everybody konws its offset and the total number of
   * particles of each
794
   * type */
795

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

800
801
  /* Open header to write simulation properties */
  /* message("Writing runtime parameters..."); */
Matthieu Schaller's avatar
Matthieu Schaller committed
802
803
  h_grp =
      H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
804
  if (h_grp < 0) error("Error while creating runtime parameters group\n");
805
806

  /* Write the relevant information */
807
  io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
808
809
810

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

812
813
  /* Open header to write simulation properties */
  /* message("Writing file header..."); */
814
  h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
815
  if (h_grp < 0) error("Error while creating file header\n");
Matthieu Schaller's avatar
Matthieu Schaller committed
816

817
  /* Print the relevant information and print status */
818
  io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
819
  double dblTime = e->time;
820
  io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
821
  int dimension = (int)hydro_dimension;
822
  io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
Matthieu Schaller's avatar
Matthieu Schaller committed
823

824
  /* GADGET-2 legacy values */
825
  /* Number of particles of each type */
826
827
828
  unsigned int numParticles[swift_type_count] = {0};
  unsigned int numParticlesHighWord[swift_type_count] = {0};
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
829
830
831
    numParticles[ptype] = (unsigned int)N_total[ptype];
    numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
  }
832
833
834
835
836
837
  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);
838
  double MassTable[6] = {0., 0., 0., 0., 0., 0.};
839
840
  io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
  unsigned int flagEntropy[swift_type_count] = {0};
841
  flagEntropy[0] = writeEntropyFlag();
842
843
844
  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
845

846
847
  /* Close header */
  H5Gclose(h_grp);
Matthieu Schaller's avatar
Matthieu Schaller committed
848

849
  /* Print the code version */
850
  io_write_code_description(h_file);
851

852
  /* Print the SPH parameters */
853
854
855
  if (e->policy & engine_policy_hydro) {
    h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
856
857
858
859
860
    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);
  }
861

862
  /* Print the gravity parameters */
863
  if (e->policy & engine_policy_self_gravity) {
864
    h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
865
                      H5P_DEFAULT);
866
867
868
869
    if (h_grp < 0) error("Error while creating gravity group");
    gravity_props_print_snapshot(h_grp, e->gravity_properties);
    H5Gclose(h_grp);
  }
870
871
872
873
874
875
876

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

878
  /* Print the system of Units used in the spashot */
879
  io_write_unit_system(h_file, snapshot_units, "Units");
880
881

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

884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
  /* 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);
    }
  }

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