common_io.c 62.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
 ******************************************************************************/

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

24
25
26
27
/* This object's header. */
#include "common_io.h"

/* Local includes. */
28
#include "chemistry_io.h"
29
#include "engine.h"
30
#include "error.h"
31
#include "gravity_io.h"
32
#include "hydro.h"
33
#include "hydro_io.h"
34
#include "io_properties.h"
35
#include "kernel_hydro.h"
36
#include "part.h"
lhausamm's avatar
lhausamm committed
37
#include "part_type.h"
38
#include "stars_io.h"
39
#include "threadpool.h"
40
#include "units.h"
41
#include "version.h"
42

43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
/* 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

59
/**
60
 * @brief Converts a C data type to the HDF5 equivalent.
61
62
 *
 * This function is a trivial wrapper around the HDF5 types but allows
63
64
 * to change the exact storage types matching the code types in a transparent
 *way.
65
 */
66
hid_t io_hdf5_type(enum IO_DATA_TYPE type) {
67

68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
  switch (type) {
    case INT:
      return H5T_NATIVE_INT;
    case UINT:
      return H5T_NATIVE_UINT;
    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:
86
      return H5T_NATIVE_CHAR;
87
88
89
90
    default:
      error("Unknown type");
      return 0;
  }
91
92
}

93
94
95
96
97
/**
 * @brief Return 1 if the type has double precision
 *
 * Returns an error if the type is not FLOAT or DOUBLE
 */
98
int io_is_double_precision(enum IO_DATA_TYPE type) {
Matthieu Schaller's avatar
Matthieu Schaller committed
99

100
101
102
103
104
105
106
107
108
109
110
  switch (type) {
    case FLOAT:
      return 0;
    case DOUBLE:
      return 1;
    default:
      error("Invalid type");
      return 0;
  }
}

111
112
113
114
115
/**
 * @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.
116
 * @param type The #IO_DATA_TYPE of the attribute.
117
118
119
120
 * @param data (output) The attribute read from the HDF5 group.
 *
 * Calls #error() if an error occurs.
 */
121
void io_read_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
122
                       void* data) {
123

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

127
128
  const hid_t h_err = H5Aread(h_attr, io_hdf5_type(type), data);
  if (h_err < 0) error("Error while reading attribute '%s'", name);
129
130
131
132
133
134
135
136
137

  H5Aclose(h_attr);
}

/**
 * @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.
138
 * @param type The #IO_DATA_TYPE of the attribute.
139
140
141
142
143
 * @param data The attribute to write.
 * @param num The number of elements to write
 *
 * Calls #error() if an error occurs.
 */
144
void io_write_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
145
                        const void* data, int num) {
146

147
148
  const hid_t h_space = H5Screate(H5S_SIMPLE);
  if (h_space < 0)
149
    error("Error while creating dataspace for attribute '%s'.", name);
150

151
152
153
  hsize_t dim[1] = {(hsize_t)num};
  const hid_t h_err = H5Sset_extent_simple(h_space, 1, dim, NULL);
  if (h_err < 0)
154
    error("Error while changing dataspace shape for attribute '%s'.", name);
155

156
157
158
  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);
159

160
161
  const hid_t h_err2 = H5Awrite(h_attr, io_hdf5_type(type), data);
  if (h_err2 < 0) error("Error while reading attribute '%s'.", name);
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176

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

180
181
  const hid_t h_space = H5Screate(H5S_SCALAR);
  if (h_space < 0)
182
    error("Error while creating dataspace for attribute '%s'.", name);
183

184
185
  const hid_t h_type = H5Tcopy(H5T_C_S1);
  if (h_type < 0) error("Error while copying datatype 'H5T_C_S1'.");
186

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

190
191
  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);
192

193
194
  const hid_t h_err2 = H5Awrite(h_attr, h_type, str);
  if (h_err2 < 0) error("Error while reading attribute '%s'.", name);
195
196
197
198
199
200
201
202

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

/**
 * @brief Writes a double value as an attribute
203
 * @param grp The group in which to write
204
205
206
 * @param name The name of the attribute
 * @param data The value to write
 */
207
208
void io_write_attribute_d(hid_t grp, const char* name, double data) {
  io_write_attribute(grp, name, DOUBLE, &data, 1);
209
210
211
212
}

/**
 * @brief Writes a float value as an attribute
213
 * @param grp The group in which to write
214
215
216
 * @param name The name of the attribute
 * @param data The value to write
 */
217
218
void io_write_attribute_f(hid_t grp, const char* name, float data) {
  io_write_attribute(grp, name, FLOAT, &data, 1);
219
220
221
222
}

/**
 * @brief Writes an int value as an attribute
223
 * @param grp The group in which to write
224
225
226
 * @param name The name of the attribute
 * @param data The value to write
 */
227
228
void io_write_attribute_i(hid_t grp, const char* name, int data) {
  io_write_attribute(grp, name, INT, &data, 1);
229
230
231
232
}

/**
 * @brief Writes a long value as an attribute
233
 * @param grp The group in which to write
234
235
236
 * @param name The name of the attribute
 * @param data The value to write
 */
237
238
void io_write_attribute_l(hid_t grp, const char* name, long data) {
  io_write_attribute(grp, name, LONG, &data, 1);
239
240
241
242
}

/**
 * @brief Writes a string value as an attribute
243
 * @param grp The group in which to write
244
245
246
 * @param name The name of the attribute
 * @param str The string to write
 */
247
248
void io_write_attribute_s(hid_t grp, const char* name, const char* str) {
  io_writeStringAttribute(grp, name, str, strlen(str));
249
250
}

251
252
/**
 * @brief Reads the Unit System from an IC file.
253
254
255
256
 *
 * If the 'Units' group does not exist in the ICs, we will use the internal
 * system of units.
 *
257
 * @param h_file The (opened) HDF5 file from which to read.
258
259
 * @param ic_units The unit_system to fill.
 * @param internal_units The internal system of units to copy if needed.
260
 * @param mpi_rank The MPI rank we are on.
261
 */
262
263
264
void io_read_unit_system(hid_t h_file, struct unit_system* ic_units,
                         const struct unit_system* internal_units,
                         int mpi_rank) {
265
266
267

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

269
  if (exists == 0) {
270
271

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

274
    units_copy(ic_units, internal_units);
275
276

    return;
277

278
  } else if (exists < 0) {
279
    error("Serious problem with 'Units' group in ICs. H5Lexists gives %d",
280
          exists);
281
  }
282

283
  if (mpi_rank == 0) message("Reading IC units from ICs.");
284
  hid_t h_grp = H5Gopen(h_file, "/Units", H5P_DEFAULT);
285
286

  /* Ok, Read the damn thing */
287
  io_read_attribute(h_grp, "Unit length in cgs (U_L)", DOUBLE,
288
                    &ic_units->UnitLength_in_cgs);
289
  io_read_attribute(h_grp, "Unit mass in cgs (U_M)", DOUBLE,
290
                    &ic_units->UnitMass_in_cgs);
291
  io_read_attribute(h_grp, "Unit time in cgs (U_t)", DOUBLE,
292
                    &ic_units->UnitTime_in_cgs);
293
  io_read_attribute(h_grp, "Unit current in cgs (U_I)", DOUBLE,
294
                    &ic_units->UnitCurrent_in_cgs);
295
  io_read_attribute(h_grp, "Unit temperature in cgs (U_T)", DOUBLE,
296
                    &ic_units->UnitTemperature_in_cgs);
297
298
299
300
301

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

302
303
304
/**
 * @brief Writes the current Unit System
 * @param h_file The (opened) HDF5 file in which to write
305
 * @param us The unit_system to dump
306
 * @param groupName The name of the HDF5 group to write to
307
 */
308
309
void io_write_unit_system(hid_t h_file, const struct unit_system* us,
                          const char* groupName) {
310

311
  const hid_t h_grpunit = H5Gcreate1(h_file, groupName, 0);
312
313
  if (h_grpunit < 0) error("Error while creating Unit System group");

314
315
316
317
318
319
320
321
322
323
  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));
324
325
326
327

  H5Gclose(h_grpunit);
}

328
329
330
331
/**
 * @brief Writes the code version to the file
 * @param h_file The (opened) HDF5 file in which to write
 */
332
void io_write_code_description(hid_t h_file) {
333

334
  const hid_t h_grpcode = H5Gcreate1(h_file, "/Code", 0);
335
336
  if (h_grpcode < 0) error("Error while creating code group");

337
  io_write_attribute_s(h_grpcode, "Code", "SWIFT");
338
339
340
341
342
343
344
345
346
347
  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());
348
  io_write_attribute_s(h_grpcode, "Thread barriers", thread_barrier_version());
349
  io_write_attribute_s(h_grpcode, "Allocators", allocator_version());
Matthieu Schaller's avatar
Matthieu Schaller committed
350
#ifdef HAVE_FFTW
351
  io_write_attribute_s(h_grpcode, "FFTW library version", fftw3_version());
Matthieu Schaller's avatar
Matthieu Schaller committed
352
#endif
353
354
355
#ifdef HAVE_LIBGSL
  io_write_attribute_s(h_grpcode, "GSL library version", libgsl_version());
#endif
356
#ifdef WITH_MPI
357
  io_write_attribute_s(h_grpcode, "MPI library", mpi_version());
358
#ifdef HAVE_METIS
359
  io_write_attribute_s(h_grpcode, "METIS library version", metis_version());
360
#endif
361
362
363
#ifdef HAVE_PARMETIS
  io_write_attribute_s(h_grpcode, "ParMETIS library version",
                       parmetis_version());
364
#endif
365
#else
366
  io_write_attribute_s(h_grpcode, "MPI library", "Non-MPI version of SWIFT");
367
#endif
368
369
370
  H5Gclose(h_grpcode);
}

371
372
373
374
375
376
377
378
379
380
/**
 * @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");

381
  for (int i = 1; i < engine_maxpolicy; ++i)
382
383
384
385
386
387
388
389
    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);
}

390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
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) {
    if (!(parts[i].time_bin != time_bin_inhibited) &&
        !(parts[i].time_bin != time_bin_not_created)) {
      ++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) {
    if (!(gparts[i].time_bin != time_bin_inhibited) &&
        !(gparts[i].time_bin != time_bin_not_created) &&
        (gparts[i].type == swift_type_dark_matter)) {
      ++count;
    }
  }
  return count;
}

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) {
    if (!(sparts[i].time_bin != time_bin_inhibited) &&
        !(sparts[i].time_bin != time_bin_not_created)) {
      ++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) {
    if (!(bparts[i].time_bin != time_bin_inhibited) &&
        !(bparts[i].time_bin != time_bin_not_created)) {
      ++count;
    }
  }
  return count;
}

443
444
445
446
void io_write_cell_offsets(hid_t h_grp, const int cdim[3],
                           const struct cell* cells_top, const int nr_cells,
                           const double width[3], const int nodeID,
                           const long long global_counts[swift_type_count],
447
448
449
450
451
                           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]};
452
453
454

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

457
  /* Count of particles in each cell */
458
459
  long long *count_part = NULL, *count_gpart = NULL, *count_spart = NULL,
            *count_bpart = NULL;
460
461
462
  count_part = (long long*)malloc(nr_cells * sizeof(long long));
  count_gpart = (long long*)malloc(nr_cells * sizeof(long long));
  count_spart = (long long*)malloc(nr_cells * sizeof(long long));
463
  count_bpart = (long long*)malloc(nr_cells * sizeof(long long));
464
465

  /* Global offsets of particles in each cell */
466
467
  long long *offset_part = NULL, *offset_gpart = NULL, *offset_spart = NULL,
            *offset_bpart = NULL;
468
469
470
  offset_part = (long long*)malloc(nr_cells * sizeof(long long));
  offset_gpart = (long long*)malloc(nr_cells * sizeof(long long));
  offset_spart = (long long*)malloc(nr_cells * sizeof(long long));
471
  offset_bpart = (long long*)malloc(nr_cells * sizeof(long long));
472
473
474
475
476

  /* Offsets of the 0^th element */
  offset_part[0] = 0;
  offset_gpart[0] = 0;
  offset_spart[0] = 0;
477
  offset_bpart[0] = 0;
478
479

  /* Collect the cell information of *local* cells */
480
481
482
  long long local_offset_part = 0;
  long long local_offset_gpart = 0;
  long long local_offset_spart = 0;
483
  long long local_offset_bpart = 0;
484
485
486
487
488
  for (int i = 0; i < nr_cells; ++i) {

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

      /* Centre of each cell */
489
490
491
      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;
492
493

      /* Count real particles that will be written */
494
495
496
497
      count_part[i] = cell_count_non_inhibited_gas(&cells_top[i]);
      count_gpart[i] = cell_count_non_inhibited_dark_matter(&cells_top[i]);
      count_spart[i] = cell_count_non_inhibited_stars(&cells_top[i]);
      count_bpart[i] = cell_count_non_inhibited_black_holes(&cells_top[i]);
498
499
500

      /* Offsets including the global offset of all particles on this MPI rank
       */
501
      offset_part[i] = local_offset_part + global_offsets[swift_type_gas];
Matthieu Schaller's avatar
Matthieu Schaller committed
502
503
      offset_gpart[i] =
          local_offset_gpart + global_offsets[swift_type_dark_matter];
504
      offset_spart[i] = local_offset_spart + global_offsets[swift_type_stars];
505
506
      offset_bpart[i] =
          local_offset_bpart + global_offsets[swift_type_black_hole];
507

508
509
510
      local_offset_part += count_part[i];
      local_offset_gpart += count_gpart[i];
      local_offset_spart += count_spart[i];
511
      local_offset_bpart += count_bpart[i];
512

513
514
515
516
517
518
519
520
521
522
523
    } 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;
      count_spart[i] = 0;
524
      count_bpart[i] = 0;
525
526
527
528

      offset_part[i] = 0;
      offset_gpart[i] = 0;
      offset_spart[i] = 0;
529
      offset_bpart[i] = 0;
530
531
532
    }
  }

533
#ifdef WITH_MPI
Matthieu Schaller's avatar
Matthieu Schaller committed
534
535
  /* 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. */
536
537
  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
538
               0, MPI_COMM_WORLD);
539
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
540
541
    MPI_Reduce(count_part, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
542
  }
543

544
545
  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
546
               0, MPI_COMM_WORLD);
547
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
548
549
    MPI_Reduce(count_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
550
551
552
  }
  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
553
               0, MPI_COMM_WORLD);
554
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
555
556
    MPI_Reduce(count_spart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
557
  }
558
559
560
561
562
563
564
565
  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);
  }

566
567
  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
568
               0, MPI_COMM_WORLD);
569
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
570
571
    MPI_Reduce(offset_part, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
572
573
574
  }
  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
575
               0, MPI_COMM_WORLD);
576
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
577
578
    MPI_Reduce(offset_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
579
580
581
  }
  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
582
               0, MPI_COMM_WORLD);
583
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
584
585
    MPI_Reduce(offset_spart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
586
  }
587
588
589
590
591
592
593
  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);
  }
594

595
  /* For the centres we use a sum as MPI does not like bit-wise operations
596
597
     on floating point numbers */
  if (nodeID == 0) {
Matthieu Schaller's avatar
Matthieu Schaller committed
598
599
    MPI_Reduce(MPI_IN_PLACE, centres, 3 * nr_cells, MPI_DOUBLE, MPI_SUM, 0,
               MPI_COMM_WORLD);
600
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
601
602
    MPI_Reduce(centres, NULL, 3 * nr_cells, MPI_DOUBLE, MPI_SUM, 0,
               MPI_COMM_WORLD);
603
  }
604
605
#endif

606
  /* Only rank 0 actually writes */
607
  if (nodeID == 0) {
608

609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
    /* 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
627
628
629
630
631
    /* 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);
632
    io_write_attribute(h_subgrp, "size", DOUBLE, cell_width, 3);
Matthieu Schaller's avatar
Matthieu Schaller committed
633
634
635
636
    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
637
    hsize_t shape[2] = {(hsize_t)nr_cells, 3};
Matthieu Schaller's avatar
Matthieu Schaller committed
638
639
640
    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);
641
642
    if (h_err < 0)
      error("Error while changing shape of gas offsets data space.");
Matthieu Schaller's avatar
Matthieu Schaller committed
643
644
    hid_t h_data = H5Dcreate(h_grp, "Centres", io_hdf5_type(DOUBLE), h_space,
                             H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
645
    if (h_data < 0) error("Error while creating dataspace for gas offsets.");
Matthieu Schaller's avatar
Matthieu Schaller committed
646
647
648
    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.");
649
650
651
    H5Dclose(h_data);
    H5Sclose(h_space);

Matthieu Schaller's avatar
Matthieu Schaller committed
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
    /* 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);
    }
675

Matthieu Schaller's avatar
Matthieu Schaller committed
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
    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);
    }
694

Matthieu Schaller's avatar
Matthieu Schaller committed
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
    if (global_counts[swift_type_stars] > 0) {

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

715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
    if (global_counts[swift_type_black_hole] > 0) {

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

Matthieu Schaller's avatar
Matthieu Schaller committed
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
    H5Gclose(h_subgrp);

    /* Group containing the counts for each particle type */
    h_subgrp =
        H5Gcreate(h_grp, "Counts", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    if (h_subgrp < 0) error("Error while creating counts 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 counts");
      h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
      if (h_err < 0)
        error("Error while changing shape of gas counts 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 counts.");
      h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
                       H5P_DEFAULT, count_part);
      if (h_err < 0) error("Error while writing gas counts.");
      H5Dclose(h_data);
      H5Sclose(h_space);
    }
761

Matthieu Schaller's avatar
Matthieu Schaller committed
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
    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 counts");
      h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
      if (h_err < 0)
        error("Error while changing shape of DM counts 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 counts.");
      h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
                       H5P_DEFAULT, count_gpart);
      if (h_err < 0) error("Error while writing DM counts.");
      H5Dclose(h_data);
      H5Sclose(h_space);
    }
780

Matthieu Schaller's avatar
Matthieu Schaller committed
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
    if (global_counts[swift_type_stars] > 0) {

      shape[0] = nr_cells;
      shape[1] = 1;
      h_space = H5Screate(H5S_SIMPLE);
      if (h_space < 0)
        error("Error while creating data space for stars counts");
      h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
      if (h_err < 0)
        error("Error while changing shape of stars counts data space.");
      h_data = H5Dcreate(h_subgrp, "PartType4", io_hdf5_type(LONGLONG), h_space,
                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
      if (h_data < 0) error("Error while creating dataspace for star counts.");
      h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
                       H5P_DEFAULT, count_spart);
      if (h_err < 0) error("Error while writing star counts.");
      H5Dclose(h_data);
      H5Sclose(h_space);
    }
800

801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
    if (global_counts[swift_type_black_hole] > 0) {

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

Matthieu Schaller's avatar
Matthieu Schaller committed
822
    H5Gclose(h_subgrp);
823
824
825
826
827
828
829
  }

  /* Free everything we allocated */
  free(centres);
  free(count_part);
  free(count_gpart);
  free(count_spart);
830
  free(count_bpart);
831
832
833
  free(offset_part);
  free(offset_gpart);
  free(offset_spart);
834
  free(offset_bpart);
835
836
}

837
838
#endif /* HAVE_HDF5 */

839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
/**
 * @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);
    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);
    default:
      error("Unknown type");
      return 0;
  }
}

869
870
871
/**
 * @brief Mapper function to copy #part or #gpart fields into a buffer.
 */
872
void io_copy_mapper(void* restrict temp, int N, void* restrict extra_data) {
873

874
  const struct io_props props = *((const struct io_props*)extra_data);
875
876
877
878
  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
879
  char* restrict temp_c = (char*)temp;
880
881
  const ptrdiff_t delta = (temp_c - props.start_temp_c) / copySize;

882
883
884
  for (int k = 0; k < N; k++) {
    memcpy(&temp_c[k * copySize], props.field + (delta + k) * props.partSize,
           copySize);
885
886
887
888
  }
}

/**
889
890
 * @brief Mapper function to copy #part into a buffer of floats using a
 * conversion function.
891
 */
892
893
void io_convert_part_f_mapper(void* restrict temp, int N,
                              void* restrict extra_data) {
894

895
896
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct part* restrict parts = props.parts;
897
  const struct xpart* restrict xparts = props.xparts;
898
899
900
901
  const struct engine* e = props.e;
  const size_t dim = props.dimension;

  /* How far are we with this chunk? */
Matthieu Schaller's avatar
Matthieu Schaller committed
902
  float* restrict temp_f = (float*)temp;
903
  const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
904

905
  for (int i = 0; i < N; i++)
906
907
    props.convert_part_f(e, parts + delta + i, xparts + delta + i,
                         &temp_f[i * dim]);
908
909
}

910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
/**
 * @brief Mapper function to copy #part into a buffer of ints using a
 * conversion function.
 */
void io_convert_part_i_mapper(void* restrict temp, int N,
                              void* restrict extra_data) {

  const struct io_props props = *((const struct io_props*)extra_data);
  const struct part* restrict parts = props.parts;
  const struct xpart* restrict xparts = props.xparts;
  const struct engine* e = props.e;
  const size_t dim = props.dimension;

  /* How far are we with this chunk? */
  int* restrict temp_i = (int*)temp;
  const ptrdiff_t delta = (temp_i - props.start_temp_i) / dim;

  for (int i = 0; i < N; i++)
    props.convert_part_i(e, parts + delta + i, xparts + delta + i,
                         &temp_i[i * dim]);
}

932
/**
933
934
 * @brief Mapper function to copy #part into a buffer of doubles using a
 * conversion function.
935
 */
936
937
void io_convert_part_d_mapper(void* restrict temp, int N,
                              void* restrict extra_data) {
938

939
940
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct part* restrict parts = props.parts;
941
  const struct xpart* restrict xparts = props.xparts;
942
943
944
945
  const struct engine* e = props.e;
  const size_t dim = props.dimension;

  /* How far are we with this chunk? */
Matthieu Schaller's avatar
Matthieu Schaller committed
946
  double* restrict temp_d = (double*)temp;
947
  const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
948

949
  for (int i = 0; i < N; i++)
950
951
    props.convert_part_d(e, parts + delta + i, xparts + delta + i,
                         &temp_d[i * dim]);
952
953
}

954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
/**
 * @brief Mapper function to copy #part into a buffer of doubles using a
 * conversion function.
 */
void io_convert_part_l_mapper(void* restrict temp, int N,
                              void* restrict extra_data) {

  const struct io_props props = *((const struct io_props*)extra_data);
  const struct part* restrict parts = props.parts;
  const struct xpart* restrict xparts = props.xparts;
  const struct engine* e = props.e;
  const size_t dim = props.dimension;

  /* How far are we with this chunk? */
  long long* restrict temp_l = (long long*)temp;
  const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;

  for (int i = 0; i < N; i++)
    props.convert_part_l(e, parts + delta + i, xparts + delta + i,
                         &temp_l[i * dim]);
}

976
/**
977
978
 * @brief Mapper function to copy #gpart into a buffer of floats using a
 * conversion function.
979
 */
980
981
void io_convert_gpart_f_mapper(void* restrict temp, int N,
                               void* restrict extra_data) {
982

983
984
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct gpart* restrict gparts = props.gparts;
985
986
987
988
  const struct engine* e = props.e;
  const size_t dim = props.dimension;

  /* How far are we with this chunk? */
Matthieu Schaller's avatar
Matthieu Schaller committed
989
  float* restrict temp_f = (float*)temp;
990
  const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
991

992
993
994
995
996
  for (int i = 0; i < N; i++)
    props.convert_gpart_f(e, gparts + delta + i, &temp_f[i * dim]);
}

/**
997
998
 * @brief Mapper function to copy #gpart into a buffer of doubles using a
 * conversion function.
999
 */
1000
void io_convert_gpart_d_mapper(void* restrict temp, int N,
For faster browsing, not all history is shown. View entire blame