single_io.c 32.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
24
25
26
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

#if defined(HAVE_HDF5) && !defined(WITH_MPI)

/* Some standard headers. */
27
28
29
#include <hdf5.h>
#include <math.h>
#include <stddef.h>
30
31
32
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
33
#include <time.h>
34

35
36
37
38
/* This object's header. */
#include "single_io.h"

/* Local includes. */
39
#include "chemistry_io.h"
40
#include "common_io.h"
41
#include "cooling_io.h"
42
#include "dimension.h"
43
#include "engine.h"
44
#include "error.h"
45
#include "gravity_io.h"
46
#include "gravity_properties.h"
47
#include "hydro_io.h"
48
#include "hydro_properties.h"
49
#include "io_properties.h"
50
51
#include "kernel_hydro.h"
#include "part.h"
lhausamm's avatar
lhausamm committed
52
#include "part_type.h"
53
#include "stars_io.h"
54
#include "units.h"
55
#include "xmf.h"
56
57
58
59

/**
 * @brief Reads a data array from a given HDF5 group.
 *
60
 * @param h_grp The group from which to read.
61
 * @param prop The #io_props of the field to read
62
 * @param N The number of particles.
63
64
 * @param internal_units The #unit_system used internally
 * @param ic_units The #unit_system used in the ICs
65
 * @param cleanup_h Are we removing h-factors from the ICs?
66
67
68
69
 * @param cleanup_sqrt_a Are we cleaning-up the sqrt(a) factors in the Gadget
 * IC velocities?
 * @param h The value of the reduced Hubble constant.
 * @param a The current value of the scale-factor.
70
 *
71
 * @todo A better version using HDF5 hyper-slabs to read the file directly into
72
 * the part array will be written once the structures have been stabilized.
73
 */
74
void readArray(hid_t h_grp, const struct io_props props, size_t N,
75
               const struct unit_system* internal_units,
76
77
               const struct unit_system* ic_units, int cleanup_h,
               int cleanup_sqrt_a, double h, double a) {
78

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

  /* Check whether the dataspace exists or not */
84
  const htri_t exist = H5Lexists(h_grp, props.name, 0);
85
  if (exist < 0) {
86
    error("Error while checking the existence of data set '%s'.", props.name);
87
  } else if (exist == 0) {
88
89
    if (props.importance == COMPULSORY) {
      error("Compulsory data set '%s' not present in the file.", props.name);
90
91
    } else {
      /* message("Optional data set '%s' not present. Zeroing this particle
92
       * props...", name);	   */
93

94
      for (size_t i = 0; i < N; ++i)
95
        memset(props.field + i * props.partSize, 0, copySize);
96
97

      return;
98
    }
99
100
  }

Matthieu Schaller's avatar
Matthieu Schaller committed
101
  /* message("Reading %s '%s' array...", */
102
103
  /*         props.importance == COMPULSORY ? "compulsory" : "optional  ", */
  /*         props.name); */
104
105

  /* Open data space */
106
107
  const hid_t h_data = H5Dopen(h_grp, props.name, H5P_DEFAULT);
  if (h_data < 0) error("Error while opening data space '%s'.", props.name);
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
114
115

  /* Read HDF5 dataspace in temporary buffer */
  /* Dirty version that happens to work for vectors but should be improved */
  /* Using HDF5 dataspaces would be better */
116
117
  const hid_t h_err = H5Dread(h_data, io_hdf5_type(props.type), H5S_ALL,
                              H5S_ALL, H5P_DEFAULT, temp);
118
  if (h_err < 0) error("Error while reading data array '%s'.", props.name);
119

120
  /* Unit conversion if necessary */
121
  const double unit_factor =
122
      units_conversion_factor(ic_units, internal_units, props.units);
123
  if (unit_factor != 1. && exist != 0) {
124

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

127
    if (io_is_double_precision(props.type)) {
Matthieu Schaller's avatar
Matthieu Schaller committed
128
      double* temp_d = (double*)temp;
129
      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= unit_factor;
130
    } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
131
      float* temp_f = (float*)temp;
132
133
134
135
136
      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= unit_factor;
    }
  }

  /* Clean-up h if necessary */
137
  const float h_factor_exp = units_h_factor(internal_units, props.units);
138
  if (cleanup_h && h_factor_exp != 0.f && exist != 0) {
139

140
    /* message("Multipltying '%s' by h^%f=%f", props.name, h_factor_exp,
141
142
     * h_factor); */

143
    if (io_is_double_precision(props.type)) {
144
      double* temp_d = (double*)temp;
145
      const double h_factor = pow(h, h_factor_exp);
146
147
148
      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= h_factor;
    } else {
      float* temp_f = (float*)temp;
149
      const float h_factor = pow(h, h_factor_exp);
150
      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= h_factor;
151
152
153
    }
  }

154
  /* Clean-up a if necessary */
155
  if (cleanup_sqrt_a && a != 1. && (strcmp(props.name, "Velocities") == 0)) {
156

157
    if (io_is_double_precision(props.type)) {
158
159
160
161
162
163
164
165
166
167
      double* temp_d = (double*)temp;
      const double vel_factor = sqrt(a);
      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= vel_factor;
    } else {
      float* temp_f = (float*)temp;
      const float vel_factor = sqrt(a);
      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= vel_factor;
    }
  }

168
  /* Copy temporary buffer to particle data */
Matthieu Schaller's avatar
Matthieu Schaller committed
169
  char* temp_c = (char*)temp;
170
  for (size_t i = 0; i < N; ++i)
171
    memcpy(props.field + i * props.partSize, &temp_c[i * copySize], copySize);
172

173
174
175
176
177
178
179
180
  /* Free and close everything */
  free(temp);
  H5Dclose(h_data);
}

/**
 * @brief Writes a data array in given HDF5 group.
 *
181
 * @param e The #engine we are writing from.
182
183
184
 * @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
Matthieu Schaller's avatar
Matthieu Schaller committed
185
 * @param partTypeGroupName The name of the group containing the particles in
186
187
 * the HDF5 file.
 * @param props The #io_props of the field to read
188
 * @param N The number of particles to write.
189
190
 * @param internal_units The #unit_system used internally
 * @param snapshot_units The #unit_system used in the snapshots
191
 *
192
 * @todo A better version using HDF5 hyper-slabs to write the file directly from
193
 * the part array will be written once the structures have been stabilized.
194
 */
195
196
197
void writeArray(const struct engine* e, hid_t grp, char* fileName,
                FILE* xmfFile, char* partTypeGroupName,
                const struct io_props props, size_t N,
198
199
                const struct unit_system* internal_units,
                const struct unit_system* snapshot_units) {
200

201
  const size_t typeSize = io_sizeof_type(props.type);
202
  const size_t num_elements = N * props.dimension;
203

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

  /* Allocate temporary buffer */
207
  void* temp = NULL;
208
  if (posix_memalign((void**)&temp, IO_BUFFER_ALIGNMENT,
209
                     num_elements * typeSize) != 0)
210
    error("Unable to allocate temporary i/o buffer");
211

212
213
  /* Copy the particle data to the temporary buffer */
  io_copy_temp_buffer(temp, e, props, N, internal_units, snapshot_units);
214

215
  /* Create data space */
216
  const hid_t h_space = H5Screate(H5S_SIMPLE);
217
218
219
  if (h_space < 0)
    error("Error while creating data space for field '%s'.", props.name);

220
221
222
  int rank;
  hsize_t shape[2];
  hsize_t chunk_shape[2];
223

224
  if (props.dimension > 1) {
225
226
    rank = 2;
    shape[0] = N;
227
    shape[1] = props.dimension;
228
    chunk_shape[0] = 1 << 20; /* Just a guess...*/
229
    chunk_shape[1] = props.dimension;
230
231
232
233
  } else {
    rank = 1;
    shape[0] = N;
    shape[1] = 0;
234
    chunk_shape[0] = 1 << 20; /* Just a guess...*/
235
    chunk_shape[1] = 0;
236
237
  }

238
239
  /* Make sure the chunks are not larger than the dataset */
  if (chunk_shape[0] > N) chunk_shape[0] = N;
240

241
  /* Change shape of data space */
242
  hid_t h_err = H5Sset_extent_simple(h_space, rank, shape, shape);
243
  if (h_err < 0)
244
    error("Error while changing data space shape for field '%s'.", props.name);
245

246
  /* Dataset properties */
247
  const hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);
248
249
250

  /* Set chunk size */
  h_err = H5Pset_chunk(h_prop, rank, chunk_shape);
251
  if (h_err < 0)
252
    error("Error while setting chunk size (%llu, %llu) for field '%s'.",
253
          chunk_shape[0], chunk_shape[1], props.name);
254

255
256
257
258
259
  /* Impose check-sum to verify data corruption */
  h_err = H5Pset_fletcher32(h_prop);
  if (h_err < 0)
    error("Error while setting checksum options for field '%s'.", props.name);

260
  /* Impose data compression */
261
  if (e->snapshot_compression > 0) {
262
263
264
265
266
    h_err = H5Pset_shuffle(h_prop);
    if (h_err < 0)
      error("Error while setting shuffling options for field '%s'.",
            props.name);

267
    h_err = H5Pset_deflate(h_prop, e->snapshot_compression);
268
    if (h_err < 0)
269
270
      error("Error while setting compression options for field '%s'.",
            props.name);
271
272
  }

273
  /* Create dataset */
274
275
  const hid_t h_data = H5Dcreate(grp, props.name, io_hdf5_type(props.type),
                                 h_space, H5P_DEFAULT, h_prop, H5P_DEFAULT);
276
  if (h_data < 0) error("Error while creating dataspace '%s'.", props.name);
277

278
  /* Write temporary buffer to HDF5 dataspace */
279
280
  h_err = H5Dwrite(h_data, io_hdf5_type(props.type), h_space, H5S_ALL,
                   H5P_DEFAULT, temp);
281
  if (h_err < 0) error("Error while writing data array '%s'.", props.name);
282
283

  /* Write XMF description for this data set */
284
285
  xmf_write_line(xmfFile, fileName, partTypeGroupName, props.name, N,
                 props.dimension, props.type);
286
287

  /* Write unit conversion factors for this data set */
288
  char buffer[FIELD_BUFFER_SIZE];
289
  units_cgs_conversion_string(buffer, snapshot_units, props.units);
290
291
292
  io_write_attribute_d(
      h_data, "CGS conversion factor",
      units_cgs_conversion_factor(snapshot_units, props.units));
293
  io_write_attribute_f(h_data, "h-scale exponent", 0);
294
295
296
  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);
297

298
299
  /* Free and close everything */
  free(temp);
300
  H5Pclose(h_prop);
301
302
303
304
  H5Dclose(h_data);
  H5Sclose(h_space);
}

305
306
307
308
/**
 * @brief Reads an HDF5 initial condition file (GADGET-3 type)
 *
 * @param fileName The file to read.
309
 * @param internal_units The system units used internally
310
 * @param dim (output) The dimension of the volume.
311
 * @param parts (output) Array of #part particles.
Matthieu Schaller's avatar
Matthieu Schaller committed
312
 * @param gparts (output) Array of #gpart particles.
313
 * @param sparts (output) Array of #spart particles.
314
 * @param Ngas (output) number of Gas particles read.
Matthieu Schaller's avatar
Matthieu Schaller committed
315
 * @param Ngparts (output) The number of #gpart read.
316
 * @param Nstars (output) The number of #spart read.
317
 * @param periodic (output) 1 if the volume is periodic, 0 if not.
318
 * @param flag_entropy (output) 1 if the ICs contained Entropy in the
319
320
 * InternalEnergy field
 * @param with_hydro Are we reading gas particles ?
Matthieu Schaller's avatar
Matthieu Schaller committed
321
 * @param with_gravity Are we reading/creating #gpart arrays ?
322
 * @param with_stars Are we reading star particles ?
323
 * @param cleanup_h Are we cleaning-up h-factors from the quantities we read?
324
325
 * @param cleanup_sqrt_a Are we cleaning-up the sqrt(a) factors in the Gadget
 * IC velocities?
326
 * @param h The value of the reduced Hubble constant to use for correction.
327
 * @param a The current value of the scale-factor.
328
 * @prarm n_threads The number of threads to use for the temporary threadpool.
329
 * @param dry_run If 1, don't read the particle. Only allocates the arrays.
330
331
332
333
334
335
336
337
 *
 * 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.
 */
338
339
340
void read_ic_single(const char* fileName,
                    const struct unit_system* internal_units, double dim[3],
                    struct part** parts, struct gpart** gparts,
341
342
                    struct spart** sparts, size_t* Ngas, size_t* Ngparts,
                    size_t* Nstars, int* periodic, int* flag_entropy,
343
                    int with_hydro, int with_gravity, int with_stars,
344
345
                    int cleanup_h, int cleanup_sqrt_a, double h, double a,
                    int n_threads, int dry_run) {
346

347
348
  hid_t h_file = 0, h_grp = 0;
  /* GADGET has only cubic boxes (in cosmological mode) */
349
350
  double boxSize[3] = {0.0, -1.0, -1.0};
  /* GADGET has 6 particle types. We only keep the type 0 & 1 for now...*/
351
352
  long long numParticles[swift_type_count] = {0};
  long long numParticles_highWord[swift_type_count] = {0};
353
  size_t N[swift_type_count] = {0};
354
  int dimension = 3; /* Assume 3D if nothing is specified */
355
  size_t Ndm = 0;
356
357
358
359

  /* Open file */
  /* message("Opening file '%s' as IC.", fileName); */
  h_file = H5Fopen(fileName, H5F_ACC_RDONLY, H5P_DEFAULT);
360
  if (h_file < 0) error("Error while opening file '%s'.", fileName);
361
362
363

  /* Open header to read simulation properties */
  /* message("Reading runtime parameters..."); */
364
  h_grp = H5Gopen(h_file, "/RuntimePars", H5P_DEFAULT);
365
366
367
  if (h_grp < 0) error("Error while opening runtime parameters\n");

  /* Read the relevant information */
368
  io_read_attribute(h_grp, "PeriodicBoundariesOn", INT, periodic);
369
370
371
372
373
374

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

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

378
379
380
381
  /* 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");
382
  if (hid_dim > 0) io_read_attribute(h_grp, "Dimension", INT, &dimension);
383
384
385
386
  if (dimension != hydro_dimension)
    error("ICs dimensionality (%dD) does not match code dimensionality (%dD)",
          dimension, (int)hydro_dimension);

387
  /* Read the relevant information and print status */
388
  int flag_entropy_temp[6];
389
  io_read_attribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy_temp);
390
  *flag_entropy = flag_entropy_temp[0];
391
  io_read_attribute(h_grp, "BoxSize", DOUBLE, boxSize);
392
393
  io_read_attribute(h_grp, "NumPart_Total", LONGLONG, numParticles);
  io_read_attribute(h_grp, "NumPart_Total_HighWord", LONGLONG,
394
                    numParticles_highWord);
395

396
  for (int ptype = 0; ptype < swift_type_count; ++ptype)
397
    N[ptype] = (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);
398

399
  /* Get the box size if not cubic */
400
401
402
403
  dim[0] = boxSize[0];
  dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
  dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];

404
405
406
407
408
409
  /* 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];

410
411
412
413
414
415
416
  /* Convert the box size if we want to clean-up h-factors */
  if (cleanup_h) {
    dim[0] /= h;
    dim[1] /= h;
    dim[2] /= h;
  }

417
  /* message("Found %d particles in a %speriodic box of size [%f %f %f].",  */
418
  /* 	  *N, (periodic ? "": "non-"), dim[0], dim[1], dim[2]);  */
419
420
421
422

  /* Close header */
  H5Gclose(h_grp);

423
  /* Read the unit system used in the ICs */
Matthieu Schaller's avatar
Matthieu Schaller committed
424
425
  struct unit_system* ic_units =
      (struct unit_system*)malloc(sizeof(struct unit_system));
426
  if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
427
  io_read_unit_system(h_file, ic_units, internal_units, 0);
428
429

  /* Tell the user if a conversion will be needed */
430
431
  if (units_are_equal(ic_units, internal_units)) {

Matthieu Schaller's avatar
Matthieu Schaller committed
432
    message("IC and internal units match. No conversion needed.");
433
434

  } else {
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457

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

458
  /* Convert the dimensions of the box */
Matthieu Schaller's avatar
Matthieu Schaller committed
459
460
461
  for (int j = 0; j < 3; j++)
    dim[j] *=
        units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
462

463
  /* Allocate memory to store SPH particles */
464
  if (with_hydro) {
465
    *Ngas = N[swift_type_gas];
Matthieu Schaller's avatar
Matthieu Schaller committed
466
467
    if (posix_memalign((void**)parts, part_align,
                       *Ngas * sizeof(struct part)) != 0)
468
469
470
      error("Error while allocating memory for SPH particles");
    bzero(*parts, *Ngas * sizeof(struct part));
  }
471

472
  /* Allocate memory to store star particles */
473
  if (with_stars) {
474
    *Nstars = N[swift_type_stars];
Matthieu Schaller's avatar
Matthieu Schaller committed
475
    if (posix_memalign((void**)sparts, spart_align,
476
                       *Nstars * sizeof(struct spart)) != 0)
477
      error("Error while allocating memory for stars particles");
478
479
480
481
482
    bzero(*sparts, *Nstars * sizeof(struct spart));
  }

  /* Allocate memory to store all gravity particles */
  if (with_gravity) {
483
484
485
    Ndm = N[swift_type_dark_matter];
    *Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
               N[swift_type_dark_matter] +
486
               (with_stars ? N[swift_type_stars] : 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
487
    if (posix_memalign((void**)gparts, gpart_align,
488
489
490
491
                       *Ngparts * sizeof(struct gpart)) != 0)
      error("Error while allocating memory for gravity particles");
    bzero(*gparts, *Ngparts * sizeof(struct gpart));
  }
492
493
494
495

  /* message("Allocated %8.2f MB for particles.", *N * sizeof(struct part) /
   * (1024.*1024.)); */

496
497
  /* message("BoxSize = %lf", dim[0]); */
  /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */
498

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

502
    /* Don't do anything if no particle of this kind */
503
    if (N[ptype] == 0) continue;
504
505

    /* Open the particle group in the file */
506
    char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE];
Matthieu Schaller's avatar
Matthieu Schaller committed
507
508
    snprintf(partTypeGroupName, PARTICLE_GROUP_BUFFER_SIZE, "/PartType%d",
             ptype);
509
    h_grp = H5Gopen(h_file, partTypeGroupName, H5P_DEFAULT);
510
    if (h_grp < 0)
511
512
      error("Error while opening particle group %s.", partTypeGroupName);

513
514
    int num_fields = 0;
    struct io_props list[100];
515
    size_t Nparticles = 0;
516

517
    /* Read particle fields into the structure */
518
519
    switch (ptype) {

520
      case swift_type_gas:
521
522
523
        if (with_hydro) {
          Nparticles = *Ngas;
          hydro_read_particles(*parts, list, &num_fields);
524
          num_fields += chemistry_read_particles(*parts, list + num_fields);
525
        }
526
527
        break;

528
      case swift_type_dark_matter:
529
530
531
532
        if (with_gravity) {
          Nparticles = Ndm;
          darkmatter_read_particles(*gparts, list, &num_fields);
        }
533
534
        break;

535
      case swift_type_stars:
536
537
        if (with_stars) {
          Nparticles = *Nstars;
538
          stars_read_particles(*sparts, list, &num_fields);
539
        }
540
541
        break;

542
      default:
543
        message("Particle Type %d not yet supported. Particles ignored", ptype);
544
545
    }

546
547
548
    /* Read everything */
    if (!dry_run)
      for (int i = 0; i < num_fields; ++i)
549
        readArray(h_grp, list[i], Nparticles, internal_units, ic_units,
550
                  cleanup_h, cleanup_sqrt_a, h, a);
551

552
553
554
555
    /* Close particle group */
    H5Gclose(h_grp);
  }

556
557
  /* Duplicate the parts for gravity */
  if (!dry_run && with_gravity) {
558

559
560
561
    /* Let's initialise a bit of thread parallelism here */
    struct threadpool tp;
    threadpool_init(&tp, n_threads);
562

563
564
565
566
    /* Prepare the DM particles */
    io_prepare_dm_gparts(&tp, *gparts, Ndm);

    /* Duplicate the hydro particles into gparts */
567
    if (with_hydro) io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas, Ndm);
568
569
570

    /* Duplicate the star particles into gparts */
    if (with_stars)
571
      io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
572
573
574

    threadpool_clean(&tp);
  }
575

576
577
  /* message("Done Reading particles..."); */

578
579
580
  /* Clean up */
  free(ic_units);

581
582
583
584
  /* Close file */
  H5Fclose(h_file);
}

585
586
587
588
/**
 * @brief Writes an HDF5 output file (GADGET-3 type) with its XMF descriptor
 *
 * @param e The engine containing all the system.
589
 * @param baseName The common part of the snapshot file name.
590
591
 * @param internal_units The #unit_system used internally
 * @param snapshot_units The #unit_system used in the snapshots
592
593
594
 *
 * Creates an HDF5 output file and writes the particles contained
 * in the engine. If such a file already exists, it is erased and replaced
595
 * by the new one.
596
597
598
599
600
 * The companion XMF file is also updated accordingly.
 *
 * Calls #error() if an error occurs.
 *
 */
601
void write_output_single(struct engine* e, const char* baseName,
602
603
                         const struct unit_system* internal_units,
                         const struct unit_system* snapshot_units) {
604

605
  hid_t h_file = 0, h_grp = 0;
Matthieu Schaller's avatar
Matthieu Schaller committed
606
  const size_t Ngas = e->s->nr_parts;
607
  const size_t Nstars = e->s->nr_sparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
608
  const size_t Ntot = e->s->nr_gparts;
609
610
  int periodic = e->s->periodic;
  int numFiles = 1;
611
612
613
  const struct part* parts = e->s->parts;
  const struct xpart* xparts = e->s->xparts;
  const struct gpart* gparts = e->s->gparts;
614
  struct gpart* dmparts = NULL;
615
  const struct spart* sparts = e->s->sparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
616
  const struct cooling_function_data* cooling = e->cooling_func;
lhausamm's avatar
lhausamm committed
617
  struct swift_params* params = e->parameter_file;
618

619
  /* Number of unassociated gparts */
620
  const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
621

Matthieu Schaller's avatar
Matthieu Schaller committed
622
623
  long long N_total[swift_type_count] = {
      (long long)Ngas, (long long)Ndm, 0, 0, (long long)Nstars, 0};
624

625
  /* File name */
626
  char fileName[FILENAME_BUFFER_SIZE];
627
628
  if (e->snapshot_label_delta == 1)
    snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
629
             e->snapshot_output_count + e->snapshot_label_first);
630
631
  else
    snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%06i.hdf5", baseName,
632
633
             e->snapshot_output_count * e->snapshot_label_delta +
                 e->snapshot_label_first);
634
635

  /* First time, we need to create the XMF file */
636
  if (e->snapshot_output_count == 0) xmf_create_file(baseName);
637

638
  /* Prepare the XMF file for the new entry */
639
  FILE* xmfFile = 0;
640
  xmfFile = xmf_prepare_file(baseName);
641
642

  /* Write the part corresponding to this specific output */
643
  xmf_write_outputheader(xmfFile, fileName, e->time);
644
645
646

  /* Open file */
  /* message("Opening file '%s'.", fileName); */
647
  h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
648
  if (h_file < 0) error("Error while opening file '%s'.", fileName);
649
650
651

  /* Open header to write simulation properties */
  /* message("Writing runtime parameters..."); */
652
653
  h_grp =
      H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
654
  if (h_grp < 0) error("Error while creating runtime parameters group\n");
655
656

  /* Write the relevant information */
657
  io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
658
659
660

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

662
663
  /* Open header to write simulation properties */
  /* message("Writing file header..."); */
664
  h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
665
666
  if (h_grp < 0) error("Error while creating file header\n");

667
  /* Print the relevant information and print status */
668
  io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
669
  double dblTime = e->time;
670
  io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
671
  int dimension = (int)hydro_dimension;
672
  io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
673
674
  io_write_attribute(h_grp, "Redshift", DOUBLE, &e->cosmology->z, 1);
  io_write_attribute(h_grp, "Scale-factor", DOUBLE, &e->cosmology->a, 1);
675
  io_write_attribute_s(h_grp, "Code", "SWIFT");
676
677
  time_t tm = time(NULL);
  io_write_attribute_s(h_grp, "Snapshot date", ctime(&tm));
678
679

  /* GADGET-2 legacy values */
680
  /* Number of particles of each type */
681
682
683
  unsigned int numParticles[swift_type_count] = {0};
  unsigned int numParticlesHighWord[swift_type_count] = {0};
  for (int ptype = 0; ptype < swift_type_count; ++ptype) {
Matthieu Schaller's avatar
Matthieu Schaller committed
684
685
686
    numParticles[ptype] = (unsigned int)N_total[ptype];
    numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
  }
687
688
689
690
691
692
693
694
695
  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);
  double MassTable[swift_type_count] = {0};
  io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
  unsigned int flagEntropy[swift_type_count] = {0};
696
  flagEntropy[0] = writeEntropyFlag();
697
698
699
  io_write_attribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
                     swift_type_count);
  io_write_attribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
700
701
702
703

  /* Close header */
  H5Gclose(h_grp);

704
  /* Print the code version */
705
  io_write_code_description(h_file);
706

707
708
709
  /* Print the run's policy */
  io_write_engine_policy(h_file, e);

710
  /* Print the SPH parameters */
711
712
713
  if (e->policy & engine_policy_hydro) {
    h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
714
715
    if (h_grp < 0) error("Error while creating SPH group");
    hydro_props_print_snapshot(h_grp, e->hydro_properties);
716
    hydro_write_flavour(h_grp);
717
718
719
    H5Gclose(h_grp);
  }

720
721
722
723
  /* Print the subgrid parameters */
  h_grp = H5Gcreate(h_file, "/SubgridScheme", H5P_DEFAULT, H5P_DEFAULT,
                    H5P_DEFAULT);
  if (h_grp < 0) error("Error while creating subgrid group");
724
725
  cooling_write_flavour(h_grp);
  chemistry_write_flavour(h_grp);
726
727
  H5Gclose(h_grp);

728
  /* Print the gravity parameters */
729
  if (e->policy & engine_policy_self_gravity) {
730
    h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
731
                      H5P_DEFAULT);
732
733
734
735
    if (h_grp < 0) error("Error while creating gravity group");
    gravity_props_print_snapshot(h_grp, e->gravity_properties);
    H5Gclose(h_grp);
  }
736

737
738
739
740
741
742
743
744
745
  /* Print the stellar parameters */
  if (e->policy & engine_policy_stars) {
    h_grp = H5Gcreate(h_file, "/StarsScheme", H5P_DEFAULT, H5P_DEFAULT,
		      H5P_DEFAULT);
    if (h_grp < 0) error("Error while creating stars group");
    stars_props_print_snapshot(h_grp, e->stars_properties);
    H5Gclose(h_grp);
  }

746
  /* Print the cosmological model  */
747
748
749
750
751
752
753
754
755
  h_grp =
      H5Gcreate(h_file, "/Cosmology", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
  if (h_grp < 0) error("Error while creating cosmology group");
  if (e->policy & engine_policy_cosmology)
    io_write_attribute_i(h_grp, "Cosmological run", 1);
  else
    io_write_attribute_i(h_grp, "Cosmological run", 0);
  cosmology_write_model(h_grp, e->cosmology);
  H5Gclose(h_grp);
756

757
  /* Print the runtime parameters */
758
759
  h_grp =
      H5Gcreate(h_file, "/Parameters", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
760
  if (h_grp < 0) error("Error while creating parameters group");
761
762
763
  parser_write_params_to_hdf5(e->parameter_file, h_grp, 1);
  H5Gclose(h_grp);

764
  /* Print the runtime unused parameters */
Matthieu Schaller's avatar
Matthieu Schaller committed
765
766
  h_grp = H5Gcreate(h_file, "/UnusedParameters", H5P_DEFAULT, H5P_DEFAULT,
                    H5P_DEFAULT);
767
768
  if (h_grp < 0) error("Error while creating parameters group");
  parser_write_params_to_hdf5(e->parameter_file, h_grp, 0);
769
770
  H5Gclose(h_grp);

771
  /* Print the system of Units used in the spashot */
772
  io_write_unit_system(h_file, snapshot_units, "Units");
773
774

  /* Print the system of Units used internally */
775
  io_write_unit_system(h_file, internal_units, "InternalCodeUnits");
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
806
807
808
809
  /* Tell the user if a conversion will be needed */
  if (e->verbose) {
    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);
    }
  }

810
  /* Loop over all particle types */
811
  for (int ptype = 0; ptype < swift_type_count; ptype++) {
812

813
814
    /* Don't do anything if no particle of this kind */
    if (numParticles[ptype] == 0) continue;
Matthieu Schaller's avatar
Matthieu Schaller committed
815

816
    /* Add the global information for that particle type to the XMF meta-file */
817
818
    xmf_write_groupheader(xmfFile, fileName, numParticles[ptype],
                          (enum part_type)ptype);
819

820
    /* Open the particle group in the file */
821
    char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE];
Matthieu Schaller's avatar
Matthieu Schaller committed
822
823
    snprintf(partTypeGroupName, PARTICLE_GROUP_BUFFER_SIZE, "/PartType%d",
             ptype);
824
825
    h_grp = H5Gcreate(h_file, partTypeGroupName, H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
826
    if (h_grp < 0) error("Error while creating particle group.\n");
827

828
829
    int num_fields = 0;
    struct io_props list[100];
830
    size_t N = 0;