From 44f618ee742f08d80b51e26cd5b425ae07fef644 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Mon, 12 Nov 2018 17:33:18 +0100
Subject: [PATCH] Sort the extra particles at the end of their top-level cells.
---
src/cell.c | 29 +++++++++++++++++++
src/cell.h | 3 ++
src/space.c | 82 +++++++++++++++++++++++++++++++++++++++++++----------
src/space.h | 4 +++
4 files changed, 103 insertions(+), 15 deletions(-)
diff --git a/src/cell.c b/src/cell.c
index 1c4d1f9e32..57795efcd3 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -3794,6 +3794,35 @@ void cell_convert_spart_to_gpart(const struct engine *e, struct cell *c,
#endif
}
+void cell_reorder_extra_parts(struct cell *c) {
+
+ struct part *parts = c->hydro.parts;
+ struct xpart *xparts = c->hydro.xparts;
+ const int count_real = c->hydro.count;
+ const int count_total = count_real + space_extra_parts;
+
+ int first_not_extra = count_real;
+
+ /* Find extra particles */
+ for (int i = 0; i < count_real; ++i) {
+ if (parts[i].time_bin == time_bin_not_created) {
+
+ /* Find the first non-extra particle after the end of the
+ real particles */
+ while (parts[first_not_extra].time_bin == time_bin_not_created) {
+ ++first_not_extra;
+ }
+
+ if (first_not_extra >= count_total)
+ error("Looking for extra particles beyond this cell's range!");
+
+ 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!");
+ }
+ }
+}
+
/**
* @brief Can we use the MM interactions fo a given pair of cells?
*
diff --git a/src/cell.h b/src/cell.h
index 5d542e5f2d..009e6522db 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -663,6 +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);
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/space.c b/src/space.c
index a27280d411..6f78ed8674 100644
--- a/src/space.c
+++ b/src/space.c
@@ -767,12 +767,12 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
size_t size_gparts = s->size_gparts;
size_t size_sparts = s->size_sparts;
- /* Number of inhibited particles found on the node */
+ /* Counter for the number of inhibited particles found on the node */
size_t count_inhibited_parts = 0;
size_t count_inhibited_gparts = 0;
size_t count_inhibited_sparts = 0;
- /* Number of extra particles found on the node */
+ /* Counter for the number of extra particles found on the node */
size_t count_extra_parts = 0;
size_t count_extra_gparts = 0;
size_t count_extra_sparts = 0;
@@ -1157,22 +1157,26 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
}
#endif /* SWIFT_DEBUG_CHECKS */
- /* Extract the cell counts from the sorted indices. */
+ /* Extract the cell counts from the sorted indices. Deduct the extra
+ * particles. */
size_t last_index = 0;
h_index[nr_parts] = s->nr_cells; // sentinel.
for (size_t k = 0; k < nr_parts; k++) {
if (h_index[k] < h_index[k + 1]) {
- cells_top[h_index[k]].hydro.count = k - last_index + 1;
+ cells_top[h_index[k]].hydro.count =
+ k - last_index + 1 - space_extra_parts;
last_index = k + 1;
}
}
- /* Extract the cell counts from the sorted indices. */
+ /* Extract the cell counts from the sorted indices. Deduct the extra
+ * particles. */
size_t last_sindex = 0;
s_index[nr_sparts] = s->nr_cells; // sentinel.
for (size_t k = 0; k < nr_sparts; k++) {
if (s_index[k] < s_index[k + 1]) {
- cells_top[s_index[k]].stars.count = k - last_sindex + 1;
+ cells_top[s_index[k]].stars.count =
+ k - last_sindex + 1 - space_extra_sparts;
last_sindex = k + 1;
}
}
@@ -1252,12 +1256,14 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
}
#endif /* SWIFT_DEBUG_CHECKS */
- /* Extract the cell counts from the sorted indices. */
+ /* Extract the cell counts from the sorted indices. Deduct the extra
+ * particles. */
size_t last_gindex = 0;
g_index[nr_gparts] = s->nr_cells;
for (size_t k = 0; k < nr_gparts; k++) {
if (g_index[k] < g_index[k + 1]) {
- cells_top[g_index[k]].grav.count = k - last_gindex + 1;
+ cells_top[g_index[k]].grav.count =
+ k - last_gindex + 1 - space_extra_gparts;
last_gindex = k + 1;
}
}
@@ -1302,10 +1308,10 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
c->hydro.xparts = xfinger;
c->grav.parts = gfinger;
c->stars.parts = sfinger;
- finger = &finger[c->hydro.count];
- xfinger = &xfinger[c->hydro.count];
- gfinger = &gfinger[c->grav.count];
- sfinger = &sfinger[c->stars.count];
+ finger = &finger[c->hydro.count + space_extra_parts];
+ xfinger = &xfinger[c->hydro.count + space_extra_parts];
+ gfinger = &gfinger[c->grav.count + space_extra_gparts];
+ sfinger = &sfinger[c->stars.count + space_extra_sparts];
/* Add this cell to the list of local cells */
s->local_cells_top[s->nr_local_cells] = k;
@@ -1328,8 +1334,12 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
clocks_from_ticks(getticks() - tic2), clocks_getunit());
}
- /* At this point, we have the upper-level cells, old or new. Now make
- sure that the parts in each cell are ok. */
+ /* Re-order the extra particles such that they are at the end of their cell's
+ memory pool. */
+ space_reorder_extras(s, verbose);
+
+ /* At this point, we have the upper-level cells. Now recursively split each
+ cell to get the full AMR grid. */
space_split(s, verbose);
#ifdef SWIFT_DEBUG_CHECKS
@@ -1369,6 +1379,39 @@ void space_split(struct space *s, int verbose) {
clocks_getunit());
}
+void space_reorder_extra_parts_mapper(void *map_data, int num_cells,
+ void *extra_data) {
+
+ struct cell *cells_top = (struct cell *)map_data;
+ for (int ind = 0; ind < num_cells; ind++) {
+ struct cell *c = &cells_top[ind];
+ cell_reorder_extra_parts(c);
+ }
+}
+
+/**
+ * @brief Re-orders the particles in each cell such that the extra particles
+ * for on-the-fly creation are located at the end of their respective cells.
+ *
+ * This assumes that all the particles (real and extra) have already been sorted
+ * in their correct top-level cell.
+ *
+ * @param s The #space to act upon.
+ * @param verbose Are we talkative?
+ */
+void space_reorder_extras(struct space *s, int verbose) {
+
+#ifdef WITH_MPI
+ if (space_extra_parts || space_extra_gparts || space_extra_sparts)
+ error("Need an MPI-proof version of this.");
+#endif
+
+ /* 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);
+}
+
/**
* @brief #threadpool mapper function to sanitize the cells
*
@@ -2250,6 +2293,8 @@ void space_split_recursive(struct space *s, struct cell *c,
#ifdef SWIFT_DEBUG_CHECKS
if (parts[k].time_bin == time_bin_inhibited)
error("Inhibited particle present in space_split()");
+ if (parts[k].time_bin == time_bin_not_created)
+ error("Extra particle present in space_split()");
#endif
buff[k].x[0] = parts[k].x[0];
buff[k].x[1] = parts[k].x[1];
@@ -2264,6 +2309,8 @@ void space_split_recursive(struct space *s, struct cell *c,
#ifdef SWIFT_DEBUG_CHECKS
if (gparts[k].time_bin == time_bin_inhibited)
error("Inhibited particle present in space_split()");
+ if (gparts[k].time_bin == time_bin_not_created)
+ error("Extra particle present in space_split()");
#endif
gbuff[k].x[0] = gparts[k].x[0];
gbuff[k].x[1] = gparts[k].x[1];
@@ -2278,6 +2325,8 @@ void space_split_recursive(struct space *s, struct cell *c,
#ifdef SWIFT_DEBUG_CHECKS
if (sparts[k].time_bin == time_bin_inhibited)
error("Inhibited particle present in space_split()");
+ if (sparts[k].time_bin == time_bin_not_created)
+ error("Extra particle present in space_split()");
#endif
sbuff[k].x[0] = sparts[k].x[0];
sbuff[k].x[1] = sparts[k].x[1];
@@ -3234,8 +3283,11 @@ void space_convert_quantities_mapper(void *restrict map_data, int count,
const ptrdiff_t index = parts - s->parts;
struct xpart *restrict xparts = s->xparts + index;
+ /* Loop over all the particles ignoring the extra buffer ones for on-the-fly
+ * cretion */
for (int k = 0; k < count; k++)
- hydro_convert_quantities(&parts[k], &xparts[k], cosmo, hydro_props);
+ if (parts[k].time_bin <= num_time_bins)
+ hydro_convert_quantities(&parts[k], &xparts[k], cosmo, hydro_props);
}
/**
diff --git a/src/space.h b/src/space.h
index cdf83e4ada..4f1a39c287 100644
--- a/src/space.h
+++ b/src/space.h
@@ -72,6 +72,9 @@ extern int space_subsize_self_grav;
extern int space_subsize_pair_stars;
extern int space_subsize_self_stars;
extern int space_subdepth_grav;
+extern int space_extra_parts;
+extern int space_extra_gparts;
+extern int space_extra_sparts;
/**
* @brief The space in which the cells and particles reside.
@@ -282,6 +285,7 @@ void space_recycle_list(struct space *s, struct cell *cell_list_begin,
struct gravity_tensors *multipole_list_begin,
struct gravity_tensors *multipole_list_end);
void space_split(struct space *s, int verbose);
+void space_reorder_extras(struct space *s, int verbose);
void space_split_mapper(void *map_data, int num_elements, void *extra_data);
void space_list_useful_top_level_cells(struct space *s);
void space_parts_get_cell_index(struct space *s, int *ind, int *cell_counts,
--
GitLab