Better way to handle FPEs in the core of the gravity calculation. Vectorises better.

src/gravity/Default/gravity_iact.h

@@ -25,6 +25,9 @@
 #include "kernel_long_gravity.h"
 #include "multipole.h"
 
+/* Standard headers */
+#include <float.h>
+
 /**
  * @brief Computes the intensity of the force at a point generated by a
  * point-mass.
@@ -44,19 +47,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 + FLT_MIN);
+
   /* 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 = sqrtf(r2);
+    const float r = r2 * r_inv;
     const float ui = r * h_inv;
 
     float W_f_ij, W_pot_ij;
@@ -89,24 +92,19 @@ __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) {
 
-  float r;
+  /* Get the inverse distance */
+  const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
+  const float r = r2 * r_inv;
 
   /* 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;