diff --git a/theory/latex/Figures/Gravity/force.png b/theory/latex/Figures/Gravity/force.png new file mode 100644 index 0000000000000000000000000000000000000000..074aaf0dd82468f30f696cb89f52bfb1d3e689b1 Binary files /dev/null and b/theory/latex/Figures/Gravity/force.png differ diff --git a/theory/latex/Figures/Gravity/force.py b/theory/latex/Figures/Gravity/force.py new file mode 100644 index 0000000000000000000000000000000000000000..a829522f3b625b91528e7c78702393c56b6c4f2c --- /dev/null +++ b/theory/latex/Figures/Gravity/force.py @@ -0,0 +1,101 @@ +import matplotlib +matplotlib.use("Agg") +from pylab import * +# tex stuff +#rc('font',**{'family':'serif','serif':['Palatino']}) +params = {'axes.labelsize': 14, + 'text.fontsize': 14, + 'legend.fontsize': 14, + 'xtick.labelsize': 14, + 'ytick.labelsize': 14, + 'text.usetex': True, + + 'figure.figsize' : (9,6), + 'figure.subplot.left' : 0.08, # the left side of the subplots of the figure + 'figure.subplot.right' : 0.95 , # the right side of the subplots of the figure + 'figure.subplot.bottom' : 0.1 , # the bottom of the subplots of the figure + 'figure.subplot.top' : 0.93 , # the top of the subplots of the figure + 'figure.subplot.wspace' : 0.2 , # the amount of width reserved for blank space between subplots + 'figure.subplot.hspace' : 0.2 , # the amount of height reserved for white space between subplots + + 'lines.markersize' : 4, + + #'axes.formatter.limits' : (, 0), + + 'text.latex.unicode': True + } +rcParams.update(params) +rc('font', family='serif') +import sys +import math +import os + + + +epsilon = 1. +h = 2.8*epsilon +r_s = 20. +r_c = 4.5*r_s + +r=arange(0.1, 400, 1/1000.) + +numElem = size(r) +f = zeros(numElem) +fac = zeros(numElem) +kernel = zeros(numElem) + +for i in range(numElem): + if ( r[i] >= r_c ): + f[i] = 0. + elif ( r[i] >= h ): + f[i] = (1. / r[i]**3) + else: + u = r[i] / h + if( u < 0.5 ): + f[i] = (10.666666666667 + u * u * (32.0 * u - 38.4)) / h**3 + else: + f[i] = (21.333333333333 - 48.0 * u + 38.4 * u * u - 10.666666666667 * u * u * u - 0.066666666667 / (u * u * u)) / h**3 + + fac[i] = math.erfc( r[i] / ( 2. * r_s ) ) + ( r[i] / ( r_s * sqrt( math.pi ) ) ) * exp( - r[i] * r[i] / ( 4 * r_s * r_s ) ) + f[i] = f[i]*fac[i] + +for i in range(numElem): + u = r[i] / h + if( u < 0.5 ): + kernel[i] = (10.666666666667 + u * u * (32.0 * u - 38.4)) / h**3 + else: + kernel[i] = (21.333333333333 - 48.0 * u + 38.4 * u * u - 10.666666666667 * u * u * u - 0.066666666667 / (u * u * u)) / h**3 + +figure() +loglog(r/h, 1/r**3, 'b-', label="Newton's force") +loglog(r/h, kernel, 'r-', label="Softening kernel") +loglog(r/h, fac/r**3, 'g-', label="Tree force",) +loglog(r/h, f, 'k--', label="Total short-range force", linewidth=3) + + +plot([1, 1], [1e-10, 1e10], 'k--') +text(0.85, 7e-9, "$h=2.8\epsilon$", rotation="vertical", fontsize=14, va="bottom") + +plot([epsilon/h, epsilon/h], [1e-10, 1e10], 'k--') +text(0.85*epsilon/h, 7e-9, "$\epsilon$", rotation="vertical", fontsize=14, va="bottom") + +plot([r_s/h, r_s/h], [1e-10, 1e10], 'k--') +text(0.85*r_s/h, 7e-9, "$r_s$", rotation="vertical", fontsize=14, va="bottom") + +plot([r_c/h, r_c/h], [1e-10, 1e10], 'k--') +text(0.85*r_c/h, 7e-9, "$r_c$", rotation="vertical", fontsize=14, va="bottom") + + +legend(loc="upper right") + + +grid() + +xlim(1e-1, 200) +xlabel("r/h") + +ylim(4e-10, 1e2) +ylabel("Normalized acceleration") + + +savefig("force.png") diff --git a/theory/latex/swift.tex b/theory/latex/swift.tex index baf8cc656cd9a9996889ee0494b92eed518ce0a1..a93933665488076b1d5b9a51d3033aa5935a3c17 100755 --- a/theory/latex/swift.tex +++ b/theory/latex/swift.tex @@ -197,39 +197,29 @@ time. In cosmological simulations, the softening length is a function of time. \section{Gravity interactions} -The interaction is computed between all pairs of particles. The acceleration of particle $i$ due to one other particle -$j$ is given by: +The gravity interactions are split in two parts, the sort range interactions computed via the FMM method and the long +range interactions computed through a PM scheme. Short range interactions are also soften on a scale $\epsilon$ to +reduce the effect of spurious two-body relaxation. For implementation purposes, we define the quantity $h_c = +2.8\epsilon$, the radius at which the force becomes Newtonian and is not soften any more. + +The short range force is computed thanks to FMM. The expression for the acceleration of a particle $i$ due to a body +(particle or monopole) $j$ is given by: \begin{equation} - \vec{a_i} = \left\lbrace \begin{array}{rcl} - \Bigg. G_N\dfrac{m_j}{|\vec{r}_{ij}|^3}\vec{r}_{ij}& \mbox{if} & \Big.\dfrac{|\vec{r}_{ij}|}{\epsilon} -\geq 1 \\ - \Bigg. G_N\dfrac{m_j}{\epsilon^3}\phi\left(\frac{|\vec{r}_{ij}|}{\epsilon}\right)\vec{r}_{ij} & \mbox{if} & -\Big.\dfrac{|\vec{r}_{ij}|}{\epsilon} < 1 \\ - \end{array} -\right. + \vec{a} = f \cdot \vec{r}_{ij}, \end{equation} - -where the softened potential $\phi(u)$ is given by: + where $f$ is the norm of the acceleration given by \begin{equation} -\phi(u) = \left\lbrace \begin{array}{rcl} - \frac{32}{3} - \frac{192}{5}u^2 + 32u^3& \mbox{if} & u < \Big.\frac{1}{2}\\ - \frac{64}{3} - 48u + \frac{192}{5} u^2 - \frac{32}{3}u^3 - \frac{2}{30}\dfrac{1}{u^3}& \mbox{if} & \Big.u -\geq \frac{1}{2}\\ - \end{array} -\right. + f = G\left\lbrace \begin{array}{rcl} + \frac{1}{|\vec{r}_{ij}|^3} & \mbox{if} & |\vec{r}_{ij}| > h_c \\ + (\frac{32}{3} + \frac{|\vec{r}_{ij}|^2}{h_c^2} (32 \frac{|\vec{r}_{ij}|}{h_c} - 38.4)) / h_c^3 & +\mbox{if} & h_c \geq |\vec{r}_{ij}| > \frac{h_c}{2} \\ + \frac{1}{|\vec{r}_{ij}|^3} & \mbox{if} & |\vec{r}_{ij}| \leq \frac{h_c}{2} \\ + \end{array}\right. \end{equation} -and $G_N$ is Newton's constant for gravity. The softening length used in this equation is the maximum of the softening -length of both particles, i.e. $\epsilon = \max(\epsilon_i, \epsilon_j)$. Notice that this softened potential is -derived from the smoothing kernel used in the SPH part of the calculation with a smoothing length $h=1.4\epsilon$. This -corresponds to a Plummer sphere of size $\epsilon$. \\ - - -\textcolor{red}{This is the softening potential used in GADGET. There is something weird about this function that I want -to check but for now it should be sufficient.}\\ - + The maximal time step a particle is allowed to have is related to the acceleration by: \begin{equation}