common_io.c 83.5 KB
Newer Older
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
4
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk).
5
 *
6
7
8
9
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
10
 *
11
12
13
14
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
15
 *
16
17
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
18
 *
19
20
21
22
23
 ******************************************************************************/

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

24
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 "fof_io.h"
37
#include "gravity_io.h"
38
#include "hydro.h"
39
#include "hydro_io.h"
40
#include "io_properties.h"
41
#include "kernel_hydro.h"
42
#include "part.h"
lhausamm's avatar
lhausamm committed
43
#include "part_type.h"
44
#include "star_formation_io.h"
45
#include "stars_io.h"
46
#include "threadpool.h"
47
#include "tracers_io.h"
48
#include "units.h"
49
#include "velociraptor_io.h"
50
#include "version.h"
51

52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
/* 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

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

77
78
79
80
81
  switch (type) {
    case INT:
      return H5T_NATIVE_INT;
    case UINT:
      return H5T_NATIVE_UINT;
Loic Hausammann's avatar
Loic Hausammann committed
82
83
    case UINT64:
      return H5T_NATIVE_UINT64;
84
85
86
87
88
89
90
91
92
93
94
95
96
    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:
97
      return H5T_NATIVE_CHAR;
98
99
100
101
    default:
      error("Unknown type");
      return 0;
  }
102
103
}

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

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

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

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

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

  H5Aclose(h_attr);
}

Josh Borrow's avatar
Josh Borrow committed
144
145
146
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
210
211
212
213
214
215
216
217
218
219
220
221
222
/**
 * @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 */
  if (redshift_from_snapshot == -1.0) {
    return;
  }

  const double current_redshift = 1.0 / a - 1.0;
  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
223
224
225
226
227
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
/**
 * @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
258
259
260
261
 *
 * @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
262
263
264
265
266
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
 * @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);
}

367
368
369
370
371
/**
 * @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.
372
 * @param type The #IO_DATA_TYPE of the attribute.
373
374
375
376
377
 * @param data The attribute to write.
 * @param num The number of elements to write
 *
 * Calls #error() if an error occurs.
 */
378
void io_write_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
379
                        const void* data, int num) {
380

381
382
  const hid_t h_space = H5Screate(H5S_SIMPLE);
  if (h_space < 0)
383
    error("Error while creating dataspace for attribute '%s'.", name);
384

385
386
387
  hsize_t dim[1] = {(hsize_t)num};
  const hid_t h_err = H5Sset_extent_simple(h_space, 1, dim, NULL);
  if (h_err < 0)
388
    error("Error while changing dataspace shape for attribute '%s'.", name);
389

390
391
392
  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);
393

394
395
  const hid_t h_err2 = H5Awrite(h_attr, io_hdf5_type(type), data);
  if (h_err2 < 0) error("Error while reading attribute '%s'.", name);
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410

  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.
 */
411
412
void io_writeStringAttribute(hid_t grp, const char* name, const char* str,
                             int length) {
413

414
415
  const hid_t h_space = H5Screate(H5S_SCALAR);
  if (h_space < 0)
416
    error("Error while creating dataspace for attribute '%s'.", name);
417

418
419
  const hid_t h_type = H5Tcopy(H5T_C_S1);
  if (h_type < 0) error("Error while copying datatype 'H5T_C_S1'.");
420

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

424
425
  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);
426

427
428
  const hid_t h_err2 = H5Awrite(h_attr, h_type, str);
  if (h_err2 < 0) error("Error while reading attribute '%s'.", name);
429
430
431
432
433
434
435
436

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

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

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

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

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

/**
 * @brief Writes a string value as an attribute
477
 * @param grp The group in which to write
478
479
480
 * @param name The name of the attribute
 * @param str The string to write
 */
481
482
void io_write_attribute_s(hid_t grp, const char* name, const char* str) {
  io_writeStringAttribute(grp, name, str, strlen(str));
483
484
}

485
486
/**
 * @brief Reads the Unit System from an IC file.
487
488
489
490
 *
 * If the 'Units' group does not exist in the ICs, we will use the internal
 * system of units.
 *
491
 * @param h_file The (opened) HDF5 file from which to read.
492
493
 * @param ic_units The unit_system to fill.
 * @param internal_units The internal system of units to copy if needed.
494
 * @param mpi_rank The MPI rank we are on.
495
 */
496
497
498
void io_read_unit_system(hid_t h_file, struct unit_system* ic_units,
                         const struct unit_system* internal_units,
                         int mpi_rank) {
499
500
501

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

503
  if (exists == 0) {
504
505

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

508
    units_copy(ic_units, internal_units);
509
510

    return;
511

512
  } else if (exists < 0) {
513
    error("Serious problem with 'Units' group in ICs. H5Lexists gives %d",
514
          exists);
515
  }
516

517
  if (mpi_rank == 0) message("Reading IC units from ICs.");
518
  hid_t h_grp = H5Gopen(h_file, "/Units", H5P_DEFAULT);
519
520

  /* Ok, Read the damn thing */
521
  io_read_attribute(h_grp, "Unit length in cgs (U_L)", DOUBLE,
522
                    &ic_units->UnitLength_in_cgs);
523
  io_read_attribute(h_grp, "Unit mass in cgs (U_M)", DOUBLE,
524
                    &ic_units->UnitMass_in_cgs);
525
  io_read_attribute(h_grp, "Unit time in cgs (U_t)", DOUBLE,
526
                    &ic_units->UnitTime_in_cgs);
527
  io_read_attribute(h_grp, "Unit current in cgs (U_I)", DOUBLE,
528
                    &ic_units->UnitCurrent_in_cgs);
529
  io_read_attribute(h_grp, "Unit temperature in cgs (U_T)", DOUBLE,
530
                    &ic_units->UnitTemperature_in_cgs);
531
532
533
534
535

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

536
537
538
/**
 * @brief Writes the current Unit System
 * @param h_file The (opened) HDF5 file in which to write
539
 * @param us The unit_system to dump
540
 * @param groupName The name of the HDF5 group to write to
541
 */
542
543
void io_write_unit_system(hid_t h_file, const struct unit_system* us,
                          const char* groupName) {
544

545
  const hid_t h_grpunit = H5Gcreate1(h_file, groupName, 0);
546
547
  if (h_grpunit < 0) error("Error while creating Unit System group");

548
549
550
551
552
553
554
555
556
557
  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));
558
559
560
561

  H5Gclose(h_grpunit);
}

562
563
564
565
/**
 * @brief Writes the code version to the file
 * @param h_file The (opened) HDF5 file in which to write
 */
566
void io_write_code_description(hid_t h_file) {
567

568
  const hid_t h_grpcode = H5Gcreate1(h_file, "/Code", 0);
569
570
  if (h_grpcode < 0) error("Error while creating code group");

571
  io_write_attribute_s(h_grpcode, "Code", "SWIFT");
572
573
574
575
576
577
578
579
580
581
  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());
582
  io_write_attribute_s(h_grpcode, "Thread barriers", thread_barrier_version());
583
  io_write_attribute_s(h_grpcode, "Allocators", allocator_version());
Matthieu Schaller's avatar
Matthieu Schaller committed
584
#ifdef HAVE_FFTW
585
  io_write_attribute_s(h_grpcode, "FFTW library version", fftw3_version());
Matthieu Schaller's avatar
Matthieu Schaller committed
586
#endif
587
588
589
#ifdef HAVE_LIBGSL
  io_write_attribute_s(h_grpcode, "GSL library version", libgsl_version());
#endif
590
#ifdef WITH_MPI
591
  io_write_attribute_s(h_grpcode, "MPI library", mpi_version());
592
#ifdef HAVE_METIS
593
  io_write_attribute_s(h_grpcode, "METIS library version", metis_version());
594
#endif
595
596
597
#ifdef HAVE_PARMETIS
  io_write_attribute_s(h_grpcode, "ParMETIS library version",
                       parmetis_version());
598
#endif
599
#else
600
  io_write_attribute_s(h_grpcode, "MPI library", "Non-MPI version of SWIFT");
601
#endif
Loic Hausammann's avatar
Loic Hausammann committed
602
  io_write_attribute_i(h_grpcode, "RandomSeed", SWIFT_RANDOM_SEED_XOR);
603
604
605
  H5Gclose(h_grpcode);
}

606
607
608
609
610
611
612
613
614
615
/**
 * @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");

616
  for (int i = 1; i < engine_maxpolicy; ++i)
617
618
619
620
621
622
623
624
    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);
}

625
626
627
628
629
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) {
630
631
    if ((parts[i].time_bin != time_bin_inhibited) &&
        (parts[i].time_bin != time_bin_not_created)) {
632
633
634
635
636
637
638
639
640
641
642
      ++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) {
643
644
    if ((gparts[i].time_bin != time_bin_inhibited) &&
        (gparts[i].time_bin != time_bin_not_created) &&
645
646
647
648
649
650
651
        (gparts[i].type == swift_type_dark_matter)) {
      ++count;
    }
  }
  return count;
}

652
653
654
655
656
657
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) {
658
659
    if ((gparts[i].time_bin != time_bin_inhibited) &&
        (gparts[i].time_bin != time_bin_not_created) &&
660
661
662
663
664
665
666
        (gparts[i].type == swift_type_dark_matter_background)) {
      ++count;
    }
  }
  return count;
}

667
668
669
670
671
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) {
672
673
    if ((sparts[i].time_bin != time_bin_inhibited) &&
        (sparts[i].time_bin != time_bin_not_created)) {
674
675
676
677
678
679
680
681
682
683
684
      ++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) {
685
686
    if ((bparts[i].time_bin != time_bin_inhibited) &&
        (bparts[i].time_bin != time_bin_not_created)) {
687
688
689
690
691
692
      ++count;
    }
  }
  return count;
}

693
694
void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
                           const double pos_dithering[3],
695
696
697
                           const struct cell* cells_top, const int nr_cells,
                           const double width[3], const int nodeID,
                           const long long global_counts[swift_type_count],
698
699
700
701
702
                           const long long global_offsets[swift_type_count],
                           const struct unit_system* internal_units,
                           const struct unit_system* snapshot_units) {

  double cell_width[3] = {width[0], width[1], width[2]};
703
704
705

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

708
  /* Count of particles in each cell */
709
710
  long long *count_part = NULL, *count_gpart = NULL,
            *count_background_gpart = NULL, *count_spart = NULL,
711
            *count_bpart = NULL;
712
713
  count_part = (long long*)malloc(nr_cells * sizeof(long long));
  count_gpart = (long long*)malloc(nr_cells * sizeof(long long));
714
  count_background_gpart = (long long*)malloc(nr_cells * sizeof(long long));
715
  count_spart = (long long*)malloc(nr_cells * sizeof(long long));
716
  count_bpart = (long long*)malloc(nr_cells * sizeof(long long));
717
718

  /* Global offsets of particles in each cell */
719
720
  long long *offset_part = NULL, *offset_gpart = NULL,
            *offset_background_gpart = NULL, *offset_spart = NULL,
721
            *offset_bpart = NULL;
722
723
  offset_part = (long long*)malloc(nr_cells * sizeof(long long));
  offset_gpart = (long long*)malloc(nr_cells * sizeof(long long));
724
  offset_background_gpart = (long long*)malloc(nr_cells * sizeof(long long));
725
  offset_spart = (long long*)malloc(nr_cells * sizeof(long long));
726
  offset_bpart = (long long*)malloc(nr_cells * sizeof(long long));
727
728
729
730

  /* Offsets of the 0^th element */
  offset_part[0] = 0;
  offset_gpart[0] = 0;
731
  offset_background_gpart[0] = 0;
732
  offset_spart[0] = 0;
733
  offset_bpart[0] = 0;
734
735

  /* Collect the cell information of *local* cells */
736
737
  long long local_offset_part = 0;
  long long local_offset_gpart = 0;
738
  long long local_offset_background_gpart = 0;
739
  long long local_offset_spart = 0;
740
  long long local_offset_bpart = 0;
741
742
743
744
745
  for (int i = 0; i < nr_cells; ++i) {

    if (cells_top[i].nodeID == nodeID) {

      /* Centre of each cell */
746
747
748
      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;
749

750
751
752
753
754
755
756
757
758
759
760
      /* 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]);

761
      /* Count real particles that will be written */
762
763
      count_part[i] = cell_count_non_inhibited_gas(&cells_top[i]);
      count_gpart[i] = cell_count_non_inhibited_dark_matter(&cells_top[i]);
764
765
      count_background_gpart[i] =
          cell_count_non_inhibited_background_dark_matter(&cells_top[i]);
766
767
      count_spart[i] = cell_count_non_inhibited_stars(&cells_top[i]);
      count_bpart[i] = cell_count_non_inhibited_black_holes(&cells_top[i]);
768
769
770

      /* Offsets including the global offset of all particles on this MPI rank
       */
771
      offset_part[i] = local_offset_part + global_offsets[swift_type_gas];
Matthieu Schaller's avatar
Matthieu Schaller committed
772
773
      offset_gpart[i] =
          local_offset_gpart + global_offsets[swift_type_dark_matter];
774
775
776
      offset_background_gpart[i] =
          local_offset_background_gpart +
          global_offsets[swift_type_dark_matter_background];
777
      offset_spart[i] = local_offset_spart + global_offsets[swift_type_stars];
778
779
      offset_bpart[i] =
          local_offset_bpart + global_offsets[swift_type_black_hole];
780

781
782
      local_offset_part += count_part[i];
      local_offset_gpart += count_gpart[i];
783
      local_offset_background_gpart += count_background_gpart[i];
784
      local_offset_spart += count_spart[i];
785
      local_offset_bpart += count_bpart[i];
786

787
788
789
790
791
792
793
794
795
796
    } 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;
797
      count_background_gpart[i] = 0;
798
      count_spart[i] = 0;
799
      count_bpart[i] = 0;
800
801
802

      offset_part[i] = 0;
      offset_gpart[i] = 0;
803
      offset_background_gpart[i] = 0;
804
      offset_spart[i] = 0;
805
      offset_bpart[i] = 0;
806
807
808
    }
  }

809
#ifdef WITH_MPI
Matthieu Schaller's avatar
Matthieu Schaller committed
810
811
  /* 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. */
812
813
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, count_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
Matthieu Schaller's avatar
Matthieu Schaller committed
814
               0, MPI_COMM_WORLD);
815
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
816
817
    MPI_Reduce(count_part, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
818
819
820
  }
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, count_gpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
Matthieu Schaller's avatar
Matthieu Schaller committed
821
               0, MPI_COMM_WORLD);
822
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
823
824
    MPI_Reduce(count_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
825
  }
826
827
828
829
830
831
832
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, count_background_gpart, nr_cells,
               MPI_LONG_LONG_INT, MPI_BOR, 0, MPI_COMM_WORLD);
  } else {
    MPI_Reduce(count_background_gpart, NULL, nr_cells, MPI_LONG_LONG_INT,
               MPI_BOR, 0, MPI_COMM_WORLD);
  }
833
834
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, count_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
Matthieu Schaller's avatar
Matthieu Schaller committed
835
               0, MPI_COMM_WORLD);
836
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
837
838
    MPI_Reduce(count_spart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
839
  }
840
841
842
843
844
845
846
847
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, count_bpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
               0, MPI_COMM_WORLD);
  } else {
    MPI_Reduce(count_bpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
  }

848
849
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, offset_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
Matthieu Schaller's avatar
Matthieu Schaller committed
850
               0, MPI_COMM_WORLD);
851
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
852
853
    MPI_Reduce(offset_part, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
854
855
856
  }
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, offset_gpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
Matthieu Schaller's avatar
Matthieu Schaller committed
857
               0, MPI_COMM_WORLD);
858
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
859
860
    MPI_Reduce(offset_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
861
  }
862
863
864
865
866
867
868
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, offset_background_gpart, nr_cells,
               MPI_LONG_LONG_INT, MPI_BOR, 0, MPI_COMM_WORLD);
  } else {
    MPI_Reduce(offset_background_gpart, NULL, nr_cells, MPI_LONG_LONG_INT,
               MPI_BOR, 0, MPI_COMM_WORLD);
  }
869
870
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, offset_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
Matthieu Schaller's avatar
Matthieu Schaller committed
871
               0, MPI_COMM_WORLD);
872
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
873
874
    MPI_Reduce(offset_spart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
875
  }
876
877
878
879
880
881
882
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, offset_bpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
               0, MPI_COMM_WORLD);
  } else {
    MPI_Reduce(offset_bpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
  }
883

884
  /* For the centres we use a sum as MPI does not like bit-wise operations
885
886
     on floating point numbers */
  if (nodeID == 0) {
Matthieu Schaller's avatar
Matthieu Schaller committed
887
888
    MPI_Reduce(MPI_IN_PLACE, centres, 3 * nr_cells, MPI_DOUBLE, MPI_SUM, 0,
               MPI_COMM_WORLD);
889
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
890
891
    MPI_Reduce(centres, NULL, 3 * nr_cells, MPI_DOUBLE, MPI_SUM, 0,
               MPI_COMM_WORLD);
892
  }
893
894
#endif

895
  /* Only rank 0 actually writes */
896
  if (nodeID == 0) {
897

898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
    /* 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
916
917
918
919
920
    /* 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);
921
    io_write_attribute(h_subgrp, "size", DOUBLE, cell_width, 3);
Matthieu Schaller's avatar
Matthieu Schaller committed
922
923
924
925
    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
926
    hsize_t shape[2] = {(hsize_t)nr_cells, 3};
Matthieu Schaller's avatar
Matthieu Schaller committed
927
928
929
    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);
930
931
    if (h_err < 0)
      error("Error while changing shape of gas offsets data space.");
Matthieu Schaller's avatar
Matthieu Schaller committed
932
933
    hid_t h_data = H5Dcreate(h_grp, "Centres", io_hdf5_type(DOUBLE), h_space,
                             H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
934
    if (h_data < 0) error("Error while creating dataspace for gas offsets.");
Matthieu Schaller's avatar
Matthieu Schaller committed
935
936
937
    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.");
938
939
940
    H5Dclose(h_data);
    H5Sclose(h_space);

Matthieu Schaller's avatar
Matthieu Schaller committed
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
    /* Group containing the offsets for each particle type */
    h_subgrp =
        H5Gcreate(h_grp, "Offsets", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    if (h_subgrp < 0) error("Error while creating offsets sub-group");

    if (global_counts[swift_type_gas] > 0) {

      shape[0] = nr_cells;
      shape[1] = 1;
      h_space = H5Screate(H5S_SIMPLE);
      if (h_space < 0) error("Error while creating data space for gas offsets");
      h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
      if (h_err < 0)
        error("Error while changing shape of gas offsets data space.");
      h_data = H5Dcreate(h_subgrp, "PartType0", io_hdf5_type(LONGLONG), h_space,
                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
      if (h_data < 0) error("Error while creating dataspace for gas offsets.");
      h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
                       H5P_DEFAULT, offset_part);
      if (h_err < 0) error("Error while writing gas offsets.");
      H5Dclose(h_data);
      H5Sclose(h_space);
    }
964

Matthieu Schaller's avatar
Matthieu Schaller committed
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
    if (global_counts[swift_type_dark_matter] > 0) {

      shape[0] = nr_cells;
      shape[1] = 1;
      h_space = H5Screate(H5S_SIMPLE);
      if (h_space < 0) error("Error while creating data space for DM offsets");
      h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
      if (h_err < 0)
        error("Error while changing shape of DM offsets data space.");
      h_data = H5Dcreate(h_subgrp, "PartType1", io_hdf5_type(LONGLONG), h_space,
                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
      if (h_data < 0) error("Error while creating dataspace for DM offsets.");
      h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
                       H5P_DEFAULT, offset_gpart);
      if (h_err < 0) error("Error while writing DM offsets.");
      H5Dclose(h_data);
      H5Sclose(h_space);
    }
983

984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
    if (global_counts[swift_type_dark_matter_background] > 0) {

      shape[0] = nr_cells;
      shape[1] = 1;
      h_space = H5Screate(H5S_SIMPLE);
      if (h_space < 0)
        error("Error while creating data space for background DM offsets");
      h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
      if (h_err < 0)
        error(
            "Error while changing shape of background DM offsets data space.");
      h_data = H5Dcreate(h_subgrp, "PartType2", io_hdf5_type(LONGLONG), h_space,
                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
      if (h_data < 0)
        error("Error while creating dataspace for background DM offsets.");
      h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
                       H5P_DEFAULT, offset_background_gpart);