diff --git a/theory/kernel/spline.nb b/theory/kernel/spline_3.nb
similarity index 97%
rename from theory/kernel/spline.nb
rename to theory/kernel/spline_3.nb
index 1f7bed9f69cf5abf038de5220d92e0df0db364a0..d59c7f43846fd6217c2b98e193e410f6b6268cc5 100644
--- a/theory/kernel/spline.nb
+++ b/theory/kernel/spline_3.nb
@@ -10,10 +10,10 @@
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[ 157, 7]
-NotebookDataLength[ 36914, 860]
-NotebookOptionsPosition[ 35529, 807]
-NotebookOutlinePosition[ 35868, 822]
-CellTagsIndexPosition[ 35825, 819]
+NotebookDataLength[ 36970, 862]
+NotebookOptionsPosition[ 35595, 809]
+NotebookOutlinePosition[ 35934, 824]
+CellTagsIndexPosition[ 35891, 821]
WindowFrame->Normal*)
(* Beginning of Notebook Content *)
@@ -681,22 +681,24 @@ Cell[BoxData[
RowBox[{"N", "[",
RowBox[{"Expand", "[",
RowBox[{
- FractionBox["3", "4"], " ",
+ RowBox[{"-",
+ FractionBox["3", "4"]}], " ",
RowBox[{
SuperscriptBox[
RowBox[{"(",
RowBox[{"2", "-", "q"}], ")"}], "2"], "/", "Pi"}]}], "]"}],
"]"}]], "Input",
CellChangeTimes->{{3.560160709698295*^9, 3.560160723505558*^9}, {
- 3.560161185019305*^9, 3.560161189166279*^9}}],
+ 3.560161185019305*^9, 3.560161189166279*^9}, 3.560237508328278*^9}],
Cell[BoxData[
- RowBox[{"0.954929658551372`", "\[VeryThinSpace]", "-",
- RowBox[{"0.954929658551372`", " ", "q"}], "+",
+ RowBox[{
+ RowBox[{"-", "0.954929658551372`"}], "+",
+ RowBox[{"0.954929658551372`", " ", "q"}], "-",
RowBox[{"0.238732414637843`", " ",
SuperscriptBox["q", "2"]}]}]], "Output",
CellChangeTimes->{{3.560160720336149*^9, 3.560160724559553*^9},
- 3.5601612038241377`*^9}]
+ 3.5601612038241377`*^9, 3.5602375092839746`*^9}]
}, Open ]],
Cell[CellGroupData[{
@@ -806,7 +808,7 @@ Cell[BoxData[
}, Open ]]
},
WindowSize->{740, 867},
-WindowMargins->{{Automatic, -1398}, {57, Automatic}},
+WindowMargins->{{Automatic, -1324}, {Automatic, 61}},
FrontEndVersion->"8.0 for Linux x86 (64-bit) (November 7, 2010)",
StyleDefinitions->"Default.nb"
]
@@ -843,7 +845,7 @@ Cell[12697, 311, 265, 7, 30, "Input"],
Cell[12965, 320, 7987, 137, 221, "Output"]
}, Open ]],
Cell[CellGroupData[{
-Cell[20989, 462, 414, 10, 30, InheritFromParent],
+Cell[20989, 462, 414, 10, 30, "Input"],
Cell[21406, 474, 9020, 153, 223, "Output"]
}, Open ]],
Cell[CellGroupData[{
@@ -851,16 +853,16 @@ Cell[30463, 632, 247, 5, 30, "Input"],
Cell[30713, 639, 1094, 35, 65, "Output"]
}, Open ]],
Cell[CellGroupData[{
-Cell[31844, 679, 360, 11, 54, "Input"],
-Cell[32207, 692, 296, 6, 30, "Output"]
+Cell[31844, 679, 404, 12, 54, "Input"],
+Cell[32251, 693, 318, 7, 30, "Output"]
}, Open ]],
Cell[CellGroupData[{
-Cell[32540, 703, 541, 17, 54, "Input"],
-Cell[33084, 722, 238, 6, 30, "Output"]
+Cell[32606, 705, 541, 17, 54, "Input"],
+Cell[33150, 724, 238, 6, 30, "Output"]
}, Open ]],
Cell[CellGroupData[{
-Cell[33359, 733, 79, 2, 30, "Input"],
-Cell[33441, 737, 2072, 67, 118, "Output"]
+Cell[33425, 735, 79, 2, 30, "Input"],
+Cell[33507, 739, 2072, 67, 118, "Output"]
}, Open ]]
}
]
diff --git a/theory/latex/sph.tex b/theory/latex/sph.tex
index 5ac8b6fd33f950807d4305d9a0fdc81d432f81bb..5b24fe891840f46c83fd8034792e12c98b99e5ef 100755
--- a/theory/latex/sph.tex
+++ b/theory/latex/sph.tex
@@ -34,7 +34,7 @@ Every particle contains the following information:
\end{tabular}
\end{table}
-For optimisation purposes, any function of these qunatities could be stored. For instance, $1/h$ instead of $h$ or
+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.
\section{Kernel function}
@@ -70,7 +70,7 @@ function that goes to $0$ if $x > \zeta h$. \\
Coming back to the simplest case, the derivatives of the kernel function are given by:
\begin{eqnarray*}
- \vec\nabla W(\vec{x},h) &=& \frac{1}{h^4}f'\left(\frac{|\vec{x}|}{h}\right) \frac{\vec{x}}{|\vec{x}|} \\
+ \vec\nabla_x W(\vec{x},h) &=& \frac{1}{h^4}f'\left(\frac{|\vec{x}|}{h}\right) \frac{\vec{x}}{|\vec{x}|} \\
\frac{\partial W(\vec{x},h)}{\partial h} &=&- \frac{1}{h^4}\left[3f\left(\frac{|\vec{x}|}{h}\right) +
\frac{|\vec{x}|}{h}f'\left(\frac{|\vec{x}|}{h}\right)\right]
\end{eqnarray*}
@@ -87,7 +87,7 @@ with
\right. \\
&=&\left\lbrace \begin{array}{rcl}
-0.95493 q + 0.716197 q^2& \mbox{if} & 0 \leq q < 1 \\
- 0.95493 -0.95493 q+0.238732 q^2 & \mbox{if} & 1 \leq q < 2 \\
+ -0.95493+0.95493 q-0.238732 q^2 & \mbox{if} & 1 \leq q < 2 \\
0 & \mbox{if} & q \geq 2 \\
\end{array}
\right.
@@ -99,7 +99,7 @@ change this without really touching the code.
Having a compilation option somewhere which activates a given form of $f$ and changes the value of $\zeta$ accordingly
would be great.
-\section{SPH equations}
+\section{First SPH loop (density)}
In the first loop of the algorithm, the secondary quantities of particle $i$ are computed from the primary ones in the
following way:
@@ -109,9 +109,25 @@ following way:
h_i &=& \eta \left(\frac{m_i}{\rho_i} \right)^{1/3}
\end{eqnarray}
-where $\eta \approx 1.2$ is a constant. These two equations can be solved iterativelly using a Newton-Raphson or
-bissection scheme. In practice, the loop is performed over all particles $j$ which are at a distance $r_{ij}<\zeta
-h$ from the particle of interest. One has to iterate those two equations until their outcomes are stable. \\
+where $\eta \approx 1.2$ is a constant. These two equations can be solved iteratively using a Newton-Raphson or
+bisection scheme. In practice, the loop is performed over all particles $j$ which are at a distance $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
+
+\begin{equation}
+ N_{ngb} = \frac{4}{3}\pi \left(\zeta h\right)^3 \sum_j W(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
+
+\begin{equation}
+ N_{ngb} = \frac{4}{3}\pi\left(\zeta \eta\right)^3
+\end{equation}
+
+We usually use $N_{ngb} = 48$, which corresponds to a sub-optimal value of $\eta=1.127$. The optimal value should be
+$N_{ngb}=57.9$ but this is computationally more expensive and the improvement over $N_{ngb}=48$ is not obvious.\\
+
To increase the convergence rate, one can use the derivative of the density with respect to the smoothing length in the
Newton iterations:
@@ -127,7 +143,13 @@ This term is given by
\Omega_i = 1 + \frac{h_i}{3\rho_i}\sum_b m_b\frac{\partial W(r_{ij},h_i)}{\partial h}
\end{equation}
-Once those quantities have been obtained, the force estiamtion loop can be started.
+This concludes the first SPH loop in the standard implementation. More complicated quantities such as
+$\vec\nabla\times\vec v_i$ or $\vec\nabla\cdot\vec v_i$ are sometimes computed here if they are needed in the force
+loop.
+
+\section{Second SPH loop (forces)}
+
+Once those quantities have been obtained, the force estimation loop can be started.
First, the pressure has to be evaluated evaluated using the equation of state
\begin{equation}
@@ -138,16 +160,16 @@ where $\gamma$ is the polytropic index. Usually, $\gamma = \frac{5}{3}$.
The second loop is used to compute the accelerations (tertiary quantities). The exact expressions are
\begin{eqnarray}
- \vec{a} &=& - \sum_j m_j\left[\frac{P_i}{\Omega_i\rho_i^2}\vec{\nabla} W(r_{ij}, h_i) +
-\frac{P_j}{\Omega_j\rho_j^2}\vec{\nabla}
-W(r_{ij}, h_j) \right] \\
- \frac{du}{dt} &=& \frac{P_i}{\rho_i^2} \sum_j m_j (\vec{v_i}-\vec{v_j})\cdot\vec{\nabla} W(r_{ij}, h_i) \\
- \frac{dh}{dt} &=& \frac{h_i}{3}\sum_j \frac{m_j}{\rho_j} \left(\vec{v_j} - \vec{v_i} \right) \cdot\vec{\nabla}W(r_{ij},
+ \vec{a} &=& - \sum_j m_j\left[\frac{P_i}{\Omega_i\rho_i^2}\vec{\nabla_r} W(r_{ij}, h_i) +
+\frac{P_j}{\Omega_j\rho_j^2}\vec{\nabla_r}W(r_{ij}, h_j) \right] \\
+ \frac{du}{dt} &=& \frac{P_i}{\rho_i^2} \sum_j m_j (\vec{v_i}-\vec{v_j})\cdot\vec{\nabla_r} W(r_{ij}, h_i) \\
+ \frac{dh}{dt} &=& \frac{h_i}{3}\sum_j \frac{m_j}{\rho_j} \left(\vec{v_j} - \vec{v_i} \right)
+\cdot\vec{\nabla_r}W(r_{ij},
h_i)
\end{eqnarray}
In practice the loop is here performed over all pairs of particles such that $r_{ij} < \zeta h_i$ or $r_{ij} < \zeta
-h_j$. In general, the equations are more involved as they will contain terms to mimimc the effect of viscosity or
+h_j$. In general, the equations are more involved as they will contain terms to mimic the effect of viscosity or
thermal conduction. These terms are pure functions of the properties of particles $i$ and $j$ and are thus very simple
to insert once the code is stabilized.\\
@@ -165,14 +187,14 @@ The time step is then given by the Courant relation:
where the CFL parameter usually takes a value between $0.2$ and $0.3$. The integration in time can then take place. The
leapfrog integrator is usually used as it behaves well when coupled to gravity. \\
-In the case where only one global timestep is used for all particles, the minimal timestep of all particles is reduced
+In the case where only one global time step is used for all particles, the minimal time step of all particles is reduced
and used. \\
Notice that $h$ has to be recomputed through the iterative process
-presented in the previous section at every timestep. The time
+presented in the previous section at every time step. The time
derivative of the smoothing length only give a rough estimate of its
change. It only provides a good guess for the Newton-Raphson (or
-bissection) scheme.
+bisection) scheme.
\section{Conserved quantities}