test27cells.c 20.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/*******************************************************************************
 * 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/>.
 *
 ******************************************************************************/

Matthieu Schaller's avatar
Matthieu Schaller committed
20 21 22 23
/* Config parameters. */
#include "../config.h"

/* Some standard headers. */
24
#include <fenv.h>
25
#include <stdio.h>
26 27 28
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
Matthieu Schaller's avatar
Matthieu Schaller committed
29 30

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

33
#if defined(WITH_VECTORIZATION)
34 35
#define DOSELF1 runner_doself1_branch_density
#define DOSELF1_SUBSET runner_doself_subset_branch_density
36
#define DOPAIR1_SUBSET runner_dopair_subset_branch_density
37
#define DOPAIR1 runner_dopair1_branch_density
38
#ifdef TEST_DOSELF_SUBSET
39
#define DOSELF1_NAME "runner_doself_subset_branch_density"
40
#else
41
#define DOSELF1_NAME "runner_doself1_branch_density"
42
#endif
43
#ifdef TEST_DOPAIR_SUBSET
44
#define DOPAIR1_NAME "runner_dopair_subset_branch_density"
45
#else
46
#define DOPAIR1_NAME "runner_dopair1_branch_density"
47
#endif
48 49 50
#endif

#ifndef DOSELF1
51
#define DOSELF1 runner_doself1_branch_density
52
#define DOSELF1_SUBSET runner_doself_subset_branch_density
53
#ifdef TEST_DOSELF_SUBSET
54
#define DOSELF1_NAME "runner_doself_subset_branch_density"
55
#else
56
#define DOSELF1_NAME "runner_doself1_branch_density"
57
#endif
58
#endif
59

60
#ifndef DOPAIR1
61
#define DOPAIR1 runner_dopair1_branch_density
62
#define DOPAIR1_SUBSET runner_dopair_subset_branch_density
63
#ifdef TEST_DOPAIR_SUBSET
64
#define DOPAIR1_NAME "runner_dopair1_subset_branch_density"
65
#else
66
#define DOPAIR1_NAME "runner_dopair1_branch_density"
67
#endif
68
#endif
69

70
#define NODE_ID 0
71

72 73 74 75 76 77 78
enum velocity_types {
  velocity_zero,
  velocity_random,
  velocity_divergent,
  velocity_rotating
};

Matthieu Schaller's avatar
Matthieu Schaller committed
79 80 81 82 83 84 85
/**
 * @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.
James Willis's avatar
James Willis committed
86
 * @param h The smoothing length of the particles in units of the inter-particle
87
 * separation.
Matthieu Schaller's avatar
Matthieu Schaller committed
88 89
 * @param density The density of the fluid.
 * @param partId The running counter of IDs.
James Willis's avatar
James Willis committed
90
 * @param pert The perturbation to apply to the particles in the cell in units
91
 * of the inter-particle separation.
Matthieu Schaller's avatar
Matthieu Schaller committed
92
 * @param vel The type of velocity field (0, random, divergent, rotating)
93
 * @param h_pert The perturbation to apply to the smoothing length.
94
 */
95
struct cell *make_cell(size_t n, double *offset, double size, double h,
96
                       double density, long long *partId, double pert,
97
                       enum velocity_types vel, double h_pert) {
98 99
  const size_t count = n * n * n;
  const double volume = size * size * size;
100
  float h_max = 0.f;
101
  struct cell *cell = (struct cell *)malloc(sizeof(struct cell));
102 103
  bzero(cell, sizeof(struct cell));

104
  if (posix_memalign((void **)&cell->hydro.parts, part_align,
105
                     count * sizeof(struct part)) != 0) {
106 107
    error("couldn't allocate particles, no. of particles: %d", (int)count);
  }
108
  bzero(cell->hydro.parts, count * sizeof(struct part));
109

110
  /* Construct the parts */
111
  struct part *part = cell->hydro.parts;
112 113 114 115 116 117 118 119 120 121 122 123
  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] =
            offset[0] +
            size * (x + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
        part->x[1] =
            offset[1] +
            size * (y + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
        part->x[2] =
            offset[2] +
            size * (z + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
        switch (vel) {
          case velocity_zero:
            part->v[0] = 0.f;
            part->v[1] = 0.f;
            part->v[2] = 0.f;
            break;
          case velocity_random:
            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);
            break;
          case velocity_divergent:
            part->v[0] = part->x[0] - 1.5 * size;
            part->v[1] = part->x[1] - 1.5 * size;
            part->v[2] = part->x[2] - 1.5 * size;
            break;
          case velocity_rotating:
            part->v[0] = part->x[1];
            part->v[1] = -part->x[0];
            part->v[2] = 0.f;
            break;
        }
146
        if (h_pert)
147
          part->h = size * h * random_uniform(1.f, h_pert) / (float)n;
148 149
        else
          part->h = size * h / (float)n;
150
        h_max = fmaxf(h_max, part->h);
151
        part->id = ++(*partId);
152

153
#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
154
        part->conserved.mass = density * volume / count;
155 156 157 158 159 160 161

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

162
#else
163
        part->mass = density * volume / count;
164
#endif
165 166 167 168 169 170

#if defined(HOPKINS_PE_SPH)
        part->entropy = 1.f;
        part->entropy_one_over_gamma = 1.f;
#endif

171 172 173 174 175 176 177
        part->time_bin = 1;

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

178 179 180 181 182
        ++part;
      }
    }
  }

183
  /* Cell properties */
184
  cell->split = 0;
185 186
  cell->hydro.h_max = h_max;
  cell->hydro.count = count;
187
  cell->hydro.dx_max_part = 0.;
188
  cell->hydro.dx_max_sort = 0.;
189 190 191
  cell->width[0] = size;
  cell->width[1] = size;
  cell->width[2] = size;
192 193 194 195
  cell->loc[0] = offset[0];
  cell->loc[1] = offset[1];
  cell->loc[2] = offset[2];

196
  cell->hydro.ti_old_part = 8;
197 198
  cell->hydro.ti_end_min = 8;
  cell->hydro.ti_end_max = 8;
Matthieu Schaller's avatar
Matthieu Schaller committed
199
  cell->nodeID = NODE_ID;
200

201
  shuffle_particles(cell->hydro.parts, cell->hydro.count);
202

203
  cell->hydro.sorted = 0;
204
  cell->hydro.sort = NULL;
205 206 207 208 209

  return cell;
}

void clean_up(struct cell *ci) {
210
  free(ci->hydro.parts);
211
  free(ci->hydro.sort);
212 213 214 215 216 217 218
  free(ci);
}

/**
 * @brief Initializes all particles field to be ready for a density calculation
 */
void zero_particle_fields(struct cell *c) {
219 220 221 222 223 224 225 226 227 228 229 230
#ifdef SHADOWFAX_SPH
  struct hydro_space hs;
  hs.anchor[0] = 0.;
  hs.anchor[1] = 0.;
  hs.anchor[2] = 0.;
  hs.side[0] = 1.;
  hs.side[1] = 1.;
  hs.side[2] = 1.;
  struct hydro_space *hspointer = &hs;
#else
  struct hydro_space *hspointer = NULL;
#endif
231 232
  for (int pid = 0; pid < c->hydro.count; pid++) {
    hydro_init_part(&c->hydro.parts[pid], hspointer);
233 234 235
  }
}

236 237 238
/**
 * @brief Ends the loop by adding the appropriate coefficients
 */
239
void end_calculation(struct cell *c, const struct cosmology *cosmo) {
240 241
  for (int pid = 0; pid < c->hydro.count; pid++) {
    hydro_end_density(&c->hydro.parts[pid], cosmo);
242 243

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

249 250 251
/**
 * @brief Dump all the particles to a file
 */
252 253
void dump_particle_fields(char *fileName, struct cell *main_cell,
                          struct cell **cells) {
254 255
  FILE *file = fopen(fileName, "w");

256
  /* Write header */
257
  fprintf(file,
258 259
          "# %4s %10s %10s %10s %10s %10s %10s %13s %13s %13s %13s %13s "
          "%13s %13s %13s\n",
260 261
          "ID", "pos_x", "pos_y", "pos_z", "v_x", "v_y", "v_z", "rho", "rho_dh",
          "wcount", "wcount_dh", "div_v", "curl_vx", "curl_vy", "curl_vz");
262

263
  fprintf(file, "# Main cell --------------------------------------------\n");
264

265
  /* Write main cell */
266
  for (int pid = 0; pid < main_cell->hydro.count; pid++) {
267
    fprintf(file,
268 269
            "%6llu %10f %10f %10f %10f %10f %10f %13e %13e %13e %13e %13e "
            "%13e %13e %13e\n",
270 271 272 273 274
            main_cell->hydro.parts[pid].id, main_cell->hydro.parts[pid].x[0],
            main_cell->hydro.parts[pid].x[1], main_cell->hydro.parts[pid].x[2],
            main_cell->hydro.parts[pid].v[0], main_cell->hydro.parts[pid].v[1],
            main_cell->hydro.parts[pid].v[2],
            hydro_get_comoving_density(&main_cell->hydro.parts[pid]),
275
#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
276
            0.f,
Josh Borrow's avatar
Josh Borrow committed
277 278
#elif defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN) || \
    defined(ANARCHY_PU_SPH)
279
            main_cell->hydro.parts[pid].density.pressure_bar_dh,
280
#else
281
            main_cell->hydro.parts[pid].density.rho_dh,
282
#endif
283 284
            main_cell->hydro.parts[pid].density.wcount,
            main_cell->hydro.parts[pid].density.wcount_dh,
285
#if defined(GADGET2_SPH) || defined(HOPKINS_PE_SPH) || \
Josh Borrow's avatar
Josh Borrow committed
286
    defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN)
287 288 289 290
            main_cell->hydro.parts[pid].density.div_v,
            main_cell->hydro.parts[pid].density.rot_v[0],
            main_cell->hydro.parts[pid].density.rot_v[1],
            main_cell->hydro.parts[pid].density.rot_v[2]
Josh Borrow's avatar
Josh Borrow committed
291
#elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
Josh Borrow's avatar
Josh Borrow committed
292 293 294 295 296
            /* this is required because of the variable AV scheme */
            main_cell->hydro.parts[pid].viscosity.div_v,
            main_cell->hydro.parts[pid].density.rot_v[0],
            main_cell->hydro.parts[pid].density.rot_v[1],
            main_cell->hydro.parts[pid].density.rot_v[2]
297 298 299
#else
            0., 0., 0., 0.
#endif
300
    );
301 302
  }

303 304 305 306 307 308
  /* Write all other cells */
  for (int i = 0; i < 3; ++i) {
    for (int j = 0; j < 3; ++j) {
      for (int k = 0; k < 3; ++k) {
        struct cell *cj = cells[i * 9 + j * 3 + k];
        if (cj == main_cell) continue;
309

310 311 312
        fprintf(file,
                "# Offset: [%2d %2d %2d] -----------------------------------\n",
                i - 1, j - 1, k - 1);
313

314
        for (int pjd = 0; pjd < cj->hydro.count; pjd++) {
315 316
          fprintf(
              file,
317 318
              "%6llu %10f %10f %10f %10f %10f %10f %13e %13e %13e %13e %13e "
              "%13e %13e %13e\n",
319 320 321 322 323
              cj->hydro.parts[pjd].id, cj->hydro.parts[pjd].x[0],
              cj->hydro.parts[pjd].x[1], cj->hydro.parts[pjd].x[2],
              cj->hydro.parts[pjd].v[0], cj->hydro.parts[pjd].v[1],
              cj->hydro.parts[pjd].v[2],
              hydro_get_comoving_density(&cj->hydro.parts[pjd]),
324
#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
325 326
              0.f,
#else
327
              main_cell->hydro.parts[pjd].density.rho_dh,
328
#endif
329 330
              cj->hydro.parts[pjd].density.wcount,
              cj->hydro.parts[pjd].density.wcount_dh,
331
#if defined(GADGET2_SPH) || defined(HOPKINS_PE_SPH)
332 333 334 335
              cj->hydro.parts[pjd].density.div_v,
              cj->hydro.parts[pjd].density.rot_v[0],
              cj->hydro.parts[pjd].density.rot_v[1],
              cj->hydro.parts[pjd].density.rot_v[2]
Josh Borrow's avatar
Josh Borrow committed
336
#elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
Josh Borrow's avatar
Josh Borrow committed
337 338 339 340 341
              /* this is required because of the variable AV scheme */
              cj->hydro.parts[pjd].viscosity.div_v,
              cj->hydro.parts[pjd].density.rot_v[0],
              cj->hydro.parts[pjd].density.rot_v[1],
              cj->hydro.parts[pjd].density.rot_v[2]
342 343 344
#else
              0., 0., 0., 0.
#endif
345
          );
346 347
        }
      }
348
    }
349 350 351 352 353
  }
  fclose(file);
}

/* Just a forward declaration... */
354 355
void runner_dopair1_branch_density(struct runner *r, struct cell *ci,
                                   struct cell *cj);
356
void runner_doself1_branch_density(struct runner *r, struct cell *c);
James Willis's avatar
James Willis committed
357 358 359 360 361
void runner_dopair_subset_branch_density(struct runner *r,
                                         struct cell *restrict ci,
                                         struct part *restrict parts_i,
                                         int *restrict ind, int count,
                                         struct cell *restrict cj);
Matthieu Schaller's avatar
Matthieu Schaller committed
362 363 364 365
void runner_doself_subset_branch_density(struct runner *r,
                                         struct cell *restrict ci,
                                         struct part *restrict parts,
                                         int *restrict ind, int count);
366

367
/* And go... */
368
int main(int argc, char *argv[]) {
369

370
#ifdef HAVE_SETAFFINITY
371
  engine_pin();
372 373
#endif

374
  size_t runs = 0, particles = 0;
375
  double h = 1.23485, size = 1., rho = 1.;
376
  double perturbation = 0., h_pert = 0.;
377
  char outputFileNameExtension[100] = "";
378
  char outputFileName[200] = "";
379
  enum velocity_types vel = velocity_zero;
Matthieu Schaller's avatar
Matthieu Schaller committed
380

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

385 386
/* Choke on FP-exceptions */
#ifdef HAVE_FE_ENABLE_EXCEPT
Matthieu Schaller's avatar
Matthieu Schaller committed
387
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
388
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
389

390
  /* Get some randomness going */
391 392
  srand(0);

393
  char c;
394
  while ((c = getopt(argc, argv, "m:s:h:p:n:r:t:d:f:v:")) != -1) {
395 396 397 398
    switch (c) {
      case 'h':
        sscanf(optarg, "%lf", &h);
        break;
399 400 401
      case 'p':
        sscanf(optarg, "%lf", &h_pert);
        break;
402 403 404
      case 's':
        sscanf(optarg, "%lf", &size);
        break;
405
      case 'n':
406 407 408 409 410 411 412 413
        sscanf(optarg, "%zu", &particles);
        break;
      case 'r':
        sscanf(optarg, "%zu", &runs);
        break;
      case 'd':
        sscanf(optarg, "%lf", &perturbation);
        break;
414 415 416
      case 'm':
        sscanf(optarg, "%lf", &rho);
        break;
417 418 419
      case 'f':
        strcpy(outputFileNameExtension, optarg);
        break;
420
      case 'v':
Matthieu Schaller's avatar
Matthieu Schaller committed
421
        sscanf(optarg, "%d", (int *)&vel);
422
        break;
423 424 425 426 427 428
      case '?':
        error("Unknown option.");
        break;
    }
  }

429
  if (h < 0 || particles == 0 || runs == 0) {
430
    printf(
431 432 433 434
        "\nUsage: %s -n PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
        "\nGenerates 27 cells, filled with particles on a Cartesian grid."
        "\nThese are then interacted using runner_dopair1_density() and "
        "runner_doself1_density()."
435
        "\n\nOptions:"
436
        "\n-h DISTANCE=1.2348 - Smoothing length in units of <x>"
437
        "\n-p                 - Random fractional change in h, h=h*random(1,p)"
438 439
        "\n-m rho             - Physical density in the cell"
        "\n-s size            - Physical size of the cell"
440
        "\n-d pert            - Perturbation to apply to the particles [0,1["
441 442
        "\n-v type (0,1,2,3)  - Velocity field: (zero, random, divergent, "
        "rotating)"
443
        "\n-f fileName        - Part of the file name used to save the dumps\n",
444 445 446 447
        argv[0]);
    exit(1);
  }

448
  /* Help users... */
449 450
  message("DOSELF1 function called: %s", DOSELF1_NAME);
  message("DOPAIR1 function called: %s", DOPAIR1_NAME);
451
  message("Vector size: %d", VEC_SIZE);
452
  message("Adiabatic index: ga = %f", hydro_gamma);
453
  message("Hydro implementation: %s", SPH_IMPLEMENTATION);
454 455
  message("Smoothing length: h = %f", h * size);
  message("Kernel:               %s", kernel_name);
456
  message("Neighbour target: N = %f", pow_dimension(h) * kernel_norm);
457 458 459
  message("Density target: rho = %f", rho);
  message("div_v target:   div = %f", vel == 2 ? 3.f : 0.f);
  message("curl_v target: curl = [0., 0., %f]", vel == 3 ? -2.f : 0.f);
460

461
  printf("\n");
462

463 464
  /* Build the infrastructure */
  struct space space;
465 466 467 468
  space.periodic = 1;
  space.dim[0] = 3.;
  space.dim[1] = 3.;
  space.dim[2] = 3.;
469

470
  struct hydro_props hp;
471
  hydro_props_init_no_hydro(&hp);
472 473
  hp.eta_neighbours = h;
  hp.h_tolerance = 1e0;
474
  hp.h_max = FLT_MAX;
475 476
  hp.max_smoothing_iterations = 1;
  hp.CFL_condition = 0.1;
477

478
  struct engine engine;
479 480
  engine.s = &space;
  engine.time = 0.1f;
481
  engine.ti_current = 8;
482
  engine.max_active_bin = num_time_bins;
483
  engine.hydro_properties = &hp;
484
  engine.nodeID = NODE_ID;
485

486 487 488 489
  struct cosmology cosmo;
  cosmology_init_no_cosmo(&cosmo);
  engine.cosmology = &cosmo;

490
  struct runner runner;
491 492
  runner.e = &engine;

493 494 495 496 497 498 499 500
  /* Construct some cells */
  struct cell *cells[27];
  struct cell *main_cell;
  static long long partId = 0;
  for (int i = 0; i < 3; ++i) {
    for (int j = 0; j < 3; ++j) {
      for (int k = 0; k < 3; ++k) {
        double offset[3] = {i * size, j * size, k * size};
501 502 503
        cells[i * 9 + j * 3 + k] =
            make_cell(particles, offset, size, h, rho, &partId, perturbation,
                      vel, h_pert);
504

505
        runner_do_drift_part(&runner, cells[i * 9 + j * 3 + k], 0);
506

Loic Hausammann's avatar
Loic Hausammann committed
507
        runner_do_hydro_sort(&runner, cells[i * 9 + j * 3 + k], 0x1FFF, 0, 0);
508 509 510 511
      }
    }
  }

512
  /* Store the main cell for future use */
513 514
  main_cell = cells[13];

515 516 517
  ticks timings[27];
  for (int i = 0; i < 27; i++) timings[i] = 0;

518
  ticks time = 0;
519 520
  for (size_t i = 0; i < runs; ++i) {
    /* Zero the fields */
521
    for (int j = 0; j < 27; ++j) zero_particle_fields(cells[j]);
522

523
    const ticks tic = getticks();
524

525
#ifdef WITH_VECTORIZATION
526 527 528 529
    runner.ci_cache.count = 0;
    cache_init(&runner.ci_cache, 512);
    runner.cj_cache.count = 0;
    cache_init(&runner.cj_cache, 512);
530 531
#endif

532
#if defined(TEST_DOSELF_SUBSET) || defined(TEST_DOPAIR_SUBSET)
James Willis's avatar
James Willis committed
533 534
    int *pid = NULL;
    int count = 0;
535
    if ((pid = (int *)malloc(sizeof(int) * main_cell->hydro.count)) == NULL)
James Willis's avatar
James Willis committed
536
      error("Can't allocate memory for pid.");
537 538
    for (int k = 0; k < main_cell->hydro.count; k++)
      if (part_is_active(&main_cell->hydro.parts[k], &engine)) {
James Willis's avatar
James Willis committed
539 540 541
        pid[count] = k;
        ++count;
      }
542 543
#endif

544
    /* Run all the pairs */
545 546 547 548
    for (int j = 0; j < 27; ++j) {
      if (cells[j] != main_cell) {
        const ticks sub_tic = getticks();

549
#ifdef TEST_DOPAIR_SUBSET
550
        DOPAIR1_SUBSET(&runner, main_cell, main_cell->hydro.parts, pid, count,
James Willis's avatar
James Willis committed
551
                       cells[j]);
552 553 554
#else
        DOPAIR1(&runner, main_cell, cells[j]);
#endif
555

556
        timings[j] += getticks() - sub_tic;
557 558 559
      }
    }

James Willis's avatar
James Willis committed
560
    /* And now the self-interaction */
561 562
    const ticks self_tic = getticks();

563
#ifdef TEST_DOSELF_SUBSET
564
    DOSELF1_SUBSET(&runner, main_cell, main_cell->hydro.parts, pid, count);
565 566 567
#else
    DOSELF1(&runner, main_cell);
#endif
568

569
    timings[13] += getticks() - self_tic;
570

571
    const ticks toc = getticks();
572 573
    time += toc - tic;

574
    /* Let's get physical ! */
575
    end_calculation(main_cell, &cosmo);
576

577 578
    /* Dump if necessary */
    if (i % 50 == 0) {
lhausamm's avatar
lhausamm committed
579
      sprintf(outputFileName, "swift_dopair_27_%.150s.dat",
580 581
              outputFileNameExtension);
      dump_particle_fields(outputFileName, main_cell, cells);
582 583 584 585
    }
  }

  /* Output timing */
586 587 588 589 590 591 592 593 594 595 596 597 598 599 600
  ticks corner_time = timings[0] + timings[2] + timings[6] + timings[8] +
                      timings[18] + timings[20] + timings[24] + timings[26];

  ticks edge_time = timings[1] + timings[3] + timings[5] + timings[7] +
                    timings[9] + timings[11] + timings[15] + timings[17] +
                    timings[19] + timings[21] + timings[23] + timings[25];

  ticks face_time = timings[4] + timings[10] + timings[12] + timings[14] +
                    timings[16] + timings[22];

  message("Corner calculations took       : %15lli ticks.", corner_time / runs);
  message("Edge calculations took         : %15lli ticks.", edge_time / runs);
  message("Face calculations took         : %15lli ticks.", face_time / runs);
  message("Self calculations took         : %15lli ticks.", timings[13] / runs);
  message("SWIFT calculation took         : %15lli ticks.", time / runs);
601 602 603 604

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

  /* Zero the fields */
605 606 607 608 609 610 611
  for (int i = 0; i < 27; ++i) zero_particle_fields(cells[i]);

  const ticks tic = getticks();

  /* Run all the brute-force pairs */
  for (int j = 0; j < 27; ++j)
    if (cells[j] != main_cell) pairs_all_density(&runner, main_cell, cells[j]);
612

613 614
  /* And now the self-interaction */
  self_all_density(&runner, main_cell);
615

616
  const ticks toc = getticks();
617

618
  /* Let's get physical ! */
619
  end_calculation(main_cell, &cosmo);
620 621

  /* Dump */
lhausamm's avatar
lhausamm committed
622
  sprintf(outputFileName, "brute_force_27_%.150s.dat", outputFileNameExtension);
623
  dump_particle_fields(outputFileName, main_cell, cells);
624 625

  /* Output timing */
Matthieu Schaller's avatar
Matthieu Schaller committed
626
  message("Brute force calculation took : %15lli ticks.", toc - tic);
627 628

  /* Clean things to make the sanitizer happy ... */
629
  for (int i = 0; i < 27; ++i) clean_up(cells[i]);
630

631
#ifdef WITH_VECTORIZATION
632 633
  cache_clean(&runner.ci_cache);
  cache_clean(&runner.cj_cache);
634
#endif
635

636 637
  return 0;
}