Skip to content
Snippets Groups Projects
Commit 3f7769ae authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Revert to using the approximate exponential in the gravity long-range truncation.

parent 07f1e5b1
Branches
Tags
No related merge requests found
...@@ -49,4 +49,16 @@ __attribute__((always_inline)) INLINE static float good_approx_expf(float x) { ...@@ -49,4 +49,16 @@ __attribute__((always_inline)) INLINE static float good_approx_expf(float x) {
x * ((1.f / 120.f) + (1.f / 720.f) * x))))); x * ((1.f / 120.f) + (1.f / 720.f) * x)))));
} }
/**
* @brief Approximate version of exp(x) using a 6th order Taylor expansion
*/
__attribute__((always_inline)) INLINE static double good_approx_exp(double x) {
return 1. +
x * (1. +
x * (0.5 +
x * ((1. / 6.) +
x * ((1. / 24.) +
x * ((1. / 120.) + (1. / 720.) * x)))));
}
#endif /* SWIFT_APPROX_MATH_H */ #endif /* SWIFT_APPROX_MATH_H */
...@@ -49,7 +49,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_pot_eval( ...@@ -49,7 +49,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_pot_eval(
#else #else
const float x = 2.f * u; const float x = 2.f * u;
const float exp_x = expf(x); const float exp_x = good_approx_expf(x);
const float alpha = 1.f / (1.f + exp_x); const float alpha = 1.f / (1.f + exp_x);
/* We want 2 - 2 exp(x) * alpha */ /* We want 2 - 2 exp(x) * alpha */
...@@ -83,7 +83,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval( ...@@ -83,7 +83,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval(
#else #else
const float arg = 2.f * u; const float arg = 2.f * u;
const float exp_arg = expf(arg); const float exp_arg = good_approx_expf(arg);
const float term = 1.f / (1.f + exp_arg); const float term = 1.f / (1.f + exp_arg);
*W = arg * exp_arg * term * term - exp_arg * term + 1.f; *W = arg * exp_arg * term * term - exp_arg * term + 1.f;
...@@ -111,4 +111,4 @@ __attribute__((always_inline)) INLINE static void fourier_kernel_long_grav_eval( ...@@ -111,4 +111,4 @@ __attribute__((always_inline)) INLINE static void fourier_kernel_long_grav_eval(
#endif #endif
} }
#endif // SWIFT_KERNEL_LONG_GRAVITY_H #endif /* SWIFT_KERNEL_LONG_GRAVITY_H */
...@@ -41,15 +41,18 @@ const double eps = 0.02; ...@@ -41,15 +41,18 @@ const double eps = 0.02;
* @param s String used to identify this check in messages * @param s String used to identify this check in messages
*/ */
void check_value(double a, double b, const char *s) { 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) if (fabs(a - b) / fabs(a + b) > 1e-6 && fabs(a - b) > 1.e-6)
error("Values are inconsistent: %12.15e %12.15e (%s)!", a, b, s); error("Values are inconsistent: %12.15e %12.15e (%s)!", a, b, s);
} }
/* Definitions of the potential and force that match /* Definitions of the potential and force that match
exactly the theory document */ exactly the theory document */
double S(double x) { return exp(x) / (1. + exp(x)); } double S(double x) { return good_approx_exp(x) / (1. + good_approx_exp(x)); }
double S_prime(double x) { return exp(x) / ((1. + exp(x)) * (1. + exp(x))); } double S_prime(double x) {
return good_approx_exp(x) /
((1. + good_approx_exp(x)) * (1. + good_approx_exp(x)));
}
double potential(double mass, double r, double H, double rlr) { double potential(double mass, double r, double H, double rlr) {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment