testPotentialPair now verifies that the potential is correctly computed for...

testPotentialPair now verifies that the potential is correctly computed for the P-M interaction up to 3rd order.

src/gravity/Default/gravity_iact.h

@@ -38,6 +38,7 @@
  * @param h_inv3 Cube of the inverse of the softening length.
  * @param mass Mass of the point-mass.
  * @param f_ij (return) The force intensity.
+ * @param pot_ij (return) The potential.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
     float r2, float h2, float h_inv, float h_inv3, float mass, float *f_ij,
@@ -82,6 +83,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
  * @param mass Mass of the point-mass.
  * @param rlr_inv Inverse of the mesh smoothing scale.
  * @param f_ij (return) The force intensity.
+ * @param pot_ij (return) The potential.
  */
 __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
     float r2, float h2, float h_inv, float h_inv3, float mass, float rlr_inv,
@@ -146,9 +148,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
   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;
-  *f_y = f_ij;
-  *f_z = f_ij;
+  *f_x = -f_ij * r_x;
+  *f_y = -f_ij * r_y;
+  *f_z = -f_ij * r_z;
 #if 0  // else
 
   /* Get the inverse distance */

src/runner_doiact_grav.h

@@ -183,7 +183,7 @@ static INLINE void runner_dopair_grav_pp_full(const struct engine *e,
                                               struct gpart *restrict gparts_j) {
 
   TIMER_TIC;
-  
+
   /* Loop over all particles in ci... */
   for (int pid = 0; pid < gcount_i; pid++) {
 
@@ -277,7 +277,7 @@ static INLINE void runner_dopair_grav_pp_truncated(
     struct gpart *restrict gparts_j) {
 
   TIMER_TIC;
-   
+
   /* Loop over all particles in ci... */
   for (int pid = 0; pid < gcount_i; pid++) {
 

tests/testPotentialPair.c

@@ -96,7 +96,7 @@ int main() {
   e.time_base = 1e-10;
 
   struct space s;
-  s.periodic = 1;
+  s.periodic = 0;
   s.width[0] = 0.2;
   s.width[1] = 0.2;
   s.width[2] = 0.2;
@@ -112,7 +112,7 @@ int main() {
   bzero(&r, sizeof(struct runner));
   r.e = &e;
 
-  const double rlr = props.a_smooth * s.width[0];
+  const double rlr = props.a_smooth * s.width[0] * FLT_MAX;
 
   /* Init the cache for gravity interaction */
   gravity_cache_init(&r.ci_gravity_cache, num_tests);
@@ -123,7 +123,7 @@ int main() {
   struct cell ci, cj;
   bzero(&ci, sizeof(struct cell));
   bzero(&cj, sizeof(struct cell));
-  
+
   ci.width[0] = 1.;
   ci.width[1] = 1.;
   ci.width[2] = 1.;
@@ -157,28 +157,19 @@ int main() {
   /* Set the multipoles */
   ci.multipole->r_max = 0.1;
   cj.multipole->r_max = 0.1;
-  
+
   /* Allocate the particles */
-  posix_memalign((void **)&ci.gparts, gpart_align, ci.gcount * sizeof(struct gpart));
+  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));
 
-  posix_memalign((void **)&cj.gparts, gpart_align, cj.gcount * sizeof(struct gpart));
+  if (posix_memalign((void **)&cj.gparts, gpart_align,
+                     cj.gcount * sizeof(struct gpart)) != 0)
+    error("Error allocating gparts for cell ci");
   bzero(cj.gparts, cj.gcount * sizeof(struct gpart));
-  
-  /* Create the massive particle */
-  ci.gparts[0].x[0] = 1.;
-  ci.gparts[0].x[1] = 0.5;
-  ci.gparts[0].x[2] = 0.5;
-  ci.gparts[0].mass = 1.;
-  ci.gparts[0].epsilon = eps;
-  ci.gparts[0].time_bin = 1;
-  ci.gparts[0].type = swift_type_dark_matter;
-  ci.gparts[0].id_or_neg_offset = 1;
-#ifdef SWIFT_DEBUG_CHECKS
-  ci.gparts[0].ti_drift = 8;
-#endif
 
-  /* Create the mass-less particles */
+  /* Create the mass-less test particles */
   for (int n = 0; n < num_tests; ++n) {
 
     struct gpart *gp = &cj.gparts[n];
@@ -196,6 +187,23 @@ int main() {
 #endif
   }
 
+  /***********************************************/
+  /* Let's start by testing the P-P interactions */
+  /***********************************************/
+
+  /* Create the massive particle */
+  ci.gparts[0].x[0] = 0.;
+  ci.gparts[0].x[1] = 0.5;
+  ci.gparts[0].x[2] = 0.5;
+  ci.gparts[0].mass = 1.;
+  ci.gparts[0].epsilon = eps;
+  ci.gparts[0].time_bin = 1;
+  ci.gparts[0].type = swift_type_dark_matter;
+  ci.gparts[0].id_or_neg_offset = 1;
+#ifdef SWIFT_DEBUG_CHECKS
+  ci.gparts[0].ti_drift = 8;
+#endif
+
   /* Now compute the forces */
   runner_dopair_grav_pp(&r, &ci, &cj);
 
@@ -204,16 +212,62 @@ int main() {
     const struct gpart *gp = &cj.gparts[n];
     const struct gpart *gp2 = &ci.gparts[0];
 
-    double pot_true = potential(ci.gparts[0].mass, gp->x[0] - gp2->x[0], gp->epsilon, rlr);
+    double pot_true =
+        potential(ci.gparts[0].mass, gp->x[0] - gp2->x[0], gp->epsilon, rlr);
     double acc_true =
         acceleration(ci.gparts[0].mass, gp->x[0] - gp2->x[0], gp->epsilon, rlr);
 
-    message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - gp2->x[0], gp->a_grav[0], acc_true, gp->potential, pot_true);
+    message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - gp2->x[0],
+            gp->a_grav[0], acc_true, gp->potential, pot_true);
 
     check_value(gp->potential, pot_true, "potential");
     check_value(gp->a_grav[0], acc_true, "acceleration");
   }
 
+  message("\n\t\t P-P interactions all good\n");
+
+  /* Reset the accelerations */
+  for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]);
+
+  /*********************************/
+  /* Now, test the PM interactions */
+  /*********************************/
+
+  /* Set an opening angle that allows M-P interactions */
+  props.theta_crit2 = 1.;
+
+  ci.gparts[0].mass = 0.;
+  ci.multipole->CoM[0] = 0.;
+  ci.multipole->CoM[1] = 0.5;
+  ci.multipole->CoM[2] = 0.5;
+
+  bzero(&ci.multipole->m_pole, sizeof(struct multipole));
+  bzero(&cj.multipole->m_pole, sizeof(struct multipole));
+  ci.multipole->m_pole.M_000 = 1.;
+
+  /* Now 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 = potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0],
+                                gp->epsilon, rlr * FLT_MAX);
+    double acc_true =
+        acceleration(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], gp->epsilon,
+                     rlr * FLT_MAX);
+
+    message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - mpole->CoM[0],
+            gp->a_grav[0], acc_true, gp->potential, pot_true);
+
+    check_value(gp->potential, pot_true, "potential");
+    check_value(gp->a_grav[0], acc_true, "acceleration");
+  }
+
+  message("\n\t\t P-M interactions all good\n");
+
   free(ci.multipole);
   free(cj.multipole);
   free(ci.gparts);

tests/testPotentialSelf.c

@@ -129,8 +129,9 @@ int main() {
   c.ti_gravity_end_min = 8;
   c.ti_gravity_end_max = 8;
 
-  posix_memalign((void **)&c.gparts, gpart_align,
-                 c.gcount * sizeof(struct gpart));
+  if (posix_memalign((void **)&c.gparts, gpart_align,
+                     c.gcount * sizeof(struct gpart)) != 0)
+    error("Impossible to allocate memory for the gparts.");
   bzero(c.gparts, c.gcount * sizeof(struct gpart));
 
   /* Create the massive particle */