Small overall updates to the self-gravity scheme tex file.

theory/Multipoles/fmm_standalone.tex

@@ -13,10 +13,10 @@
 
 \maketitle
 
-We use the multi-index notation of \cite{Dehnen2014} to simplify expressions.
-
 \input{potential_softening}
+\input{fmm_summary}
 \input{gravity_derivatives}
+\input{mesh_summary}
 
 \bibliographystyle{mnras}
 \bibliography{./bibliography.bib}

theory/Multipoles/fmm_summary.tex

+\subsection{Evaluating the forces using FMM}
+\label{ssec:fmm_summary}
+
+We use the multi-index notation of \cite{Dehnen2014} to simplify expressions.
+
+

theory/Multipoles/gravity_derivatives.tex

-\subsection{Derivatives of the gravitational potential}
+\subsection{Notes on the derivatives of the gravitational potential}
+\label{ssec:grav_derivatives}
 
 The calculation of all the
 $D_\mathbf{n}(x,y,z) \equiv \nabla^{\mathbf{n}}\phi(x,y,z)$ terms up
@@ -9,7 +10,7 @@ and additions of each of the coordinates and the inverse distance. We
 achieve this by writing $\phi$ as a composition of functions
 $\phi(u(x,y,z))$ and apply the \textit{Fa\`a di Bruno}
 formula \citep[i.e. the ``chain rule'' for higher order derivatives,
-e.g.][]{Hardy2006} to construct our terms:
+ see e.g.][]{Hardy2006} to construct our terms:
 
 \begin{equation}
 \label{eq:faa_di_bruno}
@@ -18,35 +19,39 @@ e.g.][]{Hardy2006} to construct our terms:
 A} \frac{\partial^{|B|}}{\prod_{c\in B}\partial x_c} u(x,y,z),
 \end{equation}
 where $A$ is the set of all partitions of $\lbrace1,\cdots, n\rbrace$,
-$B$ is a block of a partition $A$ and $|\cdot|$ denotes the
+$B$ is a block of a partition in the set $A$ and $|\cdot|$ denotes the
 cardinality of a set. For generic functions $\phi$ and $u$ this
 formula yields an untracktable number of terms; an 8th-order
 derivative will have $4140$ (!)  terms in the sum\footnote{The number
-of terms in the sum is given by the Bell number of the same order}. \\
-We choose to write
+  of terms in the sum is given by the Bell number of the same
+  order.}. \\ For the un-softened gravitational potential, we choose to write
 \begin{align}
    \phi(x,y,z) &= 1 / \sqrt{u(x,y,z)}, \\
    u(x,y,z) &= x^2 + y^2 + z^2.
 \end{align}
 This choice allows to have derivatives of any order of $\phi(u)$ that
-only depend on powers of $u$:
+can be easily expressed and only depend on powers of $u$:
 
 \begin{equation}
-f^{(n)}(u) = \frac{\Gamma(\frac{1}{2})}{\Gamma(\frac{1}{2} -
-n)}\frac{1}{u^{n+\frac{1}{2}}}.
+\phi^{(n)}(u) = (-1)^n\cdot\frac{(2n-1)!!}{2^n}\cdot\frac{1}{u^{n+\frac{1}{2}}},
 \end{equation}
-More importantly, this choice of decomposition allows us to have
-derivatives of $u$ only up to second order in $x$, $y$ or $z$. The
-number of non-zero terms in eq. \ref{eq:faa_di_bruno} is hence
-drastically reduced. For instance, when computing
-$D_{(4,1,3)} \equiv \frac{\partial^8}{\partial x^4 \partial y \partial
-z^3} \phi$, $4100$ of the $4140$ terms will involve at least one
+where $!!$ denotes the semi-factorial. More importantly, this
+choice of decomposition allows us to have non-zero derivatives of $u$
+only up to second order in $x$, $y$ or $z$. The number of non-zero
+terms in eq. \ref{eq:faa_di_bruno} is hence drastically reduced. For
+instance, when computing $D_{(4,1,3)}(\mathbf{r}) \equiv
+\frac{\partial^8}{\partial x^4 \partial y \partial z^3}
+\phi(u(x,y,z))$, $4100$ of the $4140$ terms will involve at least one
 zero-valued derivative (e.g. $\partial^3/\partial x^3$ or
 $\partial^2/\partial x\partial y$) of $u$. Furthermore, among the 40
-remaining terms, many will involve the same derivatives and can be
-grouped together, leaving us with a sum of six products of $x$,$y$ and
-$z$. This is generally the case for most of the $D_\mathbf{n}$'s and
-figuring out which terms are identical in a given set of partitions of
-$\lbrace1,\cdots, n\rbrace$ is an interesting exercise in
-combinatorics left for the reader \citep[see also][]{Hardy2006}.
+remaining terms, many will involve the same combination of derivatives
+of $u$ and can be grouped together, leaving us with a sum of six
+products of $x$,$y$ and $z$. This is generally the case for most of
+the $D_\mathbf{n}$'s and figuring out which terms are identical in a
+given set of partitions of $\lbrace1,\cdots, n\rbrace$ is an
+interesting exercise in combinatorics left for the reader \citep[see
+  also][]{Hardy2006}. We use a \texttt{python} script based on this
+technique to generate the actual C routines used within \swift. Some
+examples of these terms are given in Appendix
+\ref{sec:pot_derivatives}.
 

theory/Multipoles/mesh_summary.tex

+\subsection{Coupling the FMM to a mesh for long-range forces}
+\label{ssec:mesh_summary}

theory/Multipoles/potential.py → theory/Multipoles/plot_potential.py

@@ -141,7 +141,7 @@ 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)
+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)
@@ -166,16 +166,3 @@ ylim(0, 0.95)
 ylabel("$|\\overrightarrow{\\nabla}\\phi(r)|$", labelpad=0)
 
 savefig("potential.pdf")
-
-
-
-
-#Construct potential
-# phi = np.zeros(np.size(r))
-# for i in range(np.size(r)):
-#     if r[i] > 2*epsilon:
-#         phi[i] = 1./ r[i]
-#     elif r[i] > epsilon:
-#         phi[i] = -(1./epsilon) * ((32./3.)*u[i]**2 - (48./3.)*u[i]**3 + (38.4/4.)*u[i]**4 - (32./15.)*u[i]**5 + (2./30.)*u[i]**(-1) - (9/5.))
-#     else:
-#         phi[i] = -(1./epsilon) * ((32./6.)*u[i]**2 - (38.4/4.)*u[i]**4 + (32./5.)*u[i]**4 - (7./5.))

theory/Multipoles/potential_derivatives.tex

-\subsection{Derivatives of the potential}
+\section{Derivatives of the potential}
+\label{sec:pot_derivatives}
 
 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
@@ -14,7 +15,8 @@ D_{000}(\mathbf{r}) = \phi (\mathbf{r},H) =
 \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.
+"chain rule". We show representative examples of the first few
+relevant ones here split by order.
 
 \begin{align}
 D_{100}(\mathbf{r}) = \frac{\partial}{\partial r_x} \phi (\mathbf{r},H) = 
@@ -25,6 +27,8 @@ D_{100}(\mathbf{r}) = \frac{\partial}{\partial r_x} \phi (\mathbf{r},H) =
 \right.\nonumber
 \end{align}
 
+\noindent\rule{6cm}{0.4pt}
+
 \begin{align}
 D_{200}(\mathbf{r}) = \frac{\partial^2}{\partial r_x^2} \phi (\mathbf{r},H) = 
 \left\lbrace\begin{array}{rcl}
@@ -44,6 +48,8 @@ D_{110}(\mathbf{r}) = \frac{\partial^2}{\partial r_x\partial r_y} \phi (\mathbf{
 \right.\nonumber
 \end{align}
 
+\noindent\rule{6cm}{0.4pt}
+
 \begin{align}
 D_{300}(\mathbf{r}) = \frac{\partial^3}{\partial r_x^3} \phi (\mathbf{r},H) = 
 \left\lbrace\begin{array}{rcl}
@@ -73,3 +79,25 @@ D_{111}(\mathbf{r}) = \frac{\partial^3}{\partial r_x\partial r_y\partial r_z} \p
 \end{array}
 \right.\nonumber
 \end{align}
+
+\noindent\rule{6cm}{0.4pt}
+
+\begin{align}
+  D_{400}(\mathbf{r}) &=
+  \nonumber
+\end{align}
+
+\begin{align}
+  D_{310}(\mathbf{r}) &=
+  \nonumber
+\end{align}
+
+\begin{align}
+  D_{220}(\mathbf{r}) &=
+  \nonumber
+\end{align}
+
+\begin{align}
+  D_{211}(\mathbf{r}) &=
+  \nonumber
+\end{align}

theory/Multipoles/potential_softening.tex

 \subsection{Gravitational softening}
+\label{ssec:potential_softening}
 
 To avoid artificial two-body relaxation, the Dirac
 $\delta$-distribution of particles is convolved with a softening
@@ -6,7 +7,10 @@ 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
+is cheaper to compute and has a very similar overall shape. The C2
+kernel has the advantage of being branch-free leading to an expression
+which is faster to evaluate using vector units available on modern
+architectures; it also does not require any division. We set
 $\tilde\delta(\mathbf{x}) = \rho(|\mathbf{x}|) = W(|\mathbf{x}|,
 3\epsilon_{\rm Plummer})$, with $W(r, H)$ given by
 
@@ -41,12 +45,13 @@ 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$.}
+\caption{The density (top), potential (middle) and forces (bottom)
+  generated py a point mass in our softened gravitational scheme.
+  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$.  %(for completeness, we chose
+  %$\epsilon=2$).
+  }
 \label{fig:fmm:softening}
 \end{figure}

theory/Multipoles/run.sh

 #!/bin/bash
-python potential.py
+echo "Generating figures..."
+python plot_potential.py
+echo "Generating PDF..."
 pdflatex -jobname=fmm fmm_standalone.tex
 bibtex fmm.aux
 pdflatex -jobname=fmm fmm_standalone.tex