From c5a703a297c6eb8fed2625be818d8c5296e96da6 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sat, 9 Jun 2018 13:53:45 +0200
Subject: [PATCH] Full expression for the Gadget2 truncation function and its
derivatives.
---
src/kernel_long_gravity.h | 37 +++++++++++++++++++++++++++++-
theory/Multipoles/mesh_summary.tex | 2 +-
2 files changed, 37 insertions(+), 2 deletions(-)
diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h
index f818d67d9f..0adfcb303d 100644
--- a/src/kernel_long_gravity.h
+++ b/src/kernel_long_gravity.h
@@ -31,7 +31,7 @@
#include <float.h>
#include <math.h>
-//#define GADGET2_LONG_RANGE_CORRECTION
+#define GADGET2_LONG_RANGE_CORRECTION
struct truncated_derivatives {
@@ -47,6 +47,40 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_derivatives(
const float r, const float rs_inv,
struct truncated_derivatives *const derivs) {
+#ifdef GADGET2_LONG_RANGE_CORRECTION
+
+ /* Powers of u=r/2r_s */
+ const float u = 0.5f * r * rs_inv;
+ const float u2 = u * u;
+ const float u3 = u2 * u;
+ const float u4 = u3 * u;
+
+ /* Powers of (1/r_s) */
+ const float rs_inv2 = rs_inv * rs_inv;
+ const float rs_inv3 = rs_inv2 * rs_inv;
+ const float rs_inv4 = rs_inv3 * rs_inv;
+ const float rs_inv5 = rs_inv4 * rs_inv;
+
+ /* Derivatives of \chi */
+ derivs->chi_0 = erfcf(u);
+ derivs->chi_1 = -rs_inv;
+ derivs->chi_2 = rs_inv2 * u;
+ derivs->chi_3 = -rs_inv3 * (u2 - 0.5f);
+ derivs->chi_4 = rs_inv4 * (u3 - 1.5f * u);
+ derivs->chi_5 = -rs_inv5 * (u4 - u2 + 0.75f);
+
+ const float one_over_sqrt_pi = ((float)(M_2_SQRTPI * 0.5));
+ const float common_factor = one_over_sqrt_pi * expf(-u2);
+
+ /* Multiply in the common factors */
+ derivs->chi_1 *= common_factor;
+ derivs->chi_2 *= common_factor;
+ derivs->chi_3 *= common_factor;
+ derivs->chi_4 *= common_factor;
+ derivs->chi_5 *= common_factor;
+
+#else
+
/* Powers of 2/r_s */
const float c0 = 1.f;
const float c1 = 2.f * rs_inv;
@@ -80,6 +114,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_derivatives(
derivs->chi_4 = -2.f * exp_x * c4 * (24.f * a5 - 36.f * a4 + 14.f * a3 - a2);
derivs->chi_5 = -2.f * exp_x * c5 *
(120.f * a6 - 240.f * a5 + 150.f * a4 - 30.f * a3 + a2);
+#endif
}
/**
diff --git a/theory/Multipoles/mesh_summary.tex b/theory/Multipoles/mesh_summary.tex
index afab390bfe..1b952f0dba 100644
--- a/theory/Multipoles/mesh_summary.tex
+++ b/theory/Multipoles/mesh_summary.tex
@@ -24,7 +24,7 @@ This function alongside the trunctation function used in
$\chi(r, r_s) = \mathrm{erfc}(\frac{1}{2}\frac{r}{r_s})$.} is shown
on Fig.~\ref{fig:fmm:potential_short}. This choice of $\sigma(w)$ can
seem rather cumbersome at first but writing
-$\alpha(w) \equiv (1+e^w)^{-1}$, one can expressed all derivatives of
+$\alpha(w) \equiv (1+e^w)^{-1}$, one can express all derivatives of
$\sigma(w)$ as simple polynomials in $\alpha(w)$ (with an identical
$e^w$ pre-factor), which are easy and cheap to evaluate (see Appendix
\ref{sec:pot_derivatives}). For instance, in the case of the direct
--
GitLab