### Re-factor a bit the theory document

parent 23ae1c60
 \subsection{The multipole acceptance criterion} The main remaining question is to decide when two cells are far enough from each others that the truncated Taylor expansion used as approximation for the potential (eq. \ref{eq:fmm:expansion}) is accurate enough. The criterion used to make that decision is called the \emph{multipole acceptance criterion} (MAC). \\ We know that (\ref{eq:fmm:expansion}) is converging towards the correct answer provided $1>|\mathbf{r}_a + \mathbf{r}_b| / |\mathbf{R}|$. This is hence the most basic (and always necessary) MAC that can be designed. If this ratio is lower, the accuracy (at a fixed expansion order) is improved and it is hence common practice to define a critical \emph{opening angle} $\theta_{\rm cr}$ and allow the use of the multipole approximation between two cells if \begin{equation} \theta_{\rm cr} > \frac{\rho_A + \rho_B} {|\mathbf{R}|}. \label{eq:fmm:angle} \end{equation} This lets users have a second handle on the accuracy on the gravity calculation besides the much more involved change in the expansion order $p$ of the FMM method. Typical values for the opening angle are in the range $[0.3, 0.7]$, with the cost of the simulation growing as $\theta_{\rm cr}$ decreases. \\ This method has the drawback of using a uniform criterion across the entire simulation volume and time evolution, which means that the chosen value of $\theta_{\rm cr}$ could be too small in some regions (leading to too many operations for the expected accuracy) and too large in some other other ones (leading to a lower level of accuracy than expected). \swift instead uses a more adaptive criterion to decide when the multipole approximation can be used. This is based on the error analysis of FMM by \cite{Dehnen2014} and is summarised below for completeness. The key idea is to exploit the additional information about the distribution of particles that is encoded in the higher-order multipole terms.\\ We start by defining the scalar quantity $P_{A,n}$, the \emph{power} of the multipole of order $n$ of the particles in cell $A$, via \begin{equation} P_{A,n}^2 = \sum_{|\mathbf{m}|=n} \frac{\mathbf{m}!}{|\mathbf{m}|!}\mathsf{M}_{A,\mathbf{m}}^2, \end{equation} where the sum runs over all the multipole terms of order $n$ in the cell\footnote{Note that $P_{0} \equiv \mathsf{M}_{(0,0,0)}$ is just the mass of the cell and since \swift uses the centre of mass as the centre of expansion of the multipoles, $P_{1} = 0$.}. This quantity is a simple upper bound for the amplitude of the multipole ($\mathsf{M}_{A, \mathbf{m}} < P_{A,|\mathbf{m}|}/|\mathbf{m}|!$) and can hence be used to estimate the importance of the terms of a given order in the Taylor series of the potential. Following \cite{Dehnen2014} we then consider a sink cell $A$ and a source cell $B$ (figure \ref{fig:fmm:cells}) for which we evaluate at order $p$ the scalar \begin{equation} E_{BA,p} = \frac{1}{M_B|\mathbf{R}|^p} \sum_{n=0}^p \binom{p}{n} P_{B,n} \rho_A^{p-n}, \label{eq:fmm:e_ab} \end{equation} with $M_B \equiv \mathsf{M}_{B,(0,0,0)}$, the sum of the mass of the particles in cell $B$. Note that since $P_{B,n} \leq M_B \rho_B^n$, we have $E_{BA, p} \leq \left((\rho_A + \rho_B)/|\mathbf{R}|\right)^p$, where the right-hand side is the expression used in the basic opening angle condition (\ref{eq:fmm:angle}). We finally scale the $E_{BA,p}$'s by the relative size of the two cells to define the error estimator $\tilde{E}_{BA,p}$: \begin{equation} \tilde{E}_{BA,p} = 8\frac{\max(\rho_A, \rho_B)}{\rho_A + \rho_B}E_{BA,p}. \label{eq:fmm:e_ab_tilde} \end{equation} As shown by \cite{Dehnen2014}, these quantities are excellent estimators of the error made in computing the accelerations between two cells using the M2L and M2P kernels at a given order. We can hence use this property to design a new MAC by demanding that the estimated acceleration error is no larger than a certain fraction of the smallest acceleration in the sink cell $A$. This means we can use the FMM approximation between to approximate the accelerations in cell $A$ due to the particles in cell $B$ if \begin{equation} \tilde{E}_{BA,p} \frac{M_B}{|\mathbf{R}|^2} < \epsilon_{\rm FMM} \min_{a\in A}\left(|\mathbf{a}_a|\right) \quad \rm{and} \quad \frac{\rho_A + \rho_B} {|\mathbf{R}|} < 1, \label{eq:fmm:mac} \end{equation} where the $\mathbf{a}_a$ are the accelerations of the particles in cell $A$ and $\epsilon_{\rm FMM}$ is a tolerance parameter. Since this is self-referencing (i.e. we need the accelerations to decide how to compute the accelerations), we need to use a an estimator of $|\mathbf{a}_a|$. In \swift, we follow the strategy used by \gadget and use the acceleration of the previous time-step\footnote{On the first time-step of a simulation this value has not been computed yet. We hence run a fake 0th time-step with the simpler MAC (eq. \ref{eq:fmm:angle}), which is good enough to obtain approximations of the accelerations.}. The minimal norm of the acceleration in a given cell can be computed at the same time as the P2M and M2M kernels are evaluated in the tree construction phase. The second condition in (\ref{eq:fmm:mac}) is necessary to ensure the convergence of the Taylor expansion.\\ One important difference between this criterion and the purely geometric one (\ref{eq:fmm:angle}) is that it is not symmetric in $A \leftrightarrow B$ (i.e. $E_{AB,p} \neq E_{BA,p}$). This implies that there are cases where a multipole in cell $A$ can be used to compute the field tensors in cell $B$ but the multipole in $B$ cannot be used to compute the $\mathsf{F}$ values of cell $A$ and vice versa. This affects the tree walk by breaking the symmetry and potentially leading to cells of different sizes interacting. \\ For the M2P kernel, the sink is a single particle $a$ and hence $\rho_A = 0$, which simplifies some of the expressions above. In this case, at order $p$, we get: \begin{equation} E_{BA,p} = \frac{P_{B,p}}{M_B |\mathbf{R}|^p}, \qquad \tilde{E}_{BA,p} = 8E_{BA,p} \nonumber \end{equation} Note that, in this case, only the power term of the order of the scheme appears; not a sum over the lower-order ones. This leads to the following MAC for the M2P kernel: \begin{equation} 8\frac{P_{B,p}}{|\mathbf{R}|^{p+2}} < \epsilon_{\rm FMM} |\mathbf{a}_a| \quad \rm{and} \quad \frac{\rho_B} {|\mathbf{R}|} < 1. \label{eq:fmm:mac_m2p} \end{equation} The value of $\epsilon_{\rm FMM}$ could in principle be different than the one used for the M2L MAC. One special case is of particular interest to link our expression to other results. Using the expression for order $2$ and the approximation $P_{B,p} \approx M_B \rho_B^p$, we get \begin{equation} 8\frac{M_B}{|\mathbf{R}|^2}\left(\frac{\rho_B}{|\mathbf{R}|}\right)^2 < \epsilon_{\rm FMM} |\mathbf{a}_a| \nonumber \end{equation} for our MAC. This is the same expression as the adaptive opening angle used by \gadget \cite[see eq.18 of][]{Springel2005} up to numerical factors and definition of the size of a multipole ($\rho$ vs. the cell edge). Note, however, that, in practice, since formally $P_{B,p} \leq M_B \rho_B^p$, the dependence is slightly different.\\ We conclude this section by noting that whilst the derivation of the FMM equations and of the simple geometric MAC (eq. \ref{eq:fmm:angle}) do not make any assumptions about the functional form of $\varphi(r)$, the more advanced MAC is valid in the specific case of the gravitational potential $\varphi(r) = m/r$ as can be inferred from the $m/r^2$ term appearing on the LHS of the criteria (\ref{eq:fmm:mac}) and (\ref{eq:fmm:mac_m2p}). \subsubsection{Modifications for softened and truncated gravity} \begin{figure} \includegraphics[width=\columnwidth]{mac_potential.pdf} \caption{The gravitational forces $f_{\rm SWIFT}$ computed by SWIFT (green line) including the force softening on the smallest scales and the long-range periodic mesh truncation on the largest scales for a simulation box of size $L$, a mesh scale-length $r_s$ and Plummer-equivalent softening $\epsilon_{\rm Plummer}$. The approximate fast estimator of the forces used in the MAC $f_{\rm MAC}$ is shown using yellow dash-dotted lines. Note that, by construction, $f_{\rm SWIFT} \leq f_{\rm MAC} \leq 1/r^2$ for all distances $r$.} \label{fig:fmm:mac_potential} \end{figure} One drawback of using expression (\ref{eq:fmm:mac}) in the case of a softened potential (or a potential truncated to apply long-range forces from a mesh (Sec. \ref{ssec:mesh_summary}) is that the $M/R^2$ term will overestimate the expected contribution from the multipole to the filed tensors, sometimes by large factors. This difference is shown on fig. \ref{fig:fmm:mac_potential}, with for instance a ratio of $3$ between the true forces and the Newtonian values reached a the scale of the Plummer softening. Using the simple expression (\ref{eq:fmm:mac}) will make the MAC too aggressive by preventing it from using a given multipole as it will be difficult to make the large term $M/R^2$ be below the fixed fraction $\epsilon_{\rm FMM}$ of the total acceleration of the receiving cell. This implies more computation as it will force the tree-walk algorithm to use more interactions by going to the daughter cells. The estimation of the contribution of the multipole in the MAC should hence be replaced by a more realistic term, closer to the one actually used in the interactions (eq. \ref{eq:fmm:force_norm}). In simulations with periodic boundary conditions, the same reasoning applies to the truncated force at the radii overlapping with the scale $r_s$ of the mesh forces. However, both the short- and long-range truncation functions are expensive to evaluate in the context of the MAC which is called a large number of times during a tree walk. We hence, construct a cheaper to evaluate estimator $f_{\rm MAC}$ that is closer to the true forces than the purely Newtonian term: \begin{align} f_{\rm MAC}(r) = \left\lbrace\begin{array}{rcl} \left(\frac{9}{5}\right)^2 H^{-2} & \mbox{if} & r < \frac{5}{9}H,\\ r^{-2} & \mbox{if} & \frac{5}{9}H \leq r < \frac{5}{3}r_s, \\ \left(\frac{9}{5}\right)^2 r_s^{-2} r^{-4} & \mbox{if} & \frac{5}{3}r_s \leq r. \\ \end{array} \right. \label{eq:fmm:f_mac} \end{align} This esimator is shown as a dot-dashed line on Fig. \ref{fig:fmm:mac_potential} and obeys the relation $f_{\rm SWIFT}(r) \leq f_{\rm MAC}(r) \leq 1/r^2$, with $f_{\rm SWIFT}(r)$ the true truncated and softened forces (green line).
 ... ... @@ -34,6 +34,7 @@ \input{fmm_summary} %\input{gravity_derivatives} \input{mesh_summary} \input{fmm_mac} \input{exact_forces} \bibliographystyle{mnras} ... ...
 ... ... @@ -224,140 +224,4 @@ any particle in that cell as \\ \textcolor{red}{MORE WORDS HERE} \subsubsection{The multipole acceptance criterion} The main remaining question is to decide when two cells are far enough from each others that the truncated Taylor expansion used as approximation for the potential (eq. \ref{eq:fmm:expansion}) is accurate enough. The criterion used to make that decision is called the \emph{multipole acceptance criterion} (MAC). \\ We know that (\ref{eq:fmm:expansion}) is converging towards the correct answer provided $1>|\mathbf{r}_a + \mathbf{r}_b| / |\mathbf{R}|$. This is hence the most basic (and always necessary) MAC that can be designed. If this ratio is lower, the accuracy (at a fixed expansion order) is improved and it is hence common practice to define a critical \emph{opening angle} $\theta_{\rm cr}$ and allow the use of the multipole approximation between two cells if \begin{equation} \theta_{\rm cr} > \frac{\rho_A + \rho_B} {|\mathbf{R}|}. \label{eq:fmm:angle} \end{equation} This lets users have a second handle on the accuracy on the gravity calculation besides the much more involved change in the expansion order $p$ of the FMM method. Typical values for the opening angle are in the range $[0.3, 0.7]$, with the cost of the simulation growing as $\theta_{\rm cr}$ decreases. \\ This method has the drawback of using a uniform criterion across the entire simulation volume and time evolution, which means that the chosen value of $\theta_{\rm cr}$ could be too small in some regions (leading to too many operations for the expected accuracy) and too large in some other other ones (leading to a lower level of accuracy than expected). \swift instead uses a more adaptive criterion to decide when the multipole approximation can be used. This is based on the error analysis of FMM by \cite{Dehnen2014} and is summarised below for completeness. The key idea is to exploit the additional information about the distribution of particles that is encoded in the higher-order multipole terms.\\ We start by defining the scalar quantity $P_{A,n}$, the \emph{power} of the multipole of order $n$ of the particles in cell $A$, via \begin{equation} P_{A,n}^2 = \sum_{|\mathbf{m}|=n} \frac{\mathbf{m}!}{|\mathbf{m}|!}\mathsf{M}_{A,\mathbf{m}}^2, \end{equation} where the sum runs over all the multipole terms of order $n$ in the cell\footnote{Note that $P_{0} \equiv \mathsf{M}_{(0,0,0)}$ is just the mass of the cell and since \swift uses the centre of mass as the centre of expansion of the multipoles, $P_{1} = 0$.}. This quantity is a simple upper bound for the amplitude of the multipole ($\mathsf{M}_{A, \mathbf{m}} < P_{A,|\mathbf{m}|}/|\mathbf{m}|!$) and can hence be used to estimate the importance of the terms of a given order in the Taylor series of the potential. Following \cite{Dehnen2014} we then consider a sink cell $A$ and a source cell $B$ (figure \ref{fig:fmm:cells}) for which we evaluate at order $p$ the scalar \begin{equation} E_{BA,p} = \frac{1}{M_B|\mathbf{R}|^p} \sum_{n=0}^p \binom{p}{n} P_{B,n} \rho_A^{p-n}, \label{eq:fmm:e_ab} \end{equation} with $M_B \equiv \mathsf{M}_{B,(0,0,0)}$, the sum of the mass of the particles in cell $B$. Note that since $P_{B,n} \leq M_B \rho_B^n$, we have $E_{BA, p} \leq \left((\rho_A + \rho_B)/|\mathbf{R}|\right)^p$, where the right-hand side is the expression used in the basic opening angle condition (\ref{eq:fmm:angle}). We finally scale the $E_{BA,p}$'s by the relative size of the two cells to define the error estimator $\tilde{E}_{BA,p}$: \begin{equation} \tilde{E}_{BA,p} = 8\frac{\max(\rho_A, \rho_B)}{\rho_A + \rho_B}E_{BA,p}. \label{eq:fmm:e_ab_tilde} \end{equation} As shown by \cite{Dehnen2014}, these quantities are excellent estimators of the error made in computing the accelerations between two cells using the M2L and M2P kernels at a given order. We can hence use this property to design a new MAC by demanding that the estimated acceleration error is no larger than a certain fraction of the smallest acceleration in the sink cell $A$. This means we can use the FMM approximation between to approximate the accelerations in cell $A$ due to the particles in cell $B$ if \begin{equation} \tilde{E}_{BA,p} \frac{M_B}{|\mathbf{R}|^2} < \epsilon_{\rm FMM} \min_{a\in A}\left(|\mathbf{a}_a|\right) \quad \rm{and} \quad \frac{\rho_A + \rho_B} {|\mathbf{R}|} < 1, \label{eq:fmm:mac} \end{equation} where the $\mathbf{a}_a$ are the accelerations of the particles in cell $A$ and $\epsilon_{\rm FMM}$ is a tolerance parameter. Since this is self-referencing (i.e. we need the accelerations to decide how to compute the accelerations), we need to use a an estimator of $|\mathbf{a}_a|$. In \swift, we follow the strategy used by \gadget and use the acceleration of the previous time-step\footnote{On the first time-step of a simulation this value has not been computed yet. We hence run a fake 0th time-step with the simpler MAC (eq. \ref{eq:fmm:angle}), which is good enough to obtain approximations of the accelerations.}. The minimal norm of the acceleration in a given cell can be computed at the same time as the P2M and M2M kernels are evaluated in the tree construction phase. The second condition in (\ref{eq:fmm:mac}) is necessary to ensure the convergence of the Taylor expansion.\\ One important difference between this criterion and the purely geometric one (\ref{eq:fmm:angle}) is that it is not symmetric in $A \leftrightarrow B$ (i.e. $E_{AB,p} \neq E_{BA,p}$). This implies that there are cases where a multipole in cell $A$ can be used to compute the field tensors in cell $B$ but the multipole in $B$ cannot be used to compute the $\mathsf{F}$ values of cell $A$ and vice versa. This affects the tree walk by breaking the symmetry and potentially leading to cells of different sizes interacting. \\ For the M2P kernel, the sink is a single particle $a$ and hence $\rho_A = 0$, which simplifies some of the expressions above. In this case, at order $p$, we get: \begin{equation} E_{BA,p} = \frac{P_{B,p}}{M_B |\mathbf{R}|^p}, \qquad \tilde{E}_{BA,p} = 8E_{BA,p} \nonumber \end{equation} Note that, in this case, only the power term of the order of the scheme appears; not a sum over the lower-order ones. This leads to the following MAC for the M2P kernel: \begin{equation} 8\frac{P_{B,p}}{|\mathbf{R}|^{p+2}} < \epsilon_{\rm FMM} |\mathbf{a}_a| \quad \rm{and} \quad \frac{\rho_B} {|\mathbf{R}|} < 1. \label{eq:fmm:mac_m2p} \end{equation} The value of $\epsilon_{\rm FMM}$ could in principle be different than the one used for the M2L MAC. One special case is of particular interest to link our expression to other results. Using the expression for order $2$ and the approximation $P_{B,p} \approx M_B \rho_B^p$, we get \begin{equation} 8\frac{M_B}{|\mathbf{R}|^2}\left(\frac{\rho_B}{|\mathbf{R}|}\right)^2 < \epsilon_{\rm FMM} |\mathbf{a}_a| \nonumber \end{equation} for our MAC. This is the same expression as the adaptive opening angle used by \gadget \cite[see eq.18 of][]{Springel2005} up to numerical factors and definition of the size of a multipole ($\rho$ vs. the cell edge). Note, however, that, in practice, since formally $P_{B,p} \leq M_B \rho_B^p$, the dependence is slightly different.\\ We conclude this section by noting that whilst the derivation of the FMM equations and of the simple geometric MAC (eq. \ref{eq:fmm:angle}) do not make any assumptions about the functional form of $\varphi(r)$, the more advanced MAC is valid in the specific case of the gravitational potential $\varphi(r) = m/r$ as can be inferred from the $m/r^2$ term appearing on the LHS of the criteria (\ref{eq:fmm:mac}) and (\ref{eq:fmm:mac_m2p}).
 ... ... @@ -62,7 +62,7 @@ f'(\frac{r}{H}) \times H^{-2} & \mbox{if} & r < H,\\ r^{-2} & \mbox{if} & r \geq H. \end{array} \right. \label{eq:fmm:force} \label{eq:fmm:force_norm} \end{align} The softened density profile, its corresponding potential and resulting forces are shown on Fig. \ref{fig:fmm:softening} (for more ... ...
 ... ... @@ -14,6 +14,11 @@ then echo "Generating derivative figures..." python plot_derivatives.py fi if [ ! -e mac_potential.pdf ] then echo "Generating derivative figures..." python3 plot_mac_potential.py fi echo "Generating PDF..." pdflatex -jobname=fmm fmm_standalone.tex bibtex fmm.aux ... ...
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment