test125cells.c 24.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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
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);
99
100
#elif defined(HOPKINS_PE_SPH)
  part->entropy = pressure / pow_gamma(density);
101
102
#elif defined(DEFAULT_SPH)
  part->u = pressure / (hydro_gamma_minus_one * density);
Matthieu Schaller's avatar
Matthieu Schaller committed
103
104
#elif defined(MINIMAL_SPH)
  part->u = pressure / (hydro_gamma_minus_one * density);
105
#elif defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
106
  part->primitives.P = pressure;
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
#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) {

135
  for (int i = 0; i < main_cell->count; ++i) {
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185

    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;
186
187
188

    solution[i].S_dt = 0.f;
    solution[i].u_dt = -(solution[i].P / solution[i].rho) * solution[i].div_v;
189
190
191
  }
}

192
193
194
void reset_particles(struct cell *c, struct hydro_space *hs,
                     enum velocity_field vel, enum pressure_field press,
                     float size, float density) {
195

196
  for (int i = 0; i < c->count; ++i) {
197
198
199
200
201

    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
202

203
204
    hydro_init_part(p, hs);

205
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
206
    float volume = p->conserved.mass / density;
207
208
209
210
211
#if defined(GIZMO_SPH)
    p->geometry.volume = volume;
#else
    p->cell.volume = volume;
#endif
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
212
213
214
215
216
217
218
219
    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 =
220
        p->primitives.P / hydro_gamma_minus_one * volume +
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
221
222
223
224
225
        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
226
227
228
  }
}

229
230
/**
 * @brief Constructs a cell and all of its particle in a valid state prior to
231
 * a SPH time-step.
232
233
234
235
236
 *
 * @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
237
 * separation.
238
239
 * @param density The density of the fluid.
 * @param partId The running counter of IDs.
240
241
 * @param pert The perturbation to apply to the particles in the cell in units
 *of the inter-particle separation.
242
243
 * @param vel The type of velocity field.
 * @param press The type of pressure field.
244
 */
245
struct cell *make_cell(size_t n, const double offset[3], double size, double h,
James Willis's avatar
James Willis committed
246
                       double density, long long *partId, double pert,
247
                       enum velocity_field vel, enum pressure_field press) {
248

249
250
251
252
253
254
  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,
255
                     count * sizeof(struct part)) != 0)
256
    error("couldn't allocate particles, no. of particles: %d", (int)count);
257
258
259
  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);
260
  bzero(cell->parts, count * sizeof(struct part));
261
  bzero(cell->xparts, count * sizeof(struct xpart));
262
263
264

  /* Construct the parts */
  struct part *part = cell->parts;
265
  struct xpart *xpart = cell->xparts;
266
267
268
  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
269
270
271
272
273
274
275
276
277
        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;
278
        part->h = size * h / (float)n;
279

280
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
281
        part->conserved.mass = density * volume / count;
282
#else
283
        part->mass = density * volume / count;
284
#endif
285
286
287
288

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

289
        hydro_first_init_part(part, xpart);
290

291
        part->id = ++(*partId);
292
        part->time_bin = 1;
293

294
#if defined(GIZMO_SPH)
295
        part->geometry.volume = part->conserved.mass / density;
296
297
298
299
300
301
302
303
        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 =
304
            part->primitives.P / hydro_gamma_minus_one * volume +
305
306
307
308
309
310
            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

311
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
312
313
        part->ti_drift = 8;
        part->ti_kick = 8;
314
315
#endif

316
        ++part;
317
        ++xpart;
318
319
320
321
322
323
324
325
      }
    }
  }

  /* Cell properties */
  cell->split = 0;
  cell->h_max = h;
  cell->count = count;
326
  cell->gcount = 0;
327
  cell->dx_max_part = 0.;
328
  cell->dx_max_sort = 0.;
Matthieu Schaller's avatar
Matthieu Schaller committed
329
330
331
  cell->width[0] = size;
  cell->width[1] = size;
  cell->width[2] = size;
332
333
334
335
  cell->loc[0] = offset[0];
  cell->loc[1] = offset[1];
  cell->loc[2] = offset[2];

336
  cell->ti_old_part = 8;
Matthieu Schaller's avatar
Matthieu Schaller committed
337
338
  cell->ti_end_min = 8;
  cell->ti_end_max = 8;
339
  cell->ti_sort = 0;
340

341
  // shuffle_particles(cell->parts, cell->count);
342
343
344
345
346
347
348
349
350
351

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

  return cell;
}

void clean_up(struct cell *ci) {
  free(ci->parts);
352
  free(ci->xparts);
353
354
355
356
357
358
359
360
  free(ci->sort);
  free(ci);
}

/**
 * @brief Dump all the particles to a file
 */
void dump_particle_fields(char *fileName, struct cell *main_cell,
361
                          struct solution_part *solution, int with_solution) {
362
363
364
365
  FILE *file = fopen(fileName, "w");

  /* Write header */
  fprintf(file,
James Willis's avatar
James Willis committed
366
          "# %4s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %13s %13s "
367
          "%13s %13s %13s %8s %8s\n",
368
369
370
          "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");
371
372
373
374

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

  /* Write main cell */
375
  for (int pid = 0; pid < main_cell->count; pid++) {
376
    fprintf(file,
377
378
            "%6llu %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f "
            "%8.5f "
379
            "%8.5f %8.5f %13e %13e %13e %13e %13e %8.5f %8.5f\n",
380
381
382
            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],
383
            main_cell->parts[pid].v[2], main_cell->parts[pid].h,
384
            hydro_get_density(&main_cell->parts[pid]),
385
#if defined(MINIMAL_SPH) || defined(SHADOWFAX_SPH)
Matthieu Schaller's avatar
Matthieu Schaller committed
386
            0.f,
Matthieu Schaller's avatar
Matthieu Schaller committed
387
#else
Matthieu Schaller's avatar
Matthieu Schaller committed
388
389
            main_cell->parts[pid].density.div_v,
#endif
390
391
392
393
            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]),
394
395
            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,
396
#if defined(GADGET2_SPH)
397
398
            main_cell->parts[pid].force.v_sig, main_cell->parts[pid].entropy_dt,
            0.f
399
#elif defined(DEFAULT_SPH)
400
401
            main_cell->parts[pid].force.v_sig, 0.f,
            main_cell->parts[pid].force.u_dt
Matthieu Schaller's avatar
Matthieu Schaller committed
402
#elif defined(MINIMAL_SPH)
403
            main_cell->parts[pid].force.v_sig, 0.f, main_cell->parts[pid].u_dt
404
#else
405
            0.f, 0.f, 0.f
406
#endif
Matthieu Schaller's avatar
Matthieu Schaller committed
407
            );
408
409
  }

410
411
412
413
  if (with_solution) {

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

414
    for (int pid = 0; pid < main_cell->count; pid++) {
415
416
417
      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 "
418
              "%8.5f %8.5f %13f %13f %13f %13f %13f %8.5f %8.5f\n",
419
420
421
422
423
424
              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],
425
426
              solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt,
              solution[pid].u_dt);
427
428
    }
  }
429

430
431
432
433
434
  fclose(file);
}

/* Just a forward declaration... */
void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
435
436
void runner_dopair1_branch_density(struct runner *r, struct cell *ci,
                                   struct cell *cj);
437
void runner_doself1_density(struct runner *r, struct cell *ci);
438
439
void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself2_force(struct runner *r, struct cell *ci);
440
441
442

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

444
  size_t runs = 0, particles = 0;
445
  double h = 1.23485, size = 1., rho = 2.5;
James Willis's avatar
James Willis committed
446
  double perturbation = 0.;
447
448
  char outputFileNameExtension[200] = "";
  char outputFileName[200] = "";
449
450
  enum velocity_field vel = velocity_zero;
  enum pressure_field press = pressure_const;
451
452
453
454
455

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

456
  /* Choke on FP-exceptions */
457
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
458

459
460
461
462
  /* Get some randomness going */
  srand(0);

  char c;
463
  while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:p:")) != -1) {
464
465
466
467
468
469
470
    switch (c) {
      case 'h':
        sscanf(optarg, "%lf", &h);
        break;
      case 's':
        sscanf(optarg, "%lf", &size);
        break;
471
      case 'n':
472
473
474
475
476
        sscanf(optarg, "%zu", &particles);
        break;
      case 'r':
        sscanf(optarg, "%zu", &runs);
        break;
James Willis's avatar
James Willis committed
477
478
479
      case 'd':
        sscanf(optarg, "%lf", &perturbation);
        break;
480
481
482
483
484
485
      case 'm':
        sscanf(optarg, "%lf", &rho);
        break;
      case 'f':
        strcpy(outputFileNameExtension, optarg);
        break;
486
487
488
489
490
491
      case 'v':
        sscanf(optarg, "%d", (int *)&vel);
        break;
      case 'p':
        sscanf(optarg, "%d", (int *)&press);
        break;
492
493
494
495
496
497
498
499
      case '?':
        error("Unknown option.");
        break;
    }
  }

  if (h < 0 || particles == 0 || runs == 0) {
    printf(
500
501
502
503
504
        "\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()"
505
506
507
508
        "\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
509
        "\n-d pert            - Perturbation to apply to the particles [0,1["
510
        "\n-v type (0,1,2,3)  - Velocity field: (zero, constant, divergent, "
511
        "rotating)"
512
        "\n-p type (0,1,2)    - Pressure field: (constant, gradient divergent)"
513
514
515
516
517
518
        "\n-f fileName        - Part of the file name used to save the dumps\n",
        argv[0]);
    exit(1);
  }

  /* Help users... */
519
  message("Adiabatic index: ga = %f", hydro_gamma);
520
  message("Hydro implementation: %s", SPH_IMPLEMENTATION);
521
522
  message("Smoothing length: h = %f", h * size);
  message("Kernel:               %s", kernel_name);
523
  message("Neighbour target: N = %f", pow_dimension(h) * kernel_norm);
524
  message("Density target: rho = %f", rho);
525
526
527
528
529
530
531
532
533
  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");

534
535
  printf("\n");

536
537
538
539
#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
540

541
542
  /* Build the infrastructure */
  struct space space;
543
  space.periodic = 1;
544
545
546
  space.dim[0] = 5.;
  space.dim[1] = 5.;
  space.dim[2] = 5.;
547
  hydro_space_init(&space.hs, &space);
548

549
550
551
  struct phys_const prog_const;
  prog_const.const_newton_G = 1.f;

552
  struct hydro_props hp;
553
  hp.target_neighbours = pow_dimension(h) * kernel_norm;
554
555
  hp.delta_neighbours = 4.;
  hp.h_max = FLT_MAX;
556
  hp.max_smoothing_iterations = 1;
557
  hp.CFL_condition = 0.1;
558

559
  struct engine engine;
560
  bzero(&engine, sizeof(struct engine));
561
  engine.hydro_properties = &hp;
562
  engine.physical_constants = &prog_const;
563
564
  engine.s = &space;
  engine.time = 0.1f;
Matthieu Schaller's avatar
Matthieu Schaller committed
565
  engine.ti_current = 8;
566
  engine.max_active_bin = num_time_bins;
567
568
569
570
571

  struct runner runner;
  runner.e = &engine;

  /* Construct some cells */
572
573
  struct cell *cells[125];
  struct cell *inner_cells[27];
574
  struct cell *main_cell;
575
  int count = 0;
576
  static long long partId = 0;
577
578
579
580
581
582
583
584
  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 */
585
586
        cells[i * 25 + j * 5 + k] = make_cell(
            particles, offset, size, h, rho, &partId, perturbation, vel, press);
587
588
589
590
591
592

        /* 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++;
        }
593
594
595
596
597
      }
    }
  }

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

600
601
602
603
604
605
  /* Construct the real solution */
  struct solution_part *solution =
      malloc(main_cell->count * sizeof(struct solution_part));
  get_solution(main_cell, solution, rho, vel, press, size);

  /* Start the test */
606
  ticks time = 0;
607
  for (size_t n = 0; n < runs; ++n) {
608
609
610

    const ticks tic = getticks();

611
    /* Initialise the particles */
612
    for (int j = 0; j < 125; ++j) runner_do_drift_part(&runner, cells[j], 0);
613

James Willis's avatar
James Willis committed
614
615
616
617
618
    /* Reset particles. */
    for (int i = 0; i < 125; ++i) {
      for (int n = 0; n < cells[i]->count; ++n)
        hydro_init_part(&cells[i]->parts[n], &space.hs);
    }
619

620
621
622
    /* First, sort stuff */
    for (int j = 0; j < 125; ++j) runner_do_sort(&runner, cells[j], 0x1FFF, 0);

623
/* Do the density calculation */
624
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
625

626
627
628
629
630
631
632
633
    /* Initialise the particle cache. */
#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

634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
    /* 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];

656
                if (cj > ci) runner_dopair1_branch_density(&runner, ci, cj);
657
658
659
660
661
662
663
664
              }
            }
          }
        }
      }
    }

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

#endif
669

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

/* Do the force calculation */
674
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
675

676
677
678
679
680
681
682
683
684
685
686
687
688
689
    /* Do the pairs (for the central 27 cells) */
    for (int i = 1; i < 4; i++) {
      for (int j = 1; j < 4; j++) {
        for (int k = 1; k < 4; k++) {

          struct cell *cj = cells[i * 25 + j * 5 + k];

          if (main_cell != cj) runner_dopair2_force(&runner, main_cell, cj);
        }
      }
    }

    /* And now the self-interaction for the main cell */
    runner_doself2_force(&runner, main_cell);
690
691
#endif

692
    /* Finally, give a gentle kick */
693
    runner_do_end_force(&runner, main_cell, 0);
694
695
696
697
    const ticks toc = getticks();
    time += toc - tic;

    /* Dump if necessary */
698
    if (n == 0) {
699
      sprintf(outputFileName, "swift_dopair_125_%s.dat",
700
              outputFileNameExtension);
701
      dump_particle_fields(outputFileName, main_cell, solution, 0);
702
    }
703
704
705
706
707
708

    /* Reset stuff */
    for (int i = 0; i < 125; ++i) {
      for (int n = 0; n < cells[i]->count; ++n)
        hydro_init_part(&cells[i]->parts[n], &space.hs);
    }
709
710
711
712
713
  }

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

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

717
718
719
720
  /* NOW BRUTE-FORCE CALCULATION */

  const ticks tic = getticks();

721
722
/* Kick the central cell */
// runner_do_kick1(&runner, main_cell, 0);
723

724
725
/* And drift it */
// runner_do_drift_particles(&runner, main_cell, 0);
726

727
728
729
/* Initialise the particles */
// for (int j = 0; j < 125; ++j) runner_do_drift_particles(&runner, cells[j],
// 0);
730
731

/* Do the density calculation */
732
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
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
  /* 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);
            }
          }
        }
      }
    }
  }
763

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

767
#endif
768

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

772
/* Do the force calculation */
773
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
774

775
776
777
778
  /* 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++) {
779

780
781
782
783
784
785
786
787
788
        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);
789

790
#endif
791

792
  /* Finally, give a gentle kick */
793
  runner_do_end_force(&runner, main_cell, 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
794
  // runner_do_kick2(&runner, main_cell, 0);
795

796
  const ticks toc = getticks();
797

798
799
  /* Output timing */
  message("Brute force calculation took : %15lli ticks.", toc - tic);
800

801
802
  sprintf(outputFileName, "brute_force_125_%s.dat", outputFileNameExtension);
  dump_particle_fields(outputFileName, main_cell, solution, 0);
803
804

  /* Clean things to make the sanitizer happy ... */
805
  for (int i = 0; i < 125; ++i) clean_up(cells[i]);
806
  free(solution);
807
808
809

  return 0;
}