diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index db9c843506b6ec4c7417dbbb07e7414a0f2fe1ab..6ab0ad1b5acf0ea68542a358e80d48c6f9e3b269 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -35,7 +35,8 @@
* @param c The #cell we are working on.
* @param timer Are we timing this ?
*/
-static INLINE void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
+static INLINE void runner_do_grav_down(struct runner *r, struct cell *c,
+ int timer) {
/* Some constants */
const struct engine *e = r->e;
@@ -126,8 +127,9 @@ static INLINE void runner_do_grav_down(struct runner *r, struct cell *c, int tim
* @param ci The #cell with field tensor to interact.
* @param cj The #cell with the multipole.
*/
-static INLINE void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
- struct cell *restrict cj) {
+static INLINE void runner_dopair_grav_mm(const struct runner *r,
+ struct cell *restrict ci,
+ struct cell *restrict cj) {
/* Some constants */
const struct engine *e = r->e;
@@ -440,7 +442,8 @@ static INLINE void runner_dopair_grav_pm(
* @param ci The first #cell.
* @param cj The other #cell.
*/
-static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
+static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
+ struct cell *cj) {
const struct engine *e = r->e;
@@ -641,7 +644,8 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci, stru
*
* @todo Use a local cache for the particles.
*/
-static INLINE void runner_doself_grav_pp_full(struct runner *r, struct cell *c) {
+static INLINE void runner_doself_grav_pp_full(struct runner *r,
+ struct cell *c) {
/* Some constants */
const struct engine *const e = r->e;
@@ -763,7 +767,8 @@ static INLINE void runner_doself_grav_pp_full(struct runner *r, struct cell *c)
*
* @todo Use a local cache for the particles.
*/
-static INLINE void runner_doself_grav_pp_truncated(struct runner *r, struct cell *c) {
+static INLINE void runner_doself_grav_pp_truncated(struct runner *r,
+ struct cell *c) {
/* Some constants */
const struct engine *const e = r->e;
@@ -945,8 +950,8 @@ static INLINE void runner_doself_grav_pp(struct runner *r, struct cell *c) {
*
* @todo Use a local cache for the particles.
*/
-static INLINE void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
- int gettimer) {
+static INLINE void runner_dopair_grav(struct runner *r, struct cell *ci,
+ struct cell *cj, int gettimer) {
/* Some constants */
const struct engine *e = r->e;
@@ -1097,7 +1102,8 @@ static INLINE void runner_dopair_grav(struct runner *r, struct cell *ci, struct
*
* @todo Use a local cache for the particles.
*/
-static INLINE void runner_doself_grav(struct runner *r, struct cell *c, int gettimer) {
+static INLINE void runner_doself_grav(struct runner *r, struct cell *c,
+ int gettimer) {
/* Some constants */
const struct engine *e = r->e;
@@ -1148,7 +1154,8 @@ static INLINE void runner_doself_grav(struct runner *r, struct cell *c, int gett
* @param ci The #cell of interest.
* @param timer Are we timing this ?
*/
-static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
+static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
+ int timer) {
/* Some constants */
const struct engine *e = r->e;
diff --git a/tests/testPotential.c b/tests/testPotential.c
index 79f62555b5d67afaac50327e27dbd8c02abe177f..285cf11ecf8d8a44bda87b2e1196aa3baac605cf 100644
--- a/tests/testPotential.c
+++ b/tests/testPotential.c
@@ -19,39 +19,69 @@
#include "../config.h"
/* Some standard headers. */
-#include <math.h>
#include <fenv.h>
+#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
/* Local headers. */
-#include "swift.h"
#include "runner_doiact_grav.h"
+#include "swift.h"
const int num_tests = 100;
-const double eps = 0.01;
+const double eps = 0.02;
+
+/**
+ * @brief Check that a and b are consistent (up to some relative error)
+ *
+ * @param a First value
+ * @param b Second value
+ * @param s String used to identify this check in messages
+ */
+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);
+}
/* Definitions of the potential and force that match
exactly the theory document */
-double S(double x) {
- return exp(x) / (1. + exp(x));
-}
+double S(double x) { return exp(x) / (1. + exp(x)); }
-double potential(double r, double H, double rlr) {
+double S_prime(double x) { return exp(x) / ((1. + exp(x)) * (1. + exp(x))); }
+
+double potential(double mass, double r, double H, double rlr) {
const double u = r / H;
const double x = r / rlr;
double pot;
if (u > 1.)
- pot = -1. / r;
+ pot = -mass / r;
else
- pot = -(-3.*u*u*u*u*u*u*u + 15.*u*u*u*u*u*u - 28.*u*u*u*u*u + 21.*u*u*u*u - 7.*u*u + 3.)/ H;
+ pot = -mass *
+ (-3. * u * u * u * u * u * u * u + 15. * u * u * u * u * u * u -
+ 28. * u * u * u * u * u + 21. * u * u * u * u - 7. * u * u + 3.) /
+ H;
- return pot * (2. - 2.*S(2.*x));
+ return pot * (2. - 2. * S(2. * x));
}
-
+
+double acceleration(double mass, double r, double H, double rlr) {
+
+ const double u = r / H;
+ const double x = r / rlr;
+ double acc;
+ if (u > 1.)
+ acc = -mass / (r * r * r);
+ else
+ acc = -mass * (21. * u * u * u * u * u - 90. * u * u * u * u +
+ 140. * u * u * u - 84. * u * u + 14.) /
+ (H * H * H);
+
+ return r * acc * (4. * x * S_prime(2 * x) - 2. * S(2. * x) + 2.);
+}
+
int main() {
/* Initialize CPU frequency, this also starts time. */
@@ -70,32 +100,37 @@ int main() {
s.width[1] = 0.2;
s.width[2] = 0.2;
e.s = &s;
-
+
struct gravity_props props;
props.a_smooth = 1.25;
e.gravity_properties = &props;
-
+
struct runner r;
bzero(&r, sizeof(struct runner));
r.e = &e;
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);
-
- /* Let's create one cell with a massive particle and a bunch of test particles */
+
+ /* 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));
+
+ posix_memalign((void **)&c.gparts, gpart_align,
+ c.gcount * sizeof(struct gpart));
bzero(c.gparts, c.gcount * sizeof(struct gpart));
/* Create the massive particle */
@@ -110,13 +145,13 @@ int main() {
#ifdef SWIFT_DEBUG_CHECKS
c.gparts[0].ti_drift = 8;
#endif
-
+
/* Create the mass-less particles */
for (int n = 1; n < num_tests + 1; ++n) {
struct gpart *gp = &c.gparts[n];
- gp->x[0] = n / ((double) num_tests);
+ gp->x[0] = n / ((double)num_tests);
gp->x[1] = 0.5;
gp->x[2] = 0.5;
gp->mass = 0.;
@@ -128,16 +163,24 @@ int main() {
gp->ti_drift = 8;
#endif
}
-
+
/* Now compute the forces */
runner_doself_grav_pp_truncated(&r, &c);
+ /* Verify everything */
for (int n = 1; n < num_tests + 1; ++n) {
const struct gpart *gp = &c.gparts[n];
- message("x=%f pot=%f true=%f",
- gp->x[0], gp->potential, potential(gp->x[0], gp->epsilon, rlr));
+
+ double pot_true = potential(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr);
+ double acc_true =
+ acceleration(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr);
+
+ // message("x=%e f=%e f_true=%e", gp->x[0], gp->a_grav[0], acc_true);
+
+ check_value(gp->potential, pot_true, "potential");
+ check_value(gp->a_grav[0], acc_true, "acceleration");
}
-
+
free(c.gparts);
return 0;
}