test125cells.c 22.4 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
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
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215

#if defined(GIZMO_SPH)
    p->geometry.volume = p->mass / density;
    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 =
        p->primitives.P / hydro_gamma_minus_one * p->geometry.volume +
        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
216
217
218
  }
}

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

237
238
239
240
241
242
  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,
243
                     count * sizeof(struct part)) != 0)
244
    error("couldn't allocate particles, no. of particles: %d", (int)count);
245
246
247
  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);
248
  bzero(cell->parts, count * sizeof(struct part));
249
  bzero(cell->xparts, count * sizeof(struct xpart));
250
251
252

  /* Construct the parts */
  struct part *part = cell->parts;
253
  struct xpart *xpart = cell->xparts;
254
255
256
  for (size_t x = 0; x < n; ++x) {
    for (size_t y = 0; y < n; ++y) {
      for (size_t z = 0; z < n; ++z) {
257
258
259
        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;
260
261
        part->h = size * h / (float)n;
        part->mass = density * volume / count;
262
263
264
265
266

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

        part->id = ++(*partId);
267
268
        part->ti_begin = 0;
        part->ti_end = 1;
269

270
271
272
        xpart->v_full[0] = part->v[0];
        xpart->v_full[1] = part->v[1];
        xpart->v_full[2] = part->v[2];
273
        hydro_first_init_part(part, xpart);
274
275

#if defined(GIZMO_SPH)
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
276
        part->geometry.volume = part->mass / density;
277
278
279
280
281
282
283
284
285
        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 =
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
286
            part->primitives.P / hydro_gamma_minus_one * part->geometry.volume +
287
288
289
290
291
292
            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

293
        ++part;
294
        ++xpart;
295
296
297
298
299
300
301
302
      }
    }
  }

  /* Cell properties */
  cell->split = 0;
  cell->h_max = h;
  cell->count = count;
303
  cell->gcount = 0;
304
  cell->dx_max = 0.;
Matthieu Schaller's avatar
Matthieu Schaller committed
305
306
307
  cell->width[0] = size;
  cell->width[1] = size;
  cell->width[2] = size;
308
309
310
311
312
313
314
  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;

315
  // shuffle_particles(cell->parts, cell->count);
316
317
318
319
320
321
322
323
324
325

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

  return cell;
}

void clean_up(struct cell *ci) {
  free(ci->parts);
326
  free(ci->xparts);
327
328
329
330
331
332
333
334
  free(ci->sort);
  free(ci);
}

/**
 * @brief Dump all the particles to a file
 */
void dump_particle_fields(char *fileName, struct cell *main_cell,
335
                          struct solution_part *solution, int with_solution) {
336
337
338
339
  FILE *file = fopen(fileName, "w");

  /* Write header */
  fprintf(file,
340
341
342
343
344
          "# %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");
345
346
347
348
349
350

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

  /* Write main cell */
  for (size_t pid = 0; pid < main_cell->count; pid++) {
    fprintf(file,
351
352
353
            "%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",
354
355
356
            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],
357
            main_cell->parts[pid].v[2], main_cell->parts[pid].h,
Matthieu Schaller's avatar
Matthieu Schaller committed
358
359
            main_cell->parts[pid].rho,
#ifdef MINIMAL_SPH
Matthieu Schaller's avatar
Matthieu Schaller committed
360
            0.f,
Matthieu Schaller's avatar
Matthieu Schaller committed
361
#else
Matthieu Schaller's avatar
Matthieu Schaller committed
362
363
            main_cell->parts[pid].density.div_v,
#endif
364
365
366
367
368
369
            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,
370
#if defined(GADGET2_SPH)
371
372
            main_cell->parts[pid].force.v_sig, main_cell->parts[pid].entropy_dt,
            0.f
373
#elif defined(DEFAULT_SPH)
374
375
            main_cell->parts[pid].force.v_sig, 0.f,
            main_cell->parts[pid].force.u_dt
Matthieu Schaller's avatar
Matthieu Schaller committed
376
#elif defined(MINIMAL_SPH)
377
            main_cell->parts[pid].force.v_sig, 0.f, main_cell->parts[pid].u_dt
378
#else
379
            0.f, 0.f, 0.f
380
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
381
            );
382
383
  }

384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
  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],
399
400
              solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt,
              solution[pid].u_dt);
401
402
    }
  }
403

404
405
406
407
408
409
  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);
410
411
void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself2_force(struct runner *r, struct cell *ci);
412
413
414

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

416
  size_t runs = 0, particles = 0;
417
  double h = 1.23485, size = 1., rho = 2.5;
418
419
  char outputFileNameExtension[200] = "";
  char outputFileName[200] = "";
420
421
  enum velocity_field vel = velocity_zero;
  enum pressure_field press = pressure_const;
422
423
424
425
426

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

427
  /* Choke on FP-exceptions */
428
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
429

430
431
432
433
  /* Get some randomness going */
  srand(0);

  char c;
434
  while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:p:")) != -1) {
435
436
437
438
439
440
441
    switch (c) {
      case 'h':
        sscanf(optarg, "%lf", &h);
        break;
      case 's':
        sscanf(optarg, "%lf", &size);
        break;
442
      case 'n':
443
444
445
446
447
448
449
450
451
452
453
        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;
454
455
456
457
458
459
      case 'v':
        sscanf(optarg, "%d", (int *)&vel);
        break;
      case 'p':
        sscanf(optarg, "%d", (int *)&press);
        break;
460
461
462
463
464
465
466
467
      case '?':
        error("Unknown option.");
        break;
    }
  }

  if (h < 0 || particles == 0 || runs == 0) {
    printf(
468
469
470
471
472
        "\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()"
473
474
475
476
        "\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"
477
        "\n-v type (0,1,2,3)  - Velocity field: (zero, constant, divergent, "
478
        "rotating)"
479
        "\n-p type (0,1,2)    - Pressure field: (constant, gradient divergent)"
480
481
482
483
484
485
        "\n-f fileName        - Part of the file name used to save the dumps\n",
        argv[0]);
    exit(1);
  }

  /* Help users... */
486
  message("Adiabatic index: ga = %f", hydro_gamma);
487
  message("Hydro implementation: %s", SPH_IMPLEMENTATION);
488
489
  message("Smoothing length: h = %f", h * size);
  message("Kernel:               %s", kernel_name);
490
  message("Neighbour target: N = %f", h * h * h * kernel_norm);
491
  message("Density target: rho = %f", rho);
492
493
494
495
496
497
498
499
500
  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");

501
502
503
504
505
506
507
  printf("\n");

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

508
509
510
  struct phys_const prog_const;
  prog_const.const_newton_G = 1.f;

511
  struct hydro_props hp;
512
513
514
  hp.target_neighbours = h * h * h * kernel_norm;
  hp.delta_neighbours = 1.;
  hp.max_smoothing_iterations = 1;
515
  hp.CFL_condition = 0.1;
516

517
  struct engine engine;
518
  engine.hydro_properties = &hp;
519
  engine.physical_constants = &prog_const;
520
521
522
523
524
525
526
527
  engine.s = &space;
  engine.time = 0.1f;
  engine.ti_current = 1;

  struct runner runner;
  runner.e = &engine;

  /* Construct some cells */
528
529
  struct cell *cells[125];
  struct cell *inner_cells[27];
530
  struct cell *main_cell;
531
  int count = 0;
532
  static long long partId = 0;
533
534
535
536
537
538
539
540
541
  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] =
542
            make_cell(particles, offset, size, h, rho, &partId, vel, press);
543
544
545
546
547
548

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

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

556
557
558
559
560
561
  /* 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 */
562
563
564
565
566
  ticks time = 0;
  for (size_t i = 0; i < runs; ++i) {

    const ticks tic = getticks();

567
568
569
570
571
572
    /* 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);

573
/* Do the density calculation */
574
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
575

576
577
578
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
    /* 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*/
607
    for (int j = 0; j < 27; ++j)
608
609
610
      runner_doself1_density(&runner, inner_cells[j]);

#endif
611

612
613
614
615
    /* 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 */
616
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
617

618
619
620
621
622
623
624
625
626
627
628
629
630
631
    /* 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);
632
633
#endif

634
635
636
    /* Finally, give a gentle kick */
    runner_do_kick(&runner, main_cell, 0);

637
638
639
640
    const ticks toc = getticks();
    time += toc - tic;

    /* Dump if necessary */
641
    if (i == 0) {
642
      sprintf(outputFileName, "swift_dopair_125_%s.dat",
643
              outputFileNameExtension);
644
      dump_particle_fields(outputFileName, main_cell, solution, 0);
645
646
647
648
649
650
    }
  }

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

651
652
653
  for (int j = 0; j < 125; ++j)
    reset_particles(cells[j], vel, press, size, rho);

654
655
656
657
658
659
660
661
  /* 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 */
662
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
663

664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
  /* 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);
            }
          }
        }
      }
    }
  }
693

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

697
#endif
698

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

702
/* Do the force calculation */
703
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
704

705
706
707
708
  /* 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++) {
709

710
711
712
713
714
715
716
717
718
        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);
719

720
#endif
721

722
723
  /* Finally, give a gentle kick */
  runner_do_kick(&runner, main_cell, 0);
724

725
  const ticks toc = getticks();
726

727
728
  /* Output timing */
  message("Brute force calculation took : %15lli ticks.", toc - tic);
729

730
731
  sprintf(outputFileName, "brute_force_125_%s.dat", outputFileNameExtension);
  dump_particle_fields(outputFileName, main_cell, solution, 0);
732
733

  /* Clean things to make the sanitizer happy ... */
734
  for (int i = 0; i < 125; ++i) clean_up(cells[i]);
735
  free(solution);
736
737
738

  return 0;
}