test125cells.c 26.4 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
#if defined(WITH_VECTORIZATION)
#define DOSELF2 runner_doself2_force_vec
36
#define DOPAIR2 runner_dopair2_branch_force
37
#define DOSELF2_NAME "runner_doself2_force_vec"
38
#define DOPAIR2_NAME "runner_dopair2_force_vec"
39
#endif
40
41
42
43

#ifndef DOSELF2
#define DOSELF2 runner_doself2_force
#define DOSELF2_NAME "runner_doself2_density"
44
45
46
47
#endif

#ifndef DOPAIR2
#define DOPAIR2 runner_dopair2_branch_force
48
49
50
#define DOPAIR2_NAME "runner_dopair2_force"
#endif

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
111
112
113
114
115
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);
116
117
#elif defined(HOPKINS_PE_SPH)
  part->entropy = pressure / pow_gamma(density);
118
119
#elif defined(DEFAULT_SPH)
  part->u = pressure / (hydro_gamma_minus_one * density);
Matthieu Schaller's avatar
Matthieu Schaller committed
120
121
#elif defined(MINIMAL_SPH)
  part->u = pressure / (hydro_gamma_minus_one * density);
122
#elif defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
123
  part->primitives.P = pressure;
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#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) {

152
  for (int i = 0; i < main_cell->count; ++i) {
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
198
199
200
201
202

    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;
203
204
205

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

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

213
  for (int i = 0; i < c->count; ++i) {
214
215
216
217
218

    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
219

220
221
    hydro_init_part(p, hs);

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

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

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

280
281
  float h_max = 0.f;

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

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

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

309
        hydro_first_init_part(part, xpart);
310

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

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

331
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
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
  cell->h_max = h_max;
345
  cell->count = count;
346
  cell->gcount = 0;
347
  cell->dx_max_part = 0.;
348
  cell->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->ti_old_part = 8;
Matthieu Schaller's avatar
Matthieu Schaller committed
357
358
  cell->ti_end_min = 8;
  cell->ti_end_max = 8;
359

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

  cell->sorted = 0;
Matthieu Schaller's avatar
Matthieu Schaller committed
363
  for (int k = 0; k < 13; k++) cell->sort[k] = NULL;
364
365
366
367
368
369

  return cell;
}

void clean_up(struct cell *ci) {
  free(ci->parts);
370
  free(ci->xparts);
371
  for (int k = 0; k < 13; k++)
Matthieu Schaller's avatar
Matthieu Schaller committed
372
    if (ci->sort[k] != NULL) free(ci->sort[k]);
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->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
            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],
402
            main_cell->parts[pid].v[2], main_cell->parts[pid].h,
403
            hydro_get_density(&main_cell->parts[pid]),
404
#if defined(MINIMAL_SPH) || defined(SHADOWFAX_SPH)
Matthieu Schaller's avatar
Matthieu Schaller committed
405
            0.f,
Matthieu Schaller's avatar
Matthieu Schaller committed
406
#else
Matthieu Schaller's avatar
Matthieu Schaller committed
407
408
            main_cell->parts[pid].density.div_v,
#endif
409
410
411
412
            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]),
413
414
            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,
415
#if defined(GADGET2_SPH)
416
417
            main_cell->parts[pid].force.v_sig, main_cell->parts[pid].entropy_dt,
            0.f
418
#elif defined(DEFAULT_SPH)
419
420
            main_cell->parts[pid].force.v_sig, 0.f,
            main_cell->parts[pid].force.u_dt
Matthieu Schaller's avatar
Matthieu Schaller committed
421
#elif defined(MINIMAL_SPH)
422
            main_cell->parts[pid].force.v_sig, 0.f, main_cell->parts[pid].u_dt
423
#else
424
            0.f, 0.f, 0.f
425
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
426
            );
427
428
  }

429
430
431
432
  if (with_solution) {

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

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

449
450
451
452
  fclose(file);
}

/* Just a forward declaration... */
453
454
void runner_dopair1_branch_density(struct runner *r, struct cell *ci,
                                   struct cell *cj);
455
void runner_doself1_density(struct runner *r, struct cell *ci);
456
457
void runner_dopair2_branch_force(struct runner *r, struct cell *ci,
                                 struct cell *cj);
458
void runner_doself2_force(struct runner *r, struct cell *ci);
459
void runner_doself2_force_vec(struct runner *r, struct cell *ci);
460
461
462

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

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

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

477
  /* Choke on FP-exceptions */
478
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
479

480
481
482
483
  /* Get some randomness going */
  srand(0);

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

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

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

557
558
  printf("\n");

559
560
561
562
#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
563

564
565
  /* Build the infrastructure */
  struct space space;
566
  space.periodic = 1;
567
568
569
  space.dim[0] = 5.;
  space.dim[1] = 5.;
  space.dim[2] = 5.;
570
  hydro_space_init(&space.hs, &space);
571

572
573
574
  struct phys_const prog_const;
  prog_const.const_newton_G = 1.f;

575
  struct hydro_props hp;
576
  hp.eta_neighbours = h;
577
  hp.h_tolerance = 1e0;
578
  hp.h_max = FLT_MAX;
579
  hp.max_smoothing_iterations = 1;
580
  hp.CFL_condition = 0.1;
581

582
  struct engine engine;
583
  bzero(&engine, sizeof(struct engine));
584
  engine.hydro_properties = &hp;
585
  engine.physical_constants = &prog_const;
586
587
  engine.s = &space;
  engine.time = 0.1f;
Matthieu Schaller's avatar
Matthieu Schaller committed
588
  engine.ti_current = 8;
589
  engine.max_active_bin = num_time_bins;
590
591
592
593
594

  struct runner runner;
  runner.e = &engine;

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

        /* 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++;
        }
616
617
618
619
620
      }
    }
  }

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

623
624
625
626
627
  /* 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);

628
629
630
  ticks timings[27];
  for (int i = 0; i < 27; i++) timings[i] = 0;

631
  /* Start the test */
632
  ticks time = 0;
633
  for (size_t n = 0; n < runs; ++n) {
634
635
636

    const ticks tic = getticks();

637
    /* Initialise the particles */
638
    for (int j = 0; j < 125; ++j) runner_do_drift_part(&runner, cells[j], 0);
639

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

646
    /* First, sort stuff */
647
648
    for (int j = 0; j < 125; ++j)
      runner_do_sort(&runner, cells[j], 0x1FFF, 0, 0);
649

650
/* Do the density calculation */
651
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
652

Matthieu Schaller's avatar
Matthieu Schaller committed
653
/* Initialise the particle cache. */
654
655
656
657
658
659
660
#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

661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
    /* 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];

683
                if (cj > ci) runner_dopair1_branch_density(&runner, ci, cj);
684
685
686
687
688
689
690
691
              }
            }
          }
        }
      }
    }

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

#endif
696

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

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

703
704
705
706
707
708
709
710
#ifdef WITH_VECTORIZATION
    /* Initialise the cache. */
    runner.ci_cache.count = 0;
    runner.cj_cache.count = 0;
    cache_init(&runner.ci_cache, 512);
    cache_init(&runner.cj_cache, 512);
#endif

711
    int ctr = 0;
712
713
714
715
716
717
718
    /* 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];

719
          if (main_cell != cj) {
James Willis's avatar
James Willis committed
720

721
722
            const ticks sub_tic = getticks();

723
            DOPAIR2(&runner, main_cell, cj);
724

725
            timings[ctr++] += getticks() - sub_tic;
726
          }
727
728
729
730
        }
      }
    }

731
732
    ticks self_tic = getticks();

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

736
    timings[26] += getticks() - self_tic;
737
738
#endif

739
    /* Finally, give a gentle kick */
740
    runner_do_end_force(&runner, main_cell, 0);
741
742
743
744
    const ticks toc = getticks();
    time += toc - tic;

    /* Dump if necessary */
745
    if (n == 0) {
746
      sprintf(outputFileName, "swift_dopair_125_%s.dat",
747
              outputFileNameExtension);
748
      dump_particle_fields(outputFileName, main_cell, solution, 0);
749
    }
750

James Willis's avatar
James Willis committed
751
    for (int i = 0; i < 125; ++i) {
752
753
      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
754
    }
755
  }
756
757

  /* Output timing */
758
759
760
761
762
763
764
765
766
767
  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];

768
769
  ticks self_time = timings[26];

770
771
772
  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);
773
  message("Self calculations took         : %15lli ticks.", self_time / runs);
774
  message("SWIFT calculation took         : %15lli ticks.", time / runs);
775

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

779
780
781
782
  /* NOW BRUTE-FORCE CALCULATION */

  const ticks tic = getticks();

783
784
/* Kick the central cell */
// runner_do_kick1(&runner, main_cell, 0);
785

786
787
/* And drift it */
// runner_do_drift_particles(&runner, main_cell, 0);
788

789
790
791
/* Initialise the particles */
// for (int j = 0; j < 125; ++j) runner_do_drift_particles(&runner, cells[j],
// 0);
792
793

/* Do the density calculation */
794
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
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
820
821
822
823
824
  /* 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);
            }
          }
        }
      }
    }
  }
825

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

829
#endif
830

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

834
/* Do the force calculation */
835
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
836

837
838
839
840
  /* 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++) {
841

842
843
844
845
846
847
848
849
850
        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);
851

852
#endif
853

854
  /* Finally, give a gentle kick */
855
  runner_do_end_force(&runner, main_cell, 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
856
  // runner_do_kick2(&runner, main_cell, 0);
857

858
  const ticks toc = getticks();
859

860
861
  /* Output timing */
  message("Brute force calculation took : %15lli ticks.", toc - tic);
862

863
864
  sprintf(outputFileName, "brute_force_125_%s.dat", outputFileNameExtension);
  dump_particle_fields(outputFileName, main_cell, solution, 0);
865
866

  /* Clean things to make the sanitizer happy ... */
867
  for (int i = 0; i < 125; ++i) clean_up(cells[i]);
868
  free(solution);
869
870
871

  return 0;
}