Use floating-point precision throughout the M2L kernel evaluation.

src/gravity_derivatives.h

@@ -85,6 +85,20 @@ struct potential_derivatives {
 #endif
 };
 
+/**
+ * @brief Compute all the relevent derivatives of the softened and truncated
+ * gravitational potential.
+ *
+ * @param r_x x-component of distance vector
+ * @param r_y y-component of distance vector
+ * @param r_z z-component of distance vector
+ * @param r2 Square norm of distance vector
+ * @param r_inv Inverse norm of distance vector
+ * @param eps Softening length.
+ * @param eps2 Square of softening length.
+ * @param eps_inv Inverse of softening length.
+ * @param pot (return) The structure containing all the derivatives.
+ */
 __attribute__((always_inline)) INLINE static void compute_potential_derivatives(
     float r_x, float r_y, float r_z, float r2, float r_inv, float eps,
     float eps2, float eps_inv, struct potential_derivatives *pot) {

src/gravity_properties.c

@@ -59,7 +59,7 @@ void gravity_props_init(struct gravity_props *p,
   /* Softening lengths */
   p->epsilon = 3. * parser_get_param_double(params, "Gravity:epsilon");
   p->epsilon2 = p->epsilon * p->epsilon;
-  p->epsilon_inv = 1. / p->epsilon;
+  p->epsilon_inv = 1.f / p->epsilon;
 }
 
 void gravity_props_print(const struct gravity_props *p) {

src/gravity_properties.h

@@ -58,13 +58,13 @@ struct gravity_props {
   double theta_crit_inv;
 
   /*! Softening length */
-  double epsilon;
+  float epsilon;
 
   /*! Square of softening length */
-  double epsilon2;
+  float epsilon2;
 
   /*! Inverse of softening length */
-  double epsilon_inv;
+  float epsilon_inv;
 };
 
 void gravity_props_print(const struct gravity_props *p);

src/multipole.h

@@ -1516,14 +1516,14 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
                                const double dim[3]) {
 
   /* Recover some constants */
-  const double eps = props->epsilon;
-  const double eps2 = props->epsilon2;
-  const double eps_inv = props->epsilon_inv;
+  const float eps = props->epsilon;
+  const float eps2 = props->epsilon2;
+  const float eps_inv = props->epsilon_inv;
 
   /* Compute distance vector */
-  double dx = pos_b[0] - pos_a[0];
-  double dy = pos_b[1] - pos_a[1];
-  double dz = pos_b[2] - pos_a[2];
+  float dx = (float)(pos_b[0] - pos_a[0]);
+  float dy = (float)(pos_b[1] - pos_a[1]);
+  float dz = (float)(pos_b[2] - pos_a[2]);
 
   /* Apply BC */
   if (periodic) {
@@ -1533,8 +1533,8 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
   }
 
   /* Compute distance */
-  const double r2 = dx * dx + dy * dy + dz * dz;
-  const double r_inv = 1. / sqrt(r2);
+  const float r2 = dx * dx + dy * dy + dz * dz;
+  const float r_inv = 1. / sqrtf(r2);
 
   /* Compute all derivatives */
   struct potential_derivatives pot;