Added more documentation of the softened FMM scheme

theory/Multipoles/bibliography.bib

@@ -24,4 +24,75 @@ archivePrefix = "arXiv",
    pages={13},
    year={2006}
  }
-	      
+
+@ARTICLE{Monaghan1985,
+   author = {{Monaghan}, J.~J. and {Lattanzio}, J.~C.},
+    title = "{A refined particle method for astrophysical problems}",
+  journal = {\aap},
+ keywords = {Computational Astrophysics, Gravitational Collapse, Gravitational Fields, Many Body Problem, Molecular Clouds, Stellar Evolution, Angular Momentum, Binary Stars, Computational Grids, Interpolation, Kernel Functions, Particle Mass, Stellar Orbits},
+     year = 1985,
+    month = aug,
+   volume = 149,
+    pages = {135-143},
+   adsurl = {http://adsabs.harvard.edu/abs/1985A%26A...149..135M},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@ARTICLE{Price2007,
+   author = {{Price}, D.~J. and {Monaghan}, J.~J.},
+    title = "{An energy-conserving formalism for adaptive gravitational force softening in smoothed particle hydrodynamics and N-body codes}",
+  journal = {\mnras},
+   eprint = {astro-ph/0610872},
+ keywords = {gravitation, hydrodynamics, methods: N-body simulations, methods: numerical},
+     year = 2007,
+    month = feb,
+   volume = 374,
+    pages = {1347-1358},
+      doi = {10.1111/j.1365-2966.2006.11241.x},
+   adsurl = {http://adsabs.harvard.edu/abs/2007MNRAS.374.1347P},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@ARTICLE{Dehnen2001,
+   author = {{Dehnen}, W.},
+    title = "{Towards optimal softening in three-dimensional N-body codes - I. Minimizing the force error}",
+  journal = {\mnras},
+   eprint = {astro-ph/0011568},
+ keywords = {STELLAR DYNAMICS, METHODS: N-BODY SIMULATIONS, METHODS: NUMERICAL},
+     year = 2001,
+    month = jun,
+   volume = 324,
+    pages = {273-291},
+      doi = {10.1046/j.1365-8711.2001.04237.x},
+   adsurl = {http://adsabs.harvard.edu/abs/2001MNRAS.324..273D},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@ARTICLE{Hernquist1989,
+   author = {{Hernquist}, L. and {Katz}, N.},
+    title = "{TREESPH - A unification of SPH with the hierarchical tree method}",
+  journal = {\apjs},
+ keywords = {Computational Fluid Dynamics, Computerized Simulation, Data Smoothing, Magnetohydrodynamics, Trees (Mathematics), Dynamical Systems, Many Body Problem, Monte Carlo Method, Spatial Resolution},
+     year = 1989,
+    month = jun,
+   volume = 70,
+    pages = {419-446},
+      doi = {10.1086/191344},
+   adsurl = {http://adsabs.harvard.edu/abs/1989ApJS...70..419H},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@Article{Wendland1995,
+author="Wendland, Holger",
+title="Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree",
+journal="Advances in Computational Mathematics",
+year="1995",
+volume="4",
+number="1",
+pages="389--396",
+abstract="We construct a new class of positive definite and compactly supported radial functions which consist of a univariate polynomial within their support. For given smoothness and space dimension it is proved that they are of minimal degree and unique up to a constant factor. Finally, we establish connections between already known functions of this kind.",
+issn="1572-9044",
+doi="10.1007/BF02123482",
+url="http://dx.doi.org/10.1007/BF02123482"
+}
+

theory/Multipoles/fmm_standalone.tex

-\documentclass[fleqn, usenatbib, useAMS,a4paper]{mnras}
+\documentclass[fleqn, usenatbib, useAMS, a4paper]{mnras}
 \usepackage{graphicx}
 \usepackage{amsmath,paralist,xcolor,xspace,amssymb}
 \usepackage{times}
@@ -13,10 +13,18 @@
 
 \maketitle
 
+We use the multi-index notation of \cite{Dehnen2014} to simplify expressions.
+
+\input{potential_softening}
 \input{gravity_derivatives}
 
 \bibliographystyle{mnras}
 \bibliography{./bibliography.bib}
 
+\appendix
+\onecolumn
+\input{potential_derivatives}
+
+
 
 \end{document}

theory/Multipoles/gravity_derivatives.tex

-We use the multi-index notation of \cite{Dehnen2014} to simplify expressions.
-
 \subsection{Derivatives of the gravitational potential}
 
 The calculation of all the

theory/Multipoles/potential.py

@@ -20,7 +20,6 @@ import matplotlib
 matplotlib.use("Agg")
 from pylab import *
 from scipy import integrate
-import distinct_colours as colours
 from scipy.optimize import curve_fit
 from scipy.optimize import fsolve
 from matplotlib.font_manager import FontProperties
@@ -29,15 +28,15 @@ import math
 
 params = {'axes.labelsize': 9,
 'axes.titlesize': 10,
-'font.size': 12,
+'font.size': 10,
 'legend.fontsize': 10,
-'xtick.labelsize': 9,
-'ytick.labelsize': 9,
+'xtick.labelsize': 8,
+'ytick.labelsize': 8,
 'text.usetex': True,
 'figure.figsize' : (3.15,3.15),
 'figure.subplot.left'    : 0.115,
 'figure.subplot.right'   : 0.99  ,
-'figure.subplot.bottom'  : 0.08  ,
+'figure.subplot.bottom'  : 0.065  ,
 'figure.subplot.top'     : 0.99  ,
 'figure.subplot.wspace'  : 0.  ,
 'figure.subplot.hspace'  : 0.  ,
@@ -112,52 +111,59 @@ for i in range(np.size(r)):
         F_gadget2[i] = u[i] * (10.666667 + u[i]**2 * (32. * u[i] - 38.4)) / epsilon_gadget**2
 
 figure()
+colors=['#4477AA', '#CC6677', '#DDCC77', '#117733']
 
-# Potential
+# Density
 subplot(311)
-plot(r, phi_newton, 'g--', lw=0.8, label="${\\rm Newtonian}$")
-plot(r, phi_plummer, 'b:', lw=0.8, label="${\\rm Plummer}$")
-plot(r, phi_gadget2, 'm-.', lw=0.8, label="${\\rm Gadget-2}$")
-plot(r, phi, 'r-', lw=0.8, label="${\\rm SWIFT}$")
-plot([epsilon, epsilon], [-10, 10], 'k-', alpha=0.5, lw=0.5)
-#text(epsilon+0.05, 2., "$\\epsilon$", color='k', alpha=0.5, rotation=90, va="top", ha="left", fontsize=8) 
-text(1.2, 1.6, "$\\epsilon_{\\rm{Plummer}} = \\frac{\\epsilon}{%d}$"%plummer_equivalent_factor, color='k', fontsize=8, backgroundcolor='w')
-ylim(0, 2.1)
-ylabel("$|\\phi(r)|$", labelpad=0)
+plot(r, W_newton, '--', lw=1.4, label="${\\rm Newtonian}$", color=colors[0])
+plot(r, W_plummer, ':', lw=1.4, label="${\\rm Plummer}$", color=colors[1])
+plot(r, W_gadget2, '-', lw=1.4, label="${\\rm Gadget}$", color=colors[2])
+plot(r, W, '-', lw=1.4, label="${\\rm SWIFT}$", color=colors[3])
+plot([epsilon, epsilon], [0, 10], 'k-', alpha=0.5, lw=0.5)
+plot([epsilon/plummer_equivalent_factor, epsilon/plummer_equivalent_factor], [0, 10], 'k-', alpha=0.5, lw=0.5)
 
+legend(loc="upper right", frameon=True, handletextpad=0.1, handlelength=3.2, fontsize=8)
 
 xlim(0,r_max_plot)
 xticks([0., 0.5, 1., 1.5, 2., 2.5], ["", "", "", "", "", ""])
 
-# Force
+ylim(0., 0.84)
+yticks([0, 0.2, 0.4, 0.6, 0.8], ["$0$", "$0.2$", "$0.4$", "$0.6$", "$0.8$"])
+ylabel("$\\rho(r)$", labelpad=2)
+
+# Potential
 subplot(312)
-plot(r, F_newton, 'g--', lw=0.8)
-plot(r, F_plummer, 'b:', lw=0.8)
-plot(r, F_gadget2, 'm-.', lw=0.8)
-plot(r, F, 'r-', lw=0.8)
-plot([epsilon, epsilon], [0, 10], 'k-', alpha=0.5, lw=0.5)
+plot(r, phi_newton, '--', lw=1.4, label="${\\rm Newtonian}$", color=colors[0])
+plot(r, phi_plummer, ':', lw=1.4, label="${\\rm Plummer}$", color=colors[1])
+plot(r, phi_gadget2, '-', lw=1.4, label="${\\rm Gadget}$", color=colors[2])
+plot(r, phi, '-', lw=1.4, label="${\\rm SWIFT}$", color=colors[3])
+plot([epsilon, epsilon], [-10, 10], 'k-', alpha=0.5, lw=0.5)
+plot([epsilon/plummer_equivalent_factor, epsilon/plummer_equivalent_factor], [0, 10], 'k-', alpha=0.5, lw=0.5)
+
+ylim(0, 2.3)
+ylabel("$|\\phi(r)|$", labelpad=1)
+#yticks([0., 0.5, 1., 1.5, 2., 2.5], ["$%.1f$"%(0.*epsilon), "$%.1f$"%(0.5*epsilon), "$%.1f$"%(1.*epsilon), "$%.1f$"%(1.5*epsilon), "$%.1f$"%(2.*epsilon)])
 
 xlim(0,r_max_plot)
 xticks([0., 0.5, 1., 1.5, 2., 2.5], ["", "", "", "", "", ""])
 
-ylim(0, 0.95)
-ylabel("$|\\overrightarrow{\\nabla}\\phi(r)|$", labelpad=0)
-
-# Density
+# Force
 subplot(313)
-plot(r, W_newton, 'g--', lw=0.8, label="${\\rm Newtonian}$")
-plot(r, W_plummer, 'b:', lw=0.8, label="${\\rm Plummer}$")
-plot(r, W_gadget2, 'm-.', lw=0.8, label="${\\rm Gadget\\textendash 2}$")
-plot(r, W, 'r-', lw=0.8, label="${\\rm SWIFT}$")
+plot(r, F_newton, '--', lw=1.4, color=colors[0])
+plot(r, F_plummer, ':', lw=1.4, color=colors[1])
+plot(r, F_gadget2, '-', lw=1.4, color=colors[2])
+plot(r, F, '-', lw=1.4, color=colors[3])
 plot([epsilon, epsilon], [0, 10], 'k-', alpha=0.5, lw=0.5)
+plot([epsilon/plummer_equivalent_factor, epsilon/plummer_equivalent_factor], [0, 10], 'k-', alpha=0.5, lw=0.5)
+text(epsilon+0.03, 0.05, "$\\epsilon$", color='k', alpha=0.5, rotation=90, va="bottom", ha="left", fontsize=8)
+text(epsilon/plummer_equivalent_factor+0.03, 0.05, "$\\epsilon_{\\rm Plummer}$", color='k', alpha=0.5, rotation=90, va="bottom", ha="left", fontsize=8) 
+
 xlim(0,r_max_plot)
 xticks([0., 0.5, 1., 1.5, 2., 2.5], ["$%.1f$"%(0./epsilon), "", "$%.1f$"%(1./epsilon), "", "$%.1f$"%(2./epsilon)])
-xlabel("$r/\\epsilon$", labelpad=-4)
-legend(loc="upper right", frameon=True, handletextpad=0.1, handlelength=3.2, fontsize=8)
+xlabel("$r/H$", labelpad=-7)
 
-ylim(0., 0.92)
-yticks([0, 0.2, 0.4, 0.6, 0.8], ["$0$", "$0.2$", "$0.4$", "$0.6$", "$0.8$"])
-ylabel("$\\rho(r)$", labelpad=0)
+ylim(0, 0.95)
+ylabel("$|\\overrightarrow{\\nabla}\\phi(r)|$", labelpad=0)
 
 savefig("potential.pdf")
 

theory/Multipoles/potential_derivatives.tex

+\subsection{Derivatives of the potential}
+
+For completeness, we give here the full expression for the first few
+derivatives of the potential that are used in our FMM scheme. We use
+the notation $\mathbf{r}=(r_x, r_y, r_z)$, $r = |\mathbf{r}|$ and
+$u=r/H$. Starting from the potential (Eq. \ref{eq:fmm:potential},
+reproduced here for clarity), 
+\begin{align}
+D_{000}(\mathbf{r}) = \phi (\mathbf{r},H) = 
+\left\lbrace\begin{array}{rcl}
+\frac{1}{H} \left(-3u^7 + 15u^6 - 28u^5 + 21u^4 - 7u^2 + 3\right) & \mbox{if} & u < 1,\\
+\frac{1}{r} & \mbox{if} & u \geq 1, 
+\end{array}
+\right.\nonumber
+\end{align}
+we can construct the higher order terms by successively applying the
+"chain rule". We show examples of the first few relevant ones here.
+
+\begin{align}
+D_{100}(\mathbf{r}) = \frac{\partial}{\partial r_x} \phi (\mathbf{r},H) = 
+\left\lbrace\begin{array}{rcl}
+\frac{r_x}{H^3} \left(-21u^5 + 90u^4 - 140u^3 + 84u^2 - 14\right) & \mbox{if} & u < 1,\\
+-\frac{r_x}{r^3} & \mbox{if} & u \geq 1, 
+\end{array}
+\right.\nonumber
+\end{align}
+
+\begin{align}
+D_{200}(\mathbf{r}) = \frac{\partial^2}{\partial r_x^2} \phi (\mathbf{r},H) = 
+\left\lbrace\begin{array}{rcl}
+\frac{r_x^2}{H^5}\left(-105u^3+360u^2-420u+168\right) +
+\frac{1}{H^3} \left(-21u^5 + 90u^4 - 140u^3 + 84u^2 - 14\right) & \mbox{if} & u < 1,\\
+3\frac{r_x^2}{r^5} - \frac{1}{r^3} & \mbox{if} & u \geq 1, 
+\end{array}
+\right.\nonumber
+\end{align}
+
+\begin{align}
+D_{110}(\mathbf{r}) = \frac{\partial^2}{\partial r_x\partial r_y} \phi (\mathbf{r},H) = 
+\left\lbrace\begin{array}{rcl}
+\frac{r_xr_y}{H^5}\left(-105u^3+360u^2-420u+168\right) & \mbox{if} & u < 1,\\
+3\frac{r_xr_y}{r^5} & \mbox{if} & u \geq 1, 
+\end{array}
+\right.\nonumber
+\end{align}

theory/Multipoles/potential_softening.tex

+\subsection{Gravitational softening}
+
+To avoid artificial two-body relaxation, the Dirac
+$\delta$-distribution of particles is convolved with a softening
+kernel of a given fixed, but time-variable, scale-length
+$\epsilon$. Instead of the commonly used spline kernel of
+\cite{Monaghan1985} (e.g. in \textsc{Gadget}), we use a C2 kernel
+\citep{Wendland1995} which leads to an expression for the force that
+is cheaper to compute and has a very similar overall shape. We set
+$\tilde\delta(\mathbf{x}) = \rho(|\mathbf{x}|) = W(|\mathbf{x}|,
+3\epsilon_{\rm Plummer})$, with $W(r, H)$ given by
+
+\begin{align}
+W(r,H) &= \frac{21}{2\pi H^3} \times \nonumber \\
+&\left\lbrace\begin{array}{rcl}
+4u^5 - 15u^4 + 20u^3 - 10u^2 + 1 & \mbox{if} & u < 1,\\
+0 & \mbox{if} & u \geq 1,
+\end{array}
+\right.
+\end{align}
+and $u = r/H$. The potential $\phi(r,H)$ corresponding to this density distribution reads
+\begin{align}
+\phi = 
+\left\lbrace\begin{array}{rcl}
+\frac{1}{H} (-3u^7 + 15u^6 - 28u^5 + 21u^4 - 7u^2 + 3) & \mbox{if} & u < 1,\\
+\frac{1}{r} & \mbox{if} & u \geq 1.
+\end{array}
+\right.
+\label{eq:fmm:potential}
+\end{align}
+
+These choices, lead to a potential at $|\mathbf{x}| = 0$ equal to the
+central potential of a Plummer sphere (i.e. $1/\epsilon_{\rm
+Plummer}$)\footnote{Note the factor $3$ in the definition of
+$\rho(|\mathbf{x}|)$ which differs from the factor $2.8$ used
+in \textsc{Gadget} as a consequence of the change of kernel
+shape.}. The softened density profile, its corresponding potential and
+resulting forces are shown on Fig. \ref{fig:fmm:softening} (for
+details see Sec. 2 of~\cite{Price2007}).
+
+
+\begin{figure}
+\includegraphics[width=\columnwidth]{potential.pdf}
+\caption{The density (top), potential (middle) and forces (bottom) of
+generated py a point mass in our softened gravitational scheme (for
+completeness, we chose $\epsilon=2$). A
+Plummer-equivalent sphere is shown for comparison. The spline kernel
+of \citet{Monaghan1985}, used in \textsc{Gadget}, is shown for
+comparison but note that it has not been re-scaled to match the
+Plummer-sphere potential at $r=0$.}
+\label{fig:fmm:softening}
+\end{figure}

theory/Multipoles/run.sh

 #!/bin/bash
-python kernels.py
+python potential.py
 pdflatex -jobname=fmm fmm_standalone.tex
 bibtex fmm.aux
 pdflatex -jobname=fmm fmm_standalone.tex