Skip to content
GitLab
Explore
Sign in
Primary navigation
Search or go to…
Project
SWIFTsim
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Deploy
Releases
Model registry
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
SWIFT
SWIFTsim
Commits
dc14c25e
Commit
dc14c25e
authored
12 years ago
by
Matthieu Schaller
Browse files
Options
Downloads
Patches
Plain Diff
Corrected some typos in the theory PDF.
Former-commit-id: 7fbf7981fa52d1c13d62cdd564c46d88f167b594
parent
d5dfe1b3
No related branches found
No related tags found
No related merge requests found
Changes
1
Hide whitespace changes
Inline
Side-by-side
Showing
1 changed file
theory/latex/sph.tex
+38
-29
38 additions, 29 deletions
theory/latex/sph.tex
with
38 additions
and
29 deletions
theory/latex/sph.tex
+
38
−
29
View file @
dc14c25e
...
...
@@ -18,7 +18,8 @@ Every particle contains the following information:
\begin{table}
[h]
\centering
\begin{tabular}
{
|l|l|c|l|
}
Quantity
&
Type
&
Symbol
&
Unit
\\
\hline
\textbf
{
Quantity
}
&
\textbf
{
Type
}
&
\textbf
{
Symbol
}
&
\textbf
{
Units
}
\\
\hline
\hline
Position
&
Primary
&
$
\vec
{
x
}$
&
$
[
m
]
$
\\
Velocity
&
Primary
&$
\vec
{
v
}$
&
$
[
m
\cdot
s
^{
-
1
}
]
$
\\
...
...
@@ -41,19 +42,21 @@ Every particle contains the following information:
\end{tabular}
\end{table}
Secondary quantities are computed from the primary one in a loop
of
particle neighbors. Tertiary
ones are computed from
se
ondary ones in another loop.
\\
Secondary quantities are computed from the primary one in a loop
(density loop) over all
particle neighbors. Tertiary
ones are computed from sec
ondary ones in another loop
(force loop)
.
\\
For optimization purposes, any function of these quantities could be stored. For instance,
$
1
/
h
$
instead of
$
h
$
or
$
\frac
{
P
}{
\rho\Omega
}$
instead of
$
\Omega
$
may be options worth exploring.
\\
The four quantities in the second half of the table are used in improved state-of-the-art implementations of SPH. In a
first approximations, they can be neglected.
The four quantities in the second part of the table are used in improved state-of-the-art implementations of SPH. In a
first approximation, they can be neglected.
\\
In what follows, we will use
$
\vec
{
r
}_{
ij
}
=
\vec
{
x
_
i
}
-
\vec
{
x
_
j
}$
and
$
\hat
{
r
}_{
ij
}
=
\vec
{
r
}_{
ij
}
/
x|
\vec
{
r
}_{
ij
}
|
$
to simplify the notation.
\section
{
Kernel function
}
In what follows, we will use
$
\vec
{
r
}_{
ij
}
=
\vec
{
x
_
i
}
-
\vec
{
x
_
j
}$
and
$
\hat
{
r
}_{
ij
}
=
|
\vec
{
r
}_{
ij
}
|
$
. The kernel
function can always be decomposed as:
The kernel function can always be decomposed as:
\begin{equation}
W(
\vec
{
x
}
, h) =
\frac
{
1
}{
h
^
3
}
f
\left
(
\frac
{
|
\vec
{
x
}
|
}{
h
}
\right
)
...
...
@@ -128,14 +131,15 @@ where $\eta \approx 1.2$ is a constant. These two equations can be solved iterat
bisection scheme. In practice, the loop is performed over all particles
$
j
$
which are at a distance
$
|
\vec
{
r
}_{
ij
}
|<
\zeta
h
$
from the particle of interest. One has to iterate those two equations until their outcomes are stable.
\\
Another measure of the accuracy of
$
h
$
is
to use
the weighted number of neighbors which (in 3D) reads
Another measure of the accuracy of
$
h
$
is the weighted number of neighbors which (in 3D) reads
\begin{equation}
N
_{
ngb
}
=
\frac
{
4
}{
3
}
\pi
\left
(
\zeta
h
\right
)
^
3
\sum
_
j W(
\vec
{
r
}_{
ij
}
,h
_
i)
\end{equation}
The (magical) value of
$
N
_{
ngb
}$
is a numerical parameter and its value can be expressed as a function of the more
physically relevant parameter
$
\eta
$
. In 3D this reads
One then change
$
h
$
until an optimal value for
$
N
_{
ngb
}$
is reached. GADGET uses a bisection algorithm to do so.
The (magical) value of
$
N
_{
ngb
}$
to obtain is a numerical parameter and its value can be expressed as a function of the
more physically relevant parameter
$
\eta
$
. In 3D the relation between those quantities is
\begin{equation}
N
_{
ngb
}
=
\frac
{
4
}{
3
}
\pi\left
(
\zeta
\eta\right
)
^
3
...
...
@@ -212,7 +216,7 @@ and used. \\
Notice that
$
h
$
has to be recomputed through the iterative process
presented in the previous section at every time step. The time
derivative of the smoothing length only give a rough estimate of its
derivative of the smoothing length only give
s
a rough estimate of its
change. It only provides a good guess for the Newton-Raphson (or
bisection) scheme.
...
...
@@ -221,21 +225,21 @@ bisection) scheme.
The usual scheme uses a kick-drift-kick leap-frog integrator. A full time step of size
$
\Delta
t
$
consists of the
following sub-steps:
\\
\textbf
{
f
irst kick
}
Compute velocity and internal energy at half step.
\textbf
{
F
irst kick
}
Compute velocity and internal energy at half step.
\begin{eqnarray*}
\tilde
{
\vec
{
v
}}_
i
&
=
&
\vec
{
v
}_
i +
\textstyle\frac
{
1
}{
2
}
\Delta
t ~
\vec
{
a
}_
i
\\
\tilde
u
_
i
&
=
&
u
_
i +
\textstyle\frac
{
1
}{
2
}
\Delta
t ~
\frac
{
du
_
i
}{
dt
}
\end{eqnarray*}
\textbf
{
d
rift
}
Advance time and position by a full step.
\textbf
{
D
rift
}
Advance time and position by a full step.
\begin{eqnarray*}
t
&
\leftarrow
&
t +
\Delta
t
\\
\vec
{
x
}_
i
&
\leftarrow
&
\vec
{
x
}_
i +
\Delta
t
\tilde
{
\vec
{
v
}}_
i
\\
\end{eqnarray*}
\textbf
{
p
rediction
}
Estimate velocity, internal energy and smoothing length at full step
\textbf
{
P
rediction
}
Estimate velocity, internal energy and smoothing length at full step
\begin{eqnarray*}
\vec
{
v
}_
i
&
\leftarrow
&
\vec
{
v
}_
i +
\Delta
t
\vec
{
a
}_
i
\\
...
...
@@ -243,15 +247,16 @@ u_i &\leftarrow& u_i + \Delta t ~\frac{du_i}{dt} \\
h
_
i
&
\leftarrow
&
h
_
i +
\Delta
t ~
\frac
{
dh
_
i
}{
dt
}
\\
\end{eqnarray*}
\textbf
{
SPH loop 1
}
Compute
$
\rho
_
i
$
, correct
$
h
_
i
$
and compute
$
\Omega
_
i
$
using the first SPH loop.
\\
\textbf
{
SPH loop 1
}
Compute
$
\rho
_
i
$
, correct
$
h
_
i
$
if needed using bisection algorithm and compute
$
\Omega
_
i
$
using the
first SPH loop.
\\
\textbf
{
SPH loop 2
}
Compute
$
\vec
{
a
_
i
}$
and
$
\frac
{
d
u
_
i
}{
dt
}$
using the second SPH loop.
\\
\textbf
{
SPH loop 2
}
Compute
$
\vec
{
a
_
i
}$
,
$
\frac
{
du
_
i
}{
dt
}$
and
$
\frac
{
d
h
_
i
}{
dt
}$
using the second SPH loop.
\\
\textbf
{
Gravity
}
Compute accelerations due to gravity.
\\
\textbf
{
Cooling
}
Compute the change in internal energy due to radiative cooling.
\\
\textbf
{
Cooling
}
Compute the change in internal energy due to radiative cooling
and heating
.
\\
\textbf
{
s
econd kick
}
Compute velocity and internal energy at end of step.
\textbf
{
S
econd kick
}
Compute velocity and internal energy at end of step.
\begin{eqnarray*}
\vec
{
v
}_
i
&
=
&
\tilde
{
\vec
{
v
}}_
i +
\textstyle\frac
{
1
}{
2
}
\Delta
t ~
\vec
{
a
}_
i
\\
...
...
@@ -272,9 +277,10 @@ E &=&\sum_i m_i\left(\frac{1}{2}|\vec{v_i}|^2+u_i\right)\\
A(s)
&
=
&
\left
(
\gamma
-1
\right
)
\sum
_
i
\frac
{
u
_
i
}{
\rho
_
i
^{
\gamma
- 1
}}
\end{eqnarray}
The conservation of those quantities in the code depends on the quality of the time integrator.
\\
The conservation of those quantities in the code depends on the quality of the time integrator. The leap-frog
integrator of the previous section should preserve these quantities to machine precision.
\\
Notice that the entropic function
$
A
(
s
)
$
is not the ``physical'' entropy
$
s
$
but is related to it through a monotonic
function. It is just a convenient way to represent entropy.
function. It is just a
more
convenient way to represent entropy.
\section
{
Improved SPH equations
}
...
...
@@ -340,7 +346,8 @@ using a shock detector and then a slow decay of the viscosity with time. Followi
\end{eqnarray}
with (usually)
$
\alpha
_{
max
}
=
2
$
,
$
\alpha
_{
min
}
=
0
.
1
$
and
$
\sigma
=
0
.
1
$
. The
$
\alpha
$
term in (
\ref
{
eq:visc
}
) is then
replaced by
$
\bar\alpha
=
\frac
{
1
}{
2
}
(
\alpha
_
i
+
\alpha
_
j
)
$
. The same applies to
$
\beta
$
.
\\
replaced by
$
\bar\alpha
=
\frac
{
1
}{
2
}
(
\alpha
_
i
+
\alpha
_
j
)
$
. The same applies to
$
\beta
$
which is now,
$
\beta
=
3
\bar\alpha
$
.
\\
The divergence of the velocity field can be computed in the density loop and the exact expression is
\begin{equation}
...
...
@@ -350,14 +357,15 @@ The divergence of the velocity field can be computed in the density loop and the
\subsection
{
Thermal conductivity
}
The thermal conductivity dissipat
ing
energy at
contact
discontinuities can be model
l
ed by adding
another term to the
internal energy equation
\ref
{
eq:dudt
}
.
The thermal conductivity
which
dissipat
es
energy at discontinuities
in the energy field
can be modeled by adding
another term to the
internal energy
evolution
equation
(
\ref
{
eq:dudt
}
)
.
\begin{equation}
\frac
{
du
_
i
}{
dt
}
\stackrel
{
cond
}{
=
}
-
\sum
_
j
\alpha
_
u v
_{
sig,u
}
\left
(u
_
i - u
_
j
\right
)
\bar
{
F
}_{
ij
}
\end{equation}
This time, the signal velocity must vanish when we are dealing with a contact discontinuity. A good choice is to used
This time, the signal velocity must vanish when we are dealing with a contact discontinuity as no energy should flow
between the two regions in this case. A good choice is to use
\begin{equation}
v
_{
sig,u
}
=
\sqrt
{
\frac
{
2|P
_
i-P
_
j|
}{
\rho
_
i+
\rho
_
j
}}
...
...
@@ -365,12 +373,13 @@ This time, the signal velocity must vanish when we are dealing with a contact di
Once again, the
$
\alpha
$
-term is made to decay far from any discontinuity. In this case, the equation reads
\begin{equation}
\frac
{
d
\alpha
_{
u,i
}}{
dt
}
=
\frac
{
h
_
i
\nabla
^
2 u
_
i
}{
10
}
-
\frac
{
\alpha
_{
u,i
}
c
_
i
\sigma
}{
h
_
i
}
\end{equation}
\begin{eqnarray}
\frac
{
d
\alpha
_{
u,i
}}{
dt
}
&
=
&
\mathcal
{
S
}_
u -
\frac
{
\alpha
_{
u,i
}
c
_
i
\sigma
}{
h
_
i
}
\\
\mathcal
{
S
}_
u
&
=
&
\frac
{
h
_
i
\nabla
^
2 u
_
i
}{
10
}
\end{eqnarray}
where a
ll constants usually take the same value than in the viscosity terms. T
he laplacian of
$
u
$
can
again be computed
in the density loop and reads
where a
gain
$
\sigma
=
0
.
1
$
as in the viscosity terms. As in the velocity divergence case, t
he laplacian of
$
u
$
can
be
computed
in the density loop and reads
\begin{equation}
\nabla
^
2 u
_
i =
\frac
{
2
}{
\rho
_
i
}
\sum
_
j m
_
j
\left
(u
_
j - u
_
i
\right
)
\frac
{
F
_{
ij
}}{
|
\vec
{
r
}_{
ij
}
|
}
...
...
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
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!
Save comment
Cancel
Please
register
or
sign in
to comment