From ad54b220ecdfc567e18a981590ce8f77e6bc0fe8 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Tue, 5 Jun 2018 18:28:40 +0200
Subject: [PATCH] Use SWIFT's long-range forces correction instead of Gadget's

---
 src/kernel_long_gravity.h | 22 ++++++++++++++++------
 1 file changed, 16 insertions(+), 6 deletions(-)

diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h
index c88cbada0c..987e20c92d 100644
--- a/src/kernel_long_gravity.h
+++ b/src/kernel_long_gravity.h
@@ -30,7 +30,7 @@
 /* Standard headers */
 #include <math.h>
 
-#define GADGET2_LONG_RANGE_CORRECTION
+//#define GADGET2_LONG_RANGE_CORRECTION
 
 /**
  * @brief Computes the long-range correction term for the potential calculation
@@ -51,7 +51,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); //good_approx_expf(x);
   const float alpha = 1.f / (1.f + exp_x);
 
   /* We want 2 - 2 exp(x) * alpha */
@@ -84,12 +84,22 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval(
   *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);
+  const float x = 2.f * u;
+  const float exp_x = expf(x); //good_approx_expf(x);
+  const float alpha = 1.f / (1.f + exp_x);
 
-  *W = arg * exp_arg * term * term - exp_arg * term + 1.f;
+  /* We want 2*(x*alpha - x*alpha^2 - exp(x)*alpha + 1) */
+  *W = 1.f - alpha;
+  *W = *W * x - exp_x;
+  *W = *W * alpha + 1.f;
   *W *= 2.f;
+
+  /* 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
 }
 
-- 
GitLab