diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 4a4d5cf7ad6bc191a07d0f3e544f3835ca9b8f5d..06dbb0398287b200a59aa2bc1086e6aa4d51a420 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/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 */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index ebdba8993d32d0c74c08a1de697527dbff5833d3..cec450d84449d2b8cbb6790fcf9a596bb24649be 100644
--- a/src/runner_doiact_grav.h
+++ b/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++) {
diff --git a/tests/testPotentialPair.c b/tests/testPotentialPair.c
index 00ca87e094a62d2e40bec17aafc2e08dfdf60144..a463f19089d976360f90e8a980386ab031dc126b 100644
--- a/tests/testPotentialPair.c
+++ b/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);
diff --git a/tests/testPotentialSelf.c b/tests/testPotentialSelf.c
index 285cf11ecf8d8a44bda87b2e1196aa3baac605cf..33fe574e34bcc49ff8bac6f19f22b55fda19f186 100644
--- a/tests/testPotentialSelf.c
+++ b/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 */