test125cells.c 24.2 KB
Newer Older
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Copyright (C) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/

20
21
22
23
/* Config parameters. */
#include "../config.h"

/* Some standard headers. */
24
25
26
27
28
29
#include <fenv.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

30
31
/* Local headers. */
#include "swift.h"
32

James Willis's avatar
James Willis committed
33
34
35
36
37
38
39
40
41
42
43
44
45
//#if defined(WITH_VECTORIZATION)
//#define DOSELF2 runner_doself2_force_vec
////#define DOPAIR2 runner_dopair2_force_vec
//#define DOSELF2_NAME "runner_doself2_force_vec"
//#define DOPAIR2_NAME "runner_dopair2_force"
//#endif

#ifndef DOSELF2
#define DOSELF2 runner_doself2_force
#define DOSELF2_NAME "runner_doself2_density"
#define DOPAIR2_NAME "runner_dopair2_force"
#endif

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

147
  for (int i = 0; i < main_cell->count; ++i) {
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
186
187
188
189
190
191
192
193
194
195
196
197

    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;
198
199
200

    solution[i].S_dt = 0.f;
    solution[i].u_dt = -(solution[i].P / solution[i].rho) * solution[i].div_v;
201
202
203
  }
}

204
205
206
void reset_particles(struct cell *c, struct hydro_space *hs,
                     enum velocity_field vel, enum pressure_field press,
                     float size, float density) {
207

208
  for (int i = 0; i < c->count; ++i) {
209
210
211
212
213

    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
214

215
216
    hydro_init_part(p, hs);

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

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

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

  /* Construct the parts */
  struct part *part = cell->parts;
275
  struct xpart *xpart = cell->xparts;
276
277
278
  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
279
280
281
282
283
284
285
286
287
        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;
288
        part->h = size * h / (float)n;
289

290
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
291
        part->conserved.mass = density * volume / count;
292
#else
293
        part->mass = density * volume / count;
294
#endif
295
296
297
298

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

299
        hydro_first_init_part(part, xpart);
300

301
        part->id = ++(*partId);
302
        part->time_bin = 1;
303

304
#if defined(GIZMO_SPH)
305
        part->geometry.volume = part->conserved.mass / density;
306
307
308
309
310
311
312
313
        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 =
314
            part->primitives.P / hydro_gamma_minus_one * volume +
315
316
317
318
319
320
            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

321
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
322
323
        part->ti_drift = 8;
        part->ti_kick = 8;
324
325
#endif

326
        ++part;
327
        ++xpart;
328
329
330
331
332
333
334
335
      }
    }
  }

  /* Cell properties */
  cell->split = 0;
  cell->h_max = h;
  cell->count = count;
336
  cell->gcount = 0;
337
  cell->dx_max = 0.;
338
  cell->dx_max_sort = 0.;
Matthieu Schaller's avatar
Matthieu Schaller committed
339
340
341
  cell->width[0] = size;
  cell->width[1] = size;
  cell->width[2] = size;
342
343
344
345
  cell->loc[0] = offset[0];
  cell->loc[1] = offset[1];
  cell->loc[2] = offset[2];

Matthieu Schaller's avatar
Matthieu Schaller committed
346
347
348
  cell->ti_old = 8;
  cell->ti_end_min = 8;
  cell->ti_end_max = 8;
349
  cell->ti_sort = 0;
350

351
  // shuffle_particles(cell->parts, cell->count);
352
353
354
355
356
357
358
359
360
361

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

  return cell;
}

void clean_up(struct cell *ci) {
  free(ci->parts);
362
  free(ci->xparts);
363
364
365
366
367
368
369
370
  free(ci->sort);
  free(ci);
}

/**
 * @brief Dump all the particles to a file
 */
void dump_particle_fields(char *fileName, struct cell *main_cell,
371
                          struct solution_part *solution, int with_solution) {
372
373
374
375
  FILE *file = fopen(fileName, "w");

  /* Write header */
  fprintf(file,
James Willis's avatar
James Willis committed
376
377
          "# %4s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %13s %13s "
          "%13s %8s %8s %8s %8s\n",
378
379
380
          "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");
381
382
383
384

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

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

420
421
422
423
  if (with_solution) {

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

424
    for (int pid = 0; pid < main_cell->count; pid++) {
425
426
427
428
429
430
431
432
433
434
      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 "
              "%8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f\n",
              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],
435
436
              solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt,
              solution[pid].u_dt);
437
438
    }
  }
439

440
441
442
443
444
445
  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);
446
447
void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself2_force(struct runner *r, struct cell *ci);
448
449
450

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

452
  size_t runs = 0, particles = 0;
453
  double h = 1.23485, size = 1., rho = 2.5;
James Willis's avatar
James Willis committed
454
  double perturbation = 0.;
455
456
  char outputFileNameExtension[200] = "";
  char outputFileName[200] = "";
457
458
  enum velocity_field vel = velocity_zero;
  enum pressure_field press = pressure_const;
459
460
461
462
463

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

464
  /* Choke on FP-exceptions */
465
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
466

467
468
469
470
  /* Get some randomness going */
  srand(0);

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

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

  /* Help users... */
527
  message("Adiabatic index: ga = %f", hydro_gamma);
528
  message("Hydro implementation: %s", SPH_IMPLEMENTATION);
529
530
  message("Smoothing length: h = %f", h * size);
  message("Kernel:               %s", kernel_name);
531
  message("Neighbour target: N = %f", pow_dimension(h) * kernel_norm);
532
  message("Density target: rho = %f", rho);
533
534
535
536
537
538
539
540
541
  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");

542
543
  printf("\n");

544
545
546
547
#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
548

549
550
  /* Build the infrastructure */
  struct space space;
551
552
553
554
555
  space.periodic = 1;
  space.dim[0] = 3.;
  space.dim[1] = 3.;
  space.dim[2] = 3.;
  hydro_space_init(&space.hs, &space);
556

557
558
559
  struct phys_const prog_const;
  prog_const.const_newton_G = 1.f;

560
  struct hydro_props hp;
561
562
  hp.target_neighbours = pow_dimension(h) * kernel_norm;
  hp.delta_neighbours = 2.;
563
  hp.max_smoothing_iterations = 1;
564
  hp.CFL_condition = 0.1;
565

566
  struct engine engine;
567
  bzero(&engine, sizeof(struct engine));
568
  engine.hydro_properties = &hp;
569
  engine.physical_constants = &prog_const;
570
571
  engine.s = &space;
  engine.time = 0.1f;
Matthieu Schaller's avatar
Matthieu Schaller committed
572
  engine.ti_current = 8;
573
  engine.max_active_bin = num_time_bins;
574
575
576
577
578

  struct runner runner;
  runner.e = &engine;

  /* Construct some cells */
579
580
  struct cell *cells[125];
  struct cell *inner_cells[27];
581
  struct cell *main_cell;
582
  int count = 0;
583
  static long long partId = 0;
584
585
586
587
588
589
590
591
592
  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 */
        cells[i * 25 + j * 5 + k] =
James Willis's avatar
James Willis committed
593
            make_cell(particles, offset, size, h, rho, &partId, perturbation, vel, press);
594
595
596
597
598
599

        /* 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++;
        }
600
601
602
603
604
      }
    }
  }

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

607
608
609
610
611
612
  /* 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 */
613
  ticks time = 0;
614
  for (size_t n = 0; n < runs; ++n) {
615
616
617

    const ticks tic = getticks();

618
619
620
    /* Initialise the particles */
    for (int j = 0; j < 125; ++j)
      runner_do_drift_particles(&runner, cells[j], 0);
621

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

625
/* Do the density calculation */
626
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
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
655
656
657
658
    /* 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*/
659
    for (int j = 0; j < 27; ++j)
660
661
662
      runner_doself1_density(&runner, inner_cells[j]);

#endif
663

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

/* Do the force calculation */
668
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
669

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

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

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

James Willis's avatar
James Willis committed
699
700
701
702
703
704
705
706
707
  //for (size_t n = 0; n < 100*runs; ++n) {
  //  ticks self_tic = getticks();

  //  DOSELF2(&runner, main_cell);

  //  self_force_time += getticks() - self_tic;
  //  
  //}

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

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

714
715
716
717
  /* NOW BRUTE-FORCE CALCULATION */

  const ticks tic = getticks();

718
719
/* Kick the central cell */
// runner_do_kick1(&runner, main_cell, 0);
720

721
722
/* And drift it */
// runner_do_drift_particles(&runner, main_cell, 0);
723

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

/* Do the density calculation */
729
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
730

731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
  /* 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);
            }
          }
        }
      }
    }
  }
760

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

764
#endif
765

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

769
/* Do the force calculation */
770
#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
771

772
773
774
775
  /* 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++) {
776

777
778
779
780
781
782
783
784
785
        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);
786

787
#endif
788

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

793
  const ticks toc = getticks();
794

795
796
  /* Output timing */
  message("Brute force calculation took : %15lli ticks.", toc - tic);
797

798
799
  sprintf(outputFileName, "brute_force_125_%s.dat", outputFileNameExtension);
  dump_particle_fields(outputFileName, main_cell, solution, 0);
800
801

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

  return 0;
}