test27cells.c 9.62 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
 *
 * 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/>.
 *
 ******************************************************************************/

#include <fenv.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <unistd.h>
#include "swift.h"

/**
 * Returns a random number (uniformly distributed) in [a,b[
 */
double random_uniform(double a, double b) {
  return (rand() / (double)RAND_MAX) * (a - b) + a;
}

/* n is both particles per axis and box size:
 * particles are generated on a mesh with unit spacing
 */
37
38
39
40
struct cell *make_cell(size_t n, double *offset, double size, double h,
                       double density, long long *partId, double pert) {
  const size_t count = n * n * n;
  const double volume = size * size * size;
41
42
43
  struct cell *cell = malloc(sizeof(struct cell));
  bzero(cell, sizeof(struct cell));

44
45
  if (posix_memalign((void **)&cell->parts, part_align,
                     count * sizeof(struct part)) != 0) {
46
47
48
49
    error("couldn't allocate particles, no. of particles: %d", (int)count);
  }
  bzero(cell->parts, count * sizeof(struct part));

50
51
52
53
54
  /* Construct the parts */
  struct part *part = cell->parts;
  for (size_t x = 0; x < n; ++x) {
    for (size_t y = 0; y < n; ++y) {
      for (size_t z = 0; z < n; ++z) {
55
        // Add .5 for symmetry: 0.5, 1.5, 2.5 vs. 0, 1, 2
56
57
58
59
60
61
62
63
64
65
66
67
68
        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;
        part->v[0] = 1. * random_uniform(-0.1, 0.1);
        part->v[1] = 1. * random_uniform(-0.1, 0.1);
        part->v[2] = 1. * random_uniform(-0.1, 0.1);
        part->h = size * h / (float)n;
69
        part->id = ++(*partId);
70
        part->mass = density * volume / count;
71
72
73
74
75
76
77
        part->ti_begin = 0;
        part->ti_end = 1;
        ++part;
      }
    }
  }

78
  /* Cell properties */
79
80
81
82
  cell->split = 0;
  cell->h_max = h;
  cell->count = count;
  cell->dx_max = 0.;
83
84
85
  cell->h[0] = size;
  cell->h[1] = size;
  cell->h[2] = size;
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
  cell->loc[0] = offset[0];
  cell->loc[1] = offset[1];
  cell->loc[2] = offset[2];

  cell->ti_end_min = 1;
  cell->ti_end_max = 1;

  cell->sorted = 0;
  cell->sort = NULL;
  cell->sortsize = 0;
  runner_dosort(NULL, cell, 0x1FFF, 0);

  return cell;
}

void clean_up(struct cell *ci) {
  free(ci->parts);
  free(ci->sort);
  free(ci);
}

/**
 * @brief Initializes all particles field to be ready for a density calculation
 */
void zero_particle_fields(struct cell *c) {

  for (size_t pid = 0; pid < c->count; pid++) {
    c->parts[pid].rho = 0.f;
    c->parts[pid].rho_dh = 0.f;
    hydro_init_part(&c->parts[pid]);
  }
}

119
120
121
122
123
124
125
126
127
128
/**
 * @brief Ends the loop by adding the appropriate coefficients
 */
void end_calculation(struct cell *c) {

  for (size_t pid = 0; pid < c->count; pid++) {
    hydro_end_density(&c->parts[pid], 1);
  }
}

129
130
131
/**
 * @brief Dump all the particles to a file
 */
132
133
void dump_particle_fields(char *fileName, struct cell *main_cell,
                          struct cell **cells) {
134
135
136
137
138
139
140

  FILE *file = fopen(fileName, "w");

  fprintf(file,
          "# ID  pos:[x y z]  rho  rho_dh  wcount  wcount_dh  div_v  curl_v:[x "
          "y z]\n");

141
142
143
144
145
146
147
148
149
150
151
152
  fprintf(file, "# -----------------------------------\n");

  for (size_t pid = 0; pid < main_cell->count; pid++) {
    fprintf(file, "%6llu %f %f %f %f %f %f %f %f %f %f %f\n",
            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].rho, main_cell->parts[pid].rho_dh,
            main_cell->parts[pid].density.wcount,
            main_cell->parts[pid].density.wcount_dh,
            main_cell->parts[pid].div_v, main_cell->parts[pid].density.rot_v[0],
            main_cell->parts[pid].density.rot_v[1],
            main_cell->parts[pid].density.rot_v[2]);
153
154
  }

155
156
157
158
159
160
  for (int j = 0; j < 27; ++j) {

    struct cell *cj = cells[j];
    if (cj == main_cell) continue;

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

162
163
164
165
166
167
168
169
    for (size_t pjd = 0; pjd < cj->count; pjd++) {
      fprintf(file, "%6llu %f %f %f %f %f %f %f %f %f %f %f\n",
              cj->parts[pjd].id, cj->parts[pjd].x[0], cj->parts[pjd].x[1],
              cj->parts[pjd].x[2], cj->parts[pjd].rho, cj->parts[pjd].rho_dh,
              cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
              cj->parts[pjd].div_v, cj->parts[pjd].density.rot_v[0],
              cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]);
    }
170
171
172
173
174
175
176
  }

  fclose(file);
}

/* Just a forward declaration... */
void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
177
void runner_doself1_density(struct runner *r, struct cell *ci);
178

179
/* And go... */
180
int main(int argc, char *argv[]) {
181
182
183
184

  size_t runs = 0, particles = 0;
  double h = 1.1255, size = 1., rho = 1.;
  double perturbation = 0.;
185
186
187
188
189
190
191
  char outputFileNameExtension[200] = "";
  char outputFileName[200] = "";

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

192
  /* Get some randomness going */
193
194
  srand(0);

195
196
  char c;
  while ((c = getopt(argc, argv, "m:s:h:p:r:t:d:f:")) != -1) {
197
198
199
200
    switch (c) {
      case 'h':
        sscanf(optarg, "%lf", &h);
        break;
201
202
203
      case 's':
        sscanf(optarg, "%lf", &size);
        break;
204
205
206
207
208
209
210
211
212
      case 'p':
        sscanf(optarg, "%zu", &particles);
        break;
      case 'r':
        sscanf(optarg, "%zu", &runs);
        break;
      case 'd':
        sscanf(optarg, "%lf", &perturbation);
        break;
213
214
215
      case 'm':
        sscanf(optarg, "%lf", &rho);
        break;
216
217
218
219
220
221
222
223
224
      case 'f':
        strcpy(outputFileNameExtension, optarg);
        break;
      case '?':
        error("Unknown option.");
        break;
    }
  }

225
  if (h < 0 || particles == 0 || runs == 0) {
226
227
228
229
230
    printf(
        "\nUsage: %s -p PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
        "\nGenerates a cell pair, filled with particles on a Cartesian grid."
        "\nThese are then interacted using runner_dopair1_density."
        "\n\nOptions:"
231
232
233
234
235
        "\n-h DISTANCE=1.1255 - Smoothing length"
	"\n-m rho             - Physical density in the cell"
	"\n-s size            - Physical size of the cell"
        "\n-d pert            - Perturbation to apply to the particles [0,1["
        "\n-f fileName        - Part of the file name used to save the dumps\n",
236
237
238
239
        argv[0]);
    exit(1);
  }

240
241
  /* Build the infrastructure */
  struct space space;
242
243
244
  space.periodic = 0;
  space.h_max = h;

245
  struct engine engine;
246
247
248
  engine.s = &space;
  engine.time = 0.1f;
  engine.ti_current = 1;
249
250

  struct runner runner;
251
252
  runner.e = &engine;

253
254
255
256
257
258
259
  /* Construct some cells */
  struct cell *cells[27];
  struct cell *main_cell;
  static long long partId = 0;
  for (int i = 0; i < 3; ++i) {
    for (int j = 0; j < 3; ++j) {
      for (int k = 0; k < 3; ++k) {
260

261
        double offset[3] = {i * size, j * size, k * size};
262

263
264
265
266
267
268
269
270
271
        cells[i * 9 + j * 3 + k] =
            make_cell(particles, offset, size, h, rho, &partId, perturbation);
      }
    }
  }

  main_cell = cells[13];

  ticks time = 0;
272
273
274
  for (size_t i = 0; i < runs; ++i) {

    /* Zero the fields */
275
    for (int j = 0; j < 27; ++j) zero_particle_fields(cells[j]);
276

277
    const ticks tic = getticks();
278

279
280
281
282
    /* Run all the pairs */
    for (int j = 0; j < 27; ++j)
      if (cells[j] != main_cell)
        runner_dopair1_density(&runner, main_cell, cells[j]);
283

284
285
286
287
    /* And now the self-interaction */
    runner_doself1_density(&runner, main_cell);

    const ticks toc = getticks();
288
289
    time += toc - tic;

290
291
292
    /* Let's get physical ! */
    end_calculation(main_cell);

293
294
    /* Dump if necessary */
    if (i % 50 == 0) {
295
296
297
      sprintf(outputFileName, "swift_dopair_27_%s.dat",
              outputFileNameExtension);
      dump_particle_fields(outputFileName, main_cell, cells);
298
299
300
301
    }
  }

  /* Output timing */
Matthieu Schaller's avatar
Matthieu Schaller committed
302
  message("SWIFT calculation took       : %15lli ticks.", time / runs);
303
304
305
306

  /* Now perform a brute-force version for accuracy tests */

  /* Zero the fields */
307
308
309
310
311
312
313
  for (int i = 0; i < 27; ++i) zero_particle_fields(cells[i]);

  const ticks tic = getticks();

  /* Run all the brute-force pairs */
  for (int j = 0; j < 27; ++j)
    if (cells[j] != main_cell) pairs_all_density(&runner, main_cell, cells[j]);
314

315
316
  /* And now the self-interaction */
  self_all_density(&runner, main_cell);
317

318
  const ticks toc = getticks();
319

320
321
  /* Let's get physical ! */
  end_calculation(main_cell);
322
323
324

  /* Dump */
  sprintf(outputFileName, "brute_force_27_%s.dat", outputFileNameExtension);
325
  dump_particle_fields(outputFileName, main_cell, cells);
326
327

  /* Output timing */
Matthieu Schaller's avatar
Matthieu Schaller committed
328
  message("Brute force calculation took : %15lli ticks.", toc - tic);
329
330

  /* Clean things to make the sanitizer happy ... */
331
  for (int i = 0; i < 27; ++i) clean_up(cells[i]);
332
333
334

  return 0;
}