test125cells.c 30.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
#if defined(WITH_VECTORIZATION)
#define DOSELF2_NAME "runner_doself2_force_vec"
35
#define DOPAIR2_NAME "runner_dopair2_force_vec"
36
#endif
37

38
#ifndef DOSELF2_NAME
39 40 41 42
#define DOSELF2_NAME "runner_doself2_density"
#define DOPAIR2_NAME "runner_dopair2_force"
#endif

43
#define NODE_ID 0
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 98 99 100 101 102 103 104 105 106 107 108 109
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);
110 111
#elif defined(HOPKINS_PE_SPH)
  part->entropy = pressure / pow_gamma(density);
Josh Borrow's avatar
Josh Borrow committed
112
#elif defined(PHANTOM_SPH)
113
  part->u = pressure / (hydro_gamma_minus_one * density);
114 115
#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) ||           \
    defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH) || \
Josh Borrow's avatar
Josh Borrow committed
116
    defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
117
  part->u = pressure / (hydro_gamma_minus_one * density);
118
#elif defined(PLANETARY_SPH)
119
  part->u = pressure / (hydro_gamma_minus_one * density);
120
#elif defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
121
  part->primitives.P = pressure;
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
#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) {

150
  for (int i = 0; i < main_cell->hydro.count; ++i) {
151

152
    solution[i].id = main_cell->hydro.parts[i].id;
153

154 155 156
    solution[i].x[0] = main_cell->hydro.parts[i].x[0];
    solution[i].x[1] = main_cell->hydro.parts[i].x[1];
    solution[i].x[2] = main_cell->hydro.parts[i].x[2];
157

158 159 160
    solution[i].v[0] = main_cell->hydro.parts[i].v[0];
    solution[i].v[1] = main_cell->hydro.parts[i].v[1];
    solution[i].v[2] = main_cell->hydro.parts[i].v[2];
161

162
    solution[i].h = main_cell->hydro.parts[i].h;
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 198 199 200

    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;
201 202 203

    solution[i].S_dt = 0.f;
    solution[i].u_dt = -(solution[i].P / solution[i].rho) * solution[i].div_v;
204 205 206
  }
}

207 208 209
void reset_particles(struct cell *c, struct hydro_space *hs,
                     enum velocity_field vel, enum pressure_field press,
                     float size, float density) {
210

211
  for (int i = 0; i < c->hydro.count; ++i) {
212

213
    struct part *p = &c->hydro.parts[i];
214 215 216

    set_velocity(p, vel, size);
    set_energy_state(p, press, size, density);
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
217

218 219
    hydro_init_part(p, hs);

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

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

265 266
  const size_t count = n * n * n;
  const double volume = size * size * size;
267
  struct cell *cell = (struct cell *)malloc(sizeof(struct cell));
268 269
  bzero(cell, sizeof(struct cell));

270
  if (posix_memalign((void **)&cell->hydro.parts, part_align,
271
                     count * sizeof(struct part)) != 0)
272
    error("couldn't allocate particles, no. of particles: %d", (int)count);
273
  if (posix_memalign((void **)&cell->hydro.xparts, xpart_align,
274 275
                     count * sizeof(struct xpart)) != 0)
    error("couldn't allocate particles, no. of x-particles: %d", (int)count);
276 277
  bzero(cell->hydro.parts, count * sizeof(struct part));
  bzero(cell->hydro.xparts, count * sizeof(struct xpart));
278

279
  float h_max = 0.f;
James Willis's avatar
James Willis committed
280

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

299
#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
300
        part->conserved.mass = density * volume / count;
301
#else
302
        part->mass = density * volume / count;
303
#endif
304 305 306 307

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

308
        hydro_first_init_part(part, xpart);
309

310
        part->id = ++(*partId);
311
        part->time_bin = 1;
312

313
#if defined(GIZMO_MFV_SPH)
314
        part->geometry.volume = part->conserved.mass / density;
315 316 317 318 319 320 321 322
        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 =
323
            part->primitives.P / hydro_gamma_minus_one * volume +
324 325 326 327
            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]) /
328 329 330
                part->conserved.mass;
#endif

331
#ifdef SWIFT_DEBUG_CHECKS
332 333
        part->ti_drift = 8;
        part->ti_kick = 8;
334 335
#endif

336
        ++part;
337
        ++xpart;
338 339 340 341 342 343
      }
    }
  }

  /* Cell properties */
  cell->split = 0;
344 345
  cell->hydro.h_max = h_max;
  cell->hydro.count = count;
346 347
  cell->grav.count = 0;
  cell->hydro.dx_max_part = 0.;
348
  cell->hydro.dx_max_sort = 0.;
Matthieu Schaller's avatar
Matthieu Schaller committed
349 350 351
  cell->width[0] = size;
  cell->width[1] = size;
  cell->width[2] = size;
352 353 354 355
  cell->loc[0] = offset[0];
  cell->loc[1] = offset[1];
  cell->loc[2] = offset[2];

356
  cell->hydro.ti_old_part = 8;
357 358
  cell->hydro.ti_end_min = 8;
  cell->hydro.ti_end_max = 8;
359
  cell->nodeID = NODE_ID;
360

361
  // shuffle_particles(cell->hydro.parts, cell->hydro.count);
362

363
  cell->hydro.sorted = 0;
364
  cell->hydro.sort = NULL;
365 366 367 368 369

  return cell;
}

void clean_up(struct cell *ci) {
370 371
  free(ci->hydro.parts);
  free(ci->hydro.xparts);
372
  free(ci->hydro.sort);
373 374 375 376 377 378 379
  free(ci);
}

/**
 * @brief Dump all the particles to a file
 */
void dump_particle_fields(char *fileName, struct cell *main_cell,
380
                          struct solution_part *solution, int with_solution) {
381 382 383 384
  FILE *file = fopen(fileName, "w");

  /* Write header */
  fprintf(file,
James Willis's avatar
James Willis committed
385
          "# %4s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %13s %13s "
386
          "%13s %13s %13s %8s %8s\n",
387 388 389
          "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");
390 391 392 393

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

  /* Write main cell */
394
  for (int pid = 0; pid < main_cell->hydro.count; pid++) {
395
    fprintf(file,
396 397
            "%6llu %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f "
            "%8.5f "
398
            "%8.5f %8.5f %13e %13e %13e %13e %13e %8.5f %8.5f\n",
399 400 401 402 403
            main_cell->hydro.parts[pid].id, main_cell->hydro.parts[pid].x[0],
            main_cell->hydro.parts[pid].x[1], main_cell->hydro.parts[pid].x[2],
            main_cell->hydro.parts[pid].v[0], main_cell->hydro.parts[pid].v[1],
            main_cell->hydro.parts[pid].v[2], main_cell->hydro.parts[pid].h,
            hydro_get_comoving_density(&main_cell->hydro.parts[pid]),
404 405 406
#if defined(MINIMAL_SPH) || defined(PLANETARY_SPH) ||   \
    defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH) || \
    defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN)
Matthieu Schaller's avatar
Matthieu Schaller committed
407
            0.f,
Josh Borrow's avatar
Josh Borrow committed
408
#elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
Josh Borrow's avatar
Josh Borrow committed
409
            main_cell->hydro.parts[pid].viscosity.div_v,
410
#else
411
            main_cell->hydro.parts[pid].density.div_v,
Matthieu Schaller's avatar
Matthieu Schaller committed
412
#endif
413
            hydro_get_drifted_comoving_entropy(&main_cell->hydro.parts[pid]),
Loic Hausammann's avatar
Loic Hausammann committed
414 415
            hydro_get_drifted_comoving_internal_energy(
                &main_cell->hydro.parts[pid]),
416 417
            hydro_get_comoving_pressure(&main_cell->hydro.parts[pid]),
            hydro_get_comoving_soundspeed(&main_cell->hydro.parts[pid]),
Loic Hausammann's avatar
Loic Hausammann committed
418 419 420 421
            main_cell->hydro.parts[pid].a_hydro[0],
            main_cell->hydro.parts[pid].a_hydro[1],
            main_cell->hydro.parts[pid].a_hydro[2],
            main_cell->hydro.parts[pid].force.h_dt,
422
#if defined(GADGET2_SPH)
423 424
            main_cell->hydro.parts[pid].force.v_sig,
            main_cell->hydro.parts[pid].entropy_dt, 0.f
Josh Borrow's avatar
Josh Borrow committed
425 426
#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
    defined(HOPKINS_PU_SPH_MONAGHAN)
427 428
            main_cell->hydro.parts[pid].force.v_sig, 0.f,
            main_cell->hydro.parts[pid].u_dt
Josh Borrow's avatar
Josh Borrow committed
429
#elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
Josh Borrow's avatar
Josh Borrow committed
430 431
            main_cell->hydro.parts[pid].viscosity.v_sig, 0.f,
            main_cell->hydro.parts[pid].u_dt
432
#else
433
            0.f, 0.f, 0.f
434
#endif
435
    );
436 437
  }

438 439 440 441
  if (with_solution) {

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

442
    for (int pid = 0; pid < main_cell->hydro.count; pid++) {
443 444 445
      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 "
446
              "%8.5f %8.5f %13f %13f %13f %13f %13f %8.5f %8.5f\n",
447 448 449 450 451 452
              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],
453 454
              solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt,
              solution[pid].u_dt);
455 456
    }
  }
457

458 459 460 461
  fclose(file);
}

/* Just a forward declaration... */
462 463
void runner_dopair1_branch_density(struct runner *r, struct cell *ci,
                                   struct cell *cj);
464
void runner_doself1_branch_density(struct runner *r, struct cell *ci);
Josh Borrow's avatar
Josh Borrow committed
465 466 467 468 469
#ifdef EXTRA_HYDRO_LOOP
void runner_dopair1_branch_gradient(struct runner *r, struct cell *ci,
                                    struct cell *cj);
void runner_doself1_branch_gradient(struct runner *r, struct cell *ci);
#endif /* EXTRA_HYDRO LOOP */
470 471
void runner_dopair2_branch_force(struct runner *r, struct cell *ci,
                                 struct cell *cj);
472
void runner_doself2_branch_force(struct runner *r, struct cell *ci);
Josh Borrow's avatar
Josh Borrow committed
473 474
void runner_doself2_force(struct runner *r, struct cell *ci);
void runner_doself2_force_vec(struct runner *r, struct cell *ci);
475 476 477

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

479
#ifdef HAVE_SETAFFINITY
James Willis's avatar
James Willis committed
480
  engine_pin();
481 482
#endif

483
  size_t runs = 0, particles = 0;
484
  double h = 1.23485, size = 1., rho = 2.5;
James Willis's avatar
James Willis committed
485
  double perturbation = 0.;
486
  char outputFileNameExtension[100] = "";
487
  char outputFileName[200] = "";
488 489
  enum velocity_field vel = velocity_zero;
  enum pressure_field press = pressure_const;
490 491 492 493 494

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

495 496
/* Choke on FP-exceptions */
#ifdef HAVE_FE_ENABLE_EXCEPT
497
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
498
#endif
499

500 501 502 503
  /* Get some randomness going */
  srand(0);

  char c;
504
  while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:p:")) != -1) {
505 506 507 508 509 510 511
    switch (c) {
      case 'h':
        sscanf(optarg, "%lf", &h);
        break;
      case 's':
        sscanf(optarg, "%lf", &size);
        break;
512
      case 'n':
513 514 515 516 517
        sscanf(optarg, "%zu", &particles);
        break;
      case 'r':
        sscanf(optarg, "%zu", &runs);
        break;
James Willis's avatar
James Willis committed
518 519 520
      case 'd':
        sscanf(optarg, "%lf", &perturbation);
        break;
521 522 523 524 525 526
      case 'm':
        sscanf(optarg, "%lf", &rho);
        break;
      case 'f':
        strcpy(outputFileNameExtension, optarg);
        break;
527 528 529 530 531 532
      case 'v':
        sscanf(optarg, "%d", (int *)&vel);
        break;
      case 'p':
        sscanf(optarg, "%d", (int *)&press);
        break;
533 534 535 536 537 538 539 540
      case '?':
        error("Unknown option.");
        break;
    }
  }

  if (h < 0 || particles == 0 || runs == 0) {
    printf(
541 542 543 544 545
        "\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()"
546 547 548 549
        "\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
550
        "\n-d pert            - Perturbation to apply to the particles [0,1["
551
        "\n-v type (0,1,2,3)  - Velocity field: (zero, constant, divergent, "
552
        "rotating)"
553
        "\n-p type (0,1,2)    - Pressure field: (constant, gradient divergent)"
554 555 556 557 558 559
        "\n-f fileName        - Part of the file name used to save the dumps\n",
        argv[0]);
    exit(1);
  }

  /* Help users... */
560 561
  message("DOSELF2 function called: %s", DOSELF2_NAME);
  message("DOPAIR2 function called: %s", DOPAIR2_NAME);
562
  message("Adiabatic index: ga = %f", hydro_gamma);
563
  message("Hydro implementation: %s", SPH_IMPLEMENTATION);
564 565
  message("Smoothing length: h = %f", h * size);
  message("Kernel:               %s", kernel_name);
566
  message("Neighbour target: N = %f", pow_dimension(h) * kernel_norm);
567
  message("Density target: rho = %f", rho);
568 569 570 571 572 573 574 575 576
  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");

577 578
  printf("\n");

579 580 581 582
#if !defined(HYDRO_DIMENSION_3D)
  message("test125cells only useful in 3D. Change parameters in const.h !");
  return 1;
#endif
583

584 585
  /* Build the infrastructure */
  struct space space;
586
  space.periodic = 1;
587 588 589
  space.dim[0] = 5.;
  space.dim[1] = 5.;
  space.dim[2] = 5.;
590
  hydro_space_init(&space.hs, &space);
591

592 593 594
  struct phys_const prog_const;
  prog_const.const_newton_G = 1.f;

595
  struct hydro_props hp;
596
  hydro_props_init_no_hydro(&hp);
597
  hp.eta_neighbours = h;
598
  hp.h_tolerance = 1e0;
599
  hp.h_max = FLT_MAX;
600
  hp.h_min = 0.f;
601
  hp.h_min_ratio = 0.f;
Josh Borrow's avatar
Josh Borrow committed
602
  hp.max_smoothing_iterations = 10;
603
  hp.CFL_condition = 0.1;
604

605
  struct engine engine;
606
  bzero(&engine, sizeof(struct engine));
607
  engine.hydro_properties = &hp;
608
  engine.physical_constants = &prog_const;
609 610
  engine.s = &space;
  engine.time = 0.1f;
611
  engine.ti_current = 8;
612
  engine.max_active_bin = num_time_bins;
613
  engine.nodeID = NODE_ID;
614

615 616 617 618
  struct cosmology cosmo;
  cosmology_init_no_cosmo(&cosmo);
  engine.cosmology = &cosmo;

619 620 621 622
  struct runner runner;
  runner.e = &engine;

  /* Construct some cells */
623 624
  struct cell *cells[125];
  struct cell *inner_cells[27];
625
  struct cell *main_cell;
626
  int count = 0;
627
  static long long partId = 0;
628 629 630 631 632 633 634 635
  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 */
636 637
        cells[i * 25 + j * 5 + k] = make_cell(
            particles, offset, size, h, rho, &partId, perturbation, vel, press);
638 639 640 641 642 643

        /* 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++;
        }
644 645 646 647 648
      }
    }
  }

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

651
  /* Construct the real solution */
652
  struct solution_part *solution = (struct solution_part *)malloc(
653
      main_cell->hydro.count * sizeof(struct solution_part));
654 655
  get_solution(main_cell, solution, rho, vel, press, size);

656 657 658
  ticks timings[27];
  for (int i = 0; i < 27; i++) timings[i] = 0;

659
  /* Start the test */
660
  ticks time = 0;
661
  for (size_t n = 0; n < runs; ++n) {
662 663 664

    const ticks tic = getticks();

665
    /* Initialise the particles */
666
    for (int j = 0; j < 125; ++j) runner_do_drift_part(&runner, cells[j], 0);
667

James Willis's avatar
James Willis committed
668 669
    /* Reset particles. */
    for (int i = 0; i < 125; ++i) {
670 671
      for (int pid = 0; pid < cells[i]->hydro.count; ++pid)
        hydro_init_part(&cells[i]->hydro.parts[pid], &space.hs);
James Willis's avatar
James Willis committed
672
    }
673

674
    /* First, sort stuff */
675
    for (int j = 0; j < 125; ++j)
Loic Hausammann's avatar
Loic Hausammann committed
676
      runner_do_hydro_sort(&runner, cells[j], 0x1FFF, 0, 0);
677

678
      /* Do the density calculation */
679

Matthieu Schaller's avatar
Matthieu Schaller committed
680
/* Initialise the particle cache. */
681 682 683
#ifdef WITH_VECTORIZATION
    runner.ci_cache.count = 0;
    runner.cj_cache.count = 0;
684
    cache_init(&runner.ci_cache, 512);
685 686 687
    cache_init(&runner.cj_cache, 512);
#endif

Josh Borrow's avatar
Josh Borrow committed
688
    /* Run all the  (only once !)*/
689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709
    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];

710
                if (cj > ci) runner_dopair1_branch_density(&runner, ci, cj);
711 712 713 714 715 716 717 718
              }
            }
          }
        }
      }
    }

    /* And now the self-interaction for the central cells*/
719
    for (int j = 0; j < 27; ++j)
720
      runner_doself1_branch_density(&runner, inner_cells[j]);
721 722

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

Josh Borrow's avatar
Josh Borrow committed
725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769
#ifdef EXTRA_HYDRO_LOOP
    /* We need to do the gradient loop and the extra ghost! */
    message(
        "Extra hydro loop detected, running gradient loop in test125cells.");

    /* 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_branch_gradient(&runner, ci, cj);
              }
            }
          }
        }
      }
    }

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

    /* Extra ghost to finish everything on the central cells */
    for (int j = 0; j < 27; ++j)
      runner_do_extra_ghost(&runner, inner_cells[j], 0);

#endif /* EXTRA_HYDRO_LOOP */

770
      /* Do the force calculation */
771

772 773
#ifdef WITH_VECTORIZATION
    /* Initialise the cache. */
774 775
    cache_clean(&runner.ci_cache);
    cache_clean(&runner.cj_cache);
776 777 778 779
    cache_init(&runner.ci_cache, 512);
    cache_init(&runner.cj_cache, 512);
#endif

780
    int ctr = 0;
781 782 783 784 785 786 787
    /* 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];

788
          if (main_cell != cj) {
James Willis's avatar
James Willis committed
789

790 791
            const ticks sub_tic = getticks();

792
            runner_dopair2_branch_force(&runner, main_cell, cj);
793

794
            timings[ctr++] += getticks() - sub_tic;
795
          }
796 797 798 799
        }
      }
    }

800 801
    ticks self_tic = getticks();

802
    /* And now the self-interaction for the main cell */
803
    runner_doself2_branch_force(&runner, main_cell);
James Willis's avatar
James Willis committed
804

805
    timings[26] += getticks() - self_tic;
806

807
    /* Finally, give a gentle kick */
808
    runner_do_end_hydro_force(&runner, main_cell, 0);
809 810 811 812
    const ticks toc = getticks();
    time += toc - tic;

    /* Dump if necessary */
813
    if (n == 0) {
lhausamm's avatar
lhausamm committed
814
      sprintf(outputFileName, "swift_dopair_125_%.150s.dat",
815
              outputFileNameExtension);
816
      dump_particle_fields(outputFileName, main_cell, solution, 0);
817
    }
James Willis's avatar
James Willis committed
818 819

    for (int i = 0; i < 125; ++i) {
820 821
      for (int pid = 0; pid < cells[i]->hydro.count; ++pid)
        hydro_init_part(&cells[i]->hydro.parts[pid], &space.hs);
James Willis's avatar
James Willis committed
822
    }
823 824
  }

825
  /* Output timing */
826 827 828 829 830 831 832 833 834 835
  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];

836 837
  ticks self_time = timings[26];

Josh Borrow's avatar
Josh Borrow committed
838 839 840 841 842
  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);
  message("Self calculations took:       %15lli ticks.", self_time / runs);
  message("SWIFT calculation took:       %15lli ticks.", time / runs);
843

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

847 848 849 850
  /* NOW BRUTE-FORCE CALCULATION */

  const ticks tic = getticks();

851 852
  /* Kick the central cell */
  // runner_do_kick1(&runner, main_cell, 0);
853

854 855
  /* And drift it */
  // runner_do_drift_particles(&runner, main_cell, 0);
856

857 858 859
  /* Initialise the particles */
  // for (int j = 0; j < 125; ++j) runner_do_drift_particles(&runner, cells[j],
  // 0);
860

861
  /* Do the density calculation */
862

863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891
  /* 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);
            }
          }
        }
      }
    }
  }
892

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

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

Josh Borrow's avatar
Josh Borrow committed
899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940
#ifdef EXTRA_HYDRO_LOOP
  /* We need to do the gradient loop and the extra ghost! */

  /* 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_gradient(&runner, ci, cj);
            }
          }
        }
      }
    }
  }

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

  /* Extra ghost to finish everything on the central cells */
  for (int j = 0; j < 27; ++j)
    runner_do_extra_ghost(&runner, inner_cells[j], 0);

#endif /* EXTRA_HYDRO_LOOP */

941
  /* Do the force calculation */
942

943 944 945 946
  /* 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++) {
947

948 949 950 951 952 953 954 955 956
        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);
957

958
  /* Finally, give a gentle kick */
959
  runner_do_end_hydro_force(&runner, main_cell, 0);
960
  // runner_do_kick2(&runner, main_cell, 0);
961

962
  const ticks toc = getticks();
963

964
  /* Output timing */
Josh Borrow's avatar
Josh Borrow committed
965
  message("Brute force calculation took: %15lli ticks.", toc - tic);
966

lhausamm's avatar
lhausamm committed
967 968
  sprintf(outputFileName, "brute_force_125_%.150s.dat",
          outputFileNameExtension);
969
  dump_particle_fields(outputFileName, main_cell, solution, 0);
970 971

  /* Clean things to make the sanitizer happy ... */
972
  for (int i = 0; i < 125; ++i) clean_up(cells[i]);
973
  free(solution);
974

975
#ifdef WITH_VECTORIZATION
976 977
  cache_clean(&runner.ci_cache);
  cache_clean(&runner.cj_cache);
978
#endif
979

980 981
  return 0;
}