space.c 237 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 "black_holes.h"
45
#include "chemistry.h"
46
#include "const.h"
47
#include "cooling.h"
48
#include "debug.h"
49
#include "engine.h"
50
#include "error.h"
51 52
#include "gravity.h"
#include "hydro.h"
53
#include "kernel_hydro.h"
54
#include "lock.h"
55
#include "memswap.h"
56
#include "memuse.h"
57
#include "minmax.h"
58
#include "multipole.h"
59
#include "pressure_floor.h"
60
#include "proxy.h"
61
#include "restart.h"
Loic Hausammann's avatar
Loic Hausammann committed
62
#include "sink.h"
63
#include "sort_part.h"
Loic Hausammann's avatar
Loic Hausammann committed
64
#include "space_unique_id.h"
65
#include "star_formation.h"
66
#include "star_formation_logger.h"
Loic Hausammann's avatar
Loic Hausammann committed
67
#include "stars.h"
68
#include "threadpool.h"
69
#include "tools.h"
70
#include "tracers.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
71 72 73

/* Split size. */
int space_splitsize = space_splitsize_default;
74 75
int space_subsize_pair_hydro = space_subsize_pair_hydro_default;
int space_subsize_self_hydro = space_subsize_self_hydro_default;
76 77
int space_subsize_pair_stars = space_subsize_pair_stars_default;
int space_subsize_self_stars = space_subsize_self_stars_default;
78
int space_subsize_pair_grav = space_subsize_pair_grav_default;
79
int space_subsize_self_grav = space_subsize_self_grav_default;
80
int space_subdepth_diff_grav = space_subdepth_diff_grav_default;
81
int space_maxsize = space_maxsize_default;
82

83 84 85
/*! Number of extra #part we allocate memory for per top-level cell */
int space_extra_parts = space_extra_parts_default;

86 87 88
/*! Number of extra #spart we allocate memory for per top-level cell */
int space_extra_sparts = space_extra_sparts_default;

89 90 91
/*! Number of extra #bpart we allocate memory for per top-level cell */
int space_extra_bparts = space_extra_bparts_default;

92 93 94
/*! Number of extra #gpart we allocate memory for per top-level cell */
int space_extra_gparts = space_extra_gparts_default;

Loic Hausammann's avatar
Loic Hausammann committed
95 96 97
/*! Number of extra #sink we allocate memory for per top-level cell */
int space_extra_sinks = space_extra_sinks_default;

98 99 100
/*! Maximum number of particles per ghost */
int engine_max_parts_per_ghost = engine_max_parts_per_ghost_default;
int engine_max_sparts_per_ghost = engine_max_sparts_per_ghost_default;
101
int engine_max_parts_per_cooling = engine_max_parts_per_cooling_default;
102

103
/*! Maximal depth at which the stars resort task can be pushed */
104
int engine_star_resort_task_depth = engine_star_resort_task_depth_default;
105

106 107
/*! Expected maximal number of strays received at a rebuild */
int space_expected_max_nr_strays = space_expected_max_nr_strays_default;
Loic Hausammann's avatar
Loic Hausammann committed
108
#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_CELL_GRAPH)
109 110
int last_cell_id;
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
111

112 113 114 115 116 117 118 119 120
/**
 * @brief Interval stack necessary for parallel particle sorting.
 */
struct qstack {
  volatile ptrdiff_t i, j;
  volatile int min, max;
  volatile int ready;
};

121 122 123 124 125 126
/**
 * @brief Information required to compute the particle cell indices.
 */
struct index_data {
  struct space *s;
  int *ind;
127
  int *cell_counts;
128 129 130
  size_t count_inhibited_part;
  size_t count_inhibited_gpart;
  size_t count_inhibited_spart;
131
  size_t count_inhibited_bpart;
Loic Hausammann's avatar
Loic Hausammann committed
132
  size_t count_inhibited_sink;
133 134 135
  size_t count_extra_part;
  size_t count_extra_gpart;
  size_t count_extra_spart;
136
  size_t count_extra_bpart;
Loic Hausammann's avatar
Loic Hausammann committed
137
  size_t count_extra_sink;
138 139
};

140
/**
141
 * @brief Recursively dismantle a cell tree.
142
 *
143 144
 * @param s The #space.
 * @param c The #cell to recycle.
Matthieu Schaller's avatar
Matthieu Schaller committed
145 146 147 148 149 150
 * @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.
151
 */
152
void space_rebuild_recycle_rec(struct space *s, struct cell *c,
153 154
                               struct cell **cell_rec_begin,
                               struct cell **cell_rec_end,
155 156
                               struct gravity_tensors **multipole_rec_begin,
                               struct gravity_tensors **multipole_rec_end) {
157
  if (c->split)
158
    for (int k = 0; k < 8; k++)
159
      if (c->progeny[k] != NULL) {
160 161 162 163 164 165
        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];
166

167
        if (s->with_self_gravity) {
168 169
          c->progeny[k]->grav.multipole->next = *multipole_rec_begin;
          *multipole_rec_begin = c->progeny[k]->grav.multipole;
170
        }
171 172

        if (*cell_rec_end == NULL) *cell_rec_end = *cell_rec_begin;
173
        if (s->with_self_gravity && *multipole_rec_end == NULL)
174 175
          *multipole_rec_end = *multipole_rec_begin;

176
        c->progeny[k]->grav.multipole = NULL;
177 178 179 180
        c->progeny[k] = NULL;
      }
}

181 182 183 184 185 186 187 188
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];
189
    struct cell *cell_rec_begin = NULL, *cell_rec_end = NULL;
190 191
    struct gravity_tensors *multipole_rec_begin = NULL,
                           *multipole_rec_end = NULL;
192 193 194 195 196
    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);
197
    c->hydro.sorts = NULL;
Loic Hausammann's avatar
Loic Hausammann committed
198
    c->stars.sorts = NULL;
199
    c->nr_tasks = 0;
200 201 202 203
    c->grav.nr_mm_tasks = 0;
    c->hydro.density = NULL;
    c->hydro.gradient = NULL;
    c->hydro.force = NULL;
204
    c->hydro.limiter = NULL;
205 206
    c->grav.grav = NULL;
    c->grav.mm = NULL;
207 208
    c->hydro.dx_max_part = 0.0f;
    c->hydro.dx_max_sort = 0.0f;
Loic Hausammann's avatar
Loic Hausammann committed
209
    c->stars.dx_max_part = 0.f;
Loic Hausammann's avatar
Loic Hausammann committed
210
    c->stars.dx_max_sort = 0.f;
211
    c->black_holes.dx_max_part = 0.f;
212
    c->hydro.sorted = 0;
Matthieu Schaller's avatar
Matthieu Schaller committed
213
    c->hydro.sort_allocated = 0;
Loic Hausammann's avatar
Loic Hausammann committed
214
    c->stars.sorted = 0;
215
    c->hydro.count = 0;
216
    c->hydro.count_total = 0;
217
    c->hydro.updated = 0;
218
    c->grav.count = 0;
219
    c->grav.count_total = 0;
220
    c->grav.updated = 0;
Loic Hausammann's avatar
Loic Hausammann committed
221
    c->sinks.count = 0;
222
    c->stars.count = 0;
223
    c->stars.count_total = 0;
224
    c->stars.updated = 0;
225 226 227
    c->black_holes.count = 0;
    c->black_holes.count_total = 0;
    c->black_holes.updated = 0;
228 229 230 231 232 233
    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;
234 235
    c->hydro.star_formation = NULL;
    c->hydro.stars_resort = NULL;
236 237
    c->stars.ghost = NULL;
    c->stars.density = NULL;
Alexei Borissov's avatar
Alexei Borissov committed
238
    c->stars.feedback = NULL;
239 240 241
    c->black_holes.density_ghost = NULL;
    c->black_holes.swallow_ghost[0] = NULL;
    c->black_holes.swallow_ghost[1] = NULL;
242
    c->black_holes.swallow_ghost[2] = NULL;
243
    c->black_holes.density = NULL;
244
    c->black_holes.swallow = NULL;
245 246
    c->black_holes.do_gas_swallow = NULL;
    c->black_holes.do_bh_swallow = NULL;
247
    c->black_holes.feedback = NULL;
248 249
    c->kick1 = NULL;
    c->kick2 = NULL;
250
    c->timestep = NULL;
251
    c->timestep_limiter = NULL;
252
    c->timestep_sync = NULL;
253
    c->hydro.end_force = NULL;
254
    c->hydro.drift = NULL;
255
    c->stars.drift = NULL;
256 257
    c->stars.stars_in = NULL;
    c->stars.stars_out = NULL;
258
    c->black_holes.drift = NULL;
259 260
    c->black_holes.black_holes_in = NULL;
    c->black_holes.black_holes_out = NULL;
261
    c->grav.drift = NULL;
262
    c->grav.drift_out = NULL;
263 264
    c->hydro.cooling_in = NULL;
    c->hydro.cooling_out = NULL;
265 266 267 268 269
    c->hydro.cooling = NULL;
    c->grav.long_range = NULL;
    c->grav.down_in = NULL;
    c->grav.down = NULL;
    c->grav.mesh = NULL;
270
    c->grav.end_force = NULL;
271
    c->top = c;
272
    c->super = c;
273 274 275 276
    c->hydro.super = c;
    c->grav.super = c;
    c->hydro.parts = NULL;
    c->hydro.xparts = NULL;
277
    c->grav.parts = NULL;
278
    c->grav.parts_rebuild = NULL;
Loic Hausammann's avatar
Loic Hausammann committed
279
    c->sinks.parts = NULL;
280
    c->stars.parts = NULL;
281
    c->stars.parts_rebuild = NULL;
282
    c->black_holes.parts = NULL;
283
    c->flags = 0;
284 285 286 287
    c->hydro.ti_end_min = -1;
    c->hydro.ti_end_max = -1;
    c->grav.ti_end_min = -1;
    c->grav.ti_end_max = -1;
Loic Hausammann's avatar
Loic Hausammann committed
288
    c->stars.ti_end_min = -1;
289
    c->stars.ti_end_max = -1;
290 291
    c->black_holes.ti_end_min = -1;
    c->black_holes.ti_end_max = -1;
292
    star_formation_logger_init(&c->stars.sfh);
Loic Hausammann's avatar
Loic Hausammann committed
293
#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_CELL_GRAPH)
294 295
    c->cellID = 0;
#endif
296 297
    if (s->with_self_gravity)
      bzero(c->grav.multipole, sizeof(struct gravity_tensors));
298

299 300
    cell_free_hydro_sorts(c);
    cell_free_stars_sorts(c);
301
#if WITH_MPI
302
    c->mpi.tag = -1;
303
    c->mpi.recv = NULL;
304
    c->mpi.send = NULL;
305 306 307 308
#endif
  }
}

309 310
/**
 * @brief Free up any allocated cells.
311 312
 *
 * @param s The #space.
313 314
 */
void space_free_cells(struct space *s) {
315 316 317

  ticks tic = getticks();

Matthieu Schaller's avatar
Matthieu Schaller committed
318
  threadpool_map(&s->e->threadpool, space_rebuild_recycle_mapper, s->cells_top,
319 320
                 s->nr_cells, sizeof(struct cell), threadpool_auto_chunk_size,
                 s);
321
  s->maxdepth = 0;
322 323 324 325

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

328 329 330 331
/**
 * @brief Free any memory in use for foreign particles.
 *
 * @param s The #space.
332 333
 * @param clear_cell_pointers Are we also setting all the foreign cell particle
 * pointers to NULL?
334
 */
335
void space_free_foreign_parts(struct space *s, const int clear_cell_pointers) {
336 337 338

#ifdef WITH_MPI
  if (s->parts_foreign != NULL) {
Peter W. Draper's avatar
Peter W. Draper committed
339
    swift_free("parts_foreign", s->parts_foreign);
340 341 342 343
    s->size_parts_foreign = 0;
    s->parts_foreign = NULL;
  }
  if (s->gparts_foreign != NULL) {
Peter W. Draper's avatar
Peter W. Draper committed
344
    swift_free("gparts_foreign", s->gparts_foreign);
345 346 347 348
    s->size_gparts_foreign = 0;
    s->gparts_foreign = NULL;
  }
  if (s->sparts_foreign != NULL) {
Peter W. Draper's avatar
Peter W. Draper committed
349
    swift_free("sparts_foreign", s->sparts_foreign);
350 351 352
    s->size_sparts_foreign = 0;
    s->sparts_foreign = NULL;
  }
353 354 355 356 357
  if (s->bparts_foreign != NULL) {
    swift_free("bparts_foreign", s->bparts_foreign);
    s->size_bparts_foreign = 0;
    s->bparts_foreign = NULL;
  }
358 359 360 361 362 363 364
  if (clear_cell_pointers) {
    for (int k = 0; k < s->e->nr_proxies; k++) {
      for (int j = 0; j < s->e->proxies[k].nr_cells_in; j++) {
        cell_unlink_foreign_particles(s->e->proxies[k].cells_in[j]);
      }
    }
  }
365 366 367
#endif
}

368
/**
369
 * @brief Re-build the top-level cell grid.
370
 *
371
 * @param s The #space.
372
 * @param verbose Print messages to stdout or not.
373
 */
374
void space_regrid(struct space *s, int verbose) {
375

376
  const size_t nr_parts = s->nr_parts;
Loic Hausammann's avatar
Loic Hausammann committed
377
  const size_t nr_sparts = s->nr_sparts;
378
  const size_t nr_bparts = s->nr_bparts;
379
  const ticks tic = getticks();
380
  const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0;
381

382
  /* Run through the cells and get the current h_max. */
383
  // tic = getticks();
384
  float h_max = s->cell_min / kernel_gamma / space_stretch;
385
  if (nr_parts > 0) {
386 387 388 389 390 391

    /* 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]];
392
        if (c->hydro.h_max > h_max) {
393
          h_max = c->hydro.h_max;
394
        }
Loic Hausammann's avatar
Loic Hausammann committed
395
        if (c->stars.h_max > h_max) {
396
          h_max = c->stars.h_max;
Loic Hausammann's avatar
Loic Hausammann committed
397
        }
398 399 400
        if (c->black_holes.h_max > h_max) {
          h_max = c->black_holes.h_max;
        }
401
      }
402 403

      /* Can we instead use all the top-level cells? */
404
    } else if (s->cells_top != NULL) {
Tom Theuns's avatar
Tom Theuns committed
405
      for (int k = 0; k < s->nr_cells; k++) {
406
        const struct cell *c = &s->cells_top[k];
407
        if (c->nodeID == engine_rank && c->hydro.h_max > h_max) {
408
          h_max = c->hydro.h_max;
409
        }
Loic Hausammann's avatar
Loic Hausammann committed
410
        if (c->nodeID == engine_rank && c->stars.h_max > h_max) {
411
          h_max = c->stars.h_max;
Loic Hausammann's avatar
Loic Hausammann committed
412
        }
413 414 415
        if (c->nodeID == engine_rank && c->black_holes.h_max > h_max) {
          h_max = c->black_holes.h_max;
        }
416
      }
417 418

      /* Last option: run through the particles */
419
    } else {
420
      for (size_t k = 0; k < nr_parts; k++) {
Matthieu Schaller's avatar
Matthieu Schaller committed
421
        if (s->parts[k].h > h_max) h_max = s->parts[k].h;
422
      }
Loic Hausammann's avatar
Loic Hausammann committed
423 424 425
      for (size_t k = 0; k < nr_sparts; k++) {
        if (s->sparts[k].h > h_max) h_max = s->sparts[k].h;
      }
426 427 428
      for (size_t k = 0; k < nr_bparts; k++) {
        if (s->bparts[k].h > h_max) h_max = s->bparts[k].h;
      }
429 430 431 432 433 434 435 436 437 438
    }
  }

/* 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)
439
      error("Failed to aggregate the rebuild flag across nodes.");
440 441 442
    h_max = buff;
  }
#endif
443
  if (verbose) message("h_max is %.3e (cell_min=%.3e).", h_max, s->cell_min);
444 445

  /* Get the new putative cell dimensions. */
446
  const int cdim[3] = {
447 448 449 450 451 452
      (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))};
453 454 455 456 457

  /* 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 "
458 459 460
        "is switched on.\nThis error is often caused by any of the "
        "followings:\n"
        " - too few particles to generate a sensible grid,\n"
461 462
        " - the initial value of 'Scheduler:max_top_level_cells' is too "
        "small,\n"
463
        " - the (minimal) time-step is too large leading to particles with "
464
        "predicted smoothing lengths too large for the box size,\n"
465
        " - particles with velocities so large that they move by more than two "
466
        "box sizes per time-step.\n");
467

468 469 470
/* 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. */
471
#ifdef WITH_MPI
Matthieu Schaller's avatar
Matthieu Schaller committed
472
  double oldwidth[3];
473 474 475 476 477 478 479 480
  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];
481 482 483
    oldwidth[0] = s->width[0];
    oldwidth[1] = s->width[1];
    oldwidth[2] = s->width[2];
484

485 486
    if ((oldnodeIDs =
             (int *)swift_malloc("nodeIDs", sizeof(int) * s->nr_cells)) == NULL)
487 488 489 490 491 492 493
      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);
494
          oldnodeIDs[cid] = s->cells_top[cid].nodeID;
495 496 497 498 499
        }
      }
    }
  }

Peter W. Draper's avatar
Peter W. Draper committed
500
  /* Are we about to allocate new top level cells without a regrid?
Peter W. Draper's avatar
Peter W. Draper committed
501
   * Can happen when restarting the application. */
502
  const int no_regrid = (s->cells_top == NULL && oldnodeIDs == NULL);
503 504 505 506
#endif

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

510 511
/* Be verbose about this. */
#ifdef SWIFT_DEBUG_CHECKS
512
    message("(re)griding space cdim=(%d %d %d)", cdim[0], cdim[1], cdim[2]);
513 514 515
    fflush(stdout);
#endif

516
    /* Free the old cells, if they were allocated. */
517
    if (s->cells_top != NULL) {
518
      space_free_cells(s);
519 520 521
      swift_free("local_cells_with_tasks_top", s->local_cells_with_tasks_top);
      swift_free("local_cells_top", s->local_cells_top);
      swift_free("cells_with_particles_top", s->cells_with_particles_top);
522 523
      swift_free("local_cells_with_particles_top",
                 s->local_cells_with_particles_top);
524 525
      swift_free("cells_top", s->cells_top);
      swift_free("multipoles_top", s->multipoles_top);
526 527
    }

528 529 530 531
    /* 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);

532
    /* Set the new cell dimensions only if smaller. */
533
    for (int k = 0; k < 3; k++) {
534
      s->cdim[k] = cdim[k];
535 536
      s->width[k] = s->dim[k] / cdim[k];
      s->iwidth[k] = 1.0 / s->width[k];
537
    }
538
    const float dmin = min3(s->width[0], s->width[1], s->width[2]);
539 540 541

    /* Allocate the highest level of cells. */
    s->tot_cells = s->nr_cells = cdim[0] * cdim[1] * cdim[2];
542

543
    if (swift_memalign("cells_top", (void **)&s->cells_top, cell_align,
544
                       s->nr_cells * sizeof(struct cell)) != 0)
545
      error("Failed to allocate top-level cells.");
546
    bzero(s->cells_top, s->nr_cells * sizeof(struct cell));
547

548
    /* Allocate the multipoles for the top-level cells. */
549
    if (s->with_self_gravity) {
550 551
      if (swift_memalign("multipoles_top", (void **)&s->multipoles_top,
                         multipole_align,
552
                         s->nr_cells * sizeof(struct gravity_tensors)) != 0)
553
        error("Failed to allocate top-level multipoles.");
554
      bzero(s->multipoles_top, s->nr_cells * sizeof(struct gravity_tensors));
555 556
    }

557
    /* Allocate the indices of local cells */
558 559
    if (swift_memalign("local_cells_top", (void **)&s->local_cells_top,
                       SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0)
560 561 562
      error("Failed to allocate indices of local top-level cells.");
    bzero(s->local_cells_top, s->nr_cells * sizeof(int));

563
    /* Allocate the indices of local cells with tasks */
564 565
    if (swift_memalign("local_cells_with_tasks_top",
                       (void **)&s->local_cells_with_tasks_top,
566
                       SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0)
567
      error("Failed to allocate indices of local top-level cells with tasks.");
568 569
    bzero(s->local_cells_with_tasks_top, s->nr_cells * sizeof(int));

570
    /* Allocate the indices of cells with particles */
571 572
    if (swift_memalign("cells_with_particles_top",
                       (void **)&s->cells_with_particles_top,
573
                       SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0)
574 575
      error("Failed to allocate indices of top-level cells with particles.");
    bzero(s->cells_with_particles_top, s->nr_cells * sizeof(int));
576

577
    /* Allocate the indices of local cells with particles */
578 579
    if (swift_memalign("local_cells_with_particles_top",
                       (void **)&s->local_cells_with_particles_top,
580 581 582 583 584 585
                       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));

586
    /* Set the cells' locks */
587
    for (int k = 0; k < s->nr_cells; k++) {
588
      if (lock_init(&s->cells_top[k].hydro.lock) != 0)
589
        error("Failed to init spinlock for hydro.");
590
      if (lock_init(&s->cells_top[k].grav.plock) != 0)
591
        error("Failed to init spinlock for gravity.");
592
      if (lock_init(&s->cells_top[k].grav.mlock) != 0)
593
        error("Failed to init spinlock for multipoles.");
Loic Hausammann's avatar
Loic Hausammann committed
594 595
      if (lock_init(&s->cells_top[k].grav.star_formation_lock) != 0)
        error("Failed to init spinlock for star formation (gpart).");
596
      if (lock_init(&s->cells_top[k].stars.lock) != 0)
597
        error("Failed to init spinlock for stars.");
Loic Hausammann's avatar
Loic Hausammann committed
598 599
      if (lock_init(&s->cells_top[k].sinks.lock) != 0)
        error("Failed to init spinlock for sinks.");
Matthieu Schaller's avatar
Matthieu Schaller committed
600 601
      if (lock_init(&s->cells_top[k].black_holes.lock) != 0)
        error("Failed to init spinlock for black holes.");
602
      if (lock_init(&s->cells_top[k].stars.star_formation_lock) != 0)
Loic Hausammann's avatar
Loic Hausammann committed
603
        error("Failed to init spinlock for star formation (spart).");
604
    }
605 606

    /* Set the cell location and sizes. */
607 608 609
    for (int i = 0; i < cdim[0]; i++)
      for (int j = 0; j < cdim[1]; j++)
        for (int k = 0; k < cdim[2]; k++) {
610 611
          const size_t cid = cell_getid(cdim, i, j, k);
          struct cell *restrict c = &s->cells_top[cid];
612 613 614 615 616 617
          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];
618 619
          c->dmin = dmin;
          c->depth = 0;
620
          c->split = 0;
621
          c->hydro.count = 0;
622 623
          c->grav.count = 0;
          c->stars.count = 0;
624
          c->top = c;
625
          c->super = c;
626 627
          c->hydro.super = c;
          c->grav.super = c;
628 629
          c->hydro.ti_old_part = ti_current;
          c->grav.ti_old_part = ti_current;
630
          c->stars.ti_old_part = ti_current;
631
          c->black_holes.ti_old_part = ti_current;
632
          c->grav.ti_old_multipole = ti_current;
633
#ifdef WITH_MPI
634
          c->mpi.tag = -1;
635
          c->mpi.recv = NULL;
636
          c->mpi.send = NULL;
637
#endif  // WITH_MPI
638
          if (s->with_self_gravity) c->grav.multipole = &s->multipoles_top[cid];
Loic Hausammann's avatar
Loic Hausammann committed
639
#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_CELL_GRAPH)
Matthieu Schaller's avatar
Matthieu Schaller committed
640 641
          c->cellID = -last_cell_id;
          last_cell_id++;
642
#endif
Pedro Gonnet's avatar
Pedro Gonnet committed
643
        }
644 645

    /* Be verbose about the change. */
646 647 648
    if (verbose)
      message("set cell dimensions to [ %i %i %i ].", cdim[0], cdim[1],
              cdim[2]);
649

650
#ifdef WITH_MPI
651 652 653 654 655
    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)
656 657 658
        message(
            "basic cell dimensions have increased - recalculating the "
            "global partition.");
659

Matthieu Schaller's avatar
Matthieu Schaller committed
660
      if (!partition_space_to_space(oldwidth, oldcdim, oldnodeIDs, s)) {
661 662 663 664

        /* Failed, try another technique that requires no settings. */
        message("Failed to get a new partition, trying less optimal method");
        struct partition initial_partition;
665
#if defined(HAVE_PARMETIS) || defined(HAVE_METIS)
666 667 668 669 670 671 672 673 674 675 676 677 678
        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);
679

680
      /* Finished with these. */
681
      swift_free("nodeIDs", oldnodeIDs);
Peter W. Draper's avatar
Peter W. Draper committed
682 683

    } else if (no_regrid && s->e != NULL) {
Peter W. Draper's avatar
Peter W. Draper committed
684 685 686 687 688 689 690 691 692 693 694 695 696
      /* 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);
697
    }
Pedro Gonnet's avatar
Pedro Gonnet committed
698
#endif /* WITH_MPI */
699 700 701 702

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

703
  }      /* re-build upper-level cells? */
704
  else { /* Otherwise, just clean up the cells. */
705 706

    /* Free the old cells, if they were allocated. */
707
    space_free_cells(s);
708
  }
709 710 711 712

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

715 716 717 718 719 720
/**
 * @brief Allocate memory for the extra particles used for on-the-fly creation.
 *
 * This rarely actually allocates memory. Most of the time, we convert
 * pre-allocated memory inot extra particles.
 *
721 722
 * This function also sets the extra particles' location to their top-level
 * cells. They can then be sorted into their correct memory position later on.
723 724 725 726
 *
 * @param s The current #space.
 * @param verbose Are we talkative?
 */
727 728 729 730
void space_allocate_extras(struct space *s, int verbose) {

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

731 732
  /* Anything to do here? (Abort if we don't want extras)*/
  if (space_extra_parts == 0 && space_extra_gparts == 0 &&
733
      space_extra_sparts == 0 && space_extra_bparts == 0)
734 735
    return;

736 737 738 739 740 741
  /* The top-level cells */
  const struct cell *cells = s->cells_top;
  const double half_cell_width[3] = {0.5 * cells[0].width[0],
                                     0.5 * cells[0].width[1],
                                     0.5 * cells[0].width[2]};

742
  /* The current number of particles (including spare ones) */
743 744 745
  size_t nr_parts = s->nr_parts;
  size_t nr_gparts = s->nr_gparts;
  size_t nr_sparts = s->nr_sparts;
746
  size_t nr_bparts = s->nr_bparts;
Loic Hausammann's avatar
Loic Hausammann committed
747
  size_t nr_sinks = s->nr_sinks;
748

749 750 751 752
  /* The current number of actual particles */
  size_t nr_actual_parts = nr_parts - s->nr_extra_parts;
  size_t nr_actual_gparts = nr_gparts - s->nr_extra_gparts;
  size_t nr_actual_sparts = nr_sparts - s->nr_extra_sparts;
753
  size_t nr_actual_bparts = nr_bparts - s->nr_extra_bparts;
Loic Hausammann's avatar
Loic Hausammann committed
754
  size_t nr_actual_sinks = nr_sinks - s->nr_extra_sinks;
755

756 757 758 759
  /* The number of particles we allocated memory for (MPI overhead) */
  size_t size_parts = s->size_parts;
  size_t size_gparts = s->size_gparts;
  size_t size_sparts = s->size_sparts;
760
  size_t size_bparts = s->size_bparts;
761

762 763 764 765 766
  int *local_cells = (int *)malloc(sizeof(int) * s->nr_cells);
  if (local_cells == NULL)
    error("Failed to allocate list of local top-level cells");

  /* List the local cells */
767
  size_t nr_local_cells = 0;
768 769 770 771 772 773
  for (int i = 0; i < s->nr_cells; ++i) {
    if (s->cells_top[i].nodeID == local_nodeID) {
      local_cells[nr_local_cells] = i;
      ++nr_local_cells;
    }
  }
774

775
  /* Number of extra particles we want for each type */
776 777 778
  const size_t expected_num_extra_parts = nr_local_cells * space_extra_parts;
  const size_t expected_num_extra_gparts = nr_local_cells * space_extra_gparts;
  const size_t expected_num_extra_sparts = nr_local_cells * space_extra_sparts;
779
  const size_t expected_num_extra_bparts = nr_local_cells * space_extra_bparts;
Loic Hausammann's avatar
Loic Hausammann committed
780
  const size_t expected_num_extra_sinks = nr_local_cells * space_extra_sinks;
781

782
  if (verbose) {
Loic Hausammann's avatar
Loic Hausammann committed
783 784 785 786 787 788
    message("Currently have %zd/%zd/%zd/%zd/%zd real particles.",
            nr_actual_parts, nr_actual_gparts, nr_actual_sinks,
            nr_actual_sparts, nr_actual_bparts);
    message("Currently have %zd/%zd/%zd/%zd/%zd spaces for extra particles.",
            s->nr_extra_parts, s->nr_extra_gparts, s->nr_extra_sinks,
            s->nr_extra_sparts, s->nr_extra_bparts);
789
    message(
Loic Hausammann's avatar
Loic Hausammann committed
790 791
        "Requesting space for future %zd/%zd/%zd/%zd/%zd "
        "part/gpart/sinks/sparts/bparts.",
792
        expected_num_extra_parts, expected_num_extra_gparts,
Loic Hausammann's avatar
Loic Hausammann committed
793 794
        expected_num_extra_sinks, expected_num_extra_sparts,
        expected_num_extra_bparts);
795
  }
796

797 798 799 800 801 802
  if (expected_num_extra_parts < s->nr_extra_parts)
    error("Reduction in top-level cells number not handled.");
  if (expected_num_extra_gparts < s->nr_extra_gparts)
    error("Reduction in top-level cells number not handled.");
  if (expected_num_extra_sparts < s->nr_extra_sparts)
    error("Reduction in top-level cells number not handled.");
803 804
  if (expected_num_extra_bparts < s->nr_extra_bparts)
    error("Reduction in top-level cells number not handled.");
Loic Hausammann's avatar
Loic Hausammann committed
805 806
  if (expected_num_extra_sinks < s->nr_extra_sinks)
    error("Reduction in top-level cells number not handled.");
807

808 809
  /* Do we have enough space for the extra gparts (i.e. we haven't used up any)
   * ? */
Loic Hausammann's avatar
Loic Hausammann committed
810
  if (nr_actual_gparts + expected_num_extra_gparts > nr_gparts) {
811 812 813 814 815 816 817 818 819 820 821 822 823 824 825

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

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

      size_gparts = (nr_actual_gparts + expected_num_extra_gparts) *
                    engine_redistribute_alloc_margin;

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

      /* Create more space for parts */
      struct gpart *gparts_new = NULL;
826
      if (swift_memalign("gparts", (void **)&gparts_new, gpart_align,
827 828 829 830
                         sizeof(struct gpart) * size_gparts) != 0)
        error("Failed to allocate new gpart data");
      const ptrdiff_t delta = gparts_new - s->gparts;
      memcpy(gparts_new, s->gparts, sizeof(struct gpart) * s->size_gparts);
831
      swift_free("gparts", s->gparts);
832 833 834 835 836 837 838 839 840 841 842 843 844 845
      s->gparts = gparts_new;

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

      /* We now need to reset all the part and spart pointers */
      for (size_t i = 0; i < nr_parts; ++i) {
        if (s->parts[i].time_bin != time_bin_not_created)
          s->parts[i].gpart += delta;
      }
      for (size_t i = 0; i < nr_sparts; ++i) {
        if (s->sparts[i].time_bin != time_bin_not_created)
          s->sparts[i].gpart += delta;
      }
Loic Hausammann's avatar
Loic Hausammann committed
846 847 848 849
      for (size_t i = 0; i < nr_bparts; ++i) {
        if (s->bparts[i].time_bin != time_bin_not_created)
          s->bparts[i].gpart += delta;
      }
850 851 852 853 854 855 856 857 858 859 860
    }

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

861
    /* Put the spare particles in their correct cell */
862
    size_t local_cell_id = 0;
863 864
    int current_cell = local_cells[local_cell_id];
    int count_in_cell = 0;
865 866 867 868 869 870 871 872 873 874 875
    size_t count_extra_gparts = 0;
    for (size_t i = 0; i < nr_actual_gparts + expected_num_extra_gparts; ++i) {

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

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

        /* We want the extra particles to be at the centre of their cell */
876 877 878
        s->gparts[i].x[0] = cells[current_cell].loc[0] + half_cell_width[0];
        s->gparts[i].x[1] = cells[current_cell].loc[1] + half_cell_width[1];
        s->gparts[i].x[2] = cells[current_cell].loc[2] + half_cell_width[2];
879 880 881 882 883 884 885
        ++count_in_cell;
        count_extra_gparts++;
      }

      /* Once we have reached the number of extra gpart per cell, we move to the
       * next */
      if (count_in_cell == space_extra_gparts) {
886 887 888 889 890
        ++local_cell_id;

        if (local_cell_id == nr_local_cells) break;

        current_cell = local_cells[local_cell_id];
891 892 893 894 895 896 897 898 899 900 901 902 903 904 905
        count_in_cell = 0;
      }
    }

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

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

906 907
  /* Do we have enough space for the extra parts (i.e. we haven't used up any) ?
   */
908
  if (nr_actual_parts + expected_num_extra_parts > nr_parts) {
909

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

912 913 914 915 916 917 918 919 920 921 922 923
    /* Do we need to reallocate? */
    if (nr_actual_parts + expected_num_extra_parts > size_parts) {

      size_parts = (nr_actual_parts + expected_num_extra_parts) *
                   engine_redistribute_alloc_margin;

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

      /* Create more space for parts */
      struct part *parts_new = NULL;
924
      if (swift_memalign("parts", (void **)&parts_new, part_align,
925 926 927
                         sizeof(struct part) * size_parts) != 0)
        error("Failed to allocate new part data");
      memcpy(parts_new, s->parts, sizeof(struct part) * s->size_parts);
928
      swift_free("parts", s->parts);
929 930 931 932
      s->parts = parts_new;

      /* Same for xparts */
      struct xpart *xparts_new = NULL;
933
      if (swift_memalign("xparts", (void **)&xparts_new, xpart_align,
934 935 936
                         sizeof(struct xpart) * size_parts) != 0)
        error("Failed to allocate new xpart data");
      memcpy(xparts_new, s->xparts, sizeof(struct xpart) * s->size_parts);
937
      swift_free("xparts", s->xparts);
938 939 940 941 942 943
      s->xparts = xparts_new;

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

944
    /* Turn some of the allocated spares into particles we can use */
945 946 947 948 949
    for (size_t i = nr_parts; i < nr_actual_parts + expected_num_extra_parts;
         ++i) {
      bzero(&s->parts[i], sizeof(struct part));
      bzero(&s->xparts[i], sizeof(struct xpart));
      s->parts[i].time_bin = time_bin_not_created;
950
      s->parts[i].id = -42;
951 952
    }

953
    /* Put the spare particles in their correct cell */
954
    size_t local_cell_id = 0;
955 956
    int current_cell = local_cells[local_cell_id];
    int count_in_cell = 0;