Implement an optional minimal smoothing length based on the softening of the particle

src/hydro_properties.c

@@ -30,12 +30,16 @@
 #include "dimension.h"
 #include "equation_of_state.h"
 #include "error.h"
+#include "gravity_properties.h"
 #include "hydro.h"
 #include "kernel_hydro.h"
+#include "parser.h"
+#include "units.h"
 
 #define hydro_props_default_max_iterations 30
 #define hydro_props_default_volume_change 1.4f
 #define hydro_props_default_h_max FLT_MAX
+#define hydro_props_default_h_min_ratio 0.f
 #define hydro_props_default_h_tolerance 1e-4
 #define hydro_props_default_init_temp 0.f
 #define hydro_props_default_min_temp 0.f
@@ -104,6 +108,13 @@ void hydro_props_init(struct hydro_props *p,
   p->h_max = parser_get_opt_param_float(params, "SPH:h_max",
                                         hydro_props_default_h_max);
 
+  /* Minimal smoothing length ratio to softening */
+  p->h_min_ratio = parser_get_opt_param_float(params, "SPH:h_min_ratio",
+                                              hydro_props_default_h_min_ratio);
+
+  /* Temporarily set the minimal softening to 0. */
+  p->h_min = 0.f;
+
   /* Number of iterations to converge h */
   p->max_smoothing_iterations = parser_get_opt_param_int(
       params, "SPH:max_ghost_iterations", hydro_props_default_max_iterations);
@@ -317,6 +328,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
  * @param p the struct
  */
 void hydro_props_init_no_hydro(struct hydro_props *p) {
+
   p->eta_neighbours = 1.2348;
   p->h_tolerance = hydro_props_default_h_tolerance;
   p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
@@ -325,6 +337,8 @@ void hydro_props_init_no_hydro(struct hydro_props *p) {
       (pow_dimension(delta_eta) - pow_dimension(p->eta_neighbours)) *
       kernel_norm;
   p->h_max = hydro_props_default_h_max;
+  p->h_min = 0.f;
+  p->h_min_ratio = hydro_props_default_h_min_ratio;
   p->max_smoothing_iterations = hydro_props_default_max_iterations;
   p->CFL_condition = 0.1;
   p->log_max_h_change = logf(powf(1.4, hydro_dimension_inv));
@@ -352,6 +366,20 @@ void hydro_props_init_no_hydro(struct hydro_props *p) {
   p->diffusion.alpha_min = hydro_props_default_diffusion_alpha_min;
 }
 
+/**
+ * @brief Update the global properties of the hydro scheme for that time-step.
+ *
+ * @param p The properties to update.
+ * @param gp The properties of the gravity scheme.
+ * @param cosmo The cosmological model.
+ */
+void hydro_props_update(struct hydro_props *p, const struct gravity_props *gp,
+                        const struct cosmology *cosmo) {
+
+  /* Update the minimal allowed smoothing length */
+  p->h_min = p->h_min_ratio * gp->epsilon_cur;
+}
+
 /**
  * @brief Write a hydro_props struct to the given FILE as a stream of bytes.
  *

src/hydro_properties.h

@@ -32,10 +32,14 @@
 #endif
 
 /* Local includes. */
-#include "parser.h"
-#include "physical_constants.h"
 #include "restart.h"
-#include "units.h"
+
+/* Forward declarations */
+struct cosmology;
+struct swift_params;
+struct gravity_props;
+struct phys_const;
+struct unit_system;
 
 /**
  * @brief Contains all the constants and parameters of the hydro scheme
@@ -57,6 +61,12 @@ struct hydro_props {
   /*! Maximal smoothing length */
   float h_max;
 
+  /*! Minimal smoothing length expressed as ratio to softening length */
+  float h_min_ratio;
+
+  /*! Minimal smoothing length */
+  float h_min;
+
   /*! Maximal number of iterations to converge h */
   int max_smoothing_iterations;
 
@@ -131,6 +141,9 @@ void hydro_props_init(struct hydro_props *p,
                       const struct unit_system *us,
                       struct swift_params *params);
 
+void hydro_props_update(struct hydro_props *p, const struct gravity_props *gp,
+                        const struct cosmology *cosmo);
+
 #if defined(HAVE_HDF5)
 void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p);
 #endif

src/runner.c

@@ -1230,6 +1230,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
   const struct cosmology *cosmo = e->cosmology;
   const struct chemistry_global_data *chemistry = e->chemistry;
   const float hydro_h_max = e->hydro_properties->h_max;
+  const float hydro_h_min = e->hydro_properties->h_min;
   const float eps = e->hydro_properties->h_tolerance;
   const float hydro_eta_dim =
       pow_dimension(e->hydro_properties->eta_neighbours);
@@ -1294,6 +1295,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
         const float h_old = p->h;
         const float h_old_dim = pow_dimension(h_old);
         const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);
+
         float h_new;
         int has_no_neighbours = 0;
 
@@ -1330,11 +1332,13 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 #endif
 
           /* Skip if h is already h_max and we don't have enough neighbours */
-          if ((p->h >= hydro_h_max) && (f < 0.f)) {
+          if (((p->h >= hydro_h_max) && (f < 0.f)) ||
+              ((p->h <= hydro_h_min) && (f > 0.f))) {
 
           /* We have a particle whose smoothing length is already set (wants
-           * to be larger but has already hit the maximum). So, just tidy up
-           * as if the smoothing length had converged correctly  */
+           * to be larger but has already hit the maximum OR wants to be smaller
+           * but has already reached the minimum). So, just tidy up as if the
+           * smoothing length had converged correctly  */
 
 #ifdef EXTRA_HYDRO_LOOP
 
@@ -1436,8 +1440,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
             p->h = h_new;
           }
 
-          /* If below the absolute maximum, try again */
-          if (p->h < hydro_h_max) {
+          /* If within the allowed range, try again */
+          if (p->h < hydro_h_max && p->h > hydro_h_min) {
 
             /* Flag for another round of fun */
             pid[redo] = pid[i];
@@ -1453,7 +1457,18 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
             /* Off we go ! */
             continue;
 
-          } else {
+          } else if (p->h <= hydro_h_min) {
+
+            /* Ok, this particle is a lost cause... */
+            p->h = hydro_h_min;
+
+            /* Do some damage control if no neighbours at all were found */
+            if (has_no_neighbours) {
+              hydro_part_has_no_neighbours(p, xp, cosmo);
+              chemistry_part_has_no_neighbours(p, xp, chemistry, cosmo);
+            }
+
+          } else if (p->h >= hydro_h_max) {
 
             /* Ok, this particle is a lost cause... */
             p->h = hydro_h_max;

src/space.c

@@ -3351,6 +3351,10 @@ void space_first_init_parts_mapper(void *restrict map_data, int count,
 
   const struct hydro_props *hydro_props = s->e->hydro_properties;
   const float u_init = hydro_props->initial_internal_energy;
+  const float hydro_h_min_ratio = e->hydro_properties->h_min_ratio;
+
+  const struct gravity_props *grav_props = s->e->gravity_properties;
+  const int with_gravity = e->policy & engine_policy_self_gravity;
 
   const struct chemistry_global_data *chemistry = e->chemistry;
   const struct cooling_function_data *cool_func = e->cooling_func;
@@ -3360,6 +3364,12 @@ void space_first_init_parts_mapper(void *restrict map_data, int count,
     if (p[k].h <= 0.)
       error("Invalid value of smoothing length for part %lld h=%e", p[k].id,
             p[k].h);
+
+    if (with_gravity) {
+      const struct gpart *gp = p[k].gpart;
+      const float softening = gravity_get_softening(gp, grav_props);
+      p->h = max(p->h, softening * hydro_h_min_ratio);
+    }
   }
 
   /* Convert velocities to internal units */