diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 8f90db42ae0f0f4e5e24c56c427f8f58963eb233..e73e1b55fec95855e9b40a4eece18465a5c88b0d 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 1 // 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,
@@ -240,9 +240,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
m->M_000, r_s_inv, &f_ij, &pot_ij);
- *f_x = -f_ij * r_x;
- *f_y = -f_ij * r_y;
- *f_z = -f_ij * r_z;
+ *f_x = f_ij * r_x;
+ *f_y = f_ij * r_y;
+ *f_z = f_ij * r_z;
*pot = pot_ij;
#else
diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h
index f11f760022d326cafd2ae5da2bc535273b8cc947..2c14c0d888552e91c22e801d88a99ff0f529f341 100644
--- a/src/kernel_long_gravity.h
+++ b/src/kernel_long_gravity.h
@@ -31,7 +31,7 @@
#include <float.h>
#include <math.h>
-#define GADGET2_LONG_RANGE_CORRECTION
+//#define GADGET2_LONG_RANGE_CORRECTION
#ifdef GADGET2_LONG_RANGE_CORRECTION
#define kernel_long_gravity_truncation_name "Gadget-like (using erfc())"
@@ -121,7 +121,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_derivatives(
const float x = c1 * r;
/* e^(2r / r_s) */
- const float exp_x = good_approx_expf(x);
+ const float exp_x = expf(x); // good_approx_expf(x);
/* 1 / alpha(w) */
const float a_inv = 1.f + exp_x;
@@ -164,7 +164,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_pot_eval(
#else
const float x = 2.f * u;
- const float exp_x = good_approx_expf(x);
+ const float exp_x = expf(x); // good_approx_expf(x);
const float alpha = 1.f / (1.f + exp_x);
/* We want 2 - 2 exp(x) * alpha */
@@ -197,7 +197,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval(
#else
const float x = 2.f * u;
- const float exp_x = good_approx_expf(x);
+ const float exp_x = expf(x); // good_approx_expf(x);
const float alpha = 1.f / (1.f + exp_x);
/* We want 2*(x*alpha - x*alpha^2 - exp(x)*alpha + 1) */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 3d5aa78bf1cb7a84ecf548dc9cc914fecf0bd54d..2c12d81223c2e399a953eeaa665a27210803a750 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -751,10 +751,9 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
ci->gparts, cj->gparts);
/* Then the M2P */
- if (0)
- runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j,
- multi_j, dim, r_s_inv, e, ci->gparts,
- gcount_i, cj);
+ runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j,
+ multi_j, dim, r_s_inv, e, ci->gparts,
+ gcount_i, cj);
}
if (cj_active && symmetric) {
@@ -764,10 +763,9 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
cj->gparts, ci->gparts);
/* Then the M2P */
- if (0)
- runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i,
- multi_i, dim, r_s_inv, e, cj->gparts,
- gcount_j, ci);
+ runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i,
+ multi_i, dim, r_s_inv, e, cj->gparts,
+ gcount_j, ci);
}
} else {
@@ -783,10 +781,8 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
ci->gparts, cj->gparts);
/* Then the M2P */
- if (0)
- runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
- periodic, dim, e, ci->gparts, gcount_i,
- cj);
+ runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
+ periodic, dim, e, ci->gparts, gcount_i, cj);
}
if (cj_active && symmetric) {
@@ -796,10 +792,8 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
cj->gparts, ci->gparts);
/* Then the M2P */
- if (0)
- runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
- periodic, dim, e, cj->gparts, gcount_j,
- ci);
+ runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
+ periodic, dim, e, cj->gparts, gcount_j, ci);
}
}
}
diff --git a/tests/testPotentialPair.c b/tests/testPotentialPair.c
index b282606ffe3829398e29c17afa15ade1e2700cf8..b180efe159aec516e2c11948bc4738d266d17d61 100644
--- a/tests/testPotentialPair.c
+++ b/tests/testPotentialPair.c
@@ -42,7 +42,7 @@ const double eps = 0.1;
*/
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)
- error("Values are inconsistent: %12.15e %12.15e (%s)!", a, b, s);
+ error("Values are inconsistent: SWIFT:%12.15e true:%12.15e (%s)!", a, b, s);
}
/* Definitions of the potential and force that match
@@ -92,25 +92,34 @@ int main(int argc, char *argv[]) {
#ifdef HAVE_FE_ENABLE_EXCEPT
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
#endif
-
+
/* Initialise a few things to get us going */
+
+ /* Non-truncated forces first */
+ double rlr = FLT_MAX;
+
struct engine e;
e.max_active_bin = num_time_bins;
e.time = 0.1f;
e.ti_current = 8;
e.time_base = 1e-10;
+ struct space s;
+ s.periodic = 0;
+ e.s = &s;
+
struct pm_mesh mesh;
mesh.periodic = 0;
mesh.dim[0] = 10.;
mesh.dim[1] = 10.;
mesh.dim[2] = 10.;
- mesh.r_s_inv = 0.;
+ mesh.r_s = rlr;
+ mesh.r_s_inv = 1. / rlr;
mesh.r_cut_min = 0.;
+ mesh.r_cut_max = FLT_MAX;
e.mesh = &mesh;
struct gravity_props props;
- props.a_smooth = 1.25;
props.theta_crit2 = 0.;
props.epsilon_cur = eps;
e.gravity_properties = &props;
@@ -119,8 +128,6 @@ int main(int argc, char *argv[]) {
bzero(&r, sizeof(struct runner));
r.e = &e;
- const double rlr = FLT_MAX;
-
/* Init the cache for gravity interaction */
gravity_cache_init(&r.ci_gravity_cache, num_tests);
gravity_cache_init(&r.cj_gravity_cache, num_tests);
@@ -225,8 +232,8 @@ int main(int argc, char *argv[]) {
double acc_true =
acceleration(ci.gparts[0].mass, gp->x[0] - gp2->x[0], 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");
@@ -262,13 +269,13 @@ int main(int argc, char *argv[]) {
const struct gravity_tensors *mpole = ci.multipole;
const double epsilon = gravity_get_softening(gp, &props);
- double pot_true = potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0],
- epsilon, rlr * FLT_MAX);
- double acc_true = acceleration(
- mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr * FLT_MAX);
+ double pot_true =
+ potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr);
+ double acc_true = acceleration(mpole->m_pole.M_000,
+ gp->x[0] - mpole->CoM[0], epsilon, rlr);
- 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);
+ /* 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");
@@ -279,11 +286,47 @@ int main(int argc, char *argv[]) {
/* Reset the accelerations */
for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]);
-/***************************************/
-/* Test the high-order PM interactions */
-/***************************************/
+ /***************************************/
+ /* Test the truncated PM interactions */
+ /***************************************/
+ rlr = 2.;
+ mesh.r_s = rlr;
+ mesh.r_s_inv = 1. / rlr;
+ mesh.periodic = 1;
+ s.periodic = 1;
+ props.epsilon_cur = FLT_MIN; /* No softening */
+
+ /* Now compute the forces */
+ runner_dopair_grav_pp(&r, &ci, &cj, 1);
+
+ /* Verify everything */
+ for (int n = 0; n < num_tests; ++n) {
+ const struct gpart *gp = &cj.gparts[n];
+ const struct gravity_tensors *mpole = ci.multipole;
+ const double epsilon = gravity_get_softening(gp, &props);
+
+ double pot_true =
+ potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr);
+ double acc_true = acceleration(mpole->m_pole.M_000,
+ gp->x[0] - mpole->CoM[0], epsilon, rlr);
+
+ /* 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 truncated P-M interactions all good\n");
+
+ /***************************************/
+ /* Test the high-order PM interactions */
+ /***************************************/
+
+ /* Reset the accelerations */
+ for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]);
-#if SELF_GRAVITY_MULTIPOLE_ORDER >= 3
+#if 0 // SELF_GRAVITY_MULTIPOLE_ORDER >= 3
/* Let's make ci more interesting */
free(ci.gparts);
@@ -346,13 +389,10 @@ int main(int argc, char *argv[]) {
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, epsilon, rlr * FLT_MAX);
- acc_true[0] -=
- acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[0] / d;
- acc_true[1] -=
- acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[1] / d;
- acc_true[2] -=
- acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[2] / d;
+ pot_true += potential(gp2->mass, d, epsilon, rlr);
+ acc_true[0] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[0] / d;
+ acc_true[1] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[1] / d;
+ 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",
diff --git a/tests/testPotentialSelf.c b/tests/testPotentialSelf.c
index 9d57770c862b69009f6e69ba4652843d013af809..6c68825689ea0786d9bfdc8b57c5c5896d1565c8 100644
--- a/tests/testPotentialSelf.c
+++ b/tests/testPotentialSelf.c
@@ -95,7 +95,7 @@ int main(int argc, char *argv[]) {
#ifdef HAVE_FE_ENABLE_EXCEPT
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
#endif
-
+
/* Initialise a few things to get us going */
struct engine e;
e.max_active_bin = num_time_bins;