diff --git a/src/cell.c b/src/cell.c
index 67d72bf6205bafdddcdd7db7b53f5e41311a13e0..4c5064e2036d2682be7f0ec52845935ae6b5d4c9 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -3864,11 +3864,12 @@ void cell_reorder_extra_sparts(struct cell *c, const ptrdiff_t sparts_offset) {
}
}
-void cell_reorder_extra_gparts(struct cell *c, const ptrdiff_t sparts_offset) {
+void cell_reorder_extra_gparts(struct cell *c, struct part *parts,
+ struct spart *sparts) {
struct gpart *gparts = c->grav.parts;
const int count_real = c->grav.count;
- const int count_total = count_real + space_extra_parts;
+ const int count_total = count_real + space_extra_gparts;
if (c->depth != 0)
error("This function should only be called on top-level cells!");
@@ -3891,11 +3892,11 @@ void cell_reorder_extra_gparts(struct cell *c, const ptrdiff_t sparts_offset) {
#endif
/* Swap everything (including pointers) */
- memswap(&gparts[i], &gparts[first_not_extra], sizeof(struct spart));
+ memswap(&gparts[i], &gparts[first_not_extra], sizeof(struct gpart));
if (gparts[i].type == swift_type_gas) {
- error("Need to handle this.");
+ parts[-gparts[i].id_or_neg_offset].gpart = &gparts[i];
} else if (gparts[i].type == swift_type_stars) {
- error("Need to handle this.");
+ sparts[-gparts[i].id_or_neg_offset].gpart = &gparts[i];
}
}
}
diff --git a/src/cell.h b/src/cell.h
index cf286f02f11179c62e4f7e6fb6afab74324da609..b6bdf1561bcc46a0aef09bd82ba019cfe3455087 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -664,7 +664,8 @@ void cell_convert_part_to_gpart(const struct engine *e, struct cell *c,
void cell_convert_spart_to_gpart(const struct engine *e, struct cell *c,
struct spart *sp);
void cell_reorder_extra_parts(struct cell *c, const ptrdiff_t parts_offset);
-void cell_reorder_extra_gparts(struct cell *c, const ptrdiff_t gparts_offset);
+void cell_reorder_extra_gparts(struct cell *c, struct part *parts,
+ struct spart *sparts);
void cell_reorder_extra_sparts(struct cell *c, const ptrdiff_t sparts_offset);
int cell_can_use_pair_mm(const struct cell *ci, const struct cell *cj,
const struct engine *e, const struct space *s);
diff --git a/src/engine.c b/src/engine.c
index 9239d00c27a4e29faecf730227741d284ef8d468..0530ae3fb9e48baaa8bd3645563e2c3b95a4d2c1 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1992,9 +1992,10 @@ void engine_rebuild(struct engine *e, int repartitioned,
space_rebuild(e->s, repartitioned, e->verbose);
/* Update the global counters of particles */
- long long num_particles[3] = {(long long)e->s->nr_parts,
- (long long)e->s->nr_gparts,
- (long long)e->s->nr_sparts};
+ long long num_particles[3] = {
+ (long long)(e->s->nr_parts - e->s->nr_extra_parts),
+ (long long)(e->s->nr_gparts - e->s->nr_extra_gparts),
+ (long long)(e->s->nr_sparts - e->s->nr_extra_sparts)};
#ifdef WITH_MPI
MPI_Allreduce(MPI_IN_PLACE, num_particles, 3, MPI_LONG_LONG, MPI_SUM,
MPI_COMM_WORLD);
diff --git a/src/part.c b/src/part.c
index 3a626e652cf28f0376cadc1d9a40ab85b752e6c1..ec3627d728f69f469cc7d75eb2beb9ae39ed107e 100644
--- a/src/part.c
+++ b/src/part.c
@@ -139,8 +139,9 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
for (size_t k = 0; k < nr_gparts; ++k) {
- /* We have a DM particle */
- if (gparts[k].type == swift_type_dark_matter) {
+ /* We have a real DM particle */
+ if (gparts[k].type == swift_type_dark_matter &&
+ gparts[k].time_bin != time_bin_not_created) {
/* Check that it's not linked */
if (gparts[k].id_or_neg_offset <= 0)
diff --git a/src/space.c b/src/space.c
index 98f93d0d0d6810e600de33690a08662547519677..88da0b5fed5511b1151d3ce7e9810c336e7bf5b2 100644
--- a/src/space.c
+++ b/src/space.c
@@ -669,10 +669,103 @@ void space_allocate_extras(struct space *s, int verbose) {
if (expected_num_extra_sparts < s->nr_extra_sparts)
error("Reduction in top-level cells number not handled.");
+ /* Do we have enough space for the extra gparts (i.e. we haven't used up any)
+ * ? */
+ if (nr_gparts + expected_num_extra_gparts > size_gparts) {
+
+ /* Ok... need to put some more in the game */
+
+ /* Do we need to reallocate? */
+ if (nr_actual_gparts + expected_num_extra_gparts > size_gparts) {
+
+ size_gparts = (nr_actual_gparts + expected_num_extra_gparts) *
+ engine_redistribute_alloc_margin;
+
+ if (verbose)
+ message("Re-allocating gparts array from %zd to %zd", s->size_gparts,
+ size_gparts);
+
+ /* Create more space for parts */
+ struct gpart *gparts_new = NULL;
+ if (posix_memalign((void **)&gparts_new, gpart_align,
+ sizeof(struct gpart) * size_gparts) != 0)
+ error("Failed to allocate new gpart data");
+ const ptrdiff_t delta = gparts_new - s->gparts;
+ memcpy(gparts_new, s->gparts, sizeof(struct gpart) * s->size_gparts);
+ free(s->gparts);
+ s->gparts = gparts_new;
+
+ /* Update the counter */
+ s->size_gparts = size_gparts;
+
+ /* We now need to reset all the part and spart pointers */
+ for (size_t i = 0; i < nr_parts; ++i) {
+ if (s->parts[i].time_bin != time_bin_not_created)
+ s->parts[i].gpart += delta;
+ }
+ for (size_t i = 0; i < nr_sparts; ++i) {
+ if (s->sparts[i].time_bin != time_bin_not_created)
+ s->sparts[i].gpart += delta;
+ }
+ }
+
+ /* Turn some of the allocated spares into particles we can use */
+ for (size_t i = nr_gparts; i < nr_actual_gparts + expected_num_extra_gparts;
+ ++i) {
+ bzero(&s->gparts[i], sizeof(struct gpart));
+ s->gparts[i].time_bin = time_bin_not_created;
+ s->gparts[i].type = swift_type_dark_matter;
+ s->gparts[i].id_or_neg_offset = -1;
+ }
+
+ /* Put the spare particles in their correct cell */
+#ifdef WITH_MPI
+ error("Need to do this correctly over MPI for only the local cells.");
+#endif
+ int count_in_cell = 0, current_cell = 0;
+ size_t count_extra_gparts = 0;
+ for (size_t i = 0; i < nr_actual_gparts + expected_num_extra_gparts; ++i) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (current_cell == s->nr_cells)
+ error("Cell counter beyond the maximal nr. cells.");
+#endif
+
+ if (s->gparts[i].time_bin == time_bin_not_created) {
+
+ /* We want the extra particles to be at the centre of their cell */
+ s->gparts[i].x[0] = cells[current_cell].loc[0] + half_cell_width[0];
+ s->gparts[i].x[1] = cells[current_cell].loc[1] + half_cell_width[1];
+ s->gparts[i].x[2] = cells[current_cell].loc[2] + half_cell_width[2];
+ ++count_in_cell;
+ count_extra_gparts++;
+ }
+
+ /* Once we have reached the number of extra gpart per cell, we move to the
+ * next */
+ if (count_in_cell == space_extra_gparts) {
+ ++current_cell;
+ count_in_cell = 0;
+ }
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (count_extra_gparts != expected_num_extra_gparts)
+ error("Constructed the wrong number of extra gparts (%zd vs. %zd)",
+ count_extra_gparts, expected_num_extra_gparts);
+#endif
+
+ /* Update the counters */
+ s->nr_gparts = nr_actual_gparts + expected_num_extra_gparts;
+ s->nr_extra_gparts = expected_num_extra_gparts;
+ }
+
/* Do we have enough space for the extra parts (i.e. we haven't used up any) ?
*/
if (expected_num_extra_parts > s->nr_extra_parts) {
+ /* Ok... need to put some more in the game */
+
/* Do we need to reallocate? */
if (nr_actual_parts + expected_num_extra_parts > size_parts) {
@@ -711,6 +804,7 @@ void space_allocate_extras(struct space *s, int verbose) {
bzero(&s->parts[i], sizeof(struct part));
bzero(&s->xparts[i], sizeof(struct xpart));
s->parts[i].time_bin = time_bin_not_created;
+ s->parts[i].id = -1;
}
/* Put the spare particles in their correct cell */
@@ -755,13 +849,12 @@ void space_allocate_extras(struct space *s, int verbose) {
s->nr_extra_parts = expected_num_extra_parts;
}
- if (nr_gparts + expected_num_extra_gparts > size_gparts) {
- }
-
/* Do we have enough space for the extra sparts (i.e. we haven't used up any)
* ? */
if (nr_sparts + expected_num_extra_sparts > size_sparts) {
+ /* Ok... need to put some more in the game */
+
/* Do we need to reallocate? */
if (nr_actual_sparts + expected_num_extra_sparts > size_sparts) {
@@ -790,6 +883,7 @@ void space_allocate_extras(struct space *s, int verbose) {
++i) {
bzero(&s->sparts[i], sizeof(struct spart));
s->sparts[i].time_bin = time_bin_not_created;
+ s->sparts[i].id = -1;
}
/* Put the spare particles in their correct cell */
@@ -833,6 +927,13 @@ void space_allocate_extras(struct space *s, int verbose) {
s->nr_sparts = nr_actual_sparts + expected_num_extra_sparts;
s->nr_extra_sparts = expected_num_extra_sparts;
}
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Verify that the links are correct */
+ if ((nr_gparts > 0 && nr_parts > 0) || (nr_gparts > 0 && nr_sparts > 0))
+ part_verify_links(s->parts, s->gparts, s->sparts, nr_parts, nr_gparts,
+ nr_sparts, verbose);
+#endif
}
/**
@@ -1504,7 +1605,7 @@ void space_reorder_extra_gparts_mapper(void *map_data, int num_cells,
for (int ind = 0; ind < num_cells; ind++) {
struct cell *c = &cells_top[ind];
- cell_reorder_extra_gparts(c, c->grav.parts - s->gparts);
+ cell_reorder_extra_gparts(c, s->parts, s->sparts);
}
}
diff --git a/src/space.h b/src/space.h
index 4f1a39c287b491a727f0775e2812a2475887cee3..efb3ce407f23463bc72c82eaf39e750683c63c9b 100644
--- a/src/space.h
+++ b/src/space.h
@@ -45,8 +45,8 @@ struct cosmology;
#define space_splitsize_default 400
#define space_maxsize_default 8000000
#define space_extra_parts_default 10
-#define space_extra_gparts_default 0
-#define space_extra_sparts_default 0
+#define space_extra_gparts_default 30
+#define space_extra_sparts_default 20
#define space_expected_max_nr_strays_default 100
#define space_subsize_pair_hydro_default 256000000
#define space_subsize_self_hydro_default 32000