testPotentialPair now also tests order 3 and more M2P interactions.

src/gravity/Default/gravity_iact.h

@@ -152,7 +152,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
 /* 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 1  // SELF_GRAVITY_MULTIPOLE_ORDER < 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
 
   float f_ij, pot_ij;
   runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
@@ -176,18 +176,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   *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;
+  *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;
+  *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;
+  *pot += m->M_110 * d.D_110 + m->M_101 * d.D_101 + m->M_011 * d.D_011;
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 3
 
@@ -195,7 +195,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   *f_x += m->M_300 * d.D_400 + m->M_030 * d.D_130 + m->M_003 * d.D_103;
   *f_y += m->M_300 * d.D_310 + m->M_030 * d.D_040 + m->M_003 * d.D_013;
   *f_z += m->M_300 * d.D_301 + m->M_030 * d.D_031 + m->M_003 * d.D_004;
-  *pot -= m->M_300 * d.D_100 + m->M_030 * d.D_030 + m->M_003 * d.D_003;
+  *pot += m->M_300 * d.D_100 + m->M_030 * d.D_030 + m->M_003 * d.D_003;
 #endif
 
   /* Take care of the the sign convention */
@@ -235,7 +235,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
 /* 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 1  // SELF_GRAVITY_MULTIPOLE_ORDER < 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
 
   float f_ij, pot_ij;
   runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
@@ -265,12 +265,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
   *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;
+  *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;
+  *pot += m->M_110 * d.D_110 + m->M_101 * d.D_101 + m->M_011 * d.D_011;
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 3
 
@@ -278,7 +278,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
   *f_x += m->M_300 * d.D_400 + m->M_030 * d.D_130 + m->M_003 * d.D_103;
   *f_y += m->M_300 * d.D_310 + m->M_030 * d.D_040 + m->M_003 * d.D_013;
   *f_z += m->M_300 * d.D_301 + m->M_030 * d.D_031 + m->M_003 * d.D_004;
-  *pot -= m->M_300 * d.D_100 + m->M_030 * d.D_030 + m->M_003 * d.D_003;
+  *pot += m->M_300 * d.D_100 + m->M_030 * d.D_030 + m->M_003 * d.D_003;
 #endif
 
   /* Take care of the the sign convention */

tests/testPotentialPair.c

@@ -39,12 +39,19 @@ const double eps = 0.1;
  * @param a First value
  * @param b Second value
  * @param s String used to identify this check in messages
+ * @param rel_tol Maximal relative error
+ * @param limit Minimal value to consider in the tests
  */
-void check_value(double a, double b, const char *s) {
-  if (fabs(a - b) / fabs(a + b) > 2e-6 && fabs(a - b) > 1.e-6)
+void check_value_backend(double a, double b, const char *s, double rel_tol,
+                         double limit) {
+  if (fabs(a - b) / fabs(a + b) > rel_tol && fabs(a - b) > limit)
     error("Values are inconsistent: SWIFT:%12.15e true:%12.15e (%s)!", a, b, s);
 }
 
+void check_value(double a, double b, const char *s) {
+  check_value_backend(a, b, s, 2e-6, 1e-6);
+}
+
 /* Definitions of the potential and force that match
    exactly the theory document */
 double S(double x) { return exp(x) / (1. + exp(x)); }
@@ -319,14 +326,14 @@ int main(int argc, char *argv[]) {
 
   message("\n\t\t truncated P-M interactions all good\n");
 
-  /***************************************/
-  /* Test the high-order PM interactions */
-  /***************************************/
+  /************************************************/
+  /* Test the high-order periodic PM interactions */
+  /************************************************/
 
   /* Reset the accelerations */
   for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]);
 
-#if 0  // SELF_GRAVITY_MULTIPOLE_ORDER >= 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER >= 3
 
   /* Let's make ci more interesting */
   free(ci.gparts);
@@ -367,8 +374,6 @@ int main(int argc, char *argv[]) {
   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 */
@@ -377,7 +382,6 @@ int main(int argc, char *argv[]) {
   /* 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.};
 
@@ -395,12 +399,14 @@ int main(int argc, char *argv[]) {
       acc_true[2] -= acceleration(gp2->mass, d, epsilon, rlr) * 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]);
+    /* const struct gravity_tensors *mpole = ci.multipole; */
+    /* 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");
+    check_value_backend(gp->potential, pot_true, "potential", 1e-2, 1e-6);
+    check_value_backend(gp->a_grav[0], acc_true[0], "acceleration", 1e-2, 1e-6);
   }
 
   message("\n\t\t high-order P-M interactions all good\n");