common_io.c 22.7 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
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
/**
 * @brief Mapper function to copy #part or #gpart fields into a buffer.
 */
void io_copy_mapper(void *restrict temp, int N,
		    void *restrict extra_data) {

  const struct io_props props = *((const struct io_props*) extra_data);
  const size_t typeSize = io_sizeof_type(props.type);
  const size_t copySize = typeSize * props.dimension;

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

  for(int k = 0; k < N; k++) {
    memcpy(&temp_c[k * copySize], props.field + (delta + k) * props.partSize, copySize);
  }
}

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

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

  /* How far are we with this chunk? */
  float *restrict temp_f = temp;
  const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
  
  for (int i = 0; i < N; i++)
    props.convert_part_f(e, parts + delta + i, &temp_f[i * dim]);
}

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

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

  /* How far are we with this chunk? */
  double *restrict temp_d = temp;
  const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
  
  for (int i = 0; i < N; i++)
    props.convert_part_d(e, parts + delta + i, &temp_d[i * dim]);
}

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

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

  /* How far are we with this chunk? */
  float *restrict temp_f = temp;
  const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
  
  for (int i = 0; i < N; i++)
    props.convert_gpart_f(e, gparts + delta + i, &temp_f[i * dim]);
}

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

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

  /* How far are we with this chunk? */
  double *restrict temp_d = temp;
  const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
  
  for (int i = 0; i < N; i++)
    props.convert_gpart_d(e, gparts + delta + i, &temp_d[i * dim]);
}


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

528
    /* Prepare some parameters */
529
    char* temp_c = temp;
530
    props.start_temp_c = temp_c;
531

532
533
534
535
    /* Copy the whole thing into a buffer */
    threadpool_map((struct threadpool*) &e->threadpool, io_copy_mapper,
    		   temp_c, N, copySize, 0, (void*) &props);
		   
536
537
538
539
  } else { /* Converting particle to data */

    if (props.convert_part_f != NULL) {

540
541
542
543
      /* Prepare some parameters */
      float *temp_f = temp;
      props.start_temp_f = temp;
      props.e = e;
544

545
546
547
      /* Copy the whole thing into a buffer */
      threadpool_map((struct threadpool*) &e->threadpool, io_convert_part_f_mapper,
		     temp_f, N, copySize, 0, (void*) &props);
548
549

    } else if (props.convert_part_d != NULL) {
550
551
552
553
554
      
      /* Prepare some parameters */
      double *temp_d = temp;
      props.start_temp_d = temp;
      props.e = e;
555

556
557
558
      /* Copy the whole thing into a buffer */
      threadpool_map((struct threadpool*) &e->threadpool, io_convert_part_d_mapper,
		     temp_d, N, copySize, 0, (void*) &props);
559
560
561

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

562
563
564
565
      /* Prepare some parameters */
      float *temp_f = temp;
      props.start_temp_f = temp;
      props.e = e;
566

567
568
569
      /* Copy the whole thing into a buffer */
      threadpool_map((struct threadpool*) &e->threadpool, io_convert_gpart_f_mapper,
		     temp_f, N, copySize, 0, (void*) &props);
570
571
572

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

573
574
575
576
      /* Prepare some parameters */
      double *temp_d = temp;
      props.start_temp_d = temp;
      props.e = e;
577

578
579
580
      /* Copy the whole thing into a buffer */
      threadpool_map((struct threadpool*) &e->threadpool, io_convert_gpart_d_mapper,
		     temp_d, N, copySize, 0, (void*) &props);
581
582
583
584
585
586
587
588
589
590
591
592
593
594

    } 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)) {
595
      swift_declare_aligned_ptr(double, temp_d, temp, IO_BUFFER_ALIGNMENT);
596
597
      for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= factor;
    } else {
598
      swift_declare_aligned_ptr(float, temp_f, temp, IO_BUFFER_ALIGNMENT);
599
600
601
602
603
      for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor;
    }
  }
}

604
605
#endif /* HAVE_HDF5 */

606
607
608
609
610
611
/* ------------------------------------------------------------------------------------------------
 * This part writes the XMF file descriptor enabling a visualisation through
 * ParaView
 * ------------------------------------------------------------------------------------------------
 */

612
613
614
615
616
617
618
619
/**
 * @brief Prepares the XMF file for the new entry
 *
 * Creates a temporary file on the disk in order to copy the right lines.
 *
 * @todo Use a proper XML library to avoid stupid copies.
 */

620
621
622
623
624
625
626
627
628
629
/**
 * @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
 *
 * @param gparts The array of #gpart freshly read in.
 * @param Ndm The number of DM particles read in.
 */
630
void io_prepare_dm_gparts(struct gpart* const gparts, size_t Ndm) {
631
632

  /* Let's give all these gparts a negative id */
633
  for (size_t i = 0; i < Ndm; ++i) {
Matthieu Schaller's avatar
typo    
Matthieu Schaller committed
634
    /* 0 or negative ids are not allowed */
635
    if (gparts[i].id_or_neg_offset <= 0)
636
      error("0 or negative ID for DM particle %zu: ID=%lld", i,
637
            gparts[i].id_or_neg_offset);
638
639

    /* Set gpart type */
640
    gparts[i].type = swift_type_dark_matter;
641
642
643
644
645
646
647
648
649
650
651
  }
}

/**
 * @brief Copy every #part into the corresponding #gpart and link them.
 *
 * This function assumes that the DM particles are all at the start of the
 * gparts array and adds the hydro particles afterwards
 *
 * @param parts The array of #part freshly read in.
 * @param gparts The array of #gpart freshly read in with all the DM particles
652
 * at the start
653
654
655
 * @param Ngas The number of gas particles read in.
 * @param Ndm The number of DM particles read in.
 */
656
657
658
void io_duplicate_hydro_gparts(struct part* const parts,
                               struct gpart* const gparts, size_t Ngas,
                               size_t Ndm) {
659

660
  for (size_t i = 0; i < Ngas; ++i) {
661
662
663
664
665
666
667
668
669
670

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

671
    gparts[i + Ndm].mass = hydro_get_mass(&parts[i]);
672

673
674
675
    /* Set gpart type */
    gparts[i + Ndm].type = swift_type_gas;

676
    /* Link the particles */
677
    gparts[i + Ndm].id_or_neg_offset = -i;
678
679
680
681
    parts[i].gpart = &gparts[i + Ndm];
  }
}

682
683
684
685
686
687
688
689
690
691
692
693
/**
 * @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
 *
 * @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.
 */
694
695
696
void io_duplicate_star_gparts(struct spart* const sparts,
                              struct gpart* const gparts, size_t Nstars,
                              size_t Ndm) {
697
698
699
700
701
702
703
704
705
706
707
708
709
710

  for (size_t i = 0; i < Nstars; ++i) {

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

711
712
713
    /* Set gpart type */
    gparts[i + Ndm].type = swift_type_star;

714
715
716
717
718
719
    /* Link the particles */
    gparts[i + Ndm].id_or_neg_offset = -i;
    sparts[i].gpart = &gparts[i + Ndm];
  }
}

720
721
722
723
724
725
726
727
/**
 * @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.
 */
728
729
void io_collect_dm_gparts(const struct gpart* const gparts, size_t Ntot,
                          struct gpart* const dmparts, size_t Ndm) {
730

731
  size_t count = 0;
732
733

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

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

739
    /* And collect the DM ones */
740
    if (gparts[i].type == swift_type_dark_matter) {
741
      dmparts[count] = gparts[i];
742
743
744
745
746
747
      count++;
    }
  }

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