test125cells.c 20.8 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
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
#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;
181
182
183

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

187
188
189
190
191
192
193
194
195
196
197
198
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);
  }
}

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

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

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

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

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

250
251
252
        xpart->v_full[0] = part->v[0];
        xpart->v_full[1] = part->v[1];
        xpart->v_full[2] = part->v[2];
253
        hydro_first_init_part(part, xpart);
254
        ++part;
255
        ++xpart;
256
257
258
259
260
261
262
263
      }
    }
  }

  /* Cell properties */
  cell->split = 0;
  cell->h_max = h;
  cell->count = count;
264
  cell->gcount = 0;
265
  cell->dx_max = 0.;
Matthieu Schaller's avatar
Matthieu Schaller committed
266
267
268
  cell->width[0] = size;
  cell->width[1] = size;
  cell->width[2] = size;
269
270
271
272
273
274
275
  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;

276
  // shuffle_particles(cell->parts, cell->count);
277
278
279
280
281
282
283
284
285
286

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

  return cell;
}

void clean_up(struct cell *ci) {
  free(ci->parts);
287
  free(ci->xparts);
288
289
290
291
292
293
294
295
  free(ci->sort);
  free(ci);
}

/**
 * @brief Dump all the particles to a file
 */
void dump_particle_fields(char *fileName, struct cell *main_cell,
296
                          struct solution_part *solution, int with_solution) {
297
298
299
300
  FILE *file = fopen(fileName, "w");

  /* Write header */
  fprintf(file,
301
302
303
304
305
          "# %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");
306
307
308
309
310
311

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

  /* Write main cell */
  for (size_t pid = 0; pid < main_cell->count; pid++) {
    fprintf(file,
312
313
314
            "%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",
315
316
317
            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],
318
            main_cell->parts[pid].v[2], main_cell->parts[pid].h,
Matthieu Schaller's avatar
Matthieu Schaller committed
319
320
            main_cell->parts[pid].rho,
#ifdef MINIMAL_SPH
Matthieu Schaller's avatar
Matthieu Schaller committed
321
            0.f,
Matthieu Schaller's avatar
Matthieu Schaller committed
322
#else
Matthieu Schaller's avatar
Matthieu Schaller committed
323
324
            main_cell->parts[pid].density.div_v,
#endif
325
326
327
328
329
330
            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,
331
#if defined(GADGET2_SPH)
332
333
            main_cell->parts[pid].force.v_sig, main_cell->parts[pid].entropy_dt,
            0.f
334
#elif defined(DEFAULT_SPH)
335
336
            main_cell->parts[pid].force.v_sig, 0.f,
            main_cell->parts[pid].force.u_dt
Matthieu Schaller's avatar
Matthieu Schaller committed
337
#elif defined(MINIMAL_SPH)
338
            main_cell->parts[pid].force.v_sig, 0.f, main_cell->parts[pid].u_dt
339
#else
340
            0.f, 0.f, 0.f
341
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
342
            );
343
344
  }

345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
  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],
360
361
              solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt,
              solution[pid].u_dt);
362
363
    }
  }
364

365
366
367
368
369
370
  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);
371
372
void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself2_force(struct runner *r, struct cell *ci);
373
374
375

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

377
  size_t runs = 0, particles = 0;
378
  double h = 1.23485, size = 1., rho = 2.5;
379
380
  char outputFileNameExtension[200] = "";
  char outputFileName[200] = "";
381
382
  enum velocity_field vel = velocity_zero;
  enum pressure_field press = pressure_const;
383
384
385
386
387

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

388
  /* Choke on FP-exceptions */
389
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
390

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

  char c;
395
  while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:p:")) != -1) {
396
397
398
399
400
401
402
    switch (c) {
      case 'h':
        sscanf(optarg, "%lf", &h);
        break;
      case 's':
        sscanf(optarg, "%lf", &size);
        break;
403
      case 'n':
404
405
406
407
408
409
410
411
412
413
414
        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;
415
416
417
418
419
420
      case 'v':
        sscanf(optarg, "%d", (int *)&vel);
        break;
      case 'p':
        sscanf(optarg, "%d", (int *)&press);
        break;
421
422
423
424
425
426
427
428
      case '?':
        error("Unknown option.");
        break;
    }
  }

  if (h < 0 || particles == 0 || runs == 0) {
    printf(
429
430
431
432
433
        "\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()"
434
435
436
437
        "\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"
438
        "\n-v type (0,1,2,3)  - Velocity field: (zero, constant, divergent, "
439
        "rotating)"
440
        "\n-p type (0,1,2)    - Pressure field: (constant, gradient divergent)"
441
442
443
444
445
446
        "\n-f fileName        - Part of the file name used to save the dumps\n",
        argv[0]);
    exit(1);
  }

  /* Help users... */
447
  message("Adiabatic index: ga = %f", hydro_gamma);
448
  message("Hydro implementation: %s", SPH_IMPLEMENTATION);
449
450
  message("Smoothing length: h = %f", h * size);
  message("Kernel:               %s", kernel_name);
451
  message("Neighbour target: N = %f", pow_dimension(h) * kernel_norm);
452
  message("Density target: rho = %f", rho);
453
454
455
456
457
458
459
460
461
  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");

462
463
  printf("\n");

464
465
466
467
#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
468

469
470
471
472
473
  /* Build the infrastructure */
  struct space space;
  space.periodic = 0;
  space.h_max = h;

474
475
476
  struct phys_const prog_const;
  prog_const.const_newton_G = 1.f;

477
  struct hydro_props hp;
478
479
480
  hp.target_neighbours = h * h * h * kernel_norm;
  hp.delta_neighbours = 1.;
  hp.max_smoothing_iterations = 1;
481
  hp.CFL_condition = 0.1;
482

483
  struct engine engine;
484
  engine.hydro_properties = &hp;
485
  engine.physical_constants = &prog_const;
486
487
488
489
490
491
492
493
  engine.s = &space;
  engine.time = 0.1f;
  engine.ti_current = 1;

  struct runner runner;
  runner.e = &engine;

  /* Construct some cells */
494
495
  struct cell *cells[125];
  struct cell *inner_cells[27];
496
  struct cell *main_cell;
497
  int count = 0;
498
  static long long partId = 0;
499
500
501
502
503
504
505
506
507
  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] =
508
            make_cell(particles, offset, size, h, rho, &partId, vel, press);
509
510
511
512
513
514

        /* 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++;
        }
515
516
517
518
519
      }
    }
  }

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

522
523
524
525
526
527
  /* 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 */
528
529
530
531
532
  ticks time = 0;
  for (size_t i = 0; i < runs; ++i) {

    const ticks tic = getticks();

533
534
535
536
537
538
    /* 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);

539
/* Do the density calculation */
540
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
541

542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
    /* 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*/
573
    for (int j = 0; j < 27; ++j)
574
575
576
      runner_doself1_density(&runner, inner_cells[j]);

#endif
577

578
579
580
581
    /* 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 */
582
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
583

584
585
586
587
588
589
590
591
592
593
594
595
596
597
    /* 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);
598
599
#endif

600
601
602
    /* Finally, give a gentle kick */
    runner_do_kick(&runner, main_cell, 0);

603
604
605
606
    const ticks toc = getticks();
    time += toc - tic;

    /* Dump if necessary */
607
    if (i == 0) {
608
      sprintf(outputFileName, "swift_dopair_125_%s.dat",
609
              outputFileNameExtension);
610
      dump_particle_fields(outputFileName, main_cell, solution, 0);
611
612
613
614
615
616
    }
  }

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

617
618
619
  for (int j = 0; j < 125; ++j)
    reset_particles(cells[j], vel, press, size, rho);

620
621
622
623
624
625
626
627
  /* 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 */
628
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
629

630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
  /* 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);
            }
          }
        }
      }
    }
  }
659

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

663
#endif
664

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

668
/* Do the force calculation */
669
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
670

671
672
673
674
  /* 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++) {
675

676
677
678
679
680
681
682
683
684
        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);
685

686
#endif
687

688
689
  /* Finally, give a gentle kick */
  runner_do_kick(&runner, main_cell, 0);
690

691
  const ticks toc = getticks();
692

693
694
  /* Output timing */
  message("Brute force calculation took : %15lli ticks.", toc - tic);
695

696
697
  sprintf(outputFileName, "brute_force_125_%s.dat", outputFileNameExtension);
  dump_particle_fields(outputFileName, main_cell, solution, 0);
698
699

  /* Clean things to make the sanitizer happy ... */
700
  for (int i = 0; i < 125; ++i) clean_up(cells[i]);
701
  free(solution);
702
703
704

  return 0;
}