diff --git a/examples/main.c b/examples/main.c
index ad28d16c25024f07f1f8d3433bc35c14fa24a426..37974c652a412d5ef5b152214c99ba4b96102640 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -418,6 +418,8 @@ int main(int argc, char *argv[]) {
parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
const int replicate =
parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
+ const int clean_h_values =
+ parser_get_opt_param_int(params, "InitialConditions:cleanup_h", 0);
if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
fflush(stdout);
@@ -619,7 +621,7 @@ int main(int argc, char *argv[]) {
#endif
/* Initialise the particles */
- engine_init_particles(&e, flag_entropy_ICs);
+ engine_init_particles(&e, flag_entropy_ICs, clean_h_values);
/* Write the state of the system before starting time integration. */
engine_dump_snapshot(&e);
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 353f0dbe34ebe14a057e877de289b54516a043fc..97e7948317387356d4eed5dbc3b9acf673486fba 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -60,6 +60,7 @@ Gravity:
# Parameters related to the initial conditions
InitialConditions:
file_name: SedovBlast/sedov.hdf5 # The file to read
+ cleanup_h: 0 # (Optional) Clean the values of h that are read in. Set to 1 to activate.
h_scaling: 1. # (Optional) A scaling factor to apply to all smoothing lengths in the ICs.
shift_x: 0. # (Optional) A shift to apply to all particles read from the ICs (in internal units).
shift_y: 0.
diff --git a/src/cell.c b/src/cell.c
index dbccfd2f42cabf38417cd87de0450489240884be..15ea31bef8e9219f47d4a060474deec55f54f0a1 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -941,53 +941,52 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_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.
+ * Each cell with <1000 part will be processed. We limit h to be the size of
+ * the cell and replace 0s with a good estimate.
*
* @param c The cell.
+ * @param treated Has the cell already been sanitized at this level ?
*/
-void cell_sanitize(struct cell *c) {
+void cell_sanitize(struct cell *c, int treated) {
const int count = c->count;
struct part *parts = c->parts;
+ float h_max = 0.f;
- /* 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) {
+ /* Treat cells will <1000 particles */
+ if (count < 1000 && !treated) {
- const float h = logf(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 = (h_var > 0.f) ? sqrtf(h_var) : 0.1f * h_mean;
-
- /* Choose a cut */
- const float h_limit = expf(h_mean + 4.f * h_std);
-
- /* Be verbose this is not innocuous */
- message("Cell properties: h_min= %f h_max= %f geometric mean= %f.",
- expf(h_min), expf(h_max), expf(h_mean));
+ /* Get an upper bound on h */
+ const float upper_h_max = c->dmin / (1.2f * kernel_gamma);
- if (c->h_max > h_limit) {
+ /* Apply it */
+ for (int i = 0; i < count; ++i) {
+ if (parts[i].h == 0.f || parts[i].h > upper_h_max)
+ parts[i].h = upper_h_max;
+ }
+ }
- message("Smoothing lengths will be limited to (mean + 4sigma)= %f.",
- h_limit);
+ /* Recurse and gather the new h_max values */
+ if (c->split) {
- /* Apply the cut */
- for (int i = 0; i < count; ++i) parts->h = min(parts[i].h, h_limit);
+ for (int k = 0; k < 8; ++k) {
+ if (c->progeny[k] != NULL) {
- c->h_max = h_limit;
+ /* Recurse */
+ cell_sanitize(c->progeny[k], (count < 1000));
+ /* And collect */
+ h_max = max(h_max, c->progeny[k]->h_max);
+ }
+ }
} else {
- message("Smoothing lengths will not be limited.");
+ /* Get the new value of h_max */
+ for (int i = 0; i < count; ++i) h_max = max(h_max, parts[i].h);
}
+
+ /* Record the change */
+ c->h_max = h_max;
}
/**
diff --git a/src/cell.h b/src/cell.h
index 2e32533402110040310be88629d0fb33f0128c62..178f77592543ed1a617d76d3d9fae91079395353 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -344,7 +344,7 @@ struct cell {
void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
struct cell_buff *buff, struct cell_buff *sbuff,
struct cell_buff *gbuff);
-void cell_sanitize(struct cell *c);
+void cell_sanitize(struct cell *c, int treated);
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 417c9f626d7e2f8d96d49d8d2bed942102b96e4f..66170b83ebfc814f562008ddcfe1eea8a98d8c0f 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2790,8 +2790,10 @@ void engine_print_task_counts(struct engine *e) {
* @brief Rebuild the space and tasks.
*
* @param e The #engine.
+ * @param clean_h_values Are we cleaning up the values of h before building
+ * the tasks ?
*/
-void engine_rebuild(struct engine *e) {
+void engine_rebuild(struct engine *e, int clean_h_values) {
const ticks tic = getticks();
@@ -2802,7 +2804,7 @@ void engine_rebuild(struct engine *e) {
space_rebuild(e->s, e->verbose);
/* Initial cleaning up session ? */
- if (e->s->sanitized == 0) space_sanitize(e->s);
+ if (clean_h_values) space_sanitize(e->s);
/* If in parallel, exchange the cell structure. */
#ifdef WITH_MPI
@@ -2856,7 +2858,7 @@ void engine_prepare(struct engine *e) {
if (e->forcerepart) engine_repartition(e);
/* Do we need rebuilding ? */
- if (e->forcerebuild) engine_rebuild(e);
+ if (e->forcerebuild) engine_rebuild(e, 0);
/* Unskip active tasks and check for rebuild */
engine_unskip(e);
@@ -3218,9 +3220,12 @@ void engine_launch(struct engine *e, int nr_runners) {
*
* @param e The #engine
* @param flag_entropy_ICs Did the 'Internal Energy' of the particles actually
- *contain entropy ?
+ * contain entropy ?
+ * @param clean_h_values Are we cleaning up the values of h before building
+ * the tasks ?
*/
-void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
+void engine_init_particles(struct engine *e, int flag_entropy_ICs,
+ int clean_h_values) {
struct space *s = e->s;
@@ -3237,7 +3242,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
}
/* Construct all cells and tasks to start everything */
- engine_rebuild(e);
+ engine_rebuild(e, clean_h_values);
/* No time integration. We just want the density and ghosts */
engine_skip_force_and_kick(e);
diff --git a/src/engine.h b/src/engine.h
index e62b12332d3ac1b985b8f6d7181ea66824ec4f13..d621f4f6fba3880811980a45c5a6b43e3009461f 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -272,7 +272,8 @@ void engine_init(struct engine *e, struct space *s,
struct sourceterms *sourceterms);
void engine_launch(struct engine *e, int nr_runners);
void engine_prepare(struct engine *e);
-void engine_init_particles(struct engine *e, int flag_entropy_ICs);
+void engine_init_particles(struct engine *e, int flag_entropy_ICs,
+ int clean_h_values);
void engine_step(struct engine *e);
void engine_maketasks(struct engine *e);
void engine_split(struct engine *e, struct partition *initial_partition);
@@ -281,7 +282,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
int *ind_gpart, size_t *Ngpart,
size_t offset_sparts, int *ind_spart,
size_t *Nspart);
-void engine_rebuild(struct engine *e);
+void engine_rebuild(struct engine *e, int clean_h_values);
void engine_repartition(struct engine *e);
void engine_repartition_trigger(struct engine *e);
void engine_makeproxies(struct engine *e);
diff --git a/src/space.c b/src/space.c
index b1612876b6339fb29648a87e9aec93a1d8f64664..5d7e306ebe331d449daf607a3405712b4f42f4b1 100644
--- a/src/space.c
+++ b/src/space.c
@@ -962,28 +962,33 @@ void space_split(struct space *s, struct cell *cells, int nr_cells,
}
/**
- * @brief Runs through the top-level cells and checks whether tasks associated
- * with them can be split. If not, try to sanitize the cells.
+ * @brief #threadpool mapper function to sanitize the cells
*
- * @param s The #space to act upon.
+ * @param map_data Pointers towards the top-level cells.
+ * @param num_cells The number of top-level cells.
+ * @param extra_data Unused parameters.
*/
-void space_sanitize(struct space *s) {
-
- s->sanitized = 1;
+void space_sanitize_mapper(void *map_data, int num_cells, void *extra_data) {
+ /* Unpack the inputs. */
+ struct cell *cells_top = (struct cell *)map_data;
- for (int k = 0; k < s->nr_cells; k++) {
+ for (int ind = 0; ind < num_cells; ind++) {
+ struct cell *c = &cells_top[ind];
+ cell_sanitize(c, 0);
+ }
+}
- struct cell *c = &s->cells_top[k];
- const double min_width = c->dmin;
+/**
+ * @brief Runs through the top-level cells and sanitize their h values
+ *
+ * @param s The #space to act upon.
+ */
+void space_sanitize(struct space *s) {
- /* Do we have a problem ? */
- if (c->h_max * kernel_gamma * space_stretch > min_width * 0.5 &&
- c->count > space_maxcount) {
+ if (s->e->nodeID == 0) message("Cleaning up unreasonable values of h");
- /* Ok, clean-up the mess */
- cell_sanitize(c);
- }
- }
+ threadpool_map(&s->e->threadpool, space_sanitize_mapper, s->cells_top,
+ s->nr_cells, sizeof(struct cell), 1, NULL);
}
/**
@@ -2630,7 +2635,6 @@ void space_init(struct space *s, const struct swift_params *params,
s->dim[0] = dim[0];
s->dim[1] = dim[1];
s->dim[2] = dim[2];
- s->sanitized = 0;
s->periodic = periodic;
s->gravity = gravity;
s->nr_parts = Npart;
diff --git a/src/space.h b/src/space.h
index e8e8600349c97ff8a60f0fcf2964d6ec514a7589..5c26f359e15b757b66519d958c595f670c5d6990 100644
--- a/src/space.h
+++ b/src/space.h
@@ -139,9 +139,6 @@ struct space {
/*! Number of queues in the system. */
int nr_queues;
- /*! Has this space already been sanitized ? */
- int sanitized;
-
/*! The associated engine. */
struct engine *e;