space.c 121 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
    if (s->local_cells_top != NULL) {
      for (int k = 0; k < s->nr_local_cells; ++k) {
        const struct cell *c = &s->cells_top[s->local_cells_top[k]];
289
290
        if (c->hydro.h_max > h_max) {
          h_max = s->cells_top[k].hydro.h_max;
291
        }
Loic Hausammann's avatar
Loic Hausammann committed
292
293
294
        if (c->stars.h_max > h_max) {
          h_max = s->cells_top[k].stars.h_max;
        }
295
296
      }
    } else if (s->cells_top != NULL) {
Tom Theuns's avatar
Tom Theuns committed
297
      for (int k = 0; k < s->nr_cells; k++) {
298
        const struct cell *c = &s->cells_top[k];
299
300
        if (c->nodeID == engine_rank && c->hydro.h_max > h_max) {
          h_max = s->cells_top[k].hydro.h_max;
301
        }
Loic Hausammann's avatar
Loic Hausammann committed
302
303
304
        if (c->nodeID == engine_rank && c->stars.h_max > h_max) {
          h_max = s->cells_top[k].stars.h_max;
        }
305
306
      }
    } else {
307
      for (size_t k = 0; k < nr_parts; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
308
        if (s->parts[k].h > h_max) h_max = s->parts[k].h;
309
      }
Loic Hausammann's avatar
Loic Hausammann committed
310
311
312
      for (size_t k = 0; k < nr_sparts; k++) {
        if (s->sparts[k].h > h_max) h_max = s->sparts[k].h;
      }
313
314
315
316
317
318
319
320
321
322
    }
  }

/* 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)
323
      error("Failed to aggregate the rebuild flag across nodes.");
324
325
326
    h_max = buff;
  }
#endif
327
  if (verbose) message("h_max is %.3e (cell_min=%.3e).", h_max, s->cell_min);
328
329

  /* Get the new putative cell dimensions. */
330
  const int cdim[3] = {
331
332
333
334
335
336
      (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))};
337
338
339
340
341

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

352
353
354
/* 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. */
355
#ifdef WITH_MPI
Matthieu Schaller's avatar
Matthieu Schaller committed
356
  double oldwidth[3];
357
358
359
360
361
362
363
364
  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];
365
366
367
    oldwidth[0] = s->width[0];
    oldwidth[1] = s->width[1];
    oldwidth[2] = s->width[2];
368
369
370
371
372
373
374
375
376

    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);
377
          oldnodeIDs[cid] = s->cells_top[cid].nodeID;
378
379
380
381
382
        }
      }
    }
  }

Peter W. Draper's avatar
Peter W. Draper committed
383
  /* Are we about to allocate new top level cells without a regrid?
Peter W. Draper's avatar
Peter W. Draper committed
384
   * Can happen when restarting the application. */
385
  const int no_regrid = (s->cells_top == NULL && oldnodeIDs == NULL);
386
387
388
389
#endif

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

393
394
/* Be verbose about this. */
#ifdef SWIFT_DEBUG_CHECKS
395
    message("(re)griding space cdim=(%d %d %d)", cdim[0], cdim[1], cdim[2]);
396
397
398
    fflush(stdout);
#endif

399
    /* Free the old cells, if they were allocated. */
400
    if (s->cells_top != NULL) {
401
      space_free_cells(s);
402
      free(s->local_cells_with_tasks_top);
403
      free(s->local_cells_top);
404
      free(s->cells_top);
405
      free(s->multipoles_top);
406
407
    }

408
409
410
411
    /* 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);

412
    /* Set the new cell dimensions only if smaller. */
413
    for (int k = 0; k < 3; k++) {
414
      s->cdim[k] = cdim[k];
415
416
      s->width[k] = s->dim[k] / cdim[k];
      s->iwidth[k] = 1.0 / s->width[k];
417
    }
418
    const float dmin = min3(s->width[0], s->width[1], s->width[2]);
419
420
421

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

427
428
    /* Allocate the multipoles for the top-level cells. */
    if (s->gravity) {
429
      if (posix_memalign((void **)&s->multipoles_top, multipole_align,
430
                         s->nr_cells * sizeof(struct gravity_tensors)) != 0)
431
        error("Failed to allocate top-level multipoles.");
432
      bzero(s->multipoles_top, s->nr_cells * sizeof(struct gravity_tensors));
433
434
    }

435
    /* Allocate the indices of local cells */
436
    if (posix_memalign((void **)&s->local_cells_top, SWIFT_STRUCT_ALIGNMENT,
437
438
439
440
                       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));

441
    /* Allocate the indices of local cells with tasks */
442
443
    if (posix_memalign((void **)&s->local_cells_with_tasks_top,
                       SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0)
444
445
446
      error("Failed to allocate indices of local top-level cells.");
    bzero(s->local_cells_with_tasks_top, s->nr_cells * sizeof(int));

447
    /* Set the cells' locks */
448
    for (int k = 0; k < s->nr_cells; k++) {
449
      if (lock_init(&s->cells_top[k].hydro.lock) != 0)
450
        error("Failed to init spinlock for hydro.");
451
      if (lock_init(&s->cells_top[k].grav.plock) != 0)
452
        error("Failed to init spinlock for gravity.");
453
      if (lock_init(&s->cells_top[k].grav.mlock) != 0)
454
        error("Failed to init spinlock for multipoles.");
455
      if (lock_init(&s->cells_top[k].stars.lock) != 0)
456
457
        error("Failed to init spinlock for stars.");
    }
458
459

    /* Set the cell location and sizes. */
460
461
462
    for (int i = 0; i < cdim[0]; i++)
      for (int j = 0; j < cdim[1]; j++)
        for (int k = 0; k < cdim[2]; k++) {
463
464
          const size_t cid = cell_getid(cdim, i, j, k);
          struct cell *restrict c = &s->cells_top[cid];
465
466
467
468
469
470
          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];
471
472
          c->dmin = dmin;
          c->depth = 0;
473
          c->hydro.count = 0;
474
475
          c->grav.count = 0;
          c->stars.count = 0;
476
          c->super = c;
477
478
          c->hydro.super = c;
          c->grav.super = c;
479
480
          c->hydro.ti_old_part = ti_current;
          c->grav.ti_old_part = ti_current;
481
          c->grav.ti_old_multipole = ti_current;
Pedro Gonnet's avatar
Pedro Gonnet committed
482
#ifdef WITH_MPI
483
          c->mpi.tag = -1;
Pedro Gonnet's avatar
Pedro Gonnet committed
484
#endif  // WITH_MPI
485
          if (s->gravity) c->grav.multipole = &s->multipoles_top[cid];
486
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
487
488
          c->cellID = -last_cell_id;
          last_cell_id++;
489
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
490
        }
491
492

    /* Be verbose about the change. */
493
494
495
    if (verbose)
      message("set cell dimensions to [ %i %i %i ].", cdim[0], cdim[1],
              cdim[2]);
496

497
#ifdef WITH_MPI
498
499
500
501
502
    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)
503
504
505
        message(
            "basic cell dimensions have increased - recalculating the "
            "global partition.");
506

Matthieu Schaller's avatar
Matthieu Schaller committed
507
      if (!partition_space_to_space(oldwidth, oldcdim, oldnodeIDs, s)) {
508
509
510
511

        /* Failed, try another technique that requires no settings. */
        message("Failed to get a new partition, trying less optimal method");
        struct partition initial_partition;
512
#if defined(HAVE_PARMETIS) || defined(HAVE_METIS)
513
514
515
516
517
518
519
520
521
522
523
524
525
        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);
526

527
528
      /* Finished with these. */
      free(oldnodeIDs);
Peter W. Draper's avatar
Peter W. Draper committed
529
530

    } else if (no_regrid && s->e != NULL) {
Peter W. Draper's avatar
Peter W. Draper committed
531
532
533
534
535
536
537
538
539
540
541
542
543
      /* 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);
544
    }
Pedro Gonnet's avatar
Pedro Gonnet committed
545
#endif /* WITH_MPI */
546
547
548
549

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

550
  }      /* re-build upper-level cells? */
551
  else { /* Otherwise, just clean up the cells. */
552
553

    /* Free the old cells, if they were allocated. */
554
    space_free_cells(s);
555
  }
556
557
558
559

  if (verbose)
    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
            clocks_getunit());
560
}
561
562
563
564
565

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

Matthieu Schaller's avatar
Matthieu Schaller committed
571
  const ticks tic = getticks();
572

573
574
/* Be verbose about this. */
#ifdef SWIFT_DEBUG_CHECKS
Pedro Gonnet's avatar
Pedro Gonnet committed
575
  if (s->e->nodeID == 0 || verbose) message("(re)building space");
576
577
  fflush(stdout);
#endif
578
579

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

Pedro Gonnet's avatar
Pedro Gonnet committed
582
583
  size_t nr_parts = s->nr_parts;
  size_t nr_gparts = s->nr_gparts;
584
  size_t nr_sparts = s->nr_sparts;
585
586
587
  int count_inhibited_parts = 0;
  int count_inhibited_gparts = 0;
  int count_inhibited_sparts = 0;
588
  struct cell *restrict cells_top = s->cells_top;
589
  const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0;
590

Pedro Gonnet's avatar
Pedro Gonnet committed
591
592
593
  /* 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. */
594
  const size_t ind_size = s->size_parts + 100;
595
596
597
598
599
600
  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)
601
602
    space_parts_get_cell_index(s, ind, cell_part_counts, &count_inhibited_parts,
                               verbose);
603

604
  /* Run through the gravity particles and get their cell index. */
605
  const size_t gind_size = s->size_gparts + 100;
606
607
608
609
610
  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
611
  if (s->size_gparts > 0)
612
613
    space_gparts_get_cell_index(s, gind, cell_gpart_counts,
                                &count_inhibited_gparts, verbose);
614

615
616
  /* Run through the star particles and get their cell index. */
  const size_t sind_size = s->size_sparts + 100;
617
618
619
620
621
  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.");
622
  if (s->size_sparts > 0)
623
624
    space_sparts_get_cell_index(s, sind, cell_spart_counts,
                                &count_inhibited_sparts, verbose);
625

626
627
628
629
630
631
632
633
634
#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

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

637
  /* Move non-local parts and inhibited parts to the end of the list. */
638
  if (!repartitioned && (s->e->nr_nodes > 1 || count_inhibited_parts > 0)) {
639
    for (size_t k = 0; k < nr_parts; /* void */) {
640
641

      /* Inhibited particle or foreign particle */
642
      if (ind[k] == -1 || cells_top[ind[k]].nodeID != local_nodeID) {
643
644

        /* One fewer particle */
645
        nr_parts -= 1;
646

647
648
        /* Swap the particle */
        memswap(&s->parts[k], &s->parts[nr_parts], sizeof(struct part));
649

650
651
652
653
654
655
656
        /* 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;
        }
657

658
659
660
661
        /* 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));
662

663
      } else {
664
665
        /* Increment when not exchanging otherwise we need to retest "k".*/
        k++;
666
      }
667
668
    }
  }
669

670
#ifdef SWIFT_DEBUG_CHECKS
Peter W. Draper's avatar
Peter W. Draper committed
671
  /* Check that all parts are in the correct places. */
672
  int check_count_inhibited_part = 0;
673
  for (size_t k = 0; k < nr_parts; k++) {
674
    if (ind[k] == -1 || cells_top[ind[k]].nodeID != local_nodeID) {
675
676
677
678
      error("Failed to move all non-local parts to send list");
    }
  }
  for (size_t k = nr_parts; k < s->nr_parts; k++) {
679
    if (ind[k] != -1 && cells_top[ind[k]].nodeID == local_nodeID) {
680
      error("Failed to remove local parts from send list");
681
    }
682
    if (ind[k] == -1) ++check_count_inhibited_part;
683
  }
684
685
  if (check_count_inhibited_part != count_inhibited_parts)
    error("Counts of inhibited particles do not match!");
686
#endif /* SWIFT_DEBUG_CHECKS */
687

688
  /* Move non-local sparts and inhibited sparts to the end of the list. */
689
  if (!repartitioned && (s->e->nr_nodes > 1 || count_inhibited_sparts > 0)) {
690
    for (size_t k = 0; k < nr_sparts; /* void */) {
691
692

      /* Inhibited particle or foreign particle */
693
      if (sind[k] == -1 || cells_top[sind[k]].nodeID != local_nodeID) {
694
695

        /* One fewer particle */
696
        nr_sparts -= 1;
697

698
699
        /* Swap the particle */
        memswap(&s->sparts[k], &s->sparts[nr_sparts], sizeof(struct spart));
700

701
702
703
704
705
706
707
        /* 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;
        }
708

709
710
        /* Swap the index */
        memswap(&sind[k], &sind[nr_sparts], sizeof(int));
711

712
      } else {
713
714
        /* Increment when not exchanging otherwise we need to retest "k".*/
        k++;
715
716
717
718
719
      }
    }
  }

#ifdef SWIFT_DEBUG_CHECKS
720
721
  /* Check that all sparts are in the correct place. */
  int check_count_inhibited_spart = 0;
722
  for (size_t k = 0; k < nr_sparts; k++) {
723
    if (sind[k] == -1 || cells_top[sind[k]].nodeID != local_nodeID) {
724
725
726
727
      error("Failed to move all non-local sparts to send list");
    }
  }
  for (size_t k = nr_sparts; k < s->nr_sparts; k++) {
728
    if (sind[k] != -1 && cells_top[sind[k]].nodeID == local_nodeID) {
729
730
      error("Failed to remove local sparts from send list");
    }
731
    if (sind[k] == -1) ++check_count_inhibited_spart;
732
  }
733
734
  if (check_count_inhibited_spart != count_inhibited_sparts)
    error("Counts of inhibited s-particles do not match!");
735
#endif /* SWIFT_DEBUG_CHECKS */
736

737
  /* Move non-local gparts and inhibited parts to the end of the list. */
738
  if (!repartitioned && (s->e->nr_nodes > 1 || count_inhibited_gparts > 0)) {
739
    for (size_t k = 0; k < nr_gparts; /* void */) {
740
741

      /* Inhibited particle or foreign particle */
742
      if (gind[k] == -1 || cells_top[gind[k]].nodeID != local_nodeID) {
743
744

        /* One fewer particle */
745
        nr_gparts -= 1;
746

747
748
        /* Swap the particle */
        memswap(&s->gparts[k], &s->gparts[nr_gparts], sizeof(struct gpart));
749

750
751
752
753
754
755
756
757
758
759
760
761
762
        /* 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];
        }
763

764
765
        /* Swap the index */
        memswap(&gind[k], &gind[nr_gparts], sizeof(int));
766
      } else {
767
768
        /* Increment when not exchanging otherwise we need to retest "k".*/
        k++;
769
      }
770
771
    }
  }
772

773
#ifdef SWIFT_DEBUG_CHECKS
774
775
  /* Check that all gparts are in the correct place. */
  int check_count_inhibited_gpart = 0;
776
  for (size_t k = 0; k < nr_gparts; k++) {
777
    if (gind[k] == -1 || cells_top[gind[k]].nodeID != local_nodeID) {
778
779
780
781
      error("Failed to move all non-local gparts to send list");
    }
  }
  for (size_t k = nr_gparts; k < s->nr_gparts; k++) {
782
    if (gind[k] != -1 && cells_top[gind[k]].nodeID == local_nodeID) {
783
      error("Failed to remove local gparts from send list");
784
    }
785
    if (gind[k] == -1) ++check_count_inhibited_gpart;
786
  }
787
788
  if (check_count_inhibited_gpart != count_inhibited_gparts)
    error("Counts of inhibited g-particles do not match!");
789
790
791
#endif /* SWIFT_DEBUG_CHECKS */

#ifdef WITH_MPI
792

793
  /* Exchange the strays, note that this potentially re-allocates
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
     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
  }
819

820
821
822
823
824
825
826
827
  /* 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;
    }
  }
828

829
  /* Re-allocate the index array for the parts if needed.. */
830
  if (s->nr_parts + 1 > ind_size) {
831
    int *ind_new;
832
    if ((ind_new = (int *)malloc(sizeof(int) * (s->nr_parts + 1))) == NULL)
833
      error("Failed to allocate temporary particle indices.");
834
    memcpy(ind_new, ind, sizeof(int) * nr_parts);
835
836
    free(ind);
    ind = ind_new;
837
838
  }

839
  /* Re-allocate the index array for the sparts if needed.. */
840
  if (s->nr_sparts + 1 > sind_size) {
841
842
    int *sind_new;
    if ((sind_new = (int *)malloc(sizeof(int) * (s->nr_sparts + 1))) == NULL)
843
      error("Failed to allocate temporary s-particle indices.");
844
845
846
847
848
    memcpy(sind_new, sind, sizeof(int) * nr_sparts);
    free(sind);
    sind = sind_new;
  }

849
850
851
  const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
  const double ih[3] = {s->iwidth[0], s->iwidth[1], s->iwidth[2]};

852
  /* Assign each received part to its cell. */
Pedro Gonnet's avatar
Pedro Gonnet committed
853
  for (size_t k = nr_parts; k < s->nr_parts; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
854
    const struct part *const p = &s->parts[k];
855
    ind[k] =
856
        cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
857
    cell_part_counts[ind[k]]++;
858
#ifdef SWIFT_DEBUG_CHECKS
859
    if (cells_top[ind[k]].nodeID != local_nodeID)
860
      error("Received part that does not belong to me (nodeID=%i).",
861
            cells_top[ind[k]].nodeID);
862
#endif
863
  }
864
  nr_parts = s->nr_parts;
865

866
  /* Assign each received spart to its cell. */
867
868
869
870
  for (size_t k = nr_sparts; k < s->nr_sparts; k++) {
    const struct spart *const sp = &s->sparts[k];
    sind[k] =
        cell_getid(cdim, sp->x[0] * ih[0], sp->x[1] * ih[1], sp->x[2] * ih[2]);
871
    cell_spart_counts[sind[k]]++;