space.c 151 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
 * This file is part of SWIFT.
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *               2015 Peter W. Draper (p.w.draper@durham.ac.uk)
 *               2016 John A. Regan (john.a.regan@durham.ac.uk)
 *                    Tom Theuns (tom.theuns@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/>.
 *
 ******************************************************************************/
Pedro Gonnet's avatar
Pedro Gonnet committed
23
24
25
26
27
28
29
30

/* Config parameters. */
#include "../config.h"

/* Some standard headers. */
#include <float.h>
#include <limits.h>
#include <math.h>
31
#include <stdlib.h>
32
#include <string.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
33

34
35
/* MPI headers. */
#ifdef WITH_MPI
36
#include <mpi.h>
37
38
#endif

39
40
41
/* This object's header. */
#include "space.h"

Pedro Gonnet's avatar
Pedro Gonnet committed
42
/* Local headers. */
43
#include "atomic.h"
44
#include "chemistry.h"
45
#include "const.h"
46
#include "cooling.h"
47
#include "debug.h"
48
#include "engine.h"
49
#include "error.h"
50
51
#include "gravity.h"
#include "hydro.h"
52
#include "kernel_hydro.h"
53
#include "lock.h"
54
#include "memswap.h"
55
#include "minmax.h"
56
#include "multipole.h"
57
#include "restart.h"
58
#include "sort_part.h"
59
#include "stars.h"
60
#include "threadpool.h"
61
#include "tools.h"
62
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
63
64
65

/* Split size. */
int space_splitsize = space_splitsize_default;
66
67
68
int space_subsize_pair_hydro = space_subsize_pair_hydro_default;
int space_subsize_self_hydro = space_subsize_self_hydro_default;
int space_subsize_pair_grav = space_subsize_pair_grav_default;
69
int space_subsize_self_grav = space_subsize_self_grav_default;
70
int space_subdepth_diff_grav = space_subdepth_diff_grav_default;
71
int space_maxsize = space_maxsize_default;
72

73
74
75
/*! Number of extra #part we allocate memory for per top-level cell */
int space_extra_parts = space_extra_parts_default;

76
77
78
/*! Number of extra #spart we allocate memory for per top-level cell */
int space_extra_sparts = space_extra_sparts_default;

79
80
81
/*! Number of extra #gpart we allocate memory for per top-level cell */
int space_extra_gparts = space_extra_gparts_default;

82
83
/*! Expected maximal number of strays received at a rebuild */
int space_expected_max_nr_strays = space_expected_max_nr_strays_default;
84
85
86
#ifdef SWIFT_DEBUG_CHECKS
int last_cell_id;
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
87

88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
/**
 * @brief Interval stack necessary for parallel particle sorting.
 */
struct qstack {
  volatile ptrdiff_t i, j;
  volatile int min, max;
  volatile int ready;
};

/**
 * @brief Parallel particle-sorting stack
 */
struct parallel_sort {
  struct part *parts;
  struct gpart *gparts;
  struct xpart *xparts;
104
  struct spart *sparts;
105
106
107
108
109
110
  int *ind;
  struct qstack *stack;
  unsigned int stack_size;
  volatile unsigned int first, last, waiting;
};

111
112
113
114
115
116
/**
 * @brief Information required to compute the particle cell indices.
 */
struct index_data {
  struct space *s;
  int *ind;
117
  int *cell_counts;
118
119
120
121
122
123
  size_t count_inhibited_part;
  size_t count_inhibited_gpart;
  size_t count_inhibited_spart;
  size_t count_extra_part;
  size_t count_extra_gpart;
  size_t count_extra_spart;
124
125
};

126
/**
127
 * @brief Recursively dismantle a cell tree.
128
 *
129
130
 * @param s The #space.
 * @param c The #cell to recycle.
Matthieu Schaller's avatar
Matthieu Schaller committed
131
132
133
134
135
136
 * @param cell_rec_begin Pointer to the start of the list of cells to recycle.
 * @param cell_rec_end Pointer to the end of the list of cells to recycle.
 * @param multipole_rec_begin Pointer to the start of the list of multipoles to
 * recycle.
 * @param multipole_rec_end Pointer to the end of the list of multipoles to
 * recycle.
137
 */
138
void space_rebuild_recycle_rec(struct space *s, struct cell *c,
139
140
                               struct cell **cell_rec_begin,
                               struct cell **cell_rec_end,
141
142
                               struct gravity_tensors **multipole_rec_begin,
                               struct gravity_tensors **multipole_rec_end) {
143
  if (c->split)
144
    for (int k = 0; k < 8; k++)
145
      if (c->progeny[k] != NULL) {
146
147
148
149
150
151
        space_rebuild_recycle_rec(s, c->progeny[k], cell_rec_begin,
                                  cell_rec_end, multipole_rec_begin,
                                  multipole_rec_end);

        c->progeny[k]->next = *cell_rec_begin;
        *cell_rec_begin = c->progeny[k];
152

153
        if (s->with_self_gravity) {
154
155
          c->progeny[k]->grav.multipole->next = *multipole_rec_begin;
          *multipole_rec_begin = c->progeny[k]->grav.multipole;
156
        }
157
158

        if (*cell_rec_end == NULL) *cell_rec_end = *cell_rec_begin;
159
        if (s->with_self_gravity && *multipole_rec_end == NULL)
160
161
          *multipole_rec_end = *multipole_rec_begin;

162
        c->progeny[k]->grav.multipole = NULL;
163
164
165
166
        c->progeny[k] = NULL;
      }
}

167
168
169
170
171
172
173
174
void space_rebuild_recycle_mapper(void *map_data, int num_elements,
                                  void *extra_data) {

  struct space *s = (struct space *)extra_data;
  struct cell *cells = (struct cell *)map_data;

  for (int k = 0; k < num_elements; k++) {
    struct cell *c = &cells[k];
175
    struct cell *cell_rec_begin = NULL, *cell_rec_end = NULL;
176
177
    struct gravity_tensors *multipole_rec_begin = NULL,
                           *multipole_rec_end = NULL;
178
179
180
181
182
    space_rebuild_recycle_rec(s, c, &cell_rec_begin, &cell_rec_end,
                              &multipole_rec_begin, &multipole_rec_end);
    if (cell_rec_begin != NULL)
      space_recycle_list(s, cell_rec_begin, cell_rec_end, multipole_rec_begin,
                         multipole_rec_end);
183
    c->hydro.sorts = NULL;
184
    c->stars.sorts = NULL;
185
    c->nr_tasks = 0;
186
187
188
189
    c->grav.nr_mm_tasks = 0;
    c->hydro.density = NULL;
    c->hydro.gradient = NULL;
    c->hydro.force = NULL;
190
    c->hydro.limiter = NULL;
191
192
    c->grav.grav = NULL;
    c->grav.mm = NULL;
193
194
    c->hydro.dx_max_part = 0.0f;
    c->hydro.dx_max_sort = 0.0f;
Loic Hausammann's avatar
Loic Hausammann committed
195
    c->stars.dx_max_part = 0.f;
Loic Hausammann's avatar
Loic Hausammann committed
196
    c->stars.dx_max_sort = 0.f;
197
    c->hydro.sorted = 0;
Loic Hausammann's avatar
Loic Hausammann committed
198
    c->stars.sorted = 0;
199
    c->hydro.count = 0;
200
    c->hydro.count_total = 0;
201
202
    c->hydro.updated = 0;
    c->hydro.inhibited = 0;
203
    c->grav.count = 0;
204
    c->grav.count_total = 0;
205
206
    c->grav.updated = 0;
    c->grav.inhibited = 0;
207
    c->stars.count = 0;
208
    c->stars.count_total = 0;
209
210
    c->stars.updated = 0;
    c->stars.inhibited = 0;
211
212
213
214
215
216
    c->grav.init = NULL;
    c->grav.init_out = NULL;
    c->hydro.extra_ghost = NULL;
    c->hydro.ghost_in = NULL;
    c->hydro.ghost_out = NULL;
    c->hydro.ghost = NULL;
217
218
    c->stars.ghost = NULL;
    c->stars.density = NULL;
Alexei Borissov's avatar
Alexei Borissov committed
219
    c->stars.feedback = NULL;
220
221
    c->kick1 = NULL;
    c->kick2 = NULL;
222
    c->timestep = NULL;
223
    c->timestep_limiter = NULL;
224
    c->hydro.end_force = NULL;
225
    c->hydro.drift = NULL;
226
    c->stars.drift = NULL;
227
228
    c->stars.stars_in = NULL;
    c->stars.stars_out = NULL;
229
    c->grav.drift = NULL;
230
    c->grav.drift_out = NULL;
231
232
233
234
235
    c->hydro.cooling = NULL;
    c->grav.long_range = NULL;
    c->grav.down_in = NULL;
    c->grav.down = NULL;
    c->grav.mesh = NULL;
236
    c->grav.end_force = NULL;
237
    c->super = c;
238
239
240
241
    c->hydro.super = c;
    c->grav.super = c;
    c->hydro.parts = NULL;
    c->hydro.xparts = NULL;
242
243
    c->grav.parts = NULL;
    c->stars.parts = NULL;
244
    c->hydro.do_sub_sort = 0;
Loic Hausammann's avatar
Loic Hausammann committed
245
    c->stars.do_sub_sort = 0;
246
    c->hydro.do_sub_drift = 0;
247
248
    c->grav.do_sub_drift = 0;
    c->stars.do_sub_drift = 0;
249
250
    c->hydro.do_sub_limiter = 0;
    c->hydro.do_limiter = 0;
251
252
253
254
    c->hydro.ti_end_min = -1;
    c->hydro.ti_end_max = -1;
    c->grav.ti_end_min = -1;
    c->grav.ti_end_max = -1;
Loic Hausammann's avatar
Loic Hausammann committed
255
    c->stars.ti_end_min = -1;
256
    c->stars.ti_end_max = -1;
257
258
259
#ifdef SWIFT_DEBUG_CHECKS
    c->cellID = 0;
#endif
260
261
    if (s->with_self_gravity)
      bzero(c->grav.multipole, sizeof(struct gravity_tensors));
Loic Hausammann's avatar
Loic Hausammann committed
262
    for (int i = 0; i < 13; i++) {
263
264
265
      if (c->hydro.sort[i] != NULL) {
        free(c->hydro.sort[i]);
        c->hydro.sort[i] = NULL;
266
      }
Loic Hausammann's avatar
Loic Hausammann committed
267
268
269
270
271
      if (c->stars.sort[i] != NULL) {
        free(c->stars.sort[i]);
        c->stars.sort[i] = NULL;
      }
    }
272
#if WITH_MPI
273
274
    c->mpi.tag = -1;

275
276
277
278
    c->mpi.hydro.recv_xv = NULL;
    c->mpi.hydro.recv_rho = NULL;
    c->mpi.hydro.recv_gradient = NULL;
    c->mpi.grav.recv = NULL;
279
    c->mpi.stars.recv = NULL;
280
    c->mpi.recv_ti = NULL;
281
    c->mpi.limiter.recv = NULL;
282

283
284
285
286
    c->mpi.hydro.send_xv = NULL;
    c->mpi.hydro.send_rho = NULL;
    c->mpi.hydro.send_gradient = NULL;
    c->mpi.grav.send = NULL;
287
    c->mpi.stars.send = NULL;
288
    c->mpi.send_ti = NULL;
289
    c->mpi.limiter.send = NULL;
290
291
292
293
#endif
  }
}

294
295
296
297
/**
 * @brief Free up any allocated cells.
 */
void space_free_cells(struct space *s) {
298
299
300

  ticks tic = getticks();

Matthieu Schaller's avatar
Matthieu Schaller committed
301
302
  threadpool_map(&s->e->threadpool, space_rebuild_recycle_mapper, s->cells_top,
                 s->nr_cells, sizeof(struct cell), 0, s);
303
  s->maxdepth = 0;
304
305
306
307

  if (s->e->verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
308
309
}

310
/**
311
 * @brief Re-build the top-level cell grid.
312
 *
313
 * @param s The #space.
314
 * @param verbose Print messages to stdout or not.
315
 */
316
void space_regrid(struct space *s, int verbose) {
317

318
  const size_t nr_parts = s->nr_parts;
Loic Hausammann's avatar
Loic Hausammann committed
319
  const size_t nr_sparts = s->nr_sparts;
320
  const ticks tic = getticks();
321
  const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0;
322

323
  /* Run through the cells and get the current h_max. */
324
  // tic = getticks();
325
  float h_max = s->cell_min / kernel_gamma / space_stretch;
326
  if (nr_parts > 0) {
327
328
329
330
331
332

    /* Can we use the list of local non-empty top-level cells? */
    if (s->local_cells_with_particles_top != NULL) {
      for (int k = 0; k < s->nr_local_cells_with_particles; ++k) {
        const struct cell *c =
            &s->cells_top[s->local_cells_with_particles_top[k]];
333
        if (c->hydro.h_max > h_max) {
334
          h_max = c->hydro.h_max;
335
        }
Loic Hausammann's avatar
Loic Hausammann committed
336
        if (c->stars.h_max > h_max) {
337
          h_max = c->stars.h_max;
Loic Hausammann's avatar
Loic Hausammann committed
338
        }
339
      }
340
341

      /* Can we instead use all the top-level cells? */
342
    } else if (s->cells_top != NULL) {
Tom Theuns's avatar
Tom Theuns committed
343
      for (int k = 0; k < s->nr_cells; k++) {
344
        const struct cell *c = &s->cells_top[k];
345
        if (c->nodeID == engine_rank && c->hydro.h_max > h_max) {
346
          h_max = c->hydro.h_max;
347
        }
Loic Hausammann's avatar
Loic Hausammann committed
348
        if (c->nodeID == engine_rank && c->stars.h_max > h_max) {
349
          h_max = c->stars.h_max;
Loic Hausammann's avatar
Loic Hausammann committed
350
        }
351
      }
352
353

      /* Last option: run through the particles */
354
    } else {
355
      for (size_t k = 0; k < nr_parts; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
356
        if (s->parts[k].h > h_max) h_max = s->parts[k].h;
357
      }
Loic Hausammann's avatar
Loic Hausammann committed
358
359
360
      for (size_t k = 0; k < nr_sparts; k++) {
        if (s->sparts[k].h > h_max) h_max = s->sparts[k].h;
      }
361
362
363
364
365
366
367
368
369
370
    }
  }

/* If we are running in parallel, make sure everybody agrees on
   how large the largest cell should be. */
#ifdef WITH_MPI
  {
    float buff;
    if (MPI_Allreduce(&h_max, &buff, 1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD) !=
        MPI_SUCCESS)
371
      error("Failed to aggregate the rebuild flag across nodes.");
372
373
374
    h_max = buff;
  }
#endif
375
  if (verbose) message("h_max is %.3e (cell_min=%.3e).", h_max, s->cell_min);
376
377

  /* Get the new putative cell dimensions. */
378
  const int cdim[3] = {
379
380
381
382
383
384
      (int)floor(s->dim[0] /
                 fmax(h_max * kernel_gamma * space_stretch, s->cell_min)),
      (int)floor(s->dim[1] /
                 fmax(h_max * kernel_gamma * space_stretch, s->cell_min)),
      (int)floor(s->dim[2] /
                 fmax(h_max * kernel_gamma * space_stretch, s->cell_min))};
385
386
387
388
389

  /* Check if we have enough cells for periodicity. */
  if (s->periodic && (cdim[0] < 3 || cdim[1] < 3 || cdim[2] < 3))
    error(
        "Must have at least 3 cells in each spatial dimension when periodicity "
390
391
392
        "is switched on.\nThis error is often caused by any of the "
        "followings:\n"
        " - too few particles to generate a sensible grid,\n"
393
394
        " - the initial value of 'Scheduler:max_top_level_cells' is too "
        "small,\n"
395
        " - the (minimal) time-step is too large leading to particles with "
396
        "predicted smoothing lengths too large for the box size,\n"
397
        " - particles with velocities so large that they move by more than two "
398
        "box sizes per time-step.\n");
399

400
401
402
/* In MPI-Land, changing the top-level cell size requires that the
 * global partition is recomputed and the particles redistributed.
 * Be prepared to do that. */
403
#ifdef WITH_MPI
Matthieu Schaller's avatar
Matthieu Schaller committed
404
  double oldwidth[3];
405
406
407
408
409
410
411
412
  double oldcdim[3];
  int *oldnodeIDs = NULL;
  if (cdim[0] < s->cdim[0] || cdim[1] < s->cdim[1] || cdim[2] < s->cdim[2]) {

    /* Capture state of current space. */
    oldcdim[0] = s->cdim[0];
    oldcdim[1] = s->cdim[1];
    oldcdim[2] = s->cdim[2];
413
414
415
    oldwidth[0] = s->width[0];
    oldwidth[1] = s->width[1];
    oldwidth[2] = s->width[2];
416
417
418
419
420
421
422
423
424

    if ((oldnodeIDs = (int *)malloc(sizeof(int) * s->nr_cells)) == NULL)
      error("Failed to allocate temporary nodeIDs.");

    int cid = 0;
    for (int i = 0; i < s->cdim[0]; i++) {
      for (int j = 0; j < s->cdim[1]; j++) {
        for (int k = 0; k < s->cdim[2]; k++) {
          cid = cell_getid(oldcdim, i, j, k);
425
          oldnodeIDs[cid] = s->cells_top[cid].nodeID;
426
427
428
429
430
        }
      }
    }
  }

Peter W. Draper's avatar
Peter W. Draper committed
431
  /* Are we about to allocate new top level cells without a regrid?
Peter W. Draper's avatar
Peter W. Draper committed
432
   * Can happen when restarting the application. */
433
  const int no_regrid = (s->cells_top == NULL && oldnodeIDs == NULL);
434
435
436
437
#endif

  /* Do we need to re-build the upper-level cells? */
  // tic = getticks();
438
  if (s->cells_top == NULL || cdim[0] < s->cdim[0] || cdim[1] < s->cdim[1] ||
439
440
      cdim[2] < s->cdim[2]) {

441
442
/* Be verbose about this. */
#ifdef SWIFT_DEBUG_CHECKS
443
    message("(re)griding space cdim=(%d %d %d)", cdim[0], cdim[1], cdim[2]);
444
445
446
    fflush(stdout);
#endif

447
    /* Free the old cells, if they were allocated. */
448
    if (s->cells_top != NULL) {
449
      space_free_cells(s);
450
      free(s->local_cells_with_tasks_top);
451
      free(s->local_cells_top);
452
      free(s->cells_with_particles_top);
453
      free(s->local_cells_with_particles_top);
454
      free(s->cells_top);
455
      free(s->multipoles_top);
456
457
    }

458
459
460
461
    /* Also free the task arrays, these will be regenerated and we can use the
     * memory while copying the particle arrays. */
    if (s->e != NULL) scheduler_free_tasks(&s->e->sched);

462
    /* Set the new cell dimensions only if smaller. */
463
    for (int k = 0; k < 3; k++) {
464
      s->cdim[k] = cdim[k];
465
466
      s->width[k] = s->dim[k] / cdim[k];
      s->iwidth[k] = 1.0 / s->width[k];
467
    }
468
    const float dmin = min3(s->width[0], s->width[1], s->width[2]);
469
470
471

    /* Allocate the highest level of cells. */
    s->tot_cells = s->nr_cells = cdim[0] * cdim[1] * cdim[2];
472
    if (posix_memalign((void **)&s->cells_top, cell_align,
473
                       s->nr_cells * sizeof(struct cell)) != 0)
474
      error("Failed to allocate top-level cells.");
475
    bzero(s->cells_top, s->nr_cells * sizeof(struct cell));
476

477
    /* Allocate the multipoles for the top-level cells. */
478
    if (s->with_self_gravity) {
479
      if (posix_memalign((void **)&s->multipoles_top, multipole_align,
480
                         s->nr_cells * sizeof(struct gravity_tensors)) != 0)
481
        error("Failed to allocate top-level multipoles.");
482
      bzero(s->multipoles_top, s->nr_cells * sizeof(struct gravity_tensors));
483
484
    }

485
    /* Allocate the indices of local cells */
486
    if (posix_memalign((void **)&s->local_cells_top, SWIFT_STRUCT_ALIGNMENT,
487
488
489
490
                       s->nr_cells * sizeof(int)) != 0)
      error("Failed to allocate indices of local top-level cells.");
    bzero(s->local_cells_top, s->nr_cells * sizeof(int));

491
    /* Allocate the indices of local cells with tasks */
492
493
    if (posix_memalign((void **)&s->local_cells_with_tasks_top,
                       SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0)
494
      error("Failed to allocate indices of local top-level cells with tasks.");
495
496
    bzero(s->local_cells_with_tasks_top, s->nr_cells * sizeof(int));

497
    /* Allocate the indices of cells with particles */
498
    if (posix_memalign((void **)&s->cells_with_particles_top,
499
                       SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0)
500
501
      error("Failed to allocate indices of top-level cells with particles.");
    bzero(s->cells_with_particles_top, s->nr_cells * sizeof(int));
502

503
504
505
506
507
508
509
510
    /* Allocate the indices of local cells with particles */
    if (posix_memalign((void **)&s->local_cells_with_particles_top,
                       SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0)
      error(
          "Failed to allocate indices of local top-level cells with "
          "particles.");
    bzero(s->local_cells_with_particles_top, s->nr_cells * sizeof(int));

511
    /* Set the cells' locks */
512
    for (int k = 0; k < s->nr_cells; k++) {
513
      if (lock_init(&s->cells_top[k].hydro.lock) != 0)
514
        error("Failed to init spinlock for hydro.");
515
      if (lock_init(&s->cells_top[k].grav.plock) != 0)
516
        error("Failed to init spinlock for gravity.");
517
      if (lock_init(&s->cells_top[k].grav.mlock) != 0)
518
        error("Failed to init spinlock for multipoles.");
519
      if (lock_init(&s->cells_top[k].stars.lock) != 0)
520
521
        error("Failed to init spinlock for stars.");
    }
522
523

    /* Set the cell location and sizes. */
524
525
526
    for (int i = 0; i < cdim[0]; i++)
      for (int j = 0; j < cdim[1]; j++)
        for (int k = 0; k < cdim[2]; k++) {
527
528
          const size_t cid = cell_getid(cdim, i, j, k);
          struct cell *restrict c = &s->cells_top[cid];
529
530
531
532
533
534
          c->loc[0] = i * s->width[0];
          c->loc[1] = j * s->width[1];
          c->loc[2] = k * s->width[2];
          c->width[0] = s->width[0];
          c->width[1] = s->width[1];
          c->width[2] = s->width[2];
535
536
          c->dmin = dmin;
          c->depth = 0;
537
          c->split = 0;
538
          c->hydro.count = 0;
539
540
          c->grav.count = 0;
          c->stars.count = 0;
541
          c->super = c;
542
543
          c->hydro.super = c;
          c->grav.super = c;
544
545
          c->hydro.ti_old_part = ti_current;
          c->grav.ti_old_part = ti_current;
546
          c->stars.ti_old_part = ti_current;
547
          c->grav.ti_old_multipole = ti_current;
Pedro Gonnet's avatar
Pedro Gonnet committed
548
#ifdef WITH_MPI
549
          c->mpi.tag = -1;
550
551
552
553
554
555
          c->mpi.hydro.recv_xv = NULL;
          c->mpi.hydro.recv_rho = NULL;
          c->mpi.hydro.recv_gradient = NULL;
          c->mpi.hydro.send_xv = NULL;
          c->mpi.hydro.send_rho = NULL;
          c->mpi.hydro.send_gradient = NULL;
556
557
          c->mpi.stars.send = NULL;
          c->mpi.stars.recv = NULL;
558
559
          c->mpi.grav.recv = NULL;
          c->mpi.grav.send = NULL;
Pedro Gonnet's avatar
Pedro Gonnet committed
560
#endif  // WITH_MPI
561
          if (s->with_self_gravity) c->grav.multipole = &s->multipoles_top[cid];
562
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
563
564
          c->cellID = -last_cell_id;
          last_cell_id++;
565
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
566
        }
567
568

    /* Be verbose about the change. */
569
570
571
    if (verbose)
      message("set cell dimensions to [ %i %i %i ].", cdim[0], cdim[1],
              cdim[2]);
572

573
#ifdef WITH_MPI
574
575
576
577
578
    if (oldnodeIDs != NULL) {
      /* We have changed the top-level cell dimension, so need to redistribute
       * cells around the nodes. We repartition using the old space node
       * positions as a grid to resample. */
      if (s->e->nodeID == 0)
579
580
581
        message(
            "basic cell dimensions have increased - recalculating the "
            "global partition.");
582

Matthieu Schaller's avatar
Matthieu Schaller committed
583
      if (!partition_space_to_space(oldwidth, oldcdim, oldnodeIDs, s)) {
584
585
586
587

        /* Failed, try another technique that requires no settings. */
        message("Failed to get a new partition, trying less optimal method");
        struct partition initial_partition;
588
#if defined(HAVE_PARMETIS) || defined(HAVE_METIS)
589
590
591
592
593
594
595
596
597
598
599
600
601
        initial_partition.type = INITPART_METIS_NOWEIGHT;
#else
        initial_partition.type = INITPART_VECTORIZE;
#endif
        partition_initial_partition(&initial_partition, s->e->nodeID,
                                    s->e->nr_nodes, s);
      }

      /* Re-distribute the particles to their new nodes. */
      engine_redistribute(s->e);

      /* Make the proxies. */
      engine_makeproxies(s->e);
602

603
604
      /* Finished with these. */
      free(oldnodeIDs);
Peter W. Draper's avatar
Peter W. Draper committed
605
606

    } else if (no_regrid && s->e != NULL) {
Peter W. Draper's avatar
Peter W. Draper committed
607
608
609
610
611
612
613
614
615
616
617
618
619
      /* If we have created the top-levels cells and not done an initial
       * partition (can happen when restarting), then the top-level cells
       * are not assigned to a node, we must do that and then associate the
       * particles with the cells. Note requires that
       * partition_store_celllist() was called once before, or just before
       * dumping the restart files.*/
      partition_restore_celllist(s, s->e->reparttype);

      /* Now re-distribute the particles, should just add to cells? */
      engine_redistribute(s->e);

      /* Make the proxies. */
      engine_makeproxies(s->e);
620
    }
Pedro Gonnet's avatar
Pedro Gonnet committed
621
#endif /* WITH_MPI */
622
623
624
625

    // message( "rebuilding upper-level cells took %.3f %s." ,
    // clocks_from_ticks(double)(getticks() - tic), clocks_getunit());

626
  }      /* re-build upper-level cells? */
627
  else { /* Otherwise, just clean up the cells. */
628
629

    /* Free the old cells, if they were allocated. */
630
    space_free_cells(s);
631
  }
632
633
634
635

  if (verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
636
}
637

638
639
640
641
642
643
/**
 * @brief Allocate memory for the extra particles used for on-the-fly creation.
 *
 * This rarely actually allocates memory. Most of the time, we convert
 * pre-allocated memory inot extra particles.
 *
644
645
 * This function also sets the extra particles' location to their top-level
 * cells. They can then be sorted into their correct memory position later on.
646
647
648
649
 *
 * @param s The current #space.
 * @param verbose Are we talkative?
 */
650
651
652
653
void space_allocate_extras(struct space *s, int verbose) {

  const int local_nodeID = s->e->nodeID;

654
655
656
657
658
  /* Anything to do here? (Abort if we don't want extras)*/
  if (space_extra_parts == 0 && space_extra_gparts == 0 &&
      space_extra_sparts == 0)
    return;

659
660
661
662
663
664
  /* The top-level cells */
  const struct cell *cells = s->cells_top;
  const double half_cell_width[3] = {0.5 * cells[0].width[0],
                                     0.5 * cells[0].width[1],
                                     0.5 * cells[0].width[2]};

665
  /* The current number of particles (including spare ones) */
666
667
668
669
  size_t nr_parts = s->nr_parts;
  size_t nr_gparts = s->nr_gparts;
  size_t nr_sparts = s->nr_sparts;

670
671
672
673
674
  /* The current number of actual particles */
  size_t nr_actual_parts = nr_parts - s->nr_extra_parts;
  size_t nr_actual_gparts = nr_gparts - s->nr_extra_gparts;
  size_t nr_actual_sparts = nr_sparts - s->nr_extra_sparts;

675
676
677
678
679
680
681
682
683
  /* The number of particles we allocated memory for (MPI overhead) */
  size_t size_parts = s->size_parts;
  size_t size_gparts = s->size_gparts;
  size_t size_sparts = s->size_sparts;

  int local_cells = 0;
  for (int i = 0; i < s->nr_cells; ++i)
    if (s->cells_top[i].nodeID == local_nodeID) local_cells++;

684
685
  /* Number of extra particles we want for each type */
  const size_t expected_num_extra_parts = local_cells * space_extra_parts;
686
687
  const size_t expected_num_extra_gparts = local_cells * space_extra_gparts;
  const size_t expected_num_extra_sparts = local_cells * space_extra_sparts;
688

689
690
691
  if (verbose) {
    message("Currently have %zd/%zd/%zd real particles.", nr_actual_parts,
            nr_actual_gparts, nr_actual_sparts);
692
    message("Currently have %zd/%zd/%zd spaces for extra particles.",
693
            s->nr_extra_parts, s->nr_extra_gparts, s->nr_extra_sparts);
694
    message("Requesting space for future %zd/%zd/%zd part/gpart/sparts.",
695
696
697
            expected_num_extra_parts, expected_num_extra_gparts,
            expected_num_extra_sparts);
  }
698

699
700
701
702
703
704
705
  if (expected_num_extra_parts < s->nr_extra_parts)
    error("Reduction in top-level cells number not handled.");
  if (expected_num_extra_gparts < s->nr_extra_gparts)
    error("Reduction in top-level cells number not handled.");
  if (expected_num_extra_sparts < s->nr_extra_sparts)
    error("Reduction in top-level cells number not handled.");

706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
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
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
  /* Do we have enough space for the extra gparts (i.e. we haven't used up any)
   * ? */
  if (nr_gparts + expected_num_extra_gparts > size_gparts) {

    /* Ok... need to put some more in the game */

    /* Do we need to reallocate? */
    if (nr_actual_gparts + expected_num_extra_gparts > size_gparts) {

      size_gparts = (nr_actual_gparts + expected_num_extra_gparts) *
                    engine_redistribute_alloc_margin;

      if (verbose)
        message("Re-allocating gparts array from %zd to %zd", s->size_gparts,
                size_gparts);

      /* Create more space for parts */
      struct gpart *gparts_new = NULL;
      if (posix_memalign((void **)&gparts_new, gpart_align,
                         sizeof(struct gpart) * size_gparts) != 0)
        error("Failed to allocate new gpart data");
      const ptrdiff_t delta = gparts_new - s->gparts;
      memcpy(gparts_new, s->gparts, sizeof(struct gpart) * s->size_gparts);
      free(s->gparts);
      s->gparts = gparts_new;

      /* Update the counter */
      s->size_gparts = size_gparts;

      /* We now need to reset all the part and spart pointers */
      for (size_t i = 0; i < nr_parts; ++i) {
        if (s->parts[i].time_bin != time_bin_not_created)
          s->parts[i].gpart += delta;
      }
      for (size_t i = 0; i < nr_sparts; ++i) {
        if (s->sparts[i].time_bin != time_bin_not_created)
          s->sparts[i].gpart += delta;
      }
    }

    /* Turn some of the allocated spares into particles we can use */
    for (size_t i = nr_gparts; i < nr_actual_gparts + expected_num_extra_gparts;
         ++i) {
      bzero(&s->gparts[i], sizeof(struct gpart));
      s->gparts[i].time_bin = time_bin_not_created;
      s->gparts[i].type = swift_type_dark_matter;
      s->gparts[i].id_or_neg_offset = -1;
    }

      /* Put the spare particles in their correct cell */
#ifdef WITH_MPI
    error("Need to do this correctly over MPI for only the local cells.");
#endif
    int count_in_cell = 0, current_cell = 0;
    size_t count_extra_gparts = 0;
    for (size_t i = 0; i < nr_actual_gparts + expected_num_extra_gparts; ++i) {

#ifdef SWIFT_DEBUG_CHECKS
      if (current_cell == s->nr_cells)
        error("Cell counter beyond the maximal nr. cells.");
#endif

      if (s->gparts[i].time_bin == time_bin_not_created) {

        /* We want the extra particles to be at the centre of their cell */
        s->gparts[i].x[0] = cells[current_cell].loc[0] + half_cell_width[0];
        s->gparts[i].x[1] = cells[current_cell].loc[1] + half_cell_width[1];
        s->gparts[i].x[2] = cells[current_cell].loc[2] + half_cell_width[2];
        ++count_in_cell;
        count_extra_gparts++;
      }

      /* Once we have reached the number of extra gpart per cell, we move to the
       * next */
      if (count_in_cell == space_extra_gparts) {
        ++current_cell;
        count_in_cell = 0;
      }
    }

#ifdef SWIFT_DEBUG_CHECKS
    if (count_extra_gparts != expected_num_extra_gparts)
      error("Constructed the wrong number of extra gparts (%zd vs. %zd)",
            count_extra_gparts, expected_num_extra_gparts);
#endif

    /* Update the counters */
    s->nr_gparts = nr_actual_gparts + expected_num_extra_gparts;
    s->nr_extra_gparts = expected_num_extra_gparts;
  }

797
798
  /* Do we have enough space for the extra parts (i.e. we haven't used up any) ?
   */
799
800
  if (expected_num_extra_parts > s->nr_extra_parts) {

801
802
    /* Ok... need to put some more in the game */

803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
    /* Do we need to reallocate? */
    if (nr_actual_parts + expected_num_extra_parts > size_parts) {

      size_parts = (nr_actual_parts + expected_num_extra_parts) *
                   engine_redistribute_alloc_margin;

      if (verbose)
        message("Re-allocating parts array from %zd to %zd", s->size_parts,
                size_parts);

      /* Create more space for parts */
      struct part *parts_new = NULL;
      if (posix_memalign((void **)&parts_new, part_align,
                         sizeof(struct part) * size_parts) != 0)
        error("Failed to allocate new part data");
      memcpy(parts_new, s->parts, sizeof(struct part) * s->size_parts);
      free(s->parts);
      s->parts = parts_new;

      /* Same for xparts */
      struct xpart *xparts_new = NULL;
      if (posix_memalign((void **)&xparts_new, xpart_align,
                         sizeof(struct xpart) * size_parts) != 0)
        error("Failed to allocate new xpart data");
      memcpy(xparts_new, s->xparts, sizeof(struct xpart) * s->size_parts);
      free(s->xparts);
      s->xparts = xparts_new;

      /* Update the counter */
      s->size_parts = size_parts;
    }

835
    /* Turn some of the allocated spares into particles we can use */
836
837
838
839
840
    for (size_t i = nr_parts; i < nr_actual_parts + expected_num_extra_parts;
         ++i) {
      bzero(&s->parts[i], sizeof(struct part));
      bzero(&s->xparts[i], sizeof(struct xpart));
      s->parts[i].time_bin = time_bin_not_created;
841
      s->parts[i].id = -1;
842
843
    }

844
845
846
847
848
849
850
      /* Put the spare particles in their correct cell */
#ifdef WITH_MPI
    error("Need to do this correctly over MPI for only the local cells.");
#endif
    int count_in_cell = 0, current_cell = 0;
    size_t count_extra_parts = 0;
    for (size_t i = 0; i < nr_actual_parts + expected_num_extra_parts; ++i) {
851
852
853
854
855
856

#ifdef SWIFT_DEBUG_CHECKS
      if (current_cell == s->nr_cells)
        error("Cell counter beyond the maximal nr. cells.");
#endif

857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
      if (s->parts[i].time_bin == time_bin_not_created) {

        /* We want the extra particles to be at the centre of their cell */
        s->parts[i].x[0] = cells[current_cell].loc[0] + half_cell_width[0];
        s->parts[i].x[1] = cells[current_cell].loc[1] + half_cell_width[1];
        s->parts[i].x[2] = cells[current_cell].loc[2] + half_cell_width[2];
        ++count_in_cell;
        count_extra_parts++;
      }

      /* Once we have reached the number of extra part per cell, we move to the
       * next */
      if (count_in_cell == space_extra_parts) {
        ++current_cell;
        count_in_cell = 0;
      }
    }

#ifdef SWIFT_DEBUG_CHECKS
    if (count_extra_parts != expected_num_extra_parts)
      error("Constructed the wrong number of extra parts (%zd vs. %zd)",
            count_extra_parts, expected_num_extra_parts);
#endif

881
882
883
    /* Update the counters */
    s->nr_parts = nr_actual_parts + expected_num_extra_parts;
    s->nr_extra_parts = expected_num_extra_parts;
884
  }
885

886
887
  /* Do we have enough space for the extra sparts (i.e. we haven't used up any)
   * ? */
888
  if (nr_actual_sparts + expected_num_extra_sparts > nr_sparts) {
889

890
891
    /* Ok... need to put some more in the game */

892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
    /* Do we need to reallocate? */
    if (nr_actual_sparts + expected_num_extra_sparts > size_sparts) {

      size_sparts = (nr_actual_sparts + expected_num_extra_sparts) *
                    engine_redistribute_alloc_margin;

      if (verbose)
        message("Re-allocating sparts array from %zd to %zd", s->size_sparts,
                size_sparts);

      /* Create more space for parts */
      struct spart *sparts_new = NULL;
      if (posix_memalign((void **)&sparts_new, spart_align,
                         sizeof(struct spart) * size_sparts) != 0)
        error("Failed to allocate new spart data");
      memcpy(sparts_new, s->sparts, sizeof(struct spart) * s->size_sparts);
      free(s->sparts);
      s->sparts = sparts_new;

      /* Update the counter */
      s->size_sparts = size_sparts;
    }

    /* Turn some of the allocated spares into particles we can use */
    for (size_t i = nr_sparts; i < nr_actual_sparts + expected_num_extra_sparts;
         ++i) {
      bzero(&s->sparts[i], sizeof(struct spart));
      s->sparts[i].time_bin = time_bin_not_created;
920
      s->sparts[i].id = -42;
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
    }

      /* Put the spare particles in their correct cell */
#ifdef WITH_MPI
    error("Need to do this correctly over MPI for only the local cells.");
#endif
    int count_in_cell = 0, current_cell = 0;
    size_t count_extra_sparts = 0;
    for (size_t i = 0; i < nr_actual_sparts + expected_num_extra_sparts; ++i) {

#ifdef SWIFT_DEBUG_CHECKS
      if (current_cell == s->nr_cells)
        error("Cell counter beyond the maximal nr. cells.");
#endif

      if (s->sparts[i].time_bin == time_bin_not_created) {

        /* We want the extra particles to be at the centre of their cell */
        s->sparts[i].x[0] = cells[current_cell].loc[0] + half_cell_width[0];
        s->sparts[i].x[1] = cells[current_cell].loc[1] + half_cell_width[1];
        s->sparts[i].x[2] = cells[current_cell].loc[2] + half_cell_width[2];
        ++count_in_cell;
        count_extra_sparts++;
      }

      /* Once we have reached the number of extra spart per cell, we move to the
       * next */
      if (count_in_cell == space_extra_sparts) {
        ++current_cell;
        count_in_cell = 0;
      }
    }

#ifdef SWIFT_DEBUG_CHECKS
    if (count_extra_sparts != expected_num_extra_sparts)
      error("Constructed the wrong number of extra sparts (%zd vs. %zd)",
            count_extra_sparts, expected_num_extra_sparts);
#endif

    /* Update the counters */
    s->nr_sparts = nr_actual_sparts + expected_num_extra_sparts;
    s->nr_extra_sparts = expected_num_extra_sparts;
963
  }
964
965
966
967
968
969
970

#ifdef SWIFT_DEBUG_CHECKS
  /* Verify that the links are correct */
  if ((nr_gparts > 0 && nr_parts > 0) || (nr_gparts > 0 && nr_sparts > 0))
    part_verify_links(s->parts, s->gparts, s->sparts, nr_parts, nr_gparts,
                      nr_sparts, verbose);
#endif
971
972
}

973
974
975
976
/**
 * @brief Re-build the cells as well as the tasks.
 *
 * @param s The #space in which to update the cells.
977
 * @param repartitioned Did we just repartition?
978
 * @param verbose Print messages to stdout or not
979
 */
980
void space_rebuild(struct space *s, int repartitioned, int verbose) {
981

Matthieu Schaller's avatar
Matthieu Schaller committed
982
  const ticks tic = getticks();
983

984
985
/* Be verbose about this. */
#ifdef SWIFT_DEBUG_CHECKS
Pedro Gonnet's avatar
Pedro Gonnet committed
986
  if (s->e->nodeID == 0 || verbose) message("(re)building space");
987
988
  fflush(stdout);
#endif
989
990

  /* Re-grid if necessary, or just re-set the cell data. */
991
  space_regrid(s, verbose);
992

993
  /* Allocate extra space for particles that will be created */
994
  if (s->with_star_formation) space_allocate_extras(s, verbose);
995

996
997
  struct cell *cells_top = s->cells_top;
  const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0;
998
  const int local_nodeID = s->e->nodeID;
999
1000

  /* The current number of particles */
For faster browsing, not all history is shown. View entire blame