Added function to compute long-range FFT correction term

src/Makefile.am

@@ -47,8 +47,8 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
-		 vector.h runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h kick.h \
-		 timestep.h drift.h \
+		 kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h units.h intrinsics.h \
+		 minmax.h kick.h timestep.h drift.h \
 		 gravity.h gravity_io.h \
 		 gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
 		 gravity/Default/gravity_debug.h gravity/Default/gravity_part.h  \

src/gravity/Default/gravity_iact.h

@@ -23,6 +23,7 @@
 /* Includes. */
 #include "const.h"
 #include "kernel_gravity.h"
+#include "kernel_long_gravity.h"
 #include "multipole.h"
 #include "vector.h"
 
@@ -49,7 +50,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
 
   if (r >= hi) {
 
-    /* Get Newtonian graavity */
+    /* Get Newtonian gravity */
     fi = mj * ir * ir * ir;
 
   } else {
@@ -61,7 +62,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
 
   if (r >= hj) {
 
-    /* Get Newtonian graavity */
+    /* Get Newtonian gravity */
     fj = mi * ir * ir * ir;
 
   } else {
@@ -102,7 +103,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
 
   if (r >= hi) {
 
-    /* Get Newtonian graavity */
+    /* Get Newtonian gravity */
     f = mj * ir * ir * ir;
 
   } else {

src/kernel_gravity.h

@@ -63,9 +63,9 @@ static const float
                                     0.f}; /* 1 < u */
 
 /**
- * @brief Computes the kernel function.
+ * @brief Computes the gravity softening function.
  *
- * @param u The ratio of the distance to the smoothing length $u = x/h$.
+ * @param u The ratio of the distance to the softening length $u = x/h$.
  * @param W (return) The value of the kernel function $W(x,h)$.
  */
 __attribute__((always_inline)) INLINE static void kernel_grav_eval(

src/kernel_long_gravity.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_KERNEL_LONG_GRAVITY_H
+#define SWIFT_KERNEL_LONG_GRAVITY_H
+
+#include <math.h>
+
+/* Includes. */
+#include "const.h"
+#include "inline.h"
+#include "vector.h"
+
+#define one_over_sqrt_pi ((float)(M_2_SQRTPI * 0.5))
+
+/**
+ * @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 W (return) The value of the kernel function.
+ */
+__attribute__((always_inline)) INLINE static void kernel_long_grav_eval(
+    float u, float *const W) {
+
+  const float arg1 = u * 0.5f;
+  const float arg2 = u * one_over_sqrt_pi;
+  const float arg3 = -arg1 * arg1;
+
+  const float term1 = erfc(arg1);
+  const float term2 = arg2 * expf(arg3);
+
+  *W = term1 + term2;
+}
+
+#endif  // SWIFT_KERNEL_LONG_GRAVITY_H