diff --git a/theory/Multipoles/mesh_summary.tex b/theory/Multipoles/mesh_summary.tex
index 5e9a0a2fed2d95474a975aed39c17c117118970f..afab390bfe9d9ae5aa416b470d8ec44d50c8d944 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 2468625d1408556f4f5fb00db17d5a331becdac4..d5a81c75530cc4d8c3bc0bfda20d06d5e0316b3b 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)