Added a parameter to the gravity properties to contain the distance below...

Added a parameter to the gravity properties to contain the distance below which all computed forces are Newtonian

examples/EAGLE_100/eagle_100.yml

@@ -28,8 +28,6 @@ Gravity:
   eta:                   0.025    # Constant dimensionless multiplier for time integration. 
   epsilon:               0.0001   # Softening length (in internal units).
   theta:                 0.7      # Opening angle (Multipole acceptance criterion)
-  a_smooth:              1.25     # PM smoothing scale in top-level cell-size units
-  r_cut:                 4.5      # Distance beyong wwhich FMM forces are zero in top-level cell-size units
   
 # Parameters for the hydrodynamics scheme
 SPH:

examples/EAGLE_12/eagle_12.yml

@@ -28,8 +28,6 @@ Gravity:
   eta:                   0.025    # Constant dimensionless multiplier for time integration.
   epsilon:               0.0001   # Softening length (in internal units).
   theta:                 0.7      # Opening angle (Multipole acceptance criterion)
-  a_smooth:              1.25     # PM smoothing scale in top-level cell-size units
-  r_cut:                 4.5      # Distance beyong wwhich FMM forces are zero in top-level cell-size un
   
 # Parameters for the hydrodynamics scheme
 SPH:

examples/EAGLE_25/eagle_25.yml

@@ -28,8 +28,6 @@ Gravity:
   eta:                   0.025    # Constant dimensionless multiplier for time integration. 
   epsilon:               0.0001   # Softening length (in internal units).
   theta:                 0.7      # Opening angle (Multipole acceptance criterion)
-  a_smooth:              1.25     # PM smoothing scale in top-level cell-size units
-  r_cut:                 4.5      # Distance beyong wwhich FMM forces are zero in top-level cell-size units
   
 # Parameters for the hydrodynamics scheme
 SPH:

examples/EAGLE_50/eagle_50.yml

@@ -23,6 +23,12 @@ Snapshots:
 Statistics:
   delta_time:          1e-2 # Time between statistics output
 
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                   0.025    # Constant dimensionless multiplier for time integration.
+  epsilon:               0.0001   # Softening length (in internal units).
+  theta:                 0.7      # Opening angle (Multipole acceptance criterion)
+
 # 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/UniformDMBox/uniformBox.yml

@@ -28,8 +28,7 @@ Gravity:
   eta:                   0.025    # Constant dimensionless multiplier for time integration. 
   theta:                 0.7      # Opening angle (Multipole acceptance criterion)
   epsilon:               0.00001  # Softening length (in internal units).
-  a_smooth:              1.25     # PM smoothing scale in top-level cell-size units
-  r_cut:                 4.5      # Distance beyong wwhich FMM forces are zero in top-level cell-size un  
+ 
 # Parameters governing the conserved quantities statistics
 Statistics:
   delta_time:          1e-2 # Time between statistics output

examples/parameter_example.yml

@@ -55,7 +55,8 @@ Gravity:
   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:                 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_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 related to the initial conditions
 InitialConditions:

src/gravity_properties.c

@@ -33,7 +33,8 @@
 #include "kernel_gravity.h"
 
 #define gravity_props_default_a_smooth 1.25f
-#define gravity_props_default_r_cut 4.5f
+#define gravity_props_default_r_cut_max 4.5f
+#define gravity_props_default_r_cut_min 0.1f
 
 void gravity_props_init(struct gravity_props *p,
                         const struct swift_params *params) {
@@ -41,8 +42,10 @@ void gravity_props_init(struct gravity_props *p,
   /* Tree-PM parameters */
   p->a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth",
                                            gravity_props_default_a_smooth);
-  p->r_cut = parser_get_opt_param_float(params, "Gravity:r_cut",
-                                        gravity_props_default_r_cut);
+  p->r_cut_max = parser_get_opt_param_float(params, "Gravity:r_cut_max",
+                                            gravity_props_default_r_cut_max);
+  p->r_cut_min = parser_get_opt_param_float(params, "Gravity:r_cut_min",
+                                            gravity_props_default_r_cut_min);
 
   /* Time integration */
   p->eta = parser_get_param_float(params, "Gravity:eta");
@@ -69,9 +72,10 @@ void gravity_props_print(const struct gravity_props *p) {
   message("Self-gravity softening:    epsilon=%.4f (Plummer equivalent: %.4f)",
           p->epsilon, p->epsilon / 3.);
 
-  message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth);
+  message("Self-gravity mesh smoothing-scale: a_smooth=%f", p->a_smooth);
 
-  message("Self-gravity MM cut-off: r_cut=%f", p->r_cut);
+  message("Self-gravity tree cut-off: r_cut_max=%f", p->r_cut_max);
+  message("Self-gravity truncation cut-off: r_cut_min=%f", p->r_cut_min);
 }
 
 #if defined(HAVE_HDF5)
@@ -84,7 +88,8 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
                        p->epsilon / 3.);
   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, "MM a_smooth", p->a_smooth);
-  io_write_attribute_f(h_grpgrav, "MM r_cut", p->r_cut);
+  io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
+  io_write_attribute_f(h_grpgrav, "Mesh r_cut_max", p->r_cut_max);
+  io_write_attribute_f(h_grpgrav, "Mesh r_cut_min", p->r_cut_min);
 }
 #endif

src/gravity_properties.h

@@ -34,9 +34,16 @@
  */
 struct gravity_props {
 
-  /* Tree-PM parameters */
+  /*! Mesh smoothing scale in units of top-level cell size */
   float a_smooth;
-  float r_cut;
+
+  /*! Distance below which the truncated mesh force is Newtonian in units of
+   * a_smooth */
+  float r_cut_min;
+
+  /*! Distance above which the truncated mesh force is negligible in units of
+   * a_smooth */
+  float r_cut_max;
 
   /*! Time integration dimensionless multiplier */
   float eta;

src/runner_doiact_grav.h

@@ -664,7 +664,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
   const double cell_width = s->width[0];
   const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
   const double theta_crit_inv = props->theta_crit_inv;
-  const double max_distance = props->a_smooth * props->r_cut * cell_width;
+  const double max_distance = props->a_smooth * props->r_cut_max * cell_width;
   const double max_distance2 = max_distance * max_distance;
   struct gravity_tensors *const mi = ci->multipole;
   const double CoM[3] = {mi->CoM[0], mi->CoM[1], mi->CoM[2]};
@@ -701,9 +701,11 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
       const double dz = nearest(CoM[2] - mj->CoM[2], dim[2]);
       const double r2 = dx * dx + dy * dy + dz * dz;
 
+      /* Are we beyond the distance where the truncated forces are 0 ?*/
       if (r2 > max_distance2) {
 
 #ifdef SWIFT_DEBUG_CHECKS
+        /* Need to account for the interactions we missed */
         mi->pot.num_interacted += mj->m_pole.num_gpart;
 #endif
         continue;

src/tools.c

@@ -755,7 +755,7 @@ void gravity_n2(struct gpart *gparts, const int gcount,
                 const struct gravity_props *gravity_properties, float rlr) {
 
   const float rlr_inv = 1. / rlr;
-  const float r_cut = gravity_properties->r_cut;
+  const float r_cut = gravity_properties->r_cut_max;
   const float max_d = r_cut * rlr;
   const float max_d2 = max_d * max_d;