test125cells.c 24.1 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
435
  fclose(file);
}

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

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

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

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

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

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

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

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

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

532
533
  printf("\n");

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

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

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

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

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

  struct runner runner;
  runner.e = &engine;

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

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

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

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

    const ticks tic = getticks();

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

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

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

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

624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
    /* 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_density(&runner, ci, cj);
              }
            }
          }
        }
      }
    }

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

#endif
659

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

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

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

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

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

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

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

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

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

  const ticks tic = getticks();

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

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

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

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

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

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

757
#endif
758

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

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

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

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

780
#endif
781

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

786
  const ticks toc = getticks();
787

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

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

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

  return 0;
}