Commit 1f81ee19 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Test the potential calculation also with the long-range truncation.

parent bab1068c
......@@ -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;
......
......@@ -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;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment