testActivePair.c 23.3 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
 *
 * 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/>.
 *
 ******************************************************************************/
#include "../config.h"

/* Some standard headers. */
#include <fenv.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
26
#include <time.h>
James Willis's avatar
James Willis committed
27
#include <unistd.h>
28 29 30 31

/* Local headers. */
#include "swift.h"

32 33
#define NODE_ID 1

34 35
/* Typdef function pointer for interaction function. */
typedef void (*interaction_func)(struct runner *, struct cell *, struct cell *);
Josh Borrow's avatar
Josh Borrow committed
36 37
typedef void (*init_func)(struct cell *, const struct cosmology *,
                          const struct hydro_props *);
38
typedef void (*finalise_func)(struct cell *, const struct cosmology *);
39

40 41 42 43 44 45 46 47 48 49 50 51 52 53
/**
 * @brief Constructs a cell and all of its particle in a valid state prior to
 * a DOPAIR or DOSELF calcuation.
 *
 * @param n The cube root of the number of particles.
 * @param offset The position of the cell offset from (0,0,0).
 * @param size The cell size.
 * @param h The smoothing length of the particles in units of the inter-particle
 * separation.
 * @param density The density of the fluid.
 * @param partId The running counter of IDs.
 * @param pert The perturbation to apply to the particles in the cell in units
 * of the inter-particle separation.
 * @param h_pert The perturbation to apply to the smoothing length.
54 55
 * @param fraction_active The fraction of particles that should be active in the
 * cell.
56 57
 */
struct cell *make_cell(size_t n, double *offset, double size, double h,
58
                       double density, long long *partId, double pert,
59
                       double h_pert, double fraction_active) {
60 61
  const size_t count = n * n * n;
  const double volume = size * size * size;
62
  float h_max = 0.f;
63
  struct cell *cell = (struct cell *)malloc(sizeof(struct cell));
64 65
  bzero(cell, sizeof(struct cell));

66
  if (posix_memalign((void **)&cell->hydro.parts, part_align,
67 68 69
                     count * sizeof(struct part)) != 0) {
    error("couldn't allocate particles, no. of particles: %d", (int)count);
  }
70
  bzero(cell->hydro.parts, count * sizeof(struct part));
71 72

  /* Construct the parts */
73
  struct part *part = cell->hydro.parts;
74 75 76 77
  for (size_t x = 0; x < n; ++x) {
    for (size_t y = 0; y < n; ++y) {
      for (size_t z = 0; z < n; ++z) {
        part->x[0] =
James Willis's avatar
James Willis committed
78 79
            offset[0] +
            size * (x + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
80
        part->x[1] =
James Willis's avatar
James Willis committed
81 82
            offset[1] +
            size * (y + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
83
        part->x[2] =
James Willis's avatar
James Willis committed
84 85
            offset[2] +
            size * (z + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
86 87 88
        part->v[0] = random_uniform(-0.05, 0.05);
        part->v[1] = random_uniform(-0.05, 0.05);
        part->v[2] = random_uniform(-0.05, 0.05);
89 90 91 92 93 94

        if (h_pert)
          part->h = size * h * random_uniform(1.f, h_pert) / (float)n;
        else
          part->h = size * h / (float)n;
        h_max = fmaxf(h_max, part->h);
95
        part->id = ++(*partId);
96

97
/* Set the mass */
98
#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
99
        part->conserved.mass = density * volume / count;
100 101 102 103 104

#ifdef SHADOWFAX_SPH
        double anchor[3] = {0., 0., 0.};
        double side[3] = {1., 1., 1.};
        voronoi_cell_init(&part->cell, part->x, anchor, side);
105
#endif /* SHADOWFAX_SPH */
106

107 108
#else
        part->mass = density * volume / count;
109
#endif /* GIZMO_MFV_SPH */
110

111 112 113
/* Set the thermodynamic variable */
#if defined(GADGET2_SPH)
        part->entropy = 1.f;
114 115
#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) ||           \
    defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH) || \
Josh Borrow's avatar
Josh Borrow committed
116
    defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
117 118
        part->u = 1.f;
#elif defined(HOPKINS_PE_SPH)
119 120 121
        part->entropy = 1.f;
        part->entropy_one_over_gamma = 1.f;
#endif
122 123

        /* Set the time-bin */
James Willis's avatar
James Willis committed
124 125
        if (random_uniform(0, 1.f) < fraction_active)
          part->time_bin = 1;
126 127
        else
          part->time_bin = num_time_bins + 1;
128 129 130 131 132 133 134 135 136 137 138 139 140

#ifdef SWIFT_DEBUG_CHECKS
        part->ti_drift = 8;
        part->ti_kick = 8;
#endif

        ++part;
      }
    }
  }

  /* Cell properties */
  cell->split = 0;
141 142
  cell->hydro.h_max = h_max;
  cell->hydro.count = count;
143
  cell->hydro.dx_max_part = 0.;
144
  cell->hydro.dx_max_sort = 0.;
145 146 147
  cell->width[0] = size;
  cell->width[1] = size;
  cell->width[2] = size;
148 149 150 151
  cell->loc[0] = offset[0];
  cell->loc[1] = offset[1];
  cell->loc[2] = offset[2];

152
  cell->hydro.ti_old_part = 8;
153 154
  cell->hydro.ti_end_min = 8;
  cell->hydro.ti_end_max = 10;
Matthieu Schaller's avatar
Matthieu Schaller committed
155
  cell->nodeID = NODE_ID;
156

157
  shuffle_particles(cell->hydro.parts, cell->hydro.count);
158

159
  cell->hydro.sorted = 0;
160
  cell->hydro.sort = NULL;
James Willis's avatar
James Willis committed
161

162 163 164 165
  return cell;
}

void clean_up(struct cell *ci) {
166
  cell_free_hydro_sorts(ci);
167 168 169 170 171 172
  free(ci);
}

/**
 * @brief Initializes all particles field to be ready for a density calculation
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
173
void zero_particle_fields_density(struct cell *c, const struct cosmology *cosmo,
Josh Borrow's avatar
Josh Borrow committed
174
                                  const struct hydro_props *hydro_props) {
175 176
  for (int pid = 0; pid < c->hydro.count; pid++) {
    hydro_init_part(&c->hydro.parts[pid], NULL);
177 178 179
  }
}

180
/**
181
 * @brief Initializes all particles field to be ready for a force calculation
182
 */
Josh Borrow's avatar
Josh Borrow committed
183 184
void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo,
                                const struct hydro_props *hydro_props) {
185 186 187
  for (int pid = 0; pid < c->hydro.count; pid++) {
    struct part *p = &c->hydro.parts[pid];
    struct xpart *xp = &c->hydro.xparts[pid];
188 189 190 191 192 193 194 195 196 197 198

/* Mimic the result of a density calculation */
#ifdef GADGET2_SPH
    p->rho = 1.f;
    p->density.rho_dh = 0.f;
    p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
    p->density.wcount_dh = 0.f;
    p->density.rot_v[0] = 0.f;
    p->density.rot_v[1] = 0.f;
    p->density.rot_v[2] = 0.f;
    p->density.div_v = 0.f;
199
#endif /* GADGET-2 */
Josh Borrow's avatar
Josh Borrow committed
200
#if defined(MINIMAL_SPH) || defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
201 202 203 204
    p->rho = 1.f;
    p->density.rho_dh = 0.f;
    p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
    p->density.wcount_dh = 0.f;
Josh Borrow's avatar
Josh Borrow committed
205
    p->viscosity.v_sig = hydro_get_comoving_soundspeed(p);
206
#endif /* MINIMAL */
207 208 209 210 211 212 213
#ifdef HOPKINS_PE_SPH
    p->rho = 1.f;
    p->rho_bar = 1.f;
    p->density.rho_dh = 0.f;
    p->density.pressure_dh = 0.f;
    p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
    p->density.wcount_dh = 0.f;
214
#endif /* PRESSURE-ENTROPY */
Josh Borrow's avatar
Josh Borrow committed
215 216
#if defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN) || \
    defined(ANARCHY_PU_SPH)
217 218 219 220 221 222 223
    p->rho = 1.f;
    p->pressure_bar = 0.6666666;
    p->density.rho_dh = 0.f;
    p->density.pressure_bar_dh = 0.f;
    p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
    p->density.wcount_dh = 0.f;
#endif /* PRESSURE-ENERGY */
Josh Borrow's avatar
Josh Borrow committed
224
#if defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH)
Josh Borrow's avatar
Josh Borrow committed
225
    /* Initialise viscosity variables */
226
    p->force.pressure = hydro_get_comoving_pressure(p);
Josh Borrow's avatar
Josh Borrow committed
227 228 229 230 231
    p->viscosity.alpha = 0.8;
    p->viscosity.div_v = 0.f;
    p->viscosity.div_v_previous_step = 0.f;
    p->viscosity.v_sig = hydro_get_comoving_soundspeed(p);
#endif /* ANARCHY_PU_SPH viscosity variables */
232 233

    /* And prepare for a round of force tasks. */
Josh Borrow's avatar
Josh Borrow committed
234
    hydro_prepare_force(p, xp, cosmo, hydro_props, 0.);
235 236 237 238 239 240 241 242
    hydro_reset_acceleration(p);
  }
}

/**
 * @brief Ends the density loop by adding the appropriate coefficients
 */
void end_calculation_density(struct cell *c, const struct cosmology *cosmo) {
243 244
  for (int pid = 0; pid < c->hydro.count; pid++) {
    hydro_end_density(&c->hydro.parts[pid], cosmo);
245 246

    /* Recover the common "Neighbour number" definition */
247 248
    c->hydro.parts[pid].density.wcount *= pow_dimension(c->hydro.parts[pid].h);
    c->hydro.parts[pid].density.wcount *= kernel_norm;
249 250 251
  }
}

252 253 254 255
/**
 * @brief Ends the force loop by adding the appropriate coefficients
 */
void end_calculation_force(struct cell *c, const struct cosmology *cosmo) {
256
  for (int pid = 0; pid < c->hydro.count; pid++) {
257 258
    struct part *volatile part = &c->hydro.parts[pid];
    hydro_end_force(part, cosmo);
259 260 261
  }
}

262 263 264 265
/**
 * @brief Dump all the particles to a file
 */
void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
266
  FILE *file = fopen(fileName, "a");
267 268

  /* Write header */
269
  fprintf(file, "# %4s %13s %13s\n", "ID", "wcount", "h_dt");
270 271 272

  fprintf(file, "# ci --------------------------------------------\n");

273 274 275 276
  for (int pid = 0; pid < ci->hydro.count; pid++) {
    fprintf(file, "%6llu %13e %13e\n", ci->hydro.parts[pid].id,
            ci->hydro.parts[pid].density.wcount,
            ci->hydro.parts[pid].force.h_dt);
277 278 279 280
  }

  fprintf(file, "# cj --------------------------------------------\n");

281 282 283 284
  for (int pjd = 0; pjd < cj->hydro.count; pjd++) {
    fprintf(file, "%6llu %13e %13e\n", cj->hydro.parts[pjd].id,
            cj->hydro.parts[pjd].density.wcount,
            cj->hydro.parts[pjd].force.h_dt);
285 286 287 288 289 290 291
  }

  fclose(file);
}

/* Just a forward declaration... */
void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
James Willis's avatar
James Willis committed
292 293
void runner_dopair2_force_vec(struct runner *r, struct cell *ci,
                              struct cell *cj);
294 295 296
void runner_doself1_density_vec(struct runner *r, struct cell *ci);
void runner_dopair1_branch_density(struct runner *r, struct cell *ci,
                                   struct cell *cj);
297
void runner_dopair2_branch_force(struct runner *r, struct cell *ci,
James Willis's avatar
James Willis committed
298
                                 struct cell *cj);
299

300
/**
James Willis's avatar
James Willis committed
301 302
 * @brief Computes the pair interactions of two cells using SWIFT and a brute
 * force implementation.
303
 */
James Willis's avatar
James Willis committed
304 305
void test_pair_interactions(struct runner *runner, struct cell **ci,
                            struct cell **cj, char *swiftOutputFileName,
James Willis's avatar
James Willis committed
306 307
                            char *bruteForceOutputFileName,
                            interaction_func serial_interaction,
308 309
                            interaction_func vec_interaction, init_func init,
                            finalise_func finalise) {
310

Loic Hausammann's avatar
Loic Hausammann committed
311 312
  runner_do_hydro_sort(runner, *ci, 0x1FFF, 0, 0);
  runner_do_hydro_sort(runner, *cj, 0x1FFF, 0, 0);
313 314

  /* Zero the fields */
Josh Borrow's avatar
Josh Borrow committed
315 316
  init(*ci, runner->e->cosmology, runner->e->hydro_properties);
  init(*cj, runner->e->cosmology, runner->e->hydro_properties);
317 318

  /* Run the test */
319
  vec_interaction(runner, *ci, *cj);
320 321

  /* Let's get physical ! */
322 323
  finalise(*ci, runner->e->cosmology);
  finalise(*cj, runner->e->cosmology);
James Willis's avatar
James Willis committed
324

325 326 327 328 329 330
  /* Dump if necessary */
  dump_particle_fields(swiftOutputFileName, *ci, *cj);

  /* Now perform a brute-force version for accuracy tests */

  /* Zero the fields */
Josh Borrow's avatar
Josh Borrow committed
331 332
  init(*ci, runner->e->cosmology, runner->e->hydro_properties);
  init(*cj, runner->e->cosmology, runner->e->hydro_properties);
333 334

  /* Run the brute-force test */
335
  serial_interaction(runner, *ci, *cj);
336 337

  /* Let's get physical ! */
338 339
  finalise(*ci, runner->e->cosmology);
  finalise(*cj, runner->e->cosmology);
340 341 342 343

  dump_particle_fields(bruteForceOutputFileName, *ci, *cj);
}

344 345 346
/**
 * @brief Computes the pair interactions of two cells in various configurations.
 */
James Willis's avatar
James Willis committed
347 348 349 350
void test_all_pair_interactions(
    struct runner *runner, double *offset2, size_t particles, double size,
    double h, double rho, long long *partId, double perturbation, double h_pert,
    char *swiftOutputFileName, char *bruteForceOutputFileName,
351 352
    interaction_func serial_interaction, interaction_func vec_interaction,
    init_func init, finalise_func finalise) {
353 354 355 356

  double offset1[3] = {0, 0, 0};
  struct cell *ci, *cj;

357
  /* Only one particle in each cell. */
James Willis's avatar
James Willis committed
358 359
  ci = make_cell(1, offset1, size, h, rho, partId, perturbation, h_pert, 1.);
  cj = make_cell(1, offset2, size, h, rho, partId, perturbation, h_pert, 1.);
360 361

  test_pair_interactions(runner, &ci, &cj, swiftOutputFileName,
James Willis's avatar
James Willis committed
362
                         bruteForceOutputFileName, serial_interaction,
363
                         vec_interaction, init, finalise);
364 365 366 367

  clean_up(ci);
  clean_up(cj);

368
  /* All active particles. */
James Willis's avatar
James Willis committed
369 370 371 372
  ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert,
                 1.);
  cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert,
                 1.);
373

James Willis's avatar
James Willis committed
374
  test_pair_interactions(runner, &ci, &cj, swiftOutputFileName,
James Willis's avatar
James Willis committed
375
                         bruteForceOutputFileName, serial_interaction,
376
                         vec_interaction, init, finalise);
377 378 379 380 381

  clean_up(ci);
  clean_up(cj);

  /* Half particles are active. */
James Willis's avatar
James Willis committed
382 383 384 385 386 387
  ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert,
                 0.5);
  cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert,
                 0.5);

  test_pair_interactions(runner, &ci, &cj, swiftOutputFileName,
James Willis's avatar
James Willis committed
388
                         bruteForceOutputFileName, serial_interaction,
389
                         vec_interaction, init, finalise);
390 391 392 393 394

  clean_up(ci);
  clean_up(cj);

  /* All particles inactive. */
James Willis's avatar
James Willis committed
395 396 397 398
  ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert,
                 0.);
  cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert,
                 0.);
399

James Willis's avatar
James Willis committed
400
  test_pair_interactions(runner, &ci, &cj, swiftOutputFileName,
James Willis's avatar
James Willis committed
401
                         bruteForceOutputFileName, serial_interaction,
402
                         vec_interaction, init, finalise);
403 404 405 406 407

  clean_up(ci);
  clean_up(cj);

  /* 10% of particles active. */
James Willis's avatar
James Willis committed
408 409 410 411 412 413
  ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert,
                 0.1);
  cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert,
                 0.1);

  test_pair_interactions(runner, &ci, &cj, swiftOutputFileName,
James Willis's avatar
James Willis committed
414
                         bruteForceOutputFileName, serial_interaction,
415
                         vec_interaction, init, finalise);
416 417 418 419 420

  clean_up(ci);
  clean_up(cj);

  /* One active cell one inactive cell. */
James Willis's avatar
James Willis committed
421 422 423 424 425 426
  ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert,
                 1.0);
  cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert,
                 0.);

  test_pair_interactions(runner, &ci, &cj, swiftOutputFileName,
James Willis's avatar
James Willis committed
427
                         bruteForceOutputFileName, serial_interaction,
428
                         vec_interaction, init, finalise);
429 430 431 432 433

  clean_up(ci);
  clean_up(cj);

  /* One active cell one inactive cell. */
James Willis's avatar
James Willis committed
434 435 436 437 438 439
  ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert,
                 0.);
  cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert,
                 1.0);

  test_pair_interactions(runner, &ci, &cj, swiftOutputFileName,
James Willis's avatar
James Willis committed
440
                         bruteForceOutputFileName, serial_interaction,
441
                         vec_interaction, init, finalise);
442 443 444 445 446 447 448 449

  clean_up(ci);
  clean_up(cj);

  /* Smaller cells, all active. */
  ci = make_cell(2, offset1, size, h, rho, partId, perturbation, h_pert, 1.0);
  cj = make_cell(2, offset2, size, h, rho, partId, perturbation, h_pert, 1.0);

James Willis's avatar
James Willis committed
450
  test_pair_interactions(runner, &ci, &cj, swiftOutputFileName,
James Willis's avatar
James Willis committed
451
                         bruteForceOutputFileName, serial_interaction,
452
                         vec_interaction, init, finalise);
James Willis's avatar
James Willis committed
453

454 455 456 457 458 459 460
  clean_up(ci);
  clean_up(cj);

  /* Different numbers of particles in each cell. */
  ci = make_cell(10, offset1, size, h, rho, partId, perturbation, h_pert, 0.5);
  cj = make_cell(3, offset2, size, h, rho, partId, perturbation, h_pert, 0.75);

James Willis's avatar
James Willis committed
461
  test_pair_interactions(runner, &ci, &cj, swiftOutputFileName,
James Willis's avatar
James Willis committed
462
                         bruteForceOutputFileName, serial_interaction,
463
                         vec_interaction, init, finalise);
James Willis's avatar
James Willis committed
464

465 466 467 468
  clean_up(ci);
  clean_up(cj);

  /* One cell inactive and the other only half active. */
James Willis's avatar
James Willis committed
469 470 471 472 473 474
  ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert,
                 0.5);
  cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert,
                 0.);

  test_pair_interactions(runner, &ci, &cj, swiftOutputFileName,
James Willis's avatar
James Willis committed
475
                         bruteForceOutputFileName, serial_interaction,
476
                         vec_interaction, init, finalise);
477 478 479

  clean_up(ci);
  clean_up(cj);
James Willis's avatar
James Willis committed
480

481
  /* One cell inactive and the other only half active. */
James Willis's avatar
James Willis committed
482 483 484 485 486 487
  ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert,
                 0.);
  cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert,
                 0.5);

  test_pair_interactions(runner, &ci, &cj, swiftOutputFileName,
James Willis's avatar
James Willis committed
488
                         bruteForceOutputFileName, serial_interaction,
489
                         vec_interaction, init, finalise);
490 491 492 493 494 495

  /* Clean things to make the sanitizer happy ... */
  clean_up(ci);
  clean_up(cj);
}

496
int main(int argc, char *argv[]) {
497
  size_t particles = 0, runs = 0, type = 0;
498
  double h = 1.23485, size = 1., rho = 1.;
499
  double perturbation = 0.1, h_pert = 1.1;
500 501
  struct space space;
  struct engine engine;
502
  struct cosmology cosmo;
Josh Borrow's avatar
Josh Borrow committed
503
  struct hydro_props hydro_props;
504
  struct runner *runner;
505
  char c;
506
  static long long partId = 0;
507
  char outputFileNameExtension[100] = "";
508 509
  char swiftOutputFileName[200] = "";
  char bruteForceOutputFileName[200] = "";
510 511 512 513 514

  /* Initialize CPU frequency, this also starts time. */
  unsigned long long cpufreq = 0;
  clocks_set_cpufreq(cpufreq);

515 516
/* Choke on FP-exceptions */
#ifdef HAVE_FE_ENABLE_EXCEPT
517
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
518
#endif
519

520 521
  /* Generate a RNG seed from time. */
  unsigned int seed = time(NULL);
522

523
  while ((c = getopt(argc, argv, "h:p:n:r:t:d:s:f:")) != -1) {
524 525 526 527 528
    switch (c) {
      case 'h':
        sscanf(optarg, "%lf", &h);
        break;
      case 'p':
529 530 531
        sscanf(optarg, "%lf", &h_pert);
        break;
      case 'n':
532 533 534 535 536 537 538 539 540 541 542
        sscanf(optarg, "%zu", &particles);
        break;
      case 'r':
        sscanf(optarg, "%zu", &runs);
        break;
      case 't':
        sscanf(optarg, "%zu", &type);
        break;
      case 'd':
        sscanf(optarg, "%lf", &perturbation);
        break;
543 544 545
      case 's':
        sscanf(optarg, "%u", &seed);
        break;
546 547 548 549 550 551 552 553 554 555 556
      case 'f':
        strcpy(outputFileNameExtension, optarg);
        break;
      case '?':
        error("Unknown option.");
        break;
    }
  }

  if (h < 0 || particles == 0 || runs == 0 || type > 2) {
    printf(
557
        "\nUsage: %s -n PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
558 559 560 561
        "\nGenerates a cell pair, filled with particles on a Cartesian grid."
        "\nThese are then interacted using runner_dopair1_density."
        "\n\nOptions:"
        "\n-t TYPE=0          - cells share face (0), edge (1) or corner (2)"
562 563
        "\n-h DISTANCE=1.2348 - smoothing length"
        "\n-p                 - Random fractional change in h, h=h*random(1,p)"
564
        "\n-d pert            - perturbation to apply to the particles [0,1["
565
        "\n-s seed            - seed for RNG"
566 567 568 569 570
        "\n-f fileName        - part of the file name used to save the dumps\n",
        argv[0]);
    exit(1);
  }

571 572 573
  /* Seed RNG. */
  message("Seed used for RNG: %d", seed);
  srand(seed);
James Willis's avatar
James Willis committed
574

575
  space.periodic = 0;
James Willis's avatar
James Willis committed
576 577 578
  space.dim[0] = 3.;
  space.dim[1] = 3.;
  space.dim[2] = 3.;
579 580 581 582 583

  engine.s = &space;
  engine.time = 0.1f;
  engine.ti_current = 8;
  engine.max_active_bin = num_time_bins;
584
  engine.nodeID = NODE_ID;
James Willis's avatar
James Willis committed
585

586 587
  cosmology_init_no_cosmo(&cosmo);
  engine.cosmology = &cosmo;
Josh Borrow's avatar
Josh Borrow committed
588 589
  hydro_props_init_no_hydro(&hydro_props);
  engine.hydro_properties = &hydro_props;
590

591 592
  if (posix_memalign((void **)&runner, SWIFT_STRUCT_ALIGNMENT,
                     sizeof(struct runner)) != 0) {
593 594
    error("couldn't allocate runner");
  }
James Willis's avatar
James Willis committed
595

596
  runner->e = &engine;
597

598
  /* Create output file names. */
lhausamm's avatar
lhausamm committed
599 600
  sprintf(swiftOutputFileName, "swift_dopair_%.150s.dat",
          outputFileNameExtension);
lhausamm's avatar
lhausamm committed
601
  sprintf(bruteForceOutputFileName, "brute_force_pair_%.150s.dat",
602
          outputFileNameExtension);
James Willis's avatar
James Willis committed
603

604 605 606 607 608 609 610 611 612 613 614
  /* Delete files if they already exist. */
  remove(swiftOutputFileName);
  remove(bruteForceOutputFileName);

#ifdef WITH_VECTORIZATION
  runner->ci_cache.count = 0;
  cache_init(&runner->ci_cache, 512);
  runner->cj_cache.count = 0;
  cache_init(&runner->cj_cache, 512);
#endif

James Willis's avatar
James Willis committed
615
  double offset[3] = {1., 0., 0.};
616

617 618 619
  /* Define which interactions to call */
  interaction_func serial_inter_func = &pairs_all_density;
  interaction_func vec_inter_func = &runner_dopair1_branch_density;
620 621
  init_func init = &zero_particle_fields_density;
  finalise_func finalise = &end_calculation_density;
622

623
  /* Test a pair of cells face-on. */
James Willis's avatar
James Willis committed
624 625
  test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId,
                             perturbation, h_pert, swiftOutputFileName,
James Willis's avatar
James Willis committed
626
                             bruteForceOutputFileName, serial_inter_func,
627
                             vec_inter_func, init, finalise);
628 629 630 631 632 633 634

  /* Test a pair of cells edge-on. */
  offset[0] = 1.;
  offset[1] = 1.;
  offset[2] = 0.;
  test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId,
                             perturbation, h_pert, swiftOutputFileName,
James Willis's avatar
James Willis committed
635
                             bruteForceOutputFileName, serial_inter_func,
636
                             vec_inter_func, init, finalise);
637 638 639 640 641 642 643

  /* Test a pair of cells corner-on. */
  offset[0] = 1.;
  offset[1] = 1.;
  offset[2] = 1.;
  test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId,
                             perturbation, h_pert, swiftOutputFileName,
James Willis's avatar
James Willis committed
644
                             bruteForceOutputFileName, serial_inter_func,
645
                             vec_inter_func, init, finalise);
James Willis's avatar
James Willis committed
646

647 648
  /* Re-assign function pointers. */
  serial_inter_func = &pairs_all_force;
649
  vec_inter_func = &runner_dopair2_branch_force;
650 651
  init = &zero_particle_fields_force;
  finalise = &end_calculation_force;
652 653

  /* Create new output file names. */
lhausamm's avatar
lhausamm committed
654
  sprintf(swiftOutputFileName, "swift_dopair2_force_%.150s.dat",
James Willis's avatar
James Willis committed
655
          outputFileNameExtension);
lhausamm's avatar
lhausamm committed
656
  sprintf(bruteForceOutputFileName, "brute_force_dopair2_%.150s.dat",
657 658 659 660 661
          outputFileNameExtension);

  /* Delete files if they already exist. */
  remove(swiftOutputFileName);
  remove(bruteForceOutputFileName);
James Willis's avatar
James Willis committed
662

663
  /* Test a pair of cells face-on. */
664 665 666
  offset[0] = 1.;
  offset[1] = 0.;
  offset[2] = 0.;
James Willis's avatar
James Willis committed
667 668
  test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId,
                             perturbation, h_pert, swiftOutputFileName,
James Willis's avatar
James Willis committed
669
                             bruteForceOutputFileName, serial_inter_func,
670
                             vec_inter_func, init, finalise);
James Willis's avatar
James Willis committed
671

672
  /* Test a pair of cells edge-on. */
James Willis's avatar
James Willis committed
673 674 675 676 677
  offset[0] = 1.;
  offset[1] = 1.;
  offset[2] = 0.;
  test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId,
                             perturbation, h_pert, swiftOutputFileName,
James Willis's avatar
James Willis committed
678
                             bruteForceOutputFileName, serial_inter_func,
679
                             vec_inter_func, init, finalise);
680

681
  /* Test a pair of cells corner-on. */
James Willis's avatar
James Willis committed
682 683 684 685 686
  offset[0] = 1.;
  offset[1] = 1.;
  offset[2] = 1.;
  test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId,
                             perturbation, h_pert, swiftOutputFileName,
James Willis's avatar
James Willis committed
687
                             bruteForceOutputFileName, serial_inter_func,
688
                             vec_inter_func, init, finalise);
689 690
  return 0;
}