Added order 4 sofetened gravity.

src/gravity_derivatives.h

@@ -32,6 +32,7 @@
 
 /* Local headers. */
 #include "inline.h"
+#include "kernel_gravity.h"
 
 struct potential_derivatives {
 
@@ -85,29 +86,73 @@ struct potential_derivatives {
 };
 
 __attribute__((always_inline)) INLINE static void compute_potential_derivatives(
-    float r_x, float r_y, float r_z, float r_inv,
-    struct potential_derivatives *pot) {
+    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) {
 
-/* Compute some (odd) powers of 1/r */
+  float Dt_1;
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
-  const float r_inv2 = r_inv * r_inv;
-  const float r_inv3 = -1.f * r_inv * r_inv2; /* -1 / r^3 */
+  float Dt_3;
 #endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 1
-  const float r_inv5 = -3.f * r_inv3 * r_inv2; /* 3 / r^5 */
+  float Dt_5;
 #endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 2
-  const float r_inv7 = -5.f * r_inv5 * r_inv2; /* -15 / r^7 */
+  float Dt_7;
 #endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 3
-  const float r_inv9 = -7.f * r_inv7 * r_inv2; /* 105 / r^9 */
+  float Dt_9;
+#endif
+
+  /* Un-softened case */
+  if (r2 > eps2) {
+
+    Dt_1 = r_inv;
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+    const float r_inv2 = r_inv * r_inv;
+    Dt_3 = -1.f * Dt_1 * r_inv2; /* -1 / r^3 */
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+    Dt_5 = -3.f * Dt_3 * r_inv2; /* 3 / r^5 */
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+    Dt_7 = -5.f * Dt_5 * r_inv2; /* -15 / r^7 */
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+    Dt_9 = -7.f * Dt_7 * r_inv2; /* 105 / r^9 */
 #endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 4
-  const float r_inv11 = -9.f * r_inv9 * r_inv2; /* -945 / r^11 */
+#error "Missing implementation for order >4"
 #endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
-#error "Missing implementation for order >5"
+
+  } else {
+    const float r = r2 * r_inv;
+    const float u = r * eps_inv;
+    const float u_inv = r_inv * eps;
+
+    Dt_1 = eps_inv * D_soft_1(u, u_inv);
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+    const float eps_inv2 = eps_inv * eps_inv;
+    const float eps_inv3 = eps_inv * eps_inv2;
+    Dt_3 = -eps_inv3 * D_soft_3(u, u_inv);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+    const float eps_inv5 = eps_inv3 * eps_inv2;
+    Dt_5 = eps_inv5 * D_soft_5(u, u_inv);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+    const float eps_inv7 = eps_inv5 * eps_inv2;
+    Dt_7 = -eps_inv7 * D_soft_7(u, u_inv);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+    const float eps_inv9 = eps_inv7 * eps_inv2;
+    Dt_9 = eps_inv9 * D_soft_9(u, u_inv);
 #endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for order >4"
+#endif
+  }
+
+/* Alright, let's get the full terms */
 
 /* Compute some powers of r_x, r_y and r_z */
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 1
@@ -135,56 +180,56 @@ __attribute__((always_inline)) INLINE static void compute_potential_derivatives(
 #endif
 
   /* Get the 0th order term */
-  pot->D_000 = r_inv;
+  pot->D_000 = Dt_1;
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
   /* 1st order derivatives */
-  pot->D_100 = r_x * r_inv3;
-  pot->D_010 = r_y * r_inv3;
-  pot->D_001 = r_z * r_inv3;
+  pot->D_100 = r_x * Dt_3;
+  pot->D_010 = r_y * Dt_3;
+  pot->D_001 = r_z * Dt_3;
 #endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 1
   /* 2nd order derivatives */
-  pot->D_200 = r_x2 * r_inv5 + r_inv3;
-  pot->D_020 = r_y2 * r_inv5 + r_inv3;
-  pot->D_002 = r_z2 * r_inv5 + r_inv3;
-  pot->D_110 = r_x * r_y * r_inv5;
-  pot->D_101 = r_x * r_z * r_inv5;
-  pot->D_011 = r_y * r_z * r_inv5;
+  pot->D_200 = r_x2 * Dt_5 + Dt_3;
+  pot->D_020 = r_y2 * Dt_5 + Dt_3;
+  pot->D_002 = r_z2 * Dt_5 + Dt_3;
+  pot->D_110 = r_x * r_y * Dt_5;
+  pot->D_101 = r_x * r_z * Dt_5;
+  pot->D_011 = r_y * r_z * Dt_5;
 #endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 2
   /* 3rd order derivatives */
-  pot->D_300 = r_x3 * r_inv7 + 3.f * r_x * r_inv5;
-  pot->D_030 = r_y3 * r_inv7 + 3.f * r_y * r_inv5;
-  pot->D_003 = r_z3 * r_inv7 + 3.f * r_z * r_inv5;
-  pot->D_210 = r_x2 * r_y * r_inv7 + r_y * r_inv5;
-  pot->D_201 = r_x2 * r_z * r_inv7 + r_z * r_inv5;
-  pot->D_120 = r_y2 * r_x * r_inv7 + r_x * r_inv5;
-  pot->D_021 = r_y2 * r_z * r_inv7 + r_z * r_inv5;
-  pot->D_102 = r_z2 * r_x * r_inv7 + r_x * r_inv5;
-  pot->D_012 = r_z2 * r_y * r_inv7 + r_y * r_inv5;
-  pot->D_111 = r_x * r_y * r_z * r_inv7;
+  pot->D_300 = r_x3 * Dt_7 + 3.f * r_x * Dt_5;
+  pot->D_030 = r_y3 * Dt_7 + 3.f * r_y * Dt_5;
+  pot->D_003 = r_z3 * Dt_7 + 3.f * r_z * Dt_5;
+  pot->D_210 = r_x2 * r_y * Dt_7 + r_y * Dt_5;
+  pot->D_201 = r_x2 * r_z * Dt_7 + r_z * Dt_5;
+  pot->D_120 = r_y2 * r_x * Dt_7 + r_x * Dt_5;
+  pot->D_021 = r_y2 * r_z * Dt_7 + r_z * Dt_5;
+  pot->D_102 = r_z2 * r_x * Dt_7 + r_x * Dt_5;
+  pot->D_012 = r_z2 * r_y * Dt_7 + r_y * Dt_5;
+  pot->D_111 = r_x * r_y * r_z * Dt_7;
 #endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 3
   /* 4th order derivatives */
-  pot->D_400 = r_x4 * r_inv9 + 6.f * r_x2 * r_inv7 + 3.f * r_inv5;
-  pot->D_040 = r_y4 * r_inv9 + 6.f * r_y2 * r_inv7 + 3.f * r_inv5;
-  pot->D_004 = r_z4 * r_inv9 + 6.f * r_z2 * r_inv7 + 3.f * r_inv5;
-  pot->D_310 = r_x3 * r_y * r_inv9 + 3.f * r_x * r_y * r_inv7;
-  pot->D_301 = r_x3 * r_z * r_inv9 + 3.f * r_x * r_z * r_inv7;
-  pot->D_130 = r_y3 * r_x * r_inv9 + 3.f * r_y * r_x * r_inv7;
-  pot->D_031 = r_y3 * r_z * r_inv9 + 3.f * r_y * r_z * r_inv7;
-  pot->D_103 = r_z3 * r_x * r_inv9 + 3.f * r_z * r_x * r_inv7;
-  pot->D_013 = r_z3 * r_y * r_inv9 + 3.f * r_z * r_y * r_inv7;
-  pot->D_220 = r_x2 * r_y2 * r_inv9 + r_x2 * r_inv7 + r_y2 * r_inv7 + r_inv5;
-  pot->D_202 = r_x2 * r_z2 * r_inv9 + r_x2 * r_inv7 + r_z2 * r_inv7 + r_inv5;
-  pot->D_022 = r_y2 * r_z2 * r_inv9 + r_y2 * r_inv7 + r_z2 * r_inv7 + r_inv5;
-  pot->D_211 = r_x2 * r_y * r_z * r_inv9 + r_y * r_z * r_inv7;
-  pot->D_121 = r_y2 * r_x * r_z * r_inv9 + r_x * r_z * r_inv7;
-  pot->D_112 = r_z2 * r_x * r_y * r_inv9 + r_x * r_y * r_inv7;
+  pot->D_400 = r_x4 * Dt_9 + 6.f * r_x2 * Dt_7 + 3.f * Dt_5;
+  pot->D_040 = r_y4 * Dt_9 + 6.f * r_y2 * Dt_7 + 3.f * Dt_5;
+  pot->D_004 = r_z4 * Dt_9 + 6.f * r_z2 * Dt_7 + 3.f * Dt_5;
+  pot->D_310 = r_x3 * r_y * Dt_9 + 3.f * r_x * r_y * Dt_7;
+  pot->D_301 = r_x3 * r_z * Dt_9 + 3.f * r_x * r_z * Dt_7;
+  pot->D_130 = r_y3 * r_x * Dt_9 + 3.f * r_y * r_x * Dt_7;
+  pot->D_031 = r_y3 * r_z * Dt_9 + 3.f * r_y * r_z * Dt_7;
+  pot->D_103 = r_z3 * r_x * Dt_9 + 3.f * r_z * r_x * Dt_7;
+  pot->D_013 = r_z3 * r_y * Dt_9 + 3.f * r_z * r_y * Dt_7;
+  pot->D_220 = r_x2 * r_y2 * Dt_9 + r_x2 * Dt_7 + r_y2 * Dt_7 + Dt_5;
+  pot->D_202 = r_x2 * r_z2 * Dt_9 + r_x2 * Dt_7 + r_z2 * Dt_7 + Dt_5;
+  pot->D_022 = r_y2 * r_z2 * Dt_9 + r_y2 * Dt_7 + r_z2 * Dt_7 + Dt_5;
+  pot->D_211 = r_x2 * r_y * r_z * Dt_9 + r_y * r_z * Dt_7;
+  pot->D_121 = r_y2 * r_x * r_z * Dt_9 + r_x * r_z * Dt_7;
+  pot->D_112 = r_z2 * r_x * r_y * Dt_9 + r_x * r_y * Dt_7;
 #endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
-#error "Missing implementation for order >5"
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+#error "Missing implementation for orders >4"
 #endif
 }
 

src/gravity_softened_derivatives.h

@@ -34,6 +34,8 @@
 #include "inline.h"
 #include "kernel_gravity.h"
 
+#if 0
+
 /*************************/
 /* 0th order derivatives */
 /*************************/
@@ -440,4 +442,6 @@ __attribute__((always_inline)) INLINE static double D_soft_111(
   return -r_x * r_y * r_z * eps_inv7 * D_soft_3(u);
 }
 
+#endif
+
 #endif /* SWIFT_GRAVITY_SOFTENED_DERIVATIVE_H */

src/kernel_gravity.h

@@ -71,46 +71,61 @@ __attribute__((always_inline)) INLINE static void kernel_grav_eval_double(
 /* Derivatives of softening kernel used for FMM */
 /************************************************/
 
-__attribute__((always_inline)) INLINE static double D_soft_0(double u) {
+__attribute__((always_inline)) INLINE static float D_soft_1(float u,
+                                                            float u_inv) {
 
   /* phi(u) = -3u^7 + 15u^6 - 28u^5 + 21u^4 - 7u^2 + 3 */
-  double phi = -3. * u + 15.;
-  phi = phi * u - 28.;
-  phi = phi * u + 21.;
+  float phi = -3.f * u + 15.f;
+  phi = phi * u - 28.f;
+  phi = phi * u + 21.f;
   phi = phi * u;
-  phi = phi * u - 7.;
+  phi = phi * u - 7.f;
   phi = phi * u;
-  phi = phi * u + 3.;
+  phi = phi * u + 3.f;
 
   return phi;
 }
 
-__attribute__((always_inline)) INLINE static double D_soft_1(double u) {
+__attribute__((always_inline)) INLINE static float D_soft_3(float u,
+                                                            float u_inv) {
 
   /* phi'(u)/u = 21u^5 - 90u^4 + 140u^3 - 84u^2 + 14 */
-  double phi = 21. * u - 90.;
-  phi = phi * u + 140.;
-  phi = phi * u - 84.;
+  float phi = 21.f * u - 90.f;
+  phi = phi * u + 140.f;
+  phi = phi * u - 84.f;
   phi = phi * u;
-  phi = phi * u + 14.;
+  phi = phi * u + 14.f;
 
   return phi;
 }
 
-__attribute__((always_inline)) INLINE static double D_soft_2(double u) {
+__attribute__((always_inline)) INLINE static float D_soft_5(float u,
+                                                            float u_inv) {
 
   /* (phi'(u)/u)'/u = -105u^3 + 360u^2 - 420u + 168 */
-  double phi = -105. * u + 360.;
-  phi = phi * u - 420.;
-  phi = phi * u + 168.;
+  float phi = -105.f * u + 360.f;
+  phi = phi * u - 420.f;
+  phi = phi * u + 168.f;
 
   return phi;
 }
 
-__attribute__((always_inline)) INLINE static double D_soft_3(double u) {
+__attribute__((always_inline)) INLINE static float D_soft_7(float u,
+                                                            float u_inv) {
 
-  /* ((phi'(u)/u)'/u)'/u = 315u - 720 + 420/u */
-  return 315. * u - 720. + 420. / u;
+  /* ((phi'(u)/u)'/u)'/u = 315u - 720 + 420u^-1 */
+  return 315.f * u - 720.f + 420.f * u_inv;
+}
+
+__attribute__((always_inline)) INLINE static float D_soft_9(float u,
+                                                            float u_inv) {
+
+  /* (((phi'(u)/u)'/u)'/u)'/u = -315u^-1 + 420u^-3 */
+  float phi = 420.f * u_inv;
+  phi = phi * u_inv - 315.f;
+  phi = phi * u_inv;
+
+  return phi;
 }
 
 #endif /* SWIFT_KERNEL_GRAVITY_H */

tests/testGravityDerivatives.c

@@ -955,7 +955,7 @@ int main() {
 
     /* Compute all derivatives */
     struct potential_derivatives pot;
-    compute_potential_derivatives(dx, dy, dz, r_inv, &pot);
+    compute_potential_derivatives(dx, dy, dz, r2, r_inv, 0., 0., FLT_MAX, &pot);
 
     /* Minimal value we care about */
     const double min = 1e-9;