test125cells.c 24.2 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
void runner_dopair1_branch_density(struct runner *r, struct cell *ci, struct cell *cj);
436
void runner_doself1_density(struct runner *r, struct cell *ci);
437
438
void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself2_force(struct runner *r, struct cell *ci);
439
440
441

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

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

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

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

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

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

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

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

533
534
  printf("\n");

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

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

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

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

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

  struct runner runner;
  runner.e = &engine;

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

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

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

599
600
601
602
603
604
  /* 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 */
605
  ticks time = 0;
606
  for (size_t n = 0; n < runs; ++n) {
607
608
609

    const ticks tic = getticks();

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

James Willis's avatar
James Willis committed
613
614
615
616
617
    /* 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);
    }
618

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

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

625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
    /* 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];

647
                if (cj > ci) runner_dopair1_branch_density(&runner, ci, cj);
648
649
650
651
652
653
654
655
              }
            }
          }
        }
      }
    }

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

#endif
660

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

/* Do the force calculation */
665
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
666

667
668
669
670
671
672
673
674
675
676
677
678
679
680
    /* 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);
681
682
#endif

683
    /* Finally, give a gentle kick */
684
    runner_do_end_force(&runner, main_cell, 0);
685
686
687
688
    const ticks toc = getticks();
    time += toc - tic;

    /* Dump if necessary */
689
    if (n == 0) {
690
      sprintf(outputFileName, "swift_dopair_125_%s.dat",
691
              outputFileNameExtension);
692
      dump_particle_fields(outputFileName, main_cell, solution, 0);
693
    }
694
695
696
697
698
699

    /* 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);
    }
700
701
702
703
704
  }

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

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

708
709
710
711
  /* NOW BRUTE-FORCE CALCULATION */

  const ticks tic = getticks();

712
713
/* Kick the central cell */
// runner_do_kick1(&runner, main_cell, 0);
714

715
716
/* And drift it */
// runner_do_drift_particles(&runner, main_cell, 0);
717

718
719
720
/* Initialise the particles */
// for (int j = 0; j < 125; ++j) runner_do_drift_particles(&runner, cells[j],
// 0);
721
722

/* Do the density calculation */
723
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
724

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
  /* 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);
            }
          }
        }
      }
    }
  }
754

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

758
#endif
759

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

763
/* Do the force calculation */
764
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
765

766
767
768
769
  /* 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++) {
770

771
772
773
774
775
776
777
778
779
        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);
780

781
#endif
782

783
  /* Finally, give a gentle kick */
784
  runner_do_end_force(&runner, main_cell, 0);
Matthieu Schaller's avatar
Matthieu Schaller committed
785
  // runner_do_kick2(&runner, main_cell, 0);
786

787
  const ticks toc = getticks();
788

789
790
  /* Output timing */
  message("Brute force calculation took : %15lli ticks.", toc - tic);
791

792
793
  sprintf(outputFileName, "brute_force_125_%s.dat", outputFileNameExtension);
  dump_particle_fields(outputFileName, main_cell, solution, 0);
794
795

  /* Clean things to make the sanitizer happy ... */
796
  for (int i = 0; i < 125; ++i) clean_up(cells[i]);
797
  free(solution);
798
799
800

  return 0;
}