common_io.c 47 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
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],
                           const long long global_offsets[swift_type_count]) {

  /* Temporary memory for the cell-by-cell information */
  double* centres = NULL;
  if (posix_memalign((void**)&centres, IO_BUFFER_ALIGNMENT,
                     nr_cells * 3 * sizeof(double)) != 0)
    error("Error allocating temporary buffer for centres");

  message("global offsets=%lld", global_offsets[0]);

404
405
406
  int* node = NULL;
  node = malloc(nr_cells * sizeof(int));

407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
  /* Count of particles in each cell */
  long long *count_part = NULL, *count_gpart = NULL, *count_spart = NULL;
  count_part = malloc(nr_cells * sizeof(long long));
  count_gpart = malloc(nr_cells * sizeof(long long));
  count_spart = malloc(nr_cells * sizeof(long long));

  /* Global offsets of particles in each cell */
  long long *offset_part = NULL, *offset_gpart = NULL, *offset_spart = NULL;
  offset_part = malloc(nr_cells * sizeof(long long));
  offset_gpart = malloc(nr_cells * sizeof(long long));
  offset_spart = malloc(nr_cells * sizeof(long long));

  /* 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 */
  int count_local_cells = 0;
426
427
428
  long long local_offset_part = 0;
  long long local_offset_gpart = 0;
  long long local_offset_spart = 0;
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
  for (int i = 0; i < nr_cells; ++i) {

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

      /* Centre of each cell */
      centres[i * 3 + 0] = cells_top[i].loc[0] + width[0] * 0.5;
      centres[i * 3 + 1] = cells_top[i].loc[1] + width[1] * 0.5;
      centres[i * 3 + 2] = cells_top[i].loc[2] + width[2] * 0.5;

      /* 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
       */
449
450
451
      offset_part[i] = local_offset_part + global_offsets[swift_type_gas];
      offset_gpart[i] = local_offset_gpart + global_offsets[swift_type_dark_matter];
      offset_spart[i] = local_offset_spart + global_offsets[swift_type_stars];
452
453

      ++count_local_cells;
454
455
456
457
458

      local_offset_part += count_part[i];
      local_offset_gpart += count_gpart[i];
      local_offset_spart += count_spart[i];

459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
    } 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;
    }
475
476

    node[i] = cells_top[i].nodeID;
477
478
  }

479
#ifdef WITH_MPI
480
481
    /* Now, reduce all the arrays. Note that we use a bit-by-bit OR here. This
       is safe as we made sure only local cells have non-zero values. */
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, count_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
	       0, MPI_COMM_WORLD);
  } else {
    MPI_Reduce(count_part, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
	       0, MPI_COMM_WORLD);
  }
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, count_gpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
	       0, MPI_COMM_WORLD);
  } else {
    MPI_Reduce(count_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
	       0, MPI_COMM_WORLD);
  }
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, count_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
	       0, MPI_COMM_WORLD);
  } else {
    MPI_Reduce(count_spart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
	       0, MPI_COMM_WORLD);
  }
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, offset_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
	       0, MPI_COMM_WORLD);
  } else {
    MPI_Reduce(offset_part, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
	       0, MPI_COMM_WORLD);
  }
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, offset_gpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
	       0, MPI_COMM_WORLD);
  } else {
    MPI_Reduce(offset_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
	       0, MPI_COMM_WORLD);
  }
  if (nodeID == 0) {
    MPI_Reduce(MPI_IN_PLACE, offset_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
	       0, MPI_COMM_WORLD);
  } else {
    MPI_Reduce(offset_spart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
	       0, MPI_COMM_WORLD);
  }


  /* For the centres we use a sum as MPI does not like bitwise operations
     on floating point numbers */
  if (nodeID == 0) {
     MPI_Reduce(MPI_IN_PLACE, centres, 3 * nr_cells, MPI_DOUBLE, MPI_SUM,
                0, MPI_COMM_WORLD);
  } else {
    MPI_Reduce(centres, NULL, 3 * nr_cells, MPI_DOUBLE, MPI_SUM,
                0, MPI_COMM_WORLD);
  }
535
536
#endif

537
538
  /* Only rank 0 actually writes */
  if (nodeID != 0) return;
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564

  /* 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);
  io_write_attribute(h_subgrp, "size", DOUBLE, width, 3);
  io_write_attribute(h_subgrp, "dimension", INT, cdim, 3);
  H5Gclose(h_subgrp);

  /* Write the centres to the group */
  hsize_t shape[2] = {nr_cells, 3};
  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);
  if (h_err < 0) error("Error while changing shape of gas offsets data space.");
  hid_t h_data = H5Dcreate(h_grp, "Centres", io_hdf5_type(DOUBLE), 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(DOUBLE), h_space, H5S_ALL, H5P_DEFAULT,
                   centres);
  if (h_err < 0) error("Error while writing centres.");
  H5Dclose(h_data);
  H5Sclose(h_space);
  free(centres);

565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
  /* Write the nodeIDs to the group */
  shape[0] = nr_cells;
  shape[1] = 1;
  h_space = H5Screate(H5S_SIMPLE);
  if (h_space < 0) error("Error while creating data space for cell centres");
  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_grp, "Nodes", io_hdf5_type(INT), 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(INT), h_space, H5S_ALL, H5P_DEFAULT,
                   node);
  if (h_err < 0) error("Error while writing centres.");
  H5Dclose(h_data);
  H5Sclose(h_space);
  free(node);


583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
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
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
  /* 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);
  }

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

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

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

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

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

  H5Gclose(h_subgrp);
}

709
710
#endif /* HAVE_HDF5 */

711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
/**
 * @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;
  }
}

741
742
743
/**
 * @brief Mapper function to copy #part or #gpart fields into a buffer.
 */
744
void io_copy_mapper(void* restrict temp, int N, void* restrict extra_data) {
745

746
  const struct io_props props = *((const struct io_props*)extra_data);
747
748
749
750
  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
751
  char* restrict temp_c = (char*)temp;
752
753
  const ptrdiff_t delta = (temp_c - props.start_temp_c) / copySize;

754
755
756
  for (int k = 0; k < N; k++) {
    memcpy(&temp_c[k * copySize], props.field + (delta + k) * props.partSize,
           copySize);
757
758
759
760
  }
}

/**
761
762
 * @brief Mapper function to copy #part into a buffer of floats using a
 * conversion function.
763
 */
764
765
void io_convert_part_f_mapper(void* restrict temp, int N,
                              void* restrict extra_data) {
766

767
768
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct part* restrict parts = props.parts;
769
  const struct xpart* restrict xparts = props.xparts;
770
771
772
773
  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
774
  float* restrict temp_f = (float*)temp;
775
  const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
776

777
  for (int i = 0; i < N; i++)
778
779
    props.convert_part_f(e, parts + delta + i, xparts + delta + i,
                         &temp_f[i * dim]);
780
781
782
}

/**
783
784
 * @brief Mapper function to copy #part into a buffer of doubles using a
 * conversion function.
785
 */
786
787
void io_convert_part_d_mapper(void* restrict temp, int N,
                              void* restrict extra_data) {
788

789
790
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct part* restrict parts = props.parts;
791
  const struct xpart* restrict xparts = props.xparts;
792
793
794
795
  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
796
  double* restrict temp_d = (double*)temp;
797
  const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
798

799
  for (int i = 0; i < N; i++)
800
801
    props.convert_part_d(e, parts + delta + i, xparts + delta + i,
                         &temp_d[i * dim]);
802
803
804
}

/**
805
806
 * @brief Mapper function to copy #gpart into a buffer of floats using a
 * conversion function.
807
 */
808
809
void io_convert_gpart_f_mapper(void* restrict temp, int N,
                               void* restrict extra_data) {
810

811
812
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct gpart* restrict gparts = props.gparts;
813
814
815
816
  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
817
  float* restrict temp_f = (float*)temp;
818
  const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
819

820
821
822
823
824
  for (int i = 0; i < N; i++)
    props.convert_gpart_f(e, gparts + delta + i, &temp_f[i * dim]);
}

/**
825
826
 * @brief Mapper function to copy #gpart into a buffer of doubles using a
 * conversion function.
827
 */
828
829
void io_convert_gpart_d_mapper(void* restrict temp, int N,
                               void* restrict extra_data) {
830

831
832
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct gpart* restrict gparts = props.gparts;
833
834
835
836
  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
837
  double* restrict temp_d = (double*)temp;
838
  const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
839

840
841
842
843
  for (int i = 0; i < N; i++)
    props.convert_gpart_d(e, gparts + delta + i, &temp_d[i * dim]);
}

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
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
/**
 * @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]);
}

884
885
886
887
888
/**
 * @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.
889
890
 * @param props The #io_props corresponding to the particle field we are
 * copying.
891
892
893
894
895
 * @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,
896
                         struct io_props props, size_t N,
897
898
899
900
901
902
903
904
905
906
                         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 */

907
    /* Prepare some parameters */
Matthieu Schaller's avatar
Matthieu Schaller committed
908
    char* temp_c = (char*)temp;
909
    props.start_temp_c = temp_c;
910

911
    /* Copy the whole thing into a buffer */
912
913
914
    threadpool_map((struct threadpool*)&e->threadpool, io_copy_mapper, temp_c,
                   N, copySize, 0, (void*)&props);

915
916
917
918
  } else { /* Converting particle to data */

    if (props.convert_part_f != NULL) {

919
      /* Prepare some parameters */
Matthieu Schaller's avatar
Matthieu Schaller committed
920
921
      float* temp_f = (float*)temp;
      props.start_temp_f = (float*)temp;
922
      props.e = e;
923

924
      /* Copy the whole thing into a buffer */
925
926
927
      threadpool_map((struct threadpool*)&e->threadpool,
                     io_convert_part_f_mapper, temp_f, N, copySize, 0,
                     (void*)&props);
928
929

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

931
      /* Prepare some parameters */
Matthieu Schaller's avatar
Matthieu Schaller committed
932
933
      double* temp_d = (double*)temp;
      props.start_temp_d = (double*)temp;
934
      props.e = e;
935

936
      /* Copy the whole thing into a buffer */
937
938
939
      threadpool_map((struct threadpool*)&e->threadpool,
                     io_convert_part_d_mapper, temp_d, N, copySize, 0,
                     (void*)&props);
940
941
942

    } else if (props.convert_gpart_f != NULL) {

943
      /* Prepare some parameters */
Matthieu Schaller's avatar
Matthieu Schaller committed
944
945
      float* temp_f = (float*)temp;
      props.start_temp_f = (float*)temp;
946
      props.e = e;
947

948
      /* Copy the whole thing into a buffer */
949
950
951
      threadpool_map((struct threadpool*)&e->threadpool,
                     io_convert_gpart_f_mapper, temp_f, N, copySize, 0,
                     (void*)&props);
952
953
954

    } else if (props.convert_gpart_d != NULL) {

955
      /* Prepare some parameters */
Matthieu Schaller's avatar
Matthieu Schaller committed
956
957
      double* temp_d = (double*)temp;
      props.start_temp_d = (double*)temp;
958
      props.e = e;
959

960
      /* Copy the whole thing into a buffer */
961
962
963
      threadpool_map((struct threadpool*)&e->threadpool,
                     io_convert_gpart_d_mapper, temp_d, N, copySize, 0,
                     (void*)&props);
964

965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
    } else if (props.convert_spart_f != NULL) {

      /* Prepare some parameters */
      float* temp_f = (float*)temp;
      props.start_temp_f = (float*)temp;
      props.e = e;

      /* Copy the whole thing into a buffer */
      threadpool_map((struct threadpool*)&e->threadpool,
                     io_convert_spart_f_mapper, temp_f, N, copySize, 0,
                     (void*)&props);

    } else if (props.convert_spart_d != NULL) {

      /* Prepare some parameters */
      double* temp_d = (double*)temp;
      props.start_temp_d = (double*)temp;
      props.e = e;

      /* Copy the whole thing into a buffer */
      threadpool_map((struct threadpool*)&e->threadpool,
                     io_convert_spart_d_mapper, temp_d, N, copySize, 0,
                     (void*)&props);

989
990
991
992
993
994
995
996
997
998
999
1000
    } else {
      error("Missing conversion function");
    }
  }

  /* Unit conversion if necessary */
  const double factor =
      units_conversion_factor(internal_units, snapshot_units, props.units);
  if (factor != 1.) {

    /* message("Converting ! factor=%e", factor); */

For faster browsing, not all history is shown. View entire blame