Prevent FPEs in the P-P gravity calculation.

src/gravity/Default/gravity_iact.h

@@ -44,19 +44,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
     float r2, float h2, float h_inv, float h_inv3, float mass, float *f_ij,
     float *pot_ij) {
 
-  /* Get the inverse distance */
-  const float r_inv = 1.f / sqrtf(r2);
-
   /* Should we soften ? */
   if (r2 >= h2) {
 
+    /* Get the inverse distance */
+    const float r_inv = 1.f / sqrtf(r2);
+
     /* Get Newtonian gravity */
     *f_ij = mass * r_inv * r_inv * r_inv;
     *pot_ij = -mass * r_inv;
 
   } else {
 
-    const float r = r2 * r_inv;
+    const float r = sqrtf(r2);
     const float ui = r * h_inv;
 
     float W_f_ij, W_pot_ij;
@@ -89,19 +89,24 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
     float r2, float h2, float h_inv, float h_inv3, float mass, float rlr_inv,
     float *f_ij, float *pot_ij) {
 
-  /* Get the inverse distance */
-  const float r_inv = 1.f / sqrtf(r2);
-  const float r = r2 * r_inv;
+  float r;
 
   /* Should we soften ? */
   if (r2 >= h2) {
 
+    /* Get the inverse distance */
+    const float r_inv = 1.f / sqrtf(r2);
+    r = r2 * r_inv;
+
     /* Get Newtonian gravity */
     *f_ij = mass * r_inv * r_inv * r_inv;
     *pot_ij = -mass * r_inv;
 
   } else {
 
+    /* Get the distance */
+    r = sqrtf(r2);
+
     const float ui = r * h_inv;
     float W_f_ij, W_pot_ij;
 

src/runner_doiact_grav.h

@@ -234,7 +234,8 @@ static INLINE void runner_dopair_grav_pp_full(const struct engine *e,
       const float r2 = dx * dx + dy * dy + dz * dz;
 
 #ifdef SWIFT_DEBUG_CHECKS
-      if (r2 == 0.f) error("Interacting particles with 0 distance");
+      if (r2 == 0.f && h_i == 0.)
+        error("Interacting particles with 0 distance and 0 softening.");
 
       /* Check that particles have been drifted to the current time */
       if (gparts_i[pid].ti_drift != e->ti_current)