Merge branch 'gravity_bf' into 'master'

Gravity brute force checks, extra options

See merge request !1014

doc/RTD/source/GettingStarted/special_modes.rst

@@ -41,3 +41,33 @@ this mode, configure the code with ``--enable-glass-making``.
 Note that this will *not* generate the initial random distribution of
 particles. An initial condition file with particles still has to be provided.
 
+Gravity force checks
+~~~~~~~~~~~~~~~~~~~~
+
+To test the accuracy of the gravitational forces approximated by the code,
+SWIFT can be configured with the option to additionally compute the "exact"
+forces for each active particle during each timestep. Here the exact forces are
+simply the Newtonian sum, i.e.,
+
+:math:`\vec{F}_{i,j} = \sum^{n}_{i \neq j} \frac{G m_i m_j}{\vec{r}_{i,j}^2}.`
+
+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 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
+another containing the ``_exact_`` forces. The total force (``a_swift``), along
+with the contributed force from each component of the tree (P2P, M2P and M2L)
+and the FFT mesh if periodic (PM) is stored (i.e., ``a_swift[0]`` = ``a_p2p[0]`` +
+``a_m2p[0]`` + ``a_m2l[0]`` + ``a_PM[0]``, for the :math:`x` component). In addition,
+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
+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

@@ -1263,3 +1263,23 @@ Showing all the parameters for a basic cosmologica test-case, one would have:
 	 snapshots. Snapshots at a given time would always have the same set of
 	 digits irrespective of the number of snapshots produced before.
 
+Gravity Force Checks
+--------------------
+
+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 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 gparts are active.
+    only_at_snapshots:      1    # Only compute exact forces during timesteps when a snapshot is being dumped.
+
+If ``only_when_all_active:1`` and ``only_at_snapshots: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

@@ -2397,6 +2397,10 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
   e->forcerebuild = 1;
   e->wallclock_time = (float)clocks_diff(&time1, &time2);
 
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  e->force_checks_snapshot_flag = 0;
+#endif
+
   if (e->verbose) message("took %.3f %s.", e->wallclock_time, clocks_getunit());
 }
 
@@ -2570,9 +2574,38 @@ void engine_step(struct engine *e) {
 #endif
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
-  /* Run the brute-force gravity calculation for some gparts */
-  if (e->policy & engine_policy_self_gravity)
-    gravity_exact_force_compute(e->s, e);
+  /* Do we need to check if all gparts are active? */
+  if (e->force_checks_only_all_active) {
+    size_t nr_gparts = e->s->nr_gparts;
+    e->all_gparts_active = 1;
+
+    /* Look for inactive gparts */
+    for (size_t i = 0; i < nr_gparts; ++i) {
+      struct gpart *gp = &e->s->gparts[i];
+
+      /* If one gpart is inactive we can stop. */
+      if (!gpart_is_active(gp, e)) {
+        e->all_gparts_active = 0;
+        break;
+      }
+    }
+  }
+
+  /* Check if we want to run force checks this timestep. */
+  if (e->policy & engine_policy_self_gravity) {
+    /* Are all gparts active (and the option is selected)? */
+    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->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);
+      }
+    }
+  }
 #endif
 
     /* Get current CPU times.*/
@@ -2609,9 +2642,24 @@ void engine_step(struct engine *e) {
   }
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
-  /* Check the accuracy of the gravity calculation */
-  if (e->policy & engine_policy_self_gravity)
-    gravity_exact_force_check(e->s, e, 1e-1);
+  /* Check if we want to run force checks this timestep. */
+  if (e->policy & engine_policy_self_gravity) {
+    /* Are all gparts active (and the option is selected)? */
+    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->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;
+      }
+    }
+  }
 #endif
 
 #ifdef SWIFT_DEBUG_CHECKS
@@ -2740,6 +2788,11 @@ void engine_check_for_dumps(struct engine *e) {
 
       case output_snapshot:
 
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+        /* Indicate we are allowed to do a brute force calculation now */
+        e->force_checks_snapshot_flag = 1;
+#endif
+
         /* Do we want a corresponding VELOCIraptor output? */
         if (with_stf && e->snapshot_invoke_stf && !e->stf_this_timestep) {
 
@@ -3841,6 +3894,13 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
   }
 #endif
 
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  e->force_checks_only_all_active =
+      parser_get_opt_param_int(params, "ForceChecks:only_when_all_active", 0);
+  e->force_checks_only_at_snapshots =
+      parser_get_opt_param_int(params, "ForceChecks:only_at_snapshots", 0);
+#endif
+
   /* Make the space link back to the engine. */
   s->e = e;
 

src/engine.h

@@ -486,6 +486,20 @@ struct engine {
 
   /* Has there been an stf this timestep? */
   char stf_this_timestep;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  /* Run brute force checks only on steps when all gparts active? */
+  int force_checks_only_all_active;
+
+  /* Run brute force checks only during snapshot timesteps? */
+  int force_checks_only_at_snapshots;
+
+  /* Are all gparts active this timestep? */
+  int all_gparts_active;
+
+  /* Flag to tell brute force checks a snapshot was recently written. */
+  int force_checks_snapshot_flag;
+#endif
 };
 
 /* Function prototypes, engine.c. */

src/gravity/Default/gravity.h

@@ -152,9 +152,20 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   gp->potential_PM = 0.f;
-  gp->a_grav_PM[0] = 0.f;
-  gp->a_grav_PM[1] = 0.f;
-  gp->a_grav_PM[2] = 0.f;
+
+  /* Track accelerations of each component. */
+  for (int i = 0; i < 3; i++) {
+    gp->a_grav_PM[i] = 0.f;
+    gp->a_grav_p2p[i] = 0.f;
+    gp->a_grav_m2p[i] = 0.f;
+    gp->a_grav_m2l[i] = 0.f;
+  }
+
+  /* Interaction counters. */
+  gp->num_interacted_m2p = 0;
+  gp->num_interacted_m2l = 0;
+  gp->num_interacted_p2p = 0;
+  gp->num_interacted_pm = 0;
 #endif
 
 #ifdef SWIFT_DEBUG_CHECKS
@@ -188,9 +199,12 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   gp->potential_PM *= const_G;
-  gp->a_grav_PM[0] *= const_G;
-  gp->a_grav_PM[1] *= const_G;
-  gp->a_grav_PM[2] *= const_G;
+  for (int i = 0; i < 3; i++) {
+    gp->a_grav_PM[i] *= const_G;
+    gp->a_grav_p2p[i] *= const_G;
+    gp->a_grav_m2p[i] *= const_G;
+    gp->a_grav_m2l[i] *= const_G;
+  }
 #endif
 
 #ifdef SWIFT_DEBUG_CHECKS

src/gravity/Default/gravity_part.h

@@ -58,9 +58,6 @@ struct gpart {
 
 #ifdef SWIFT_DEBUG_CHECKS
 
-  /* Numer of gparts this gpart interacted with */
-  long long num_interacted;
-
   /* Time of the last drift */
   integertime_t ti_drift;
 
@@ -70,6 +67,8 @@ struct gpart {
   /* Has this particle been initialised? */
   int initialised;
 
+  /* Total number of gparts this gpart interacted with */
+  long long num_interacted;
 #endif
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
@@ -80,11 +79,24 @@ struct gpart {
   /*! Potential taken from the mesh only */
   float potential_PM;
 
-  /* Brute-force particle acceleration. */
+  /* Acceleration taken from each component of the tree */
+  float a_grav_p2p[3];
+  float a_grav_m2p[3];
+  float a_grav_m2l[3];
+
+  /* Brute-force particle accelerations */
   double a_grav_exact[3];
+  double a_grav_exact_long[3];
+  double a_grav_exact_short[3];
 
   /* Brute-force particle potential. */
   double potential_exact;
+
+  /* Type specific interaction counters */
+  long long num_interacted_m2p;
+  long long num_interacted_m2l;
+  long long num_interacted_p2p;
+  long long num_interacted_pm;
 #endif
 
 } SWIFT_STRUCT_ALIGN;

src/gravity/MultiSoftening/gravity.h

@@ -162,9 +162,20 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   gp->potential_PM = 0.f;
-  gp->a_grav_PM[0] = 0.f;
-  gp->a_grav_PM[1] = 0.f;
-  gp->a_grav_PM[2] = 0.f;
+
+  /* Track accelerations of each component. */
+  for (int i = 0; i < 3; i++) {
+    gp->a_grav_PM[i] = 0.f;
+    gp->a_grav_p2p[i] = 0.f;
+    gp->a_grav_m2p[i] = 0.f;
+    gp->a_grav_m2l[i] = 0.f;
+  }
+
+  /* Interaction counters. */
+  gp->num_interacted_m2p = 0;
+  gp->num_interacted_m2l = 0;
+  gp->num_interacted_p2p = 0;
+  gp->num_interacted_pm = 0;
 #endif
 
 #ifdef SWIFT_DEBUG_CHECKS
@@ -202,9 +213,12 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   gp->potential_PM *= const_G;
-  gp->a_grav_PM[0] *= const_G;
-  gp->a_grav_PM[1] *= const_G;
-  gp->a_grav_PM[2] *= const_G;
+  for (int i = 0; i < 3; i++) {
+    gp->a_grav_PM[i] *= const_G;
+    gp->a_grav_p2p[i] *= const_G;
+    gp->a_grav_m2p[i] *= const_G;
+    gp->a_grav_m2l[i] *= const_G;
+  }
 #endif
 
 #ifdef SWIFT_DEBUG_CHECKS

src/gravity/MultiSoftening/gravity_part.h

@@ -59,9 +59,6 @@ struct gpart {
 
 #ifdef SWIFT_DEBUG_CHECKS
 
-  /* Numer of gparts this gpart interacted with */
-  long long num_interacted;
-
   /* Time of the last drift */
   integertime_t ti_drift;
 
@@ -71,6 +68,8 @@ struct gpart {
   /* Has this particle been initialised? */
   int initialised;
 
+  /* Total number of gparts this gpart interacted with */
+  long long num_interacted;
 #endif
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
@@ -81,11 +80,24 @@ struct gpart {
   /*! Potential taken from the mesh only */
   float potential_PM;
 
-  /* Brute-force particle acceleration. */
+  /* Acceleration taken from each component of the tree */
+  float a_grav_p2p[3];
+  float a_grav_m2p[3];
+  float a_grav_m2l[3];
+
+  /* Brute-force particle accelerations */
   double a_grav_exact[3];
+  double a_grav_exact_long[3];
+  double a_grav_exact_short[3];
 
   /* Brute-force particle potential. */
   double potential_exact;
+
+  /* Type specific interaction counters */
+  long long num_interacted_m2l;
+  long long num_interacted_m2p;
+  long long num_interacted_p2p;
+  long long num_interacted_pm;
 #endif
 
 } SWIFT_STRUCT_ALIGN;

src/kernel_long_gravity.h

@@ -208,6 +208,42 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval(
 #endif
 }
 
+/**
+ * @brief Computes the long-range correction term for the force calculation
+ * coming from FFT in double precision.
+ *
+ * @param u The ratio of the distance to the FFT cell scale \f$u = r/r_s\f$.
+ * @param W (return) The value of the kernel function.
+ */
+__attribute__((always_inline)) INLINE static void
+kernel_long_grav_force_eval_double(const double u, double *const W) {
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+#ifdef GADGET2_LONG_RANGE_CORRECTION
+
+  const double one_over_sqrt_pi = M_2_SQRTPI * 0.5;
+
+  const double arg1 = u * 0.5;
+  const double arg2 = -arg1 * arg1;
+
+  const double term1 = erfc(arg1);
+  const double term2 = u * one_over_sqrt_pi * exp(arg2);
+
+  *W = term1 + term2;
+#else
+
+  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. - alpha;
+  *W = *W * x - exp_x;
+  *W = *W * alpha + 1.;
+  *W *= 2.;
+#endif
+#endif
+}
+
 /**
  * @brief Returns the long-range truncation of the Poisson potential in Fourier
  * space.

src/multipole.h

@@ -92,8 +92,16 @@ struct grav_tensor {
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 5
 #error "Missing implementation for order >5"
 #endif
-#ifdef SWIFT_DEBUG_CHECKS
 
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  /* Number of gparts interacted through the tree. */
+  long long num_interacted_tree;
+
+  /* Number of gparts interacted through the FFT mesh */
+  long long num_interacted_pm;
+#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
   /* Total number of gpart this field tensor interacted with */
   long long num_interacted;
 
@@ -168,7 +176,7 @@ struct multipole {
 #error "Missing implementation for order >5"
 #endif
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
 
   /* Total number of gpart in this multipole */
   long long num_gpart;
@@ -308,8 +316,13 @@ INLINE static void gravity_field_tensors_add(
     struct grav_tensor *restrict la, const struct grav_tensor *restrict lb) {
 #ifdef SWIFT_DEBUG_CHECKS
   if (lb->num_interacted == 0) error("Adding tensors that did not interact");
+
   la->num_interacted += lb->num_interacted;
 #endif
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  la->num_interacted_tree += lb->num_interacted_tree;
+  la->num_interacted_pm += lb->num_interacted_pm;
+#endif
 
   la->interacted = 1;
 
@@ -597,7 +610,7 @@ INLINE static void gravity_multipole_add(struct multipole *restrict ma,
 #error "Missing implementation for order >5"
 #endif
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
   ma->num_gpart += mb->num_gpart;
 #endif
 }
@@ -1319,7 +1332,7 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
 #error "Missing implementation for order >5"
 #endif
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
   multi->m_pole.num_gpart = gcount;
 #endif
 }
@@ -1576,7 +1589,7 @@ INLINE static void gravity_M2M(struct multipole *restrict m_a,
 #error "Missing implementation for order >5"
 #endif
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
   m_a->num_gpart = m_b->num_gpart;
 #endif
 }
@@ -1599,6 +1612,10 @@ INLINE static void gravity_M2L_apply(
   l_b->num_interacted += m_a->num_gpart;
 #endif
 
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  l_b->num_interacted_tree += m_a->num_gpart;
+#endif
+
   /* Record that this tensor has received contributions */
   l_b->interacted = 1;
 
@@ -2091,8 +2108,13 @@ INLINE static void gravity_L2L(struct grav_tensor *restrict la,
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (lb->num_interacted == 0) error("Shifting tensors that did not interact");
+
   la->num_interacted = lb->num_interacted;
 #endif
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  la->num_interacted_tree = lb->num_interacted_tree;
+  la->num_interacted_pm = lb->num_interacted_pm;
+#endif
 
   /* Distance to shift by */
   const double dx[3] = {pos_a[0] - pos_b[0], pos_a[1] - pos_b[1],
@@ -2444,9 +2466,15 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (lb->num_interacted == 0) error("Interacting with empty field tensor");
+
   gp->num_interacted += lb->num_interacted;
 #endif
 
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  gp->num_interacted_m2l += lb->num_interacted_tree;
+  gp->num_interacted_pm += lb->num_interacted_pm;
+#endif
+
   /* Local accumulator */
   double a_grav[3] = {0., 0., 0.};
   double pot = 0.;
@@ -2565,6 +2593,12 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
   gp->a_grav[1] += a_grav[1];
   gp->a_grav[2] += a_grav[2];
   gravity_add_comoving_potential(gp, pot);
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  gp->a_grav_m2l[0] += a_grav[0];
+  gp->a_grav_m2l[1] += a_grav[1];
+  gp->a_grav_m2l[2] += a_grav[2];
+#endif
 }
 
 /**

src/runner_doiact_grav.c

@@ -266,6 +266,12 @@ static INLINE void runner_dopair_grav_pp_full(
       if (pjd < gcount_j && !gpart_is_inhibited(&gparts_j[pjd], e))
         gparts_i[pid].num_interacted++;
 #endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+      /* Update the p2p interaction counter if it's not a padded gpart */
+      if (pjd < gcount_j && !gpart_is_inhibited(&gparts_j[pjd], e))
+        gparts_i[pid].num_interacted_p2p++;
+#endif
     }
 
     /* Store everything back in cache */
@@ -273,6 +279,12 @@ static INLINE void runner_dopair_grav_pp_full(
     ci_cache->a_y[pid] += a_y;
     ci_cache->a_z[pid] += a_z;
     ci_cache->pot[pid] += pot;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    gparts_i[pid].a_grav_p2p[0] += a_x;
+    gparts_i[pid].a_grav_p2p[1] += a_y;
+    gparts_i[pid].a_grav_p2p[2] += a_z;
+#endif
   }
 }
 
@@ -410,6 +422,12 @@ static INLINE void runner_dopair_grav_pp_truncated(
       if (pjd < gcount_j && !gpart_is_inhibited(&gparts_j[pjd], e))
         gparts_i[pid].num_interacted++;
 #endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+      /* Update the p2p interaction counter if it's not a padded gpart */
+      if (pjd < gcount_j && !gpart_is_inhibited(&gparts_j[pjd], e))
+        gparts_i[pid].num_interacted_p2p++;
+#endif
     }
 
     /* Store everything back in cache */
@@ -417,6 +435,12 @@ static INLINE void runner_dopair_grav_pp_truncated(
     ci_cache->a_y[pid] += a_y;
     ci_cache->a_z[pid] += a_z;
     ci_cache->pot[pid] += pot;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    gparts_i[pid].a_grav_p2p[0] += a_x;
+    gparts_i[pid].a_grav_p2p[1] += a_y;
+    gparts_i[pid].a_grav_p2p[2] += a_z;
+#endif
   }
 }
 
@@ -543,6 +567,16 @@ static INLINE void runner_dopair_grav_pm_full(
     if (pid < gcount_i)
       gparts_i[pid].num_interacted += cj->grav.multipole->m_pole.num_gpart;
 #endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    /* Update the M2P interaction counter and forces. */
+    if (pid < gcount_i) {
+      gparts_i[pid].num_interacted_m2p += cj->grav.multipole->m_pole.num_gpart;
+      gparts_i[pid].a_grav_m2p[0] += f_x;
+      gparts_i[pid].a_grav_m2p[1] += f_y;
+      gparts_i[pid].a_grav_m2p[2] += f_z;
+    }
+#endif
   }
 }
 
@@ -674,6 +708,16 @@ static INLINE void runner_dopair_grav_pm_truncated(
     if (pid < gcount_i)
       gparts_i[pid].num_interacted += cj->grav.multipole->m_pole.num_gpart;
 #endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    /* Update the M2P interaction counter and forces. */
+    if (pid < gcount_i) {
+      gparts_i[pid].num_interacted_m2p += cj->grav.multipole->m_pole.num_gpart;
+      gparts_i[pid].a_grav_m2p[0] += f_x;
+      gparts_i[pid].a_grav_m2p[1] += f_y;
+      gparts_i[pid].a_grav_m2p[2] += f_z;
+    }
+#endif
   }
 }
 
@@ -999,6 +1043,12 @@ static INLINE void runner_doself_grav_pp_full(
       if (pjd < gcount && !gpart_is_inhibited(&gparts[pjd], e))
         gparts[pid].num_interacted++;
 #endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+      /* Update the P2P interaction counter if it's not a padded gpart */
+      if (pjd < gcount && !gpart_is_inhibited(&gparts[pjd], e))
+        gparts[pid].num_interacted_p2p++;
+#endif
     }
 
     /* Store everything back in cache */
@@ -1006,6 +1056,12 @@ static INLINE void runner_doself_grav_pp_full(
     ci_cache->a_y[pid] += a_y;
     ci_cache->a_z[pid] += a_z;
     ci_cache->pot[pid] += pot;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    gparts[pid].a_grav_p2p[0] += a_x;
+    gparts[pid].a_grav_p2p[1] += a_y;
+    gparts[pid].a_grav_p2p[2] += a_z;
+#endif
   }
 }
 
@@ -1126,6 +1182,12 @@ static INLINE void runner_doself_grav_pp_truncated(
       if (pjd < gcount && !gpart_is_inhibited(&gparts[pjd], e))
         gparts[pid].num_interacted++;
 #endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+      /* Update the P2P interaction counter if it's not a padded gpart */
+      if (pjd < gcount && !gpart_is_inhibited(&gparts[pjd], e))
+        gparts[pid].num_interacted_p2p++;
+#endif
     }
 
     /* Store everything back in cache */
@@ -1133,6 +1195,12 @@ static INLINE void runner_doself_grav_pp_truncated(
     ci_cache->a_y[pid] += a_y;
     ci_cache->a_z[pid] += a_z;
     ci_cache->pot[pid] += pot;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    gparts[pid].a_grav_p2p[0] += a_x;
+    gparts[pid].a_grav_p2p[1] += a_y;
+    gparts[pid].a_grav_p2p[2] += a_z;
+#endif
   }
 }
 
@@ -1569,12 +1637,19 @@ void runner_dopair_recursive_grav(struct runner *r, struct cell *ci,
   if (periodic && r_lr_check > max_distance) {
 
 #ifdef SWIFT_DEBUG_CHECKS
-    /* Need to account for the interactions we missed */
     if (cell_is_active_gravity(ci, e))
       multi_i->pot.num_interacted += multi_j->m_pole.num_gpart;
     if (cell_is_active_gravity(cj, e))
       multi_j->pot.num_interacted += multi_i->m_pole.num_gpart;
 #endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    /* Need to account for the interactions we missed */
+    if (cell_is_active_gravity(ci, e))
+      multi_i->pot.num_interacted_pm += multi_j->m_pole.num_gpart;
+    if (cell_is_active_gravity(cj, e))
+      multi_j->pot.num_interacted_pm += multi_i->m_pole.num_gpart;
+#endif
     return;
   }
 
@@ -1782,6 +1857,11 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
         multi_i->pot.num_interacted += multi_j->m_pole.num_gpart;
 #endif
 
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+        /* Need to account for the interactions we missed */
+        multi_i->pot.num_interacted_pm += multi_j->m_pole.num_gpart;
+#endif
+
         /* Record that this multipole received a contribution */
         multi_i->pot.interacted = 1;