From 52edb2698cedb901bfaef228f352145db6fc5ff1 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Wed, 16 Aug 2017 23:31:15 +0100
Subject: [PATCH] Added 5th order derivatives
---
src/gravity_derivatives.h | 49 +++++++++++++++++++++++++++++++---
src/kernel_gravity.h | 13 +++++++++
tests/testGravityDerivatives.c | 26 ++++++++++++++++++
3 files changed, 85 insertions(+), 3 deletions(-)
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index e808f53454..27aa4f12ae 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -116,6 +116,9 @@ __attribute__((always_inline)) INLINE static void compute_potential_derivatives(
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
float Dt_9;
#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+ float Dt_11;
+#endif
/* Un-softened case */
if (r2 > eps2) {
@@ -135,7 +138,10 @@ __attribute__((always_inline)) INLINE static void compute_potential_derivatives(
Dt_9 = -7.f * Dt_7 * r_inv2; /* 105 / r^9 */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
-#error "Missing implementation for order >4"
+ Dt_11 = -9.f * Dt_9 * r_inv2; /* -945 / r^11 */
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
+#error "Missing implementation for order >5"
#endif
} else {
@@ -162,7 +168,11 @@ __attribute__((always_inline)) INLINE static void compute_potential_derivatives(
Dt_9 = eps_inv9 * D_soft_9(u, u_inv);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
-#error "Missing implementation for order >4"
+ const float eps_inv11 = eps_inv9 * eps_inv2;
+ Dt_11 = -eps_inv11 * D_soft_11(u, u_inv);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
+#error "Missing implementation for order >5"
#endif
}
@@ -243,7 +253,40 @@ __attribute__((always_inline)) INLINE static void compute_potential_derivatives(
pot->D_112 = r_z2 * r_x * r_y * Dt_9 + r_x * r_y * Dt_7;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
-#error "Missing implementation for orders >4"
+ /* 5th order derivatives */
+ pot->D_500 = r_x5 * Dt_11 + 10.f * r_x3 * Dt_9 + 15.f * r_x * Dt_7;
+ pot->D_050 = r_y5 * Dt_11 + 10.f * r_y3 * Dt_9 + 15.f * r_y * Dt_7;
+ pot->D_005 = r_z5 * Dt_11 + 10.f * r_z3 * Dt_9 + 15.f * r_z * Dt_7;
+ pot->D_410 = r_x4 * r_y * Dt_11 + 6.f * r_x2 * r_y * Dt_9 + 3.f * r_y * Dt_7;
+ pot->D_401 = r_x4 * r_z * Dt_11 + 6.f * r_x2 * r_z * Dt_9 + 3.f * r_z * Dt_7;
+ pot->D_140 = r_y4 * r_x * Dt_11 + 6.f * r_y2 * r_x * Dt_9 + 3.f * r_x * Dt_7;
+ pot->D_041 = r_y4 * r_z * Dt_11 + 6.f * r_y2 * r_z * Dt_9 + 3.f * r_z * Dt_7;
+ pot->D_104 = r_z4 * r_x * Dt_11 + 6.f * r_z2 * r_x * Dt_9 + 3.f * r_x * Dt_7;
+ pot->D_014 = r_z4 * r_y * Dt_11 + 6.f * r_z2 * r_y * Dt_9 + 3.f * r_y * Dt_7;
+ pot->D_320 = r_x3 * r_y2 * Dt_11 + r_x3 * Dt_9 + 3.f * r_x * r_y2 * Dt_9 +
+ 3.f * r_x * Dt_7;
+ pot->D_302 = r_x3 * r_z2 * Dt_11 + r_x3 * Dt_9 + 3.f * r_x * r_z2 * Dt_9 +
+ 3.f * r_x * Dt_7;
+ pot->D_230 = r_y3 * r_x2 * Dt_11 + r_y3 * Dt_9 + 3.f * r_y * r_x2 * Dt_9 +
+ 3.f * r_y * Dt_7;
+ pot->D_032 = r_y3 * r_z2 * Dt_11 + r_y3 * Dt_9 + 3.f * r_y * r_z2 * Dt_9 +
+ 3.f * r_y * Dt_7;
+ pot->D_203 = r_z3 * r_x2 * Dt_11 + r_z3 * Dt_9 + 3.f * r_z * r_x2 * Dt_9 +
+ 3.f * r_z * Dt_7;
+ pot->D_023 = r_z3 * r_y2 * Dt_11 + r_z3 * Dt_9 + 3.f * r_z * r_y2 * Dt_9 +
+ 3.f * r_z * Dt_7;
+ pot->D_311 = r_x3 * r_y * r_z * Dt_11 + 3.f * r_x * r_y * r_z * Dt_9;
+ pot->D_131 = r_y3 * r_x * r_z * Dt_11 + 3.f * r_x * r_y * r_z * Dt_9;
+ pot->D_113 = r_z3 * r_x * r_y * Dt_11 + 3.f * r_x * r_y * r_z * Dt_9;
+ pot->D_122 = r_x * r_y2 * r_z2 * Dt_11 + r_x * r_y2 * Dt_9 +
+ r_x * r_z2 * Dt_9 + r_x * Dt_7;
+ pot->D_212 = r_y * r_x2 * r_z2 * Dt_11 + r_y * r_x2 * Dt_9 +
+ r_y * r_z2 * Dt_9 + r_y * Dt_7;
+ pot->D_221 = r_z * r_x2 * r_y2 * Dt_11 + r_z * r_x2 * Dt_9 +
+ r_z * r_y2 * Dt_9 + r_z * Dt_7;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
+#error "Missing implementation for orders >5"
#endif
}
diff --git a/src/kernel_gravity.h b/src/kernel_gravity.h
index 21f65e8a87..799bda85b0 100644
--- a/src/kernel_gravity.h
+++ b/src/kernel_gravity.h
@@ -128,4 +128,17 @@ __attribute__((always_inline)) INLINE static float D_soft_9(float u,
return phi;
}
+__attribute__((always_inline)) INLINE static float D_soft_11(float u,
+ float u_inv) {
+
+ /* ((((phi'(u)/u)'/u)'/u)'/u)'/u = 315u^-3 - 1260u^-5 */
+ float phi = -1260.f * u_inv;
+ phi = phi * u_inv + 315.f;
+ phi = phi * u_inv;
+ phi = phi * u_inv;
+ phi = phi * u_inv;
+
+ return phi;
+}
+
#endif /* SWIFT_KERNEL_GRAVITY_H */
diff --git a/tests/testGravityDerivatives.c b/tests/testGravityDerivatives.c
index 772fec451d..36921d402c 100644
--- a/tests/testGravityDerivatives.c
+++ b/tests/testGravityDerivatives.c
@@ -1016,6 +1016,32 @@ int main() {
test(pot.D_112, D_112(dx, dy, dz, r_inv), tol, min, "D_112");
#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+
+ /* 5th order terms */
+ test(pot.D_500, D_500(dx, dy, dz, r_inv), tol, min, "D_500");
+ test(pot.D_050, D_050(dx, dy, dz, r_inv), tol, min, "D_050");
+ test(pot.D_005, D_005(dx, dy, dz, r_inv), tol, min, "D_005");
+ test(pot.D_410, D_410(dx, dy, dz, r_inv), tol, min, "D_410");
+ test(pot.D_401, D_401(dx, dy, dz, r_inv), tol, min, "D_401");
+ test(pot.D_140, D_140(dx, dy, dz, r_inv), tol, min, "D_140");
+ test(pot.D_041, D_041(dx, dy, dz, r_inv), tol, min, "D_041");
+ test(pot.D_104, D_104(dx, dy, dz, r_inv), tol, min, "D_104");
+ test(pot.D_014, D_014(dx, dy, dz, r_inv), tol, min, "D_014");
+ test(pot.D_320, D_320(dx, dy, dz, r_inv), tol, min, "D_320");
+ test(pot.D_302, D_302(dx, dy, dz, r_inv), tol, min, "D_302");
+ test(pot.D_230, D_230(dx, dy, dz, r_inv), tol, min, "D_230");
+ test(pot.D_032, D_032(dx, dy, dz, r_inv), tol, min, "D_032");
+ test(pot.D_203, D_203(dx, dy, dz, r_inv), tol, min, "D_203");
+ test(pot.D_023, D_023(dx, dy, dz, r_inv), tol, min, "D_023");
+ test(pot.D_311, D_311(dx, dy, dz, r_inv), tol, min, "D_311");
+ test(pot.D_131, D_131(dx, dy, dz, r_inv), tol, min, "D_131");
+ test(pot.D_113, D_113(dx, dy, dz, r_inv), tol, min, "D_113");
+ test(pot.D_122, D_122(dx, dy, dz, r_inv), tol, min, "D_122");
+ test(pot.D_212, D_212(dx, dy, dz, r_inv), tol, min, "D_212");
+ test(pot.D_221, D_221(dx, dy, dz, r_inv), tol, min, "D_221");
+
+#endif
message("All good!");
}
return 0;
--
GitLab