parallel_io.c 31.5 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
55
56
57
58
59
60
61

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  /* Close header */
  H5Gclose(h_grp);

461
  /* Read the unit system used in the ICs */
462
  struct unit_system* ic_units = malloc(sizeof(struct unit_system));
463
  if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
464
  io_read_unit_system(h_file, ic_units);
465
466
467
468
469

  /* 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
470
      message("IC and internal units match. No conversion needed.");
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
496

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

753
  /* Print the SPH parameters */
754
755
756
  if (e->policy & engine_policy_hydro) {
    h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
757
758
759
760
761
    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);
  }
762

763
  /* Print the gravity parameters */
764
  if (e->policy & engine_policy_self_gravity) {
765
    h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
766
                      H5P_DEFAULT);
767
768
769
770
    if (h_grp < 0) error("Error while creating gravity group");
    gravity_props_print_snapshot(h_grp, e->gravity_properties);
    H5Gclose(h_grp);
  }
771
772
773
774
775
776
777

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

779
  /* Print the system of Units used in the spashot */
780
  io_write_unit_system(h_file, snapshot_units, "Units");
781
782

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

785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
  /* 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);
    }
  }

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

821
822
    /* Don't do anything if no particle of this kind */
    if (N_total[ptype] == 0) continue;
Matthieu Schaller's avatar
Matthieu Schaller committed
823
824
825

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

830
831
832
    /* 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
833
834
835
             ptype);
    h_grp = H5Gcreate(h_file, partTypeGroupName, H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
836
837
838
    if (h_grp < 0) {
      error("Error while opening particle group %s.", partTypeGroupName);
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
839

840
841
    int num_fields = 0;
    struct io_props list[100];
842
    size_t Nparticles = 0;
843
844

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

847
      case swift_type_gas:
848
        Nparticles = Ngas;
849
        hydro_write_particles(parts, list, &num_fields);
Matthieu Schaller's avatar
Matthieu Schaller committed
850
851
        break;

852
      case swift_type_dark_matter:
Matthieu Schaller's avatar
Matthieu Schaller committed
853
854
855
856
857
858
859
860
861
        /* 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 */
862
        io_collect_dm_gparts(gparts, Ntot, dmparts, Ndm);
Matthieu Schaller's avatar
Matthieu Schaller committed
863
864

        /* Write DM particles */
865
        Nparticles = Ndm;
866
        darkmatter_write_particles(dmparts, list, &num_fields);
Matthieu Schaller's avatar
Matthieu Schaller committed
867
868
        break;

869
      case swift_type_star:
Matthieu Schaller's avatar
Matthieu Schaller committed
870
871
872
        Nparticles = Nstars;
        star_write_particles(sparts, list, &num_fields);
        break;
873

Matthieu Schaller's avatar
Matthieu Schaller committed
874
875
      default:
        error("Particle Type %d not yet supported. Aborting", ptype);
876
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
877

878
879
    /* Write everything */
    for (int i = 0; i < num_fields; ++i)
880
881
882
      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
883

884
    /* Free temporary array */
885
886
887
888
    if (dmparts) {
      free(dmparts);
      dmparts = 0;
    }
889

890
891
    /* Close particle group */
    H5Gclose(h_grp);
Matthieu Schaller's avatar
Matthieu Schaller committed
892

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

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

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

902
903
  /* Close property descriptor */
  H5Pclose(plist_id);
Matthieu Schaller's avatar
Matthieu Schaller committed
904