test125cells.c 26.3 KB
Newer Older
1

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

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

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

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

34
35
36
37
38
#if defined(WITH_VECTORIZATION)
#define DOSELF2 runner_doself2_force_vec
#define DOSELF2_NAME "runner_doself2_force_vec"
#define DOPAIR2_NAME "runner_dopair2_force"
#endif
39
40
41
42
43
44
45

#ifndef DOSELF2
#define DOSELF2 runner_doself2_force
#define DOSELF2_NAME "runner_doself2_density"
#define DOPAIR2_NAME "runner_dopair2_force"
#endif

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
98
99
100
101
102
103
104
105
106
107
108
109
110
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);
111
112
#elif defined(HOPKINS_PE_SPH)
  part->entropy = pressure / pow_gamma(density);
113
114
#elif defined(DEFAULT_SPH)
  part->u = pressure / (hydro_gamma_minus_one * density);
Matthieu Schaller's avatar
Matthieu Schaller committed
115
116
#elif defined(MINIMAL_SPH)
  part->u = pressure / (hydro_gamma_minus_one * density);
117
#elif defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
118
  part->primitives.P = pressure;
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
#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) {

147
  for (int i = 0; i < main_cell->count; ++i) {
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
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197

    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;
198
199
200

    solution[i].S_dt = 0.f;
    solution[i].u_dt = -(solution[i].P / solution[i].rho) * solution[i].div_v;
201
202
203
  }
}

204
205
206
void reset_particles(struct cell *c, struct hydro_space *hs,
                     enum velocity_field vel, enum pressure_field press,
                     float size, float density) {
207

208
  for (int i = 0; i < c->count; ++i) {
209
210
211
212
213

    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
214

215
216
    hydro_init_part(p, hs);

217
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
218
    float volume = p->conserved.mass / density;
219
220
221
222
223
#if defined(GIZMO_SPH)
    p->geometry.volume = volume;
#else
    p->cell.volume = volume;
#endif
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
224
225
226
227
228
229
230
231
    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.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 =
232
        p->primitives.P / hydro_gamma_minus_one * volume +
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
233
234
235
236
237
        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
238
239
240
  }
}

241
242
/**
 * @brief Constructs a cell and all of its particle in a valid state prior to
243
 * a SPH time-step.
244
245
246
247
248
 *
 * @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
249
 * separation.
250
251
 * @param density The density of the fluid.
 * @param partId The running counter of IDs.
252
253
 * @param pert The perturbation to apply to the particles in the cell in units
 *of the inter-particle separation.
254
255
 * @param vel The type of velocity field.
 * @param press The type of pressure field.
256
 */
257
struct cell *make_cell(size_t n, const double offset[3], double size, double h,
James Willis's avatar
James Willis committed
258
                       double density, long long *partId, double pert,
259
                       enum velocity_field vel, enum pressure_field press) {
260

261
262
263
264
265
266
  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,
267
                     count * sizeof(struct part)) != 0)
268
    error("couldn't allocate particles, no. of particles: %d", (int)count);
269
270
271
  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);
272
  bzero(cell->parts, count * sizeof(struct part));
273
  bzero(cell->xparts, count * sizeof(struct xpart));
274
275
276

  /* Construct the parts */
  struct part *part = cell->parts;
277
  struct xpart *xpart = cell->xparts;
278
279
280
  for (size_t x = 0; x < n; ++x) {
    for (size_t y = 0; y < n; ++y) {
      for (size_t z = 0; z < n; ++z) {
James Willis's avatar
James Willis committed
281
282
283
284
285
286
287
288
289
        part->x[0] =
            offset[0] +
            size * (x + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
        part->x[1] =
            offset[1] +
            size * (y + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
        part->x[2] =
            offset[2] +
            size * (z + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
290
        part->h = size * h / (float)n;
291

292
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
293
        part->conserved.mass = density * volume / count;
294
#else
295
        part->mass = density * volume / count;
296
#endif
297
298
299
300

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

301
        hydro_first_init_part(part, xpart);
302

303
        part->id = ++(*partId);
304
        part->time_bin = 1;
305

306
#if defined(GIZMO_SPH)
307
        part->geometry.volume = part->conserved.mass / density;
308
309
310
311
312
313
314
315
        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.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 =
316
            part->primitives.P / hydro_gamma_minus_one * volume +
317
318
319
320
321
322
            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

323
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
324
325
        part->ti_drift = 8;
        part->ti_kick = 8;
326
327
#endif

328
        ++part;
329
        ++xpart;
330
331
332
333
334
335
336
337
      }
    }
  }

  /* Cell properties */
  cell->split = 0;
  cell->h_max = h;
  cell->count = count;
338
  cell->gcount = 0;
339
  cell->dx_max_part = 0.;
340
  cell->dx_max_sort = 0.;
Matthieu Schaller's avatar
Matthieu Schaller committed
341
342
343
  cell->width[0] = size;
  cell->width[1] = size;
  cell->width[2] = size;
344
345
346
347
  cell->loc[0] = offset[0];
  cell->loc[1] = offset[1];
  cell->loc[2] = offset[2];

348
  cell->ti_old_part = 8;
Matthieu Schaller's avatar
Matthieu Schaller committed
349
350
  cell->ti_end_min = 8;
  cell->ti_end_max = 8;
351

352
  // shuffle_particles(cell->parts, cell->count);
353
354

  cell->sorted = 0;
Matthieu Schaller's avatar
Matthieu Schaller committed
355
  for (int k = 0; k < 13; k++) cell->sort[k] = NULL;
356
357
358
359
360
361

  return cell;
}

void clean_up(struct cell *ci) {
  free(ci->parts);
362
  free(ci->xparts);
363
  for (int k = 0; k < 13; k++)
Matthieu Schaller's avatar
Matthieu Schaller committed
364
    if (ci->sort[k] != NULL) free(ci->sort[k]);
365
366
367
368
369
370
371
  free(ci);
}

/**
 * @brief Dump all the particles to a file
 */
void dump_particle_fields(char *fileName, struct cell *main_cell,
372
                          struct solution_part *solution, int with_solution) {
373
374
375
376
  FILE *file = fopen(fileName, "w");

  /* Write header */
  fprintf(file,
James Willis's avatar
James Willis committed
377
          "# %4s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %13s %13s "
378
          "%13s %13s %13s %8s %8s\n",
379
380
381
          "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");
382
383
384
385

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

  /* Write main cell */
386
  for (int pid = 0; pid < main_cell->count; pid++) {
387
    fprintf(file,
388
389
            "%6llu %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f "
            "%8.5f "
390
            "%8.5f %8.5f %13e %13e %13e %13e %13e %8.5f %8.5f\n",
391
392
393
            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],
394
            main_cell->parts[pid].v[2], main_cell->parts[pid].h,
395
            hydro_get_density(&main_cell->parts[pid]),
396
#if defined(MINIMAL_SPH) || defined(SHADOWFAX_SPH)
Matthieu Schaller's avatar
Matthieu Schaller committed
397
            0.f,
Matthieu Schaller's avatar
Matthieu Schaller committed
398
#else
Matthieu Schaller's avatar
Matthieu Schaller committed
399
400
            main_cell->parts[pid].density.div_v,
#endif
401
402
403
404
            hydro_get_entropy(&main_cell->parts[pid]),
            hydro_get_internal_energy(&main_cell->parts[pid]),
            hydro_get_pressure(&main_cell->parts[pid]),
            hydro_get_soundspeed(&main_cell->parts[pid]),
405
406
            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,
407
#if defined(GADGET2_SPH)
408
409
            main_cell->parts[pid].force.v_sig, main_cell->parts[pid].entropy_dt,
            0.f
410
#elif defined(DEFAULT_SPH)
411
412
            main_cell->parts[pid].force.v_sig, 0.f,
            main_cell->parts[pid].force.u_dt
Matthieu Schaller's avatar
Matthieu Schaller committed
413
#elif defined(MINIMAL_SPH)
414
            main_cell->parts[pid].force.v_sig, 0.f, main_cell->parts[pid].u_dt
415
#else
416
            0.f, 0.f, 0.f
417
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
418
            );
419
420
  }

421
422
423
424
  if (with_solution) {

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

425
    for (int pid = 0; pid < main_cell->count; pid++) {
426
427
428
      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 "
429
              "%8.5f %8.5f %13f %13f %13f %13f %13f %8.5f %8.5f\n",
430
431
432
433
434
435
              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],
436
437
              solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt,
              solution[pid].u_dt);
438
439
    }
  }
440

441
442
443
444
445
  fclose(file);
}

/* Just a forward declaration... */
void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
446
447
void runner_dopair1_branch_density(struct runner *r, struct cell *ci,
                                   struct cell *cj);
448
void runner_doself1_density(struct runner *r, struct cell *ci);
449
450
void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself2_force(struct runner *r, struct cell *ci);
451
void runner_doself2_force_vec(struct runner *r, struct cell *ci);
452
453
454

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

James Willis's avatar
James Willis committed
456
  engine_pin();
457
  size_t runs = 0, particles = 0;
458
  double h = 1.23485, size = 1., rho = 2.5;
James Willis's avatar
James Willis committed
459
  double perturbation = 0.;
460
461
  char outputFileNameExtension[200] = "";
  char outputFileName[200] = "";
462
463
  enum velocity_field vel = velocity_zero;
  enum pressure_field press = pressure_const;
464
465
466
467
468

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

469
  /* Choke on FP-exceptions */
470
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
471

472
473
474
475
  /* Get some randomness going */
  srand(0);

  char c;
476
  while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:p:")) != -1) {
477
478
479
480
481
482
483
    switch (c) {
      case 'h':
        sscanf(optarg, "%lf", &h);
        break;
      case 's':
        sscanf(optarg, "%lf", &size);
        break;
484
      case 'n':
485
486
487
488
489
        sscanf(optarg, "%zu", &particles);
        break;
      case 'r':
        sscanf(optarg, "%zu", &runs);
        break;
James Willis's avatar
James Willis committed
490
491
492
      case 'd':
        sscanf(optarg, "%lf", &perturbation);
        break;
493
494
495
496
497
498
      case 'm':
        sscanf(optarg, "%lf", &rho);
        break;
      case 'f':
        strcpy(outputFileNameExtension, optarg);
        break;
499
500
501
502
503
504
      case 'v':
        sscanf(optarg, "%d", (int *)&vel);
        break;
      case 'p':
        sscanf(optarg, "%d", (int *)&press);
        break;
505
506
507
508
509
510
511
512
      case '?':
        error("Unknown option.");
        break;
    }
  }

  if (h < 0 || particles == 0 || runs == 0) {
    printf(
513
514
515
516
517
        "\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()"
518
519
520
521
        "\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"
James Willis's avatar
James Willis committed
522
        "\n-d pert            - Perturbation to apply to the particles [0,1["
523
        "\n-v type (0,1,2,3)  - Velocity field: (zero, constant, divergent, "
524
        "rotating)"
525
        "\n-p type (0,1,2)    - Pressure field: (constant, gradient divergent)"
526
527
528
529
530
531
        "\n-f fileName        - Part of the file name used to save the dumps\n",
        argv[0]);
    exit(1);
  }

  /* Help users... */
532
533
  message("DOSELF2 function called: %s", DOSELF2_NAME);
  message("DOPAIR2 function called: %s", DOPAIR2_NAME);
534
  message("Adiabatic index: ga = %f", hydro_gamma);
535
  message("Hydro implementation: %s", SPH_IMPLEMENTATION);
536
537
  message("Smoothing length: h = %f", h * size);
  message("Kernel:               %s", kernel_name);
538
  message("Neighbour target: N = %f", pow_dimension(h) * kernel_norm);
539
  message("Density target: rho = %f", rho);
540
541
542
543
544
545
546
547
548
  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");

549
550
  printf("\n");

551
552
553
554
#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
555

556
557
  /* Build the infrastructure */
  struct space space;
558
  space.periodic = 1;
559
560
561
  space.dim[0] = 5.;
  space.dim[1] = 5.;
  space.dim[2] = 5.;
562
  hydro_space_init(&space.hs, &space);
563

564
565
566
  struct phys_const prog_const;
  prog_const.const_newton_G = 1.f;

567
  struct hydro_props hp;
568
  hp.eta_neighbours = h;
569
  hp.h_tolerance = 1e0;
570
  hp.h_max = FLT_MAX;
571
  hp.max_smoothing_iterations = 1;
572
  hp.CFL_condition = 0.1;
573

574
  struct engine engine;
575
  bzero(&engine, sizeof(struct engine));
576
  engine.hydro_properties = &hp;
577
  engine.physical_constants = &prog_const;
578
579
  engine.s = &space;
  engine.time = 0.1f;
Matthieu Schaller's avatar
Matthieu Schaller committed
580
  engine.ti_current = 8;
581
  engine.max_active_bin = num_time_bins;
582
583
584
585
586

  struct runner runner;
  runner.e = &engine;

  /* Construct some cells */
587
588
  struct cell *cells[125];
  struct cell *inner_cells[27];
589
  struct cell *main_cell;
590
  int count = 0;
591
  static long long partId = 0;
592
593
594
595
596
597
598
599
  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 */
600
601
        cells[i * 25 + j * 5 + k] = make_cell(
            particles, offset, size, h, rho, &partId, perturbation, vel, press);
602
603
604
605
606
607

        /* 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++;
        }
608
609
610
611
612
      }
    }
  }

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

615
616
617
618
619
  /* 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);

620
621
622
  ticks timings[27];
  for (int i = 0; i < 27; i++) timings[i] = 0;

623
  /* Start the test */
624
  ticks time = 0;
625
  ticks self_force_time = 0;
626
  for (size_t n = 0; n < runs; ++n) {
627
628
629

    const ticks tic = getticks();

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

James Willis's avatar
James Willis committed
633
634
    /* Reset particles. */
    for (int i = 0; i < 125; ++i) {
635
636
      for (int pid = 0; pid < cells[i]->count; ++pid)
        hydro_init_part(&cells[i]->parts[pid], &space.hs);
James Willis's avatar
James Willis committed
637
    }
638

639
    /* First, sort stuff */
640
641
    for (int j = 0; j < 125; ++j)
      runner_do_sort(&runner, cells[j], 0x1FFF, 0, 0);
642

643
/* Do the density calculation */
644
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
645

Matthieu Schaller's avatar
Matthieu Schaller committed
646
/* Initialise the particle cache. */
647
648
649
650
651
652
653
#ifdef WITH_VECTORIZATION
    runner.ci_cache.count = 0;
    cache_init(&runner.ci_cache, 512);
    runner.cj_cache.count = 0;
    cache_init(&runner.cj_cache, 512);
#endif

654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
    /* 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];

676
                if (cj > ci) runner_dopair1_branch_density(&runner, ci, cj);
677
678
679
680
681
682
683
684
              }
            }
          }
        }
      }
    }

    /* And now the self-interaction for the central cells*/
685
    for (int j = 0; j < 27; ++j)
686
687
688
      runner_doself1_density(&runner, inner_cells[j]);

#endif
689

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

/* Do the force calculation */
694
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
695

696
    int ctr = 0;
697
698
699
700
701
702
703
    /* 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];

704
          if (main_cell != cj) {
James Willis's avatar
James Willis committed
705

706
707
708
709
710
711
712
            const ticks sub_tic = getticks();

            runner_dopair2_force(&runner, main_cell, cj);

            const ticks sub_toc = getticks();
            timings[ctr++] += sub_toc - sub_tic;
          }
713
714
715
716
        }
      }
    }

717
#ifdef WITH_VECTORIZATION
James Willis's avatar
James Willis committed
718
    /* Initialise the cache. */
719
720
    runner.ci_cache.count = 0;
    cache_init(&runner.ci_cache, 512);
721
#endif
722
723
724

    ticks self_tic = getticks();

725
    /* And now the self-interaction for the main cell */
726
    DOSELF2(&runner, main_cell);
James Willis's avatar
James Willis committed
727

James Willis's avatar
James Willis committed
728
    self_force_time += getticks() - self_tic;
729
    timings[26] += getticks() - self_tic;
730
731
#endif

732
    /* Finally, give a gentle kick */
733
    runner_do_end_force(&runner, main_cell, 0);
734
735
736
737
    const ticks toc = getticks();
    time += toc - tic;

    /* Dump if necessary */
738
    if (n == 0) {
739
      sprintf(outputFileName, "swift_dopair_125_%s.dat",
740
              outputFileNameExtension);
741
      dump_particle_fields(outputFileName, main_cell, solution, 0);
742
    }
James Willis's avatar
James Willis committed
743
744

    for (int i = 0; i < 125; ++i) {
745
746
      for (int pid = 0; pid < cells[i]->count; ++pid)
        hydro_init_part(&cells[i]->parts[pid], &space.hs);
James Willis's avatar
James Willis committed
747
    }
748
749
  }

750
  /* Output timing */
751
752
753
754
755
756
757
758
759
760
761
762
763
  ticks corner_time = timings[0] + timings[2] + timings[6] + timings[8] +
                      timings[17] + timings[19] + timings[23] + timings[25];

  ticks edge_time = timings[1] + timings[3] + timings[5] + timings[7] +
                    timings[9] + timings[11] + timings[14] + timings[16] +
                    timings[18] + timings[20] + timings[22] + timings[24];

  ticks face_time = timings[4] + timings[10] + timings[12] + timings[13] +
                    timings[15] + timings[21];

  message("Corner calculations took       : %15lli ticks.", corner_time / runs);
  message("Edge calculations took         : %15lli ticks.", edge_time / runs);
  message("Face calculations took         : %15lli ticks.", face_time / runs);
James Willis's avatar
James Willis committed
764
  message("Self calculations took         : %15lli ticks.",
James Willis's avatar
James Willis committed
765
          self_force_time / runs);
766
  message("SWIFT calculation took         : %15lli ticks.", time / runs);
767

768
  for (int j = 0; j < 125; ++j)
769
    reset_particles(cells[j], &space.hs, vel, press, size, rho);
770

771
772
773
774
  /* NOW BRUTE-FORCE CALCULATION */

  const ticks tic = getticks();

775
776
/* Kick the central cell */
// runner_do_kick1(&runner, main_cell, 0);
777

778
779
/* And drift it */
// runner_do_drift_particles(&runner, main_cell, 0);
780

781
782
783
/* Initialise the particles */
// for (int j = 0; j < 125; ++j) runner_do_drift_particles(&runner, cells[j],
// 0);
784
785

/* Do the density calculation */
786
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
787

788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
  /* 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);
            }
          }
        }
      }
    }
  }
817

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

821
#endif
822

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

826
/* Do the force calculation */
827
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
828

829
830
831
832
  /* 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++) {
833

834
835
836
837
838
839
840
841
842
        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);
843

844
#endif
845

846
  /* Finally, give a gentle kick */
847
  runner_do_end_force(&runner, main_cell, 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
848
  // runner_do_kick2(&runner, main_cell, 0);
849

850
  const ticks toc = getticks();
851

852
853
  /* Output timing */
  message("Brute force calculation took : %15lli ticks.", toc - tic);
854

855
856
  sprintf(outputFileName, "brute_force_125_%s.dat", outputFileNameExtension);
  dump_particle_fields(outputFileName, main_cell, solution, 0);
857
858

  /* Clean things to make the sanitizer happy ... */
859
  for (int i = 0; i < 125; ++i) clean_up(cells[i]);
860
  free(solution);
861
862
863

  return 0;
}