diff --git a/theory/Multipoles/bibliography.bib b/theory/Multipoles/bibliography.bib index abdb5855f52aefd9b4562851946d16bd402db68a..12e274dd63093ba1e14750249f2538c092e5268a 100644 --- a/theory/Multipoles/bibliography.bib +++ b/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" +} + diff --git a/theory/Multipoles/fmm_standalone.tex b/theory/Multipoles/fmm_standalone.tex index 6948708622dbf0cfc05daa9d63917f7a3017bfc2..fcd727a89abe95bba69b23c58ce5067c8cc53211 100644 --- a/theory/Multipoles/fmm_standalone.tex +++ b/theory/Multipoles/fmm_standalone.tex @@ -1,4 +1,4 @@ -\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} diff --git a/theory/Multipoles/gravity_derivatives.tex b/theory/Multipoles/gravity_derivatives.tex index 913bdf752c666aac5913c6e6452bdb06f3fcdc86..e4c7b1565ab6c82de5623d5a643c3a8bd1fa513f 100644 --- a/theory/Multipoles/gravity_derivatives.tex +++ b/theory/Multipoles/gravity_derivatives.tex @@ -1,5 +1,3 @@ -We use the multi-index notation of \cite{Dehnen2014} to simplify expressions. - \subsection{Derivatives of the gravitational potential} The calculation of all the diff --git a/theory/Multipoles/potential.py b/theory/Multipoles/potential.py index 9374dd4351cefb978ff8b5053563b5d3bb04c8cf..559f590762a3cbef171c5dd584cbc517879a2cec 100644 --- a/theory/Multipoles/potential.py +++ b/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") diff --git a/theory/Multipoles/potential_derivatives.tex b/theory/Multipoles/potential_derivatives.tex new file mode 100644 index 0000000000000000000000000000000000000000..84e7f9db02f6bf80ca066a2e29d2d2eb2b782569 --- /dev/null +++ b/theory/Multipoles/potential_derivatives.tex @@ -0,0 +1,45 @@ +\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} diff --git a/theory/Multipoles/potential_softening.tex b/theory/Multipoles/potential_softening.tex new file mode 100644 index 0000000000000000000000000000000000000000..1186a9cec377fd8daa94e14d024115f95ecfdc99 --- /dev/null +++ b/theory/Multipoles/potential_softening.tex @@ -0,0 +1,52 @@ +\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} diff --git a/theory/Multipoles/run.sh b/theory/Multipoles/run.sh index 34304a37dc1cb1fad3a7442544a8df793ea94d35..f25d407cd4ffe679a272f352798817f7c0c4e55a 100755 --- a/theory/Multipoles/run.sh +++ b/theory/Multipoles/run.sh @@ -1,5 +1,5 @@ #!/bin/bash -python kernels.py +python potential.py pdflatex -jobname=fmm fmm_standalone.tex bibtex fmm.aux pdflatex -jobname=fmm fmm_standalone.tex