From 33ca157160db89dc2ab3faeb61afa478f8b15cb6 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Fri, 8 Jun 2018 19:39:08 +0200
Subject: [PATCH] Some updates to the truncated M2L documentation.

---
 theory/Multipoles/mesh_summary.tex | 69 +++++++++++++++++-------------
 theory/Multipoles/plot_mesh.py     | 25 ++++++++---
 2 files changed, 58 insertions(+), 36 deletions(-)

diff --git a/theory/Multipoles/mesh_summary.tex b/theory/Multipoles/mesh_summary.tex
index 5e9a0a2fed..afab390bfe 100644
--- a/theory/Multipoles/mesh_summary.tex
+++ b/theory/Multipoles/mesh_summary.tex
@@ -5,44 +5,55 @@ We truncate the potential and forces computed via the FMM using a
 smooth function that drops quickly to zero at some scale $r_s$ set by
 the top-level mesh. Traditionally, implementations have used
 expressions which are cheap to evaluate in Fourier space
-\citep[e.g.][]{Bagla2003,
-  Springel2005}. This, however, implies a large cost for each
-interaction computed within the tree as the real-space truncation
-function won't have a simple analytic form that can be evaluated
-efficiently by computers. Since the FMM scheme involves to not only
-evaluate the forces but higher-order derivatives, a more appropiate
-choice is necessary. We use the sigmoid $S(x) \equiv \frac{e^x}{1 + e^x}$
-as the basis of our truncation function write for the potential:
+\citep[e.g.][]{Bagla2003, Springel2005}. This, however, implies a
+large cost for each interaction computed within the tree as the
+real-space truncation function won't have a simple analytic form that
+can be evaluated efficiently even by modern architectures (typically
+an $\mathrm{erf}()$ function). Since the FMM scheme involves to not
+only evaluate the forces but higher-order derivatives, a more
+appropiate choice is necessary. We use the sigmoid
+$\sigma(w) \equiv \frac{e^w}{1 + e^w}$ as the basis of our truncation
+function and write for the potential:
 
 \begin{align}
-  \varphi_s(r) &= \frac{1}{r} \chi(r, r_s) = \frac{1}{r}\times\left[2 - 2S\left(\frac{2r}{r_s}\right)\right].% \nonumber\\
+  \varphi_s(r) &= \frac{1}{r} \times \chi(r, r_s) = \frac{1}{r}\times\left[2 - 2\sigma\left(\frac{2r}{r_s}\right)\right].% \nonumber\\
   %&= \frac{1}{r}\left[2 - \frac{2e^{\frac{2r}{r_s}}}{1+e^{\frac{2r}{r_s}}}.\right] 
 \end{align}
-This function alongside the trunctation function used in \gadget is
-shown on Fig.~\ref{fig:fmm:potential_short}. This choice of $S(x)$ can
-seem rather cumbersome at first but writing $\alpha(x) \equiv (1+e^x)^{-1}$,
-one can expressed all derivatives of $S(x)$ as simple polynomials in
-$\alpha(x)$, which are easy to evaluate. For instance, in the case of
-the direct force evaluation between two particles, we obtain 
-
-
+This function alongside the trunctation function used in
+\gadget\footnote{For completeness, the \gadget expression reads:\\
+  $\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
+$\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
+force evaluation between two particles, we obtain
 \begin{align}
-  |\mathbf{f}_s(r)| &=
-  \frac{1}{r^2}\times\left[\frac{4r}{r_s}S'\left(\frac{2r}{r_s}\right) -
-    2S\left(\frac{2r}{r_s}\right) + 2\right] \nonumber \\
-  %&=
-  %\frac{1}{r^2}\left[\frac{4r}{r_s}\frac{e^{\frac{2r}{r_s}}}{(1+e^{\frac{2r}{r_s}})^2}
-  %- \frac{2e^{\frac{2r}{r_s}}}{1+e^{\frac{2r}{r_s}}} + 2\right]
+  |\mathbf{f}_s(r)| &= \left|\frac{\partial}{\partial r}\varphi_s(r)\right|
+                      = \left|\frac{\partial}{\partial r}\left(\frac{1}{r} \chi(r, r_s)\right) \right|\nonumber \\
+  &=  \frac{1}{r^2}\times\left[-\frac{4r}{r_s}\sigma'\left(\frac{2r}{r_s}\right) -
+    2\sigma\left(\frac{2r}{r_s}\right) + 2\right] \nonumber \\
   &=
-  \frac{1}{r^2}\times 2 \left[x\alpha(x) - x\alpha(x)^2 - e^x\alpha(x) + 1\right],
+  % \frac{1}{r^2}\times 2 \left[x\alpha(x) - x\alpha(x)^2 - e^x\alpha(x) + 1\right],
+  % \frac{1}{r^2}\times 2 \left[1 - e^x\alpha(x) - xe^x\alpha^2(x)\right],
+    \frac{1}{r^2}\times 2 \left[1 - e^x\left(\alpha(x) - x\alpha(x)^2\right) \right]
 \end{align}
 with $x\equiv2r/r_s$. The truncated force is compared to the Newtonian
-force on Fig.~\ref{fig:fmm:force_short}. At distance $r<r_s/10$, the
+force and to the \gadget truncated forces\footnote{For completeness,
+  the \gadget expression for the norm of the truncated forces is:
+  $|\mathbf{f}_s(r)| = \frac{1}{r^2} \times
+  \left[\mathrm{erfc}\left(\frac{1}{2}\frac{r}{r_s}\right) +
+    \frac{1}{\sqrt{\upi}}\frac{r}{r_s}\exp\left(-\frac{1}{4}\frac{r^2}{r_s^2}\right)\right]$.}
+on Fig.~\ref{fig:fmm:force_short}. At distance $r<r_s/10$, the
 truncation term is negligibly close to one and the truncated forces
 can be replaced by their Newtonian equivalent. We use this
 optimization in \swift and only compute truncated forces between pairs
 of particles that are in tree-leaves larger than $1/10$ of the mesh
-size or between two tree-leaves distant by more than that amount.
+size or between
+two tree-leaves distant by more than that amount.\\
+
+MORE WORDS HERE.\\
 
 The truncation function in Fourier space reads
 
@@ -67,9 +78,9 @@ The truncation function in Fourier space reads
 
 \begin{figure}
 \includegraphics[width=\columnwidth]{force_short.pdf}
-\caption{used in \swift (green line) and \gadget
-  (yellow line) alongside the full Newtonian force term (blue dasheed
-  line). The green dash-dotted line corresponds to the same
+\caption{Norm of the truncated forces used in \swift (green line) and
+  \gadget (yellow line) alongside the full Newtonian force term (blue
+  dasheed line). The green dash-dotted line corresponds to the same
   trunctation function where the exponential in the sigmoid is
   replaced by a sixth order Taylor expansion. At $r<r_s/10$, the
   truncated forces becomes almost equal to the Newtonian ones and can
diff --git a/theory/Multipoles/plot_mesh.py b/theory/Multipoles/plot_mesh.py
index 2468625d14..d5a81c7553 100644
--- a/theory/Multipoles/plot_mesh.py
+++ b/theory/Multipoles/plot_mesh.py
@@ -35,9 +35,9 @@ params = {'axes.labelsize': 9,
 'ytick.labelsize': 8,
 'text.usetex': True,
 'figure.figsize' : (3.15,3.15),
-'figure.subplot.left'    : 0.12,
+'figure.subplot.left'    : 0.14,
 'figure.subplot.right'   : 0.99  ,
-'figure.subplot.bottom'  : 0.09  ,
+'figure.subplot.bottom'  : 0.1  ,
 'figure.subplot.top'     : 0.99  ,
 'figure.subplot.wspace'  : 0.  ,
 'figure.subplot.hspace'  : 0.  ,
@@ -111,6 +111,9 @@ ylim(-1.1, 1.1)
 xlim(-4.1, 4.1)
 savefig("temp.pdf")
 
+def alpha(x):
+    return 1. / (1. + exp(x))
+
 # Correction in real space
 corr_short_gadget2 = special.erf(r / (2.*r_s))
 corr_short_swift = swift_corr(r / (2.*r_s))
@@ -119,6 +122,13 @@ eta_short_gadget2 = special.erfc(r / (2.*r_s)) + (r / (r_s * math.sqrt(math.pi))
 eta_short_swift = 4. * (r / r_s) * d_sigmoid(2. * r / r_s) - 2. * sigmoid(2 * r / r_s) + 2.
 eta_short_swift2 = 4. * (r / r_s) * my_d_sigmoid(2. * r / r_s) - 2. * my_sigmoid(2 * r / r_s) + 2.
 
+#x = 2. * r / r_s
+#force_corr = 2. * (1. - exp(x) * (alpha(x) - x * alpha(x)**2))
+#force_corr = 2. * (1.- x*exp(x)*alpha(x)**2 - exp(x)*alpha(x))
+#force_corr = 2. * (x*alpha(x) - x*alpha(x)**2 -exp(x)*alpha(x) + 1)
+#force_corr = abs(2 * (1. - exp(x) * alpha(x) + x * exp(2*x)*alpha(x)**2 - x*exp(x)*alpha(x)))
+#force_corr = abs(force_corr)
+
 # Corection in Fourier space
 corr_long_gadget2 = exp(-k**2*r_s**2)
 corr_long_swift = math.pi * k * r_s * csch(0.5 * math.pi * r_s * k) / 2.
@@ -155,7 +165,7 @@ xlim(1.1*r_min/ r_s, 0.9*r_max / r_s)
 ylim(1.1/r_max, 0.9/r_min)
 ylabel("$\\varphi_s(r)$", labelpad=-3)
 
-legend(loc="upper right", frameon=True, handletextpad=0.1, handlelength=3.2, fontsize=8)
+legend(loc="upper right", frameon=True, handletextpad=0.3, handlelength=1.6, fontsize=8)
 
 # Correction
 subplot(312, xscale="log", yscale="log")
@@ -174,7 +184,7 @@ ylim(3e-3, 1.5)
 #ylabel("$\\chi_s(r)$", labelpad=-3)
 ylabel("$\\varphi_s(r) \\times r$", labelpad=-2)
 
-legend(loc="center left", frameon=False, handletextpad=0.1, handlelength=2.2, fontsize=7)
+legend(loc="center left", frameon=False, handletextpad=0.3, handlelength=2.2, fontsize=7)
 
 # 1 - Correction
 subplot(313, xscale="log", yscale="log")
@@ -205,6 +215,7 @@ subplot(311, xscale="log", yscale="log")
 plot(r_rs, force_newton, '--', lw=1.4, label="${\\rm Newtonian}$", color=colors[0])
 plot(r_rs, force_short_gadget2, '-', lw=1.4, label="${\\rm Gadget}$", color=colors[2])
 plot(r_rs, force_short_swift, '-', lw=1.4, label="${\\rm SWIFT}$", color=colors[3])
+#plot(r_rs, (1./r**2) * force_corr, '-', lw=1.2, color='r')
 plot(r_rs, force_short_swift2, '-.', lw=1.4, color=colors[3])
 plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
 
@@ -213,7 +224,7 @@ ylim(1.1/r_max**2, 0.9/r_min**2)
 ylabel("$|\\mathbf{f}_s(r)|$", labelpad=-3)
 yticks([1e-4, 1e-2, 1e0, 1e2], ["$10^{-4}$", "$10^{-2}$", "$10^{0}$", "$10^{2}$"])
 
-legend(loc="upper right", frameon=True, handletextpad=0.1, handlelength=3.2, fontsize=8)
+legend(loc="upper right", frameon=True, handletextpad=0.3, handlelength=1.6, fontsize=8)
 
 # Correction
 subplot(312, xscale="log", yscale="log")
@@ -232,7 +243,7 @@ ylim(3e-3, 1.5)
 #ylabel("$\\eta_s(r)$", labelpad=-3)
 ylabel("$|\\mathbf{f}_s(r)|\\times r^2$", labelpad=-2)
 
-legend(loc="center left", frameon=False, handletextpad=0.1, handlelength=2.2, fontsize=7)
+legend(loc="center left", frameon=False, handletextpad=0.3, handlelength=2.2, fontsize=7)
 
 # 1 - Correction
 subplot(313, xscale="log", yscale="log")
@@ -264,7 +275,7 @@ plot(k_rs, phit_long_gadget2, '-', lw=1.4, label="${\\rm Gadget}$", color=colors
 plot(k_rs, phit_long_swift, '-', lw=1.4, label="${\\rm SWIFT}$", color=colors[3])
 plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
 
-legend(loc="lower left", frameon=True, handletextpad=0.1, handlelength=3.2, fontsize=8)
+legend(loc="lower left", frameon=True, handletextpad=0.3, handlelength=1.6, fontsize=8)
 
 xlim(1.1*r_min/ r_s, 0.9*r_max / r_s)
 ylim(1.1/r_max**2, 0.9/r_min**2)
-- 
GitLab