From 8f4f81dc6227166c7e6ce35bbb6e4565d5909ad8 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Mon, 25 Sep 2017 17:17:41 +0100
Subject: [PATCH] Fixed FPE problem due to vectorization of the FFT task.
---
src/kernel_long_gravity.h | 2 +-
src/runner_doiact_fft.c | 4 ++--
2 files changed, 3 insertions(+), 3 deletions(-)
diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h
index ec31c27430..2e0097ded2 100644
--- a/src/kernel_long_gravity.h
+++ b/src/kernel_long_gravity.h
@@ -78,7 +78,7 @@ __attribute__((always_inline)) INLINE static void fourier_kernel_long_grav_eval(
#else
const double u = sqrt(u2);
const double arg = M_PI_2 * u;
- *W = arg / sinh(arg);
+ *W = arg / (sinh(arg) + FLT_MIN);
#endif
}
diff --git a/src/runner_doiact_fft.c b/src/runner_doiact_fft.c
index 26b59f9f6b..a8a882e4e7 100644
--- a/src/runner_doiact_fft.c
+++ b/src/runner_doiact_fft.c
@@ -231,7 +231,7 @@ void runner_do_grav_fft(struct runner* r, int timer) {
const int kz = (k > N_half ? k - N : k);
const double kz_d = (double)kz;
const double fz = k_fac * kz_d;
- const double sinc_kz_inv = (kz != 0) ? fz / sin(fz) : 1.;
+ const double sinc_kz_inv = (kz != 0) ? fz / (sin(fz) + FLT_MIN) : 1.;
/* Norm of vector in Fourier space */
const double k2 = (kx_d * kx_d + ky_d * ky_d + kz_d * kz_d);
@@ -242,7 +242,7 @@ void runner_do_grav_fft(struct runner* r, int timer) {
/* Green function */
double W;
fourier_kernel_long_grav_eval(k2 * a_smooth2, &W);
- const double green_cor = green_fac * W / k2;
+ const double green_cor = green_fac * W / (k2 + FLT_MIN);
/* Deconvolution of CIC */
const double CIC_cor = sinc_kx_inv * sinc_ky_inv * sinc_kz_inv;
--
GitLab