Merge branch 'error_messages' into 'master'

Sanitize smoothing lengths

Here is a proposal to address problems such as #211. 

At the start, after having read in the particles, constructed the top-level cells and all their progenitors we do a quick pass to sanitize the whole thing. If a top-level cell has more than 100000 particles (say) but has an `h_max` that prevents tasks associated with this cell from being split we do the following:

 - Compute the geometric mean and standard deviation of the smoothing lengths in that cell.
 - Limit the smoothing lengths of all particles in this cell to mean+4sigma (say).

Recall that the cells are always split so that the scheduler will be able to push down the tasks at a later point. 
This procedure is applied to the top-level cells only.

Does that make sense to you ?

See merge request !258

examples/parameter_example.yml

@@ -12,6 +12,7 @@ Scheduler:
   cell_max_size:    8000000  # (Optional) Maximal number of interactions per task if we force the split (this is the default value).
   cell_sub_size:    64000000 # (Optional) Maximal number of interactions per sub-task  (this is the default value).
   cell_split_size:  400      # (Optional) Maximal number of particles per cell (this is the default value).
+  cell_max_count:   10000    # (Optional) Maximal number of particles per cell allowed before triggering a sanitizing (this is the default value).
 
 # Parameters governing the time integration
 TimeIntegration:

src/cell.c

@@ -679,6 +679,59 @@ 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 = 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));
+
+  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;
+
+  } else {
+
+    message("Smoothing lengths will not be limited.");
+  }
+}
+
 /**
  * @brief Converts hydro quantities to a valid state after the initial density
  * calculation

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);

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);

src/scheduler.c

@@ -141,6 +141,7 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
 
       /* Get a handle on the cell involved. */
       struct cell *ci = t->ci;
+      const double hi = ci->dmin;
 
       /* Foreign task? */
       if (ci->nodeID != s->nodeID) {
@@ -149,7 +150,7 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
       }
 
       /* Is this cell even split? */
-      if (ci->split) {
+      if (ci->split && ci->h_max * kernel_gamma * space_stretch < hi / 2) {
 
         /* Make a sub? */
         if (scheduler_dosub &&

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] = {
@@ -228,7 +229,14 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
   if (s->periodic && (cdim[0] < 3 || cdim[1] < 3 || cdim[2] < 3))
     error(
         "Must have at least 3 cells in each spatial dimension when periodicity "
-        "is switched on.");
+        "is switched on.\nThis error is often caused by any of the "
+        "followings:\n"
+        " - too few particles to generate a sensible grid,\n"
+        " - the initial value of 'SPH:max_smoothing_length' is too large,\n"
+        " - the (minimal) time-step is too large leading to particles with "
+        "predicted smoothing lengths too large for the box size,\n"
+        " - particle with velocities so large that they move by more than two "
+        "box sizes per time-step.\n");
 
   /* Check if we have enough cells for gravity. */
   if (s->gravity && (cdim[0] < 8 || cdim[1] < 8 || cdim[2] < 8))
@@ -691,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),
@@ -704,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)
@@ -719,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.
@@ -1240,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];
 
@@ -1571,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);

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 10000
 #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();