common_io.c 101 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 25 26
/* This object's header. */
#include "common_io.h"

27 28 29
/* Pre-inclusion as needed in other headers */
#include "engine.h"

30
/* Local includes. */
31
#include "black_holes_io.h"
32
#include "chemistry_io.h"
Josh Borrow's avatar
Josh Borrow committed
33
#include "const.h"
34
#include "cooling_io.h"
35
#include "error.h"
36
#include "feedback.h"
37
#include "fof_io.h"
38
#include "gravity_io.h"
39
#include "hydro.h"
40
#include "hydro_io.h"
41
#include "io_properties.h"
42
#include "kernel_hydro.h"
43
#include "part.h"
lhausamm's avatar
lhausamm committed
44
#include "part_type.h"
Loic Hausammann's avatar
Loic Hausammann committed
45
#include "sink_io.h"
46
#include "star_formation_io.h"
47
#include "stars_io.h"
48
#include "threadpool.h"
49
#include "tools.h"
50
#include "tracers_io.h"
51
#include "units.h"
52
#include "velociraptor_io.h"
53
#include "version.h"
54

55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
/* Some standard headers. */
#include <math.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#if defined(HAVE_HDF5)

#include <hdf5.h>

/* MPI headers. */
#ifdef WITH_MPI
#include <mpi.h>
#endif

71
/**
72
 * @brief Converts a C data type to the HDF5 equivalent.
73 74
 *
 * This function is a trivial wrapper around the HDF5 types but allows
75 76
 * to change the exact storage types matching the code types in a transparent
 *way.
77
 */
78
hid_t io_hdf5_type(enum IO_DATA_TYPE type) {
79

80 81 82 83 84
  switch (type) {
    case INT:
      return H5T_NATIVE_INT;
    case UINT:
      return H5T_NATIVE_UINT;
Loic Hausammann's avatar
Loic Hausammann committed
85 86
    case UINT64:
      return H5T_NATIVE_UINT64;
87 88 89 90 91 92 93 94 95 96 97 98 99
    case LONG:
      return H5T_NATIVE_LONG;
    case ULONG:
      return H5T_NATIVE_ULONG;
    case LONGLONG:
      return H5T_NATIVE_LLONG;
    case ULONGLONG:
      return H5T_NATIVE_ULLONG;
    case FLOAT:
      return H5T_NATIVE_FLOAT;
    case DOUBLE:
      return H5T_NATIVE_DOUBLE;
    case CHAR:
100
      return H5T_NATIVE_CHAR;
101 102 103 104
    default:
      error("Unknown type");
      return 0;
  }
105 106
}

107 108 109 110 111
/**
 * @brief Return 1 if the type has double precision
 *
 * Returns an error if the type is not FLOAT or DOUBLE
 */
112
int io_is_double_precision(enum IO_DATA_TYPE type) {
Matthieu Schaller's avatar
Matthieu Schaller committed
113

114 115 116 117 118 119 120 121 122 123 124
  switch (type) {
    case FLOAT:
      return 0;
    case DOUBLE:
      return 1;
    default:
      error("Invalid type");
      return 0;
  }
}

125
/**
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
126
 * @brief Reads an attribute (scalar) from a given HDF5 group.
127 128 129
 *
 * @param grp The group from which to read.
 * @param name The name of the attribute to read.
130
 * @param type The #IO_DATA_TYPE of the attribute.
131 132 133 134
 * @param data (output) The attribute read from the HDF5 group.
 *
 * Calls #error() if an error occurs.
 */
135
void io_read_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
136
                       void* data) {
137

138 139
  const hid_t h_attr = H5Aopen(grp, name, H5P_DEFAULT);
  if (h_attr < 0) error("Error while opening attribute '%s'", name);
140

141 142
  const hid_t h_err = H5Aread(h_attr, io_hdf5_type(type), data);
  if (h_err < 0) error("Error while reading attribute '%s'", name);
143 144 145 146

  H5Aclose(h_attr);
}

Josh Borrow's avatar
Josh Borrow committed
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209
/**
 * @brief Reads an attribute from a given HDF5 group.
 *
 * @param grp The group from which to read.
 * @param name The name of the attribute to read.
 * @param type The #IO_DATA_TYPE of the attribute.
 * @param data (output) The attribute read from the HDF5 group.
 *
 * Exits gracefully (i.e. does not read the attribute at all) if
 * it is not present, unless debugging checks are activated. If they are,
 * and the read fails, we print a warning.
 */
void io_read_attribute_graceful(hid_t grp, const char* name,
                                enum IO_DATA_TYPE type, void* data) {

  /* First, we need to check if this attribute exists to avoid raising errors
   * within the HDF5 library if we attempt to access an attribute that does
   * not exist. */
  const htri_t h_exists = H5Aexists(grp, name);

  if (h_exists <= 0) {
  /* Attribute either does not exist (0) or function failed (-ve) */
#ifdef SWIFT_DEBUG_CHECKS
    message("WARNING: attribute '%s' does not exist.", name);
#endif
  } else {
    /* Ok, now we know that it exists we can read it. */
    const hid_t h_attr = H5Aopen(grp, name, H5P_DEFAULT);

    if (h_attr >= 0) {
      const hid_t h_err = H5Aread(h_attr, io_hdf5_type(type), data);
      if (h_err < 0) {
      /* Explicitly do nothing unless debugging checks are activated */
#ifdef SWIFT_DEBUG_CHECKS
        message("WARNING: unable to read attribute '%s'", name);
#endif
      }
    } else {
#ifdef SWIFT_DEBUG_CHECKS
      if (h_attr < 0) {
        message("WARNING: was unable to open attribute '%s'", name);
      }
#endif
    }

    H5Aclose(h_attr);
  }
}

/**
 * @brief Asserts that the redshift in the initial conditions and the one
 *        specified by the parameter file match.
 *
 * @param h_grp The Header group from the ICs
 * @param a Current scale factor as specified by parameter file
 */
void io_assert_valid_header_cosmology(hid_t h_grp, double a) {

  double redshift_from_snapshot = -1.0;
  io_read_attribute_graceful(h_grp, "Redshift", DOUBLE,
                             &redshift_from_snapshot);

  /* If the Header/Redshift value is not present, then we skip this check */
210
  if (redshift_from_snapshot == -1.0) return;
Josh Borrow's avatar
Josh Borrow committed
211 212

  const double current_redshift = 1.0 / a - 1.0;
213 214 215 216

  /* Escape if non-cosmological */
  if (current_redshift == 0.) return;

Josh Borrow's avatar
Josh Borrow committed
217 218 219 220 221 222 223 224 225 226 227
  const double redshift_fractional_difference =
      fabs(redshift_from_snapshot - current_redshift) / current_redshift;

  if (redshift_fractional_difference >= io_redshift_tolerance) {
    error(
        "Initial redshift specified in parameter file (%lf) and redshift "
        "read from initial conditions (%lf) are inconsistent.",
        current_redshift, redshift_from_snapshot);
  }
}

Loic Hausammann's avatar
Gear  
Loic Hausammann committed
228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262
/**
 * @brief Reads the number of elements in a HDF5 attribute.
 *
 * @param attr The attribute from which to read.
 *
 * @return The number of elements.
 *
 * Calls #error() if an error occurs.
 */
hsize_t io_get_number_element_in_attribute(hid_t attr) {
  /* Get the dataspace */
  hid_t space = H5Aget_space(attr);
  if (space < 0) error("Failed to get data space");

  /* Read the number of dimensions */
  const int ndims = H5Sget_simple_extent_ndims(space);

  /* Read the dimensions */
  hsize_t* dims = (hsize_t*)malloc(sizeof(hsize_t) * ndims);
  H5Sget_simple_extent_dims(space, dims, NULL);

  /* Compute number of elements */
  hsize_t count = 1;
  for (int i = 0; i < ndims; i++) {
    count *= dims[i];
  }

  /* Cleanup */
  free(dims);
  H5Sclose(space);
  return count;
};

/**
 * @brief Reads an attribute (array) from a given HDF5 group.
Matthieu Schaller's avatar
Matthieu Schaller committed
263 264 265 266
 *
 * @param grp The group from which to read.
 * @param name The name of the dataset to read.
 * @param type The #IO_DATA_TYPE of the attribute.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371
 * @param data (output) The attribute read from the HDF5 group (need to be
 * already allocated).
 * @param number_element Number of elements in the attribute.
 *
 * Calls #error() if an error occurs.
 */
void io_read_array_attribute(hid_t grp, const char* name,
                             enum IO_DATA_TYPE type, void* data,
                             hsize_t number_element) {

  /* Open attribute */
  const hid_t h_attr = H5Aopen(grp, name, H5P_DEFAULT);
  if (h_attr < 0) error("Error while opening attribute '%s'", name);

  /* Get the number of elements */
  hsize_t count = io_get_number_element_in_attribute(h_attr);

  /* Check if correct number of element */
  if (count != number_element) {
    error(
        "Error found a different number of elements than expected (%lli != "
        "%lli) in attribute %s",
        count, number_element, name);
  }

  /* Read attribute */
  const hid_t h_err = H5Aread(h_attr, io_hdf5_type(type), data);
  if (h_err < 0) error("Error while reading attribute '%s'", name);

  /* Cleanup */
  H5Aclose(h_attr);
}

/**
 * @brief Reads the number of elements in a HDF5 dataset.
 *
 * @param dataset The dataset from which to read.
 *
 * @return The number of elements.
 *
 * Calls #error() if an error occurs.
 */
hsize_t io_get_number_element_in_dataset(hid_t dataset) {
  /* Get the dataspace */
  hid_t space = H5Dget_space(dataset);
  if (space < 0) error("Failed to get data space");

  /* Read the number of dimensions */
  const int ndims = H5Sget_simple_extent_ndims(space);

  /* Read the dimensions */
  hsize_t* dims = (hsize_t*)malloc(sizeof(hsize_t) * ndims);
  H5Sget_simple_extent_dims(space, dims, NULL);

  /* Compute number of elements */
  hsize_t count = 1;
  for (int i = 0; i < ndims; i++) {
    count *= dims[i];
  }

  /* Cleanup */
  free(dims);
  H5Sclose(space);
  return count;
};

/**
 * @brief Reads a dataset (array) from a given HDF5 group.
 *
 * @param grp The group from which to read.
 * @param name The name of the dataset to read.
 * @param type The #IO_DATA_TYPE of the attribute.
 * @param data (output) The attribute read from the HDF5 group (need to be
 * already allocated).
 * @param number_element Number of elements in the attribute.
 *
 * Calls #error() if an error occurs.
 */
void io_read_array_dataset(hid_t grp, const char* name, enum IO_DATA_TYPE type,
                           void* data, hsize_t number_element) {

  /* Open dataset */
  const hid_t h_dataset = H5Dopen(grp, name, H5P_DEFAULT);
  if (h_dataset < 0) error("Error while opening attribute '%s'", name);

  /* Get the number of elements */
  hsize_t count = io_get_number_element_in_dataset(h_dataset);

  /* Check if correct number of element */
  if (count != number_element) {
    error(
        "Error found a different number of elements than expected (%lli != "
        "%lli) in dataset %s",
        count, number_element, name);
  }

  /* Read dataset */
  const hid_t h_err = H5Dread(h_dataset, io_hdf5_type(type), H5S_ALL, H5S_ALL,
                              H5P_DEFAULT, data);
  if (h_err < 0) error("Error while reading dataset '%s'", name);

  /* Cleanup */
  H5Dclose(h_dataset);
}

372 373 374 375 376
/**
 * @brief Write an attribute to a given HDF5 group.
 *
 * @param grp The group in which to write.
 * @param name The name of the attribute to write.
377
 * @param type The #IO_DATA_TYPE of the attribute.
378 379 380 381 382
 * @param data The attribute to write.
 * @param num The number of elements to write
 *
 * Calls #error() if an error occurs.
 */
383
void io_write_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
384
                        const void* data, int num) {
385

386 387
  const hid_t h_space = H5Screate(H5S_SIMPLE);
  if (h_space < 0)
388
    error("Error while creating dataspace for attribute '%s'.", name);
389

390 391 392
  hsize_t dim[1] = {(hsize_t)num};
  const hid_t h_err = H5Sset_extent_simple(h_space, 1, dim, NULL);
  if (h_err < 0)
393
    error("Error while changing dataspace shape for attribute '%s'.", name);
394

395 396 397
  const hid_t h_attr =
      H5Acreate1(grp, name, io_hdf5_type(type), h_space, H5P_DEFAULT);
  if (h_attr < 0) error("Error while creating attribute '%s'.", name);
398

399 400
  const hid_t h_err2 = H5Awrite(h_attr, io_hdf5_type(type), data);
  if (h_err2 < 0) error("Error while reading attribute '%s'.", name);
401 402 403 404 405 406 407 408 409 410 411 412 413 414 415

  H5Sclose(h_space);
  H5Aclose(h_attr);
}

/**
 * @brief Write a string as an attribute to a given HDF5 group.
 *
 * @param grp The group in which to write.
 * @param name The name of the attribute to write.
 * @param str The string to write.
 * @param length The length of the string
 *
 * Calls #error() if an error occurs.
 */
416 417
void io_writeStringAttribute(hid_t grp, const char* name, const char* str,
                             int length) {
418

419 420
  const hid_t h_space = H5Screate(H5S_SCALAR);
  if (h_space < 0)
421
    error("Error while creating dataspace for attribute '%s'.", name);
422

423 424
  const hid_t h_type = H5Tcopy(H5T_C_S1);
  if (h_type < 0) error("Error while copying datatype 'H5T_C_S1'.");
425

426 427
  const hid_t h_err = H5Tset_size(h_type, length);
  if (h_err < 0) error("Error while resizing attribute type to '%i'.", length);
428

429 430
  const hid_t h_attr = H5Acreate1(grp, name, h_type, h_space, H5P_DEFAULT);
  if (h_attr < 0) error("Error while creating attribute '%s'.", name);
431

432 433
  const hid_t h_err2 = H5Awrite(h_attr, h_type, str);
  if (h_err2 < 0) error("Error while reading attribute '%s'.", name);
434 435 436 437 438 439 440 441

  H5Tclose(h_type);
  H5Sclose(h_space);
  H5Aclose(h_attr);
}

/**
 * @brief Writes a double value as an attribute
442
 * @param grp The group in which to write
443 444 445
 * @param name The name of the attribute
 * @param data The value to write
 */
446 447
void io_write_attribute_d(hid_t grp, const char* name, double data) {
  io_write_attribute(grp, name, DOUBLE, &data, 1);
448 449 450 451
}

/**
 * @brief Writes a float value as an attribute
452
 * @param grp The group in which to write
453 454 455
 * @param name The name of the attribute
 * @param data The value to write
 */
456 457
void io_write_attribute_f(hid_t grp, const char* name, float data) {
  io_write_attribute(grp, name, FLOAT, &data, 1);
458 459 460 461
}

/**
 * @brief Writes an int value as an attribute
462
 * @param grp The group in which to write
463 464 465
 * @param name The name of the attribute
 * @param data The value to write
 */
466 467
void io_write_attribute_i(hid_t grp, const char* name, int data) {
  io_write_attribute(grp, name, INT, &data, 1);
468 469 470 471
}

/**
 * @brief Writes a long value as an attribute
472
 * @param grp The group in which to write
473 474 475
 * @param name The name of the attribute
 * @param data The value to write
 */
476 477
void io_write_attribute_l(hid_t grp, const char* name, long data) {
  io_write_attribute(grp, name, LONG, &data, 1);
478 479 480 481
}

/**
 * @brief Writes a string value as an attribute
482
 * @param grp The group in which to write
483 484 485
 * @param name The name of the attribute
 * @param str The string to write
 */
486 487
void io_write_attribute_s(hid_t grp, const char* name, const char* str) {
  io_writeStringAttribute(grp, name, str, strlen(str));
488 489
}

490
/**
491
 * @brief Writes the meta-data of the run to an open hdf5 snapshot file.
492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620
 *
 * @param h_file The opened hdf5 file.
 * @param e The #engine containing the meta-data.
 * @param internal_units The system of units used internally.
 * @param snapshot_units The system of units used in snapshots.
 */
void io_write_meta_data(hid_t h_file, const struct engine* e,
                        const struct unit_system* internal_units,
                        const struct unit_system* snapshot_units) {

  hid_t h_grp;

  /* Print the code version */
  io_write_code_description(h_file);

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

  /* Print the physical constants */
  phys_const_print_snapshot(h_file, e->physical_constants);

  /* Print the SPH parameters */
  if (e->policy & engine_policy_hydro) {
    h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
                      H5P_DEFAULT);
    if (h_grp < 0) error("Error while creating SPH group");
    hydro_props_print_snapshot(h_grp, e->hydro_properties);
    hydro_write_flavour(h_grp);
    H5Gclose(h_grp);
  }

  /* 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");
  hid_t h_grp_columns =
      H5Gcreate(h_grp, "NamedColumns", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
  if (h_grp_columns < 0) error("Error while creating named columns group");
  entropy_floor_write_flavour(h_grp);
  cooling_write_flavour(h_grp, h_grp_columns, e->cooling_func);
  chemistry_write_flavour(h_grp, h_grp_columns, e);
  tracers_write_flavour(h_grp);
  feedback_write_flavour(e->feedback_props, h_grp);
  H5Gclose(h_grp_columns);
  H5Gclose(h_grp);

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

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

  /* Print the cosmological model  */
  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);

  /* 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, /*write_used=*/1);
  H5Gclose(h_grp);

  /* Print the runtime unused parameters */
  h_grp = H5Gcreate(h_file, "/UnusedParameters", 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, /*write_used=*/0);
  H5Gclose(h_grp);

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

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

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

621 622
/**
 * @brief Reads the Unit System from an IC file.
623 624 625 626
 *
 * If the 'Units' group does not exist in the ICs, we will use the internal
 * system of units.
 *
627
 * @param h_file The (opened) HDF5 file from which to read.
628 629
 * @param ic_units The unit_system to fill.
 * @param internal_units The internal system of units to copy if needed.
630
 * @param mpi_rank The MPI rank we are on.
631
 */
632 633 634
void io_read_unit_system(hid_t h_file, struct unit_system* ic_units,
                         const struct unit_system* internal_units,
                         int mpi_rank) {
635 636 637

  /* First check if it exists as this is *not* required. */
  const htri_t exists = H5Lexists(h_file, "/Units", H5P_DEFAULT);
638

639
  if (exists == 0) {
640 641

    if (mpi_rank == 0)
642
      message("'Units' group not found in ICs. Assuming internal unit system.");
643

644
    units_copy(ic_units, internal_units);
645 646

    return;
647

648
  } else if (exists < 0) {
649
    error("Serious problem with 'Units' group in ICs. H5Lexists gives %d",
650
          exists);
651
  }
652

653
  if (mpi_rank == 0) message("Reading IC units from ICs.");
654
  hid_t h_grp = H5Gopen(h_file, "/Units", H5P_DEFAULT);
655 656

  /* Ok, Read the damn thing */
657
  io_read_attribute(h_grp, "Unit length in cgs (U_L)", DOUBLE,
658
                    &ic_units->UnitLength_in_cgs);
659
  io_read_attribute(h_grp, "Unit mass in cgs (U_M)", DOUBLE,
660
                    &ic_units->UnitMass_in_cgs);
661
  io_read_attribute(h_grp, "Unit time in cgs (U_t)", DOUBLE,
662
                    &ic_units->UnitTime_in_cgs);
663
  io_read_attribute(h_grp, "Unit current in cgs (U_I)", DOUBLE,
664
                    &ic_units->UnitCurrent_in_cgs);
665
  io_read_attribute(h_grp, "Unit temperature in cgs (U_T)", DOUBLE,
666
                    &ic_units->UnitTemperature_in_cgs);
667 668 669 670 671

  /* Clean up */
  H5Gclose(h_grp);
}

672 673 674
/**
 * @brief Writes the current Unit System
 * @param h_file The (opened) HDF5 file in which to write
675
 * @param us The unit_system to dump
676
 * @param groupName The name of the HDF5 group to write to
677
 */
678 679
void io_write_unit_system(hid_t h_file, const struct unit_system* us,
                          const char* groupName) {
680

681
  const hid_t h_grpunit = H5Gcreate1(h_file, groupName, 0);
682 683
  if (h_grpunit < 0) error("Error while creating Unit System group");

684 685 686 687 688 689 690 691 692 693
  io_write_attribute_d(h_grpunit, "Unit mass in cgs (U_M)",
                       units_get_base_unit(us, UNIT_MASS));
  io_write_attribute_d(h_grpunit, "Unit length in cgs (U_L)",
                       units_get_base_unit(us, UNIT_LENGTH));
  io_write_attribute_d(h_grpunit, "Unit time in cgs (U_t)",
                       units_get_base_unit(us, UNIT_TIME));
  io_write_attribute_d(h_grpunit, "Unit current in cgs (U_I)",
                       units_get_base_unit(us, UNIT_CURRENT));
  io_write_attribute_d(h_grpunit, "Unit temperature in cgs (U_T)",
                       units_get_base_unit(us, UNIT_TEMPERATURE));
694 695 696 697

  H5Gclose(h_grpunit);
}

698 699 700 701
/**
 * @brief Writes the code version to the file
 * @param h_file The (opened) HDF5 file in which to write
 */
702
void io_write_code_description(hid_t h_file) {
703

704
  const hid_t h_grpcode = H5Gcreate1(h_file, "/Code", 0);
705 706
  if (h_grpcode < 0) error("Error while creating code group");

707
  io_write_attribute_s(h_grpcode, "Code", "SWIFT");
708 709 710 711 712 713 714 715 716 717
  io_write_attribute_s(h_grpcode, "Code Version", package_version());
  io_write_attribute_s(h_grpcode, "Compiler Name", compiler_name());
  io_write_attribute_s(h_grpcode, "Compiler Version", compiler_version());
  io_write_attribute_s(h_grpcode, "Git Branch", git_branch());
  io_write_attribute_s(h_grpcode, "Git Revision", git_revision());
  io_write_attribute_s(h_grpcode, "Git Date", git_date());
  io_write_attribute_s(h_grpcode, "Configuration options",
                       configuration_options());
  io_write_attribute_s(h_grpcode, "CFLAGS", compilation_cflags());
  io_write_attribute_s(h_grpcode, "HDF5 library version", hdf5_version());
718
  io_write_attribute_s(h_grpcode, "Thread barriers", thread_barrier_version());
719
  io_write_attribute_s(h_grpcode, "Allocators", allocator_version());
Matthieu Schaller's avatar
Matthieu Schaller committed
720
#ifdef HAVE_FFTW
721
  io_write_attribute_s(h_grpcode, "FFTW library version", fftw3_version());
Matthieu Schaller's avatar
Matthieu Schaller committed
722
#endif
723 724 725
#ifdef HAVE_LIBGSL
  io_write_attribute_s(h_grpcode, "GSL library version", libgsl_version());
#endif
726
#ifdef WITH_MPI
727
  io_write_attribute_s(h_grpcode, "MPI library", mpi_version());
728
#ifdef HAVE_METIS
729
  io_write_attribute_s(h_grpcode, "METIS library version", metis_version());
730
#endif
731 732 733
#ifdef HAVE_PARMETIS
  io_write_attribute_s(h_grpcode, "ParMETIS library version",
                       parmetis_version());
734
#endif
735
#else
736
  io_write_attribute_s(h_grpcode, "MPI library", "Non-MPI version of SWIFT");
737
#endif
Loic Hausammann's avatar
Loic Hausammann committed
738
  io_write_attribute_i(h_grpcode, "RandomSeed", SWIFT_RANDOM_SEED_XOR);
739 740 741
  H5Gclose(h_grpcode);
}

742 743 744 745 746 747 748 749 750 751
/**
 * @brief Write the #engine policy to the file.
 * @param h_file File to write to.
 * @param e The #engine to read the policy from.
 */
void io_write_engine_policy(hid_t h_file, const struct engine* e) {

  const hid_t h_grp = H5Gcreate1(h_file, "/Policy", 0);
  if (h_grp < 0) error("Error while creating policy group");

752
  for (int i = 1; i < engine_maxpolicy; ++i)
753 754 755 756 757 758 759 760
    if (e->policy & (1 << i))
      io_write_attribute_i(h_grp, engine_policy_names[i + 1], 1);
    else
      io_write_attribute_i(h_grp, engine_policy_names[i + 1], 0);

  H5Gclose(h_grp);
}

761 762 763 764 765
static long long cell_count_non_inhibited_gas(const struct cell* c) {
  const int total_count = c->hydro.count;
  struct part* parts = c->hydro.parts;
  long long count = 0;
  for (int i = 0; i < total_count; ++i) {
766 767
    if ((parts[i].time_bin != time_bin_inhibited) &&
        (parts[i].time_bin != time_bin_not_created)) {
768 769 770 771 772 773 774 775 776 777 778
      ++count;
    }
  }
  return count;
}

static long long cell_count_non_inhibited_dark_matter(const struct cell* c) {
  const int total_count = c->grav.count;
  struct gpart* gparts = c->grav.parts;
  long long count = 0;
  for (int i = 0; i < total_count; ++i) {
779 780
    if ((gparts[i].time_bin != time_bin_inhibited) &&
        (gparts[i].time_bin != time_bin_not_created) &&
781 782 783 784 785 786 787
        (gparts[i].type == swift_type_dark_matter)) {
      ++count;
    }
  }
  return count;
}

788 789 790 791 792 793
static long long cell_count_non_inhibited_background_dark_matter(
    const struct cell* c) {
  const int total_count = c->grav.count;
  struct gpart* gparts = c->grav.parts;
  long long count = 0;
  for (int i = 0; i < total_count; ++i) {
794 795
    if ((gparts[i].time_bin != time_bin_inhibited) &&
        (gparts[i].time_bin != time_bin_not_created) &&
796 797 798 799 800 801 802
        (gparts[i].type == swift_type_dark_matter_background)) {
      ++count;
    }
  }
  return count;
}

803 804 805 806 807
static long long cell_count_non_inhibited_stars(const struct cell* c) {
  const int total_count = c->stars.count;
  struct spart* sparts = c->stars.parts;
  long long count = 0;
  for (int i = 0; i < total_count; ++i) {
808 809
    if ((sparts[i].time_bin != time_bin_inhibited) &&
        (sparts[i].time_bin != time_bin_not_created)) {
810 811 812 813 814 815 816 817 818 819 820
      ++count;
    }
  }
  return count;
}

static long long cell_count_non_inhibited_black_holes(const struct cell* c) {
  const int total_count = c->black_holes.count;
  struct bpart* bparts = c->black_holes.parts;
  long long count = 0;
  for (int i = 0; i < total_count; ++i) {
821 822
    if ((bparts[i].time_bin != time_bin_inhibited) &&
        (bparts[i].time_bin != time_bin_not_created)) {
823 824 825 826 827 828
      ++count;
    }
  }
  return count;
}

Loic Hausammann's avatar
Loic Hausammann committed
829 830 831 832 833 834 835 836 837 838 839 840 841
static long long cell_count_non_inhibited_sinks(const struct cell* c) {
  const int total_count = c->sinks.count;
  struct sink* sinks = c->sinks.parts;
  long long count = 0;
  for (int i = 0; i < total_count; ++i) {
    if ((sinks[i].time_bin != time_bin_inhibited) &&
        (sinks[i].time_bin != time_bin_not_created)) {
      ++count;
    }
  }
  return count;
}

842 843 844 845 846 847 848 849 850 851 852 853 854 855 856
/**
 * @brief Write a single 1D array to a hdf5 group.
 *
 * This creates a simple Nx1 array with a chunk size of 1024x1.
 * The Fletcher-32 filter is applied to the array.
 *
 * @param h_grp The open hdf5 group.
 * @param n The number of elements in the array.
 * @param array The data to write.
 * @param type The type of the data to write.
 * @param name The name of the array.
 * @param array_content The name of the parent group (only used for error
 * messages).
 */
void io_write_array(hid_t h_grp, const int n, const void* array,
857 858
                    const enum IO_DATA_TYPE type, const char* name,
                    const char* array_content) {
859

860 861
  /* Create memory space */
  const hsize_t shape[2] = {(hsize_t)n, 1};
862 863 864 865 866 867 868
  hid_t h_space = H5Screate(H5S_SIMPLE);
  if (h_space < 0)
    error("Error while creating data space for %s %s", name, array_content);
  hid_t h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
  if (h_err < 0)
    error("Error while changing shape of %s %s data space.", name,
          array_content);
869

Matthieu Schaller's avatar
Matthieu Schaller committed
870 871 872
  /* Dataset type */
  hid_t h_type = H5Tcopy(io_hdf5_type(type));

873 874 875 876 877 878 879 880 881 882 883 884 885 886
  const hsize_t chunk[2] = {(1024 > n ? n : 1024), 1};
  hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);
  h_err = H5Pset_chunk(h_prop, 1, chunk);
  if (h_err < 0)
    error("Error while setting chunk shapes of %s %s data space.", name,
          array_content);

  /* Impose check-sum to verify data corruption */
  h_err = H5Pset_fletcher32(h_prop);
  if (h_err < 0)
    error("Error while setting check-sum filter on %s %s data space.", name,
          array_content);

  /* Write */
Matthieu Schaller's avatar
Matthieu Schaller committed
887 888
  hid_t h_data =
      H5Dcreate(h_grp, name, h_type, h_space, H5P_DEFAULT, h_prop, H5P_DEFAULT);
889 890
  if (h_data < 0)
    error("Error while creating dataspace for %s %s.", name, array_content);
891 892
  h_err = H5Dwrite(h_data, io_hdf5_type(type), h_space, H5S_ALL, H5P_DEFAULT,
                   array);
893
  if (h_err < 0) error("Error while writing %s %s.", name, array_content);
Matthieu Schaller's avatar
Matthieu Schaller committed
894
  H5Tclose(h_type);
895
  H5Dclose(h_data);
896
  H5Pclose(h_prop);
897 898 899
  H5Sclose(h_space);
}

900 901
void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
                           const double pos_dithering[3],
902 903
                           const struct cell* cells_top, const int nr_cells,
                           const double width[3], const int nodeID,
904
                           const int distributed,
905
                           const long long global_counts[swift_type_count],
906
                           const long long global_offsets[swift_type_count],
907
                           const int num_fields[swift_type_count],
908 909 910
                           const struct unit_system* internal_units,
                           const struct unit_system* snapshot_units) {

911 912 913 914 915 916 917 918 919
#ifdef SWIFT_DEBUG_CHECKS
  if (distributed) {
    if (global_offsets[0] != 0 || global_offsets[1] != 0 ||
        global_offsets[2] != 0 || global_offsets[3] != 0 ||
        global_offsets[4] != 0 || global_offsets[5] != 0)
      error("Global offset non-zero in the distributed io case!");
  }
#endif

920 921 922 923
  /* Abort if we don't have any cells yet (i.e. haven't constructed the space)
   */
  if (nr_cells == 0) return;

924
  double cell_width[3] = {width[0], width[1], width[2]};
925 926 927

  /* Temporary memory for the cell-by-cell information */
  double* centres = NULL;
928
  centres = (double*)malloc(3 * nr_cells * sizeof(double));
929

930 931 932 933
  /* Temporary memory for the cell files ID */
  int* files = NULL;
  files = (int*)malloc(nr_cells * sizeof(int));

934
  /* Count of particles in each cell */
935 936
  long long *count_part = NULL, *count_gpart = NULL,
            *count_background_gpart = NULL, *count_spart = NULL,
Loic Hausammann's avatar
Loic Hausammann committed
937
            *count_bpart = NULL, *count_sink = NULL;
938 939
  count_part = (long long*)malloc(nr_cells * sizeof(long long));
  count_gpart = (long long*)malloc(nr_cells * sizeof(long long));
940
  count_background_gpart = (long long*)malloc(nr_cells * sizeof(long long));
941
  count_spart = (long long*)malloc(nr_cells * sizeof(long long));
942
  count_bpart = (long long*)malloc(nr_cells * sizeof(long long));
Loic Hausammann's avatar
Loic Hausammann committed
943
  count_sink = (long long*)malloc(nr_cells * sizeof(long long));
944 945

  /* Global offsets of particles in each cell */
946 947
  long long *offset_part = NULL, *offset_gpart = NULL,
            *offset_background_gpart = NULL, *offset_spart = NULL,
Loic Hausammann's avatar
Loic Hausammann committed
948
            *offset_bpart = NULL, *offset_sink = NULL;
949 950
  offset_part = (long long*)malloc(nr_cells * sizeof(long long));
  offset_gpart = (long long*)malloc(nr_cells * sizeof(long long));
951
  offset_background_gpart = (long long*)malloc(nr_cells * sizeof(long long));
952
  offset_spart = (long long*)malloc(nr_cells * sizeof(long long));
953
  offset_bpart = (long long*)malloc(nr_cells * sizeof(long long));
Loic Hausammann's avatar
Loic Hausammann committed
954
  offset_sink = (long long*)malloc(nr_cells * sizeof(long long));
955 956 957 958

  /* Offsets of the 0^th element */
  offset_part[0] = 0;
  offset_gpart[0] = 0;
959
  offset_background_gpart[0] = 0;
960
  offset_spart[0] = 0;
961
  offset_bpart[0] = 0;
Loic Hausammann's avatar
Loic Hausammann committed
962
  offset_sink[0] = 0;
963 964

  /* Collect the cell information of *local* cells */
965 966
  long long local_offset_part = 0;
  long long local_offset_gpart = 0;
967
  long long local_offset_background_gpart = 0;
968
  long long local_offset_spart = 0;
969
  long long local_offset_bpart = 0;
Loic Hausammann's avatar
Loic Hausammann committed
970
  long long local_offset_sink = 0;
971 972
  for (int i = 0; i < nr_cells; ++i) {

973 974 975 976 977 978 979 980
    /* Store in which file this cell will be found */
    if (distributed) {
      files[i] = cells_top[i].nodeID;
    } else {
      files[i] = 0;
    }

    /* Is the cell on this node (i.e. we have full information */
981 982 983
    if (cells_top[i].nodeID == nodeID) {

      /* Centre of each cell */
984 985 986
      centres[i * 3 + 0] = cells_top[i].loc[0] + cell_width[0] * 0.5;
      centres[i * 3 + 1] = cells_top[i].loc[1] + cell_width[1] * 0.5;
      centres[i * 3 + 2] = cells_top[i].loc[2] + cell_width[2] * 0.5;
987

988 989 990 991 992 993 994 995 996 997 998
      /* Undo the dithering since the particles will have this vector applied to
       * them */
      centres[i * 3 + 0] = centres[i * 3 + 0] - pos_dithering[0];
      centres[i * 3 + 1] = centres[i * 3 + 1] - pos_dithering[1];
      centres[i * 3 + 2] = centres[i * 3 + 2] - pos_dithering[2];

      /* Finish by box wrapping to match what is done to the particles */
      centres[i * 3 + 0] = box_wrap(centres[i * 3 + 0], 0.0, dim[0]);
      centres[i * 3 + 1] = box_wrap(centres[i * 3 + 1], 0.0, dim[1]);
      centres[i * 3 + 2] = box_wrap(centres[i * 3 + 2], 0.0, dim[2]);

999
      /* Count real particles that will be written */
1000 1001
      count_part[i] = cell_count_non_inhibited_gas(&cells_top[i]);
      count_gpart[i] = cell_count_non_inhibited_dark_matter(&cells_top[i]);
1002 1003
      count_background_gpart[i] =
          cell_count_non_inhibited_background_dark_matter(&cells_top[i]);
1004 1005
      count_spart[i] = cell_count_non_inhibited_stars(&cells_top[i]);
      count_bpart[i] = cell_count_non_inhibited_black_holes(&cells_top[i]);
Loic Hausammann's avatar
Loic Hausammann committed
1006
      count_sink[i] = cell_count_non_inhibited_sinks(&cells_top[i]);
1007 1008

      /* Offsets including the global offset of all particles on this MPI rank
1009 1010
       * Note that in the distributed case, the global offsets are 0 such that
       * we actually compute the offset in the file written by this rank. */
1011
      offset_part[i] = local_offset_part + global_offsets[swift_type_gas];
Matthieu Schaller's avatar
Matthieu Schaller committed
1012 1013
      offset_gpart[i] =
          local_offset_gpart + global_offsets[swift_type_dark_matter];
1014 1015 1016
      offset_background_gpart[i] =
          local_offset_background_gpart +
          global_offsets[swift_type_dark_matter_background];
1017
      offset_spart[i] = local_offset_spart + global_offsets[swift_type_stars];
1018 1019
      offset_bpart[i] =
          local_offset_bpart + global_offsets[swift_type_black_hole];
Loic Hausammann's avatar
Loic Hausammann committed
1020
      offset_sink[i] = local_offset_sink + global_offsets[swift_type_sink];
1021

1022 1023
      local_offset_part += count_part[i];
      local_offset_gpart += count_gpart[i];
1024
      local_offset_background_gpart += count_background_gpart[i];
1025
      local_offset_spart += count_spart[i];
1026
      local_offset_bpart += count_bpart[i];
Loic Hausammann's avatar
Loic Hausammann committed
1027
      local_offset_sink += count_sink[i];
1028

1029 1030 1031 1032 1033 1034 1035 1036 1037 1038
    } else {

      /* Just zero everything for the foregin cells */

      centres[i * 3 + 0] = 0.;
      centres[i * 3 + 1] = 0.;
      centres[i * 3 + 2] = 0.;

      count_part[i] = 0;
      count_gpart[i] = 0;
1039
      count_background_gpart[i] = 0;
1040
      count_spart[i] = 0;
1041
      count_bpart[i] = 0;
Loic Hausammann's avatar
Loic Hausammann committed
1042
      count_sink[i] = 0;
1043 1044 1045

      offset_part[i] = 0;
      offset_gpart[i] = 0;
1046
      offset_background_gpart[i] = 0;
1047
      offset_spart[i] = 0;
1048
      offset_bpart[i] = 0;
Loic Hausammann's avatar
Loic Hausammann committed
1049
      offset_sink[i] = 0;
1050 1051 1052
    }
  }

1053
#ifdef WITH_MPI
Matthieu Schaller's avatar
Matthieu Schaller committed
1054 1055
  /* Now, reduce all the arrays. Note that we use a bit-wise OR here. This
     is safe as we made sure only local cells have non-zero values. */
1056
  MPI_Allreduce(MPI_IN_PLACE, count_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
1057
                MPI_COMM_WORLD);
1058
  MPI_Allreduce(MPI_IN_PLACE, count_gpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
1059
                MPI_COMM_WORLD);
1060
  MPI_Allreduce(MPI_IN_PLACE, count_background_gpart, nr_cells,
1061
                MPI_LONG_LONG_INT, MPI_BOR, MPI_COMM_WORLD);
Loic Hausammann's avatar
Loic Hausammann committed
1062 1063
  MPI_Allreduce(MPI_IN_PLACE, count_sink, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
                MPI_COMM_WORLD);
1064
  MPI_Allreduce(MPI_IN_PLACE, count_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
1065
                MPI_COMM_WORLD);
1066
  MPI_Allreduce(MPI_IN_PLACE, count_bpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
1067
                MPI_COMM_WORLD);
1068 1069

  MPI_Allreduce(MPI_IN_PLACE, offset_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
1070
                MPI_COMM_WORLD);
1071
  MPI_Allreduce(MPI_IN_PLACE, offset_gpart, nr_cells, MPI_LONG_LONG_INT,
1072
                MPI_BOR, MPI_COMM_WORLD);
1073
  MPI_Allreduce(MPI_IN_PLACE, offset_background_gpart, nr_cells,
1074
                MPI_LONG_LONG_INT, MPI_BOR, MPI_COMM_WORLD);
Loic Hausammann's avatar
Loic Hausammann committed
1075 1076
  MPI_Allreduce(MPI_IN_PLACE, offset_sink, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
                MPI_COMM_WORLD);
1077
  MPI_Allreduce(MPI_IN_PLACE, offset_spart, nr_cells, MPI_LONG_LONG_INT,
1078
                MPI_BOR, MPI_COMM_WORLD);
1079
  MPI_Allreduce(MPI_IN_PLACE, offset_bpart, nr_cells, MPI_LONG_LONG_INT,
1080
                MPI_BOR, MPI_COMM_WORLD);
1081

1082
  /* For the centres we use a sum as MPI does not like bit-wise operations
1083
     on floating point numbers */
1084
  MPI_Allreduce(MPI_IN_PLACE, centres, 3 * nr_cells, MPI_DOUBLE, MPI_SUM,
1085
                MPI_COMM_WORLD);
1086 1087
#endif

1088
  /* When writing a single file, only rank 0 writes the meta-data */
1089
  if ((distributed) || (!distributed && nodeID == 0)) {
1090

1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108
    /* Unit conversion if necessary */
    const double factor = units_conversion_factor(
        internal_units, snapshot_units, UNIT_CONV_LENGTH);
    if (factor != 1.) {

      /* Convert the cell centres */
      for (int i = 0; i < nr_cells; ++i) {
        centres[i * 3 + 0] *= factor;
        centres[i * 3 + 1] *= factor;
        centres[i * 3 + 2] *= factor;
      }

      /* Convert the cell widths */
      cell_width[0] *= factor;
      cell_width[1] *= factor;
      cell_width[2] *= factor;
    }

Matthieu Schaller's avatar
Matthieu Schaller committed
1109 1110 1111 1112 1113
    /* Write some meta-information first */
    hid_t h_subgrp =
        H5Gcreate(h_grp, "Meta-data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    if (h_subgrp < 0) error("Error while creating meta-data sub-group");
    io_write_attribute(h_subgrp, "nr_cells", INT, &nr_cells, 1);
1114
    io_write_attribute(h_subgrp, "size", DOUBLE, cell_width, 3);
Matthieu Schaller's avatar
Matthieu Schaller committed
1115 1116 1117 1118
    io_write_attribute(h_subgrp, "dimension", INT, cdim, 3);
    H5Gclose(h_subgrp);

    /* Write the centres to the group */
Peter W. Draper's avatar
Peter W. Draper committed
1119
    hsize_t shape[2] = {(hsize_t)nr_cells, 3};
Matthieu Schaller's avatar
Matthieu Schaller committed
1120 1121 1122
    hid_t h_space = H5Screate(H5S_SIMPLE);
    if (h_space < 0) error("Error while creating data space for cell centres");
    hid_t h_err = H5Sset_extent_simple(h_space, 2, shape, shape);
1123 1124
    if (h_err < 0)
      error("Error while changing shape of gas offsets data space.");
Matthieu Schaller's avatar
Matthieu Schaller committed
1125 1126
    hid_t h_data = H5Dcreate(h_grp, "Centres", io_hdf5_type(DOUBLE), h_space,
                             H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1127
    if (h_data < 0) error("Error while creating dataspace for gas offsets.");
Matthieu Schaller's avatar
Matthieu Schaller committed
1128 1129 1130
    h_err = H5Dwrite(h_data, io_hdf5_type(DOUBLE), h_space, H5S_ALL,
                     H5P_DEFAULT, centres);
    if (h_err < 0) error("Error while writing centres.");
1131 1132 1133
    H5Dclose(h_data);
    H5Sclose(h_space);

1134 1135 1136 1137 1138 1139 1140 1141
    /* Group containing the offsets and counts for each particle type */
    hid_t h_grp_offsets = H5Gcreate(h_grp, "OffsetsInFile", H5P_DEFAULT,
                                    H5P_DEFAULT, H5P_DEFAULT);
    if (h_grp_offsets < 0) error("Error while creating offsets sub-group");
    hid_t h_grp_files =
        H5Gcreate(h_grp, "Files", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    if (h_grp_files < 0) error("Error while creating filess sub-group");
    hid_t h_grp_counts =
Matthieu Schaller's avatar
Matthieu Schaller committed
1142
        H5Gcreate(h_grp, "Counts", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
1143
    if (h_grp_counts < 0) error("Error while creating counts sub-group");
1144

1145
    if (global_counts[swift_type_gas] > 0 && num_fields[swift_type_gas] > 0) {
1146 1147 1148 1149 1150 1151
      io_write_array(h_grp_files, nr_cells, files, INT, "PartType0", "files");
      io_write_array(h_grp_offsets, nr_cells, offset_part, LONGLONG,
                     "PartType0", "offsets");
      io_write_array(h_grp_counts, nr_cells, count_part, LONGLONG, "PartType0",
                     "counts");
    }
1152

1153 1154
    if (global_counts[swift_type_dark_matter] > 0 &&
        num_fields[swift_type_dark_matter] > 0) {
1155 1156 1157 1158
      io_write_array(h_grp_files, nr_cells, files, INT, "PartType1", "files");
      io_write_array(h_grp_offsets, nr_cells, offset_gpart, LONGLONG,
                     "PartType1", "offsets");
      io_write_array(h_grp_counts, nr_cells, count_gpart, LONGLONG, "PartType1",
1159
                     "counts");
1160 1161
    }

1162 1163
    if (global_counts[swift_type_dark_matter_background] > 0 &&
        num_fields[swift_type_dark_matter_background] > 0) {
1164 1165 1166 1167 1168 1169
      io_write_array(h_grp_files, nr_cells, files, INT, "PartType2", "files");
      io_write_array(h_grp_offsets, nr_cells, offset_background_gpart, LONGLONG,
                     "PartType2", "offsets");
      io_write_array(h_grp_counts, nr_cells, count_background_gpart, LONGLONG,
                     "PartType2", "counts");
    }
1170

Loic Hausammann's avatar
Loic Hausammann committed
1171 1172 1173 1174 1175 1176 1177 1178
    if (global_counts[swift_type_sink] > 0 && num_fields[swift_type_sink] > 0) {
      io_write_array(h_grp_files, nr_cells, files, INT, "PartType3", "files");
      io_write_array(h_grp_offsets, nr_cells, offset_sink, LONGLONG,
                     "PartType3", "offsets");
      io_write_array(h_grp_counts, nr_cells, count_sink, LONGLONG, "PartType3",
                     "counts");
    }

1179 1180
    if (global_counts[swift_type_stars] > 0 &&
        num_fields[swift_type_stars] > 0) {
1181 1182 1183 1184 1185 1186
      io_write_array(h_grp_files, nr_cells, files, INT, "PartType4", "files");
      io_write_array(h_grp_offsets, nr_cells, offset_spart, LONGLONG,
                     "PartType4", "offsets");
      io_write_array(h_grp_counts, nr_cells, count_spart, LONGLONG, "PartType4",
                     "counts");
    }
1187

1188 1189
    if (global_counts[swift_type_black_hole] > 0 &&
        num_fields[swift_type_black_hole] > 0) {
1190 1191 1192 1193 1194 1195
      io_write_array(h_grp_files, nr_cells, files, INT, "PartType5", "files");
      io_write_array(h_grp_offsets, nr_cells, offset_bpart, LONGLONG,
                     "PartType5", "offsets");
      io_write_array(h_grp_counts, nr_cells, count_bpart, LONGLONG, "PartType5",
                     "counts");
    }
1196

1197 1198 1199
    H5Gclose(h_grp_offsets);
    H5Gclose(h_grp_files);
    H5Gclose(h_grp_counts);
1200 1201 1202 1203
  }

  /* Free everything we allocated */
  free(centres);
1204
  free(files);
1205 1206
  free(count_part);
  free(count_gpart);
1207
  free(count_background_gpart);
1208
  free(count_spart);
1209
  free(count_bpart);
Loic Hausammann's avatar
Loic Hausammann committed
1210
  free(count_sink);
1211 1212
  free(offset_part);
  free(offset_gpart);
1213
  free(offset_background_gpart);
1214
  free(offset_spart);
1215
  free(offset_bpart);
Loic Hausammann's avatar
Loic Hausammann committed
1216
  free(offset_sink);
1217 1218
}

1219 1220
#endif /* HAVE_HDF5 */

1221 1222 1223 1224 1225 1226 1227 1228 1229 1230
/**
 * @brief Returns the memory size of the data type
 */
size_t io_sizeof_type(enum IO_DATA_TYPE type) {

  switch (type) {
    case INT:
      return sizeof(int);
    case UINT:
      return sizeof(unsigned int);
Loic Hausammann's avatar
Loic Hausammann committed
1231 1232
    case UINT64:
      return sizeof(uint64_t);
1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246
    case LONG:
      return sizeof(long);
    case ULONG:
      return sizeof(unsigned long);
    case LONGLONG:
      return sizeof(long long);
    case ULONGLONG:
      return sizeof(unsigned long long);
    case FLOAT:
      return sizeof(float);
    case DOUBLE:
      return sizeof(double);
    case CHAR:
      return sizeof(char);
Loic Hausammann's avatar
Loic Hausammann committed
1247 1248
    case SIZE_T:
      return sizeof(size_t);
1249 1250 1251 1252 1253 1254
    default:
      error("Unknown type");
      return 0;
  }
}

1255 1256 1257
/**
 * @brief Mapper function to copy #part or #gpart fields into a buffer.
 */
1258
void io_copy_mapper(void* restrict temp, int N, void* restrict extra_data) {
1259

1260
  const struct io_props props = *((const struct io_props*)extra_data);
1261 1262 1263 1264
  const size_t typeSize = io_sizeof_type(props.type);
  const size_t copySize = typeSize * props.dimension;

  /* How far are we with this chunk? */
Matthieu Schaller's avatar
Matthieu Schaller committed
1265
  char* restrict temp_c = (char*)temp;
1266 1267
  const ptrdiff_t delta = (temp_c - props.start_temp_c) / copySize;

1268 1269 1270
  for (int k = 0; k < N; k++) {
    memcpy(&temp_c[k * copySize], props.field + (delta + k) * props.partSize,
           copySize);