diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index e73e1b55fec95855e9b40a4eece18465a5c88b0d..b1073e55e735f8425736148f1e00435d7248761c 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/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 */
diff --git a/tests/testPotentialPair.c b/tests/testPotentialPair.c
index b180efe159aec516e2c11948bc4738d266d17d61..57d68264990bf6d3da2ea385126ef5daa341b787 100644
--- a/tests/testPotentialPair.c
+++ b/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");