velociraptor_interface.c 25.3 KB
Newer Older
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Copyright (c) 2018 James Willis (james.s.willis@durham.ac.uk)
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 *
 * 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.
 *
 * 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.
 *
 * 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/>.
 *
 ******************************************************************************/

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

23
/* Some standard headers. */
24
25
26
27
#include <stdio.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/types.h>
28
29
#include <unistd.h>

30
31
32
/* This object's header. */
#include "velociraptor_interface.h"

33
/* Local includes. */
34
#include "cooling.h"
35
#include "engine.h"
James Willis's avatar
James Willis committed
36
#include "hydro.h"
37
#include "swift_velociraptor_part.h"
38
#include "threadpool.h"
39
#include "velociraptor_struct.h"
40
41
42

#ifdef HAVE_VELOCIRAPTOR

43
44
/**
 * @brief Structure for passing cosmological information to VELOCIraptor.
45
46
47
 *
 * This should match the structure cosmoinfo in the file src/swiftinterface.h
 * in the VELOCIraptor code.
48
 */
49
50
51
52
53
54
55
56
57
58
59
struct cosmoinfo {

  /*! Current expansion factor of the Universe. (cosmology.a) */
  double atime;

  /*! Reduced Hubble constant (H0 / (100km/s/Mpc) (cosmology.h) */
  double littleh;

  /*! Matter density parameter (cosmology.Omega_m) */
  double Omega_m;

60
61
62
63
64
65
66
67
68
  /*! Radiation density parameter (cosmology.Omega_r) */
  double Omega_r;

  /*! Neutrino density parameter (0 in SWIFT) */
  double Omega_nu;

  /*! Neutrino density parameter (cosmology.Omega_k) */
  double Omega_k;

69
70
71
72
73
74
75
76
77
78
79
80
81
  /*! Baryon density parameter (cosmology.Omega_b) */
  double Omega_b;

  /*! Radiation constant density parameter (cosmology.Omega_lambda) */
  double Omega_Lambda;

  /*! Dark matter density parameter (cosmology.Omega_m - cosmology.Omega_b) */
  double Omega_cdm;

  /*! Dark-energy equation of state at the current time (cosmology.w)*/
  double w_de;
};

82
83
/**
 * @brief Structure for passing unit information to VELOCIraptor.
84
85
86
 *
 * This should match the structure unitinfo in the file src/swiftinterface.h
 * in the VELOCIraptor code.
87
 */
88
89
struct unitinfo {

90
  /*! Length conversion factor to kpc. */
91
92
  double lengthtokpc;

93
  /*! Velocity conversion factor to km/s. */
94
95
  double velocitytokms;

96
  /*! Mass conversion factor to solar masses. */
97
98
  double masstosolarmass;

99
  /*! Potential conversion factor to (km/s)^2. */
100
101
102
103
104
105
106
107
108
  double energyperunitmass;

  /*! Newton's gravitationl constant (phys_const.const_newton_G)*/
  double gravity;

  /*! Hubble constant at the current redshift (cosmology.H) */
  double hubbleunit;
};

109
110
111
/**
 * @brief Structure to hold the location of a top-level cell.
 */
112
113
struct cell_loc {

114
  /*! Coordinates x,y,z */
115
116
117
  double loc[3];
};

118
119
120
/**
 * @brief Structure for passing simulation information to VELOCIraptor for a
 * given call.
121
122
123
 *
 * This should match the structure siminfo in the file src/swiftinterface.h
 * in the VELOCIraptor code.
124
 */
125
126
struct siminfo {

127
128
129
130
131
132
133
134
135
136
137
138
139
  /*! Size of periodic replications */
  double period;

  /*! Mass of the high-resolution DM particles in a zoom-in run. */
  double zoomhigresolutionmass;

  /*! Mean inter-particle separation of the DM particles */
  double interparticlespacing;

  /*! Spacial extent of the simulation volume */
  double spacedimension[3];

  /*! Number of top-level cells. */
140
141
  int numcells;

142
143
144
  /*! Number of top-level cells. */
  int numcellsperdim;

145
146
147
148
149
150
151
152
153
  /*! Locations of top-level cells. */
  struct cell_loc *cell_loc;

  /*! Top-level cell width. */
  double cellwidth[3];

  /*! Inverse of the top-level cell width. */
  double icellwidth[3];

154
155
156
  /*! Holds the node ID of each top-level cell. */
  int *cellnodeids;

157
  /*! Is this a cosmological simulation? */
158
  int icosmologicalsim;
159
160

  /*! Is this a zoom-in simulation? */
161
  int izoomsim;
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176

  /*! Do we have DM particles? */
  int idarkmatter;

  /*! Do we have gas particles? */
  int igas;

  /*! Do we have star particles? */
  int istar;

  /*! Do we have BH particles? */
  int ibh;

  /*! Do we have other particles? */
  int iother;
177
178
179
180
181

#ifdef HAVE_VELOCIRAPTOR_WITH_NOMASS
  /*! Mass of the DM particles */
  double mass_uniform_box;
#endif
182
183
};

184
185
186
/**
 * @brief Structure for group information back to swift
 */
187
struct groupinfo {
188
189

  /*! Index of a #gpart in the global array on this MPI rank */
190
  int index;
191
192

  /*! Group number of the #gpart. */
193
  long long groupID;
194
195
};

pelahi's avatar
pelahi committed
196
197
198
int InitVelociraptor(char *config_name, struct unitinfo unit_info,
                     struct siminfo sim_info, const int numthreads);

199
struct groupinfo *InvokeVelociraptor(
200
201
    const int snapnum, char *output_name, struct cosmoinfo cosmo_info,
    struct siminfo sim_info, const size_t num_gravity_parts,
202
203
204
205
    const size_t num_hydro_parts, const size_t num_star_parts,
    struct swift_vel_part *swift_parts, const int *cell_node_ids,
    const int numthreads, const int return_group_flags,
    int *const num_in_groups);
206
207

#endif /* HAVE_VELOCIRAPTOR */
208

209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
/**
 * @brief Temporary structure used for the data copy mapper.
 */
struct velociraptor_copy_data {
  const struct engine *e;
  struct swift_vel_part *swift_parts;
};

/**
 * @brief Mapper function to conver the #gpart into VELOCIraptor Particles.
 *
 * @param map_data The array of #gpart.
 * @param nr_gparts The number of #gpart.
 * @param extra_data Pointer to the #engine and to the array to fill.
 */
void velociraptor_convert_particles_mapper(void *map_data, int nr_gparts,
                                           void *extra_data) {

  /* Unpack the data */
  struct gpart *restrict gparts = (struct gpart *)map_data;
  struct velociraptor_copy_data *data =
      (struct velociraptor_copy_data *)extra_data;
  const struct engine *e = data->e;
  const struct space *s = e->s;
233
234
  struct swift_vel_part *swift_parts =
      data->swift_parts + (ptrdiff_t)(gparts - s->gparts);
235
236
237
238
239

  /* Handle on the other particle types */
  const struct part *parts = s->parts;
  const struct xpart *xparts = s->xparts;
  const struct spart *sparts = s->sparts;
240
  const struct bpart *bparts = s->bparts;
241
242
243
244
245
246
247
248
249

  /* Handle on the physics modules */
  const struct cosmology *cosmo = e->cosmology;
  const struct hydro_props *hydro_props = e->hydro_properties;
  const struct unit_system *us = e->internal_units;
  const struct phys_const *phys_const = e->physical_constants;
  const struct cooling_function_data *cool_func = e->cooling_func;

  const float a_inv = e->cosmology->a_inv;
Matthieu Schaller's avatar
Matthieu Schaller committed
250
251
252
253
  const int periodic = s->periodic;
  const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
  const double pos_dithering[3] = {s->pos_dithering[0], s->pos_dithering[1],
                                   s->pos_dithering[2]};
254
255
256

  /* Convert particle properties into VELOCIraptor units.
   * VELOCIraptor wants:
Matthieu Schaller's avatar
Matthieu Schaller committed
257
   * - Un-dithered co-moving positions,
258
259
260
261
262
263
264
   * - Peculiar velocities,
   * - Co-moving potential,
   * - Physical internal energy (for the gas),
   * - Temperatures (for the gas).
   */
  for (int i = 0; i < nr_gparts; i++) {

Matthieu Schaller's avatar
Matthieu Schaller committed
265
266
267
268
269
270
271
272
273
274
275
276
    if (periodic) {
      swift_parts[i].x[0] =
          box_wrap(gparts[i].x[0] - pos_dithering[0], 0.0, dim[0]);
      swift_parts[i].x[1] =
          box_wrap(gparts[i].x[1] - pos_dithering[1], 0.0, dim[1]);
      swift_parts[i].x[2] =
          box_wrap(gparts[i].x[2] - pos_dithering[2], 0.0, dim[2]);
    } else {
      swift_parts[i].x[0] = gparts[i].x[0];
      swift_parts[i].x[1] = gparts[i].x[1];
      swift_parts[i].x[2] = gparts[i].x[2];
    }
277
278
279
280
281

    swift_parts[i].v[0] = gparts[i].v_full[0] * a_inv;
    swift_parts[i].v[1] = gparts[i].v_full[1] * a_inv;
    swift_parts[i].v[2] = gparts[i].v_full[2] * a_inv;

282
#ifndef HAVE_VELOCIRAPTOR_WITH_NOMASS
283
    swift_parts[i].mass = gravity_get_mass(&gparts[i]);
284
285
#endif

286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
    swift_parts[i].potential = gravity_get_comoving_potential(&gparts[i]);

    swift_parts[i].type = gparts[i].type;

    swift_parts[i].index = i;
#ifdef WITH_MPI
    swift_parts[i].task = e->nodeID;
#else
    swift_parts[i].task = 0;
#endif

    /* Set gas particle IDs from their hydro counterparts and set internal
     * energies. */
    switch (gparts[i].type) {

      case swift_type_gas: {
        const struct part *p = &parts[-gparts[i].id_or_neg_offset];
        const struct xpart *xp = &xparts[-gparts[i].id_or_neg_offset];

        swift_parts[i].id = parts[-gparts[i].id_or_neg_offset].id;
        swift_parts[i].u = hydro_get_drifted_physical_internal_energy(p, cosmo);
        swift_parts[i].T = cooling_get_temperature(phys_const, hydro_props, us,
                                                   cosmo, cool_func, p, xp);
      } break;

      case swift_type_stars:

        swift_parts[i].id = sparts[-gparts[i].id_or_neg_offset].id;
        swift_parts[i].u = 0.f;
        swift_parts[i].T = 0.f;
        break;

318
319
320
321
322
323
324
      case swift_type_black_hole:

        swift_parts[i].id = bparts[-gparts[i].id_or_neg_offset].id;
        swift_parts[i].u = 0.f;
        swift_parts[i].T = 0.f;
        break;

325
326
327
328
329
330
331
      case swift_type_dark_matter:

        swift_parts[i].id = gparts[i].id_or_neg_offset;
        swift_parts[i].u = 0.f;
        swift_parts[i].T = 0.f;
        break;

332
333
334
335
336
337
338
      case swift_type_dark_matter_background:

        swift_parts[i].id = gparts[i].id_or_neg_offset;
        swift_parts[i].u = 0.f;
        swift_parts[i].T = 0.f;
        break;

339
340
341
342
343
344
      default:
        error("Particle type not handled by VELOCIraptor.");
    }
  }
}

345
/**
346
347
 * @brief Initialise VELOCIraptor with configuration, units,
 * simulation info needed to run.
348
349
350
 *
 * @param e The #engine.
 */
351
void velociraptor_init(struct engine *e) {
352

353
#ifdef HAVE_VELOCIRAPTOR
354
355
  const ticks tic = getticks();

356
  /* Internal SWIFT units */
357
358
  const struct unit_system *swift_us = e->internal_units;

359
  /* CGS units and physical constants in CGS */
360
361
362
363
364
  struct unit_system cgs_us;
  units_init_cgs(&cgs_us);
  struct phys_const cgs_pc;
  phys_const_init(&cgs_us, /*params=*/NULL, &cgs_pc);

Matthieu Schaller's avatar
Matthieu Schaller committed
365
  /* Set unit conversions. */
366
  struct unitinfo unit_info;
367
368
369
370
371
372
373
374
375
376
377
  unit_info.lengthtokpc =
      units_cgs_conversion_factor(swift_us, UNIT_CONV_LENGTH) /
      (1000. * cgs_pc.const_parsec);
  unit_info.velocitytokms =
      units_cgs_conversion_factor(swift_us, UNIT_CONV_VELOCITY) / 1.0e5;
  unit_info.masstosolarmass =
      units_cgs_conversion_factor(swift_us, UNIT_CONV_MASS) /
      cgs_pc.const_solar_mass;
  unit_info.energyperunitmass =
      units_cgs_conversion_factor(swift_us, UNIT_CONV_ENERGY_PER_UNIT_MASS) /
      (1.0e10);
Matthieu Schaller's avatar
Matthieu Schaller committed
378
  unit_info.gravity = e->physical_constants->const_newton_G;
James Willis's avatar
James Willis committed
379
  unit_info.hubbleunit = e->cosmology->H0 / e->cosmology->h;
Matthieu Schaller's avatar
Matthieu Schaller committed
380

381
382
  /* Gather some information about the simulation */
  struct siminfo sim_info;
Matthieu Schaller's avatar
Matthieu Schaller committed
383

384
385
  /* Are we running with cosmology? */
  if (e->policy & engine_policy_cosmology) {
Matthieu Schaller's avatar
Matthieu Schaller committed
386
    sim_info.icosmologicalsim = 1;
387
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
388
    sim_info.icosmologicalsim = 0;
389
  }
390
391
392
393
394
395
396

  /* Are we running a zoom? */
  if (e->s->with_DM_background) {
    sim_info.izoomsim = 1;
  } else {
    sim_info.izoomsim = 0;
  }
397
398
399
400
401

  /* Tell VELOCIraptor what we have in the simulation */
  sim_info.idarkmatter = (e->total_nr_gparts - e->total_nr_parts > 0);
  sim_info.igas = (e->policy & engine_policy_hydro);
  sim_info.istar = (e->policy & engine_policy_stars);
402
  sim_info.ibh = (e->policy & engine_policy_black_holes);
403
  sim_info.iother = 0;
404
405
406
407
408
409
410
411
412

  /* Be nice, talk! */
  if (e->verbose) {
    message("VELOCIraptor conf: Length conversion factor: %e",
            unit_info.lengthtokpc);
    message("VELOCIraptor conf: Velocity conversion factor: %e",
            unit_info.velocitytokms);
    message("VELOCIraptor conf: Mass conversion factor: %e",
            unit_info.masstosolarmass);
413
    message("VELOCIraptor conf: Internal energy conversion factor: %e",
414
415
416
417
418
419
420
            unit_info.energyperunitmass);
    message("VELOCIraptor conf: G: %e", unit_info.gravity);
    message("VELOCIraptor conf: H0/h: %e", unit_info.hubbleunit);
    message("VELOCIraptor conf: Config file name: %s", e->stf_config_file_name);
    message("VELOCIraptor conf: Cosmological Simulation: %d",
            sim_info.icosmologicalsim);
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
421
422

  /* Initialise VELOCIraptor. */
423
424
  if (InitVelociraptor(e->stf_config_file_name, unit_info, sim_info,
                       e->nr_threads) != 1)
425
    error("VELOCIraptor initialisation failed.");
426
427
428
429

  if (e->verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
430
#else
431
  error("SWIFT not configured to run with VELOCIraptor.");
432
#endif /* HAVE_VELOCIRAPTOR */
433
434
435
436
437
438
}

/**
 * @brief Run VELOCIraptor with current particle data.
 *
 * @param e The #engine.
439
 * @param linked_with_snap Are we running at the same time as a snapshot dump?
440
 */
441
void velociraptor_invoke(struct engine *e, const int linked_with_snap) {
442

443
#ifdef HAVE_VELOCIRAPTOR
444

445
  /* Handle on the particles */
446
  const struct space *s = e->s;
Matthieu Schaller's avatar
Matthieu Schaller committed
447
  const size_t nr_gparts = s->nr_gparts;
448
449
  const size_t nr_parts = s->nr_parts;
  const size_t nr_sparts = s->nr_sparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
450
  const int nr_cells = s->nr_cells;
451
  const struct cell *cells_top = s->cells_top;
452

453
454
455
456
457
458
  /* Start by freeing some of the unnecessary memory to give VR some breathing
     space */
#ifdef WITH_MPI
  space_free_foreign_parts(e->s, /*clear_cell_pointers=*/1);
#endif

Matthieu Schaller's avatar
Matthieu Schaller committed
459
  /* Allow thread to run on any core for the duration of the call to
460
461
   * VELOCIraptor so that  when OpenMP threads are spawned
   * they can run on any core on the processor. */
Matthieu Schaller's avatar
Matthieu Schaller committed
462
463
464
465
  const int nr_cores = sysconf(_SC_NPROCESSORS_ONLN);
  pthread_t thread = pthread_self();

  /* Set affinity mask to include all cores on the CPU for VELOCIraptor. */
466
  cpu_set_t cpuset;
Matthieu Schaller's avatar
Matthieu Schaller committed
467
468
469
470
  CPU_ZERO(&cpuset);
  for (int j = 0; j < nr_cores; j++) CPU_SET(j, &cpuset);
  pthread_setaffinity_np(thread, sizeof(cpu_set_t), &cpuset);

471
472
  /* Set cosmology information for this point in time */
  struct cosmoinfo cosmo_info;
473
474
475
476
  cosmo_info.atime = e->cosmology->a;
  cosmo_info.littleh = e->cosmology->h;
  cosmo_info.Omega_m = e->cosmology->Omega_m;
  cosmo_info.Omega_b = e->cosmology->Omega_b;
477
478
479
  cosmo_info.Omega_r = e->cosmology->Omega_r;
  cosmo_info.Omega_k = e->cosmology->Omega_k;
  cosmo_info.Omega_nu = 0.;
480
481
482
483
  cosmo_info.Omega_Lambda = e->cosmology->Omega_lambda;
  cosmo_info.Omega_cdm = e->cosmology->Omega_m - e->cosmology->Omega_b;
  cosmo_info.w_de = e->cosmology->w;

484
485
486
487
488
489
490
491
492
493
494
495
496
  /* Report the cosmo info we use */
  if (e->verbose) {
    message("VELOCIraptor conf: Scale factor: %e", cosmo_info.atime);
    message("VELOCIraptor conf: Little h: %e", cosmo_info.littleh);
    message("VELOCIraptor conf: Omega_m: %e", cosmo_info.Omega_m);
    message("VELOCIraptor conf: Omega_b: %e", cosmo_info.Omega_b);
    message("VELOCIraptor conf: Omega_Lambda: %e", cosmo_info.Omega_Lambda);
    message("VELOCIraptor conf: Omega_cdm: %e", cosmo_info.Omega_cdm);
    message("VELOCIraptor conf: w_de: %e", cosmo_info.w_de);
  }

  /* Update the simulation information */
  struct siminfo sim_info;
497

498
  /* Period of the box (Note we assume a cubic box!) */
499
  if (e->s->periodic) {
500
    sim_info.period = s->dim[0];
501
  } else {
502
503
    sim_info.period = 0.0;
  }
504
505

  /* Tell VELOCIraptor this is not a zoom-in simulation */
506
507
  sim_info.zoomhigresolutionmass = -1.0;

508
  /* Are we running with cosmology? */
509
510
  if (e->policy & engine_policy_cosmology) {
    sim_info.icosmologicalsim = 1;
511
512
513
514
515
516
517
518

    /* Are we running a zoom? */
    if (e->s->with_DM_background) {
      sim_info.izoomsim = 1;
    } else {
      sim_info.izoomsim = 0;
    }

519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
    /* Collect the mass of the non-background gpart */
    double high_res_DM_mass = 0.;
    for (size_t i = 0; i < e->s->nr_gparts; ++i) {
      const struct gpart *gp = &e->s->gparts[i];
      if (gp->type == swift_type_dark_matter &&
          gp->time_bin != time_bin_inhibited &&
          gp->time_bin != time_bin_not_created) {
        high_res_DM_mass = gp->mass;
        break;
      }
    }

#ifdef WITH_MPI
    /* We need to all-reduce this in case one of the nodes had 0 DM particles.
     */
    MPI_Allreduce(MPI_IN_PLACE, &high_res_DM_mass, 1, MPI_DOUBLE, MPI_MAX,
                  MPI_COMM_WORLD);
#endif

538
539
    const double Omega_m = e->cosmology->Omega_m;
    const double Omega_b = e->cosmology->Omega_b;
540
    const double critical_density_0 = e->cosmology->critical_density_0;
541

542
543
544
    /* Linking length based on the mean DM inter-particle separation
     * in the zoom region and assuming the mean density of the Universe
     * is used in the zoom region. */
545
    double mean_matter_density;
546
547
548
549
550
    if (s->with_hydro)
      mean_matter_density = (Omega_m - Omega_b) * critical_density_0;
    else
      mean_matter_density = Omega_m * critical_density_0;

551
552
553
    sim_info.interparticlespacing =
        cbrt(high_res_DM_mass / mean_matter_density);

554
  } else {
555
    sim_info.izoomsim = 0;
556
557
    sim_info.icosmologicalsim = 0;
    sim_info.interparticlespacing = -1.;
558
559
  }

560
561
562
563
564
565
566
567
568
569
570
571
572
573
#ifdef HAVE_VELOCIRAPTOR_WITH_NOMASS
  /* Assume all particles have the same mass */
  double DM_mass = 0.;
  for (size_t i = 0; i < e->s->nr_gparts; ++i) {
    const struct gpart *gp = &e->s->gparts[i];
    if (gp->time_bin != time_bin_inhibited &&
        gp->time_bin != time_bin_not_created) {
      DM_mass = gp->mass;
      break;
    }
  }
  sim_info.mass_uniform_box = DM_mass;
#endif

574
  /* Set the spatial extent of the simulation volume */
575
576
577
  sim_info.spacedimension[0] = s->dim[0];
  sim_info.spacedimension[1] = s->dim[1];
  sim_info.spacedimension[2] = s->dim[2];
578

579
  /* Store number of top-level cells */
580
  sim_info.numcells = s->nr_cells;
581
  sim_info.numcellsperdim = s->cdim[0]; /* We assume a cubic box! */
582
583
  if (s->cdim[0] != s->cdim[1] || s->cdim[0] != s->cdim[2])
    error("Trying to run VR on a non-cubic number of top-level cells");
584
585

  /* Size and inverse size of the top-level cells in VELOCIraptor units */
586
587
588
589
590
591
  sim_info.cellwidth[0] = s->cells_top[0].width[0];
  sim_info.cellwidth[1] = s->cells_top[0].width[1];
  sim_info.cellwidth[2] = s->cells_top[0].width[2];
  sim_info.icellwidth[0] = s->iwidth[0];
  sim_info.icellwidth[1] = s->iwidth[1];
  sim_info.icellwidth[2] = s->iwidth[2];
592

593
594
  ticks tic = getticks();

595
596
  /* Allocate and populate array of cell node IDs and positions. */
  int *cell_node_ids = NULL;
597
598
  if (swift_memalign("VR.cell_loc", (void **)&sim_info.cell_loc,
                     SWIFT_STRUCT_ALIGNMENT,
599
600
                     s->nr_cells * sizeof(struct cell_loc)) != 0)
    error("Failed to allocate top-level cell locations for VELOCIraptor.");
601
602
  if (swift_memalign("VR.cell_nodeID", (void **)&cell_node_ids,
                     SWIFT_STRUCT_ALIGNMENT, nr_cells * sizeof(int)) != 0)
603
604
    error("Failed to allocate list of cells node IDs for VELOCIraptor.");

605
  for (int i = 0; i < s->nr_cells; i++) {
606
607
    cell_node_ids[i] = cells_top[i].nodeID;

Matthieu Schaller's avatar
Matthieu Schaller committed
608
609
610
611
612
613
614
615
616
617
618
619
    if (s->periodic) {
      sim_info.cell_loc[i].loc[0] =
          box_wrap(cells_top[i].loc[0] - s->pos_dithering[0], 0.0, s->dim[0]);
      sim_info.cell_loc[i].loc[1] =
          box_wrap(cells_top[i].loc[1] - s->pos_dithering[1], 0.0, s->dim[1]);
      sim_info.cell_loc[i].loc[2] =
          box_wrap(cells_top[i].loc[2] - s->pos_dithering[2], 0.0, s->dim[2]);
    } else {
      sim_info.cell_loc[i].loc[0] = cells_top[i].loc[0];
      sim_info.cell_loc[i].loc[1] = cells_top[i].loc[1];
      sim_info.cell_loc[i].loc[2] = cells_top[i].loc[2];
    }
620
  }
621

622
623
624
625
626
627
628
629
630
631
632
633
634
  if (e->verbose) {
    message("VELOCIraptor conf: Space dimensions: (%e,%e,%e)",
            sim_info.spacedimension[0], sim_info.spacedimension[1],
            sim_info.spacedimension[2]);
    message("VELOCIraptor conf: No. of top-level cells: %d", sim_info.numcells);
    message(
        "VELOCIraptor conf: Top-level cell locations range: (%e,%e,%e) -> "
        "(%e,%e,%e)",
        sim_info.cell_loc[0].loc[0], sim_info.cell_loc[0].loc[1],
        sim_info.cell_loc[0].loc[2],
        sim_info.cell_loc[sim_info.numcells - 1].loc[0],
        sim_info.cell_loc[sim_info.numcells - 1].loc[1],
        sim_info.cell_loc[sim_info.numcells - 1].loc[2]);
635
636
  }

637
638
639
640
641
  /* Report timing */
  if (e->verbose)
    message("VR Collecting top-level cell info took %.3f %s.",
            clocks_from_ticks(getticks() - tic), clocks_getunit());

642
  /* Mention the number of particles being sent */
643
  if (e->verbose)
644
645
646
    message(
        "VELOCIraptor conf: MPI rank %d sending %zu gparts to VELOCIraptor.",
        engine_rank, nr_gparts);
Matthieu Schaller's avatar
Matthieu Schaller committed
647

648
649
650
  /* Generate directory name for this output - start with snapshot directory, if
   * specified */
  char outputDirName[FILENAME_BUFFER_SIZE] = "";
651
  if (strcmp(e->snapshot_subdir, engine_default_snapshot_subdir) != 0) {
652
653
654
655
    if (snprintf(outputDirName, FILENAME_BUFFER_SIZE, "%s/",
                 e->snapshot_subdir) >= FILENAME_BUFFER_SIZE) {
      error("FILENAME_BUFFER_SIZE is to small for snapshot directory name!");
    }
656
    if (engine_rank == 0) io_make_snapshot_subdir(e->snapshot_subdir);
657
658
659
660
661
662
663
#ifdef WITH_MPI
    MPI_Barrier(MPI_COMM_WORLD);
#endif
  }

  /* Then create output-specific subdirectory if necessary */
  char subDirName[FILENAME_BUFFER_SIZE] = "";
664
  if (strcmp(e->stf_subdir_per_output, engine_default_stf_subdir_per_output) != 0) {
665
666
    if (snprintf(subDirName, FILENAME_BUFFER_SIZE, "%s%s_%04i/", outputDirName,
                 e->stf_subdir_per_output,
667
                 e->stf_output_count) >= FILENAME_BUFFER_SIZE) {
668
669
670
      error(
          "FILENAME_BUFFER_SIZE is to small for Velociraptor directory name!");
    }
671
    if (engine_rank == 0) io_make_snapshot_subdir(subDirName);
672
673
674
675
676
677
678
#ifdef WITH_MPI
    MPI_Barrier(MPI_COMM_WORLD);
#endif
  } else {
    /* Not making separate directories so subDirName=outputDirName */
    strncpy(subDirName, outputDirName, FILENAME_BUFFER_SIZE);
  }
679

680
  /* What should the filename be? */
681
682
683
684
685
686
  char outputFileName[FILENAME_BUFFER_SIZE];
  if (snprintf(outputFileName, FILENAME_BUFFER_SIZE, "%s%s_%04i.VELOCIraptor",
               subDirName, e->stf_base_name,
               e->stf_output_count) >= FILENAME_BUFFER_SIZE) {
    error("FILENAME_BUFFER_SIZE is too small for Velociraptor file name!");
  }
687

688
689
  tic = getticks();

Matthieu Schaller's avatar
Matthieu Schaller committed
690
691
692
  /* Allocate and populate an array of swift_vel_parts to be passed to
   * VELOCIraptor. */
  struct swift_vel_part *swift_parts = NULL;
693
  if (swift_memalign("VR.parts", (void **)&swift_parts, part_align,
Matthieu Schaller's avatar
Matthieu Schaller committed
694
695
696
                     nr_gparts * sizeof(struct swift_vel_part)) != 0)
    error("Failed to allocate array of particles for VELOCIraptor.");

697
698
  struct velociraptor_copy_data copy_data = {e, swift_parts};
  threadpool_map(&e->threadpool, velociraptor_convert_particles_mapper,
699
700
                 s->gparts, nr_gparts, sizeof(struct gpart),
                 threadpool_auto_chunk_size, &copy_data);
701

702
703
704
705
706
707
708
  /* Report timing */
  if (e->verbose)
    message("VR Collecting particle info took %.3f %s.",
            clocks_from_ticks(getticks() - tic), clocks_getunit());

  tic = getticks();

709
710
711
712
  /* Values returned by VELOCIRaptor */
  int num_gparts_in_groups = -1;
  struct groupinfo *group_info = NULL;

713
714
715
716
717
718
719
#ifdef SWIFT_MEMUSE_REPORTS
  char report_filename[60];
  sprintf(report_filename, "memuse-VR-report-rank%d-step%d.txt", e->nodeID,
          e->step);
  memuse_log_dump(report_filename);
#endif

Matthieu Schaller's avatar
Matthieu Schaller committed
720
  /* Call VELOCIraptor. */
721
  group_info = (struct groupinfo *)InvokeVelociraptor(
722
      e->stf_output_count, outputFileName, cosmo_info, sim_info, nr_gparts, nr_parts,
723
      nr_sparts, swift_parts, cell_node_ids, e->nr_threads, linked_with_snap,
724
      &num_gparts_in_groups);
725

726
727
728
729
730
  /* Report that the memory was freed */
  memuse_log_allocation("VR.cell_loc", sim_info.cell_loc, 0, 0);
  memuse_log_allocation("VR.cell_nodeID", cell_node_ids, 0, 0);
  memuse_log_allocation("VR.parts", swift_parts, 0, 0);

731
  /* Check that the ouput is valid */
732
  if (linked_with_snap && group_info == NULL && num_gparts_in_groups < 0) {
Matthieu Schaller's avatar
Matthieu Schaller committed
733
    error("Exiting. Call to VELOCIraptor failed on rank: %d.", e->nodeID);
734
  }
735
736
737
  if (!linked_with_snap && group_info != NULL) {
    error("VELOCIraptor returned an array whilst it should not have.");
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
738

739
740
  /* Report timing */
  if (e->verbose)
741
    message("VR Invocation of velociraptor took %.3f %s.",
742
743
744
745
            clocks_from_ticks(getticks() - tic), clocks_getunit());

  tic = getticks();

746
  /* Assign the group IDs back to the gparts */
747
  if (linked_with_snap) {
748

749
750
    if (swift_memalign("VR.group_data", (void **)&s->gpart_group_data,
                       part_align,
751
752
753
754
755
                       nr_gparts * sizeof(struct velociraptor_gpart_data)) != 0)
      error("Failed to allocate array of gpart data for VELOCIraptor i/o.");

    struct velociraptor_gpart_data *data = s->gpart_group_data;

756
    /* Zero the array (gparts not in groups have an ID of 0) */
757
758
759
    bzero(data, nr_gparts * sizeof(struct velociraptor_gpart_data));

    /* Copy the data at the right place */
760
    for (int i = 0; i < num_gparts_in_groups; i++) {
761
      data[group_info[i].index].groupID = group_info[i].groupID;
762
763
    }

764
765
766
767
768
    /* Report timing */
    if (e->verbose)
      message("VR Copying group information back took %.3f %s.",
              clocks_from_ticks(getticks() - tic), clocks_getunit());

769
    /* Free the array returned by VELOCIraptor */
770
    swift_free("VR.group_data", group_info);
771
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
772

773
774
775
  /* Reset the pthread affinity mask after VELOCIraptor returns. */
  pthread_setaffinity_np(thread, sizeof(cpu_set_t), engine_entry_affinity());

776
  /* Increase output counter (if not linked with snapshot) */
777
  if (!linked_with_snap) e->stf_output_count++;
778

779
780
  /* Record we have ran stf this timestep */
  e->stf_this_timestep = 1;
781

782
783
784
785
786
787
  /* Reallocate the memory that was freed earlier */
#ifdef WITH_MPI

  engine_allocate_foreign_particles(e);
#endif

788
#else
789
  error("SWIFT not configured to run with VELOCIraptor.");
790
#endif /* HAVE_VELOCIRAPTOR */
791
}