common_io.c 20.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
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 "error.h"
44
#include "hydro.h"
45
#include "kernel_hydro.h"
46
47
#include "part.h"
#include "units.h"
48
#include "version.h"
49

Matthieu Schaller's avatar
Matthieu Schaller committed
50
51
const char* particle_type_names[NUM_PARTICLE_TYPES] = {
    "Gas", "DM", "Boundary", "Dummy", "Star", "BH"};
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 hdf5Type(enum DATA_TYPE type) {
61

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

/**
 * @brief Returns the memory size of the data type
 */
90
size_t sizeOfType(enum 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
122
/**
 * @brief Return 1 if the type has double precision
 *
 * Returns an error if the type is not FLOAT or DOUBLE
 */
int isDoublePrecision(enum 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
140
141
142
143
144
/**
 * @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.
 * @param type The #DATA_TYPE of the attribute.
 * @param data (output) The attribute read from the HDF5 group.
 *
 * Calls #error() if an error occurs.
 */
145
146
void readAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data) {
  hid_t h_attr = 0, h_err = 0;
147
148

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

/**
 * @brief Writes an int value as an attribute
267
 * @param grp The group in which to write
268
269
270
271
 * @param name The name of the attribute
 * @param data The value to write
 */

272
void writeAttribute_i(hid_t grp, const char* name, int data) {
273
274
275
276
277
  writeAttribute(grp, name, INT, &data, 1);
}

/**
 * @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
void writeAttribute_l(hid_t grp, const char* name, long data) {
283
284
285
286
287
  writeAttribute(grp, name, LONG, &data, 1);
}

/**
 * @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
void writeAttribute_s(hid_t grp, const char* name, const char* str) {
293
294
295
  writeStringAttribute(grp, name, str, strlen(str));
}

296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
/**
 * @brief Reads the Unit System from an IC file.
 * @param h_file The (opened) HDF5 file from which to read.
 * @param us The UnitSystem to fill.
 *
 * If the 'Units' group does not exist in the ICs, cgs units will be assumed
 */
void readUnitSystem(hid_t h_file, struct UnitSystem* us) {

  hid_t h_grp = H5Gopen(h_file, "/Units", H5P_DEFAULT);

  if (h_grp < 0) {
    message("'Units' group not found in ICs. Assuming CGS unit system.");

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

  /* Ok, Read the damn thing */
  readAttribute(h_grp, "Unit length in cgs (U_L)", DOUBLE,
                &us->UnitLength_in_cgs);
  readAttribute(h_grp, "Unit mass in cgs (U_M)", DOUBLE, &us->UnitMass_in_cgs);
  readAttribute(h_grp, "Unit time in cgs (U_t)", DOUBLE, &us->UnitTime_in_cgs);
  readAttribute(h_grp, "Unit current in cgs (U_I)", DOUBLE,
                &us->UnitCurrent_in_cgs);
  readAttribute(h_grp, "Unit temperature in cgs (U_T)", DOUBLE,
                &us->UnitTemperature_in_cgs);

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

334
335
336
/**
 * @brief Writes the current Unit System
 * @param h_file The (opened) HDF5 file in which to write
337
338
 * @param us The UnitSystem to dump
 * @param groupName The name of the HDF5 group to write to
339
 */
340
341
void writeUnitSystem(hid_t h_file, const struct UnitSystem* us,
                     const char* groupName) {
342

343
344
  hid_t h_grpunit = 0;
  h_grpunit = H5Gcreate1(h_file, groupName, 0);
345
346
347
  if (h_grpunit < 0) error("Error while creating Unit System group");

  writeAttribute_d(h_grpunit, "Unit mass in cgs (U_M)",
348
                   units_get_base_unit(us, UNIT_MASS));
349
  writeAttribute_d(h_grpunit, "Unit length in cgs (U_L)",
350
                   units_get_base_unit(us, UNIT_LENGTH));
351
  writeAttribute_d(h_grpunit, "Unit time in cgs (U_t)",
352
                   units_get_base_unit(us, UNIT_TIME));
353
  writeAttribute_d(h_grpunit, "Unit current in cgs (U_I)",
354
                   units_get_base_unit(us, UNIT_CURRENT));
355
  writeAttribute_d(h_grpunit, "Unit temperature in cgs (U_T)",
356
                   units_get_base_unit(us, UNIT_TEMPERATURE));
357
358
359
360

  H5Gclose(h_grpunit);
}

361
362
363
364
365
366
367
368
369
370
371
/**
 * @brief Writes the code version to the file
 * @param h_file The (opened) HDF5 file in which to write
 */
void writeCodeDescription(hid_t h_file) {
  hid_t h_grpcode = 0;

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

  writeAttribute_s(h_grpcode, "Code Version", package_version());
372
373
  writeAttribute_s(h_grpcode, "Compiler Name", compiler_name());
  writeAttribute_s(h_grpcode, "Compiler Version", compiler_version());
374
375
  writeAttribute_s(h_grpcode, "Git Branch", git_branch());
  writeAttribute_s(h_grpcode, "Git Revision", git_revision());
376
377
  writeAttribute_s(h_grpcode, "Configuration options", configuration_options());
  writeAttribute_s(h_grpcode, "CFLAGS", compilation_cflags());
378
  writeAttribute_s(h_grpcode, "HDF5 library version", hdf5_version());
Matthieu Schaller's avatar
Matthieu Schaller committed
379
380
381
#ifdef HAVE_FFTW
  writeAttribute_s(h_grpcode, "FFTW library version", fftw3_version());
#endif
382
383
384
#ifdef WITH_MPI
  writeAttribute_s(h_grpcode, "MPI library", mpi_version());
#ifdef HAVE_METIS
385
386
  writeAttribute_s(h_grpcode, "METIS library version", metis_version());
#endif
387
388
389
#else
  writeAttribute_s(h_grpcode, "MPI library", "Non-MPI version of SWIFT");
#endif
390
391
392
  H5Gclose(h_grpcode);
}

393
394
395
396
397
398
/* ------------------------------------------------------------------------------------------------
 * This part writes the XMF file descriptor enabling a visualisation through
 * ParaView
 * ------------------------------------------------------------------------------------------------
 */

399
400
401
402
403
404
405
/**
 * @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.
 */
406
FILE* prepareXMFfile(const char* baseName) {
407
408
  char buffer[1024];

409
410
411
412
413
414
  char fileName[FILENAME_BUFFER_SIZE];
  char tempFileName[FILENAME_BUFFER_SIZE];
  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s.xmf", baseName);
  snprintf(tempFileName, FILENAME_BUFFER_SIZE, "%s_temp.xmf", baseName);
  FILE* xmfFile = fopen(fileName, "r");
  FILE* tempFile = fopen(tempFileName, "w");
415

416
  if (xmfFile == NULL) error("Unable to open current XMF file.");
417

418
  if (tempFile == NULL) error("Unable to open temporary file.");
419
420
421

  /* First we make a temporary copy of the XMF file and count the lines */
  int counter = 0;
422
423
424
425
  while (fgets(buffer, 1024, xmfFile) != NULL) {
    counter++;
    fprintf(tempFile, "%s", buffer);
  }
426
427
  fclose(tempFile);
  fclose(xmfFile);
428

429
  /* We then copy the XMF file back up to the closing lines */
430
431
  xmfFile = fopen(fileName, "w");
  tempFile = fopen(tempFileName, "r");
432

433
  if (xmfFile == NULL) error("Unable to open current XMF file.");
434

435
  if (tempFile == NULL) error("Unable to open temporary file.");
436

437
438
439
  int i = 0;
  while (fgets(buffer, 1024, tempFile) != NULL && i < counter - 3) {
    i++;
440
441
    fprintf(xmfFile, "%s", buffer);
  }
442
443
  fprintf(xmfFile, "\n");
  fclose(tempFile);
444
  remove(tempFileName);
445

446
447
448
449
450
451
  return xmfFile;
}

/**
 * @brief Writes the begin of the XMF file
 *
452
453
 * @todo Exploit the XML nature of the XMF format to write a proper XML writer
 *and simplify all the XMF-related stuff.
454
 */
455
456
457
458
459
void createXMFfile(const char* baseName) {

  char fileName[FILENAME_BUFFER_SIZE];
  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s.xmf", baseName);
  FILE* xmfFile = fopen(fileName, "w");
460
461
462

  fprintf(xmfFile, "<?xml version=\"1.0\" ?> \n");
  fprintf(xmfFile, "<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []> \n");
463
464
465
  fprintf(
      xmfFile,
      "<Xdmf xmlns:xi=\"http://www.w3.org/2003/XInclude\" Version=\"2.1\">\n");
466
  fprintf(xmfFile, "<Domain>\n");
467
468
469
  fprintf(xmfFile,
          "<Grid Name=\"TimeSeries\" GridType=\"Collection\" "
          "CollectionType=\"Temporal\">\n\n");
470
471
472
473
474
475
476
477
478

  fprintf(xmfFile, "</Grid>\n");
  fprintf(xmfFile, "</Domain>\n");
  fprintf(xmfFile, "</Xdmf>\n");

  fclose(xmfFile);
}

/**
479
480
 * @brief Writes the part of the XMF entry presenting the geometry of the
 *snapshot
481
482
483
484
485
 *
 * @param xmfFile The file to write in.
 * @param hdfFileName The name of the HDF5 file corresponding to this output.
 * @param time The current simulation time.
 */
486
void writeXMFoutputheader(FILE* xmfFile, char* hdfFileName, float time) {
487
  /* Write end of file */
488

489
  fprintf(xmfFile, "<!-- XMF description for file: %s -->\n", hdfFileName);
490
491
  fprintf(xmfFile,
          "<Grid GridType=\"Collection\" CollectionType=\"Spatial\">\n");
492
493
494
495
496
497
498
  fprintf(xmfFile, "<Time Type=\"Single\" Value=\"%f\"/>\n", time);
}

/**
 * @brief Writes the end of the XMF file (closes all open markups)
 *
 * @param xmfFile The file to write in.
499
500
 * @param output The number of this output.
 * @param time The current simulation time.
501
 */
502
void writeXMFoutputfooter(FILE* xmfFile, int output, float time) {
503
  /* Write end of the section of this time step */
504

505
506
507
  fprintf(xmfFile,
          "\n</Grid> <!-- End of meta-data for output=%03i, time=%f -->\n",
          output, time);
508
  fprintf(xmfFile, "\n</Grid> <!-- timeSeries -->\n");
509
510
  fprintf(xmfFile, "</Domain>\n");
  fprintf(xmfFile, "</Xdmf>\n");
511

512
513
514
  fclose(xmfFile);
}

515
void writeXMFgroupheader(FILE* xmfFile, char* hdfFileName, size_t N,
Matthieu Schaller's avatar
Matthieu Schaller committed
516
                         enum PARTICLE_TYPE ptype) {
517
  fprintf(xmfFile, "\n<Grid Name=\"%s\" GridType=\"Uniform\">\n",
Matthieu Schaller's avatar
Matthieu Schaller committed
518
          particle_type_names[ptype]);
519
  fprintf(xmfFile,
520
          "<Topology TopologyType=\"Polyvertex\" Dimensions=\"%zu\"/>\n", N);
521
522
  fprintf(xmfFile, "<Geometry GeometryType=\"XYZ\">\n");
  fprintf(xmfFile,
523
          "<DataItem Dimensions=\"%zu 3\" NumberType=\"Double\" "
524
525
          "Precision=\"8\" "
          "Format=\"HDF\">%s:/PartType%d/Coordinates</DataItem>\n",
526
          N, hdfFileName, (int)ptype);
527
528
529
530
  fprintf(xmfFile,
          "</Geometry>\n <!-- Done geometry for %s, start of particle fields "
          "list -->\n",
          particle_type_names[ptype]);
531
532
533
}

void writeXMFgroupfooter(FILE* xmfFile, enum PARTICLE_TYPE ptype) {
534
  fprintf(xmfFile, "</Grid> <!-- End of meta-data for parttype=%s -->\n",
Matthieu Schaller's avatar
Matthieu Schaller committed
535
          particle_type_names[ptype]);
536
537
}

538
539
540
541
542
/**
 * @brief Writes the lines corresponding to an array of the HDF5 output
 *
 * @param xmfFile The file in which to write
 * @param fileName The name of the HDF5 file associated to this XMF descriptor.
543
544
 * @param partTypeGroupName The name of the group containing the particles in
 *the HDF5 file.
545
546
547
548
549
550
551
 * @param name The name of the array in the HDF5 file.
 * @param N The number of particles.
 * @param dim The dimension of the quantity (1 for scalars, 3 for vectors).
 * @param type The type of the data to write.
 *
 * @todo Treat the types in a better way.
 */
552
553
554
void writeXMFline(FILE* xmfFile, const char* fileName,
                  const char* partTypeGroupName, const char* name, size_t N,
                  int dim, enum DATA_TYPE type) {
555
556
557
558
559
  fprintf(xmfFile,
          "<Attribute Name=\"%s\" AttributeType=\"%s\" Center=\"Node\">\n",
          name, dim == 1 ? "Scalar" : "Vector");
  if (dim == 1)
    fprintf(xmfFile,
560
            "<DataItem Dimensions=\"%zu\" NumberType=\"Double\" "
561
562
            "Precision=\"%d\" Format=\"HDF\">%s:%s/%s</DataItem>\n",
            N, type == FLOAT ? 4 : 8, fileName, partTypeGroupName, name);
563
  else
564
    fprintf(xmfFile,
565
            "<DataItem Dimensions=\"%zu %d\" NumberType=\"Double\" "
566
567
            "Precision=\"%d\" Format=\"HDF\">%s:%s/%s</DataItem>\n",
            N, dim, type == FLOAT ? 4 : 8, fileName, partTypeGroupName, name);
568
569
570
  fprintf(xmfFile, "</Attribute>\n");
}

571
572
573
574
575
576
577
578
579
580
/**
 * @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.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
581
void prepare_dm_gparts(struct gpart* const gparts, size_t Ndm) {
582
583

  /* Let's give all these gparts a negative id */
584
  for (size_t i = 0; i < Ndm; ++i) {
Matthieu Schaller's avatar
typo    
Matthieu Schaller committed
585
    /* 0 or negative ids are not allowed */
586
    if (gparts[i].id_or_neg_offset <= 0)
587
      error("0 or negative ID for DM particle %zu: ID=%lld", i,
588
            gparts[i].id_or_neg_offset);
589
590
591
592
593
594
595
596
597
598
599
  }
}

/**
 * @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
600
 * at the start
601
602
603
 * @param Ngas The number of gas particles read in.
 * @param Ndm The number of DM particles read in.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
604
605
606
void duplicate_hydro_gparts(struct part* const parts,
                            struct gpart* const gparts, size_t Ngas,
                            size_t Ndm) {
607

608
  for (size_t i = 0; i < Ngas; ++i) {
609
610
611
612
613
614
615
616
617
618

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

619
    gparts[i + Ndm].mass = hydro_get_mass(&parts[i]);
620
621

    /* Link the particles */
622
    gparts[i + Ndm].id_or_neg_offset = -i;
623
624
625
626
    parts[i].gpart = &gparts[i + Ndm];
  }
}

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
/**
 * @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.
 */
void duplicate_star_gparts(struct spart* const sparts,
                           struct gpart* const gparts, size_t Nstars,
                           size_t Ndm) {

  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;

    /* Link the particles */
    gparts[i + Ndm].id_or_neg_offset = -i;
    sparts[i].gpart = &gparts[i + Ndm];
  }
}

662
663
664
665
666
667
668
669
/**
 * @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.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
670
671
void collect_dm_gparts(const struct gpart* const gparts, size_t Ntot,
                       struct gpart* const dmparts, size_t Ndm) {
672

673
  size_t count = 0;
674
675

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

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

681
    /* And collect the DM ones */
682
    if (gparts[i].id_or_neg_offset > 0) {
683
      dmparts[count] = gparts[i];
684
685
686
687
688
689
      count++;
    }
  }

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

694
#endif