diff --git a/.gitignore b/.gitignore
index b7d64d4fbd54e018ebfdf9d27902eeaa5cf45e19..9f1af75edb069b8b4ebef81b1557403ff06d480b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -108,6 +108,7 @@ tests/testLogger
tests/benchmarkInteractions
tests/testGravityDerivatives
tests/testPotentialSelf
+tests/testPotentialPair
theory/latex/swift.pdf
theory/SPH/Kernels/kernels.pdf
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index d0003da9a8a311d6f413440f50fb5a00a7b68eeb..4a4d5cf7ad6bc191a07d0f3e544f3835ca9b8f5d 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -136,14 +136,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
* @param f_x (return) The x-component of the acceleration.
* @param f_y (return) The y-component of the acceleration.
* @param f_z (return) The z-component of the acceleration.
+ * @param pot (return) The potential.
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
float r_x, float r_y, float r_z, float r2, float h, float h_inv,
- const struct multipole *m, float *f_x, float *f_y, float *f_z) {
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
- runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv3, m->M_000, f_ij);
-#else
+ const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) {
+
+ //#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
+ 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;
+#if 0 // else
/* Get the inverse distance */
const float r_inv = 1.f / sqrtf(r2);
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 6ab0ad1b5acf0ea68542a358e80d48c6f9e3b269..ebdba8993d32d0c74c08a1de697527dbff5833d3 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -21,6 +21,7 @@
#define SWIFT_RUNNER_DOIACT_GRAV_H
/* Includes. */
+#include "active.h"
#include "cell.h"
#include "gravity.h"
#include "inline.h"
@@ -182,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++) {
@@ -276,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++) {
@@ -380,6 +381,7 @@ static INLINE void runner_dopair_grav_pm(
swift_declare_aligned_ptr(float, a_x, ci_cache->a_x, SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(float, a_y, ci_cache->a_y, SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(float, a_z, ci_cache->a_z, SWIFT_CACHE_ALIGNMENT);
+ swift_declare_aligned_ptr(float, pot, ci_cache->pot, SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(int, active, ci_cache->active,
SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(int, use_mpole, ci_cache->use_mpole,
@@ -415,14 +417,15 @@ static INLINE void runner_dopair_grav_pm(
const float r2 = dx * dx + dy * dy + dz * dz;
/* Interact! */
- float f_x, f_y, f_z;
- runner_iact_grav_pm(dx, dy, dz, r2, h_i, h_inv_i, multi_j, &f_x, &f_y,
- &f_z);
+ float f_x, f_y, f_z, pot_ij;
+ runner_iact_grav_pm(dx, dy, dz, r2, h_i, h_inv_i, multi_j, &f_x, &f_y, &f_z,
+ &pot_ij);
/* Store it back */
a_x[pid] = f_x;
a_y[pid] = f_y;
a_z[pid] = f_z;
+ pot[pid] = pot_ij;
#ifdef SWIFT_DEBUG_CHECKS
/* Update the interaction counter */
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 3a4f451a4a3f69907772ce63d0155fa18ba60444..209a31f5b2d9e4e0947a814539e0aebf18421950 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -26,7 +26,8 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
testAdiabaticIndex testRiemannExact testRiemannTRRS testRiemannHLLC \
testMatrixInversion testThreadpool testDump testLogger testInteractions.sh \
testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
- testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf
+ testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf \
+ testPotentialPair
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
@@ -36,7 +37,7 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testAdiabaticIndex testRiemannExact testRiemannTRRS \
testRiemannHLLC testMatrixInversion testDump testLogger \
testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
- testGravityDerivatives testPotentialSelf
+ testGravityDerivatives testPotentialSelf testPotentialPair
# Rebuild tests when SWIFT is updated.
$(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -102,6 +103,8 @@ testGravityDerivatives_SOURCES = testGravityDerivatives.c
testPotentialSelf_SOURCES = testPotentialSelf.c
+testPotentialPair_SOURCES = testPotentialPair.c
+
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testActivePair.sh \
test27cells.sh test27cellsPerturbed.sh testParser.sh testPeriodicBC.sh \
diff --git a/tests/testPotential.c b/tests/testPotentialPair.c
similarity index 62%
rename from tests/testPotential.c
rename to tests/testPotentialPair.c
index 285cf11ecf8d8a44bda87b2e1196aa3baac605cf..00ca87e094a62d2e40bec17aafc2e08dfdf60144 100644
--- a/tests/testPotential.c
+++ b/tests/testPotentialPair.c
@@ -31,7 +31,7 @@
#include "swift.h"
const int num_tests = 100;
-const double eps = 0.02;
+const double eps = 0.1;
/**
* @brief Check that a and b are consistent (up to some relative error)
@@ -96,6 +96,7 @@ int main() {
e.time_base = 1e-10;
struct space s;
+ s.periodic = 1;
s.width[0] = 0.2;
s.width[1] = 0.2;
s.width[2] = 0.2;
@@ -103,6 +104,8 @@ int main() {
struct gravity_props props;
props.a_smooth = 1.25;
+ props.r_cut_min = 0.;
+ props.theta_crit2 = 0.;
e.gravity_properties = &props;
struct runner r;
@@ -112,46 +115,75 @@ int main() {
const double rlr = props.a_smooth * s.width[0];
/* Init the cache for gravity interaction */
- gravity_cache_init(&r.ci_gravity_cache, num_tests * 2);
+ gravity_cache_init(&r.ci_gravity_cache, num_tests);
+ gravity_cache_init(&r.cj_gravity_cache, num_tests);
/* Let's create one cell with a massive particle and a bunch of test particles
*/
- struct cell c;
- bzero(&c, sizeof(struct cell));
- c.width[0] = 1.;
- c.width[1] = 1.;
- c.width[2] = 1.;
- c.loc[0] = 0.;
- c.loc[1] = 0.;
- c.loc[2] = 0.;
- c.gcount = 1 + num_tests;
- c.ti_old_gpart = 8;
- c.ti_gravity_end_min = 8;
- c.ti_gravity_end_max = 8;
-
- posix_memalign((void **)&c.gparts, gpart_align,
- c.gcount * sizeof(struct gpart));
- bzero(c.gparts, c.gcount * sizeof(struct gpart));
-
+ 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.;
+ ci.loc[0] = 0.;
+ ci.loc[1] = 0.;
+ ci.loc[2] = 0.;
+ ci.gcount = 1;
+ ci.ti_old_gpart = 8;
+ ci.ti_old_multipole = 8;
+ ci.ti_gravity_end_min = 8;
+ ci.ti_gravity_end_max = 8;
+
+ cj.width[0] = 1.;
+ cj.width[1] = 1.;
+ cj.width[2] = 1.;
+ cj.loc[0] = 1.;
+ cj.loc[1] = 0.;
+ cj.loc[2] = 0.;
+ cj.gcount = num_tests;
+ cj.ti_old_gpart = 8;
+ cj.ti_old_multipole = 8;
+ cj.ti_gravity_end_min = 8;
+ cj.ti_gravity_end_max = 8;
+
+ /* Allocate multipoles */
+ ci.multipole = malloc(sizeof(struct gravity_tensors));
+ cj.multipole = malloc(sizeof(struct gravity_tensors));
+ bzero(ci.multipole, sizeof(struct gravity_tensors));
+ bzero(cj.multipole, sizeof(struct gravity_tensors));
+
+ /* 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));
+ bzero(ci.gparts, ci.gcount * sizeof(struct gpart));
+
+ posix_memalign((void **)&cj.gparts, gpart_align, cj.gcount * sizeof(struct gpart));
+ bzero(cj.gparts, cj.gcount * sizeof(struct gpart));
+
/* Create the massive particle */
- c.gparts[0].x[0] = 0.;
- c.gparts[0].x[1] = 0.5;
- c.gparts[0].x[2] = 0.5;
- c.gparts[0].mass = 1.;
- c.gparts[0].epsilon = eps;
- c.gparts[0].time_bin = 1;
- c.gparts[0].type = swift_type_dark_matter;
- c.gparts[0].id_or_neg_offset = 1;
+ 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
- c.gparts[0].ti_drift = 8;
+ ci.gparts[0].ti_drift = 8;
#endif
/* Create the mass-less particles */
- for (int n = 1; n < num_tests + 1; ++n) {
+ for (int n = 0; n < num_tests; ++n) {
- struct gpart *gp = &c.gparts[n];
+ struct gpart *gp = &cj.gparts[n];
- gp->x[0] = n / ((double)num_tests);
+ gp->x[0] = 1. + (n + 1) / ((double)num_tests);
gp->x[1] = 0.5;
gp->x[2] = 0.5;
gp->mass = 0.;
@@ -165,22 +197,26 @@ int main() {
}
/* Now compute the forces */
- runner_doself_grav_pp_truncated(&r, &c);
+ runner_dopair_grav_pp(&r, &ci, &cj);
/* Verify everything */
- for (int n = 1; n < num_tests + 1; ++n) {
- const struct gpart *gp = &c.gparts[n];
+ for (int n = 0; n < num_tests; ++n) {
+ const struct gpart *gp = &cj.gparts[n];
+ const struct gpart *gp2 = &ci.gparts[0];
- double pot_true = potential(c.gparts[0].mass, gp->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(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr);
+ acceleration(ci.gparts[0].mass, gp->x[0] - gp2->x[0], gp->epsilon, rlr);
- // message("x=%e f=%e f_true=%e", gp->x[0], gp->a_grav[0], acc_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");
}
- free(c.gparts);
+ free(ci.multipole);
+ free(cj.multipole);
+ free(ci.gparts);
+ free(cj.gparts);
return 0;
}