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
8b385175
Commit
8b385175
authored
Jul 31, 2017
by
Matthieu Schaller
Browse files
Added the short-range correction equations to the FMM pdf documentation.
parent
1ed63fd2
Changes
2
Hide whitespace changes
Inline
Side-by-side
theory/Multipoles/mesh_summary.tex
View file @
8b385175
\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
}
...
...
@@ -9,7 +26,14 @@
\begin{figure}
\includegraphics
[width=\columnwidth]
{
potential
_
long
.pdf
}
\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
View file @
8b385175
...
...
@@ -67,56 +67,42 @@ phi_newton = 1. / r
phit_newton
=
1.
/
k
**
2
force_newton
=
1.
/
r
**
2
def
smoothstep
(
x
):
#S_2(x)
ret
=
6
*
x
**
5
-
15
*
x
**
4
+
10
*
x
**
3
#ret = 3*x**2 - 2*x**3
ret
[
x
<
0
]
=
0.
ret
[
x
>
1
]
=
1.
return
ret
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
)
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
1.
/
(
1.
+
1.
/
my_exp
(
x
))
#return x / sqrt(1. + x**2)
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. * smoothstep(x/4. + 1./2.) - 1.
#return sigmoid(4. * x)
return
2
*
sigmoid
(
4
*
x
)
-
1
def
d_swift_corr
(
x
):
return
2
*
d_sigmoid
(
4
*
x
)
def
csch
(
x
):
# hyperbolic cosecant
return
1.
/
sinh
(
x
)
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
)
#plot(x, exp(x), '-', color=colors[0])
#plot(x, my_exp(x), '-', color=colors[1])
savefig
(
"temp.pdf"
)
#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_long_gadget2
=
exp
(
-
k
**
2
*
r_s
**
2
)
corr_short_swift
=
swift_corr
(
r
/
(
2.
*
r_s
))
#corr_long_swift = k * r_s * math.pi / (2. * sinh(0.5 * math.pi * r_s * k))
corr_long_swift
=
math
.
pi
*
k
*
r_s
*
csch
(
0.5
*
math
.
pi
*
r_s
*
k
)
/
2.
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
=
(
2.
*
r
*
d_swift_corr
(
r
/
(
2.
*
r_s
))
/
r_s
)
-
2
*
sigmoid
(
2
*
r
/
(
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
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
...
...
@@ -149,18 +135,17 @@ legend(loc="upper right", frameon=True, handletextpad=0.1, handlelength=3.2, fon
# Correction
subplot
(
312
,
xscale
=
"log"
,
yscale
=
"log"
)
#
plot(r_rs, np.ones(np.size(r)), '--', lw=1.4, color=colors[0])
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.zeros(np.size(r)), '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
)
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("$\\chi_s(r)$", labelpad=-3)
ylabel
(
"$
\\
varphi_s(r)
\\
times r$"
,
labelpad
=-
2
)
# 1 - Correction
subplot
(
313
,
xscale
=
"log"
,
yscale
=
"log"
)
...
...
@@ -168,12 +153,13 @@ 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
)
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 - \\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
)
...
...
@@ -200,16 +186,17 @@ legend(loc="upper right", frameon=True, handletextpad=0.1, handlelength=3.2, fon
# 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
)),
'k--'
,
alpha
=
0.5
,
lw
=
0.5
)
plot
(
r_rs
,
np
.
ones
(
np
.
size
(
r
))
*
0.01
,
'k--'
,
alpha
=
0.5
,
lw
=
0.5
)
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("$\\eta_s(r)$", labelpad=-3)
ylabel
(
"$|
\\
mathbf{f}_s(r)|
\\
times r^2$"
,
labelpad
=-
2
)
# 1 - Correction
subplot
(
313
,
xscale
=
"log"
,
yscale
=
"log"
)
...
...
@@ -217,12 +204,13 @@ 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
)
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 - \\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
)
...
...
@@ -231,7 +219,7 @@ savefig("force_short.pdf")
##################################################################################################
figure
()
subplot
(
2
11
,
xscale
=
"log"
,
yscale
=
"log"
)
subplot
(
3
11
,
xscale
=
"log"
,
yscale
=
"log"
)
# Potential
plot
(
k_rs
,
phit_newton
,
'--'
,
lw
=
1.4
,
label
=
"${
\\
rm Newtonian}$"
,
color
=
colors
[
0
])
...
...
@@ -245,22 +233,35 @@ legend(loc="lower left", frameon=True, handletextpad=0.1, handlelength=3.2, font
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
(
212
,
xscale
=
"log"
,
yscale
=
"log"
)
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
(
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
)
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"
)
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment