test125cells.c 22.6 KB
Newer Older
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Copyright (C) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 *
 * 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/>.
 *
 ******************************************************************************/

20
21
22
23
/* Config parameters. */
#include "../config.h"

/* Some standard headers. */
24
25
26
27
28
29
#include <fenv.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

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

33
34
35
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
enum velocity_field {
  velocity_zero,
  velocity_const,
  velocity_divergent,
  velocity_rotating
};

enum pressure_field { pressure_const, pressure_gradient, pressure_divergent };

void set_velocity(struct part *part, enum velocity_field vel, float size) {

  switch (vel) {
    case velocity_zero:
      part->v[0] = 0.f;
      part->v[1] = 0.f;
      part->v[2] = 0.f;
      break;
    case velocity_const:
      part->v[0] = 1.f;
      part->v[1] = 0.f;
      part->v[2] = 0.f;
      break;
    case velocity_divergent:
      part->v[0] = part->x[0] - 2.5 * size;
      part->v[1] = part->x[1] - 2.5 * size;
      part->v[2] = part->x[2] - 2.5 * size;
      break;
    case velocity_rotating:
      part->v[0] = part->x[1];
      part->v[1] = -part->x[0];
      part->v[2] = 0.f;
      break;
  }
}

float get_pressure(double x[3], enum pressure_field press, float size) {

  float r2 = 0.;
  float dx[3] = {0.f};

  switch (press) {
    case pressure_const:
      return 1.5f;
      break;
    case pressure_gradient:
      return 1.5f * x[0]; /* gradient along x */
      break;
    case pressure_divergent:
      dx[0] = x[0] - 2.5 * size;
      dx[1] = x[1] - 2.5 * size;
      dx[2] = x[2] - 2.5 * size;
      r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
      return sqrt(r2) + 1.5f;
      break;
  }
  return 0.f;
}

void set_energy_state(struct part *part, enum pressure_field press, float size,
                      float density) {

  const float pressure = get_pressure(part->x, press, size);

#if defined(GADGET2_SPH)
  part->entropy = pressure / pow_gamma(density);
98
99
#elif defined(DEFAULT_SPH)
  part->u = pressure / (hydro_gamma_minus_one * density);
Matthieu Schaller's avatar
Matthieu Schaller committed
100
101
#elif defined(MINIMAL_SPH)
  part->u = pressure / (hydro_gamma_minus_one * density);
102
#elif defined(GIZMO_SPH) || defined(SHADOWSWIFT)
103
  part->primitives.P = pressure;
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
#else
  error("Need to define pressure here !");
#endif
}

struct solution_part {

  long long id;
  double x[3];
  float v[3];
  float a_hydro[3];
  float h;
  float rho;
  float div_v;
  float S;
  float u;
  float P;
  float c;
  float h_dt;
  float v_sig;
  float S_dt;
  float u_dt;
};

void get_solution(const struct cell *main_cell, struct solution_part *solution,
                  float density, enum velocity_field vel,
                  enum pressure_field press, float size) {

  for (size_t i = 0; i < main_cell->count; ++i) {

    solution[i].id = main_cell->parts[i].id;

    solution[i].x[0] = main_cell->parts[i].x[0];
    solution[i].x[1] = main_cell->parts[i].x[1];
    solution[i].x[2] = main_cell->parts[i].x[2];

    solution[i].v[0] = main_cell->parts[i].v[0];
    solution[i].v[1] = main_cell->parts[i].v[1];
    solution[i].v[2] = main_cell->parts[i].v[2];

    solution[i].h = main_cell->parts[i].h;

    solution[i].rho = density;

    solution[i].P = get_pressure(solution[i].x, press, size);
    solution[i].u = solution[i].P / (solution[i].rho * hydro_gamma_minus_one);
    solution[i].S = solution[i].P / pow_gamma(solution[i].rho);
    solution[i].c = sqrt(hydro_gamma * solution[i].P / solution[i].rho);

    if (vel == velocity_divergent)
      solution[i].div_v = 3.f;
    else
      solution[i].div_v = 0.f;

    solution[i].h_dt = solution[i].h * solution[i].div_v / 3.;

    float gradP[3] = {0.f};
    if (press == pressure_gradient) {
      gradP[0] = 1.5f;
      gradP[1] = 0.f;
      gradP[2] = 0.f;
    } else if (press == pressure_divergent) {
      float dx[3];
      dx[0] = solution[i].x[0] - 2.5 * size;
      dx[1] = solution[i].x[1] - 2.5 * size;
      dx[2] = solution[i].x[2] - 2.5 * size;
      float r = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
      if (r > 0.) {
        gradP[0] = dx[0] / r;
        gradP[1] = dx[1] / r;
        gradP[2] = dx[2] / r;
      }
    }

    solution[i].a_hydro[0] = -gradP[0] / solution[i].rho;
    solution[i].a_hydro[1] = -gradP[1] / solution[i].rho;
    solution[i].a_hydro[2] = -gradP[2] / solution[i].rho;

    solution[i].v_sig = 2.f * solution[i].c;
183
184
185

    solution[i].S_dt = 0.f;
    solution[i].u_dt = -(solution[i].P / solution[i].rho) * solution[i].div_v;
186
187
188
  }
}

189
190
191
192
193
194
195
196
197
void reset_particles(struct cell *c, enum velocity_field vel,
                     enum pressure_field press, float size, float density) {

  for (size_t i = 0; i < c->count; ++i) {

    struct part *p = &c->parts[i];

    set_velocity(p, vel, size);
    set_energy_state(p, press, size, density);
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
198

199
#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
200
201
202
203
204
205
    float volume = p->mass / density;
#if defined(GIZMO_SPH)
    p->geometry.volume = volume;
#else
    p->cell.volume = volume;
#endif
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
206
207
208
209
210
211
212
213
214
    p->primitives.rho = density;
    p->primitives.v[0] = p->v[0];
    p->primitives.v[1] = p->v[1];
    p->primitives.v[2] = p->v[2];
    p->conserved.mass = p->mass;
    p->conserved.momentum[0] = p->conserved.mass * p->v[0];
    p->conserved.momentum[1] = p->conserved.mass * p->v[1];
    p->conserved.momentum[2] = p->conserved.mass * p->v[2];
    p->conserved.energy =
215
        p->primitives.P / hydro_gamma_minus_one * volume +
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
216
217
218
219
220
        0.5f * (p->conserved.momentum[0] * p->conserved.momentum[0] +
                p->conserved.momentum[1] * p->conserved.momentum[1] +
                p->conserved.momentum[2] * p->conserved.momentum[2]) /
            p->conserved.mass;
#endif
221
222
223
  }
}

224
225
/**
 * @brief Constructs a cell and all of its particle in a valid state prior to
226
 * a SPH time-step.
227
228
229
230
231
 *
 * @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
232
 * separation.
233
234
 * @param density The density of the fluid.
 * @param partId The running counter of IDs.
235
236
 * @param vel The type of velocity field.
 * @param press The type of pressure field.
237
 */
238
struct cell *make_cell(size_t n, const double offset[3], double size, double h,
239
240
                       double density, long long *partId,
                       enum velocity_field vel, enum pressure_field press) {
241

242
243
244
245
246
247
  const size_t count = n * n * n;
  const double volume = size * size * size;
  struct cell *cell = malloc(sizeof(struct cell));
  bzero(cell, sizeof(struct cell));

  if (posix_memalign((void **)&cell->parts, part_align,
248
                     count * sizeof(struct part)) != 0)
249
    error("couldn't allocate particles, no. of particles: %d", (int)count);
250
251
252
  if (posix_memalign((void **)&cell->xparts, xpart_align,
                     count * sizeof(struct xpart)) != 0)
    error("couldn't allocate particles, no. of x-particles: %d", (int)count);
253
  bzero(cell->parts, count * sizeof(struct part));
254
  bzero(cell->xparts, count * sizeof(struct xpart));
255
256
257

  /* Construct the parts */
  struct part *part = cell->parts;
258
  struct xpart *xpart = cell->xparts;
259
260
261
  for (size_t x = 0; x < n; ++x) {
    for (size_t y = 0; y < n; ++y) {
      for (size_t z = 0; z < n; ++z) {
262
263
264
        part->x[0] = offset[0] + size * (x + 0.5) / (float)n;
        part->x[1] = offset[1] + size * (y + 0.5) / (float)n;
        part->x[2] = offset[2] + size * (z + 0.5) / (float)n;
265
266
        part->h = size * h / (float)n;
        part->mass = density * volume / count;
267
268
269
270
271

        set_velocity(part, vel, size);
        set_energy_state(part, press, size, density);

        part->id = ++(*partId);
272
273
        part->ti_begin = 0;
        part->ti_end = 1;
274

275
        hydro_first_init_part(part, xpart);
276

277
#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
278
279
280
281
282
283
        float volume = part->mass / density;
#ifdef GIZMO_SPH
        part->geometry.volume = volume;
#else
        part->cell.volume = volume;
#endif
284
285
286
287
288
289
290
291
292
        part->primitives.rho = density;
        part->primitives.v[0] = part->v[0];
        part->primitives.v[1] = part->v[1];
        part->primitives.v[2] = part->v[2];
        part->conserved.mass = part->mass;
        part->conserved.momentum[0] = part->conserved.mass * part->v[0];
        part->conserved.momentum[1] = part->conserved.mass * part->v[1];
        part->conserved.momentum[2] = part->conserved.mass * part->v[2];
        part->conserved.energy =
293
            part->primitives.P / hydro_gamma_minus_one * volume +
294
295
296
297
298
299
            0.5f * (part->conserved.momentum[0] * part->conserved.momentum[0] +
                    part->conserved.momentum[1] * part->conserved.momentum[1] +
                    part->conserved.momentum[2] * part->conserved.momentum[2]) /
                part->conserved.mass;
#endif

300
        ++part;
301
        ++xpart;
302
303
304
305
306
307
308
309
      }
    }
  }

  /* Cell properties */
  cell->split = 0;
  cell->h_max = h;
  cell->count = count;
310
  cell->gcount = 0;
311
  cell->dx_max = 0.;
Matthieu Schaller's avatar
Matthieu Schaller committed
312
313
314
  cell->width[0] = size;
  cell->width[1] = size;
  cell->width[2] = size;
315
316
317
318
319
320
321
  cell->loc[0] = offset[0];
  cell->loc[1] = offset[1];
  cell->loc[2] = offset[2];

  cell->ti_end_min = 1;
  cell->ti_end_max = 1;

322
  // shuffle_particles(cell->parts, cell->count);
323
324
325
326
327
328
329
330
331
332

  cell->sorted = 0;
  cell->sort = NULL;
  cell->sortsize = 0;

  return cell;
}

void clean_up(struct cell *ci) {
  free(ci->parts);
333
  free(ci->xparts);
334
335
336
337
338
339
340
341
  free(ci->sort);
  free(ci);
}

/**
 * @brief Dump all the particles to a file
 */
void dump_particle_fields(char *fileName, struct cell *main_cell,
342
                          struct solution_part *solution, int with_solution) {
343
344
345
346
  FILE *file = fopen(fileName, "w");

  /* Write header */
  fprintf(file,
347
348
349
350
351
          "# %4s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s "
          "%8s %8s %8s %8s %8s\n",
          "ID", "pos_x", "pos_y", "pos_z", "v_x", "v_y", "v_z", "h", "rho",
          "div_v", "S", "u", "P", "c", "a_x", "a_y", "a_z", "h_dt", "v_sig",
          "dS/dt", "du/dt");
352
353
354
355
356
357

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

  /* Write main cell */
  for (size_t pid = 0; pid < main_cell->count; pid++) {
    fprintf(file,
358
359
360
            "%6llu %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f "
            "%8.5f "
            "%8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f\n",
361
362
363
            main_cell->parts[pid].id, main_cell->parts[pid].x[0],
            main_cell->parts[pid].x[1], main_cell->parts[pid].x[2],
            main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
364
            main_cell->parts[pid].v[2], main_cell->parts[pid].h,
Matthieu Schaller's avatar
Matthieu Schaller committed
365
            main_cell->parts[pid].rho,
366
#if defined(MINIMAL_SPH) || defined(SHADOWSWIFT)
Matthieu Schaller's avatar
Matthieu Schaller committed
367
            0.f,
Matthieu Schaller's avatar
Matthieu Schaller committed
368
#else
Matthieu Schaller's avatar
Matthieu Schaller committed
369
370
            main_cell->parts[pid].density.div_v,
#endif
371
372
373
374
375
376
            hydro_get_entropy(&main_cell->parts[pid], 0.f),
            hydro_get_internal_energy(&main_cell->parts[pid], 0.f),
            hydro_get_pressure(&main_cell->parts[pid], 0.f),
            hydro_get_soundspeed(&main_cell->parts[pid], 0.f),
            main_cell->parts[pid].a_hydro[0], main_cell->parts[pid].a_hydro[1],
            main_cell->parts[pid].a_hydro[2], main_cell->parts[pid].force.h_dt,
377
#if defined(GADGET2_SPH)
378
379
            main_cell->parts[pid].force.v_sig, main_cell->parts[pid].entropy_dt,
            0.f
380
#elif defined(DEFAULT_SPH)
381
382
            main_cell->parts[pid].force.v_sig, 0.f,
            main_cell->parts[pid].force.u_dt
Matthieu Schaller's avatar
Matthieu Schaller committed
383
#elif defined(MINIMAL_SPH)
384
            main_cell->parts[pid].force.v_sig, 0.f, main_cell->parts[pid].u_dt
385
#else
386
            0.f, 0.f, 0.f
387
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
388
            );
389
390
  }

391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
  if (with_solution) {

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

    for (size_t pid = 0; pid < main_cell->count; pid++) {
      fprintf(file,
              "%6llu %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f "
              "%8.5f %8.5f "
              "%8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f\n",
              solution[pid].id, solution[pid].x[0], solution[pid].x[1],
              solution[pid].x[2], solution[pid].v[0], solution[pid].v[1],
              solution[pid].v[2], solution[pid].h, solution[pid].rho,
              solution[pid].div_v, solution[pid].S, solution[pid].u,
              solution[pid].P, solution[pid].c, solution[pid].a_hydro[0],
              solution[pid].a_hydro[1], solution[pid].a_hydro[2],
406
407
              solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt,
              solution[pid].u_dt);
408
409
    }
  }
410

411
412
413
414
415
416
  fclose(file);
}

/* Just a forward declaration... */
void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself1_density(struct runner *r, struct cell *ci);
417
418
void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself2_force(struct runner *r, struct cell *ci);
419
420
421

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

423
  size_t runs = 0, particles = 0;
424
  double h = 1.23485, size = 1., rho = 2.5;
425
426
  char outputFileNameExtension[200] = "";
  char outputFileName[200] = "";
427
428
  enum velocity_field vel = velocity_zero;
  enum pressure_field press = pressure_const;
429
430
431
432
433

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

434
  /* Choke on FP-exceptions */
435
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
436

437
438
439
440
  /* Get some randomness going */
  srand(0);

  char c;
441
  while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:p:")) != -1) {
442
443
444
445
446
447
448
    switch (c) {
      case 'h':
        sscanf(optarg, "%lf", &h);
        break;
      case 's':
        sscanf(optarg, "%lf", &size);
        break;
449
      case 'n':
450
451
452
453
454
455
456
457
458
459
460
        sscanf(optarg, "%zu", &particles);
        break;
      case 'r':
        sscanf(optarg, "%zu", &runs);
        break;
      case 'm':
        sscanf(optarg, "%lf", &rho);
        break;
      case 'f':
        strcpy(outputFileNameExtension, optarg);
        break;
461
462
463
464
465
466
      case 'v':
        sscanf(optarg, "%d", (int *)&vel);
        break;
      case 'p':
        sscanf(optarg, "%d", (int *)&press);
        break;
467
468
469
470
471
472
473
474
      case '?':
        error("Unknown option.");
        break;
    }
  }

  if (h < 0 || particles == 0 || runs == 0) {
    printf(
475
476
477
478
479
        "\nUsage: %s -n PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
        "\nGenerates 125 cells, filled with particles on a Cartesian grid."
        "\nThese are then interacted using runner_dopair1_density() and "
        "runner_doself1_density() followed by runner_dopair2_force() and "
        "runner_doself2_force()"
480
481
482
483
        "\n\nOptions:"
        "\n-h DISTANCE=1.2348 - Smoothing length in units of <x>"
        "\n-m rho             - Physical density in the cell"
        "\n-s size            - Physical size of the cell"
484
        "\n-v type (0,1,2,3)  - Velocity field: (zero, constant, divergent, "
485
        "rotating)"
486
        "\n-p type (0,1,2)    - Pressure field: (constant, gradient divergent)"
487
488
489
490
491
492
        "\n-f fileName        - Part of the file name used to save the dumps\n",
        argv[0]);
    exit(1);
  }

  /* Help users... */
493
  message("Adiabatic index: ga = %f", hydro_gamma);
494
  message("Hydro implementation: %s", SPH_IMPLEMENTATION);
495
496
  message("Smoothing length: h = %f", h * size);
  message("Kernel:               %s", kernel_name);
497
  message("Neighbour target: N = %f", pow_dimension(h) * kernel_norm);
498
  message("Density target: rho = %f", rho);
499
500
501
502
503
504
505
506
507
  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);
  if (press == pressure_const)
    message("P field constant");
  else if (press == pressure_gradient)
    message("P field gradient");
  else
    message("P field divergent");

508
509
  printf("\n");

510
511
512
513
#if !defined(HYDRO_DIMENSION_3D)
  message("test125cells only useful in 3D. Change parameters in const.h !");
  return 1;
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
514

515
516
517
518
  /* Build the infrastructure */
  struct space space;
  space.periodic = 0;

519
520
521
  struct phys_const prog_const;
  prog_const.const_newton_G = 1.f;

522
  struct hydro_props hp;
523
524
  hp.target_neighbours = pow_dimension(h) * kernel_norm;
  hp.delta_neighbours = 2.;
525
  hp.max_smoothing_iterations = 1;
526
  hp.CFL_condition = 0.1;
527

528
  struct engine engine;
529
  engine.hydro_properties = &hp;
530
  engine.physical_constants = &prog_const;
531
532
533
534
535
536
537
538
  engine.s = &space;
  engine.time = 0.1f;
  engine.ti_current = 1;

  struct runner runner;
  runner.e = &engine;

  /* Construct some cells */
539
540
  struct cell *cells[125];
  struct cell *inner_cells[27];
541
  struct cell *main_cell;
542
  int count = 0;
543
  static long long partId = 0;
544
545
546
547
548
549
550
551
552
  for (int i = 0; i < 5; ++i) {
    for (int j = 0; j < 5; ++j) {
      for (int k = 0; k < 5; ++k) {

        /* Position of the cell */
        const double offset[3] = {i * size, j * size, k * size};

        /* Construct it */
        cells[i * 25 + j * 5 + k] =
553
            make_cell(particles, offset, size, h, rho, &partId, vel, press);
554
555
556
557
558
559

        /* Store the inner cells */
        if (i > 0 && i < 4 && j > 0 && j < 4 && k > 0 && k < 4) {
          inner_cells[count] = cells[i * 25 + j * 5 + k];
          count++;
        }
560
561
562
563
564
      }
    }
  }

  /* Store the main cell for future use */
565
  main_cell = cells[62];
566

567
568
569
570
571
572
  /* Construct the real solution */
  struct solution_part *solution =
      malloc(main_cell->count * sizeof(struct solution_part));
  get_solution(main_cell, solution, rho, vel, press, size);

  /* Start the test */
573
574
575
576
577
  ticks time = 0;
  for (size_t i = 0; i < runs; ++i) {

    const ticks tic = getticks();

578
579
580
581
582
583
    /* First, sort stuff */
    for (int j = 0; j < 125; ++j) runner_do_sort(&runner, cells[j], 0x1FFF, 0);

    /* Initialise the particles */
    for (int j = 0; j < 125; ++j) runner_do_init(&runner, cells[j], 0);

584
/* Do the density calculation */
585
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
586

587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
    /* Run all the pairs (only once !)*/
    for (int i = 0; i < 5; i++) {
      for (int j = 0; j < 5; j++) {
        for (int k = 0; k < 5; k++) {

          struct cell *ci = cells[i * 25 + j * 5 + k];

          for (int ii = -1; ii < 2; ii++) {
            int iii = i + ii;
            if (iii < 0 || iii >= 5) continue;
            iii = (iii + 5) % 5;
            for (int jj = -1; jj < 2; jj++) {
              int jjj = j + jj;
              if (jjj < 0 || jjj >= 5) continue;
              jjj = (jjj + 5) % 5;
              for (int kk = -1; kk < 2; kk++) {
                int kkk = k + kk;
                if (kkk < 0 || kkk >= 5) continue;
                kkk = (kkk + 5) % 5;

                struct cell *cj = cells[iii * 25 + jjj * 5 + kkk];

                if (cj > ci) runner_dopair1_density(&runner, ci, cj);
              }
            }
          }
        }
      }
    }

    /* And now the self-interaction for the central cells*/
618
    for (int j = 0; j < 27; ++j)
619
620
621
      runner_doself1_density(&runner, inner_cells[j]);

#endif
622

623
624
625
626
    /* Ghost to finish everything on the central cells */
    for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j]);

/* Do the force calculation */
627
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
628

629
630
631
632
633
634
635
636
637
638
639
640
641
642
    /* Do the pairs (for the central 27 cells) */
    for (int i = 1; i < 4; i++) {
      for (int j = 1; j < 4; j++) {
        for (int k = 1; k < 4; k++) {

          struct cell *cj = cells[i * 25 + j * 5 + k];

          if (main_cell != cj) runner_dopair2_force(&runner, main_cell, cj);
        }
      }
    }

    /* And now the self-interaction for the main cell */
    runner_doself2_force(&runner, main_cell);
643
644
#endif

645
646
647
    /* Finally, give a gentle kick */
    runner_do_kick(&runner, main_cell, 0);

648
649
650
651
    const ticks toc = getticks();
    time += toc - tic;

    /* Dump if necessary */
652
    if (i == 0) {
653
      sprintf(outputFileName, "swift_dopair_125_%s.dat",
654
              outputFileNameExtension);
655
      dump_particle_fields(outputFileName, main_cell, solution, 0);
656
657
658
659
660
661
    }
  }

  /* Output timing */
  message("SWIFT calculation took       : %15lli ticks.", time / runs);

662
663
664
  for (int j = 0; j < 125; ++j)
    reset_particles(cells[j], vel, press, size, rho);

665
666
667
668
669
670
671
672
  /* NOW BRUTE-FORCE CALCULATION */

  const ticks tic = getticks();

  /* Initialise the particles */
  for (int j = 0; j < 125; ++j) runner_do_init(&runner, cells[j], 0);

/* Do the density calculation */
673
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
674

675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
  /* Run all the pairs (only once !)*/
  for (int i = 0; i < 5; i++) {
    for (int j = 0; j < 5; j++) {
      for (int k = 0; k < 5; k++) {

        struct cell *ci = cells[i * 25 + j * 5 + k];

        for (int ii = -1; ii < 2; ii++) {
          int iii = i + ii;
          if (iii < 0 || iii >= 5) continue;
          iii = (iii + 5) % 5;
          for (int jj = -1; jj < 2; jj++) {
            int jjj = j + jj;
            if (jjj < 0 || jjj >= 5) continue;
            jjj = (jjj + 5) % 5;
            for (int kk = -1; kk < 2; kk++) {
              int kkk = k + kk;
              if (kkk < 0 || kkk >= 5) continue;
              kkk = (kkk + 5) % 5;

              struct cell *cj = cells[iii * 25 + jjj * 5 + kkk];

              if (cj > ci) pairs_all_density(&runner, ci, cj);
            }
          }
        }
      }
    }
  }
704

705
706
  /* And now the self-interaction for the central cells*/
  for (int j = 0; j < 27; ++j) self_all_density(&runner, inner_cells[j]);
707

708
#endif
709

710
711
  /* Ghost to finish everything on the central cells */
  for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j]);
712

713
/* Do the force calculation */
714
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
715

716
717
718
719
  /* Do the pairs (for the central 27 cells) */
  for (int i = 1; i < 4; i++) {
    for (int j = 1; j < 4; j++) {
      for (int k = 1; k < 4; k++) {
720

721
722
723
724
725
726
727
728
729
        struct cell *cj = cells[i * 25 + j * 5 + k];

        if (main_cell != cj) pairs_all_force(&runner, main_cell, cj);
      }
    }
  }

  /* And now the self-interaction for the main cell */
  self_all_force(&runner, main_cell);
730

731
#endif
732

733
734
  /* Finally, give a gentle kick */
  runner_do_kick(&runner, main_cell, 0);
735

736
  const ticks toc = getticks();
737

738
739
  /* Output timing */
  message("Brute force calculation took : %15lli ticks.", toc - tic);
740

741
742
  sprintf(outputFileName, "brute_force_125_%s.dat", outputFileNameExtension);
  dump_particle_fields(outputFileName, main_cell, solution, 0);
743
744

  /* Clean things to make the sanitizer happy ... */
745
  for (int i = 0; i < 125; ++i) clean_up(cells[i]);
746
  free(solution);
747
748
749

  return 0;
}