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 "black_holes_io.h"
35
#include "cooling.h"
36
#include "engine.h"
37
#include "gravity_io.h"
James Willis's avatar
James Willis committed
38
#include "hydro.h"
39
40
#include "hydro_io.h"
#include "stars_io.h"
41
#include "swift_velociraptor_part.h"
42
#include "threadpool.h"
43
#include "velociraptor_struct.h"
44
45
46

#ifdef HAVE_VELOCIRAPTOR

47
48
/**
 * @brief Structure for passing cosmological information to VELOCIraptor.
49
50
51
 *
 * This should match the structure cosmoinfo in the file src/swiftinterface.h
 * in the VELOCIraptor code.
52
 */
53
54
55
56
57
58
59
60
61
62
63
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;

64
65
66
67
68
69
70
71
72
  /*! 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;

73
74
75
76
77
78
79
80
81
82
83
84
85
  /*! 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;
};

86
87
/**
 * @brief Structure for passing unit information to VELOCIraptor.
88
89
90
 *
 * This should match the structure unitinfo in the file src/swiftinterface.h
 * in the VELOCIraptor code.
91
 */
92
93
struct unitinfo {

94
  /*! Length conversion factor to kpc. */
95
96
  double lengthtokpc;

97
  /*! Velocity conversion factor to km/s. */
98
99
  double velocitytokms;

100
  /*! Mass conversion factor to solar masses. */
101
102
  double masstosolarmass;

103
  /*! Potential conversion factor to (km/s)^2. */
104
105
106
107
108
109
110
111
112
  double energyperunitmass;

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

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

113
114
115
/**
 * @brief Structure to hold the location of a top-level cell.
 */
116
117
struct cell_loc {

118
  /*! Coordinates x,y,z */
119
120
121
  double loc[3];
};

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

131
132
133
134
135
136
137
138
139
140
141
142
143
  /*! 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. */
144
145
  int numcells;

146
147
148
  /*! Number of top-level cells. */
  int numcellsperdim;

149
150
151
152
153
154
155
156
157
  /*! 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];

158
159
160
  /*! Holds the node ID of each top-level cell. */
  int *cellnodeids;

161
  /*! Is this a cosmological simulation? */
162
  int icosmologicalsim;
163
164

  /*! Is this a zoom-in simulation? */
165
  int izoomsim;
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180

  /*! 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;
181
182
183
184
185

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

188
189
190
/**
 * @brief Structure for group information back to swift
 */
191
struct groupinfo {
192
193

  /*! Index of a #gpart in the global array on this MPI rank */
194
  int index;
195
196

  /*! Group number of the #gpart. */
197
  long long groupID;
198
199
};

pelahi's avatar
pelahi committed
200
201
202
int InitVelociraptor(char *config_name, struct unitinfo unit_info,
                     struct siminfo sim_info, const int numthreads);

203
struct groupinfo *InvokeVelociraptor(
204
205
    const int snapnum, char *output_name, struct cosmoinfo cosmo_info,
    struct siminfo sim_info, const size_t num_gravity_parts,
206
207
208
209
    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);
210
211

#endif /* HAVE_VELOCIRAPTOR */
212

213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
/**
 * @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;
237
238
  struct swift_vel_part *swift_parts =
      data->swift_parts + (ptrdiff_t)(gparts - s->gparts);
239
  const ptrdiff_t index_offset = gparts - s->gparts;
240
241
242
243
244

  /* Handle on the other particle types */
  const struct part *parts = s->parts;
  const struct xpart *xparts = s->xparts;
  const struct spart *sparts = s->sparts;
245
  const struct bpart *bparts = s->bparts;
246
247
248
249
250
251
252
253
254
255

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

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

264
#ifndef HAVE_VELOCIRAPTOR_WITH_NOMASS
265
    swift_parts[i].mass = gravity_get_mass(&gparts[i]);
266
267
#endif

268
269
270
271
    swift_parts[i].potential = gravity_get_comoving_potential(&gparts[i]);

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

272
    swift_parts[i].index = i + index_offset;
273
274
275
276
277
278
279
280
281
282
283
284
285
286
#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];

287
288
        convert_part_pos(e, p, xp, swift_parts[i].x);
        convert_part_vel(e, p, xp, swift_parts[i].v);
289
290
291
292
293
294
        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;

295
296
      case swift_type_stars: {
        const struct spart *sp = &sparts[-gparts[i].id_or_neg_offset];
297

298
299
        convert_spart_pos(e, sp, swift_parts[i].x);
        convert_spart_vel(e, sp, swift_parts[i].v);
300
301
302
        swift_parts[i].id = sparts[-gparts[i].id_or_neg_offset].id;
        swift_parts[i].u = 0.f;
        swift_parts[i].T = 0.f;
303
      } break;
304

305
306
      case swift_type_black_hole: {
        const struct bpart *bp = &bparts[-gparts[i].id_or_neg_offset];
307

308
309
        convert_bpart_pos(e, bp, swift_parts[i].x);
        convert_bpart_vel(e, bp, swift_parts[i].v);
310
311
312
        swift_parts[i].id = bparts[-gparts[i].id_or_neg_offset].id;
        swift_parts[i].u = 0.f;
        swift_parts[i].T = 0.f;
313
      } break;
314

315
316
      case swift_type_dark_matter:

317
318
        convert_gpart_pos(e, &(gparts[i]), swift_parts[i].x);
        convert_gpart_vel(e, &(gparts[i]), swift_parts[i].v);
319
320
321
322
323
        swift_parts[i].id = gparts[i].id_or_neg_offset;
        swift_parts[i].u = 0.f;
        swift_parts[i].T = 0.f;
        break;

324
325
      case swift_type_dark_matter_background:

326
327
        convert_gpart_pos(e, &(gparts[i]), swift_parts[i].x);
        convert_gpart_vel(e, &(gparts[i]), swift_parts[i].v);
328
329
330
331
332
        swift_parts[i].id = gparts[i].id_or_neg_offset;
        swift_parts[i].u = 0.f;
        swift_parts[i].T = 0.f;
        break;

333
334
335
336
337
338
      default:
        error("Particle type not handled by VELOCIraptor.");
    }
  }
}

339
/**
340
341
 * @brief Initialise VELOCIraptor with configuration, units,
 * simulation info needed to run.
342
343
344
 *
 * @param e The #engine.
 */
345
void velociraptor_init(struct engine *e) {
346

347
#ifdef HAVE_VELOCIRAPTOR
348
349
  const ticks tic = getticks();

350
  /* Internal SWIFT units */
351
352
  const struct unit_system *swift_us = e->internal_units;

353
  /* CGS units and physical constants in CGS */
354
355
356
357
358
  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
359
  /* Set unit conversions. */
360
  struct unitinfo unit_info;
361
362
363
364
365
366
367
368
369
370
371
  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
372
  unit_info.gravity = e->physical_constants->const_newton_G;
James Willis's avatar
James Willis committed
373
  unit_info.hubbleunit = e->cosmology->H0 / e->cosmology->h;
Matthieu Schaller's avatar
Matthieu Schaller committed
374

375
376
  /* Gather some information about the simulation */
  struct siminfo sim_info;
Matthieu Schaller's avatar
Matthieu Schaller committed
377

378
379
  /* Are we running with cosmology? */
  if (e->policy & engine_policy_cosmology) {
Matthieu Schaller's avatar
Matthieu Schaller committed
380
    sim_info.icosmologicalsim = 1;
381
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
382
    sim_info.icosmologicalsim = 0;
383
  }
384
385
386
387
388
389
390

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

  /* 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);
396
  sim_info.ibh = (e->policy & engine_policy_black_holes);
397
  sim_info.iother = 0;
398
399
400
401
402
403
404
405
406

  /* 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);
407
    message("VELOCIraptor conf: Internal energy conversion factor: %e",
408
409
410
411
412
413
414
            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
415
416

  /* Initialise VELOCIraptor. */
417
418
  if (InitVelociraptor(e->stf_config_file_name, unit_info, sim_info,
                       e->nr_threads) != 1)
419
    error("VELOCIraptor initialisation failed.");
420
421
422
423

  if (e->verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
424
#else
425
  error("SWIFT not configured to run with VELOCIraptor.");
426
#endif /* HAVE_VELOCIRAPTOR */
427
428
429
430
431
432
}

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

437
#ifdef HAVE_VELOCIRAPTOR
438

439
  /* Handle on the particles */
440
  const struct space *s = e->s;
Matthieu Schaller's avatar
Matthieu Schaller committed
441
  const size_t nr_gparts = s->nr_gparts;
442
443
  const size_t nr_parts = s->nr_parts;
  const size_t nr_sparts = s->nr_sparts;
Matthieu Schaller's avatar
Matthieu Schaller committed
444
  const int nr_cells = s->nr_cells;
445
  const struct cell *cells_top = s->cells_top;
446

447
448
449
450
451
452
  /* 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
453
  /* Allow thread to run on any core for the duration of the call to
454
455
   * VELOCIraptor so that  when OpenMP threads are spawned
   * they can run on any core on the processor. */
Matthieu Schaller's avatar
Matthieu Schaller committed
456
457
458
459
  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. */
460
  cpu_set_t cpuset;
Matthieu Schaller's avatar
Matthieu Schaller committed
461
462
463
464
  CPU_ZERO(&cpuset);
  for (int j = 0; j < nr_cores; j++) CPU_SET(j, &cpuset);
  pthread_setaffinity_np(thread, sizeof(cpu_set_t), &cpuset);

465
466
  /* Set cosmology information for this point in time */
  struct cosmoinfo cosmo_info;
467
468
469
470
  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;
471
472
473
  cosmo_info.Omega_r = e->cosmology->Omega_r;
  cosmo_info.Omega_k = e->cosmology->Omega_k;
  cosmo_info.Omega_nu = 0.;
474
475
476
477
  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;

478
479
480
481
482
483
484
485
486
487
488
489
490
  /* 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;
491

492
  /* Period of the box (Note we assume a cubic box!) */
493
  if (e->s->periodic) {
494
    sim_info.period = s->dim[0];
495
  } else {
496
497
    sim_info.period = 0.0;
  }
498
499

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

502
  /* Are we running with cosmology? */
503
504
  if (e->policy & engine_policy_cosmology) {
    sim_info.icosmologicalsim = 1;
505
506
507
508
509
510
511
512

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

513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
    /* 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

532
533
    const double Omega_m = e->cosmology->Omega_m;
    const double Omega_b = e->cosmology->Omega_b;
534
    const double critical_density_0 = e->cosmology->critical_density_0;
535

536
537
538
    /* 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. */
539
    double mean_matter_density;
540
541
542
543
544
    if (s->with_hydro)
      mean_matter_density = (Omega_m - Omega_b) * critical_density_0;
    else
      mean_matter_density = Omega_m * critical_density_0;

545
546
547
    sim_info.interparticlespacing =
        cbrt(high_res_DM_mass / mean_matter_density);

548
  } else {
549
    sim_info.izoomsim = 0;
550
551
    sim_info.icosmologicalsim = 0;
    sim_info.interparticlespacing = -1.;
552
553
  }

554
555
556
557
558
559
560
561
562
563
564
565
566
567
#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

568
  /* Set the spatial extent of the simulation volume */
569
570
571
  sim_info.spacedimension[0] = s->dim[0];
  sim_info.spacedimension[1] = s->dim[1];
  sim_info.spacedimension[2] = s->dim[2];
572

573
  /* Store number of top-level cells */
574
  sim_info.numcells = s->nr_cells;
575
  sim_info.numcellsperdim = s->cdim[0]; /* We assume a cubic box! */
576
577
  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");
578
579

  /* Size and inverse size of the top-level cells in VELOCIraptor units */
580
581
582
583
584
585
  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];
586

587
588
  ticks tic = getticks();

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

599
  for (int i = 0; i < s->nr_cells; i++) {
600
601
    cell_node_ids[i] = cells_top[i].nodeID;

Matthieu Schaller's avatar
Matthieu Schaller committed
602
603
604
605
606
607
608
609
610
611
612
613
    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];
    }
614
  }
615

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

631
632
633
634
635
  /* Report timing */
  if (e->verbose)
    message("VR Collecting top-level cell info took %.3f %s.",
            clocks_from_ticks(getticks() - tic), clocks_getunit());

636
  /* Mention the number of particles being sent */
637
  if (e->verbose)
638
639
640
    message(
        "VELOCIraptor conf: MPI rank %d sending %zu gparts to VELOCIraptor.",
        engine_rank, nr_gparts);
Matthieu Schaller's avatar
Matthieu Schaller committed
641

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

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

675
  /* What should the filename be? */
676
677
678
679
680
681
  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!");
  }
682

683
684
  tic = getticks();

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

692
693
  struct velociraptor_copy_data copy_data = {e, swift_parts};
  threadpool_map(&e->threadpool, velociraptor_convert_particles_mapper,
694
695
                 s->gparts, nr_gparts, sizeof(struct gpart),
                 threadpool_auto_chunk_size, &copy_data);
696

697
698
699
700
701
702
703
  /* Report timing */
  if (e->verbose)
    message("VR Collecting particle info took %.3f %s.",
            clocks_from_ticks(getticks() - tic), clocks_getunit());

  tic = getticks();

704
705
706
707
  /* Values returned by VELOCIRaptor */
  int num_gparts_in_groups = -1;
  struct groupinfo *group_info = NULL;

708
709
710
711
712
713
714
#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
715
  /* Call VELOCIraptor. */
716
  group_info = (struct groupinfo *)InvokeVelociraptor(
717
718
719
      e->stf_output_count, outputFileName, cosmo_info, sim_info, nr_gparts,
      nr_parts, nr_sparts, swift_parts, cell_node_ids, e->nr_threads,
      linked_with_snap, &num_gparts_in_groups);
720

721
722
723
724
725
  /* 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);

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

734
735
  /* Report timing */
  if (e->verbose)
736
    message("VR Invocation of velociraptor took %.3f %s.",
737
738
739
740
            clocks_from_ticks(getticks() - tic), clocks_getunit());

  tic = getticks();

741
  /* Assign the group IDs back to the gparts */
742
  if (linked_with_snap) {
743

744
745
    if (swift_memalign("VR.group_data", (void **)&s->gpart_group_data,
                       part_align,
746
747
748
749
750
                       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;

751
    /* Zero the array (gparts not in groups have an ID of 0) */
752
753
754
    bzero(data, nr_gparts * sizeof(struct velociraptor_gpart_data));

    /* Copy the data at the right place */
755
    for (int i = 0; i < num_gparts_in_groups; i++) {
756
      data[group_info[i].index].groupID = group_info[i].groupID;
757
758
    }

759
760
761
762
763
    /* Report timing */
    if (e->verbose)
      message("VR Copying group information back took %.3f %s.",
              clocks_from_ticks(getticks() - tic), clocks_getunit());

764
    /* Free the array returned by VELOCIraptor */
765
    swift_free("VR.group_data", group_info);
766
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
767

768
769
770
  /* Reset the pthread affinity mask after VELOCIraptor returns. */
  pthread_setaffinity_np(thread, sizeof(cpu_set_t), engine_entry_affinity());

771
  /* Increase output counter (if not linked with snapshot) */
772
  if (!linked_with_snap) e->stf_output_count++;
773

774
775
  /* Record we have ran stf this timestep */
  e->stf_this_timestep = 1;
776

777
778
779
780
781
782
  /* Reallocate the memory that was freed earlier */
#ifdef WITH_MPI

  engine_allocate_foreign_particles(e);
#endif

783
#else
784
  error("SWIFT not configured to run with VELOCIraptor.");
785
#endif /* HAVE_VELOCIRAPTOR */
786
}