velociraptor_interface.c 13.4 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
#include <unistd.h>

26
27
28
/* This object's header. */
#include "velociraptor_interface.h"

29
/* Local includes. */
30
#include "engine.h"
James Willis's avatar
James Willis committed
31
#include "hydro.h"
32
33
34
35
#include "swift_velociraptor_part.h"

#ifdef HAVE_VELOCIRAPTOR

36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
/* Structure for passing cosmological information to VELOCIraptor. */
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;

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

/* Structure for passing unit information to VELOCIraptor. */
struct unitinfo {

  /* Length conversion factor to kpc. */
  double lengthtokpc;

  /* Velocity conversion factor to km/s. */
  double velocitytokms;

  /* Mass conversion factor to solar masses. */
  double masstosolarmass;

  /* Potential conversion factor. */
  double energyperunitmass;

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

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

/* Structure to hold the location of a top-level cell. */
struct cell_loc {

  /* Coordinates x,y,z */
  double loc[3];
};

/* Structure for passing simulation information to VELOCIraptor. */
struct siminfo {
  double period, zoomhigresolutionmass, interparticlespacing, spacedimension[3];

  /* Number of top-cells. */
  int numcells;

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

  int icosmologicalsim;
};

109
110
111
112
/* VELOCIraptor interface. */
int InitVelociraptor(char *config_name, char *output_name,
                     struct cosmoinfo cosmo_info, struct unitinfo unit_info,
                     struct siminfo sim_info);
113

114
int InvokeVelociraptor(const size_t num_gravity_parts,
115
                       const size_t num_hydro_parts, const int snapnum,
116
117
118
119
                       struct swift_vel_part *swift_parts,
                       const int *cell_node_ids, char *output_name);

#endif /* HAVE_VELOCIRAPTOR */
120

121
/**
Matthieu Schaller's avatar
Matthieu Schaller committed
122
123
 * @brief Initialise VELOCIraptor with input and output file names along with
 * cosmological info needed to run.
124
125
 *
 * @param e The #engine.
126
 * @param linked_with_snap Are we running at the same time as a snapshot dump?
127
 */
128
void velociraptor_init(struct engine *e, const int linked_with_snap) {
129

130
#ifdef HAVE_VELOCIRAPTOR
131

Matthieu Schaller's avatar
Matthieu Schaller committed
132
133
134
135
  struct space *s = e->s;
  struct cosmoinfo cosmo_info;
  struct unitinfo unit_info;
  struct siminfo sim_info;
136
137
  const ticks tic = getticks();

Matthieu Schaller's avatar
Matthieu Schaller committed
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
  /* Set cosmological constants. */
  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;
  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;

  message("Scale factor: %e", cosmo_info.atime);
  message("Little h: %e", cosmo_info.littleh);
  message("Omega_m: %e", cosmo_info.Omega_m);
  message("Omega_b: %e", cosmo_info.Omega_b);
  message("Omega_Lambda: %e", cosmo_info.Omega_Lambda);
  message("Omega_cdm: %e", cosmo_info.Omega_cdm);
  message("w_de: %e", cosmo_info.w_de);

James Willis's avatar
James Willis committed
155
156
157
  if (e->cosmology->w != -1.)
    error("w_de is not 1. It is: %lf", e->cosmology->w);

Matthieu Schaller's avatar
Matthieu Schaller committed
158
159
160
161
162
163
  /* Set unit conversions. */
  unit_info.lengthtokpc = 1.0;
  unit_info.velocitytokms = 1.0;
  unit_info.masstosolarmass = 1.0;
  unit_info.energyperunitmass = 1.0;
  unit_info.gravity = e->physical_constants->const_newton_G;
James Willis's avatar
James Willis committed
164
  unit_info.hubbleunit = e->cosmology->H0 / e->cosmology->h;
Matthieu Schaller's avatar
Matthieu Schaller committed
165
166
167
168
169
170
171
172
173
174
175
176

  message("Length conversion factor: %e", unit_info.lengthtokpc);
  message("Velocity conversion factor: %e", unit_info.velocitytokms);
  message("Mass conversion factor: %e", unit_info.masstosolarmass);
  message("Potential conversion factor: %e", unit_info.energyperunitmass);
  message("G: %e", unit_info.gravity);
  message("H: %e", unit_info.hubbleunit);

  /* TODO: Find the total number of DM particles when running with star
   * particles and BHs. */
  const int total_nr_dmparts = e->total_nr_gparts - e->total_nr_parts;

177
  /* Set periodicity information. */
Matthieu Schaller's avatar
Matthieu Schaller committed
178
  if (e->s->periodic) {
179
180
181
    /* Physical size of box in VELOCIraptor units (kpc). */
    sim_info.period = unit_info.lengthtokpc * s->dim[0];
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
182
    sim_info.period = 0.0;
183
184
185
186
  }

  /* Are we running with cosmology? */
  if (e->policy & engine_policy_cosmology) {
Matthieu Schaller's avatar
Matthieu Schaller committed
187
    sim_info.icosmologicalsim = 1;
188
  } else {
Matthieu Schaller's avatar
Matthieu Schaller committed
189
    sim_info.icosmologicalsim = 0;
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
  }

  /* Append base name with the current output number */
  char outputFileName[PARSER_MAX_LINE_SIZE + 128];

  /* What should the filename be? */
  if (linked_with_snap) {
    snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128,
             "stf_%s_%04i.VELOCIraptor", e->snapshot_base_name,
             e->snapshot_output_count);
  } else {
    snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128, "%s_%04i.VELOCIraptor",
             e->stf_base_name, e->stf_output_count);
  }

  /* Gather the rest of the information */
  sim_info.zoomhigresolutionmass = -1.0; /* Placeholder. */
  sim_info.interparticlespacing = sim_info.period / cbrt(total_nr_dmparts);
Matthieu Schaller's avatar
Matthieu Schaller committed
208
209
210
211
212
213
214
215
  sim_info.spacedimension[0] = unit_info.lengthtokpc * s->dim[0];
  sim_info.spacedimension[1] = unit_info.lengthtokpc * s->dim[1];
  sim_info.spacedimension[2] = unit_info.lengthtokpc * s->dim[2];
  sim_info.numcells = s->nr_cells;

  sim_info.cellwidth[0] = unit_info.lengthtokpc * s->cells_top[0].width[0];
  sim_info.cellwidth[1] = unit_info.lengthtokpc * s->cells_top[0].width[1];
  sim_info.cellwidth[2] = unit_info.lengthtokpc * s->cells_top[0].width[2];
216

Matthieu Schaller's avatar
Matthieu Schaller committed
217
218
219
220
  sim_info.icellwidth[0] = s->iwidth[0] / unit_info.lengthtokpc;
  sim_info.icellwidth[1] = s->iwidth[1] / unit_info.lengthtokpc;
  sim_info.icellwidth[2] = s->iwidth[2] / unit_info.lengthtokpc;

James Willis's avatar
James Willis committed
221
222
  /* Only allocate cell location array on first call to velociraptor_init(). */
  if (e->cell_loc == NULL) {
223

James Willis's avatar
James Willis committed
224
225
226
227
228
229
230
231
232
233
    /* Allocate and populate top-level cell locations. */
    if (posix_memalign((void **)&(e->cell_loc), 32,
                       s->nr_cells * sizeof(struct cell_loc)) != 0)
      error("Failed to allocate top-level cell locations for VELOCIraptor.");

    for (int i = 0; i < s->nr_cells; i++) {
      e->cell_loc[i].loc[0] = unit_info.lengthtokpc * s->cells_top[i].loc[0];
      e->cell_loc[i].loc[1] = unit_info.lengthtokpc * s->cells_top[i].loc[1];
      e->cell_loc[i].loc[2] = unit_info.lengthtokpc * s->cells_top[i].loc[2];
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
234
235
236
237
  }

  sim_info.cell_loc = e->cell_loc;

238
  message("Config file name: %s", e->stf_config_file_name);
Matthieu Schaller's avatar
Matthieu Schaller committed
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
  message("Period: %e", sim_info.period);
  message("Zoom high res mass: %e", sim_info.zoomhigresolutionmass);
  message("Inter-particle spacing: %e", sim_info.interparticlespacing);
  message("Cosmological: %d", sim_info.icosmologicalsim);
  message("Space dimensions: (%e,%e,%e)", sim_info.spacedimension[0],
          sim_info.spacedimension[1], sim_info.spacedimension[2]);
  message("No. of top-level cells: %d", sim_info.numcells);
  message("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]);

  /* Initialise VELOCIraptor. */
254
255
  if (!InitVelociraptor(e->stf_config_file_name, outputFileName, cosmo_info,
                        unit_info, sim_info)) {
Matthieu Schaller's avatar
Matthieu Schaller committed
256
    error("Exiting. VELOCIraptor initialisation failed.");
257
  }
258
259
260
261

  if (e->verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
262
263
#else
  error("SWIFT not configure to run with VELOCIraptor.");
264
#endif /* HAVE_VELOCIRAPTOR */
265
266
267
268
269
270
}

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

275
#ifdef HAVE_VELOCIRAPTOR
276
277

  const struct space *s = e->s;
Matthieu Schaller's avatar
Matthieu Schaller committed
278
279
280
281
282
  struct gpart *gparts = s->gparts;
  struct part *parts = s->parts;
  const size_t nr_gparts = s->nr_gparts;
  const size_t nr_hydro_parts = s->nr_parts;
  const int nr_cells = s->nr_cells;
283
284

  const ticks tic = getticks();
285

Matthieu Schaller's avatar
Matthieu Schaller committed
286
  /* Allow thread to run on any core for the duration of the call to
287
288
   * VELOCIraptor so that  when OpenMP threads are spawned
   * they can run on any core on the processor. */
Matthieu Schaller's avatar
Matthieu Schaller committed
289
290
291
292
  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. */
293
  cpu_set_t cpuset;
Matthieu Schaller's avatar
Matthieu Schaller committed
294
295
296
297
298
299
  CPU_ZERO(&cpuset);
  for (int j = 0; j < nr_cores; j++) CPU_SET(j, &cpuset);

  pthread_setaffinity_np(thread, sizeof(cpu_set_t), &cpuset);

  /* Allocate and populate array of cell node IDs. */
300
  int *cell_node_ids = NULL;
Matthieu Schaller's avatar
Matthieu Schaller committed
301
302
303
304
305
  if (posix_memalign((void **)&cell_node_ids, 32, nr_cells * sizeof(int)) != 0)
    error("Failed to allocate list of cells node IDs for VELOCIraptor.");

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

306
  /* Mention the number of particles being sent */
307
308
309
  if (e->verbose)
    message("MPI rank %d sending %zu gparts to VELOCIraptor.", engine_rank,
            nr_gparts);
Matthieu Schaller's avatar
Matthieu Schaller committed
310

311
  /* Append base name with the current output number */
Matthieu Schaller's avatar
Matthieu Schaller committed
312
  char outputFileName[PARSER_MAX_LINE_SIZE + 128];
313

314
315
316
317
318
319
320
321
322
  /* What should the filename be? */
  if (linked_with_snap) {
    snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128,
             "stf_%s_%04i.VELOCIraptor", e->snapshot_base_name,
             e->snapshot_output_count);
  } else {
    snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128, "%s_%04i.VELOCIraptor",
             e->stf_base_name, e->stf_output_count);
  }
Matthieu Schaller's avatar
Matthieu Schaller committed
323
324
325
326
327
328
329
330
331
332

  /* Allocate and populate an array of swift_vel_parts to be passed to
   * VELOCIraptor. */
  struct swift_vel_part *swift_parts = NULL;
  if (posix_memalign((void **)&swift_parts, part_align,
                     nr_gparts * sizeof(struct swift_vel_part)) != 0)
    error("Failed to allocate array of particles for VELOCIraptor.");

  bzero(swift_parts, nr_gparts * sizeof(struct swift_vel_part));

333
  const float a_inv = e->cosmology->a_inv;
Matthieu Schaller's avatar
Matthieu Schaller committed
334
335
336

  /* Convert particle properties into VELOCIraptor units */
  for (size_t i = 0; i < nr_gparts; i++) {
337

Matthieu Schaller's avatar
Matthieu Schaller committed
338
339
340
    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];
341
342
343
344
345

    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;

Matthieu Schaller's avatar
Matthieu Schaller committed
346
    swift_parts[i].mass = gravity_get_mass(&gparts[i]);
James Willis's avatar
James Willis committed
347
    swift_parts[i].potential = gravity_get_comoving_potential(&gparts[i]);
348

Matthieu Schaller's avatar
Matthieu Schaller committed
349
350
351
352
    swift_parts[i].type = gparts[i].type;

    /* Set gas particle IDs from their hydro counterparts and set internal
     * energies. */
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
    switch (gparts[i].type) {

      case swift_type_gas:

        swift_parts[i].id = parts[-gparts[i].id_or_neg_offset].id;
        swift_parts[i].u = hydro_get_drifted_physical_internal_energy(
            &parts[-gparts[i].id_or_neg_offset], e->cosmology);
        break;

      case swift_type_dark_matter:

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

      default:
        error("Particle type not handled by VELOCIraptor.");
370
    }
Matthieu Schaller's avatar
Matthieu Schaller committed
371
  }
372

Matthieu Schaller's avatar
Matthieu Schaller committed
373
  /* Call VELOCIraptor. */
374
  if (!InvokeVelociraptor(nr_gparts, nr_hydro_parts, e->stf_output_count,
375
                          swift_parts, cell_node_ids, outputFileName))
Matthieu Schaller's avatar
Matthieu Schaller committed
376
377
378
379
380
381
382
383
384
    error("Exiting. Call to VELOCIraptor failed on rank: %d.", e->nodeID);

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

  /* Free cell node ids after VELOCIraptor has copied them. */
  free(cell_node_ids);
  free(swift_parts);

385
  /* Increase output counter (if not linked with snapshots) */
386
  if (!linked_with_snap) e->stf_output_count++;
387

388
389
390
  if (e->verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
391
#else
Matthieu Schaller's avatar
Matthieu Schaller committed
392
  error("SWIFT not configure to run with VELOCIraptor.");
393
#endif /* HAVE_VELOCIRAPTOR */
394
}