space.c 125 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
      error("Failed to allocate indices of local top-level cells with tasks.");
454
455
    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->split = 0;
497
          c->hydro.count = 0;
498
499
          c->grav.count = 0;
          c->stars.count = 0;
500
          c->super = c;
501
502
          c->hydro.super = c;
          c->grav.super = c;
503
504
          c->hydro.ti_old_part = ti_current;
          c->grav.ti_old_part = ti_current;
505
          c->grav.ti_old_multipole = ti_current;
Pedro Gonnet's avatar
Pedro Gonnet committed
506
#ifdef WITH_MPI
507
          c->mpi.tag = -1;
Pedro Gonnet's avatar
Pedro Gonnet committed
508
#endif  // WITH_MPI
509
          if (s->gravity) c->grav.multipole = &s->multipoles_top[cid];
510
#ifdef SWIFT_DEBUG_CHECKS
Matthieu Schaller's avatar
Matthieu Schaller committed
511
512
          c->cellID = -last_cell_id;
          last_cell_id++;
513
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
514
        }
515
516

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Pedro Gonnet's avatar
Pedro Gonnet committed
615
616
617
  /* 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. */
618
  const size_t ind_size = s->size_parts + 100;
619
620
621
622
623
624
  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)
625
626
    space_parts_get_cell_index(s, ind, cell_part_counts, &count_inhibited_parts,
                               verbose);
627

628
  /* Run through the gravity particles and get their cell index. */
629
  const size_t gind_size = s->size_gparts + 100;
630
631
632
633
634
  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
635
  if (s->size_gparts > 0)
636
637
    space_gparts_get_cell_index(s, gind, cell_gpart_counts,
                                &count_inhibited_gparts, verbose);
638

639
640
  /* Run through the star particles and get their cell index. */
  const size_t sind_size = s->size_sparts + 100;
641
642
643
644
645
  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.");
646
  if (s->size_sparts > 0)
647
648
    space_sparts_get_cell_index(s, sind, cell_spart_counts,
                                &count_inhibited_sparts, verbose);
649

650
651
652
653
654
655
656
657
#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
658

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

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

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

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

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

674
675
676
677
678
679
680
        /* 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;
        }
681

682
683
684
685
        /* 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));
686

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

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

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

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

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

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

725
726
727
728
729
730
731
        /* 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;
        }
732

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

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

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

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

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

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

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

774
775
776
777
778
779
780
781
782
783
784
785
786
        /* 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];
        }
787

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

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

#ifdef WITH_MPI
816

817
  /* Exchange the strays, note that this potentially re-allocates
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
     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
  }
843

844
845
846
847
848
849
850
851
  /* 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;
    }
  }
852

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

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