Sorry I messed up the git so it's one big commit.

Updated RTD as per your recomendations.

Added ForceChecks to paramater_example.yml

Changed brute_force_gravity_flag to force_checks_snapshot_flag

Removed epsilon from gravity_exact_force_file_exits

Made functions in src/kernel_long_gravity.h double

doc/RTD/source/GettingStarted/special_modes.rst

@@ -53,8 +53,8 @@ simply the Newtonian sum, i.e.,
 
 To run SWIFT in this mode, configure the code with
 ``--enable-gravity-force-checks=N``, which means that the exact forces will be
-computed for every :math:`N^{th}` particle (i.e., to compute the exact forces
-for all particles set ``N=1``).
+computed for every :math:`N^{th}` particle based on their ID (i.e., to compute
+the exact forces for all particles set ``N=1``).
 
 Two `.dat` files will be output during each timestep, one containing the forces
 (really it is the accelerations that are stored) as computed by ``_swift_``, and
@@ -66,5 +66,8 @@ the number of particles contributing to each force component is also stored
 (these numbers will add up to :math:`n-1`).   
 
 This mode will slow down the code *considerably*, and it is not recommended to
-run on problems with more than :math:`10^{5}` particles. This mode must be run
-on a single node/rank and is only designed for pure gravity tests (i.e., DMO).
+be run on problems with more than :math:`10^{5}` particles when
+``--enable-gravity-force-checks=1``. For larger runs, sampling a sub-set of
+particles via the argument ``N`` of the configuration option is recommended.
+This mode must be run on a single node/rank, and is primarily designed for pure
+gravity tests (i.e., DMO).

doc/RTD/source/ParameterFiles/parameter_description.rst

@@ -1258,19 +1258,20 @@ Showing all the parameters for a basic cosmologica test-case, one would have:
 Gravity Force Checks
 --------------------
 
-By default, when configured with ``--enable-gravity-force-checks`` the "exact"
-forces will be computed for all active particles during each timestep.
+By default, when the code is configured with ``--enable-gravity-force-checks``,
+the "exact" forces of all active gparts are computed during each timestep.
 
-To give a bit more control over this, you can select to only compute the exact
-forces during the timesteps that all the particles are active, and/or only at
-the timesteps when a snapshot is being dumped, i.e.,
+To give a bit more control over this, you can select to only perform the exact
+force computation during the timesteps that all gparts are active, and/or only
+at the timesteps when a snapshot is being dumped, i.e.,
 
 .. code:: YAML
 
   ForceChecks:
-    only_when_all_active:   1    # Only compute exact forces during timesteps when all particles are active.
+    only_when_all_active:   1    # Only compute exact forces during timesteps when all gparts are active.
     only_at_snapshots:      1    # Only compute exact forces during timesteps when a snapshot is being dumped.
 
-Note that if ``only_when_all_active:1`` and ``only_when_all_active:1``, the
-exact forces will be computed during the first timestep after the snapshot
-output time when all particles are active.
+If ``only_when_all_active:1`` and ``only_when_all_active:1`` are enabled
+together, and all the gparts are not active during the timestep of the snapshot
+dump, the exact forces computation is performed on the first timestep at which
+all the gparts are active after that snapshot output timestep.

examples/parameter_example.yml

@@ -77,7 +77,11 @@ Gravity:
   dithering:                     1         # (Optional) Activate the dithering of the gravity mesh at every rebuild (this is the default value).
   dithering_ratio:               1.0       # (Optional) Magnitude of each component of the dithering vector in units of the top-level cell sizes (this is the default value).
 
-# Parameters for the Friends-Of-Friends algorithm
+# Parameters when running with SWIFT_GRAVITY_FORCE_CHECKS 
+ForceChecks:
+  only_when_all_active: 1  # (Optional) Only compute exact forces during timesteps when all gparts are active (default: 0).
+  only_at_snapshots:    1  # (Optional) Only compute exact forces during timesteps when a snapshot is being dumped (default: 0).
+
 FOF:
   basename:                        fof_output  # Filename for the FOF outputs (Unused when FoF is only run to seed BHs).
   scale_factor_first:              0.91        # Scale-factor of first FoF black hole seeding calls (needed for cosmological runs).

src/engine.c

@@ -2165,7 +2165,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
   e->wallclock_time = (float)clocks_diff(&time1, &time2);
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
-  e->brute_force_gravity_flag = 0;
+  e->force_checks_snapshot_flag = 0;
 #endif
 
   if (e->verbose) message("took %.3f %s.", e->wallclock_time, clocks_getunit());
@@ -2362,12 +2362,13 @@ void engine_step(struct engine *e) {
     if ((e->all_gparts_active && e->force_checks_only_all_active) ||
         !e->force_checks_only_all_active) {
       /* Is this a snapshot timestep (and the option is selected)? */
-      if ((e->brute_force_gravity_flag == 1 &&
-           e->force_checks_only_at_snapshots) ||
-          !e->force_checks_only_at_snapshots) {
-
-        /* Do checks */
-        gravity_exact_force_compute(e->s, e);
+      if ((e->force_checks_snapshot_flag &&
+        e->force_checks_only_at_snapshots) ||
+        !e->force_checks_only_at_snapshots) {
+        
+          /* Do checks */
+          gravity_exact_force_compute(e->s, e);
+        }
       }
     }
   }
@@ -2395,15 +2396,16 @@ void engine_step(struct engine *e) {
     if ((e->all_gparts_active && e->force_checks_only_all_active) ||
         !e->force_checks_only_all_active) {
       /* Is this a snapshot timestep (and the option is selected)? */
-      if ((e->brute_force_gravity_flag == 1 &&
-           e->force_checks_only_at_snapshots) ||
-          !e->force_checks_only_at_snapshots) {
-
-        /* Do checks */
-        gravity_exact_force_check(e->s, e, 1e-1);
-
-        /* Reset flag waiting for next output time */
-        e->brute_force_gravity_flag = 0;
+      if ((e->force_checks_snapshot_flag &&
+        e->force_checks_only_at_snapshots) ||
+        !e->force_checks_only_at_snapshots) {
+        
+          /* Do checks */
+          gravity_exact_force_check(e->s, e, 1e-1);
+
+          /* Reset flag waiting for next output time */
+          e->force_checks_snapshot_flag = 0;
+        }
       }
     }
   }
@@ -2535,7 +2537,7 @@ void engine_check_for_dumps(struct engine *e) {
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
         /* Indicate we are allowed to do a brute force calculation now */
-        e->brute_force_gravity_flag = 1;
+        e->force_checks_snapshot_flag = 1;
 #endif
 
         /* Do we want a corresponding VELOCIraptor output? */

src/engine.h

@@ -495,8 +495,8 @@ struct engine {
   /* Are all gparts active this timestep? */
   int all_gparts_active;
 
-  /* Are we doing brute force calculations this timestep? */
-  int brute_force_gravity_flag;
+  /* Flag to tell brute force checks a snapshot was recently written. */
+  int force_checks_snapshot_flag;
 #endif
 };
 

src/gravity.c

@@ -430,7 +430,7 @@ int gravity_exact_force_file_exits(const struct engine *e) {
 
     char line[100];
     char dummy1[10], dummy2[10];
-    double epsilon, newton_G;
+    double newton_G;
     int N, periodic;
     /* Reads file header */
     if (fgets(line, 100, file) != line) error("Problem reading title");
@@ -438,17 +438,12 @@ int gravity_exact_force_file_exits(const struct engine *e) {
     sscanf(line, "%s %s %le", dummy1, dummy2, &newton_G);
     if (fgets(line, 100, file) != line) error("Problem reading N");
     sscanf(line, "%s %s %d", dummy1, dummy2, &N);
-    if (fgets(line, 100, file) != line) error("Problem reading epsilon");
-    sscanf(line, "%s %s %le", dummy1, dummy2, &epsilon);
     if (fgets(line, 100, file) != line) error("Problem reading BC");
     sscanf(line, "%s %s %d", dummy1, dummy2, &periodic);
     fclose(file);
 
     /* Check whether it matches the current parameters */
     if (N == SWIFT_GRAVITY_FORCE_CHECKS && periodic == e->s->periodic &&
-        (fabs(epsilon - gravity_get_softening(0, e->gravity_properties)) /
-             epsilon <
-         1e-5) &&
         (fabs(newton_G - e->physical_constants->const_newton_G) / newton_G <
          1e-5)) {
       return 1;

src/kernel_long_gravity.h

@@ -222,25 +222,25 @@ kernel_long_grav_force_eval_double(const double u, double *const W) {
 
   const double one_over_sqrt_pi = ((double)(M_2_SQRTPI * 0.5));
 
-  const double arg1 = u * 0.5f;
-  const double arg2 = -arg1 * arg1;
+    const double arg1 = u * 0.5;
+    const double arg2 = -arg1 * arg1;
 
-  const double term1 = approx_erfcf(arg1);
-  const double term2 = u * one_over_sqrt_pi * expf(arg2);
+    const double term1 = erfc(arg1);
+    const double term2 = u * one_over_sqrt_pi * exp(arg2);
 
   *W = term1 + term2;
 #else
 
-  const double x = 2.f * u;
-  const double exp_x = expf(x);  // good_approx_expf(x);
-  const double alpha = 1.f / (1.f + exp_x);
+    const double x = 2. * u;
+    const double exp_x = exp(x); 
+    const double alpha = 1. / (1. + exp_x);
 
-  /* We want 2*(x*alpha - x*alpha^2 - exp(x)*alpha + 1) */
-  *W = 1.f - alpha;
-  *W = *W * x - exp_x;
-  *W = *W * alpha + 1.f;
-  *W *= 2.f;
-#endif
+    /* We want 2*(x*alpha - x*alpha^2 - exp(x)*alpha + 1) */
+    *W = 1. - alpha;
+    *W = *W * x - exp_x;
+    *W = *W * alpha + 1.;
+    *W *= 2.;
+  #endif
 #endif
 }