diff --git a/src/cell.c b/src/cell.c
index 7ce6fb81a8fa6875884d3f5c840c36e5177cdf6b..2c57f09c46f563fab91ba3c69b4a6c27e4ee5478 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -679,6 +679,55 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset) {
part_relink_parts(gparts, gcount, parts - parts_offset);
}
+/**
+ * @brief Sanitizes the smoothing length values of cells by setting large
+ * outliers to more sensible values.
+ *
+ * We compute the mean and standard deviation of the smoothing lengths in
+ * logarithmic space and limit values to mean + 4 sigma.
+ *
+ * @param c The cell.
+ */
+void cell_sanitize(struct cell *c) {
+
+ const int count = c->count;
+ struct part *parts = c->parts;
+
+ /* First collect some statistics */
+ float h_mean = 0.f, h_mean2 = 0.f;
+ float h_min = FLT_MAX, h_max = 0.f;
+ for (int i = 0; i < count; ++i) {
+
+ const float h = log10f(parts[i].h);
+ h_mean += h;
+ h_mean2 += h * h;
+ h_max = max(h_max, h);
+ h_min = min(h_min, h);
+ }
+ h_mean /= count;
+ h_mean2 /= count;
+ const float h_var = h_mean2 - h_mean * h_mean;
+ const float h_std = sqrtf(h_var);
+
+ /* Choose a cut */
+ const float h_limit = pow(10.f, h_mean + 4.f * h_std);
+
+ /* Be verbose this is not innocuous */
+ message("Cell properties: h_min= %f h_max= %f geometric mean= %f",
+ powf(10.f, h_min), powf(10.f, h_max), powf(10.f, h_mean));
+
+ if (c->h_max > h_limit) {
+
+ message("Smoothing lengths will be limited to (mean + 4sigma)= %f",
+ h_limit);
+
+ /* Apply the cut */
+ for (int i = 0; i < count; ++i) parts->h = min(parts[i].h, h_limit);
+
+ c->h_max = h_limit;
+ }
+}
+
/**
* @brief Converts hydro quantities to a valid state after the initial density
* calculation
diff --git a/src/cell.h b/src/cell.h
index b78cc0a8f842770f60777e3986616a175d2f33ca..f86ab765e75052c9231fdfebe225da0a5c72be68 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -218,6 +218,7 @@ struct cell {
/* Function prototypes. */
void cell_split(struct cell *c, ptrdiff_t parts_offset);
+void cell_sanitize(struct cell *c);
int cell_locktree(struct cell *c);
void cell_unlocktree(struct cell *c);
int cell_glocktree(struct cell *c);
diff --git a/src/engine.c b/src/engine.c
index b73ad2e24f91426b87c0933f7987de9ff88991c9..38e655d0a820a046a1821f323d191db672d3875c 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2243,6 +2243,8 @@ void engine_rebuild(struct engine *e) {
/* Re-build the space. */
space_rebuild(e->s, 0.0, e->verbose);
+ if (e->ti_current == 0) space_sanitize(e->s);
+
/* If in parallel, exchange the cell structure. */
#ifdef WITH_MPI
engine_exchange_cells(e);
diff --git a/src/space.c b/src/space.c
index 6a5e9af4811029390605f6da5050879c0c64690c..6baaee5200eec81f43375f7296723d8919eeb1eb 100644
--- a/src/space.c
+++ b/src/space.c
@@ -58,6 +58,7 @@
int space_splitsize = space_splitsize_default;
int space_subsize = space_subsize_default;
int space_maxsize = space_maxsize_default;
+int space_maxcount = space_maxcount_default;
/* Map shift vector to sortlist. */
const int sortlistID[27] = {
@@ -698,7 +699,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
/* At this point, we have the upper-level cells, old or new. Now make
sure that the parts in each cell are ok. */
- space_split(s, cells_top, verbose);
+ space_split(s, cells_top, s->nr_cells, verbose);
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -711,14 +712,16 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
* This is done in parallel using threads in the #threadpool.
*
* @param s The #space.
- * @param cells The cell hierarchy
+ * @param cells The cell hierarchy.
+ * @param nr_cells The number of cells.
* @param verbose Are we talkative ?
*/
-void space_split(struct space *s, struct cell *cells, int verbose) {
+void space_split(struct space *s, struct cell *cells, int nr_cells,
+ int verbose) {
const ticks tic = getticks();
- threadpool_map(&s->e->threadpool, space_split_mapper, cells, s->nr_cells,
+ threadpool_map(&s->e->threadpool, space_split_mapper, cells, nr_cells,
sizeof(struct cell), 1, s);
if (verbose)
@@ -726,6 +729,29 @@ void space_split(struct space *s, struct cell *cells, int verbose) {
clocks_getunit());
}
+/**
+ * @brief Runs through the top-level cells and checks whether tasks associated
+ * with them can be split. If not, try to sanitize the cells.
+ *
+ * @param s The #space to act upon.
+ */
+void space_sanitize(struct space *s) {
+
+ for (int k = 0; k < s->nr_cells; k++) {
+
+ struct cell *c = &s->cells_top[k];
+ const double min_width = c->dmin;
+
+ /* Do we have a problem ? */
+ if (c->h_max * kernel_gamma * space_stretch > min_width * 0.5 &&
+ c->count > space_maxcount) {
+
+ /* Ok, clean-up the mess */
+ cell_sanitize(c);
+ }
+ }
+}
+
/**
* @brief Sort the particles and condensed particles according to the given
* indices.
@@ -1247,17 +1273,17 @@ void space_map_cells_pre(struct space *s, int full,
* too many particles.
*
* @param map_data Pointer towards the top-cells.
- * @param num_elements The number of cells to treat.
+ * @param num_cells The number of cells to treat.
* @param extra_data Pointers to the #space.
*/
-void space_split_mapper(void *map_data, int num_elements, void *extra_data) {
+void space_split_mapper(void *map_data, int num_cells, void *extra_data) {
/* Unpack the inputs. */
struct space *s = (struct space *)extra_data;
struct cell *restrict cells_top = (struct cell *)map_data;
struct engine *e = s->e;
- for (int ind = 0; ind < num_elements; ind++) {
+ for (int ind = 0; ind < num_cells; ind++) {
struct cell *c = &cells_top[ind];
@@ -1578,6 +1604,8 @@ void space_init(struct space *s, const struct swift_params *params,
space_subsize_default);
space_splitsize = parser_get_opt_param_int(
params, "Scheduler:cell_split_size", space_splitsize_default);
+ space_maxcount = parser_get_opt_param_int(params, "Scheduler:cell_max_count",
+ space_maxcount_default);
if (verbose)
message("max_size set to %d, sub_size set to %d, split_size set to %d",
space_maxsize, space_subsize, space_splitsize);
diff --git a/src/space.h b/src/space.h
index 5332bc9fb5a7142f251d3ffd28dde6fac5211ea5..9cc8ea56140941241f4df2e55562289ba519f56b 100644
--- a/src/space.h
+++ b/src/space.h
@@ -41,6 +41,7 @@
#define space_splitsize_default 400
#define space_maxsize_default 8000000
#define space_subsize_default 64000000
+#define space_maxcount_default 100000
#define space_stretch 1.10f
#define space_maxreldx 0.25f
@@ -51,6 +52,7 @@
extern int space_splitsize;
extern int space_maxsize;
extern int space_subsize;
+extern int space_maxcount;
/* Map shift vector to sortlist. */
extern const int sortlistID[27];
@@ -148,6 +150,7 @@ void space_init(struct space *s, const struct swift_params *params,
double dim[3], struct part *parts, struct gpart *gparts,
size_t Npart, size_t Ngpart, int periodic, int gravity,
int verbose, int dry_run);
+void space_sanitize(struct space *s);
void space_map_cells_pre(struct space *s, int full,
void (*fun)(struct cell *c, void *data), void *data);
void space_map_parts(struct space *s,
@@ -164,7 +167,8 @@ void space_gparts_sort_mapper(void *map_data, int num_elements,
void *extra_data);
void space_rebuild(struct space *s, double h_max, int verbose);
void space_recycle(struct space *s, struct cell *c);
-void space_split(struct space *s, struct cell *cells, int verbose);
+void space_split(struct space *s, struct cell *cells, int nr_cells,
+ int verbose);
void space_split_mapper(void *map_data, int num_elements, void *extra_data);
void space_do_parts_sort();
void space_do_gparts_sort();