Commit 9b541470 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added a figure to illustrate the contributions to gravity


Former-commit-id: 9666e58c8ffadff83d316fec3c68600fe9b97477
parent a09f2c4a
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")
......@@ -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}
......
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