diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h index fab633053a468d6e5ba291b2b782d4766a912f14..77bc7c17e2ab8038a4860c74daebd19a2d785070 100644 --- a/src/kernel_long_gravity.h +++ b/src/kernel_long_gravity.h @@ -49,7 +49,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); const float alpha = 1.f / (1.f + exp_x); /* We want 2 - 2 exp(x) * alpha */ @@ -83,7 +83,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval( #else const float arg = 2.f * u; - const float exp_arg = good_approx_expf(arg); + const float exp_arg = expf(arg); const float term = 1.f / (1.f + exp_arg); *W = arg * exp_arg * term * term - exp_arg * term + 1.f; diff --git a/tests/testPotential.c b/tests/testPotential.c index 2bd642c732b67a4096c5f2c0f8eed9f7dbea2280..79f62555b5d67afaac50327e27dbd8c02abe177f 100644 --- a/tests/testPotential.c +++ b/tests/testPotential.c @@ -19,6 +19,7 @@ #include "../config.h" /* Some standard headers. */ +#include <math.h> #include <fenv.h> #include <stdio.h> #include <stdlib.h> @@ -30,25 +31,58 @@ #include "runner_doiact_grav.h" const int num_tests = 100; +const double eps = 0.01; +/* Definitions of the potential and force that match + exactly the theory document */ +double S(double x) { + return exp(x) / (1. + exp(x)); +} + +double potential(double r, double H, double rlr) { + + const double u = r / H; + const double x = r / rlr; + double pot; + if (u > 1.) + pot = -1. / 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; + + return pot * (2. - 2.*S(2.*x)); +} + int main() { /* Initialize CPU frequency, this also starts time. */ unsigned long long cpufreq = 0; clocks_set_cpufreq(cpufreq); - + + /* Initialise a few things to get us going */ struct engine e; e.max_active_bin = num_time_bins; e.time = 0.1f; e.ti_current = 8; - e.timeBase = 1e-10; + e.time_base = 1e-10; + struct space s; + s.width[0] = 0.2; + 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; - gravity_cache_init(&r.ci_gravity_cache, num_tests * 2); + 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 */ struct cell c; @@ -69,7 +103,7 @@ int main() { c.gparts[0].x[1] = 0.5; c.gparts[0].x[2] = 0.5; c.gparts[0].mass = 1.; - c.gparts[0].epsilon = 0.; + 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; @@ -86,19 +120,24 @@ int main() { gp->x[1] = 0.5; gp->x[2] = 0.5; gp->mass = 0.; - gp->epsilon = 0.; + gp->epsilon = eps; gp->time_bin = 1; gp->type = swift_type_dark_matter; gp->id_or_neg_offset = n + 1; #ifdef SWIFT_DEBUG_CHECKS gp->ti_drift = 8; #endif - } /* Now compute the forces */ - runner_doself_grav_pp_full(&r, &c); - + runner_doself_grav_pp_truncated(&r, &c); + + 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)); + } + free(c.gparts); return 0; }