Implemented potential calculation for M2P interaction up to quadrupole order....

Implemented potential calculation for M2P interaction up to quadrupole order. Still need to check accuracy.

src/gravity/Default/gravity_iact.h

@@ -144,33 +144,37 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
     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) {
 
-  //#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
+#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,
                            &f_ij, pot);
   *f_x = -f_ij * r_x;
   *f_y = -f_ij * r_y;
   *f_z = -f_ij * r_z;
-#if 0  // else
+#else
 
   /* Get the inverse distance */
   const float r_inv = 1.f / sqrtf(r2);
 
-  struct potential_derivatives_M2P pot;
-  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, &pot);
+  struct potential_derivatives_M2P d;
+  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, &d);
 
   /* 1st order terms (monopole) */
-  *f_x = m->M_000 * pot.D_100;
-  *f_y = m->M_000 * pot.D_010;
-  *f_z = m->M_000 * pot.D_001;
+  *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 * pot.D_300 + m->M_020 * pot.D_120 + m->M_002 * pot.D_102;
-  *f_y += m->M_200 * pot.D_210 + m->M_020 * pot.D_030 + m->M_002 * pot.D_012;
-  *f_z += m->M_200 * pot.D_201 + m->M_020 * pot.D_021 + m->M_002 * pot.D_003;
-  *f_x += m->M_110 * pot.D_210 + m->M_101 * pot.D_201 + m->M_011 * pot.D_111;
-  *f_y += m->M_110 * pot.D_120 + m->M_101 * pot.D_111 + m->M_011 * pot.D_021;
-  *f_z += m->M_110 * pot.D_111 + m->M_101 * pot.D_102 + m->M_011 * pot.D_012;
+  *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;
 
 #endif
 }

src/gravity_derivatives.h

@@ -95,9 +95,16 @@ struct potential_derivatives_M2L {
  */
 struct potential_derivatives_M2P {
 
+  /* 0th order terms */
+  float D_000;
+
   /* 1st order terms */
   float D_100, D_010, D_001;
 
+  /* 2nd order terms */
+  float D_200, D_020, D_002;
+  float D_110, D_101, D_011;
+
   /* 3rd order terms */
   float D_300, D_030, D_003;
   float D_210, D_201;
@@ -368,11 +375,22 @@ compute_potential_derivatives_M2P(float r_x, float r_y, float r_z, float r2,
   const float r_y3 = r_y2 * r_y;
   const float r_z3 = r_z2 * r_z;
 
+  /* 0th order derivative */
+  pot->D_000 = Dt_1;
+
   /* 1st order derivatives */
   pot->D_100 = r_x * Dt_3;
   pot->D_010 = r_y * Dt_3;
   pot->D_001 = r_z * Dt_3;
 
+  /* 2nd order derivatives */
+  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;
+
   /* 3rd order derivatives */
   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;

tests/testPotentialPair.c

@@ -229,11 +229,11 @@ int main() {
   /* Reset the accelerations */
   for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]);
 
-  /*********************************/
-  /* Now, test the PM interactions */
-  /*********************************/
+  /**********************************/
+  /* Test the basic PM interactions */
+  /**********************************/
 
-  /* Set an opening angle that allows M-P interactions */
+  /* Set an opening angle that allows P-M interactions */
   props.theta_crit2 = 1.;
 
   ci.gparts[0].mass = 0.;
@@ -266,7 +266,98 @@ int main() {
     check_value(gp->a_grav[0], acc_true, "acceleration");
   }
 
-  message("\n\t\t P-M interactions all good\n");
+  message("\n\t\t basic P-M interactions all good\n");
+
+  /* Reset the accelerations */
+  for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]);
+
+/***************************************/
+/* Test the high-order PM interactions */
+/***************************************/
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER >= 3
+
+  /* Let's make ci more interesting */
+  free(ci.gparts);
+  ci.gcount = 8;
+  if (posix_memalign((void **)&ci.gparts, gpart_align,
+                     ci.gcount * sizeof(struct gpart)) != 0)
+    error("Error allocating gparts for cell ci");
+  bzero(ci.gparts, ci.gcount * sizeof(struct gpart));
+
+  /* Place particles on a simple cube of side-length 0.2 */
+  for (int n = 0; n < 8; ++n) {
+    if (n & 1)
+      ci.gparts[n].x[0] = 0.0 - 0.1;
+    else
+      ci.gparts[n].x[0] = 0.0 + 0.1;
+
+    if (n & 2)
+      ci.gparts[n].x[1] = 0.5 - 0.1;
+    else
+      ci.gparts[n].x[1] = 0.5 + 0.1;
+
+    if (n & 2)
+      ci.gparts[n].x[2] = 0.5 - 0.1;
+    else
+      ci.gparts[n].x[2] = 0.5 + 0.1;
+
+    ci.gparts[n].mass = 1. / 8.;
+
+    ci.gparts[n].epsilon = eps;
+    ci.gparts[n].time_bin = 1;
+    ci.gparts[n].type = swift_type_dark_matter;
+    ci.gparts[n].id_or_neg_offset = 1;
+#ifdef SWIFT_DEBUG_CHECKS
+    ci.gparts[n].ti_drift = 8;
+#endif
+  }
+
+  /* Now let's make a multipole out of it. */
+  gravity_reset(ci.multipole);
+  gravity_P2M(ci.multipole, ci.gparts, ci.gcount);
+
+  // message("CoM=[%e %e %e]", ci.multipole->CoM[0], ci.multipole->CoM[1],
+  // ci.multipole->CoM[2]);
+  gravity_multipole_print(&ci.multipole->m_pole);
+
+  /* Compute the forces */
+  runner_dopair_grav_pp(&r, &ci, &cj);
+
+  /* Verify everything */
+  for (int n = 0; n < num_tests; ++n) {
+    const struct gpart *gp = &cj.gparts[n];
+    const struct gravity_tensors *mpole = ci.multipole;
+
+    double pot_true = 0, acc_true[3] = {0., 0., 0.};
+
+    for (int i = 0; i < 8; ++i) {
+      const struct gpart *gp2 = &ci.gparts[i];
+
+      const double dx[3] = {gp2->x[0] - gp->x[0], gp2->x[1] - gp->x[1],
+                            gp2->x[2] - gp->x[2]};
+      const double d = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
+
+      pot_true += potential(gp2->mass, d, gp->epsilon, rlr * FLT_MAX);
+      acc_true[0] -=
+          acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[0] / d;
+      acc_true[1] -=
+          acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[1] / d;
+      acc_true[2] -=
+          acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[2] / d;
+    }
+
+    message("x=%e f=%e f_true=%e pot=%e pot_true=%e %e %e",
+            gp->x[0] - mpole->CoM[0], gp->a_grav[0], acc_true[0], gp->potential,
+            pot_true, acc_true[1], acc_true[2]);
+
+    // check_value(gp->potential, pot_true, "potential");
+    // check_value(gp->a_grav[0], acc_true[0], "acceleration");
+  }
+
+  message("\n\t\t high-order P-M interactions all good\n");
+
+#endif
 
   free(ci.multipole);
   free(cj.multipole);