Add cosmological evolution of the gravitational softening.

examples/EAGLE_100/eagle_100.yml

@@ -37,10 +37,11 @@ Statistics:
 
 # Parameters for the self-gravity scheme
 Gravity:
-  eta:                   0.025    # Constant dimensionless multiplier for time integration. 
-  epsilon:               0.001   # Softening length (in internal units).
-  theta:                 0.85      # Opening angle (Multipole acceptance criterion)
-  
+  eta:                    0.025    # Constant dimensionless multiplier for time integration. 
+  theta:                  0.85      # Opening angle (Multipole acceptance criterion)
+  comoving_softening:     0.0026994 # Comoving softening length (in internal units).
+  max_physical_softening: 0.0007    # Physical softening length (in internal units).
+
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).

examples/EAGLE_12/eagle_12.yml

@@ -38,9 +38,10 @@ Statistics:
 
 # Parameters for the self-gravity scheme
 Gravity:
-  eta:                   0.025    # Constant dimensionless multiplier for time integration.
-  epsilon:               0.001    # Softening length (in internal units).
-  theta:                 0.85     # Opening angle (Multipole acceptance criterion)
+  eta:                    0.025     # Constant dimensionless multiplier for time integration.
+  theta:                  0.85      # Opening angle (Multipole acceptance criterion)
+  comoving_softening:     0.0026994 # Comoving softening length (in internal units).
+  max_physical_softening: 0.0007    # Physical softening length (in internal units).
   
 # Parameters for the hydrodynamics scheme
 SPH:

examples/EAGLE_25/eagle_25.yml

@@ -38,10 +38,11 @@ Statistics:
 
 # Parameters for the self-gravity scheme
 Gravity:
-  eta:                   0.025    # Constant dimensionless multiplier for time integration. 
-  epsilon:               0.001    # Softening length (in internal units).
-  theta:                 0.85     # Opening angle (Multipole acceptance criterion)
-  
+  eta:                    0.025    # Constant dimensionless multiplier for time integration. 
+  theta:                  0.85     # Opening angle (Multipole acceptance criterion)
+  comoving_softening:     0.0026994 # Comoving softening length (in internal units).
+  max_physical_softening: 0.0007    # Physical softening length (in internal units).
+
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).

examples/EAGLE_50/eagle_50.yml

@@ -37,9 +37,10 @@ Statistics:
 
 # Parameters for the self-gravity scheme
 Gravity:
-  eta:                   0.025    # Constant dimensionless multiplier for time integration.
-  epsilon:               0.001    # Softening length (in internal units).
-  theta:                 0.85     # Opening angle (Multipole acceptance criterion)
+  eta:                    0.025    # Constant dimensionless multiplier for time integration.
+  theta:                  0.85     # Opening angle (Multipole acceptance criterion)
+  comoving_softening:     0.0026994 # Comoving softening length (in internal units).
+  max_physical_softening: 0.0007    # Physical softening length (in internal units).
 
 # Parameters for the hydrodynamics scheme
 SPH:

examples/EAGLE_6/eagle_6.yml

@@ -38,9 +38,10 @@ Statistics:
 
 # Parameters for the self-gravity scheme
 Gravity:
-  eta:                   0.025    # Constant dimensionless multiplier for time integration.
-  epsilon:               0.001    # Softening length (in internal units).
-  theta:                 0.85     # Opening angle (Multipole acceptance criterion)
+  eta:                    0.025    # Constant dimensionless multiplier for time integration.
+  theta:                  0.85     # Opening angle (Multipole acceptance criterion)
+  comoving_softening:     0.0026994 # Comoving softening length (in internal units).
+  max_physical_softening: 0.0007    # Physical softening length (in internal units).
   
 # Parameters for the hydrodynamics scheme
 SPH:

examples/EvrardCollapse_3D/evrard.yml

@@ -31,11 +31,10 @@ SPH:
 
 # Parameters for the self-gravity scheme
 Gravity:
-  eta:                   0.025    # Constant dimensionless multiplier for time integration.
-  epsilon:               0.001      # Softening length (in internal units).
-  theta:                 0.9
-  a_smooth:              1.25     # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
-  r_cut:                 4.5      # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
+  eta:                    0.025    # Constant dimensionless multiplier for time integration.
+  theta:                  0.9
+  comoving_softening:     0.001    # Comoving softening length (in internal units).
+  max_physical_softening: 0.001    # Physical softening length (in internal units).
 
 # Parameters related to the initial conditions
 InitialConditions:

examples/main.c

@@ -614,7 +614,8 @@ int main(int argc, char *argv[]) {
     if (with_hydro) eos_init(&eos, params);
 
     /* Initialise the gravity properties */
-    if (with_self_gravity) gravity_props_init(&gravity_properties, params);
+    if (with_self_gravity)
+      gravity_props_init(&gravity_properties, params, &cosmo);
 
     /* Read particles and space information from (GADGET) ICs */
     char ICfileName[200] = "";

examples/parameter_example.yml

@@ -35,12 +35,13 @@ SPH:
   
 # Parameters for the self-gravity scheme
 Gravity:
-  eta:          0.025    # Constant dimensionless multiplier for time integration.
-  theta:        0.7      # Opening angle (Multipole acceptance criterion)
-  epsilon:      0.1      # Softening length (in internal units).
-  a_smooth:     1.25     # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
-  r_cut_max:    4.5      # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
-  r_cut_min:    0.1      # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
+  eta:          0.025               # Constant dimensionless multiplier for time integration.
+  theta:        0.7                 # Opening angle (Multipole acceptance criterion)
+  comoving_softening:     0.0026994 # Comoving softening length (in internal units).
+  max_physical_softening: 0.0007    # Physical softening length (in internal units).
+  a_smooth:     1.25                # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
+  r_cut_max:    4.5                 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
+  r_cut_min:    0.1                 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
 
 # Parameters for the task scheduling
 Scheduler:

src/engine.c

@@ -4396,6 +4396,7 @@ void engine_step(struct engine *e) {
   if (e->policy & engine_policy_cosmology) {
     e->time_old = e->time;
     cosmology_update(e->cosmology, e->physical_constants, e->ti_current);
+    gravity_update(e->gravity_properties, e->cosmology);
     e->time = e->cosmology->time;
     e->time_step = e->time - e->time_old;
   } else {
@@ -5193,7 +5194,7 @@ void engine_init(
     long long Ngas, long long Ndm, int policy, int verbose,
     struct repartition *reparttype, const struct unit_system *internal_units,
     const struct phys_const *physical_constants, struct cosmology *cosmo,
-    const struct hydro_props *hydro, const struct gravity_props *gravity,
+    const struct hydro_props *hydro, struct gravity_props *gravity,
     const struct external_potential *potential,
     const struct cooling_function_data *cooling_func,
     const struct chemistry_data *chemistry, struct sourceterms *sourceterms) {

src/engine.h

@@ -278,7 +278,7 @@ struct engine {
   const struct hydro_props *hydro_properties;
 
   /* Properties of the self-gravity scheme */
-  const struct gravity_props *gravity_properties;
+  struct gravity_props *gravity_properties;
 
   /* Properties of external gravitational potential */
   const struct external_potential *external_potential;
@@ -335,7 +335,7 @@ void engine_init(
     long long Ngas, long long Ndm, int policy, int verbose,
     struct repartition *reparttype, const struct unit_system *internal_units,
     const struct phys_const *physical_constants, struct cosmology *cosmo,
-    const struct hydro_props *hydro, const struct gravity_props *gravity,
+    const struct hydro_props *hydro, struct gravity_props *gravity,
     const struct external_potential *potential,
     const struct cooling_function_data *cooling_func,
     const struct chemistry_data *chemistry, struct sourceterms *sourceterms);

src/gravity/Default/gravity.h

@@ -45,7 +45,7 @@ __attribute__((always_inline)) INLINE static float gravity_get_mass(
 __attribute__((always_inline)) INLINE static float gravity_get_softening(
     const struct gpart* gp, const struct gravity_props* restrict grav_props) {
 
-  return grav_props->epsilon;
+  return grav_props->epsilon_cur;
 }
 
 /**

src/gravity_properties.c

@@ -37,7 +37,8 @@
 #define gravity_props_default_r_cut_min 0.1f
 
 void gravity_props_init(struct gravity_props *p,
-                        const struct swift_params *params) {
+                        const struct swift_params *params,
+                        const struct cosmology *cosmo) {
 
   /* Tree-PM parameters */
   p->a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth",
@@ -56,11 +57,35 @@ void gravity_props_init(struct gravity_props *p,
   p->theta_crit2 = p->theta_crit * p->theta_crit;
   p->theta_crit_inv = 1. / p->theta_crit;
 
-  /* Softening lengths */
-  p->epsilon = 3. * parser_get_param_double(params, "Gravity:epsilon");
-  p->epsilon2 = p->epsilon * p->epsilon;
-  p->epsilon_inv = 1.f / p->epsilon;
-  p->epsilon_inv3 = p->epsilon_inv * p->epsilon_inv * p->epsilon_inv;
+  /* Softening parameters */
+  p->epsilon_comoving =
+      parser_get_param_double(params, "Gravity:comoving_softening");
+  p->epsilon_max_physical =
+      parser_get_param_double(params, "Gravity:max_physical_softening");
+
+  /* Set the softening to the current time */
+  gravity_update(p, cosmo);
+}
+
+void gravity_update(struct gravity_props *p, const struct cosmology *cosmo) {
+
+  /* Current softening lengths */
+  double softening;
+  if (p->epsilon_comoving * cosmo->a > p->epsilon_max_physical)
+    softening = p->epsilon_max_physical / cosmo->a;
+  else
+    softening = p->epsilon_comoving;
+
+  /* Plummer equivalent -> internal */
+  softening *= kernel_gravity_softening_plummer_equivalent;
+
+  /* Store things */
+  p->epsilon_cur = softening;
+
+  /* Other factors */
+  p->epsilon_cur2 = softening * softening;
+  p->epsilon_cur_inv = 1. / softening;
+  p->epsilon_cur_inv3 = 1. / (softening * softening * softening);
 }
 
 void gravity_props_print(const struct gravity_props *p) {
@@ -72,8 +97,17 @@ void gravity_props_print(const struct gravity_props *p) {
 
   message("Self-gravity opening angle:  theta=%.4f", p->theta_crit);
 
-  message("Self-gravity softening:    epsilon=%.4f (Plummer equivalent: %.4f)",
-          p->epsilon, p->epsilon / 3.);
+  message(
+      "Self-gravity comoving softening:    epsilon=%.4f (Plummer equivalent: "
+      "%.4f)",
+      p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent,
+      p->epsilon_comoving);
+
+  message(
+      "Self-gravity maximal physical softening:    epsilon=%.4f (Plummer "
+      "equivalent: %.4f)",
+      p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent,
+      p->epsilon_max_physical);
 
   message("Self-gravity mesh smoothing-scale: a_smooth=%f", p->a_smooth);
 
@@ -86,9 +120,18 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
                                   const struct gravity_props *p) {
 
   io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta);
-  io_write_attribute_f(h_grpgrav, "Softening length", p->epsilon);
-  io_write_attribute_f(h_grpgrav, "Softening length (Plummer equivalent)",
-                       p->epsilon / 3.);
+  io_write_attribute_f(
+      h_grpgrav, "Comoving softening length",
+      p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent);
+  io_write_attribute_f(h_grpgrav,
+                       "Comoving Softening length (Plummer equivalent)",
+                       p->epsilon_comoving);
+  io_write_attribute_f(
+      h_grpgrav, "Maximal physical softening length",
+      p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent);
+  io_write_attribute_f(h_grpgrav,
+                       "Maximal physical softening length (Plummer equivalent)",
+                       p->epsilon_max_physical);
   io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
   io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
   io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);

src/gravity_properties.h

@@ -27,6 +27,7 @@
 #endif
 
 /* Local includes. */
+#include "cosmology.h"
 #include "parser.h"
 #include "restart.h"
 
@@ -58,22 +59,30 @@ struct gravity_props {
   /*! Inverse of opening angle */
   double theta_crit_inv;
 
-  /*! Softening length */
-  float epsilon;
+  /*! Comoving softening */
+  double epsilon_comoving;
 
-  /*! Square of softening length */
-  float epsilon2;
+  /*! Maxium physical softening */
+  double epsilon_max_physical;
 
-  /*! Inverse of softening length */
-  float epsilon_inv;
+  /*! Current sftening length */
+  float epsilon_cur;
 
-  /*! Cube of the inverse of softening length */
-  float epsilon_inv3;
+  /*! Square of current softening length */
+  float epsilon_cur2;
+
+  /*! Inverse of current softening length */
+  float epsilon_cur_inv;
+
+  /*! Cube of the inverse of current oftening length */
+  float epsilon_cur_inv3;
 };
 
 void gravity_props_print(const struct gravity_props *p);
 void gravity_props_init(struct gravity_props *p,
-                        const struct swift_params *params);
+                        const struct swift_params *params,
+                        const struct cosmology *cosmo);
+void gravity_update(struct gravity_props *p, const struct cosmology *cosmo);
 
 #if defined(HAVE_HDF5)
 void gravity_props_print_snapshot(hid_t h_grpsph,

src/kernel_gravity.h

@@ -26,6 +26,9 @@
 #include "inline.h"
 #include "minmax.h"
 
+/*! Conversion factor between Plummer softening and internal softening */
+#define kernel_gravity_softening_plummer_equivalent 3.
+
 /**
  * @brief Computes the gravity softening function for potential.
  *

src/multipole.h

@@ -1524,8 +1524,8 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
                                const double dim[3]) {
 
   /* Recover some constants */
-  const float eps = props->epsilon;
-  const float eps_inv = props->epsilon_inv;
+  const float eps = props->epsilon_cur;
+  const float eps_inv = props->epsilon_cur_inv;
 
   /* Compute distance vector */
   float dx = (float)(pos_b[0] - pos_a[0]);
@@ -2380,6 +2380,7 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
   gp->a_grav[2] += a_grav[2];
   gp->potential += pot;
 }
+
 /**
  * @brief Checks whether a cell-cell interaction can be appromixated by a M-M
  * interaction using the distance and cell radius.