common_io.c 29.8 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
145
void io_write_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
                        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
#else
362
  io_write_attribute_s(h_grpcode, "MPI library", "Non-MPI version of SWIFT");
363
#endif
364
365
366
  H5Gclose(h_grpcode);
}

367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
/**
 * @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");

  for (int i = 1; i <= engine_maxpolicy; ++i)
    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);
}

386
387
#endif /* HAVE_HDF5 */

388
389
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
/**
 * @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;
  }
}

418
419
420
/**
 * @brief Mapper function to copy #part or #gpart fields into a buffer.
 */
421
void io_copy_mapper(void* restrict temp, int N, void* restrict extra_data) {
422

423
  const struct io_props props = *((const struct io_props*)extra_data);
424
425
426
427
  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
428
  char* restrict temp_c = (char*)temp;
429
430
  const ptrdiff_t delta = (temp_c - props.start_temp_c) / copySize;

431
432
433
  for (int k = 0; k < N; k++) {
    memcpy(&temp_c[k * copySize], props.field + (delta + k) * props.partSize,
           copySize);
434
435
436
437
  }
}

/**
438
439
 * @brief Mapper function to copy #part into a buffer of floats using a
 * conversion function.
440
 */
441
442
void io_convert_part_f_mapper(void* restrict temp, int N,
                              void* restrict extra_data) {
443

444
445
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct part* restrict parts = props.parts;
446
  const struct xpart* restrict xparts = props.xparts;
447
448
449
450
  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
451
  float* restrict temp_f = (float*)temp;
452
  const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
453

454
  for (int i = 0; i < N; i++)
455
456
    props.convert_part_f(e, parts + delta + i, xparts + delta + i,
                         &temp_f[i * dim]);
457
458
459
}

/**
460
461
 * @brief Mapper function to copy #part into a buffer of doubles using a
 * conversion function.
462
 */
463
464
void io_convert_part_d_mapper(void* restrict temp, int N,
                              void* restrict extra_data) {
465

466
467
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct part* restrict parts = props.parts;
468
  const struct xpart* restrict xparts = props.xparts;
469
470
471
472
  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
473
  double* restrict temp_d = (double*)temp;
474
  const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
475

476
  for (int i = 0; i < N; i++)
477
478
    props.convert_part_d(e, parts + delta + i, xparts + delta + i,
                         &temp_d[i * dim]);
479
480
481
}

/**
482
483
 * @brief Mapper function to copy #gpart into a buffer of floats using a
 * conversion function.
484
 */
485
486
void io_convert_gpart_f_mapper(void* restrict temp, int N,
                               void* restrict extra_data) {
487

488
489
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct gpart* restrict gparts = props.gparts;
490
491
492
493
  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
494
  float* restrict temp_f = (float*)temp;
495
  const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
496

497
498
499
500
501
  for (int i = 0; i < N; i++)
    props.convert_gpart_f(e, gparts + delta + i, &temp_f[i * dim]);
}

/**
502
503
 * @brief Mapper function to copy #gpart into a buffer of doubles using a
 * conversion function.
504
 */
505
506
void io_convert_gpart_d_mapper(void* restrict temp, int N,
                               void* restrict extra_data) {
507

508
509
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct gpart* restrict gparts = props.gparts;
510
511
512
513
  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
514
  double* restrict temp_d = (double*)temp;
515
  const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
516

517
518
519
520
  for (int i = 0; i < N; i++)
    props.convert_gpart_d(e, gparts + delta + i, &temp_d[i * dim]);
}

521
522
523
524
525
/**
 * @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.
526
527
 * @param props The #io_props corresponding to the particle field we are
 * copying.
528
529
530
531
532
 * @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,
533
                         struct io_props props, size_t N,
534
535
536
537
538
539
540
541
542
543
                         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 */

544
    /* Prepare some parameters */
Matthieu Schaller's avatar
Matthieu Schaller committed
545
    char* temp_c = (char*)temp;
546
    props.start_temp_c = temp_c;
547

548
    /* Copy the whole thing into a buffer */
549
550
551
    threadpool_map((struct threadpool*)&e->threadpool, io_copy_mapper, temp_c,
                   N, copySize, 0, (void*)&props);

552
553
554
555
  } else { /* Converting particle to data */

    if (props.convert_part_f != NULL) {

556
      /* Prepare some parameters */
Matthieu Schaller's avatar
Matthieu Schaller committed
557
558
      float* temp_f = (float*)temp;
      props.start_temp_f = (float*)temp;
559
      props.e = e;
560

561
      /* Copy the whole thing into a buffer */
562
563
564
      threadpool_map((struct threadpool*)&e->threadpool,
                     io_convert_part_f_mapper, temp_f, N, copySize, 0,
                     (void*)&props);
565
566

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

568
      /* Prepare some parameters */
Matthieu Schaller's avatar
Matthieu Schaller committed
569
570
      double* temp_d = (double*)temp;
      props.start_temp_d = (double*)temp;
571
      props.e = e;
572

573
      /* Copy the whole thing into a buffer */
574
575
576
      threadpool_map((struct threadpool*)&e->threadpool,
                     io_convert_part_d_mapper, temp_d, N, copySize, 0,
                     (void*)&props);
577
578
579

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

580
      /* Prepare some parameters */
Matthieu Schaller's avatar
Matthieu Schaller committed
581
582
      float* temp_f = (float*)temp;
      props.start_temp_f = (float*)temp;
583
      props.e = e;
584

585
      /* Copy the whole thing into a buffer */
586
587
588
      threadpool_map((struct threadpool*)&e->threadpool,
                     io_convert_gpart_f_mapper, temp_f, N, copySize, 0,
                     (void*)&props);
589
590
591

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

592
      /* Prepare some parameters */
Matthieu Schaller's avatar
Matthieu Schaller committed
593
594
      double* temp_d = (double*)temp;
      props.start_temp_d = (double*)temp;
595
      props.e = e;
596

597
      /* Copy the whole thing into a buffer */
598
599
600
      threadpool_map((struct threadpool*)&e->threadpool,
                     io_convert_gpart_d_mapper, temp_d, N, copySize, 0,
                     (void*)&props);
601
602
603
604
605
606
607
608
609
610
611
612
613
614

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

    if (io_is_double_precision(props.type)) {
Matthieu Schaller's avatar
Matthieu Schaller committed
615
616
      swift_declare_aligned_ptr(double, temp_d, (double*)temp,
                                IO_BUFFER_ALIGNMENT);
617
618
      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= factor;
    } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
619
620
      swift_declare_aligned_ptr(float, temp_f, (float*)temp,
                                IO_BUFFER_ALIGNMENT);
621
622
623
624
625
      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor;
    }
  }
}

626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
void io_prepare_dm_gparts_mapper(void* restrict data, int Ndm, void* dummy) {

  struct gpart* restrict gparts = (struct gpart*)data;

  /* Let's give all these gparts a negative id */
  for (int i = 0; i < Ndm; ++i) {

    /* 0 or negative ids are not allowed */
    if (gparts[i].id_or_neg_offset <= 0)
      error("0 or negative ID for DM particle %i: ID=%lld", i,
            gparts[i].id_or_neg_offset);

    /* Set gpart type */
    gparts[i].type = swift_type_dark_matter;
  }
}

643
644
645
646
647
648
649
/**
 * @brief Prepare the DM particles (in gparts) read in for the addition of the
 * other particle types
 *
 * This function assumes that the DM particles are all at the start of the
 * gparts array
 *
650
 * @param tp The current #threadpool.
651
652
653
 * @param gparts The array of #gpart freshly read in.
 * @param Ndm The number of DM particles read in.
 */
654
655
void io_prepare_dm_gparts(struct threadpool* tp, struct gpart* const gparts,
                          size_t Ndm) {
656

657
658
  threadpool_map(tp, io_prepare_dm_gparts_mapper, gparts, Ndm,
                 sizeof(struct gpart), 0, NULL);
659
660
}

661
662
663
664
665
666
667
668
669
struct duplication_data {

  struct part* parts;
  struct gpart* gparts;
  struct spart* sparts;
  int Ndm;
  int Ngas;
  int Nstars;
};
670

671
672
673
674
675
676
677
678
679
680
void io_duplicate_hydro_gparts_mapper(void* restrict data, int Ngas,
                                      void* restrict extra_data) {

  struct duplication_data* temp = (struct duplication_data*)extra_data;
  const int Ndm = temp->Ndm;
  struct part* parts = (struct part*)data;
  const ptrdiff_t offset = parts - temp->parts;
  struct gpart* gparts = temp->gparts + offset;

  for (int i = 0; i < Ngas; ++i) {
681
682
683
684
685
686
687
688
689
690

    /* Duplicate the crucial information */
    gparts[i + Ndm].x[0] = parts[i].x[0];
    gparts[i + Ndm].x[1] = parts[i].x[1];
    gparts[i + Ndm].x[2] = parts[i].x[2];

    gparts[i + Ndm].v_full[0] = parts[i].v[0];
    gparts[i + Ndm].v_full[1] = parts[i].v[1];
    gparts[i + Ndm].v_full[2] = parts[i].v[2];

691
    gparts[i + Ndm].mass = hydro_get_mass(&parts[i]);
692

693
694
695
    /* Set gpart type */
    gparts[i + Ndm].type = swift_type_gas;

696
    /* Link the particles */
697
    gparts[i + Ndm].id_or_neg_offset = -(long long)(offset + i);
698
699
700
701
    parts[i].gpart = &gparts[i + Ndm];
  }
}

702
/**
703
 * @brief Copy every #part into the corresponding #gpart and link them.
704
 *
705
706
 * This function assumes that the DM particles are all at the start of the
 * gparts array and adds the hydro particles afterwards
707
 *
708
 * @param tp The current #threadpool.
709
710
711
712
713
 * @param parts The array of #part freshly read in.
 * @param gparts The array of #gpart freshly read in with all the DM particles
 * at the start
 * @param Ngas The number of gas particles read in.
 * @param Ndm The number of DM particles read in.
714
 */
715
716
717
718
719
720
721
void io_duplicate_hydro_gparts(struct threadpool* tp, struct part* const parts,
                               struct gpart* const gparts, size_t Ngas,
                               size_t Ndm) {
  struct duplication_data data;
  data.parts = parts;
  data.gparts = gparts;
  data.Ndm = Ndm;
722

723
724
725
726
727
728
729
730
731
732
733
734
735
736
  threadpool_map(tp, io_duplicate_hydro_gparts_mapper, parts, Ngas,
                 sizeof(struct part), 0, &data);
}

void io_duplicate_hydro_sparts_mapper(void* restrict data, int Nstars,
                                      void* restrict extra_data) {

  struct duplication_data* temp = (struct duplication_data*)extra_data;
  const int Ndm = temp->Ndm;
  struct spart* sparts = (struct spart*)data;
  const ptrdiff_t offset = sparts - temp->sparts;
  struct gpart* gparts = temp->gparts + offset;

  for (int i = 0; i < Nstars; ++i) {
737
738
739
740
741
742
743
744
745
746
747
748

    /* Duplicate the crucial information */
    gparts[i + Ndm].x[0] = sparts[i].x[0];
    gparts[i + Ndm].x[1] = sparts[i].x[1];
    gparts[i + Ndm].x[2] = sparts[i].x[2];

    gparts[i + Ndm].v_full[0] = sparts[i].v[0];
    gparts[i + Ndm].v_full[1] = sparts[i].v[1];
    gparts[i + Ndm].v_full[2] = sparts[i].v[2];

    gparts[i + Ndm].mass = sparts[i].mass;

749
    /* Set gpart type */
750
    gparts[i + Ndm].type = swift_type_stars;
751

752
    /* Link the particles */
753
    gparts[i + Ndm].id_or_neg_offset = -(long long)(offset + i);
754
755
756
757
    sparts[i].gpart = &gparts[i + Ndm];
  }
}

758
759
760
761
762
763
/**
 * @brief Copy every #spart into the corresponding #gpart and link them.
 *
 * This function assumes that the DM particles and gas particles are all at
 * the start of the gparts array and adds the star particles afterwards
 *
764
 * @param tp The current #threadpool.
765
766
767
768
769
770
 * @param sparts The array of #spart freshly read in.
 * @param gparts The array of #gpart freshly read in with all the DM and gas
 * particles at the start.
 * @param Nstars The number of stars particles read in.
 * @param Ndm The number of DM and gas particles read in.
 */
771
void io_duplicate_stars_gparts(struct threadpool* tp, struct spart* const sparts,
772
773
774
775
776
777
778
779
780
781
782
783
                              struct gpart* const gparts, size_t Nstars,
                              size_t Ndm) {

  struct duplication_data data;
  data.gparts = gparts;
  data.sparts = sparts;
  data.Ndm = Ndm;

  threadpool_map(tp, io_duplicate_hydro_sparts_mapper, sparts, Nstars,
                 sizeof(struct spart), 0, &data);
}

784
785
786
787
788
789
790
791
/**
 * @brief Copy every DM #gpart into the dmparts array.
 *
 * @param gparts The array of #gpart containing all particles.
 * @param Ntot The number of #gpart.
 * @param dmparts The array of #gpart containg DM particles to be filled.
 * @param Ndm The number of DM particles.
 */
792
793
void io_collect_dm_gparts(const struct gpart* const gparts, size_t Ntot,
                          struct gpart* const dmparts, size_t Ndm) {
794

795
  size_t count = 0;
796
797

  /* Loop over all gparts */
798
  for (size_t i = 0; i < Ntot; ++i) {
799

Matthieu Schaller's avatar
Matthieu Schaller committed
800
801
    /* message("i=%zd count=%zd id=%lld part=%p", i, count, gparts[i].id,
     * gparts[i].part); */
Matthieu Schaller's avatar
Matthieu Schaller committed
802

803
    /* And collect the DM ones */
804
    if (gparts[i].type == swift_type_dark_matter) {
805
      dmparts[count] = gparts[i];
806
807
808
809
810
811
      count++;
    }
  }

  /* Check that everything is fine */
  if (count != Ndm)
812
    error("Collected the wrong number of dm particles (%zu vs. %zu expected)",
813
814
          count, Ndm);
}
815
816
817
818

/**
 * @brief Verify the io parameter file
 *
819
 * @param params The #swift_params
820
 * @param N_total The total number of each particle type.
821
 */
822
void io_check_output_fields(const struct swift_params* params,
823
                            const long long N_total[3]) {
824

825
826
827
828
829
  /* Create some fake particles as arguments for the writing routines */
  struct part p;
  struct xpart xp;
  struct spart sp;
  struct gpart gp;
830

831
  /* Copy N_total to array with length == 6 */
832
833
  const long long nr_total[swift_type_count] = {N_total[0], N_total[1], 0,
                                                0,          N_total[2], 0};
834

835
  /* Loop over all particle types to check the fields */
836
  for (int ptype = 0; ptype < swift_type_count; ptype++) {
837

838
839
840
841
    int num_fields = 0;
    struct io_props list[100];

    /* Don't do anything if no particle of this kind */
lhausamm's avatar
lhausamm committed
842
    if (nr_total[ptype] == 0) continue;
843

844
    /* Gather particle fields from the particle structures */
845
846
    switch (ptype) {

lhausamm's avatar
lhausamm committed
847
      case swift_type_gas:
848
849
        hydro_write_particles(&p, &xp, list, &num_fields);
        num_fields += chemistry_write_particles(&p, list + num_fields);
lhausamm's avatar
lhausamm committed
850
        break;
851

lhausamm's avatar
lhausamm committed
852
      case swift_type_dark_matter:
853
        darkmatter_write_particles(&gp, list, &num_fields);
lhausamm's avatar
lhausamm committed
854
        break;
855

856
857
      case swift_type_stars:
        stars_write_particles(&sp, list, &num_fields);
lhausamm's avatar
lhausamm committed
858
        break;
859

lhausamm's avatar
lhausamm committed
860
861
862
      default:
        error("Particle Type %d not yet supported. Aborting", ptype);
    }
863

864
    /* loop over each parameter */
865
866
    for (int param_id = 0; param_id < params->paramCount; param_id++) {
      const char* param_name = params->data[param_id].name;
lhausamm's avatar
lhausamm committed
867

lhausamm's avatar
lhausamm committed
868
      char section_name[PARSER_MAX_LINE_SIZE];
869
870

      /* Skip if wrong section */
871
      sprintf(section_name, "SelectOutput:");
lhausamm's avatar
lhausamm committed
872
      if (strstr(param_name, section_name) == NULL) continue;
873

874
      /* Skip if wrong particle type */
875
      sprintf(section_name, "_%s", part_type_names[ptype]);
876
877
      if (strstr(param_name, section_name) == NULL) continue;

878
      int found = 0;
lhausamm's avatar
lhausamm committed
879

880
881
      /* loop over each possible output field */
      for (int field_id = 0; field_id < num_fields; field_id++) {
lhausamm's avatar
lhausamm committed
882
        char field_name[PARSER_MAX_LINE_SIZE];
lhausamm's avatar
Format    
lhausamm committed
883
884
        sprintf(field_name, "SelectOutput:%s_%s", list[field_id].name,
                part_type_names[ptype]);
lhausamm's avatar
lhausamm committed
885

lhausamm's avatar
lhausamm committed
886
887
        if (strcmp(param_name, field_name) == 0) {
          found = 1;
lhausamm's avatar
Format    
lhausamm committed
888
889
890
891
892
          /* check if correct input */
          int retParam = 0;
          char str[PARSER_MAX_LINE_SIZE];
          sscanf(params->data[param_id].value, "%d%s", &retParam, str);

893
          /* Check that we have a 0 or 1 */
lhausamm's avatar
Format    
lhausamm committed
894
          if (retParam != 0 && retParam != 1)
Matthieu Schaller's avatar
Matthieu Schaller committed
895
896
897
898
            message(
                "WARNING: Unexpected input for %s. Received %i but expect 0 or "
                "1. ",
                field_name, retParam);
899
900
901

          /* Found it, so move to the next one. */
          break;
lhausamm's avatar
lhausamm committed
902
        }
903
      }
lhausamm's avatar
lhausamm committed
904
      if (!found)
905
906
907
908
        message(
            "WARNING: Trying to dump particle field '%s' (read from '%s') that "
            "does not exist.",
            param_name, params->fileName);
909
910
911
    }
  }
}
912
913
914
915
916
917
918

/**
 * @brief Write the output field parameters file
 *
 * @param filename The file to write
 */
void io_write_output_field_parameter(const char* filename) {
919

lhausamm's avatar
Format    
lhausamm committed
920
  FILE* file = fopen(filename, "w");
921
  if (file == NULL) error("Error opening file '%s'", filename);
922
923

  /* Loop over all particle types */
924
  fprintf(file, "SelectOutput:\n");
925
  for (int ptype = 0; ptype < swift_type_count; ptype++) {
926

927
928
929
930
931
932
    int num_fields = 0;
    struct io_props list[100];

    /* Write particle fields from the particle structure */
    switch (ptype) {

lhausamm's avatar
Format    
lhausamm committed
933
934
935
936
      case swift_type_gas:
        hydro_write_particles(NULL, NULL, list, &num_fields);
        num_fields += chemistry_write_particles(NULL, list + num_fields);
        break;
937

lhausamm's avatar
Format    
lhausamm committed
938
939
940
      case swift_type_dark_matter:
        darkmatter_write_particles(NULL, list, &num_fields);
        break;
941

942
943
      case swift_type_stars:
        stars_write_particles(NULL, list, &num_fields);
lhausamm's avatar
Format    
lhausamm committed
944
        break;
945
946
947

      default:
        break;
948
949
    }

lhausamm's avatar
Format    
lhausamm committed
950
    if (num_fields == 0) continue;
951

952
    /* Output a header for that particle type */
lhausamm's avatar
lhausamm committed
953
    fprintf(file, "  # Particle Type %s\n", part_type_names[ptype]);
954
955
956

    /* Write all the fields of this particle type */
    for (int i = 0; i < num_fields; ++i)
lhausamm's avatar
Format    
lhausamm committed
957
      fprintf(file, "  %s_%s: 1\n", list[i].name, part_type_names[ptype]);
lhausamm's avatar
Format    
lhausamm committed
958

959
960
961
962
    fprintf(file, "\n");
  }

  fclose(file);
963
964
965
966
967

  printf(
      "List of valid ouput fields for the particle in snapshots dumped in "
      "'%s'.\n",
      filename);
968
}