test125cells.c 21.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
103
#elif defined(GIZMO_SPH)
  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
198
199
200
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);
  }
}

201
202
/**
 * @brief Constructs a cell and all of its particle in a valid state prior to
203
 * a SPH time-step.
204
205
206
207
208
 *
 * @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
209
 * separation.
210
211
 * @param density The density of the fluid.
 * @param partId The running counter of IDs.
212
213
 * @param vel The type of velocity field.
 * @param press The type of pressure field.
214
 */
215
struct cell *make_cell(size_t n, const double offset[3], double size, double h,
216
217
                       double density, long long *partId,
                       enum velocity_field vel, enum pressure_field press) {
218

219
220
221
222
223
224
  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,
225
                     count * sizeof(struct part)) != 0)
226
    error("couldn't allocate particles, no. of particles: %d", (int)count);
227
228
229
  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);
230
  bzero(cell->parts, count * sizeof(struct part));
231
  bzero(cell->xparts, count * sizeof(struct xpart));
232
233
234

  /* Construct the parts */
  struct part *part = cell->parts;
235
  struct xpart *xpart = cell->xparts;
236
237
238
  for (size_t x = 0; x < n; ++x) {
    for (size_t y = 0; y < n; ++y) {
      for (size_t z = 0; z < n; ++z) {
239
240
241
        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;
242
243
        part->h = size * h / (float)n;
        part->mass = density * volume / count;
244
245
246
247
248

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

        part->id = ++(*partId);
249
250
        part->ti_begin = 0;
        part->ti_end = 1;
251

252
253
254
        xpart->v_full[0] = part->v[0];
        xpart->v_full[1] = part->v[1];
        xpart->v_full[2] = part->v[2];
255
        hydro_first_init_part(part, xpart);
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274

#if defined(GIZMO_SPH)
        part->geometry.volume = volume;
        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 =
            part->primitives.P / hydro_gamma_minus_one * volume +
            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

275
        ++part;
276
        ++xpart;
277
278
279
280
281
282
283
284
      }
    }
  }

  /* Cell properties */
  cell->split = 0;
  cell->h_max = h;
  cell->count = count;
285
  cell->gcount = 0;
286
  cell->dx_max = 0.;
Matthieu Schaller's avatar
Matthieu Schaller committed
287
288
289
  cell->width[0] = size;
  cell->width[1] = size;
  cell->width[2] = size;
290
291
292
293
294
295
296
  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;

297
  // shuffle_particles(cell->parts, cell->count);
298
299
300
301
302
303
304
305
306
307

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

  return cell;
}

void clean_up(struct cell *ci) {
  free(ci->parts);
308
  free(ci->xparts);
309
310
311
312
313
314
315
316
  free(ci->sort);
  free(ci);
}

/**
 * @brief Dump all the particles to a file
 */
void dump_particle_fields(char *fileName, struct cell *main_cell,
317
                          struct solution_part *solution, int with_solution) {
318
319
320
321
  FILE *file = fopen(fileName, "w");

  /* Write header */
  fprintf(file,
322
323
324
325
326
          "# %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");
327
328
329
330
331
332

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

  /* Write main cell */
  for (size_t pid = 0; pid < main_cell->count; pid++) {
    fprintf(file,
333
334
335
            "%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",
336
337
338
            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],
339
            main_cell->parts[pid].v[2], main_cell->parts[pid].h,
Matthieu Schaller's avatar
Matthieu Schaller committed
340
341
            main_cell->parts[pid].rho,
#ifdef MINIMAL_SPH
Matthieu Schaller's avatar
Matthieu Schaller committed
342
            0.f,
Matthieu Schaller's avatar
Matthieu Schaller committed
343
#else
Matthieu Schaller's avatar
Matthieu Schaller committed
344
345
            main_cell->parts[pid].density.div_v,
#endif
346
347
348
349
350
351
            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,
352
#if defined(GADGET2_SPH)
353
354
            main_cell->parts[pid].force.v_sig, main_cell->parts[pid].entropy_dt,
            0.f
355
#elif defined(DEFAULT_SPH)
356
357
            main_cell->parts[pid].force.v_sig, 0.f,
            main_cell->parts[pid].force.u_dt
Matthieu Schaller's avatar
Matthieu Schaller committed
358
#elif defined(MINIMAL_SPH)
359
            main_cell->parts[pid].force.v_sig, 0.f, main_cell->parts[pid].u_dt
360
#else
361
            0.f, 0.f, 0.f
362
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
363
            );
364
365
  }

366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
  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],
381
382
              solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt,
              solution[pid].u_dt);
383
384
    }
  }
385

386
387
388
389
390
391
  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);
392
393
void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself2_force(struct runner *r, struct cell *ci);
394
395
396

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

398
  size_t runs = 0, particles = 0;
399
  double h = 1.23485, size = 1., rho = 2.5;
400
401
  char outputFileNameExtension[200] = "";
  char outputFileName[200] = "";
402
403
  enum velocity_field vel = velocity_zero;
  enum pressure_field press = pressure_const;
404
405
406
407
408

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

409
  /* Choke on FP-exceptions */
410
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
411

412
413
414
415
  /* Get some randomness going */
  srand(0);

  char c;
416
  while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:p:")) != -1) {
417
418
419
420
421
422
423
    switch (c) {
      case 'h':
        sscanf(optarg, "%lf", &h);
        break;
      case 's':
        sscanf(optarg, "%lf", &size);
        break;
424
      case 'n':
425
426
427
428
429
430
431
432
433
434
435
        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;
436
437
438
439
440
441
      case 'v':
        sscanf(optarg, "%d", (int *)&vel);
        break;
      case 'p':
        sscanf(optarg, "%d", (int *)&press);
        break;
442
443
444
445
446
447
448
449
      case '?':
        error("Unknown option.");
        break;
    }
  }

  if (h < 0 || particles == 0 || runs == 0) {
    printf(
450
451
452
453
454
        "\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()"
455
456
457
458
        "\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"
459
        "\n-v type (0,1,2,3)  - Velocity field: (zero, constant, divergent, "
460
        "rotating)"
461
        "\n-p type (0,1,2)    - Pressure field: (constant, gradient divergent)"
462
463
464
465
466
467
        "\n-f fileName        - Part of the file name used to save the dumps\n",
        argv[0]);
    exit(1);
  }

  /* Help users... */
468
  message("Adiabatic index: ga = %f", hydro_gamma);
469
  message("Hydro implementation: %s", SPH_IMPLEMENTATION);
470
471
  message("Smoothing length: h = %f", h * size);
  message("Kernel:               %s", kernel_name);
472
  message("Neighbour target: N = %f", h * h * h * kernel_norm);
473
  message("Density target: rho = %f", rho);
474
475
476
477
478
479
480
481
482
  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");

483
484
485
486
487
488
489
  printf("\n");

  /* Build the infrastructure */
  struct space space;
  space.periodic = 0;
  space.h_max = h;

490
491
492
  struct phys_const prog_const;
  prog_const.const_newton_G = 1.f;

493
  struct hydro_props hp;
494
495
496
  hp.target_neighbours = h * h * h * kernel_norm;
  hp.delta_neighbours = 1.;
  hp.max_smoothing_iterations = 1;
497
  hp.CFL_condition = 0.1;
498

499
  struct engine engine;
500
  engine.hydro_properties = &hp;
501
  engine.physical_constants = &prog_const;
502
503
504
505
506
507
508
509
  engine.s = &space;
  engine.time = 0.1f;
  engine.ti_current = 1;

  struct runner runner;
  runner.e = &engine;

  /* Construct some cells */
510
511
  struct cell *cells[125];
  struct cell *inner_cells[27];
512
  struct cell *main_cell;
513
  int count = 0;
514
  static long long partId = 0;
515
516
517
518
519
520
521
522
523
  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] =
524
            make_cell(particles, offset, size, h, rho, &partId, vel, press);
525
526
527
528
529
530

        /* 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++;
        }
531
532
533
534
535
      }
    }
  }

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

538
539
540
541
542
543
  /* 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 */
544
545
546
547
548
  ticks time = 0;
  for (size_t i = 0; i < runs; ++i) {

    const ticks tic = getticks();

549
550
551
552
553
554
    /* 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);

555
/* Do the density calculation */
556
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
557

558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
    /* 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*/
589
    for (int j = 0; j < 27; ++j)
590
591
592
      runner_doself1_density(&runner, inner_cells[j]);

#endif
593

594
595
596
597
    /* 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 */
598
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
599

600
601
602
603
604
605
606
607
608
609
610
611
612
613
    /* 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);
614
615
#endif

616
617
618
    /* Finally, give a gentle kick */
    runner_do_kick(&runner, main_cell, 0);

619
620
621
622
    const ticks toc = getticks();
    time += toc - tic;

    /* Dump if necessary */
623
    if (i == 0) {
624
      sprintf(outputFileName, "swift_dopair_125_%s.dat",
625
              outputFileNameExtension);
626
      dump_particle_fields(outputFileName, main_cell, solution, 0);
627
628
629
630
631
632
    }
  }

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

633
634
635
  for (int j = 0; j < 125; ++j)
    reset_particles(cells[j], vel, press, size, rho);

636
637
638
639
640
641
642
643
  /* 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 */
644
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
645

646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
  /* 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);
            }
          }
        }
      }
    }
  }
675

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

679
#endif
680

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

684
/* Do the force calculation */
685
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
686

687
688
689
690
  /* 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++) {
691

692
693
694
695
696
697
698
699
700
        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);
701

702
#endif
703

704
705
  /* Finally, give a gentle kick */
  runner_do_kick(&runner, main_cell, 0);
706

707
  const ticks toc = getticks();
708

709
710
  /* Output timing */
  message("Brute force calculation took : %15lli ticks.", toc - tic);
711

712
713
  sprintf(outputFileName, "brute_force_125_%s.dat", outputFileNameExtension);
  dump_particle_fields(outputFileName, main_cell, solution, 0);
714
715

  /* Clean things to make the sanitizer happy ... */
716
  for (int i = 0; i < 125; ++i) clean_up(cells[i]);
717
  free(solution);
718
719
720

  return 0;
}