parallel_io.c 31.1 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
#include "xmf.h"
52
53
54
55
56
57
58
59
60

/**
 * @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)
61
 * @param part_c A (char*) pointer on the first occurrence of the field of
62
63
64
 *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.
65
 *
66
 * @todo A better version using HDF5 hyper-slabs to read the file directly into
67
 *the part array
68
 * will be written once the structures have been stabilized.
69
 *
70
71
 * Calls #error() if an error occurs.
 */
72
73
74
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,
75
76
               const struct UnitSystem* ic_units) {

77
  const size_t typeSize = io_sizeof_type(prop.type);
78
79
  const size_t copySize = typeSize * prop.dimension;
  const size_t num_elements = N * prop.dimension;
80
81

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

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

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

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

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

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

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

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

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

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

157
    if (io_is_double_precision(prop.type)) {
158
159
160
161
162
163
      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;
    }
164
  }
165
166

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

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

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

/**
 * @brief Writes a data array in given HDF5 group.
 *
187
 * @param e The #engine we are writing from.
188
189
190
191
192
193
 * @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
194
195
 * @param internal_units The #UnitSystem used internally
 * @param snapshot_units The #UnitSystem used in the snapshots
196
 *
197
 * @todo A better version using HDF5 hyper-slabs to write the file directly from
198
 * the part array will be written once the structures have been stabilized.
199
200
 *
 */
201
void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
202
203
204
205
206
                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) {

207
  const size_t typeSize = io_sizeof_type(props.type);
208
209
210
  const size_t copySize = typeSize * props.dimension;
  const size_t num_elements = N * props.dimension;

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

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

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

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

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

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

245
    if (io_is_double_precision(props.type)) {
246
247
248
249
250
251
252
      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;
    }
  }
253
254

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

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

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

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

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

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

310
311
312
313
314
  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.    */
315
  const hid_t h_plist_id = H5Pcreate(H5P_DATASET_XFER);
316
317
318
  H5Pset_dxpl_mpio(h_plist_id, H5FD_MPIO_COLLECTIVE);

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

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

330
  /* Write unit conversion factors for this data set */
331
332
  char buffer[FIELD_BUFFER_SIZE];
  units_cgs_conversion_string(buffer, snapshot_units, props.units);
333
334
335
336
337
338
339
340
  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);
341

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

350
351
352
353
/**
 * @brief Reads an HDF5 initial condition file (GADGET-3 type) in parallel
 *
 * @param fileName The file to read.
354
 * @param internal_units The system units used internally
355
356
357
358
 * @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.
359
 * @param flag_entropy (output) 1 if the ICs contained Entropy in the
360
361
362
363
364
 * 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
365
 * @param dry_run If 1, don't read the particle. Only allocates the arrays.
366
367
368
369
370
371
372
373
374
375
376
 *
 * 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.
 *
 */
377
378
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
379
380
381
382
383
                      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) {
384

385
  hid_t h_file = 0, h_grp = 0;
386
  /* GADGET has only cubic boxes (in cosmological mode) */
Matthieu Schaller's avatar
Matthieu Schaller committed
387
  double boxSize[3] = {0.0, -1.0, -1.0};
388
389
390
391
392
  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};
393
  int dimension = 3; /* Assume 3D if nothing is specified */
394
  size_t Ndm = 0;
395
396
397
398
399
400
401
402
403
404
405
406

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

  /* Read the relevant information */
411
  io_read_attribute(h_grp, "PeriodicBoundariesOn", INT, periodic);
412
413
414
415
416
417

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

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

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

430
  /* Read the relevant information and print status */
431
  int flag_entropy_temp[6];
432
  io_read_attribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy_temp);
433
  *flag_entropy = flag_entropy_temp[0];
434
435
436
437
  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);
438

439
  for (int ptype = 0; ptype < swift_type_count; ++ptype)
Matthieu Schaller's avatar
Matthieu Schaller committed
440
441
    N_total[ptype] =
        (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);
442

443
444
445
446
  dim[0] = boxSize[0];
  dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
  dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];

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

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

  /* Close header */
  H5Gclose(h_grp);

460
461
462
  /* 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");
463
  io_read_UnitSystem(h_file, ic_units);
464
465
466
467
468

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

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

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

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

  /* Allocate memory to store star particles */
  if (with_stars) {
512
    *Nstars = N[swift_type_star];
513
514
515
516
517
518
519
520
521
    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];
522
523
524
    *Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
               N[swift_type_dark_matter] +
               (with_stars ? N[swift_type_star] : 0);
525
    if (posix_memalign((void*)gparts, gpart_align,
Matthieu Schaller's avatar
Matthieu Schaller committed
526
                       *Ngparts * sizeof(struct gpart)) != 0)
527
528
529
530
531
      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) /
532
533
   * (1024.*1024.)); */

534
  /* message("BoxSize = %lf", dim[0]); */
535
  /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */
536
537

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

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

543
544
545
    /* 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
546
             ptype);
547
548
549
550
    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
551

552
553
    int num_fields = 0;
    struct io_props list[100];
554
    size_t Nparticles = 0;
555

556
557
    /* Read particle fields into the particle structure */
    switch (ptype) {
Matthieu Schaller's avatar
Matthieu Schaller committed
558

559
      case swift_type_gas:
Matthieu Schaller's avatar
Matthieu Schaller committed
560
561
562
563
        if (with_hydro) {
          Nparticles = *Ngas;
          hydro_read_particles(*parts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
564
565
        break;

566
      case swift_type_dark_matter:
Matthieu Schaller's avatar
Matthieu Schaller committed
567
568
569
570
        if (with_gravity) {
          Nparticles = Ndm;
          darkmatter_read_particles(*gparts, list, &num_fields);
        }
Matthieu Schaller's avatar
Matthieu Schaller committed
571
572
        break;

573
      case swift_type_star:
574
        if (with_stars) {
Matthieu Schaller's avatar
Matthieu Schaller committed
575
576
577
578
          Nparticles = *Nstars;
          star_read_particles(*sparts, list, &num_fields);
        }
        break;
579

Matthieu Schaller's avatar
Matthieu Schaller committed
580
      default:
581
        message("Particle Type %d not yet supported. Particles ignored", ptype);
582
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
583

584
585
586
    /* Read everything */
    if (!dry_run)
      for (int i = 0; i < num_fields; ++i)
587
        readArray(h_grp, list[i], Nparticles, N_total[ptype], offset[ptype],
588
589
                  internal_units, ic_units);

590
591
592
    /* Close particle group */
    H5Gclose(h_grp);
  }
593

594
  /* Prepare the DM particles */
595
  if (!dry_run && with_gravity) io_prepare_dm_gparts(*gparts, Ndm);
596

597
598
  /* Duplicate the hydro particles into gparts */
  if (!dry_run && with_gravity && with_hydro)
599
    io_duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
600
601
602

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

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

607
608
609
  /* Clean up */
  free(ic_units);

610
611
612
613
614
615
616
  /* Close property handler */
  H5Pclose(h_plist_id);

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

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

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

657
  /* Number of unassociated gparts */
658
  const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
659

660
  /* File name */
661
  char fileName[FILENAME_BUFFER_SIZE];
662
663
  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%03i.hdf5", baseName,
           outputCount);
664
665

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

668
  /* Prepare the XMF file for the new entry */
669
  if (mpi_rank == 0) xmfFile = xmf_prepare_file(baseName);
670

671
  /* Open HDF5 file */
672
673
674
  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);
675
676
677
  if (h_file < 0) {
    error("Error while opening file '%s'.", fileName);
  }
678

Matthieu Schaller's avatar
Matthieu Schaller committed
679
680
  /* Compute offset in the file and total number of
   * particles */
681
682
683
684
685
  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)
686
687
    N_total[ptype] = offset[ptype] + N[ptype];

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

Matthieu Schaller's avatar
Matthieu Schaller committed
692
693
  /* Now everybody konws its offset and the total number of
   * particles of each
694
   * type */
695

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

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

  /* Write the relevant information */
707
  io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
708
709
710

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

712
713
  /* Open header to write simulation properties */
  /* message("Writing file header..."); */
714
  h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
715
  if (h_grp < 0) error("Error while creating file header\n");
Matthieu Schaller's avatar
Matthieu Schaller committed
716

717
  /* Print the relevant information and print status */
718
  io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
719
  double dblTime = e->time;
720
  io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
721
  int dimension = (int)hydro_dimension;
722
  io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
Matthieu Schaller's avatar
Matthieu Schaller committed
723

724
  /* GADGET-2 legacy values */
725
  /* Number of particles of each type */
726
727
728
  unsigned int numParticles[swift_type_count] = {0};
  unsigned int numParticlesHighWord[swift_type_count] = {0};
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
729
730
731
    numParticles[ptype] = (unsigned int)N_total[ptype];
    numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
  }
732
733
734
735
736
737
  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);
738
  double MassTable[6] = {0., 0., 0., 0., 0., 0.};
739
740
  io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
  unsigned int flagEntropy[swift_type_count] = {0};
741
  flagEntropy[0] = writeEntropyFlag();
742
743
744
  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
745

746
747
  /* Close header */
  H5Gclose(h_grp);
Matthieu Schaller's avatar
Matthieu Schaller committed
748

749
  /* Print the code version */
750
  io_write_code_description(h_file);
751

752
  /* Print the SPH parameters */
753
754
  h_grp =
      H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
755
  if (h_grp < 0) error("Error while creating SPH group");
756
  hydro_props_print_snapshot(h_grp, e->hydro_properties);
757
758
759
760
761
762
763
764
765
  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);
766

767
  /* Print the system of Units used in the spashot */
768
  io_write_UnitSystem(h_file, snapshot_units, "Units");
769
770

  /* Print the system of Units used internally */
771
  io_write_UnitSystem(h_file, internal_units, "InternalCodeUnits");
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
801
802
803
804
805
  /* 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);
    }
  }

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

809
810
    /* Don't do anything if no particle of this kind */
    if (N_total[ptype] == 0) continue;
Matthieu Schaller's avatar
Matthieu Schaller committed
811
812
813

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

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

828
829
    int num_fields = 0;
    struct io_props list[100];
830
    size_t Nparticles = 0;
831
832

    /* Write particle fields from the particle structure */
833
    switch (ptype) {
Matthieu Schaller's avatar
Matthieu Schaller committed
834

835
      case swift_type_gas:
836
        Nparticles = Ngas;
837
        hydro_write_particles(parts, list, &num_fields);
Matthieu Schaller's avatar
Matthieu Schaller committed
838
839
        break;

840
      case swift_type_dark_matter:
Matthieu Schaller's avatar
Matthieu Schaller committed
841
842
843
844
845
846
847
848
849
        /* 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 */
850
        io_collect_dm_gparts(gparts, Ntot, dmparts, Ndm);
Matthieu Schaller's avatar
Matthieu Schaller committed
851
852

        /* Write DM particles */
853
        Nparticles = Ndm;
854
        darkmatter_write_particles(dmparts, list, &num_fields);
Matthieu Schaller's avatar
Matthieu Schaller committed
855
856
        break;

857
      case swift_type_star:
Matthieu Schaller's avatar
Matthieu Schaller committed
858
859
860
        Nparticles = Nstars;
        star_write_particles(sparts, list, &num_fields);
        break;
861

Matthieu Schaller's avatar
Matthieu Schaller committed
862
863
      default:
        error("Particle Type %d not yet supported. Aborting", ptype);
864
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
865

866
867
    /* Write everything */
    for (int i = 0; i < num_fields; ++i)
868
869
870
      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
871

872
    /* Free temporary array */
873
874
875
876
    if (dmparts) {
      free(dmparts);
      dmparts = 0;
    }
877

878
879
    /* Close particle group */
    H5Gclose(h_grp);
Matthieu Schaller's avatar
Matthieu Schaller committed
880

881
    /* Close this particle group in the XMF file as well */
882
    if (mpi_rank == 0) xmf_write_groupfooter(xmfFile, (enum part_type)ptype);
883
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
884

885
  /* Write LXMF file descriptor */
886
  if (mpi_rank == 0) xmf_write_outputfooter(xmfFile, outputCount, e->time);
Matthieu Schaller's avatar
Matthieu Schaller committed
887

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

890
891
  /* Close property descriptor */
  H5Pclose(plist_id);
Matthieu Schaller's avatar
Matthieu Schaller committed
892

893
894
  /* Close file */
  H5Fclose(h_file);
Matthieu Schaller's avatar
Matthieu Schaller committed
895

896
897
898
  ++outputCount;
}

899
#endif /* HAVE_HDF5 */