Added an explicit derivation of the kernels used in SWIFT with all the...

Added an explicit derivation of the kernels used in SWIFT with all the relevant 'magical' constants.


Former-commit-id: 45efdc473a2fbae9e1d871cc1d4d6335f9d75730

theory/kernel/kernel.tex

+\documentclass[a4paper,10pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage{amsmath}
+
+%opening
+\title{SPH kernels in SWIFT}
+\author{Matthieu Schaller}
+
+\begin{document}
+
+\maketitle
+
+\section{General Definitions}
+
+The smoothing kernels used in SPH are almost always isotropic and can hence be written in 3D as
+
+\begin{equation}
+ W(\vec{x},h) = \frac{1}{h^3}f\left(\frac{|\vec{x}|}{h}\right),
+\end{equation}
+
+where $f(q)$ is a dimensionless function, usually a low-order polynomial, normalized to unity. For computational 
+reasons, this kernel 
+usually has a finite support of radius $H$. In other words,
+
+\begin{equation}
+ W(\vec{x},h) = 0\quad \forall\quad |\vec{x}| > H.
+\end{equation}
+ One can then define the weighted number of neighbours within $H$ as
+
+\begin{equation}
+ N_{ngb} = \frac{4}{3}\pi H^3 \sum_j W(\vec{x}_i - \vec{x}_j,h).
+\end{equation}
+
+The value of $N_{ngb}$ is often used in the codes to find the smoothing length of each particle via Newton iterations 
+or a bissection algorithm. $H$ is defined as \emph{the smoothing length} in the GADGET code. This definition is useful 
+for implementation reasons but does not really correspond to a true physical quantity. \\
+The main question is the definition of the smoothing length. The function $W(\vec{x},h)$ is invariant under the 
+rescaling $h\rightarrow \alpha h,~f(q)\rightarrow\alpha^{-3}f(\alpha q)$, which makes the definition of $h$ difficult. 
+This ambiguity is present in the litterature with authors using different definition of the \emph{physical} smoothing 
+length, $h=\frac{1}{2}H$ or $h=H$ for instance. \\
+A more physically motivated estimate is the standard deviation of the kernel:
+
+\begin{equation}
+ \sigma^2 = \frac{1}{3} \int \vec{x}^2~W(\vec{x},h)~d^3\vec{x}
+\end{equation}
+ 
+which then allows us to set $h=2\sigma$. This definition of the smoothing length is more physical as one can 
+demonstrate that the reconstruction of any smooth field $A(\vec{x})$ using interpolation of particles at the point 
+$\vec{x}_i$ can be expanded as
+
+\begin{equation}
+A_i  \approx A(\vec{x}_i) + \frac{1}{2}\sigma^2 \nabla^2A(\vec{x}_i) + \mathcal{O}\left(\sigma^4\right).
+\end{equation}
+
+The quantity $H/\sigma$ is independant of the choice of $h$ made and is purely a functional of $f(q)$. The number of 
+neighbours (used in the code to construct the neighborhood of a given particle) can then be expressed as a function of 
+this \emph{physical} $h$ (or $\sigma$). Or to relate it even more to the particle distribution, we can write 
+$h=\eta\Delta x$, with $\Delta x$ the mean inter-particle separation:
+
+\begin{equation}
+ N_{ngb} = \frac{4}{3}\pi \left(\frac{1}{2}\eta\frac{H}{\sigma}\right)^3 = \frac{4}{3}\pi 
+\left(\eta\zeta\right)^3
+\end{equation}
+
+This definition of the number of neighbours only depends on $f(q)$ (via $\zeta$) and on the mean inter-particle 
+separation. The problem is then fully specified by specifying a form for $f(q)$ and $\eta$. \\
+Experiments suggest that $\eta \approx 1.2 - 1.3$ is a reasonnable choice. The bigger $\eta$, the better the smoothing 
+and hence the better the reconstruction of the field. This, however, comes at a higher computational cost as more 
+interactions between neighbours will have to be computed. Also, spline kernels become instable when $\eta>1.5$. 
+
+\section{Kernels available in SWIFT}
+
+The different kernels available are listed below.
+\paragraph{Cubic Spline Kernel}
+\begin{equation*}
+ f(q) = \frac{1}{\pi}\left\lbrace \begin{array}{rcl}
+                      \frac{3}{4}q^3 - 15q^2 + 1 & \mbox{if} & 0 \leq q < 1 \\
+		      -\frac{1}{4}q^3 + \frac{3}{2}q^2-3q+2 & \mbox{if} & 1 \leq q < 2 \\
+		      0 & \mbox{if} & q \geq 2 \\
+                     \end{array}\right.
+\end{equation*}
+with $\zeta = \frac{1}{2}\sqrt{\frac{40}{3}} \approx 1.825742$. Thus, for a resolution of $\eta = 1.235$, this kernel 
+uses $N_{ngb} \approx 48$. The code uses $h = \frac{1}{2}H = \zeta \sigma$ internally.
+
+\paragraph{Quartic Spline Kernel}
+\begin{equation*}
+ f(q) = \frac{1}{20\pi}\left\lbrace \begin{array}{rcl}
+                      6q^4 - 15q^2 + \frac{115}{8} & \mbox{if} & 0 \leq q < \frac{1}{2} \\
+		      -4q^4 + 20q^3-30q^2 + 5q + \frac{55}{4} & \mbox{if} & \frac{1}{2} \leq q < \frac{3}{2} \\
+		      q^4-10q^3+\frac{75}{2}q^2-\frac{125}{2}q+\frac{625}{16} & \mbox{if} & \frac{3}{2} \leq q < 
+\frac{5}{2} \\
+		      0 & \mbox{if} & q \geq \frac{5}{2} \\
+                     \end{array}\right.
+\end{equation*}
+with $\zeta = \frac{1}{2}\sqrt{\frac{375}{23}} \approx 2.018932$. Thus, for a resolution of $\eta = 1.235$, this kernel 
+uses $N_{ngb} \approx 64.9$. The code uses $h = \frac{2}{5}H =\frac{4}{5}\zeta \sigma$ internally.
+
+\paragraph{Quintic Spline Kernel}
+\begin{equation*}
+ f(q) = \frac{1}{120\pi}\left\lbrace \begin{array}{rcl}
+                      -10q^5 + 30q^4 - 60q^2 + 66 & \mbox{if} & 0 \leq q < 1 \\
+		      5q^5 - 45q^4 + 150q^3 - 210q^2 + 75q + 51 & \mbox{if} & 1 \leq q < 2 \\
+		      -q^5 + 15q^4 - 90q^3 + 270q^2 - 405q + 243 & \mbox{if} & 2 \leq q < 3 \\
+		      0 & \mbox{if} & q \geq 3 \\
+                     \end{array}\right.
+\end{equation*}
+with $\zeta = \frac{1}{2}\sqrt{\frac{135}{7}} \approx 2.195775$. Thus, for a resolution of $\eta = 1.235$, this kernel 
+uses $N_{ngb} \approx 83.5$. The code uses $h = \frac{1}{3}H = \frac{2}{3}\zeta \sigma$ internally.
+
+\paragraph{Wendland $C2$ Kernel}
+\begin{equation*}
+ f(q) = \frac{21}{2\pi}\left\lbrace \begin{array}{rcl}
+                      4 q^5-15 q^4+20 q^3-10 q^2+1 & \mbox{if} & 0 \leq q < 1 \\
+		      0 & \mbox{if} & q \geq 1 \\
+                     \end{array}\right.
+\end{equation*}
+ with $\zeta = \frac{1}{2}\sqrt{15} \approx 1.93649$. Thus, for a resolution of $\eta = 1.235$, this kernel 
+uses $N_{ngb} \approx 57.3$. The code uses $h = H = 2\zeta \sigma$ internally.
+
+
+\paragraph{Wendland $C4$ Kernel}
+\begin{equation*}
+ f(q) = \frac{495}{32\pi}\left\lbrace \begin{array}{rcl}
+                      \frac{35}{3} q^8-64 q^7+ 140 q^6-\frac{448}{3} q^5+70 q^4-\frac{28}{3} q^2+1 & \mbox{if} & 0 
+\leq q < 1 \\
+		      0 & \mbox{if} & q \geq 1 \\
+                     \end{array}\right.
+\end{equation*}
+ with $\zeta = \frac{1}{2}\sqrt{\frac{39}{2}} \approx 2.20794$. Thus, for a resolution of $\eta = 1.235$, this kernel 
+uses $N_{ngb} \approx 84.9$. The code uses $h = H = 2\zeta \sigma$ internally.
+
+\paragraph{Wendland $C6$ Kernel}
+\begin{equation*}
+ f(q) = \frac{1365}{64\pi}\left\lbrace \begin{array}{rcl}
+                      32 q^{11}-231 q^{10}+704 q^9-1155 q^8+1056 q^7-462 q^6+66 q^4-11 q^2+1 & \mbox{if} & 0 
+\leq q < 1 \\
+		      0 & \mbox{if} & q \geq 1 \\
+                     \end{array}\right.
+\end{equation*}
+ with $\zeta = \frac{1}{2}\sqrt{24} \approx 2.44949$. Thus, for a resolution of $\eta = 1.235$, this kernel 
+uses $N_{ngb} \approx 116$. The code uses $h = H = 2\zeta \sigma$ internally.
+
+\section{Kernel Derivatives}
+
+The derivatives of the kernel function have relatively simple expressions:
+
+\begin{eqnarray*}
+ \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*}
+
+Note that for all the kernels listed above, $f'(0) = 0$. 
+
+\end{document}