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
  float h_max = 0.f;
James Willis's avatar
James Willis committed
276

277
278
  /* Construct the parts */
  struct part *part = cell->parts;
279
  struct xpart *xpart = cell->xparts;
280
281
282
  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
283
284
285
286
287
288
289
290
291
        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;
292
        part->h = size * h / (float)n;
293
        h_max = fmax(h_max, part->h);
294

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

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

304
        hydro_first_init_part(part, xpart);
305

306
        part->id = ++(*partId);
307
        part->time_bin = 1;
308

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

326
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
327
328
        part->ti_drift = 8;
        part->ti_kick = 8;
329
330
#endif

331
        ++part;
332
        ++xpart;
333
334
335
336
337
338
      }
    }
  }

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

351
  cell->ti_old_part = 8;
Matthieu Schaller's avatar
Matthieu Schaller committed
352
353
  cell->ti_end_min = 8;
  cell->ti_end_max = 8;
354

355
  // shuffle_particles(cell->parts, cell->count);
356
357

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

  return cell;
}

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

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

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

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

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

424
425
426
427
  if (with_solution) {

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

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

444
445
446
447
448
  fclose(file);
}

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

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

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

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

472
  /* Choke on FP-exceptions */
473
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
474

475
476
477
478
  /* Get some randomness going */
  srand(0);

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

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

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

552
553
  printf("\n");

554
555
556
557
#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
558

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

567
568
569
  struct phys_const prog_const;
  prog_const.const_newton_G = 1.f;

570
  struct hydro_props hp;
571
  hp.eta_neighbours = h;
572
  hp.h_tolerance = 1e0;
573
  hp.h_max = FLT_MAX;
574
  hp.max_smoothing_iterations = 1;
575
  hp.CFL_condition = 0.1;
576

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

  struct runner runner;
  runner.e = &engine;

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

        /* 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++;
        }
611
612
613
614
615
      }
    }
  }

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

618
619
620
621
622
  /* 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);

623
624
625
  ticks timings[27];
  for (int i = 0; i < 27; i++) timings[i] = 0;

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

    const ticks tic = getticks();

633
    /* Initialise the particles */
634
    for (int j = 0; j < 125; ++j) runner_do_drift_part(&runner, cells[j], 0);
635

James Willis's avatar
James Willis committed
636
637
    /* Reset particles. */
    for (int i = 0; i < 125; ++i) {
638
639
      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
640
    }
641

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

646
/* Do the density calculation */
647
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
648

Matthieu Schaller's avatar
Matthieu Schaller committed
649
/* Initialise the particle cache. */
650
651
652
653
654
655
656
#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

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

679
                if (cj > ci) runner_dopair1_branch_density(&runner, ci, cj);
680
681
682
683
684
685
686
687
              }
            }
          }
        }
      }
    }

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

#endif
692

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

/* Do the force calculation */
697
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
698

699
    int ctr = 0;
700
701
702
703
704
705
706
    /* 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];

707
          if (main_cell != cj) {
James Willis's avatar
James Willis committed
708

709
710
711
712
713
714
715
            const ticks sub_tic = getticks();

            runner_dopair2_force(&runner, main_cell, cj);

            const ticks sub_toc = getticks();
            timings[ctr++] += sub_toc - sub_tic;
          }
716
717
718
719
        }
      }
    }

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

    ticks self_tic = getticks();

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

James Willis's avatar
James Willis committed
731
    self_force_time += getticks() - self_tic;
732
    timings[26] += getticks() - self_tic;
733
734
#endif

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

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

    for (int i = 0; i < 125; ++i) {
748
749
      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
750
    }
751
752
  }

753
  /* Output timing */
754
755
756
757
758
759
760
761
762
763
764
765
766
  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
767
  message("Self calculations took         : %15lli ticks.",
James Willis's avatar
James Willis committed
768
          self_force_time / runs);
769
  message("SWIFT calculation took         : %15lli ticks.", time / runs);
770

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

774
775
776
777
  /* NOW BRUTE-FORCE CALCULATION */

  const ticks tic = getticks();

778
779
/* Kick the central cell */
// runner_do_kick1(&runner, main_cell, 0);
780

781
782
/* And drift it */
// runner_do_drift_particles(&runner, main_cell, 0);
783

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

/* Do the density calculation */
789
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
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
817
818
819
  /* 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);
            }
          }
        }
      }
    }
  }
820

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

824
#endif
825

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

829
/* Do the force calculation */
830
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
831

832
833
834
835
  /* 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++) {
836

837
838
839
840
841
842
843
844
845
        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);
846

847
#endif
848

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

853
  const ticks toc = getticks();
854

855
856
  /* Output timing */
  message("Brute force calculation took : %15lli ticks.", toc - tic);
857

858
859
  sprintf(outputFileName, "brute_force_125_%s.dat", outputFileNameExtension);
  dump_particle_fields(outputFileName, main_cell, solution, 0);
860
861

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

  return 0;
}