From c5197e0804cfe32033ff03baa4adb0fba5db924d Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sun, 12 Jun 2016 20:17:22 +0200
Subject: [PATCH] Added function to compute long-range FFT correction term
---
src/Makefile.am | 4 +--
src/gravity/Default/gravity_iact.h | 7 +++--
src/kernel_gravity.h | 4 +--
src/kernel_long_gravity.h | 50 ++++++++++++++++++++++++++++++
4 files changed, 58 insertions(+), 7 deletions(-)
create mode 100644 src/kernel_long_gravity.h
diff --git a/src/Makefile.am b/src/Makefile.am
index cd5a47ea8d..311c49570a 100644
--- a/src/Makefile.am
+++ b/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 \
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index d0624ab2b2..997bff0fb1 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/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 {
diff --git a/src/kernel_gravity.h b/src/kernel_gravity.h
index 1e39f96854..b38feb5758 100644
--- a/src/kernel_gravity.h
+++ b/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(
diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h
new file mode 100644
index 0000000000..d247c7a461
--- /dev/null
+++ b/src/kernel_long_gravity.h
@@ -0,0 +1,50 @@
+/*******************************************************************************
+ * 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
--
GitLab