test125cells.c 19.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);
100
101
102
103
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
#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;
  }
}

182
183
/**
 * @brief Constructs a cell and all of its particle in a valid state prior to
184
 * a SPH time-step.
185
186
187
188
189
 *
 * @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
190
 * separation.
191
192
 * @param density The density of the fluid.
 * @param partId The running counter of IDs.
193
194
 * @param vel The type of velocity field.
 * @param press The type of pressure field.
195
 */
196
struct cell *make_cell(size_t n, const double offset[3], double size, double h,
197
198
                       double density, long long *partId,
                       enum velocity_field vel, enum pressure_field press) {
199

200
201
202
203
204
205
  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,
206
                     count * sizeof(struct part)) != 0)
207
    error("couldn't allocate particles, no. of particles: %d", (int)count);
208
209
210
  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);
211
  bzero(cell->parts, count * sizeof(struct part));
212
  bzero(cell->xparts, count * sizeof(struct xpart));
213
214
215

  /* Construct the parts */
  struct part *part = cell->parts;
216
  struct xpart *xpart = cell->xparts;
217
218
219
  for (size_t x = 0; x < n; ++x) {
    for (size_t y = 0; y < n; ++y) {
      for (size_t z = 0; z < n; ++z) {
220
221
222
        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;
223
224
        part->h = size * h / (float)n;
        part->mass = density * volume / count;
225
226
227
228
229

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

        part->id = ++(*partId);
230
231
        part->ti_begin = 0;
        part->ti_end = 1;
232

233
234
235
        xpart->v_full[0] = part->v[0];
        xpart->v_full[1] = part->v[1];
        xpart->v_full[2] = part->v[2];
236
        ++part;
237
        ++xpart;
238
239
240
241
242
243
244
245
      }
    }
  }

  /* Cell properties */
  cell->split = 0;
  cell->h_max = h;
  cell->count = count;
246
  cell->gcount = 0;
247
  cell->dx_max = 0.;
Matthieu Schaller's avatar
Matthieu Schaller committed
248
249
250
  cell->width[0] = size;
  cell->width[1] = size;
  cell->width[2] = size;
251
252
253
254
255
256
257
  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;

258
  // shuffle_particles(cell->parts, cell->count);
259
260
261
262
263
264
265
266
267
268

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

  return cell;
}

void clean_up(struct cell *ci) {
  free(ci->parts);
269
  free(ci->xparts);
270
271
272
273
274
275
276
277
  free(ci->sort);
  free(ci);
}

/**
 * @brief Dump all the particles to a file
 */
void dump_particle_fields(char *fileName, struct cell *main_cell,
278
                          struct solution_part *solution, int with_solution) {
279
280
281
282
  FILE *file = fopen(fileName, "w");

  /* Write header */
  fprintf(file,
283
284
285
286
287
          "# %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");
288
289
290
291
292
293

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

  /* Write main cell */
  for (size_t pid = 0; pid < main_cell->count; pid++) {
    fprintf(file,
294
295
296
            "%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",
297
298
299
            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],
300
301
302
303
304
305
306
307
308
            main_cell->parts[pid].v[2], main_cell->parts[pid].h,
            main_cell->parts[pid].rho, main_cell->parts[pid].density.div_v,
            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,
            main_cell->parts[pid].force.v_sig,
309
310
311
#if defined(GADGET2_SPH)
            main_cell->parts[pid].entropy_dt, 0.f
#elif defined(DEFAULT_SPH)
Matthieu Schaller's avatar
Matthieu Schaller committed
312
            0.f, main_cell->parts[pid].force.u_dt
313
#else
Matthieu Schaller's avatar
Matthieu Schaller committed
314
            0.f, 0.f
315
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
316
            );
317
318
  }

319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
  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],
              solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt, 0.f);
335
336
    }
  }
337

338
339
340
341
342
343
  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);
344
345
void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself2_force(struct runner *r, struct cell *ci);
346
347
348

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

350
  size_t runs = 0, particles = 0;
351
  double h = 1.23485, size = 1., rho = 2.5;
352
353
  char outputFileNameExtension[200] = "";
  char outputFileName[200] = "";
354
355
  enum velocity_field vel = velocity_zero;
  enum pressure_field press = pressure_const;
356
357
358
359
360

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

361
362
363
  /* Choke on FP-exceptions */
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);

364
365
366
367
  /* Get some randomness going */
  srand(0);

  char c;
368
  while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:p:")) != -1) {
369
370
371
372
373
374
375
    switch (c) {
      case 'h':
        sscanf(optarg, "%lf", &h);
        break;
      case 's':
        sscanf(optarg, "%lf", &size);
        break;
376
      case 'n':
377
378
379
380
381
382
383
384
385
386
387
        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;
388
389
390
391
392
393
      case 'v':
        sscanf(optarg, "%d", (int *)&vel);
        break;
      case 'p':
        sscanf(optarg, "%d", (int *)&press);
        break;
394
395
396
397
398
399
400
401
      case '?':
        error("Unknown option.");
        break;
    }
  }

  if (h < 0 || particles == 0 || runs == 0) {
    printf(
402
403
404
405
406
        "\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()"
407
408
409
410
        "\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"
411
        "\n-v type (0,1,2,3)  - Velocity field: (zero, constant, divergent, "
412
        "rotating)"
413
        "\n-p type (0,1,2)    - Pressure field: (constant, gradient divergent)"
414
415
416
417
418
419
        "\n-f fileName        - Part of the file name used to save the dumps\n",
        argv[0]);
    exit(1);
  }

  /* Help users... */
420
  message("Adiabatic index: ga = %f", hydro_gamma);
421
422
  message("Smoothing length: h = %f", h * size);
  message("Kernel:               %s", kernel_name);
423
  message("Neighbour target: N = %f", h * h * h * kernel_norm);
424
  message("Density target: rho = %f", rho);
425
426
427
428
429
430
431
432
433
  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");

434
435
436
437
438
439
440
  printf("\n");

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

441
  struct hydro_props hp;
442
443
444
  hp.target_neighbours = h * h * h * kernel_norm;
  hp.delta_neighbours = 1.;
  hp.max_smoothing_iterations = 1;
445

446
  struct engine engine;
447
  engine.hydro_properties = &hp;
448
449
450
451
452
453
454
455
  engine.s = &space;
  engine.time = 0.1f;
  engine.ti_current = 1;

  struct runner runner;
  runner.e = &engine;

  /* Construct some cells */
456
457
  struct cell *cells[125];
  struct cell *inner_cells[27];
458
  struct cell *main_cell;
459
  int count = 0;
460
  static long long partId = 0;
461
462
463
464
465
466
467
468
469
  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] =
470
            make_cell(particles, offset, size, h, rho, &partId, vel, press);
471
472
473
474
475
476

        /* 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++;
        }
477
478
479
480
481
      }
    }
  }

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

484
485
486
487
488
489
  /* 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 */
490
491
492
493
494
  ticks time = 0;
  for (size_t i = 0; i < runs; ++i) {

    const ticks tic = getticks();

495
496
497
498
499
500
501
    /* 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);

/* Do the density calculation */
502
503
#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)

504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
    /* 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*/
535
    for (int j = 0; j < 27; ++j)
536
537
538
      runner_doself1_density(&runner, inner_cells[j]);

#endif
539

540
541
542
543
544
    /* 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 */
#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
545

546
547
548
549
550
551
552
553
554
555
556
557
558
559
    /* 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);
560
561
#endif

562
563
564
    /* Finally, give a gentle kick */
    runner_do_kick(&runner, main_cell, 0);

565
566
567
568
    const ticks toc = getticks();
    time += toc - tic;

    /* Dump if necessary */
569
    if (i == 0) {
570
      sprintf(outputFileName, "swift_dopair_125_%s.dat",
571
              outputFileNameExtension);
572
      dump_particle_fields(outputFileName, main_cell, solution, 0);
573
574
575
576
577
578
    }
  }

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

579
580
581
582
583
584
585
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
  /* 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 */
#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)

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

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

622
#endif
623

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

627
628
/* Do the force calculation */
#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
629

630
631
632
633
  /* 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++) {
634

635
636
637
638
639
640
641
642
643
644
        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);
#endif
645

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

649
  const ticks toc = getticks();
650

651
652
  /* Output timing */
  message("Brute force calculation took : %15lli ticks.", toc - tic);
653

654
655
  sprintf(outputFileName, "brute_force_125_%s.dat", outputFileNameExtension);
  dump_particle_fields(outputFileName, main_cell, solution, 0);
656
657

  /* Clean things to make the sanitizer happy ... */
658
  for (int i = 0; i < 125; ++i) clean_up(cells[i]);
659
  free(solution);
660
661
662

  return 0;
}