From 5e01d104c96a83bd60caad18adb0ce4ec45a7ca1 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Mon, 11 Jun 2018 17:48:04 +0200
Subject: [PATCH] Add quadrupole term in the M2P kernel.

---
 src/gravity/Default/gravity_iact.h | 48 +++++++++++++++++++++++++++++-
 src/gravity_derivatives.h          | 38 ++++++++++++++++++-----
 2 files changed, 77 insertions(+), 9 deletions(-)

diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 167b4fe3de..28a6855450 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -218,6 +218,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
     float r_x, float r_y, float r_z, float r2, float h, float h_inv,
     const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) {
 
+  /* In the case where the order is < 3, then there is only a monopole term left. */
+  /* We can default to the normal P-P interaction with the mass of the multipole */
+  /* and its CoM as the "particle" property */
 #if SELF_GRAVITY_MULTIPOLE_ORDER < 3
   float f_ij;
   runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
@@ -230,8 +233,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   /* Get the inverse distance */
   const float r_inv = 1.f / sqrtf(r2);
 
+  /* Compute the derivatives of the potential */
   struct potential_derivatives_M2P d;
-  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, &d);
+  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, FLT_MAX, &d);
 
   /* 1st order terms (monopole) */
   *f_x = m->M_000 * d.D_100;
@@ -250,6 +254,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   *f_z += m->M_110 * d.D_111 + m->M_101 * d.D_102 + m->M_011 * d.D_012;
   *pot -= m->M_110 * d.D_110 + m->M_101 * d.D_101 + m->M_011 * d.D_011;
 
+  /* Take care of the the sign convention */
+  *f_x *= -1.f;
+  *f_y *= -1.f;
+  *f_z *= -1.f;
+  *pot *= -1.f;
 #endif
 }
 
@@ -258,12 +267,49 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
     float rlr_inv, const struct multipole *m, float *f_x, float *f_y,
     float *f_z, float *pot) {
 
+  /* In the case where the order is < 3, then there is only a monopole term left. */
+  /* We can default to the normal P-P interaction with the mass of the multipole */
+  /* and its CoM as the "particle" property */
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
   float f_ij;
   runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
                                 m->M_000, rlr_inv, &f_ij, pot);
   *f_x = f_ij * r_x;
   *f_y = f_ij * r_y;
   *f_z = f_ij * r_z;
+
+#else
+
+  /* Get the inverse distance */
+  const float r_inv = 1.f / sqrtf(r2);
+
+  /* Compute the derivatives of the potential */
+  struct potential_derivatives_M2P d;
+  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1, rlr_inv, &d);
+
+  /* 1st order terms (monopole) */
+  *f_x = m->M_000 * d.D_100;
+  *f_y = m->M_000 * d.D_010;
+  *f_z = m->M_000 * d.D_001;
+  *pot = m->M_000 * d.D_000;
+
+  /* 3rd order terms (quadrupole) */
+  *f_x += m->M_200 * d.D_300 + m->M_020 * d.D_120 + m->M_002 * d.D_102;
+  *f_y += m->M_200 * d.D_210 + m->M_020 * d.D_030 + m->M_002 * d.D_012;
+  *f_z += m->M_200 * d.D_201 + m->M_020 * d.D_021 + m->M_002 * d.D_003;
+  *pot -= m->M_200 * d.D_100 + m->M_020 * d.D_020 + m->M_002 * d.D_002;
+
+  *f_x += m->M_110 * d.D_210 + m->M_101 * d.D_201 + m->M_011 * d.D_111;
+  *f_y += m->M_110 * d.D_120 + m->M_101 * d.D_111 + m->M_011 * d.D_021;
+  *f_z += m->M_110 * d.D_111 + m->M_101 * d.D_102 + m->M_011 * d.D_012;
+  *pot -= m->M_110 * d.D_110 + m->M_101 * d.D_101 + m->M_011 * d.D_011;
+
+  /* Take care of the the sign convention */
+  *f_x *= -1.f;
+  *f_y *= -1.f;
+  *f_z *= -1.f;
+  *pot *= -1.f;
+#endif
 }
 
 #endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index b2c27b095e..291b437a23 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -96,7 +96,7 @@ struct potential_derivatives_M2L {
  */
 struct potential_derivatives_M2P {
 
-  /* 0th order terms */
+  /* 0th order term */
   float D_000;
 
   /* 1st order terms */
@@ -130,9 +130,9 @@ struct potential_derivatives_M2P {
  * @param pot (return) The structure containing all the derivatives.
  */
 __attribute__((always_inline)) INLINE static void
-compute_potential_derivatives_M2L(float r_x, float r_y, float r_z, float r2,
-                                  float r_inv, float eps, float eps_inv,
-                                  int periodic, float rs_inv,
+compute_potential_derivatives_M2L(const float r_x, const float r_y, const float r_z, const float r2,
+                                  const float r_inv, const float eps, const float eps_inv,
+                                  const int periodic, const float rs_inv,
                                   struct potential_derivatives_M2L *pot) {
 
   float Dt_1;
@@ -379,11 +379,14 @@ compute_potential_derivatives_M2L(float r_x, float r_y, float r_z, float r2,
  * @param r_inv Inverse norm of distance vector
  * @param eps Softening length.
  * @param eps_inv Inverse of softening length.
+ * @param periodic Is the calculation using periodic BCs?
+ * @param rs_inv The inverse of the gravity mesh-smoothing scale.
  * @param pot (return) The structure containing all the derivatives.
  */
 __attribute__((always_inline)) INLINE static void
-compute_potential_derivatives_M2P(float r_x, float r_y, float r_z, float r2,
-                                  float r_inv, float eps, float eps_inv,
+compute_potential_derivatives_M2P(const float r_x, const float r_y, const float r_z, const float r2,
+                                  const float r_inv, const float eps, const float eps_inv,
+				  const int periodic, const float rs_inv,
                                   struct potential_derivatives_M2P *pot) {
 
   float Dt_1;
@@ -391,8 +394,8 @@ compute_potential_derivatives_M2P(float r_x, float r_y, float r_z, float r2,
   float Dt_5;
   float Dt_7;
 
-  /* Un-softened case */
-  if (r2 > eps * eps) {
+  /* Un-softened un-truncated case (Newtonian potential) */
+  if (!periodic && r2 > eps * eps) {
 
     const float r_inv2 = r_inv * r_inv;
 
@@ -401,6 +404,25 @@ compute_potential_derivatives_M2P(float r_x, float r_y, float r_z, float r2,
     Dt_5 = -3.f * Dt_3 * r_inv2; /* 3 / r^5 */
     Dt_7 = -5.f * Dt_5 * r_inv2; /* -15 / r^7 */
 
+    /* Un-softened truncated case */
+  } else if (periodic && r2 > eps * eps) {
+
+    /* Get the derivatives of the truncated potential */
+    const float r = r2 * r_inv;
+    struct truncated_derivatives d;
+    kernel_long_grav_derivatives(r, rs_inv, &d);
+
+    const float r_inv2 = r_inv * r_inv;
+    const float r_inv3 = r_inv2 * r_inv;
+    const float r_inv5 = r_inv2 * r_inv3;
+    const float r_inv7 = r_inv2 * r_inv5;
+
+    Dt_1 = d.chi_0 * r_inv;
+    Dt_3 = (r * d.chi_1 - d.chi_0) * r_inv3;
+    Dt_5 = (r * r * d.chi_2 - 3.f * r * d.chi_1 + 3.f * d.chi_0) * r_inv5;
+    Dt_7 = (r * r * r * d.chi_3 - 6.f * r * r * d.chi_2 + 15.f * r * d.chi_1 - 15.f * d.chi_0) * r_inv7;
+
+    /* Softened case */
   } else {
 
     const float r = r2 * r_inv;
-- 
GitLab