space.c 124 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"
Pedro Gonnet's avatar
Pedro Gonnet committed
62
63
64

/* Split size. */
int space_splitsize = space_splitsize_default;
65
66
67
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;
68
int space_subsize_self_grav = space_subsize_self_grav_default;
69
70
int space_subsize_pair_stars = space_subsize_pair_stars_default;
int space_subsize_self_stars = space_subsize_self_stars_default;
71
int space_subdepth_grav = space_subdepth_grav_default;
72
int space_maxsize = space_maxsize_default;
73
74
75
#ifdef SWIFT_DEBUG_CHECKS
int last_cell_id;
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
76

77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
/**
 * @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;
93
  struct spart *sparts;
94
95
96
97
98
99
  int *ind;
  struct qstack *stack;
  unsigned int stack_size;
  volatile unsigned int first, last, waiting;
};

100
101
102
103
104
105
/**
 * @brief Information required to compute the particle cell indices.
 */
struct index_data {
  struct space *s;
  int *ind;
106
  int *cell_counts;
107
108
109
  int count_inhibited_part;
  int count_inhibited_gpart;
  int count_inhibited_spart;
110
111
};

112
/**
113
 * @brief Recursively dismantle a cell tree.
114
 *
115
116
 * @param s The #space.
 * @param c The #cell to recycle.
Matthieu Schaller's avatar
Matthieu Schaller committed
117
118
119
120
121
122
 * @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.
123
 */
124
void space_rebuild_recycle_rec(struct space *s, struct cell *c,
125
126
                               struct cell **cell_rec_begin,
                               struct cell **cell_rec_end,
127
128
                               struct gravity_tensors **multipole_rec_begin,
                               struct gravity_tensors **multipole_rec_end) {
129
  if (c->split)
130
    for (int k = 0; k < 8; k++)
131
      if (c->progeny[k] != NULL) {
132
133
134
135
136
137
        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];
138
139

        if (s->gravity) {
140
141
          c->progeny[k]->grav.multipole->next = *multipole_rec_begin;
          *multipole_rec_begin = c->progeny[k]->grav.multipole;
142
        }
143
144

        if (*cell_rec_end == NULL) *cell_rec_end = *cell_rec_begin;
145
        if (s->gravity && *multipole_rec_end == NULL)
146
147
          *multipole_rec_end = *multipole_rec_begin;

148
        c->progeny[k]->grav.multipole = NULL;
149
150
151
152
        c->progeny[k] = NULL;
      }
}

153
154
155
156
157
158
159
160
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];
161
    struct cell *cell_rec_begin = NULL, *cell_rec_end = NULL;
162
163
    struct gravity_tensors *multipole_rec_begin = NULL,
                           *multipole_rec_end = NULL;
164
165
166
167
168
    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);
169
    c->hydro.sorts = NULL;
170
    c->nr_tasks = 0;
171
172
173
174
175
176
    c->grav.nr_mm_tasks = 0;
    c->hydro.density = NULL;
    c->hydro.gradient = NULL;
    c->hydro.force = NULL;
    c->grav.grav = NULL;
    c->grav.mm = NULL;
177
178
    c->hydro.dx_max_part = 0.0f;
    c->hydro.dx_max_sort = 0.0f;
Loic Hausammann's avatar
Loic Hausammann committed
179
    c->stars.dx_max_part = 0.f;
180
181
    c->hydro.sorted = 0;
    c->hydro.count = 0;
182
183
    c->hydro.updated = 0;
    c->hydro.inhibited = 0;
184
    c->grav.count = 0;
185
186
    c->grav.updated = 0;
    c->grav.inhibited = 0;
187
    c->stars.count = 0;
188
189
    c->stars.updated = 0;
    c->stars.inhibited = 0;
190
191
192
193
194
195
    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;
196
197
198
199
    c->stars.ghost_in = NULL;
    c->stars.ghost_out = NULL;
    c->stars.ghost = NULL;
    c->stars.density = NULL;
200
201
    c->kick1 = NULL;
    c->kick2 = NULL;
202
    c->timestep = NULL;
203
    c->end_force = NULL;
204
    c->hydro.drift = NULL;
205
    c->grav.drift = NULL;
206
    c->hydro.cooling = NULL;
207
    c->sourceterms = NULL;
208
209
210
211
    c->grav.long_range = NULL;
    c->grav.down_in = NULL;
    c->grav.down = NULL;
    c->grav.mesh = NULL;
212
    c->super = c;
213
214
215
216
    c->hydro.super = c;
    c->grav.super = c;
    c->hydro.parts = NULL;
    c->hydro.xparts = NULL;
217
218
    c->grav.parts = NULL;
    c->stars.parts = NULL;
219
220
221
222
223
224
225
    c->hydro.do_sub_sort = 0;
    c->grav.do_sub_drift = 0;
    c->hydro.do_sub_drift = 0;
    c->hydro.ti_end_min = -1;
    c->hydro.ti_end_max = -1;
    c->grav.ti_end_min = -1;
    c->grav.ti_end_max = -1;
226
227
228
#ifdef SWIFT_DEBUG_CHECKS
    c->cellID = 0;
#endif
229
    if (s->gravity) bzero(c->grav.multipole, sizeof(struct gravity_tensors));
230
    for (int i = 0; i < 13; i++)
231
232
233
      if (c->hydro.sort[i] != NULL) {
        free(c->hydro.sort[i]);
        c->hydro.sort[i] = NULL;
234
      }
235
#if WITH_MPI
236
237
    c->mpi.tag = -1;

238
239
240
241
    c->mpi.hydro.recv_xv = NULL;
    c->mpi.hydro.recv_rho = NULL;
    c->mpi.hydro.recv_gradient = NULL;
    c->mpi.grav.recv = NULL;
242
243
    c->mpi.recv_ti = NULL;

244
245
246
247
    c->mpi.hydro.send_xv = NULL;
    c->mpi.hydro.send_rho = NULL;
    c->mpi.hydro.send_gradient = NULL;
    c->mpi.grav.send = NULL;
248
    c->mpi.send_ti = NULL;
249
250
251
252
#endif
  }
}

253
254
255
256
/**
 * @brief Free up any allocated cells.
 */
void space_free_cells(struct space *s) {
257
258
259

  ticks tic = getticks();

Matthieu Schaller's avatar
Matthieu Schaller committed
260
261
  threadpool_map(&s->e->threadpool, space_rebuild_recycle_mapper, s->cells_top,
                 s->nr_cells, sizeof(struct cell), 0, s);
262
  s->maxdepth = 0;
263
264
265
266

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

269
/**
270
 * @brief Re-build the top-level cell grid.
271
 *
272
 * @param s The #space.
273
 * @param verbose Print messages to stdout or not.
274
 */
275
void space_regrid(struct space *s, int verbose) {
276

277
  const size_t nr_parts = s->nr_parts;
Loic Hausammann's avatar
Loic Hausammann committed
278
  const size_t nr_sparts = s->nr_sparts;
279
  const ticks tic = getticks();
280
  const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0;
281

282
  /* Run through the cells and get the current h_max. */
283
  // tic = getticks();
284
  float h_max = s->cell_min / kernel_gamma / space_stretch;
285
  if (nr_parts > 0) {
286
287
288
289
290
291

    /* 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]];
292
        if (c->hydro.h_max > h_max) {
293
          h_max = c->hydro.h_max;
294
        }
Loic Hausammann's avatar
Loic Hausammann committed
295
        if (c->stars.h_max > h_max) {
296
          h_max = c->stars.h_max;
Loic Hausammann's avatar
Loic Hausammann committed
297
        }
298
      }
299
300

      /* Can we instead use all the top-level cells? */
301
    } else if (s->cells_top != NULL) {
Tom Theuns's avatar
Tom Theuns committed
302
      for (int k = 0; k < s->nr_cells; k++) {
303
        const struct cell *c = &s->cells_top[k];
304
        if (c->nodeID == engine_rank && c->hydro.h_max > h_max) {
305
          h_max = c->hydro.h_max;
306
        }
Loic Hausammann's avatar
Loic Hausammann committed
307
        if (c->nodeID == engine_rank && c->stars.h_max > h_max) {
308
          h_max = c->stars.h_max;
Loic Hausammann's avatar
Loic Hausammann committed
309
        }
310
      }
311
312

      /* Last option: run through the particles */
313
    } else {
314
      for (size_t k = 0; k < nr_parts; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
315
        if (s->parts[k].h > h_max) h_max = s->parts[k].h;
316
      }
Loic Hausammann's avatar
Loic Hausammann committed
317
318
319
      for (size_t k = 0; k < nr_sparts; k++) {
        if (s->sparts[k].h > h_max) h_max = s->sparts[k].h;
      }
320
321
322
323
324
325
326
327
328
329
    }
  }

/* 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)
330
      error("Failed to aggregate the rebuild flag across nodes.");
331
332
333
    h_max = buff;
  }
#endif
334
  if (verbose) message("h_max is %.3e (cell_min=%.3e).", h_max, s->cell_min);
335
336

  /* Get the new putative cell dimensions. */
337
  const int cdim[3] = {
338
339
340
341
342
343
      (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))};
344
345
346
347
348

  /* 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 "
349
350
351
        "is switched on.\nThis error is often caused by any of the "
        "followings:\n"
        " - too few particles to generate a sensible grid,\n"
352
353
        " - the initial value of 'Scheduler:max_top_level_cells' is too "
        "small,\n"
354
        " - the (minimal) time-step is too large leading to particles with "
355
        "predicted smoothing lengths too large for the box size,\n"
356
        " - particles with velocities so large that they move by more than two "
357
        "box sizes per time-step.\n");
358

359
360
361
/* 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. */
362
#ifdef WITH_MPI
Matthieu Schaller's avatar
Matthieu Schaller committed
363
  double oldwidth[3];
364
365
366
367
368
369
370
371
  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];
372
373
374
    oldwidth[0] = s->width[0];
    oldwidth[1] = s->width[1];
    oldwidth[2] = s->width[2];
375
376
377
378
379
380
381
382
383

    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);
384
          oldnodeIDs[cid] = s->cells_top[cid].nodeID;
385
386
387
388
389
        }
      }
    }
  }

Peter W. Draper's avatar
Peter W. Draper committed
390
  /* Are we about to allocate new top level cells without a regrid?
Peter W. Draper's avatar
Peter W. Draper committed
391
   * Can happen when restarting the application. */
392
  const int no_regrid = (s->cells_top == NULL && oldnodeIDs == NULL);
393
394
395
396
#endif

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

400
401
/* Be verbose about this. */
#ifdef SWIFT_DEBUG_CHECKS
402
    message("(re)griding space cdim=(%d %d %d)", cdim[0], cdim[1], cdim[2]);
403
404
405
    fflush(stdout);
#endif

406
    /* Free the old cells, if they were allocated. */
407
    if (s->cells_top != NULL) {
408
      space_free_cells(s);
409
      free(s->local_cells_with_tasks_top);
410
      free(s->local_cells_top);
411
      free(s->cells_with_particles_top);
412
      free(s->local_cells_with_particles_top);
413
      free(s->cells_top);
414
      free(s->multipoles_top);
415
416
    }

417
418
419
420
    /* 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);

421
    /* Set the new cell dimensions only if smaller. */
422
    for (int k = 0; k < 3; k++) {
423
      s->cdim[k] = cdim[k];
424
425
      s->width[k] = s->dim[k] / cdim[k];
      s->iwidth[k] = 1.0 / s->width[k];
426
    }
427
    const float dmin = min3(s->width[0], s->width[1], s->width[2]);
428
429
430

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

436
437
    /* Allocate the multipoles for the top-level cells. */
    if (s->gravity) {
438
      if (posix_memalign((void **)&s->multipoles_top, multipole_align,
439
                         s->nr_cells * sizeof(struct gravity_tensors)) != 0)
440
        error("Failed to allocate top-level multipoles.");
441
      bzero(s->multipoles_top, s->nr_cells * sizeof(struct gravity_tensors));
442
443
    }

444
    /* Allocate the indices of local cells */
445
    if (posix_memalign((void **)&s->local_cells_top, SWIFT_STRUCT_ALIGNMENT,
446
447
448
449
                       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));

450
    /* Allocate the indices of local cells with tasks */
451
452
    if (posix_memalign((void **)&s->local_cells_with_tasks_top,
                       SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0)
453
454
455
      error("Failed to allocate indices of local top-level cells with tasks.");
    bzero(s->local_cells_with_tasks_top, s->nr_cells * sizeof(int));

456
    /* Allocate the indices of cells with particles */
457
    if (posix_memalign((void **)&s->cells_with_particles_top,
458
                       SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0)
459
460
      error("Failed to allocate indices of top-level cells with particles.");
    bzero(s->cells_with_particles_top, s->nr_cells * sizeof(int));
461

462
463
464
465
466
467
468
469
    /* 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));

470
    /* Set the cells' locks */
471
    for (int k = 0; k < s->nr_cells; k++) {
472
      if (lock_init(&s->cells_top[k].hydro.lock) != 0)
473
        error("Failed to init spinlock for hydro.");
474
      if (lock_init(&s->cells_top[k].grav.plock) != 0)
475
        error("Failed to init spinlock for gravity.");
476
      if (lock_init(&s->cells_top[k].grav.mlock) != 0)
477
        error("Failed to init spinlock for multipoles.");
478
      if (lock_init(&s->cells_top[k].stars.lock) != 0)
479
480
        error("Failed to init spinlock for stars.");
    }
481
482

    /* Set the cell location and sizes. */
483
484
485
    for (int i = 0; i < cdim[0]; i++)
      for (int j = 0; j < cdim[1]; j++)
        for (int k = 0; k < cdim[2]; k++) {
486
487
          const size_t cid = cell_getid(cdim, i, j, k);
          struct cell *restrict c = &s->cells_top[cid];
488
489
490
491
492
493
          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];
494
495
          c->dmin = dmin;
          c->depth = 0;
496
          c->hydro.count = 0;
497
498
          c->grav.count = 0;
          c->stars.count = 0;
499
          c->super = c;
500
501
          c->hydro.super = c;
          c->grav.super = c;
502
503
          c->hydro.ti_old_part = ti_current;
          c->grav.ti_old_part = ti_current;
504
          c->grav.ti_old_multipole = ti_current;
Pedro Gonnet's avatar
Pedro Gonnet committed
505
#ifdef WITH_MPI
506
          c->mpi.tag = -1;
Pedro Gonnet's avatar
Pedro Gonnet committed
507
#endif  // WITH_MPI
508
          if (s->gravity) c->grav.multipole = &s->multipoles_top[cid];
509
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
510
511
          c->cellID = -last_cell_id;
          last_cell_id++;
512
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
513
        }
514
515

    /* Be verbose about the change. */
516
517
518
    if (verbose)
      message("set cell dimensions to [ %i %i %i ].", cdim[0], cdim[1],
              cdim[2]);
519

520
#ifdef WITH_MPI
521
522
523
524
525
    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)
526
527
528
        message(
            "basic cell dimensions have increased - recalculating the "
            "global partition.");
529

Matthieu Schaller's avatar
Matthieu Schaller committed
530
      if (!partition_space_to_space(oldwidth, oldcdim, oldnodeIDs, s)) {
531
532
533
534

        /* Failed, try another technique that requires no settings. */
        message("Failed to get a new partition, trying less optimal method");
        struct partition initial_partition;
535
#if defined(HAVE_PARMETIS) || defined(HAVE_METIS)
536
537
538
539
540
541
542
543
544
545
546
547
548
        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);
549

550
551
      /* Finished with these. */
      free(oldnodeIDs);
Peter W. Draper's avatar
Peter W. Draper committed
552
553

    } else if (no_regrid && s->e != NULL) {
Peter W. Draper's avatar
Peter W. Draper committed
554
555
556
557
558
559
560
561
562
563
564
565
566
      /* 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);
567
    }
Pedro Gonnet's avatar
Pedro Gonnet committed
568
#endif /* WITH_MPI */
569
570
571
572

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

573
  }      /* re-build upper-level cells? */
574
  else { /* Otherwise, just clean up the cells. */
575
576

    /* Free the old cells, if they were allocated. */
577
    space_free_cells(s);
578
  }
579
580
581
582

  if (verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
583
}
584
585
586
587
588

/**
 * @brief Re-build the cells as well as the tasks.
 *
 * @param s The #space in which to update the cells.
589
 * @param repartitioned Did we just repartition?
590
 * @param verbose Print messages to stdout or not
591
 */
592
void space_rebuild(struct space *s, int repartitioned, int verbose) {
593

Matthieu Schaller's avatar
Matthieu Schaller committed
594
  const ticks tic = getticks();
595

596
597
/* Be verbose about this. */
#ifdef SWIFT_DEBUG_CHECKS
Pedro Gonnet's avatar
Pedro Gonnet committed
598
  if (s->e->nodeID == 0 || verbose) message("(re)building space");
599
600
  fflush(stdout);
#endif
601
602

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

Pedro Gonnet's avatar
Pedro Gonnet committed
605
606
  size_t nr_parts = s->nr_parts;
  size_t nr_gparts = s->nr_gparts;
607
  size_t nr_sparts = s->nr_sparts;
608
609
610
  int count_inhibited_parts = 0;
  int count_inhibited_gparts = 0;
  int count_inhibited_sparts = 0;
611
  struct cell *restrict cells_top = s->cells_top;
612
  const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0;
613

Pedro Gonnet's avatar
Pedro Gonnet committed
614
615
616
  /* Run through the particles and get their cell index. Allocates
     an index that is larger than the number of particles to avoid
     re-allocating after shuffling. */
617
  const size_t ind_size = s->size_parts + 100;
618
619
620
621
622
623
  int *ind = (int *)malloc(sizeof(int) * ind_size);
  if (ind == NULL) error("Failed to allocate temporary particle indices.");
  int *cell_part_counts = (int *)calloc(sizeof(int), s->nr_cells);
  if (cell_part_counts == NULL)
    error("Failed to allocate cell part count buffer.");
  if (s->size_parts > 0)
624
625
    space_parts_get_cell_index(s, ind, cell_part_counts, &count_inhibited_parts,
                               verbose);
626

627
  /* Run through the gravity particles and get their cell index. */
628
  const size_t gind_size = s->size_gparts + 100;
629
630
631
632
633
  int *gind = (int *)malloc(sizeof(int) * gind_size);
  if (gind == NULL) error("Failed to allocate temporary g-particle indices.");
  int *cell_gpart_counts = (int *)calloc(sizeof(int), s->nr_cells);
  if (cell_gpart_counts == NULL)
    error("Failed to allocate cell gpart count buffer.");
Pedro Gonnet's avatar
Pedro Gonnet committed
634
  if (s->size_gparts > 0)
635
636
    space_gparts_get_cell_index(s, gind, cell_gpart_counts,
                                &count_inhibited_gparts, verbose);
637

638
639
  /* Run through the star particles and get their cell index. */
  const size_t sind_size = s->size_sparts + 100;
640
641
642
643
644
  int *sind = (int *)malloc(sizeof(int) * sind_size);
  if (sind == NULL) error("Failed to allocate temporary s-particle indices.");
  int *cell_spart_counts = (int *)calloc(sizeof(int), s->nr_cells);
  if (cell_spart_counts == NULL)
    error("Failed to allocate cell gpart count buffer.");
645
  if (s->size_sparts > 0)
646
647
    space_sparts_get_cell_index(s, sind, cell_spart_counts,
                                &count_inhibited_sparts, verbose);
648

649
650
651
652
653
654
655
656
#ifdef SWIFT_DEBUG_CHECKS
  if (repartitioned && count_inhibited_parts)
    error("We just repartitioned but still found inhibited parts.");
  if (repartitioned && count_inhibited_sparts)
    error("We just repartitioned but still found inhibited sparts.");
  if (repartitioned && count_inhibited_gparts)
    error("We just repartitioned but still found inhibited gparts.");
#endif
657

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

660
  /* Move non-local parts and inhibited parts to the end of the list. */
661
  if (!repartitioned && (s->e->nr_nodes > 1 || count_inhibited_parts > 0)) {
662
    for (size_t k = 0; k < nr_parts; /* void */) {
663
664

      /* Inhibited particle or foreign particle */
665
      if (ind[k] == -1 || cells_top[ind[k]].nodeID != local_nodeID) {
666
667

        /* One fewer particle */
668
        nr_parts -= 1;
669

670
671
        /* Swap the particle */
        memswap(&s->parts[k], &s->parts[nr_parts], sizeof(struct part));
672

673
674
675
676
677
678
679
        /* Swap the link with the gpart */
        if (s->parts[k].gpart != NULL) {
          s->parts[k].gpart->id_or_neg_offset = -k;
        }
        if (s->parts[nr_parts].gpart != NULL) {
          s->parts[nr_parts].gpart->id_or_neg_offset = -nr_parts;
        }
680

681
682
683
684
        /* Swap the xpart */
        memswap(&s->xparts[k], &s->xparts[nr_parts], sizeof(struct xpart));
        /* Swap the index */
        memswap(&ind[k], &ind[nr_parts], sizeof(int));
685

686
      } else {
687
688
        /* Increment when not exchanging otherwise we need to retest "k".*/
        k++;
689
      }
690
691
    }
  }
692

693
#ifdef SWIFT_DEBUG_CHECKS
Peter W. Draper's avatar
Peter W. Draper committed
694
  /* Check that all parts are in the correct places. */
695
  int check_count_inhibited_part = 0;
696
  for (size_t k = 0; k < nr_parts; k++) {
697
    if (ind[k] == -1 || cells_top[ind[k]].nodeID != local_nodeID) {
698
699
700
701
      error("Failed to move all non-local parts to send list");
    }
  }
  for (size_t k = nr_parts; k < s->nr_parts; k++) {
702
    if (ind[k] != -1 && cells_top[ind[k]].nodeID == local_nodeID) {
703
      error("Failed to remove local parts from send list");
704
    }
705
    if (ind[k] == -1) ++check_count_inhibited_part;
706
  }
707
708
  if (check_count_inhibited_part != count_inhibited_parts)
    error("Counts of inhibited particles do not match!");
709
#endif /* SWIFT_DEBUG_CHECKS */
710

711
  /* Move non-local sparts and inhibited sparts to the end of the list. */
712
  if (!repartitioned && (s->e->nr_nodes > 1 || count_inhibited_sparts > 0)) {
713
    for (size_t k = 0; k < nr_sparts; /* void */) {
714
715

      /* Inhibited particle or foreign particle */
716
      if (sind[k] == -1 || cells_top[sind[k]].nodeID != local_nodeID) {
717
718

        /* One fewer particle */
719
        nr_sparts -= 1;
720

721
722
        /* Swap the particle */
        memswap(&s->sparts[k], &s->sparts[nr_sparts], sizeof(struct spart));
723

724
725
726
727
728
729
730
        /* Swap the link with the gpart */
        if (s->sparts[k].gpart != NULL) {
          s->sparts[k].gpart->id_or_neg_offset = -k;
        }
        if (s->sparts[nr_sparts].gpart != NULL) {
          s->sparts[nr_sparts].gpart->id_or_neg_offset = -nr_sparts;
        }
731

732
733
        /* Swap the index */
        memswap(&sind[k], &sind[nr_sparts], sizeof(int));
734

735
      } else {
736
737
        /* Increment when not exchanging otherwise we need to retest "k".*/
        k++;
738
739
740
741
742
      }
    }
  }

#ifdef SWIFT_DEBUG_CHECKS
743
744
  /* Check that all sparts are in the correct place. */
  int check_count_inhibited_spart = 0;
745
  for (size_t k = 0; k < nr_sparts; k++) {
746
    if (sind[k] == -1 || cells_top[sind[k]].nodeID != local_nodeID) {
747
748
749
750
      error("Failed to move all non-local sparts to send list");
    }
  }
  for (size_t k = nr_sparts; k < s->nr_sparts; k++) {
751
    if (sind[k] != -1 && cells_top[sind[k]].nodeID == local_nodeID) {
752
753
      error("Failed to remove local sparts from send list");
    }
754
    if (sind[k] == -1) ++check_count_inhibited_spart;
755
  }
756
757
  if (check_count_inhibited_spart != count_inhibited_sparts)
    error("Counts of inhibited s-particles do not match!");
758
#endif /* SWIFT_DEBUG_CHECKS */
759

760
  /* Move non-local gparts and inhibited parts to the end of the list. */
761
  if (!repartitioned && (s->e->nr_nodes > 1 || count_inhibited_gparts > 0)) {
762
    for (size_t k = 0; k < nr_gparts; /* void */) {
763
764

      /* Inhibited particle or foreign particle */
765
      if (gind[k] == -1 || cells_top[gind[k]].nodeID != local_nodeID) {
766
767

        /* One fewer particle */
768
        nr_gparts -= 1;
769

770
771
        /* Swap the particle */
        memswap(&s->gparts[k], &s->gparts[nr_gparts], sizeof(struct gpart));
772

773
774
775
776
777
778
779
780
781
782
783
784
785
        /* Swap the link with part/spart */
        if (s->gparts[k].type == swift_type_gas) {
          s->parts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
        } else if (s->gparts[k].type == swift_type_stars) {
          s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
        }
        if (s->gparts[nr_gparts].type == swift_type_gas) {
          s->parts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
              &s->gparts[nr_gparts];
        } else if (s->gparts[nr_gparts].type == swift_type_stars) {
          s->sparts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
              &s->gparts[nr_gparts];
        }
786

787
788
        /* Swap the index */
        memswap(&gind[k], &gind[nr_gparts], sizeof(int));
789
      } else {
790
791
        /* Increment when not exchanging otherwise we need to retest "k".*/
        k++;
792
      }
793
794
    }
  }
795

796
#ifdef SWIFT_DEBUG_CHECKS
797
798
  /* Check that all gparts are in the correct place. */
  int check_count_inhibited_gpart = 0;
799
  for (size_t k = 0; k < nr_gparts; k++) {
800
    if (gind[k] == -1 || cells_top[gind[k]].nodeID != local_nodeID) {
801
802
803
804
      error("Failed to move all non-local gparts to send list");
    }
  }
  for (size_t k = nr_gparts; k < s->nr_gparts; k++) {
805
    if (gind[k] != -1 && cells_top[gind[k]].nodeID == local_nodeID) {
806
      error("Failed to remove local gparts from send list");
807
    }
808
    if (gind[k] == -1) ++check_count_inhibited_gpart;
809
  }
810
811
  if (check_count_inhibited_gpart != count_inhibited_gparts)
    error("Counts of inhibited g-particles do not match!");
812
813
814
#endif /* SWIFT_DEBUG_CHECKS */

#ifdef WITH_MPI
815

816
  /* Exchange the strays, note that this potentially re-allocates
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
     the parts arrays. This can be skipped if we just repartitioned aspace
     there should be no strays */
  if (!repartitioned) {

    size_t nr_parts_exchanged = s->nr_parts - nr_parts;
    size_t nr_gparts_exchanged = s->nr_gparts - nr_gparts;
    size_t nr_sparts_exchanged = s->nr_sparts - nr_sparts;
    engine_exchange_strays(s->e, nr_parts, &ind[nr_parts], &nr_parts_exchanged,
                           nr_gparts, &gind[nr_gparts], &nr_gparts_exchanged,
                           nr_sparts, &sind[nr_sparts], &nr_sparts_exchanged);

    /* Set the new particle counts. */
    s->nr_parts = nr_parts + nr_parts_exchanged;
    s->nr_gparts = nr_gparts + nr_gparts_exchanged;
    s->nr_sparts = nr_sparts + nr_sparts_exchanged;
  } else {
#ifdef SWIFT_DEBUG_CHECKS
    if (s->nr_parts != nr_parts)
      error("Number of parts changing after repartition");
    if (s->nr_sparts != nr_sparts)
      error("Number of sparts changing after repartition");
    if (s->nr_gparts != nr_gparts)
      error("Number of gparts changing after repartition");
#endif
  }
842

843
844
845
846
847
848
849
850
  /* Clear non-local cell counts. */
  for (int k = 0; k < s->nr_cells; k++) {
    if (s->cells_top[k].nodeID != local_nodeID) {
      cell_part_counts[k] = 0;
      cell_spart_counts[k] = 0;
      cell_gpart_counts[k] = 0;
    }
  }
851

852
  /* Re-allocate the index array for the parts if needed.. */
853
  if (s->nr_parts + 1 > ind_size) {
854
    int *ind_new;
855
    if ((ind_new = (int *)malloc(sizeof(int) * (s->nr_parts + 1))) == NULL)
856
      error("Failed to allocate temporary particle indices.");
857
    memcpy(ind_new, ind, sizeof(int) * nr_parts);
858
859
    free(ind);
    ind = ind_new;
860
861
  }

862
  /* Re-allocate the index array for the sparts if needed.. */
863
  if (s->nr_sparts + 1 > sind_size) {
864
865
    int *sind_new;
    if ((sind_new = (int *)malloc(sizeof(int) * (s->nr_sparts + 1))) == NULL)
866
      error("Failed to allocate temporary s-particle indices.");
867
868
869
870
871
    memcpy(sind_new, sind, sizeof(int) * nr_sparts);
    free(sind);
    sind = sind_new;
  }

872
873
874
  const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
  const double ih[3] = {s->iwidth[0],