diff --git a/src/approx_math.h b/src/approx_math.h index ad07adeb4f3b1b54ca5f33d80eabb6a004d2a3aa..48319ddfd7a86c132a1cd18b4a08fa849a36a15a 100644 --- a/src/approx_math.h +++ b/src/approx_math.h @@ -36,4 +36,17 @@ __attribute__((always_inline)) INLINE static float approx_expf(float x) { return 1.f + x * (1.f + x * (0.5f + x * (1.f / 6.f + 1.f / 24.f * x))); } +/** + * @brief Approximate version of expf(x) using a 6th order Taylor expansion + * + */ +__attribute__((always_inline)) INLINE static float good_approx_expf(float x) { + return 1.f + + x * (1.f + + x * (0.5f + + x * ((1.f / 6.f) + + x * ((1.f / 24.f) + + x * ((1.f / 120.f) + (1.f / 720.f) * x))))); +} + #endif /* SWIFT_APPROX_MATH_H */ diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h index 6952681999f833bce7755a72aaee742a7fa0ed22..ec31c2743079da22d1f3dd0c8683adf674aca1e3 100644 --- a/src/kernel_long_gravity.h +++ b/src/kernel_long_gravity.h @@ -19,24 +19,30 @@ #ifndef SWIFT_KERNEL_LONG_GRAVITY_H #define SWIFT_KERNEL_LONG_GRAVITY_H -#include <math.h> +/* Config parameters. */ +#include "../config.h" -/* Includes. */ +/* Local headers. */ +#include "approx_math.h" #include "const.h" #include "inline.h" -#include "vector.h" -#define one_over_sqrt_pi ((float)(M_2_SQRTPI * 0.5)) +/* Standard headers */ +#include <math.h> /** * @brief Computes the long-range correction term for the FFT calculation. * - * @param u The ratio of the distance to the FFT cell scale $u = x/A$. + * @param u The ratio of the distance to the FFT cell scale \f$u = r/r_s\f$. * @param W (return) The value of the kernel function. */ __attribute__((always_inline)) INLINE static void kernel_long_grav_eval( float u, float *const W) { +#ifdef GADGET2_LONG_RANGE_CORRECTION + + const float one_over_sqrt_pi = ((float)(M_2_SQRTPI * 0.5)); + const float arg1 = u * 0.5f; const float arg2 = u * one_over_sqrt_pi; const float arg3 = -arg1 * arg1; @@ -45,6 +51,35 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_eval( const float term2 = arg2 * expf(arg3); *W = term1 + term2; +#else + + const float arg = 2.f * u; + const float exp_arg = good_approx_expf(arg); + const float term = 1.f / (1.f + exp_arg); + + *W = arg * exp_arg * term * term - exp_arg * term + 1.f; + *W *= 2.f; +#endif +} + +/** + * @brief Returns the long-range truncation of the Poisson potential in Fourier + * space. + * + * @param u2 The square of the Fourier mode times the cell scale + * \f$u^2 = k^2r_s^2\f$. + * @param W (return) The value of the kernel function. + */ +__attribute__((always_inline)) INLINE static void fourier_kernel_long_grav_eval( + double u2, double *const W) { + +#ifdef GADGET2_LONG_RANGE_CORRECTION + *W = exp(-u2); +#else + const double u = sqrt(u2); + const double arg = M_PI_2 * u; + *W = arg / sinh(arg); +#endif } #endif // SWIFT_KERNEL_LONG_GRAVITY_H diff --git a/src/runner_doiact_fft.c b/src/runner_doiact_fft.c index 266cae61f4c28ba6271b1ef88f6a53bcd129bfc4..26b59f9f6b864445df9190c6041ee684c456ba22 100644 --- a/src/runner_doiact_fft.c +++ b/src/runner_doiact_fft.c @@ -20,9 +20,6 @@ /* Config parameters. */ #include "../config.h" -/* Some standard headers. */ -#include <pthread.h> - #ifdef HAVE_FFTW #include <fftw3.h> #endif @@ -33,6 +30,7 @@ /* Local includes. */ #include "engine.h" #include "error.h" +#include "kernel_long_gravity.h" #include "runner.h" #include "space.h" #include "timers.h" @@ -242,7 +240,9 @@ void runner_do_grav_fft(struct runner* r, int timer) { if (k2 == 0.) continue; /* Green function */ - const double green_cor = green_fac * exp(-k2 * a_smooth2) / k2; + double W; + fourier_kernel_long_grav_eval(k2 * a_smooth2, &W); + const double green_cor = green_fac * W / k2; /* Deconvolution of CIC */ const double CIC_cor = sinc_kx_inv * sinc_ky_inv * sinc_kz_inv;