common_io.c 24.1 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
24
25
26
 ******************************************************************************/

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

#if defined(HAVE_HDF5)

/* Some standard headers. */
27
28
29
#include <hdf5.h>
#include <math.h>
#include <stddef.h>
30
31
32
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
33
34

/* MPI headers. */
35
36
37
#ifdef WITH_MPI
#include <mpi.h>
#endif
38

39
40
41
42
/* This object's header. */
#include "common_io.h"

/* Local includes. */
43
#include "engine.h"
44
#include "error.h"
45
#include "hydro.h"
46
#include "io_properties.h"
47
#include "kernel_hydro.h"
48
#include "part.h"
49
#include "threadpool.h"
50
#include "units.h"
51
#include "version.h"
52
53

/**
54
 * @brief Converts a C data type to the HDF5 equivalent.
55
56
 *
 * This function is a trivial wrapper around the HDF5 types but allows
57
58
 * to change the exact storage types matching the code types in a transparent
 *way.
59
 */
60
hid_t io_hdf5_type(enum IO_DATA_TYPE type) {
61

62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
  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:
80
      return H5T_NATIVE_CHAR;
81
82
83
84
    default:
      error("Unknown type");
      return 0;
  }
85
86
87
88
89
}

/**
 * @brief Returns the memory size of the data type
 */
90
size_t io_sizeof_type(enum IO_DATA_TYPE type) {
91

92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
  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;
  }
115
116
}

117
118
119
120
121
/**
 * @brief Return 1 if the type has double precision
 *
 * Returns an error if the type is not FLOAT or DOUBLE
 */
122
int io_is_double_precision(enum IO_DATA_TYPE type) {
Matthieu Schaller's avatar
Matthieu Schaller committed
123

124
125
126
127
128
129
130
131
132
133
134
  switch (type) {
    case FLOAT:
      return 0;
    case DOUBLE:
      return 1;
    default:
      error("Invalid type");
      return 0;
  }
}

135
136
137
138
139
/**
 * @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.
140
 * @param type The #IO_DATA_TYPE of the attribute.
141
142
143
144
 * @param data (output) The attribute read from the HDF5 group.
 *
 * Calls #error() if an error occurs.
 */
145
146
void io_read_attribute(hid_t grp, char* name, enum IO_DATA_TYPE type,
                       void* data) {
147
  hid_t h_attr = 0, h_err = 0;
148
149

  h_attr = H5Aopen(grp, name, H5P_DEFAULT);
150
151
152
  if (h_attr < 0) {
    error("Error while opening attribute '%s'", name);
  }
153

154
  h_err = H5Aread(h_attr, io_hdf5_type(type), data);
155
156
157
  if (h_err < 0) {
    error("Error while reading attribute '%s'", name);
  }
158
159
160
161
162
163
164
165
166

  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.
167
 * @param type The #IO_DATA_TYPE of the attribute.
168
169
170
171
172
 * @param data The attribute to write.
 * @param num The number of elements to write
 *
 * Calls #error() if an error occurs.
 */
173
174
void io_write_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
                        void* data, int num) {
175
176
  hid_t h_space = 0, h_attr = 0, h_err = 0;
  hsize_t dim[1] = {num};
177
178

  h_space = H5Screate(H5S_SIMPLE);
179
180
181
  if (h_space < 0) {
    error("Error while creating dataspace for attribute '%s'.", name);
  }
182
183

  h_err = H5Sset_extent_simple(h_space, 1, dim, NULL);
184
185
186
  if (h_err < 0) {
    error("Error while changing dataspace shape for attribute '%s'.", name);
  }
187

188
  h_attr = H5Acreate1(grp, name, io_hdf5_type(type), h_space, H5P_DEFAULT);
189
190
191
  if (h_attr < 0) {
    error("Error while creating attribute '%s'.", name);
  }
192

193
  h_err = H5Awrite(h_attr, io_hdf5_type(type), data);
194
195
196
  if (h_err < 0) {
    error("Error while reading attribute '%s'.", name);
  }
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211

  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.
 */
212
213
void io_writeStringAttribute(hid_t grp, const char* name, const char* str,
                             int length) {
214
  hid_t h_space = 0, h_attr = 0, h_err = 0, h_type = 0;
215
216

  h_space = H5Screate(H5S_SCALAR);
217
218
219
  if (h_space < 0) {
    error("Error while creating dataspace for attribute '%s'.", name);
  }
220
221

  h_type = H5Tcopy(H5T_C_S1);
222
223
224
  if (h_type < 0) {
    error("Error while copying datatype 'H5T_C_S1'.");
  }
225
226

  h_err = H5Tset_size(h_type, length);
227
  if (h_err < 0) {
Matthieu Schaller's avatar
Matthieu Schaller committed
228
    error("Error while resizing attribute type to '%i'.", length);
229
  }
230
231

  h_attr = H5Acreate1(grp, name, h_type, h_space, H5P_DEFAULT);
232
233
234
  if (h_attr < 0) {
    error("Error while creating attribute '%s'.", name);
  }
235

236
237
238
239
  h_err = H5Awrite(h_attr, h_type, str);
  if (h_err < 0) {
    error("Error while reading attribute '%s'.", name);
  }
240
241
242
243
244
245
246
247

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

/**
 * @brief Writes a double value as an attribute
248
 * @param grp The group in which to write
249
250
251
 * @param name The name of the attribute
 * @param data The value to write
 */
252
253
void io_write_attribute_d(hid_t grp, const char* name, double data) {
  io_write_attribute(grp, name, DOUBLE, &data, 1);
254
255
256
257
}

/**
 * @brief Writes a float value as an attribute
258
 * @param grp The group in which to write
259
260
261
 * @param name The name of the attribute
 * @param data The value to write
 */
262
263
void io_write_attribute_f(hid_t grp, const char* name, float data) {
  io_write_attribute(grp, name, FLOAT, &data, 1);
264
265
266
267
}

/**
 * @brief Writes an int value as an attribute
268
 * @param grp The group in which to write
269
270
271
 * @param name The name of the attribute
 * @param data The value to write
 */
272
273
void io_write_attribute_i(hid_t grp, const char* name, int data) {
  io_write_attribute(grp, name, INT, &data, 1);
274
275
276
277
}

/**
 * @brief Writes a long value as an attribute
278
 * @param grp The group in which to write
279
280
281
 * @param name The name of the attribute
 * @param data The value to write
 */
282
283
void io_write_attribute_l(hid_t grp, const char* name, long data) {
  io_write_attribute(grp, name, LONG, &data, 1);
284
285
286
287
}

/**
 * @brief Writes a string value as an attribute
288
 * @param grp The group in which to write
289
290
291
 * @param name The name of the attribute
 * @param str The string to write
 */
292
293
void io_write_attribute_s(hid_t grp, const char* name, const char* str) {
  io_writeStringAttribute(grp, name, str, strlen(str));
294
295
}

296
297
298
/**
 * @brief Reads the Unit System from an IC file.
 * @param h_file The (opened) HDF5 file from which to read.
299
 * @param us The unit_system to fill.
300
 * @param mpi_rank The MPI rank we are on.
301
302
303
 *
 * If the 'Units' group does not exist in the ICs, cgs units will be assumed
 */
304
void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank) {
305
306
307

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

309
  if (exists == 0) {
310
311
312

    if (mpi_rank == 0)
      message("'Units' group not found in ICs. Assuming CGS unit system.");
313
314
315
316
317
318
319
320
321

    /* Default to CGS */
    us->UnitMass_in_cgs = 1.;
    us->UnitLength_in_cgs = 1.;
    us->UnitTime_in_cgs = 1.;
    us->UnitCurrent_in_cgs = 1.;
    us->UnitTemperature_in_cgs = 1.;

    return;
322
  } else if (exists < 0) {
323
    error("Serious problem with 'Units' group in ICs. H5Lexists gives %d",
324
          exists);
325
  }
326

327
  if (mpi_rank == 0) message("Reading IC units from ICs.");
328
  hid_t h_grp = H5Gopen(h_file, "/Units", H5P_DEFAULT);
329
330

  /* Ok, Read the damn thing */
331
332
333
334
335
336
337
338
339
340
  io_read_attribute(h_grp, "Unit length in cgs (U_L)", DOUBLE,
                    &us->UnitLength_in_cgs);
  io_read_attribute(h_grp, "Unit mass in cgs (U_M)", DOUBLE,
                    &us->UnitMass_in_cgs);
  io_read_attribute(h_grp, "Unit time in cgs (U_t)", DOUBLE,
                    &us->UnitTime_in_cgs);
  io_read_attribute(h_grp, "Unit current in cgs (U_I)", DOUBLE,
                    &us->UnitCurrent_in_cgs);
  io_read_attribute(h_grp, "Unit temperature in cgs (U_T)", DOUBLE,
                    &us->UnitTemperature_in_cgs);
341
342
343
344
345

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

346
347
348
/**
 * @brief Writes the current Unit System
 * @param h_file The (opened) HDF5 file in which to write
349
 * @param us The unit_system to dump
350
 * @param groupName The name of the HDF5 group to write to
351
 */
352
353
void io_write_unit_system(hid_t h_file, const struct unit_system* us,
                          const char* groupName) {
354

355
356
  hid_t h_grpunit = 0;
  h_grpunit = H5Gcreate1(h_file, groupName, 0);
357
358
  if (h_grpunit < 0) error("Error while creating Unit System group");

359
360
361
362
363
364
365
366
367
368
  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));
369
370
371
372

  H5Gclose(h_grpunit);
}

373
374
375
376
/**
 * @brief Writes the code version to the file
 * @param h_file The (opened) HDF5 file in which to write
 */
377
void io_write_code_description(hid_t h_file) {
378
379
380
381
382
  hid_t h_grpcode = 0;

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

383
  io_write_attribute_s(h_grpcode, "Code", "SWIFT");
384
385
386
387
388
389
390
391
392
393
  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());
394
  io_write_attribute_s(h_grpcode, "Thread barriers", thread_barrier_version());
Matthieu Schaller's avatar
Matthieu Schaller committed
395
#ifdef HAVE_FFTW
396
  io_write_attribute_s(h_grpcode, "FFTW library version", fftw3_version());
Matthieu Schaller's avatar
Matthieu Schaller committed
397
#endif
398
#ifdef WITH_MPI
399
  io_write_attribute_s(h_grpcode, "MPI library", mpi_version());
400
#ifdef HAVE_METIS
401
  io_write_attribute_s(h_grpcode, "METIS library version", metis_version());
402
#endif
403
#else
404
  io_write_attribute_s(h_grpcode, "MPI library", "Non-MPI version of SWIFT");
405
#endif
406
407
408
  H5Gclose(h_grpcode);
}

409
410
#endif /* HAVE_HDF5 */

411
412
413
/**
 * @brief Mapper function to copy #part or #gpart fields into a buffer.
 */
414
void io_copy_mapper(void* restrict temp, int N, void* restrict extra_data) {
415

416
  const struct io_props props = *((const struct io_props*)extra_data);
417
418
419
420
  const size_t typeSize = io_sizeof_type(props.type);
  const size_t copySize = typeSize * props.dimension;

  /* How far are we with this chunk? */
421
  char* restrict temp_c = temp;
422
423
  const ptrdiff_t delta = (temp_c - props.start_temp_c) / copySize;

424
425
426
  for (int k = 0; k < N; k++) {
    memcpy(&temp_c[k * copySize], props.field + (delta + k) * props.partSize,
           copySize);
427
428
429
430
  }
}

/**
431
432
 * @brief Mapper function to copy #part into a buffer of floats using a
 * conversion function.
433
 */
434
435
void io_convert_part_f_mapper(void* restrict temp, int N,
                              void* restrict extra_data) {
436

437
438
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct part* restrict parts = props.parts;
439
440
441
442
  const struct engine* e = props.e;
  const size_t dim = props.dimension;

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

446
447
448
449
450
  for (int i = 0; i < N; i++)
    props.convert_part_f(e, parts + delta + i, &temp_f[i * dim]);
}

/**
451
452
 * @brief Mapper function to copy #part into a buffer of doubles using a
 * conversion function.
453
 */
454
455
void io_convert_part_d_mapper(void* restrict temp, int N,
                              void* restrict extra_data) {
456

457
458
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct part* restrict parts = props.parts;
459
460
461
462
  const struct engine* e = props.e;
  const size_t dim = props.dimension;

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

466
467
468
469
470
  for (int i = 0; i < N; i++)
    props.convert_part_d(e, parts + delta + i, &temp_d[i * dim]);
}

/**
471
472
 * @brief Mapper function to copy #gpart into a buffer of floats using a
 * conversion function.
473
 */
474
475
void io_convert_gpart_f_mapper(void* restrict temp, int N,
                               void* restrict extra_data) {
476

477
478
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct gpart* restrict gparts = props.gparts;
479
480
481
482
  const struct engine* e = props.e;
  const size_t dim = props.dimension;

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

486
487
488
489
490
  for (int i = 0; i < N; i++)
    props.convert_gpart_f(e, gparts + delta + i, &temp_f[i * dim]);
}

/**
491
492
 * @brief Mapper function to copy #gpart into a buffer of doubles using a
 * conversion function.
493
 */
494
495
void io_convert_gpart_d_mapper(void* restrict temp, int N,
                               void* restrict extra_data) {
496

497
498
  const struct io_props props = *((const struct io_props*)extra_data);
  const struct gpart* restrict gparts = props.gparts;
499
500
501
502
  const struct engine* e = props.e;
  const size_t dim = props.dimension;

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

506
507
508
509
  for (int i = 0; i < N; i++)
    props.convert_gpart_d(e, gparts + delta + i, &temp_d[i * dim]);
}

510
511
512
513
514
/**
 * @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.
515
516
 * @param props The #io_props corresponding to the particle field we are
 * copying.
517
518
519
520
521
 * @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,
522
                         struct io_props props, size_t N,
523
524
525
526
527
528
529
530
531
532
                         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 */

533
    /* Prepare some parameters */
534
    char* temp_c = temp;
535
    props.start_temp_c = temp_c;
536

537
    /* Copy the whole thing into a buffer */
538
539
540
    threadpool_map((struct threadpool*)&e->threadpool, io_copy_mapper, temp_c,
                   N, copySize, 0, (void*)&props);

541
542
543
544
  } else { /* Converting particle to data */

    if (props.convert_part_f != NULL) {

545
      /* Prepare some parameters */
546
      float* temp_f = temp;
547
548
      props.start_temp_f = temp;
      props.e = e;
549

550
      /* Copy the whole thing into a buffer */
551
552
553
      threadpool_map((struct threadpool*)&e->threadpool,
                     io_convert_part_f_mapper, temp_f, N, copySize, 0,
                     (void*)&props);
554
555

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

557
      /* Prepare some parameters */
558
      double* temp_d = temp;
559
560
      props.start_temp_d = temp;
      props.e = e;
561

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

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

569
      /* Prepare some parameters */
570
      float* temp_f = temp;
571
572
      props.start_temp_f = temp;
      props.e = e;
573

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

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

581
      /* Prepare some parameters */
582
      double* temp_d = temp;
583
584
      props.start_temp_d = temp;
      props.e = e;
585

586
      /* Copy the whole thing into a buffer */
587
588
589
      threadpool_map((struct threadpool*)&e->threadpool,
                     io_convert_gpart_d_mapper, temp_d, N, copySize, 0,
                     (void*)&props);
590
591
592
593
594
595
596
597
598
599
600
601
602
603

    } 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)) {
604
      swift_declare_aligned_ptr(double, temp_d, temp, IO_BUFFER_ALIGNMENT);
605
606
      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= factor;
    } else {
607
      swift_declare_aligned_ptr(float, temp_f, temp, IO_BUFFER_ALIGNMENT);
608
609
610
611
612
      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor;
    }
  }
}

613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
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;
  }
}

630
631
632
633
634
635
636
/**
 * @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
 *
637
 * @param tp The current #threadpool.
638
639
640
 * @param gparts The array of #gpart freshly read in.
 * @param Ndm The number of DM particles read in.
 */
641
642
void io_prepare_dm_gparts(struct threadpool* tp, struct gpart* const gparts,
                          size_t Ndm) {
643

644
645
  threadpool_map(tp, io_prepare_dm_gparts_mapper, gparts, Ndm,
                 sizeof(struct gpart), 0, NULL);
646
647
}

648
649
650
651
652
653
654
655
656
struct duplication_data {

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

658
659
660
661
662
663
664
665
666
667
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) {
668
669
670
671
672
673
674
675
676
677

    /* 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];

678
    gparts[i + Ndm].mass = hydro_get_mass(&parts[i]);
679

680
681
682
    /* Set gpart type */
    gparts[i + Ndm].type = swift_type_gas;

683
    /* Link the particles */
684
    gparts[i + Ndm].id_or_neg_offset = -i;
685
686
687
688
    parts[i].gpart = &gparts[i + Ndm];
  }
}

689
/**
690
 * @brief Copy every #part into the corresponding #gpart and link them.
691
 *
692
693
 * This function assumes that the DM particles are all at the start of the
 * gparts array and adds the hydro particles afterwards
694
 *
695
 * @param tp The current #threadpool.
696
697
698
699
700
 * @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.
701
 */
702
703
704
705
706
707
708
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;
709

710
711
712
713
714
715
716
717
718
719
720
721
722
723
  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) {
724
725
726
727
728
729
730
731
732
733
734
735

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

736
737
738
    /* Set gpart type */
    gparts[i + Ndm].type = swift_type_star;

739
740
741
742
743
744
    /* Link the particles */
    gparts[i + Ndm].id_or_neg_offset = -i;
    sparts[i].gpart = &gparts[i + Ndm];
  }
}

745
746
747
748
749
750
/**
 * @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
 *
751
 * @param tp The current #threadpool.
752
753
754
755
756
757
758
759
760
761
762
763
764
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.
 */
void io_duplicate_star_gparts(struct threadpool* tp, struct spart* const sparts,
                              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);
}

771
772
773
774
775
776
777
778
/**
 * @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.
 */
779
780
void io_collect_dm_gparts(const struct gpart* const gparts, size_t Ntot,
                          struct gpart* const dmparts, size_t Ndm) {
781

782
  size_t count = 0;
783
784

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

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

790
    /* And collect the DM ones */
791
    if (gparts[i].type == swift_type_dark_matter) {
792
      dmparts[count] = gparts[i];
793
794
795
796
797
798
      count++;
    }
  }

  /* Check that everything is fine */
  if (count != Ndm)
799
    error("Collected the wrong number of dm particles (%zu vs. %zu expected)",
800
801
          count, Ndm);
}