single_io.c 30.6 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 33
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

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

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

/*-----------------------------------------------------------------------------
 * Routines reading an IC file
 *-----------------------------------------------------------------------------*/

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

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(h_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 89
    } else {
      /* message("Optional data set '%s' not present. Zeroing this particle
90
       * prop...", name);	   */
91

92
      for (size_t i = 0; i < N; ++i)
93
        memset(prop.field + i * prop.partSize, 0, copySize);
94 95

      return;
96
    }
97 98
  }

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

  /* Open data space */
104
  const hid_t h_data = H5Dopen(h_grp, prop.name, H5P_DEFAULT);
105
  if (h_data < 0) {
106
    error("Error while opening data space '%s'.", prop.name);
107
  }
108 109

  /* Check data type */
110
  const hid_t h_type = H5Dget_type(h_data);
111
  if (h_type < 0) error("Unable to retrieve data type from the file");
112
  // if (!H5Tequal(h_type, hdf5_type(type)))
113
  //  error("Non-matching types between the code and the file");
114

115
  /* Allocate temporary buffer */
116
  void* temp = malloc(num_elements * typeSize);
117
  if (temp == NULL) error("Unable to allocate memory for temporary buffer");
118 119 120 121

  /* Read HDF5 dataspace in temporary buffer */
  /* Dirty version that happens to work for vectors but should be improved */
  /* Using HDF5 dataspaces would be better */
122 123
  const hid_t h_err = H5Dread(h_data, io_hdf5_type(prop.type), H5S_ALL, H5S_ALL,
                              H5P_DEFAULT, temp);
124
  if (h_err < 0) {
125
    error("Error while reading data array '%s'.", prop.name);
126
  }
127

128
  /* Unit conversion if necessary */
129
  const double unit_factor =
130
      units_conversion_factor(ic_units, internal_units, prop.units);
131
  if (unit_factor != 1. && exist != 0) {
132

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

135
    if (io_is_double_precision(prop.type)) {
Matthieu Schaller's avatar
Matthieu Schaller committed
136
      double* temp_d = (double*)temp;
137
      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= unit_factor;
138
    } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
139
      float* temp_f = (float*)temp;
140 141 142 143 144 145
      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= unit_factor;
    }
  }

  /* Clean-up h if necessary */
  const float h_factor_exp = units_h_factor(internal_units, prop.units);
146
  if (cleanup_h && h_factor_exp != 0.f && exist != 0) {
147 148 149 150 151 152 153 154 155 156 157
    const double h_factor = pow(h, h_factor_exp);

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

    if (io_is_double_precision(prop.type)) {
      double* temp_d = (double*)temp;
      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= h_factor;
    } else {
      float* temp_f = (float*)temp;
      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= h_factor;
158 159 160
    }
  }

161
  /* Copy temporary buffer to particle data */
Matthieu Schaller's avatar
Matthieu Schaller committed
162
  char* temp_c = (char*)temp;
163
  for (size_t i = 0; i < N; ++i)
164
    memcpy(prop.field + i * prop.partSize, &temp_c[i * copySize], copySize);
165

166 167 168 169 170 171 172 173 174 175 176 177 178
  /* Free and close everything */
  free(temp);
  H5Tclose(h_type);
  H5Dclose(h_data);
}

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

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

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

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

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

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

213
  /* Create data space */
214 215 216 217
  const hid_t h_space = H5Screate(H5S_SIMPLE);
  int rank;
  hsize_t shape[2];
  hsize_t chunk_shape[2];
218
  if (h_space < 0) {
219
    error("Error while creating data space for field '%s'.", props.name);
220 221
  }

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

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

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

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

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

  /* Impose data compression */
256
  if (e->snapshotCompression > 0) {
257 258 259 260 261
    h_err = H5Pset_shuffle(h_prop);
    if (h_err < 0)
      error("Error while setting shuffling options for field '%s'.",
            props.name);

262
    h_err = H5Pset_deflate(h_prop, e->snapshotCompression);
263
    if (h_err < 0)
264 265
      error("Error while setting compression options for field '%s'.",
            props.name);
266 267
  }

268
  /* Create dataset */
269 270
  const hid_t h_data = H5Dcreate(grp, props.name, io_hdf5_type(props.type),
                                 h_space, H5P_DEFAULT, h_prop, H5P_DEFAULT);
271
  if (h_data < 0) {
272
    error("Error while creating dataspace '%s'.", props.name);
273 274
  }

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

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

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

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

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

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

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

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

  /* Read the relevant information */
365
  io_read_attribute(h_grp, "PeriodicBoundariesOn", INT, periodic);
366 367 368 369 370 371

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

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

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

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

393
  for (int ptype = 0; ptype < swift_type_count; ++ptype)
394 395 396
    N[ptype] = ((long long)numParticles[ptype]) +
               ((long long)numParticles_highWord[ptype] << 32);

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

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

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

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

  /* Close header */
  H5Gclose(h_grp);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    threadpool_clean(&tp);
  }
574

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

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

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

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

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

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

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

623
  /* File name */
624
  char fileName[FILENAME_BUFFER_SIZE];
625
  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
626
           e->snapshotOutputCount);
627 628

  /* First time, we need to create the XMF file */
629
  if (e->snapshotOutputCount == 0) xmf_create_file(baseName);
630

631
  /* Prepare the XMF file for the new entry */
632
  FILE* xmfFile = 0;
633
  xmfFile = xmf_prepare_file(baseName);
634 635

  /* Write the part corresponding to this specific output */
636
  xmf_write_outputheader(xmfFile, fileName, e->time);
637 638 639

  /* Open file */
  /* message("Opening file '%s'.", fileName); */
640 641 642 643
  h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
  if (h_file < 0) {
    error("Error while opening file '%s'.", fileName);
  }
644 645 646

  /* Open header to write simulation properties */
  /* message("Writing runtime parameters..."); */
647 648
  h_grp =
      H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
649
  if (h_grp < 0) error("Error while creating runtime parameters group\n");
650 651

  /* Write the relevant information */
652
  io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
653 654 655

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

657 658
  /* Open header to write simulation properties */
  /* message("Writing file header..."); */
659
  h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
660 661
  if (h_grp < 0) error("Error while creating file header\n");

662
  /* Print the relevant information and print status */
663
  io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
664
  double dblTime = e->time;
665
  io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
666
  int dimension = (int)hydro_dimension;
667
  io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
668 669
  io_write_attribute(h_grp, "Redshift", DOUBLE, &e->cosmology->z, 1);
  io_write_attribute(h_grp, "Scale-factor", DOUBLE, &e->cosmology->a, 1);
670 671

  /* GADGET-2 legacy values */
672
  /* Number of particles of each type */
673 674 675
  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
676 677 678
    numParticles[ptype] = (unsigned int)N_total[ptype];
    numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
  }
679 680 681 682 683 684 685 686 687
  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};
688
  flagEntropy[0] = writeEntropyFlag();
689 690 691
  io_write_attribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
                     swift_type_count);
  io_write_attribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
692 693 694 695

  /* Close header */
  H5Gclose(h_grp);

696
  /* Print the code version */
697
  io_write_code_description(h_file);
698

699 700 701
  /* Print the run's policy */
  io_write_engine_policy(h_file, e);

702
  /* Print the SPH parameters */
703 704 705
  if (e->policy & engine_policy_hydro) {
    h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
706 707
    if (h_grp < 0) error("Error while creating SPH group");
    hydro_props_print_snapshot(h_grp, e->hydro_properties);
708
    hydro_write_flavour(h_grp);
709 710 711
    H5Gclose(h_grp);
  }

712 713 714 715
  /* 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");
716 717
  cooling_write_flavour(h_grp);
  chemistry_write_flavour(h_grp);
718 719
  H5Gclose(h_grp);

720
  /* Print the gravity parameters */
721
  if (e->policy & engine_policy_self_gravity) {
722
    h_grp = H5Gcreate(h_file, "/GravityScheme", H5P_DEFAULT, H5P_DEFAULT,
723
                      H5P_DEFAULT);
724 725 726 727
    if (h_grp < 0) error("Error while creating gravity group");
    gravity_props_print_snapshot(h_grp, e->gravity_properties);
    H5Gclose(h_grp);
  }
728

729 730 731 732 733 734 735 736 737
  /* Print the cosmological model  */
  if (e->policy & engine_policy_cosmology) {
    h_grp =
        H5Gcreate(h_file, "/Cosmology", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    if (h_grp < 0) error("Error while creating cosmology group");
    cosmology_write_model(h_grp, e->cosmology);
    H5Gclose(h_grp);
  }

738
  /* Print the runtime parameters */
739 740
  h_grp =
      H5Gcreate(h_file, "/Parameters", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
741 742 743 744
  if (h_grp < 0) error("Error while creating parameters group");
  parser_write_params_to_hdf5(e->parameter_file, h_grp);
  H5Gclose(h_grp);

745
  /* Print the system of Units used in the spashot */
746
  io_write_unit_system(h_file, snapshot_units, "Units");
747 748

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

751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783
  /* 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);
    }
  }

784
  /* Loop over all particle types */
785
  for (int ptype = 0; ptype < swift_type_count; ptype++) {
786

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