diff --git a/src/cell.c b/src/cell.c
index bcf8b064b0f214b5b91fe6f454489c9b5ca8cf2a..141c9a4079a3f00d5deb951ebacf161676baf1de 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -3794,7 +3794,7 @@ void cell_convert_spart_to_gpart(const struct engine *e, struct cell *c,
#endif
}
-void cell_reorder_extra_parts(struct cell *c) {
+void cell_reorder_extra_parts(struct cell *c, const ptrdiff_t parts_offset) {
struct part *parts = c->hydro.parts;
struct xpart *xparts = c->hydro.xparts;
@@ -3816,12 +3816,87 @@ void cell_reorder_extra_parts(struct cell *c) {
++first_not_extra;
}
+#ifdef SWIFT_DEBUG_CHECKS
if (first_not_extra >= count_total)
error("Looking for extra particles beyond this cell's range!");
+#endif
+ /* Swap everything, including g-part pointer */
memswap(&parts[i], &parts[first_not_extra], sizeof(struct part));
memswap(&xparts[i], &xparts[first_not_extra], sizeof(struct xpart));
- if (parts[i].gpart) error("Need to handle gpart pointer!");
+ if (parts[i].gpart)
+ parts[i].gpart->id_or_neg_offset = -(i + parts_offset);
+ }
+ }
+}
+
+void cell_reorder_extra_sparts(struct cell *c, const ptrdiff_t sparts_offset) {
+
+ struct spart *sparts = c->stars.parts;
+ const int count_real = c->stars.count;
+ const int count_total = count_real + space_extra_parts;
+
+ if (c->depth != 0)
+ error("This function should only be called on top-level cells!");
+
+ int first_not_extra = count_real;
+
+ /* Find extra particles */
+ for (int i = 0; i < count_real; ++i) {
+ if (sparts[i].time_bin == time_bin_not_created) {
+
+ /* Find the first non-extra particle after the end of the
+ real particles */
+ while (sparts[first_not_extra].time_bin == time_bin_not_created) {
+ ++first_not_extra;
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (first_not_extra >= count_total)
+ error("Looking for extra particles beyond this cell's range!");
+#endif
+
+ /* Swap everything, including g-part pointer */
+ memswap(&sparts[i], &sparts[first_not_extra], sizeof(struct spart));
+ if (sparts[i].gpart)
+ sparts[i].gpart->id_or_neg_offset = -(i + sparts_offset);
+ }
+ }
+}
+
+void cell_reorder_extra_gparts(struct cell *c, const ptrdiff_t sparts_offset) {
+
+ struct gpart *gparts = c->grav.parts;
+ const int count_real = c->grav.count;
+ const int count_total = count_real + space_extra_parts;
+
+ if (c->depth != 0)
+ error("This function should only be called on top-level cells!");
+
+ int first_not_extra = count_real;
+
+ /* Find extra particles */
+ for (int i = 0; i < count_real; ++i) {
+ if (gparts[i].time_bin == time_bin_not_created) {
+
+ /* Find the first non-extra particle after the end of the
+ real particles */
+ while (gparts[first_not_extra].time_bin == time_bin_not_created) {
+ ++first_not_extra;
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (first_not_extra >= count_total)
+ error("Looking for extra particles beyond this cell's range!");
+#endif
+
+ /* Swap everything (including pointers) */
+ memswap(&gparts[i], &gparts[first_not_extra], sizeof(struct spart));
+ if (gparts[i].type == swift_type_gas) {
+ error("Need to handle this.");
+ } else if (gparts[i].type == swift_type_stars) {
+ error("Need to handle this.");
+ }
}
}
}
diff --git a/src/cell.h b/src/cell.h
index 009e6522db503a1c629e3cbefddbcf5b8a7220c9..cf286f02f11179c62e4f7e6fb6afab74324da609 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -663,9 +663,9 @@ void cell_convert_part_to_gpart(const struct engine *e, struct cell *c,
struct part *p, struct xpart *xp);
void cell_convert_spart_to_gpart(const struct engine *e, struct cell *c,
struct spart *sp);
-void cell_reorder_extra_parts(struct cell *c);
-void cell_reorder_extra_gparts(struct cell *c);
-void cell_reorder_extra_sparts(struct cell *c);
+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_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);
int cell_can_use_pair_mm_rebuild(const struct cell *ci, const struct cell *cj,
diff --git a/src/engine.c b/src/engine.c
index 9f2a994ea487e30c23df3afc9b983e4f590636a9..9239d00c27a4e29faecf730227741d284ef8d468 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -5143,9 +5143,11 @@ void engine_recompute_displacement_constraint(struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check that the minimal mass collection worked */
float min_part_mass_check = FLT_MAX;
- for (size_t i = 0; i < e->s->nr_parts; ++i)
+ for (size_t i = 0; i < e->s->nr_parts; ++i) {
+ if (e->s->parts[i].time_bin >= num_time_bins) continue;
min_part_mass_check =
min(min_part_mass_check, hydro_get_mass(&e->s->parts[i]));
+ }
if (min_part_mass_check != min_mass[swift_type_gas])
error("Error collecting minimal mass of gas particles.");
#endif
diff --git a/src/space.c b/src/space.c
index d6895fb0fcf83e00bd34534d2e59d36f01d4ba60..9f313a36176ebd610d3c5f112881d5fb5f4c82c7 100644
--- a/src/space.c
+++ b/src/space.c
@@ -607,8 +607,8 @@ void space_regrid(struct space *s, int verbose) {
* This rarely actually allocates memory. Most of the time, we convert
* pre-allocated memory inot extra particles.
*
- * 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.
+ * 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.
*
* @param s The current #space.
* @param verbose Are we talkative?
@@ -617,6 +617,11 @@ void space_allocate_extras(struct space *s, int verbose) {
const int local_nodeID = s->e->nodeID;
+ /* Anything to do here? (Abort if we don't want extras)*/
+ if (space_extra_parts == 0 && space_extra_gparts == 0 &&
+ space_extra_sparts == 0)
+ return;
+
/* The top-level cells */
const struct cell *cells = s->cells_top;
const double half_cell_width[3] = {0.5 * cells[0].width[0],
@@ -650,7 +655,7 @@ void space_allocate_extras(struct space *s, int verbose) {
if (verbose) {
message("Currently have %zd/%zd/%zd real particles.", nr_actual_parts,
nr_actual_gparts, nr_actual_sparts);
- message("Currently have %zd/%zd/%zd spaces for future particles.",
+ message("Currently have %zd/%zd/%zd spaces for extra particles.",
s->nr_extra_parts, s->nr_extra_gparts, s->nr_extra_sparts);
message("Requesting space for future %zd/%zd/%zd part/gpart/sparts.",
expected_num_extra_parts, expected_num_extra_gparts,
@@ -707,6 +712,12 @@ void space_allocate_extras(struct space *s, int verbose) {
int count_in_cell = 0, current_cell = 0;
size_t count_extra_parts = 0;
for (size_t i = 0; i < nr_actual_parts + expected_num_extra_parts; ++i) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (current_cell == s->nr_cells)
+ error("Cell counter beyond the maximal nr. cells.");
+#endif
+
if (s->parts[i].time_bin == time_bin_not_created) {
/* We want the extra particles to be at the centre of their cell */
@@ -736,9 +747,8 @@ void space_allocate_extras(struct space *s, int verbose) {
s->nr_extra_parts = expected_num_extra_parts;
}
- /* Did we reduce the number of top-level cells ? */
if (expected_num_extra_parts < s->nr_extra_parts) {
- error("Un-handled case!");
+ error("Reduction in top-level cells number not handled.");
}
if (nr_gparts + expected_num_extra_gparts > size_gparts) {
@@ -1400,9 +1410,35 @@ void space_reorder_extra_parts_mapper(void *map_data, int num_cells,
void *extra_data) {
struct cell *cells_top = (struct cell *)map_data;
+ struct space *s = (struct space *)extra_data;
+
for (int ind = 0; ind < num_cells; ind++) {
struct cell *c = &cells_top[ind];
- cell_reorder_extra_parts(c);
+ cell_reorder_extra_parts(c, c->hydro.parts - s->parts);
+ }
+}
+
+void space_reorder_extra_gparts_mapper(void *map_data, int num_cells,
+ void *extra_data) {
+
+ struct cell *cells_top = (struct cell *)map_data;
+ struct space *s = (struct space *)extra_data;
+
+ for (int ind = 0; ind < num_cells; ind++) {
+ struct cell *c = &cells_top[ind];
+ cell_reorder_extra_gparts(c, c->grav.parts - s->gparts);
+ }
+}
+
+void space_reorder_extra_sparts_mapper(void *map_data, int num_cells,
+ void *extra_data) {
+
+ struct cell *cells_top = (struct cell *)map_data;
+ struct space *s = (struct space *)extra_data;
+
+ for (int ind = 0; ind < num_cells; ind++) {
+ struct cell *c = &cells_top[ind];
+ cell_reorder_extra_sparts(c, c->stars.parts - s->sparts);
}
}
@@ -1426,7 +1462,17 @@ void space_reorder_extras(struct space *s, int verbose) {
/* Re-order the gas particles */
if (space_extra_parts)
threadpool_map(&s->e->threadpool, space_reorder_extra_parts_mapper,
- s->cells_top, s->nr_cells, sizeof(struct cell), 0, NULL);
+ s->cells_top, s->nr_cells, sizeof(struct cell), 0, s);
+
+ /* Re-order the gravity particles */
+ if (space_extra_gparts)
+ threadpool_map(&s->e->threadpool, space_reorder_extra_gparts_mapper,
+ s->cells_top, s->nr_cells, sizeof(struct cell), 0, s);
+
+ /* Re-order the star particles */
+ if (space_extra_sparts)
+ threadpool_map(&s->e->threadpool, space_reorder_extra_sparts_mapper,
+ s->cells_top, s->nr_cells, sizeof(struct cell), 0, s);
}
/**
@@ -1538,13 +1584,13 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts,
/* Normal case: list its top-level cell index */
ind[k] = index;
cell_counts[index]++;
- }
- /* Compute minimal mass */
- min_mass = min(min_mass, hydro_get_mass(p));
+ /* Compute minimal mass */
+ min_mass = min(min_mass, hydro_get_mass(p));
- /* Compute sum of velocity norm */
- sum_vel_norm += p->v[0] * p->v[0] + p->v[1] * p->v[1] + p->v[2] * p->v[2];
+ /* Compute sum of velocity norm */
+ sum_vel_norm += p->v[0] * p->v[0] + p->v[1] * p->v[1] + p->v[2] * p->v[2];
+ }
/* Update the position */
p->x[0] = pos_x;
@@ -1645,17 +1691,17 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts,
/* List its top-level cell index */
ind[k] = index;
cell_counts[index]++;
- }
- if (gp->type == swift_type_dark_matter) {
+ if (gp->type == swift_type_dark_matter) {
- /* Compute minimal mass */
- min_mass = min(min_mass, gp->mass);
+ /* Compute minimal mass */
+ min_mass = min(min_mass, gp->mass);
- /* Compute sum of velocity norm */
- sum_vel_norm += gp->v_full[0] * gp->v_full[0] +
- gp->v_full[1] * gp->v_full[1] +
- gp->v_full[2] * gp->v_full[2];
+ /* Compute sum of velocity norm */
+ sum_vel_norm += gp->v_full[0] * gp->v_full[0] +
+ gp->v_full[1] * gp->v_full[1] +
+ gp->v_full[2] * gp->v_full[2];
+ }
}
/* Update the position */
@@ -1753,14 +1799,14 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts,
/* List its top-level cell index */
ind[k] = index;
cell_counts[index]++;
- }
- /* Compute minimal mass */
- min_mass = min(min_mass, sp->mass);
+ /* Compute minimal mass */
+ min_mass = min(min_mass, sp->mass);
- /* Compute sum of velocity norm */
- sum_vel_norm +=
- sp->v[0] * sp->v[0] + sp->v[1] * sp->v[1] + sp->v[2] * sp->v[2];
+ /* Compute sum of velocity norm */
+ sum_vel_norm +=
+ sp->v[0] * sp->v[0] + sp->v[1] * sp->v[1] + sp->v[2] * sp->v[2];
+ }
/* Update the position */
sp->x[0] = pos_x;