Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
SWIFT
SWIFTsim
Commits
3196f542
Commit
3196f542
authored
Jul 31, 2017
by
Matthieu Schaller
Browse files
Merge branch 'gravity_documentation' into 'master'
Documentation of the FMM-MM code See merge request
!387
parents
6775e2f3
56d84a63
Changes
10
Hide whitespace changes
Inline
Side-by-side
.gitignore
View file @
3196f542
...
...
@@ -102,6 +102,9 @@ theory/paper_pasc/pasc_paper.pdf
theory/Multipoles/fmm.pdf
theory/Multipoles/fmm_standalone.pdf
theory/Multipoles/potential.pdf
theory/Multipoles/potential_long.pdf
theory/Multipoles/potential_short.pdf
theory/Multipoles/force_short.pdf
m4/libtool.m4
m4/ltoptions.m4
...
...
theory/Multipoles/cells.pdf
View file @
3196f542
No preview for this file type
theory/Multipoles/fmm_summary.tex
View file @
3196f542
...
...
@@ -6,7 +6,7 @@ evaluate for each particle in a system the potential and associated
forces generated by all the other particles. Mathematically, this means
evaluate
\begin{equation}
\phi
(
\mathbf
{
x
}_
a) =
\sum
_{
b
\neq
a
}
G m
_
b
\phi
(
\mathbf
{
x
}_
a -
\phi
(
\mathbf
{
x
}_
a) =
\sum
_{
b
\neq
a
}
G m
_
b
\
var
phi
(
\mathbf
{
x
}_
a -
\mathbf
{
x
}_
b)
\qquad
\forall
~a
\in
N
\label
{
eq:fmm:n
_
body
}
\end{equation}
...
...
@@ -23,10 +23,11 @@ arising from nearby particles. \\
In what follows, we use the compact multi-index notation of
\cite
{
Dehnen2014
}
(repeated in appendix
\ref
{
sec:multi
_
index
_
notation
}
for completeness) to simplify expressions.
$
\mathbf
{
k
}$
,
$
\mathbf
{
m
}$
and
$
\mathbf
{
n
}$
are multi-indices and
$
\mathbf
{
r
}$
,
$
\mathbf
{
R
}$
,
$
\mathbf
{
x
}$
,
$
\mathbf
{
y
}$
and
$
\mathbf
{
z
}$
are vectors, whilst
$
a
$
and
$
b
$
are particle indices.
\\
for completeness) to simplify expressions and ease
comparisons.
$
\mathbf
{
k
}$
,
$
\mathbf
{
m
}$
and
$
\mathbf
{
n
}$
are
multi-indices and
$
\mathbf
{
r
}$
,
$
\mathbf
{
R
}$
,
$
\mathbf
{
x
}$
,
$
\mathbf
{
y
}$
and
$
\mathbf
{
z
}$
are vectors, whilst
$
a
$
and
$
b
$
are
particle indices.
\\
\begin{figure}
\includegraphics
[width=\columnwidth]
{
cells.pdf
}
...
...
@@ -46,34 +47,34 @@ with centres of mass $\mathbf{z}_A$ and $\mathbf{z}_B$
respectively, as shown on Fig.~
\ref
{
fig:fmm:cells
}
, the potential
generated by
$
b
$
at the location of
$
a
$
can be rewritten as
\begin{align}
\phi
(
\mathbf
{
x
}_
a -
\mathbf
{
x
}_
b)
&
=
\phi\left
(
\mathbf
{
x
}_
a -
\mathbf
{
z
}_
A -
\mathbf
{
x
}_
b +
\
var
phi
(
\mathbf
{
x
}_
a -
\mathbf
{
x
}_
b)
&
=
\
var
phi\left
(
\mathbf
{
x
}_
a -
\mathbf
{
z
}_
A -
\mathbf
{
x
}_
b +
\mathbf
{
z
}_
B +
\mathbf
{
z
}_
A -
\mathbf
{
z
}_
B
\right
)
\nonumber
\\
&
=
\phi\left
(
\mathbf
{
r
}_
a -
\mathbf
{
r
}_
b +
\mathbf
{
R
}
\right
)
&
=
\
var
phi\left
(
\mathbf
{
r
}_
a -
\mathbf
{
r
}_
b +
\mathbf
{
R
}
\right
)
\nonumber
\\
&
=
\sum
_
\mathbf
{
k
}
\frac
{
1
}{
\mathbf
{
k
}
!
}
\left
(
\mathbf
{
r
}_
a -
\mathbf
{
r
}_
b
\right
)
^{
\mathbf
{
k
}}
\nabla
^{
\mathbf
{
k
}}
\phi
(
\mathbf
{
R
}
)
\mathbf
{
r
}_
b
\right
)
^{
\mathbf
{
k
}}
\nabla
^{
\mathbf
{
k
}}
\
var
phi
(
\mathbf
{
R
}
)
\nonumber
\\
&
=
\sum
_
\mathbf
{
k
}
\frac
{
1
}{
\mathbf
{
k
}
!
}
\sum
_{
\mathbf
{
n
}
<
\mathbf
{
k
}}
\binom
{
\mathbf
{
k
}}{
\mathbf
{
n
}}
\mathbf
{
r
}_
a
^{
\mathbf
{
n
}}
\left
(-
\mathbf
{
r
}_
b
\right
)
^{
\mathbf
{
k
}
-
\mathbf
{
n
}}
\nabla
^{
\mathbf
{
k
}}
\phi
(
\mathbf
{
R
}
)
\nonumber
\\
\nabla
^{
\mathbf
{
k
}}
\
var
phi
(
\mathbf
{
R
}
)
\nonumber
\\
&
=
\sum
_
\mathbf
{
n
}
\frac
{
1
}{
\mathbf
{
n
}
!
}
\mathbf
{
r
}_
a
^{
\mathbf
{
n
}}
\sum
_
\mathbf
{
m
}
\frac
{
1
}{
\mathbf
{
m
}
!
}
\left
(-
\mathbf
{
r
}_
b
\right
)
^
\mathbf
{
m
}
\nabla
^{
\mathbf
{
n
}
+
\mathbf
{
m
}}
\phi
(
\mathbf
{
R
}
),
\left
(-
\mathbf
{
r
}_
b
\right
)
^
\mathbf
{
m
}
\nabla
^{
\mathbf
{
n
}
+
\mathbf
{
m
}}
\
var
phi
(
\mathbf
{
R
}
),
\end{align}
where we used the Taylor expansion of
$
\phi
$
around
$
\mathbf
{
R
}
\equiv
where we used the Taylor expansion of
$
\
var
phi
$
around
$
\mathbf
{
R
}
\equiv
\mathbf
{
z
}_
A
-
\mathbf
{
z
}_
B
$
on the third line, used
$
\mathbf
{
r
}_
a
\equiv
\mathbf
{
x
}_
a
-
\mathbf
{
z
}_
A
$
,
$
\mathbf
{
r
}_
b
\equiv
\mathbf
{
x
}_
b
-
\mathbf
{
z
}_
B
$
throughout and defined
$
\mathbf
{
m
}
\equiv
\mathbf
{
k
}
-
\mathbf
{
n
}$
on the last line. Expanding the series only up
to order
$
p
$
, we get
\begin{equation}
\phi
(
\mathbf
{
x
}_
a -
\mathbf
{
x
}_
b)
\approx
\sum
_{
\mathbf
{
n
}}^{
p
}
\
var
phi
(
\mathbf
{
x
}_
a -
\mathbf
{
x
}_
b)
\approx
\sum
_{
\mathbf
{
n
}}^{
p
}
\frac
{
1
}{
\mathbf
{
n
}
!
}
\mathbf
{
r
}_
a
^{
\mathbf
{
n
}}
\sum
_{
\mathbf
{
m
}}^{
p
-|
\mathbf
{
n
}
|
}
\frac
{
1
}{
\mathbf
{
m
}
!
}
\left
(-
\mathbf
{
r
}_
b
\right
)
^
\mathbf
{
m
}
\nabla
^{
\mathbf
{
n
}
+
\mathbf
{
m
}}
\phi
(
\mathbf
{
R
}
),
\nabla
^{
\mathbf
{
n
}
+
\mathbf
{
m
}}
\
var
phi
(
\mathbf
{
R
}
),
\label
{
eq:fmm:fmm
_
one
_
part
}
\end{equation}
with the approximation converging as
$
p
\rightarrow\infty
$
towards the
...
...
@@ -82,13 +83,13 @@ correct value provided $|\mathbf{R}|<|\mathbf{r}_a +
combine their contributions to the potential at location
$
\mathbf
{
x
}_
a
$
in cell
$
A
$
, we get
\begin{align}
\phi
_{
BA
}
(
\mathbf
{
x
}_
a)
&
=
\sum
_{
b
\in
B
}
m
_
b
\phi
(
\mathbf
{
x
}_
a -
\phi
_{
BA
}
(
\mathbf
{
x
}_
a)
&
=
\sum
_{
b
\in
B
}
G
m
_
b
\
var
phi
(
\mathbf
{
x
}_
a -
\mathbf
{
x
}_
b)
\label
{
eq:fmm:fmm
_
one
_
cell
}
\\
&
\approx
\sum
_{
\mathbf
{
n
}}^{
p
}
&
\approx
G
\sum
_{
\mathbf
{
n
}}^{
p
}
\frac
{
1
}{
\mathbf
{
n
}
!
}
\mathbf
{
r
}_
a
^{
\mathbf
{
n
}}
\sum
_{
\mathbf
{
m
}}
^{
p -|
\mathbf
{
n
}
|
}
\frac
{
1
}{
\mathbf
{
m
}
!
}
\sum
_{
b
\in
B
}
m
_
b
\left
(-
\mathbf
{
r
}_
b
\right
)
^
\mathbf
{
m
}
\nabla
^{
\mathbf
{
n
}
+
\mathbf
{
m
}}
\phi
(
\mathbf
{
R
}
)
\nonumber
.
\nabla
^{
\mathbf
{
n
}
+
\mathbf
{
m
}}
\
var
phi
(
\mathbf
{
R
}
)
\nonumber
.
\end{align}
This last equation forms the basis of the FMM. The algorithm
decomposes the equation into three separated sums evaluated at
...
...
@@ -106,12 +107,12 @@ compute the second kernel, \textsc{M2L} (multipole to local
expansion), which corresponds to the interaction of a cell with
another one:
\begin{equation}
\mathsf
{
F
}_{
\mathbf
{
n
}}
(
\mathbf
{
z
}_
A) =
\sum
_{
\mathbf
{
m
}}^{
p -|
\mathbf
{
n
}
|
}
\mathsf
{
F
}_{
\mathbf
{
n
}}
(
\mathbf
{
z
}_
A) =
G
\sum
_{
\mathbf
{
m
}}^{
p -|
\mathbf
{
n
}
|
}
\mathsf
{
M
}_{
\mathbf
{
m
}}
(
\mathbf
{
z
}_
B)
\mathsf
{
D
}_{
\mathbf
{
n
}
+
\mathbf
{
m
}}
(
\mathbf
{
R
}
),
\label
{
eq:fmm:M2L
}
\end{equation}
where
$
\mathsf
{
D
}_{
\mathbf
{
n
}
+
\mathbf
{
m
}}
(
\mathbf
{
R
}
)
\equiv
\nabla
^{
\mathbf
{
n
}
+
\mathbf
{
m
}}
\phi
(
\mathbf
{
R
}
)
$
is an order
$
n
+
m
$
\nabla
^{
\mathbf
{
n
}
+
\mathbf
{
m
}}
\
var
phi
(
\mathbf
{
R
}
)
$
is an order
$
n
+
m
$
derivative of the potential. This is the computationally expensive
step of the FMM algorithm as the number of operations in a naive
implementation using cartesian coordinates scales as
...
...
@@ -165,8 +166,8 @@ multiplications (provided $\mathsf{D}$ can be evaluated quickly, see
Sec.~
\ref
{
ssec:grav
_
derivatives
}
), which are extremely efficient
instructions on modern architectures. However, the fully expanded sums
can lead to rather large and prone to typo expressions. To avoid any
mis
s
hap, we use a
\texttt
{
python
}
to generate C code in which
all the
sums are unrolled and correct by construction. In
\swift
, we
mishap
s
, we use a
\texttt
{
python
}
script
to generate C code in which
all the
sums are unrolled and correct by construction. In
\swift
, we
implemented the kernels up to order
$
p
=
5
$
, as it proved to be accurate
enough for our purpose, but this could be extended to higher order
easily. This implies storing
$
56
$
numbers per cell for each
...
...
theory/Multipoles/gravity_derivatives.tex
View file @
3196f542
...
...
@@ -2,36 +2,36 @@
\label
{
ssec:grav
_
derivatives
}
The calculation of all the
$
\mathsf
{
D
}_
\mathbf
{
n
}
(
x,y,z
)
\equiv
\nabla
^{
\mathbf
{
n
}}
\phi
(
x,y,z
)
$
terms up
$
\mathsf
{
D
}_
\mathbf
{
n
}
(
x,y,z
)
\equiv
\nabla
^{
\mathbf
{
n
}}
\
var
phi
(
x,y,z
)
$
terms up
to the relevent order can be quite tedious and it is beneficial to
automatize the whole setup. Ideally, one would like to have an
expression for each of these terms that is only made of multiplications
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
}
achieve this by writing
$
\
var
phi
$
as a composition of functions
$
\
var
phi
(
u
(
x,y,z
))
$
and apply the
\textit
{
Fa
\`
a di Bruno
}
formula
\citep
[i.e. the ``chain rule'' for higher order derivatives,
see e.g.][]
{
Hardy2006
}
to construct our terms:
\begin{equation}
\label
{
eq:faa
_
di
_
bruno
}
\frac
{
\partial
^
n
}{
\partial
x
_
1
\cdots
\partial
x
_
n
}
\phi
(u)
=
\sum
_{
A
}
\phi
^{
(|A|)
}
(u)
\prod
_{
B
\in
\frac
{
\partial
^
n
}{
\partial
x
_
1
\cdots
\partial
x
_
n
}
\
var
phi
(u)
=
\sum
_{
A
}
\
var
phi
^{
(|A|)
}
(u)
\prod
_{
B
\in
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
$
\lbrace
1
,
\cdots
, n
\rbrace
$
,
$
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
cardinality of a set. For generic functions
$
\
var
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.
}
.
\\
For the un-softened gravitational potential, we choose to write
\begin{align}
\phi
(x,y,z)
&
= 1 /
\sqrt
{
u(x,y,z)
}
,
\\
\
var
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
This choice allows to have derivatives of any order of
$
\
var
phi
(
u
)
$
that
can be easily expressed and only depend on powers of
$
u
$
:
\begin{equation}
\phi
^{
(n)
}
(u) = (-1)
^
n
\cdot\frac
{
(2n-1)!!
}{
2
^
n
}
\cdot\frac
{
1
}{
u
^{
n+
\frac
{
1
}{
2
}}}
,
\
var
phi
^{
(n)
}
(u) = (-1)
^
n
\cdot\frac
{
(2n-1)!!
}{
2
^
n
}
\cdot\frac
{
1
}{
u
^{
n+
\frac
{
1
}{
2
}}}
,
\end{equation}
where
$
!!
$
denotes the semi-factorial. More importantly, this
choice of decomposition allows us to have non-zero derivatives of
$
u
$
...
...
@@ -39,7 +39,7 @@ 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
$
\mathsf
{
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
\
var
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 combination of derivatives
...
...
theory/Multipoles/mesh_summary.tex
View file @
3196f542
\subsection
{
Coupling the FMM to a mesh for periodic long-range forces
}
\label
{
ssec:mesh
_
summary
}
\begin{equation}
S(x) =
\frac
{
e
^
x
}{
1 + e
^
x
}
\end{equation}
\begin{align}
\varphi
_
s(r)
&
=
\frac
{
1
}{
r
}
\left
[2 - 2S\left(\frac{2r}{r_s}\right)\right]
\nonumber\\
&
=
\frac
{
1
}{
r
}
\left
[2 - \frac{2e^{\frac{2r}{r_s}}}{1+e^{\frac{2r}{r_s}}}\right]
\end{align}
\begin{align}
|
\mathbf
{
f
}_
s(r)|
&
=
\frac
{
1
}{
r
^
2
}
\left
[\frac{4r}{r_s}S'\left(\frac{2r}{r_s}\right) - 2S\left(\frac{2r}{r_s}\right) + 2\right]
\nonumber
\\
&
=
\frac
{
1
}{
r
^
2
}
\left
[\frac{4r}{r_s}\frac{e^{\frac{2r}{r_s}}}{(1+e^{\frac{2r}{r_s}})^2} - \frac{2e^{\frac{2r}{r_s}}}{1+e^{\frac{2r}{r_s}}} + 2\right]
\end{align}
\begin{equation}
\tilde\varphi
_
l(k) =
\frac
{
1
}{
k
^
2
}
\left
[\frac{\upi}{2}kr_s\textrm{csch}\left(\frac{\upi}{2}kr_s\right) \right]
\end{equation}
\begin{figure}
\includegraphics
[width=\columnwidth]
{
potential
_
short.pdf
}
\caption
{
aa
}
\label
{
fig:fmm:potential
_
short
}
\end{figure}
\begin{figure}
\includegraphics
[width=\columnwidth]
{
force
_
short.pdf
}
\caption
{
bb
}
\label
{
fig:fmm:force
_
short
}
\end{figure}
\begin{figure}
\includegraphics
[width=\columnwidth]
{
potential
_
long.pdf
}
\caption
{
cc
}
\label
{
fig:fmm:potential
_
long
}
\end{figure}
theory/Multipoles/plot_mesh.py
0 → 100644
View file @
3196f542
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
##############################################################################
import
matplotlib
matplotlib
.
use
(
"Agg"
)
from
pylab
import
*
from
scipy
import
integrate
from
scipy
import
special
from
scipy.optimize
import
curve_fit
from
scipy.optimize
import
fsolve
from
matplotlib.font_manager
import
FontProperties
import
numpy
import
math
params
=
{
'axes.labelsize'
:
9
,
'axes.titlesize'
:
10
,
'font.size'
:
10
,
'legend.fontsize'
:
10
,
'xtick.labelsize'
:
8
,
'ytick.labelsize'
:
8
,
'text.usetex'
:
True
,
'figure.figsize'
:
(
3.15
,
3.15
),
'figure.subplot.left'
:
0.12
,
'figure.subplot.right'
:
0.99
,
'figure.subplot.bottom'
:
0.09
,
'figure.subplot.top'
:
0.99
,
'figure.subplot.wspace'
:
0.
,
'figure.subplot.hspace'
:
0.
,
'lines.markersize'
:
6
,
'lines.linewidth'
:
3.
,
'text.latex.unicode'
:
True
}
rcParams
.
update
(
params
)
rc
(
'font'
,
**
{
'family'
:
'sans-serif'
,
'sans-serif'
:[
'Times'
]})
colors
=
[
'#4477AA'
,
'#CC6677'
,
'#DDCC77'
,
'#117733'
]
# Parameters
r_s
=
2.
r_min
=
1e-2
r_max
=
1.5e2
# Radius
r
=
logspace
(
log10
(
r_min
),
log10
(
r_max
),
401
)
r_rs
=
r
/
r_s
k
=
logspace
(
log10
(
r_min
/
r_s
**
2
),
log10
(
r_max
/
r_s
**
2
),
401
)
k_rs
=
k
*
r_s
# Newtonian solution
phi_newton
=
1.
/
r
phit_newton
=
1.
/
k
**
2
force_newton
=
1.
/
r
**
2
def
my_exp
(
x
):
return
1.
+
x
+
(
x
**
2
/
2.
)
+
(
x
**
3
/
6.
)
+
(
x
**
4
/
24.
)
+
(
x
**
5
/
120.
)
# + (x**6 / 720.)
#return exp(x)
def
csch
(
x
):
# hyperbolic cosecant
return
1.
/
sinh
(
x
)
def
sigmoid
(
x
):
return
my_exp
(
x
)
/
(
my_exp
(
x
)
+
1.
)
def
d_sigmoid
(
x
):
return
my_exp
(
x
)
/
((
my_exp
(
x
)
+
1
)
**
2
)
def
swift_corr
(
x
):
return
2
*
sigmoid
(
4
*
x
)
-
1
#figure()
#x = linspace(-4, 4, 100)
#plot(x, special.erf(x), '-', color=colors[0])
#plot(x, swift_corr(x), '-', color=colors[1])
#plot(x, x, '-', color=colors[2])
#ylim(-1.1, 1.1)
#xlim(-4.1, 4.1)
#savefig("temp.pdf")
# Correction in real space
corr_short_gadget2
=
special
.
erf
(
r
/
(
2.
*
r_s
))
corr_short_swift
=
swift_corr
(
r
/
(
2.
*
r_s
))
eta_short_gadget2
=
special
.
erfc
(
r
/
2.
*
r_s
)
+
(
r
/
(
r_s
*
math
.
sqrt
(
math
.
pi
)))
*
exp
(
-
r
**
2
/
(
4.
*
r_s
**
2
))
eta_short_swift
=
4.
*
(
r
/
r_s
)
*
d_sigmoid
(
2.
*
r
/
r_s
)
-
2.
*
sigmoid
(
2
*
r
/
r_s
)
+
2.
# Corection in Fourier space
corr_long_gadget2
=
exp
(
-
k
**
2
*
r_s
**
2
)
corr_long_swift
=
math
.
pi
*
k
*
r_s
*
csch
(
0.5
*
math
.
pi
*
r_s
*
k
)
/
2.
# Shortrange term
phi_short_gadget2
=
(
1.
/
r
)
*
(
1.
-
corr_short_gadget2
)
phi_short_swift
=
(
1.
/
r
)
*
(
1.
-
corr_short_swift
)
force_short_gadget2
=
(
1.
/
r
**
2
)
*
eta_short_gadget2
force_short_swift
=
(
1.
/
r
**
2
)
*
eta_short_swift
# Long-range term
phi_long_gadget2
=
(
1.
/
r
)
*
corr_short_gadget2
phi_long_swift
=
(
1.
/
r
)
*
corr_short_swift
phit_long_gadget2
=
corr_long_gadget2
/
k
**
2
phit_long_swift
=
corr_long_swift
/
k
**
2
figure
()
# Potential
subplot
(
311
,
xscale
=
"log"
,
yscale
=
"log"
)
plot
(
r_rs
,
phi_newton
,
'--'
,
lw
=
1.4
,
label
=
"${
\\
rm Newtonian}$"
,
color
=
colors
[
0
])
plot
(
r_rs
,
phi_short_gadget2
,
'-'
,
lw
=
1.4
,
label
=
"${
\\
rm Gadget}$"
,
color
=
colors
[
2
])
plot
(
r_rs
,
phi_short_swift
,
'-'
,
lw
=
1.4
,
label
=
"${
\\
rm SWIFT}$"
,
color
=
colors
[
3
])
plot
([
1.
,
1.
],
[
1e-5
,
1e5
],
'k-'
,
alpha
=
0.5
,
lw
=
0.5
)
xlim
(
1.1
*
r_min
/
r_s
,
0.9
*
r_max
/
r_s
)
ylim
(
1.1
/
r_max
,
0.9
/
r_min
)
ylabel
(
"$
\\
varphi_s(r)$"
,
labelpad
=-
3
)
legend
(
loc
=
"upper right"
,
frameon
=
True
,
handletextpad
=
0.1
,
handlelength
=
3.2
,
fontsize
=
8
)
# Correction
subplot
(
312
,
xscale
=
"log"
,
yscale
=
"log"
)
plot
(
r_rs
,
np
.
ones
(
np
.
size
(
r
)),
'--'
,
lw
=
1.4
,
color
=
colors
[
0
])
plot
(
r_rs
,
1.
-
corr_short_gadget2
,
'-'
,
lw
=
1.4
,
color
=
colors
[
2
])
plot
(
r_rs
,
1.
-
corr_short_swift
,
'-'
,
lw
=
1.4
,
color
=
colors
[
3
])
plot
(
r_rs
,
np
.
ones
(
np
.
size
(
r
))
*
0.01
,
'k:'
,
alpha
=
0.5
,
lw
=
0.5
)
plot
([
1.
,
1.
],
[
-
1e5
,
1e5
],
'k-'
,
alpha
=
0.5
,
lw
=
0.5
)
yticks
([
1e-2
,
1e-1
,
1
],
[
"$0.01$"
,
"$0.1$"
,
"$1$"
])
xlim
(
1.1
*
r_min
/
r_s
,
0.9
*
r_max
/
r_s
)
ylim
(
3e-3
,
1.5
)
#ylabel("$\\chi_s(r)$", labelpad=-3)
ylabel
(
"$
\\
varphi_s(r)
\\
times r$"
,
labelpad
=-
2
)
# 1 - Correction
subplot
(
313
,
xscale
=
"log"
,
yscale
=
"log"
)
plot
(
r_rs
,
corr_short_gadget2
,
'-'
,
lw
=
1.4
,
color
=
colors
[
2
])
plot
(
r_rs
,
corr_short_swift
,
'-'
,
lw
=
1.4
,
color
=
colors
[
3
])
plot
([
1.
,
1.
],
[
1e-5
,
1e5
],
'k-'
,
alpha
=
0.5
,
lw
=
0.5
)
plot
(
r_rs
,
np
.
ones
(
np
.
size
(
r
)),
'k:'
,
alpha
=
0.5
,
lw
=
0.5
)
plot
(
r_rs
,
np
.
ones
(
np
.
size
(
r
))
*
0.01
,
'k:'
,
alpha
=
0.5
,
lw
=
0.5
)
xlim
(
1.1
*
r_min
/
r_s
,
0.9
*
r_max
/
r_s
)
ylim
(
3e-3
,
1.5
)
#ylabel("$1 - \\chi_s(r)$", labelpad=-2)
ylabel
(
"$1 -
\\
varphi_s(r)
\\
times r$"
,
labelpad
=-
2
)
yticks
([
1e-2
,
1e-1
,
1
],
[
"$0.01$"
,
"$0.1$"
,
"$1$"
])
xlabel
(
"$r / r_s$"
,
labelpad
=-
3
)
savefig
(
"potential_short.pdf"
)
##################################################################################################
# Force
figure
()
subplot
(
311
,
xscale
=
"log"
,
yscale
=
"log"
)
plot
(
r_rs
,
force_newton
,
'--'
,
lw
=
1.4
,
label
=
"${
\\
rm Newtonian}$"
,
color
=
colors
[
0
])
plot
(
r_rs
,
force_short_gadget2
,
'-'
,
lw
=
1.4
,
label
=
"${
\\
rm Gadget}$"
,
color
=
colors
[
2
])
plot
(
r_rs
,
force_short_swift
,
'-'
,
lw
=
1.4
,
label
=
"${
\\
rm SWIFT}$"
,
color
=
colors
[
3
])
plot
([
1.
,
1.
],
[
1e-5
,
1e5
],
'k-'
,
alpha
=
0.5
,
lw
=
0.5
)
xlim
(
1.1
*
r_min
/
r_s
,
0.9
*
r_max
/
r_s
)
ylim
(
1.1
/
r_max
**
2
,
0.9
/
r_min
**
2
)
ylabel
(
"$|
\\
mathbf{f}_s(r)|$"
,
labelpad
=-
3
)
yticks
([
1e-4
,
1e-2
,
1e0
,
1e2
],
[
"$10^{-4}$"
,
"$10^{-2}$"
,
"$10^{0}$"
,
"$10^{2}$"
])
legend
(
loc
=
"upper right"
,
frameon
=
True
,
handletextpad
=
0.1
,
handlelength
=
3.2
,
fontsize
=
8
)
# Correction
subplot
(
312
,
xscale
=
"log"
,
yscale
=
"log"
)
plot
(
r_rs
,
np
.
ones
(
np
.
size
(
r
)),
'--'
,
lw
=
1.4
,
color
=
colors
[
0
])
plot
(
r_rs
,
eta_short_gadget2
,
'-'
,
lw
=
1.4
,
color
=
colors
[
2
])
plot
(
r_rs
,
eta_short_swift
,
'-'
,
lw
=
1.4
,
color
=
colors
[
3
])
plot
(
r_rs
,
np
.
ones
(
np
.
size
(
r
))
*
0.01
,
'k:'
,
alpha
=
0.5
,
lw
=
0.5
)
plot
([
1.
,
1.
],
[
-
1e5
,
1e5
],
'k-'
,
alpha
=
0.5
,
lw
=
0.5
)
yticks
([
1e-2
,
1e-1
,
1
],
[
"$0.01$"
,
"$0.1$"
,
"$1$"
])
xlim
(
1.1
*
r_min
/
r_s
,
0.9
*
r_max
/
r_s
)
ylim
(
3e-3
,
1.5
)
#ylabel("$\\eta_s(r)$", labelpad=-3)
ylabel
(
"$|
\\
mathbf{f}_s(r)|
\\
times r^2$"
,
labelpad
=-
2
)
# 1 - Correction
subplot
(
313
,
xscale
=
"log"
,
yscale
=
"log"
)
plot
(
r_rs
,
1.
-
eta_short_gadget2
,
'-'
,
lw
=
1.4
,
color
=
colors
[
2
])
plot
(
r_rs
,
1.
-
eta_short_swift
,
'-'
,
lw
=
1.4
,
color
=
colors
[
3
])
plot
([
1.
,
1.
],
[
1e-5
,
1e5
],
'k-'
,
alpha
=
0.5
,
lw
=
0.5
)
plot
(
r_rs
,
np
.
ones
(
np
.
size
(
r
)),
'k:'
,
alpha
=
0.5
,
lw
=
0.5
)
plot
(
r_rs
,
np
.
ones
(
np
.
size
(
r
))
*
0.01
,
'k:'
,
alpha
=
0.5
,
lw
=
0.5
)
xlim
(
1.1
*
r_min
/
r_s
,
0.9
*
r_max
/
r_s
)
ylim
(
3e-3
,
1.5
)
#ylabel("$1 - \\eta_s(r)$", labelpad=-2)
ylabel
(
"$1 - |
\\
mathbf{f}_s(r)|
\\
times r^2$"
,
labelpad
=-
3
)
yticks
([
1e-2
,
1e-1
,
1
],
[
"$0.01$"
,
"$0.1$"
,
"$1$"
])
xlabel
(
"$r / r_s$"
,
labelpad
=-
3
)
savefig
(
"force_short.pdf"
)
##################################################################################################
figure
()
subplot
(
311
,
xscale
=
"log"
,
yscale
=
"log"
)
# Potential
plot
(
k_rs
,
phit_newton
,
'--'
,
lw
=
1.4
,
label
=
"${
\\
rm Newtonian}$"
,
color
=
colors
[
0
])
plot
(
k_rs
,
phit_long_gadget2
,
'-'
,
lw
=
1.4
,
label
=
"${
\\
rm Gadget}$"
,
color
=
colors
[
2
])
plot
(
k_rs
,
phit_long_swift
,
'-'
,
lw
=
1.4
,
label
=
"${
\\
rm SWIFT}$"
,
color
=
colors
[
3
])
plot
(
k_rs
,
-
phit_long_swift
,
':'
,
lw
=
1.4
,
color
=
colors
[
3
])
plot
([
1.
,
1.
],
[
1e-5
,
1e5
],
'k-'
,
alpha
=
0.5
,
lw
=
0.5
)
legend
(
loc
=
"lower left"
,
frameon
=
True
,
handletextpad
=
0.1
,
handlelength
=
3.2
,
fontsize
=
8
)
xlim
(
1.1
*
r_min
/
r_s
,
0.9
*
r_max
/
r_s
)
ylim
(
1.1
/
r_max
**
2
,
0.9
/
r_min
**
2
)
ylabel
(
"$
\\
tilde{
\\
varphi_l}(k)$"
,
labelpad
=-
3
)
yticks
([
1e-4
,
1e-2
,
1e0
,
1e2
],
[
"$10^{-4}$"
,
"$10^{-2}$"
,
"$10^{0}$"
,
"$10^{2}$"
])
subplot
(
312
,
xscale
=
"log"
,
yscale
=
"log"
)
# Potential normalized
plot
(
k_rs
,
phit_newton
*
k
**
2
,
'--'
,
lw
=
1.4
,
label
=
"${
\\
rm Newtonian}$"
,
color
=
colors
[
0
])
plot
(
k_rs
,
phit_long_gadget2
*
k
**
2
,
'-'
,
lw
=
1.4
,
label
=
"${
\\
rm Gadget}$"
,
color
=
colors
[
2
])
plot
(
k_rs
,
phit_long_swift
*
k
**
2
,
'-'
,
lw
=
1.4
,
label
=
"${
\\
rm SWIFT}$"
,
color
=
colors
[
3
])
plot
([
1.
,
1.
],
[
1e-5
,
1e5
],
'k-'
,
alpha
=
0.5
,
lw
=
0.5
)
plot
(
r_rs
,
np
.
ones
(
np
.
size
(
r
))
*
0.01
,
'k:'
,
alpha
=
0.5
,
lw
=
0.5
)
xlim
(
1.1
*
r_min
/
r_s
,
0.9
*
r_max
/
r_s
)
ylim
(
3e-3
,
1.5
)
ylabel
(
"$k^2
\\
times
\\
tilde{
\\
varphi_l}(k)$"
,
labelpad
=-
3
)
yticks
([
1e-2
,
1e-1
,
1
],
[
"$0.01$"
,
"$0.1$"
,
"$1$"
])
subplot
(
313
,
xscale
=
"log"
,
yscale
=
"log"
)
plot
(
k_rs
,
1.
-
phit_long_gadget2
*
k
**
2
,
'-'
,
lw
=
1.4
,
label
=
"${
\\
rm Gadget}$"
,
color
=
colors
[
2
])
plot
(
k_rs
,
1.
-
phit_long_swift
*
k
**
2
,
'-'
,
lw
=
1.4
,
label
=
"${
\\
rm SWIFT}$"
,
color
=
colors
[
3
])
plot
([
1.
,
1.
],
[
1e-5
,
1e5
],
'k-'
,
alpha
=
0.5
,
lw
=
0.5
)
plot
(
r_rs
,
np
.
ones
(
np
.
size
(
r
)),
'k:'
,
alpha
=
0.5
,
lw
=
0.5
)
plot
(
r_rs
,
np
.
ones
(
np
.
size
(
r
))
*
0.01
,
'k:'
,
alpha
=
0.5
,
lw
=
0.5
)
xlim
(
1.1
*
r_min
/
r_s
,
0.9
*
r_max
/
r_s
)
ylim
(
3e-3
,
1.5
)
ylabel
(
"$1 - k^2
\\
times
\\
tilde{
\\
varphi_l}(k)$"
,
labelpad
=-
3
)
yticks
([
1e-2
,
1e-1
,
1
],
[
"$0.01$"
,
"$0.1$"
,
"$1$"
])
xlabel
(
"$k
\\
times r_s$"
,
labelpad
=
0
)
savefig
(
"potential_long.pdf"
)
theory/Multipoles/plot_potential.py
View file @
3196f542
...
...
@@ -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
(
"$
\\
var
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
)
...
...
@@ -163,6 +163,6 @@ xticks([0., 0.5, 1., 1.5, 2., 2.5], ["$%.1f$"%(0./epsilon), "", "$%.1f$"%(1./eps
xlabel
(
"$r/H$"
,
labelpad
=-
7
)
ylim
(
0
,
0.95
)
ylabel
(
"$|
\\
overrightarrow{
\\
nabla}
\\
phi(r)|$"
,
labelpad
=
0
)
ylabel
(
"$|
\\
overrightarrow{
\\
nabla}
\\
var
phi(r)|$"
,
labelpad
=
0
)
savefig
(
"potential.pdf"
)
theory/Multipoles/potential_derivatives.tex
View file @
3196f542
...
...
@@ -7,7 +7,7 @@ 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}
\mathsf
{
D
}_{
000
}
(
\mathbf
{
r
}
) =
\phi
(
\mathbf
{
r
}
,H) =
\mathsf
{
D
}_{
000
}
(
\mathbf
{
r
}
) =
\
var
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,
...
...
@@ -19,7 +19,7 @@ we can construct the higher order terms by successively applying the
relevant ones here split by order.
\begin{align}
\mathsf
{
D
}_{
100
}
(
\mathbf
{
r
}
) =
\frac
{
\partial
}{
\partial
r
_
x
}
\phi
(
\mathbf
{
r
}
,H) =
\mathsf
{
D
}_{
100
}
(
\mathbf
{
r
}
) =
\frac
{
\partial
}{
\partial
r
_
x
}
\
var
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,
...
...
@@ -30,7 +30,7 @@ relevant ones here split by order.
\noindent\rule
{
6cm
}{
0.4pt
}
\begin{align}
\mathsf
{
D
}_{
200
}
(
\mathbf
{
r
}
) =
\frac
{
\partial
^
2
}{
\partial
r
_
x
^
2
}
\phi
(
\mathbf
{
r
}
,H) =
\mathsf
{
D
}_{
200
}
(
\mathbf
{
r
}
) =
\frac
{
\partial
^
2
}{
\partial
r
_
x
^
2
}
\
var
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,
\\
...
...
@@ -40,7 +40,7 @@ relevant ones here split by order.
\end{align}
\begin{align}
\mathsf
{
D
}_{
110
}
(
\mathbf
{
r
}
) =
\frac
{
\partial
^
2
}{
\partial
r
_
x
\partial
r
_
y
}
\phi
(
\mathbf
{
r
}
,H) =
\mathsf
{
D
}_{
110
}
(
\mathbf
{
r
}
) =
\frac
{
\partial
^
2
}{
\partial
r
_
x
\partial
r
_
y
}
\
var
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,
...
...
@@ -51,7 +51,7 @@ relevant ones here split by order.
\noindent\rule
{
6cm
}{
0.4pt
}
\begin{align}
\mathsf
{
D
}_{
300
}
(
\mathbf
{
r
}
) =
\frac
{
\partial
^
3
}{
\partial
r
_
x
^
3
}
\phi
(
\mathbf
{
r
}
,H) =
\mathsf
{
D
}_{
300
}
(
\mathbf
{
r
}
) =
\frac
{
\partial
^
3
}{
\partial
r
_
x
^
3
}
\
var
phi
(
\mathbf
{
r
}
,H) =
\left\lbrace
\begin{array}
{
rcl
}
-
\frac
{
r
_
x
^
3
}{
H
^
7
}
\left
(315u - 720 + 420u
^{
-1
}
\right
) +
\frac
{
3r
_
x
}{
H
^
5
}
\left
(-105u
^
3+360u
^
2-420u+168
\right
)
&
\mbox
{
if
}
&
u < 1,
\\
...
...
@@ -61,7 +61,7 @@ relevant ones here split by order.
\end{align}
\begin{align}
\mathsf
{
D
}_{
210
}
(
\mathbf
{
r
}
) =
\frac
{
\partial
^
3
}{
\partial
r
_
x
^
3
}
\phi
(
\mathbf
{
r
}
,H) =
\mathsf
{
D
}_{
210
}
(
\mathbf
{
r
}
) =
\frac
{
\partial
^
3
}{
\partial
r
_
x
^
3
}
\
var
phi
(
\mathbf
{
r
}
,H) =
\left\lbrace
\begin{array}
{
rcl
}
-
\frac
{
r
_
x
^
2r
_
y
}{
H
^
7
}
\left
(315u - 720 + 420u
^{
-1
}
\right
) +
\frac
{
r
_
y
}{
H
^
5
}
\left
(-105u
^
3+360u
^
2-420u+168
\right
)
&
\mbox
{
if
}
&
u < 1,
\\
...
...
@@ -72,7 +72,7 @@ relevant ones here split by order.
\begin{align}
\mathsf
{
D
}_{
111
}
(
\mathbf
{
r
}
) =
\frac
{
\partial
^
3
}{
\partial
r
_
x
\partial
r
_
y
\partial
r
_
z
}
\phi
(
\mathbf
{
r
}
,H) =
\mathsf
{
D
}_{
111
}
(
\mathbf
{
r
}
) =
\frac
{
\partial
^
3
}{
\partial
r
_
x
\partial
r
_
y
\partial
r
_
z
}
\
var
phi
(
\mathbf
{
r
}
,H) =
\left\lbrace
\begin{array}
{
rcl
}
-
\frac
{
r
_
xr
_
yr
_
z
}{
H
^
7
}
\left
(315u - 720 + 420u
^{
-1
}
\right
)
&
\mbox
{
if
}
&
u < 1,
\\
-15
\frac
{
r
_
xr
_
yr
_
z
}{
r
^
7
}
&
\mbox
{
if
}
&
u
\geq
1,
...
...
theory/Multipoles/potential_softening.tex
View file @
3196f542
...
...
@@ -10,9 +10,10 @@ $\epsilon$. Instead of the commonly used spline kernel of
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
architectures; it also does not require any divisions to evaluate the
softened forces. 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
\\
...
...
@@ -22,9 +23,9 @@ W(r,H) &= \frac{21}{2\pi H^3} \times \nonumber \\
\end{array}
\right
.
\end{align}