common_io.c 50.9 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
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],
394
395
396
397
398
                           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]};
399
400
401

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

404
405
  /* Count of particles in each cell */
  long long *count_part = NULL, *count_gpart = NULL, *count_spart = NULL;
406
407
408
  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));
409
410
411

  /* Global offsets of particles in each cell */
  long long *offset_part = NULL, *offset_gpart = NULL, *offset_spart = NULL;
412
413
414
  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));
415
416
417
418
419
420
421

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

  /* Collect the cell information of *local* cells */
422
423
424
  long long local_offset_part = 0;
  long long local_offset_gpart = 0;
  long long local_offset_spart = 0;
425
426
427
428
429
  for (int i = 0; i < nr_cells; ++i) {

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

      /* Centre of each cell */
430
431
432
      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;
433
434
435
436
437
438
439
440
441
442
443
444

      /* Count real particles that will be written */
      count_part[i] = cells_top[i].hydro.count - cells_top[i].hydro.inhibited;
      count_gpart[i] = cells_top[i].grav.count - cells_top[i].grav.inhibited;
      count_spart[i] = cells_top[i].stars.count - cells_top[i].stars.inhibited;

      /* Only count DM gpart (gpart without friends) */
      count_gpart[i] -= count_part[i];
      count_gpart[i] -= count_spart[i];

      /* Offsets including the global offset of all particles on this MPI rank
       */
445
      offset_part[i] = local_offset_part + global_offsets[swift_type_gas];
Matthieu Schaller's avatar
Matthieu Schaller committed
446
447
      offset_gpart[i] =
          local_offset_gpart + global_offsets[swift_type_dark_matter];
448
      offset_spart[i] = local_offset_spart + global_offsets[swift_type_stars];
449

450
451
452
453
      local_offset_part += count_part[i];
      local_offset_gpart += count_gpart[i];
      local_offset_spart += count_spart[i];

454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
    } 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;

      offset_part[i] = 0;
      offset_gpart[i] = 0;
      offset_spart[i] = 0;
    }
  }

472
#ifdef WITH_MPI
Matthieu Schaller's avatar
Matthieu Schaller committed
473
474
  /* 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. */
475
476
  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
477
               0, MPI_COMM_WORLD);
478
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
479
480
    MPI_Reduce(count_part, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
481
482
483
  }
  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
484
               0, MPI_COMM_WORLD);
485
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
486
487
    MPI_Reduce(count_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
488
489
490
  }
  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
491
               0, MPI_COMM_WORLD);
492
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
493
494
    MPI_Reduce(count_spart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
495
496
497
  }
  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
498
               0, MPI_COMM_WORLD);
499
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
500
501
    MPI_Reduce(offset_part, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
502
503
504
  }
  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
505
               0, MPI_COMM_WORLD);
506
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
507
508
    MPI_Reduce(offset_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
509
510
511
  }
  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
512
               0, MPI_COMM_WORLD);
513
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
514
515
    MPI_Reduce(offset_spart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
               MPI_COMM_WORLD);
516
517
  }

518
  /* For the centres we use a sum as MPI does not like bit-wise operations
519
520
     on floating point numbers */
  if (nodeID == 0) {
Matthieu Schaller's avatar
Matthieu Schaller committed
521
522
    MPI_Reduce(MPI_IN_PLACE, centres, 3 * nr_cells, MPI_DOUBLE, MPI_SUM, 0,
               MPI_COMM_WORLD);
523
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
524
525
    MPI_Reduce(centres, NULL, 3 * nr_cells, MPI_DOUBLE, MPI_SUM, 0,
               MPI_COMM_WORLD);
526
  }
527
528
#endif

529
  /* Only rank 0 actually writes */
530
  if (nodeID == 0) {
531

532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
    /* 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
550
551
552
553
554
    /* 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);
555
    io_write_attribute(h_subgrp, "size", DOUBLE, cell_width, 3);
Matthieu Schaller's avatar
Matthieu Schaller committed
556
557
558
559
    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
560
    hsize_t shape[2] = {(hsize_t)nr_cells, 3};
Matthieu Schaller's avatar
Matthieu Schaller committed
561
562
563
    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);
564
565
    if (h_err < 0)
      error("Error while changing shape of gas offsets data space.");
Matthieu Schaller's avatar
Matthieu Schaller committed
566
567
    hid_t h_data = H5Dcreate(h_grp, "Centres", io_hdf5_type(DOUBLE), h_space,
                             H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
568
    if (h_data < 0) error("Error while creating dataspace for gas offsets.");
Matthieu Schaller's avatar
Matthieu Schaller committed
569
570
571
    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.");
572
573
574
    H5Dclose(h_data);
    H5Sclose(h_space);

Matthieu Schaller's avatar
Matthieu Schaller committed
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
    /* 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);
    }
598

Matthieu Schaller's avatar
Matthieu Schaller committed
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
    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);
    }
617

Matthieu Schaller's avatar
Matthieu Schaller committed
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
    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);
    }
637

Matthieu Schaller's avatar
Matthieu Schaller committed
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
    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);
    }
663

Matthieu Schaller's avatar
Matthieu Schaller committed
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
    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);
    }
682

Matthieu Schaller's avatar
Matthieu Schaller committed
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
    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);
    }
702

Matthieu Schaller's avatar
Matthieu Schaller committed
703
    H5Gclose(h_subgrp);
704
705
706
707
708
709
710
711
712
713
  }

  /* Free everything we allocated */
  free(centres);
  free(count_part);
  free(count_gpart);
  free(count_spart);
  free(offset_part);
  free(offset_gpart);
  free(offset_spart);
714
715
}

716
717
#endif /* HAVE_HDF5 */

718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
/**
 * @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;
  }
}

748
749
750
/**
 * @brief Mapper function to copy #part or #gpart fields into a buffer.
 */
751
void io_copy_mapper(void* restrict temp, int N, void* restrict extra_data) {
752

753
  const struct io_props props = *((const struct io_props*)extra_data);
754
755
756
757
  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
758
  char* restrict temp_c = (char*)temp;
759
760
  const ptrdiff_t delta = (temp_c - props.start_temp_c) / copySize;

761
762
763
  for (int k = 0; k < N; k++) {
    memcpy(&temp_c[k * copySize], props.field + (delta + k) * props.partSize,
           copySize);
764
765
766
767
  }
}

/**
768
769
 * @brief Mapper function to copy #part into a buffer of floats using a
 * conversion function.
770
 */
771
772
void io_convert_part_f_mapper(void* restrict temp, int N,
                              void* restrict extra_data) {
773

774
775
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct part* restrict parts = props.parts;
776
  const struct xpart* restrict xparts = props.xparts;
777
778
779
780
  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
781
  float* restrict temp_f = (float*)temp;
782
  const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
783

784
  for (int i = 0; i < N; i++)
785
786
    props.convert_part_f(e, parts + delta + i, xparts + delta + i,
                         &temp_f[i * dim]);
787
788
789
}

/**
790
791
 * @brief Mapper function to copy #part into a buffer of doubles using a
 * conversion function.
792
 */
793
794
void io_convert_part_d_mapper(void* restrict temp, int N,
                              void* restrict extra_data) {
795

796
797
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct part* restrict parts = props.parts;
798
  const struct xpart* restrict xparts = props.xparts;
799
800
801
802
  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
803
  double* restrict temp_d = (double*)temp;
804
  const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
805

806
  for (int i = 0; i < N; i++)
807
808
    props.convert_part_d(e, parts + delta + i, xparts + delta + i,
                         &temp_d[i * dim]);
809
810
}

811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
/**
 * @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]);
}

833
/**
834
835
 * @brief Mapper function to copy #gpart into a buffer of floats using a
 * conversion function.
836
 */
837
838
void io_convert_gpart_f_mapper(void* restrict temp, int N,
                               void* restrict extra_data) {
839

840
841
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct gpart* restrict gparts = props.gparts;
842
843
844
845
  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
846
  float* restrict temp_f = (float*)temp;
847
  const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
848

849
850
851
852
853
  for (int i = 0; i < N; i++)
    props.convert_gpart_f(e, gparts + delta + i, &temp_f[i * dim]);
}

/**
854
855
 * @brief Mapper function to copy #gpart into a buffer of doubles using a
 * conversion function.
856
 */
857
858
void io_convert_gpart_d_mapper(void* restrict temp, int N,
                               void* restrict extra_data) {
859

860
861
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct gpart* restrict gparts = props.gparts;
862
863
864
865
  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
866
  double* restrict temp_d = (double*)temp;
867
  const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
868

869
870
871
872
  for (int i = 0; i < N; i++)
    props.convert_gpart_d(e, gparts + delta + i, &temp_d[i * dim]);
}

873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
/**
 * @brief Mapper function to copy #gpart into a buffer of doubles using a
 * conversion function.
 */
void io_convert_gpart_l_mapper(void* restrict temp, int N,
                               void* restrict extra_data) {

  const struct io_props props = *((const struct io_props*)extra_data);
  const struct gpart* restrict gparts = props.gparts;
  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_gpart_l(e, gparts + delta + i, &temp_l[i * dim]);
}

893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
/**
 * @brief Mapper function to copy #spart into a buffer of floats using a
 * conversion function.
 */
void io_convert_spart_f_mapper(void* restrict temp, int N,
                               void* restrict extra_data) {

  const struct io_props props = *((const struct io_props*)extra_data);
  const struct spart* restrict sparts = props.sparts;
  const struct engine* e = props.e;
  const size_t dim = props.dimension;

  /* How far are we with this chunk? */
  float* restrict temp_f = (float*)temp;
  const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;

  for (int i = 0; i < N; i++)
    props.convert_spart_f(e, sparts + delta + i, &temp_f[i * dim]);
}

/**
 * @brief Mapper function to copy #spart into a buffer of doubles using a
 * conversion function.
 */
void io_convert_spart_d_mapper(void* restrict temp, int N,
                               void* restrict extra_data) {

  const struct io_props props = *((const struct io_props*)extra_data);
  const struct spart* restrict sparts = props.sparts;
  const struct engine* e = props.e;
  const size_t dim = props.dimension;

  /* How far are we with this chunk? */
  double* restrict temp_d = (double*)temp;
  const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;

  for (int i = 0; i < N; i++)
    props.convert_spart_d(e, sparts + delta + i, &temp_d[i * dim]);
}

933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
/**
 * @brief Mapper function to copy #spart into a buffer of doubles using a
 * conversion function.
 */
void io_convert_spart_l_mapper(void* restrict temp, int N,
                               void* restrict extra_data) {

  const struct io_props props = *((const struct io_props*)extra_data);
  const struct spart* restrict sparts = props.sparts;
  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_spart_l(e, sparts + delta + i, &temp_l[i * dim]);
}

953
954
955
956
957
/**
 * @brief Copy the particle data into a temporary buffer ready for i/o.
 *
 * @param temp The buffer to be filled. Must be allocated and aligned properly.
 * @param e The #engine.
958
959
 * @param props The #io_props corresponding to the particle field we are
 * copying.
960
961
962
963
964
 * @param N The number of particles to copy
 * @param internal_units The system of units used internally.
 * @param snapshot_units The system of units used for the snapshots.
 */
void io_copy_temp_buffer(void* temp, const struct engine* e,
965
                         struct io_props props, size_t N,
966
967
968
969
970
971
972
973
974
975
                         const struct unit_system* internal_units,
                         const struct unit_system* snapshot_units) {

  const size_t typeSize = io_sizeof_type(props.type);
  const size_t copySize = typeSize * props.dimension;
  const size_t num_elements = N * props.dimension;

  /* Copy particle data to temporary buffer */
  if (props.conversion == 0) { /* No conversion */

976
    /* Prepare some parameters */
Matthieu Schaller's avatar
Matthieu Schaller committed
977
    char* temp_c = (char*)temp;
978
    props.start_temp_c = temp_c;
979

980
    /* Copy the whole thing into a buffer */
981
982
983
    threadpool_map((struct threadpool*)&e->threadpool, io_copy_mapper, temp_c,
                   N, copySize, 0, (void*)&props);

984
985
986
987
  } else { /* Converting particle to data */

    if (props.convert_part_f != NULL) {

988
      /* Prepare some parameters */
Matthieu Schaller's avatar
Matthieu Schaller committed
989
990
      float* temp_f = (float*)temp;
      props.start_temp_f = (float*)temp;
991
      props.e = e;
992

993
      /* Copy the whole thing into a buffer */
994
995
996
      threadpool_map((struct threadpool*)&e->threadpool,
                     io_convert_part_f_mapper, temp_f, N, copySize, 0,
                     (void*)&props);
997
998

    } else if (props.convert_part_d != NULL) {
999

1000
      /* Prepare some parameters */
For faster browsing, not all history is shown. View entire blame