Merge branch 'master' into StefanCooling

.gitignore

@@ -83,6 +83,7 @@ theory/latex/swift.pdf
 theory/SPH/kernel/kernels.pdf
 theory/SPH/kernel/kernel_derivatives.pdf
 theory/SPH/kernel/kernel_definitions.pdf
+theory/SPH/flavours/sph_flavours.pdf
 theory/paper_pasc/pasc_paper.pdf
 
 m4/libtool.m4

examples/main.c

@@ -268,7 +268,7 @@ int main(int argc, char *argv[]) {
     message("CPU frequency used for tick conversion: %llu Hz", cpufreq);
   }
 
-  /* Report host name(s). */
+/* Report host name(s). */
 #ifdef WITH_MPI
   if (myrank == 0 || verbose > 1) {
     message("Rank %d running on: %s", myrank, hostname());

theory/SPH/flavours/bibliography.bib

+@ARTICLE{Price2012,
+   author = {{Price}, D.~J.},
+    title = "{Smoothed particle hydrodynamics and magnetohydrodynamics}",
+  journal = {Journal of Computational Physics},
+archivePrefix = "arXiv",
+   eprint = {1012.1885},
+ primaryClass = "astro-ph.IM",
+     year = 2012,
+    month = feb,
+   volume = 231,
+    pages = {759-794},
+      doi = {10.1016/j.jcp.2010.12.011},
+   adsurl = {http://adsabs.harvard.edu/abs/2012JCoPh.231..759P},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+
+@ARTICLE{Springel2005,
+   author = {{Springel}, V.},
+    title = "{The cosmological simulation code GADGET-2}",
+  journal = {\mnras},
+   eprint = {astro-ph/0505010},
+ keywords = {methods: numerical, galaxies: interactions, dark matter},
+     year = 2005,
+    month = dec,
+   volume = 364,
+    pages = {1105-1134},
+      doi = {10.1111/j.1365-2966.2005.09655.x},
+   adsurl = {http://adsabs.harvard.edu/abs/2005MNRAS.364.1105S},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+
+
+@ARTICLE{Hopkins2013,
+   author = {{Hopkins}, P.~F.},
+    title = "{A general class of Lagrangian smoothed particle hydrodynamics methods and implications for fluid mixing problems}",
+  journal = {\mnras},
+archivePrefix = "arXiv",
+   eprint = {1206.5006},
+ primaryClass = "astro-ph.IM",
+ keywords = {hydrodynamics, instabilities, turbulence, methods: numerical, cosmology: theory},
+     year = 2013,
+    month = feb,
+   volume = 428,
+    pages = {2840-2856},
+      doi = {10.1093/mnras/sts210},
+   adsurl = {http://adsabs.harvard.edu/abs/2013MNRAS.428.2840H},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@ARTICLE{Springel2002,
+   author = {{Springel}, V. and {Hernquist}, L.},
+    title = "{Cosmological smoothed particle hydrodynamics simulations: the entropy equation}",
+  journal = {\mnras},
+   eprint = {astro-ph/0111016},
+ keywords = {methods: numerical, galaxies: evolution, galaxies: starburst, methods: numerical, galaxies: evolution, galaxies: starburst},
+     year = 2002,
+    month = jul,
+   volume = 333,
+    pages = {649-664},
+      doi = {10.1046/j.1365-8711.2002.05445.x},
+   adsurl = {http://adsabs.harvard.edu/abs/2002MNRAS.333..649S},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@ARTICLE{Balsara1995,
+   author = {{Balsara}, D.~S.},
+    title = "{von Neumann stability analysis of smooth particle hydrodynamics--suggestions for optimal algorithms}",
+  journal = {Journal of Computational Physics},
+     year = 1995,
+   volume = 121,
+    pages = {357-372},
+      doi = {10.1016/S0021-9991(95)90221-X},
+   adsurl = {http://adsabs.harvard.edu/abs/1995JCoPh.121..357B},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+
+
+
+

theory/SPH/flavours/run.sh

+#!/bin/bash
+pdflatex sph_flavours.tex
+bibtex sph_flavours.aux 
+pdflatex sph_flavours.tex
+pdflatex sph_flavours.tex

theory/SPH/flavours/sph_flavours.tex

+\documentclass[fleqn, usenatbib, useAMS, a4paper]{mnras}
+\usepackage{graphicx}
+\usepackage{amsmath,paralist,xcolor,xspace,amssymb}
+\usepackage{times}
+\usepackage[super]{nth}
+
+\newcommand{\swift}{{\sc swift}\xspace}
+\newcommand{\gadget}{{\sc gadget}\xspace}
+\newcommand{\dd}[2]{\frac{\partial #1}{\partial #2}}
+\renewcommand{\vec}[1]{{\mathbf{#1}}}
+\newcommand{\Wij}{\overline{\nabla_xW_{ij}}}
+
+%\setlength{\mathindent}{0pt}
+
+%opening
+\title{SPH flavours in SWIFT}
+\author{Matthieu Schaller}
+
+\begin{document}
+
+\maketitle
+
+For vector quantities, we define $\vec{a}_{ij} \equiv (\vec{a}_i -
+\vec{a}_j)$. Barred quantities ($\bar a_{ij}$) correspond to the
+average over two particles of a quantity: $\bar a_{ij} \equiv
+\frac{1}{2}(a_i + a_j)$. To simplify notations, we also define the
+vector quantity $\Wij \equiv \frac{1}{2}\left(W(\vec{x}_{ij}, h_i) +
+\nabla_x W(\vec{x}_{ij},h_j)\right)$.
+
+\subsection{Minimal SPH}
+\label{sec:sph:minimal}
+
+This is the simplest fully-conservative version of SPH using the
+internal energy $u$ as a thermal variable that can be written
+down. Its implementation in \swift should be understood as a text-book
+example and template for more advanced implementation. A full
+derivation of the equations can be found in the review of
+\cite{Price2012}. Our implementation follows their equations (27), (43),
+(44), (45), (101), (103) and (104) with $\beta=3$ and $\alpha_u=0$. We
+summarize the equations here. 
+
+\subsubsection{Density and other fluid properties (\nth{1} neighbour loop)}
+
+For a set of particles $i$ with positions $\vec{x}_i$ with velocities
+$\vec{v}_i$, masses $m_i$, sthermal energy per unit mass $u_i$ and
+smoothing length $h_i$, we compute the density for each particle:
+
+\begin{equation}
+  \rho_i \equiv \rho(\vec{x}_i) = \sum_j m_j W(\vec{x}_{ij}, h_i),
+  \label{eq:sph:minimal:rho}
+\end{equation}
+and the derivative of its density with respect to $h$:
+
+\begin{equation}
+    \label{eq:sph:minimal:rho_dh}
+  \rho_{\partial h_i} \equiv \dd{\rho}{h}(\vec{x}_i) = \sum_j m_j \dd{W}{h}(\vec{x}_{ij}
+  , h_i).
+\end{equation}
+The gradient terms (``h-terms'') can then be computed from the density
+and its derivative:
+
+\begin{equation}
+  f_i \equiv \left(1 + \frac{h_i}{3\rho_i}\rho_{\partial h_i}
+  \right)^{-1}.
+  \label{eq:sph:minimal:f_i}
+\end{equation}
+Using the pre-defined equation of state, the pressure $P_i$ and the sound
+speed $c_i$ at the location of particle $i$ can now be computed from
+$\rho_i$ and $u_i$:
+
+\begin{align}
+  P_i &= P_{\rm eos}(\rho_i, u_i),   \label{eq:sph:minimal:P}\\
+  c_i &= c_{\rm eos}(\rho_i, u_i).   \label{eq:sph:minimal:c}
+\end{align}
+
+\subsubsection{Hydrodynamical accelerations (\nth{2} neighbour loop)}
+
+We can then proceed with the second loop over
+neighbours. The signal velocity $v_{{\rm sig},ij}$ between two particles is given by
+
+\begin{align}
+  \mu_{ij} &=
+  \begin{cases}
+  \frac{\vec{v}_{ij} \cdot \vec{x}_{ij}}{|\vec{x}_{ij}|}  & \rm{if}~
+  \vec{v}_{ij} \cdot \vec{x}_{ij} < 0,\\
+    0 &\rm{otherwise}, \\
+  \end{cases}\nonumber\\
+  v_{{\rm sig},ij} &= c_i + c_j - 3\mu_{ij}.   \label{eq:sph:minimal:v_sig}
+\end{align}
+We also use these two quantities for the calculation of the viscosity term:
+
+\begin{equation}
+\nu_{ij} = -\frac{1}{2}\frac{\alpha \mu_{ij} v_{{\rm
+      sig},ij}}{\bar\rho_{ij}}
+  \label{eq:sph:minimal:nu_ij}
+\end{equation}
+The fluid accelerations are then given by
+
+\begin{align}
+  \frac{d\vec{v}_i}{dt} = -\sum_j m_j &\left[\frac{f_iP_i}{\rho_i^2}
+  \nabla_x W(\vec{x}_{ij}, h_i)   \nonumber
+  +\frac{f_jP_j}{\rho_j^2}\nabla_x W(\vec{x}_{ij},h_j)\right. \\
+  &+ \left. \bigg.\nu_{ij} \Wij \right], \label{eq:sph:minimal:dv_dt}
+\end{align}
+and the change in internal energy,
+
+\begin{align}
+  \frac{du_i}{dt} = \sum_j m_j &\left[\frac{f_iP_i}{\rho_i^2}  \vec{v}_{ij}
+    \cdot \nabla_x W(\vec{x}_{ij}, h_i) \right. \label{eq:sph:minimal:du_dt}\\
+    &+\left. \frac{1}{2}\nu_{ij}\vec{v}_{ij}\cdot\Big. \Wij\right], \nonumber
+\end{align}
+where in both cases the first line corresponds to the standard SPH
+term and the second line to the viscuous accelerations.
+
+We also compute an estimator of the change in smoothing length to be
+used in the prediction step. This is an estimate of the local
+divergence of the velocity field compatible with the accelerations
+computed above:
+
+\begin{equation}
+  \frac{dh_i}{dt} = -\frac{1}{3}h_i \sum_j \frac{m_j}{\rho_j}
+  \vec{v}_{ij}\cdot \nabla_x W(\vec{x}_{ij}, h_i).
+  \label{eq:sph:minimal:dh_dt}
+\end{equation}
+and update the signal velocity of the particles:
+
+\begin{equation}
+  v_{{\rm sig},i} = \max_j \left( v_{{\rm sig},ij} \right).
+  \label{eq:sph:minimal:v_sig_update}
+\end{equation}
+All the quantities required for time integration have now been obtained.
+
+\subsubsection{Time integration}
+
+For each particle, we compute a time-step given by the CFL condition:
+
+\begin{equation}
+  \Delta t = 2 C_{\rm CFL} \frac{H_i}{v_{{\rm sig},i}},
+    \label{eq:sph:minimal:dt}
+\end{equation}
+where $C_{\rm CFL}$ is a free dimensionless parameter and $H_i = \gamma h_i$ is the
+kernel support size. Particles can then be ``kicked'' forward in time:
+\begin{align}
+  \vec{v}_i &\rightarrow \vec{v}_i + \frac{d\vec{v}_i}{dt} \Delta t  \label{eq:sph:minimal:kick_v}\\
+  u_i &\rightarrow u_i + \frac{du_i}{dt} \Delta t \label{eq:sph:minimal:kick_u}\\
+  P_i &\rightarrow P_{\rm eos}\left(\rho_i, u_i\right) \label{eq:sph:minimal:kick_P}, \\
+  c_i &\rightarrow c_{\rm eos}\left(\rho_i, u_i\right) \label{eq:sph:minimal:kick_c},
+\end{align}
+where we used the pre-defined equation of state to compute the new
+value of the pressure and sound-speed.
+
+\subsubsection{Particle properties prediction}
+
+Inactive particles need to have their quantities predicted forward in
+time in the ``drift'' operation. We update them as follows:
+
+\begin{align}
+  h_i &\rightarrow h_i \exp\left(\frac{1}{h_i} \frac{dh_i}{dt}
+  \Delta t\right), \label{eq:sph:minimal:drift_h}\\
+  \rho_i &\rightarrow \rho_i \exp\left(-\frac{3}{h_i} \frac{dh_i}{dt}
+  \Delta t\right), \label{eq:sph:minimal:drift_rho} \\
+  P_i &\rightarrow P_{\rm eos}\left(\rho_i, u_i + \frac{du_i}{dt} \Delta t\right), \label{eq:sph:minimal:drift_P}\\
+  c_i &\rightarrow c_{\rm eos}\left(\rho_i, u_i + \frac{du_i}{dt}
+  \Delta t\right) \label{eq:sph:minimal:drift_c},
+\end{align}
+where, as above, the last two updated quantities are obtained using
+the pre-defined equation of state. Note that the thermal energy $u_i$
+itself is \emph{not} updated.
+
+\subsection{Gadget-2 SPH}
+\label{sec:sph:gadget2}
+
+This flavour of SPH is the one implemented in the \gadget-2 code
+\citep{Springel2005}. The basic equations were derived by
+\cite{Springel2002} and also includes a \cite{Balsara1995} switch for
+the suppression of viscosity. The implementation here follows closely the
+presentation of \cite{Springel2005}. Specifically, we use their equations (5), (7),
+(8), (9), (10), (13), (14) and (17). We summarize them here for completeness.
+
+\subsubsection{Density and other fluid properties (\nth{1} neighbour loop)}
+
+For a set of particles $i$ with positions $\vec{x}_i$ with velocities
+$\vec{v}_i$, masses $m_i$, entropic function per unit mass $A_i$ and
+smoothing length $h_i$, we compute the density, derivative of the density with respect
+to $h$ and the ``f-terms'' in a similar way to the minimal-SPH case
+(Eq. \ref{eq:sph:minimal:rho}, \ref{eq:sph:minimal:rho_dh} and
+\ref{eq:sph:minimal:f_i}). From these the pressure and sound-speed can
+be computed using the predefined equation of state:
+
+\begin{align}
+  P_i &= P_{\rm eos}(\rho_i, A_i),   \label{eq:sph:gadget2:P}\\
+  c_i &= c_{\rm eos}(\rho_i, A_i).   \label{eq:sph:gadget2:c}
+\end{align}
+We additionally compute the divergence and
+curl of the velocity field using standard SPH expressions:
+
+\begin{align}
+  \nabla\cdot\vec{v}_i \equiv\nabla\cdot \vec{v}(\vec{x}_i) &= \frac{1}{\rho_i} \sum_j m_j
+  \vec{v}_{ij}\cdot\nabla_x W(\vec{x}_{ij}, h_i) \label{eq:sph:gadget2:div_v},\\ 
+    \nabla\times\vec{v}_i \equiv \nabla\times \vec{v}(\vec{x}_i) &= \frac{1}{\rho_i} \sum_j m_j
+  \vec{v}_{ij}\times\nabla_x W(\vec{x}_{ij}, h_i) \label{eq:sph:gadget2:rot_v}.
+\end{align}
+These are used to construct the \cite{Balsara1995} switch $B_i$:
+
+\begin{equation}
+  B_i = \frac{|\nabla\cdot\vec{v}_i|}{|\nabla\cdot\vec{v}_i| +
+    |\nabla\times\vec{v}_i| + 10^{-4}c_i / h_i}, \label{eq:sph:gadget2:balsara}
+\end{equation}
+where the last term in the denominator is added to prevent numerical instabilities.
+
+\subsubsection{Hydrodynamical accelerations (\nth{2} neighbour loop)}
+
+The accelerations are computed in a similar fashion to the minimal-SPH
+case with the exception of the viscosity term which is now modified to
+include the switch. Instead of Eq. \ref{eq:sph:minimal:nu_ij}, we get:
+
+\begin{equation}
+\nu_{ij} = -\frac{1}{2}\frac{\alpha \bar B_{ij} \mu_{ij} v_{{\rm sig},ij}}{\bar\rho_{ij}},
+  \label{eq:sph:gadget2:nu_ij}  
+\end{equation}
+whilst equations \ref{eq:sph:minimal:v_sig},
+\ref{eq:sph:minimal:dv_dt}, \ref{eq:sph:minimal:dh_dt} and
+\ref{eq:sph:minimal:v_sig_update} remain unchanged. The only other
+change is the equation of motion for the thermodynamic variable which
+now has to be describing the evolution of the entropic function and
+not the evolution of the thermal energy. It reads:
+
+\begin{equation}
+\frac{dA_i}{dt} = \frac{1}{2} \frac{\gamma-1}{\rho^{\gamma-1}} \sum_j
+m_j \nu_{ij}\vec{v}_{ij}\cdot \Wij
+\end{equation}
+
+\subsubsection{Time integration}
+
+The time-step condition is identical to the Minimal-SPH case
+(Eq. \ref{eq:sph:minimal:dt}). The same applies to the integration
+forward in time (Eq. \ref{eq:sph:minimal:kick_v} to
+\ref{eq:sph:minimal:kick_c}) with the exception of the change in
+internal energy (Eq. \ref{eq:sph:minimal:kick_u}) which gets replaced
+by an integration for the the entropy:
+
+
+\begin{align}
+  \vec{v}_i &\rightarrow \vec{v}_i + \frac{d\vec{v}_i}{dt} \Delta t  \label{eq:sph:gadget2:kick_v}\\
+  A_i &\rightarrow A_i + \frac{dA_i}{dt} \Delta t \label{eq:sph:gadget2:kick_A}\\
+  P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i\right) \label{eq:sph:gadget2:kick_P}, \\
+  c_i &\rightarrow c_{\rm eos}\left(\rho_i, A_i\right) \label{eq:sph:gadget2:kick_c},
+\end{align}
+where, once again, we made use of the equation of state relating
+thermodynaimcal quantities.
+
+\subsubsection{Particle properties prediction}
+
+The prediction step is also identical to the Minimal-SPH case with the
+entropic function replacing the thermal energy.
+
+\begin{align}
+  h_i &\rightarrow h_i \exp\left(\frac{1}{h_i} \frac{dh_i}{dt}
+  \Delta t\right), \label{eq:sph:gadget2:drift_h}\\
+  \rho_i &\rightarrow \rho_i \exp\left(-\frac{3}{h_i} \frac{dh_i}{dt}
+  \Delta t\right), \label{eq:sph:gadget2:drift_rho} \\
+  P_i &\rightarrow P_{\rm eos}\left(\rho_i, A_i + \frac{dA_i}{dt} \Delta t\right), \label{eq:sph:gadget2:drift_P}\\
+  c_i &\rightarrow c_{\rm eos}\left(\rho_i, A_i + \frac{dA_i}{dt}
+  \Delta t\right) \label{eq:sph:gadget2:drift_c},
+\end{align}
+where, as above, the last two updated quantities are obtained using
+the pre-defined equation of state. Note that the entropic function $A_i$
+itself is \emph{not} updated.
+
+
+
+\subsection{Pressure-Entropy SPH}
+
+\cite{Hopkins2013}
+\label{sec:sph:pe}
+
+\subsubsection{Hydrodynamical accelerations}
+
+\subsubsection{Particle properties prediction}
+
+\subsection{Pressure-Energy SPH}
+\label{sec:sph:pu}
+
+\cite{Hopkins2013}\\
+to be done
+
+\subsection{Anarchy SPH}
+\label{sec:sph:anarchy}
+
+to be done
+
+\bibliographystyle{mnras}
+\bibliography{./bibliography.bib}
+
+\end{document}

theory/SPH/kernel/bibliography.bib

+@ARTICLE{Dehnen2012,
+   author = {{Dehnen}, W. and {Aly}, H.},
+    title = "{Improving convergence in smoothed particle hydrodynamics simulations without pairing instability}",
+  journal = {\mnras},
+archivePrefix = "arXiv",
+   eprint = {1204.2471},
+ primaryClass = "astro-ph.IM",
+ keywords = {hydrodynamics, methods: numerical },
+     year = 2012,
+    month = sep,
+   volume = 425,
+    pages = {1068-1082},
+      doi = {10.1111/j.1365-2966.2012.21439.x},
+   adsurl = {http://adsabs.harvard.edu/abs/2012MNRAS.425.1068D},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+

theory/SPH/kernel/kernel_definitions.tex

-\documentclass[a4paper]{mnras}
-\usepackage[utf8]{inputenc}
-\usepackage{amsmath}
+\documentclass[fleqn, usenatbib, useAMS,a4paper]{mnras}
 \usepackage{graphicx}
-\usepackage{xspace}
+\usepackage{amsmath,paralist,xcolor,xspace,amssymb}
+\usepackage{times}
 
 \newcommand{\swift}{{\sc Swift}\xspace}
 
-
-
 %opening
 \title{SPH kernels in SWIFT}
 \author{Matthieu Schaller}
@@ -16,7 +13,9 @@
 
 \maketitle
 
-In here we follow the definitions of Dehnen \& Aly 2012.
+We follow the definitions and notation of
+\cite{Dehnen2012}. Motivation for all the material given here can be
+found in their paper.
 
 \section{General Definitions}
 
@@ -28,14 +27,12 @@ The desirable properties of an SPH kernels $W(\vec{x},h)$ are:
 \item $W(\vec{x},h)$ should have a finite support and be cheap to
   compute.
 \end{enumerate}
-
 As a consequence, the smoothing kernels used in SPH can
 hence be written (in 3D) as
 
 \begin{equation}
  W(\vec{x},h) \equiv \frac{1}{H^3}f\left(\frac{|\vec{x}|}{H}\right),
 \end{equation}
-
 where $H=\gamma h$ is defined below and $f(u)$ is a dimensionless
 function, usually a low-order polynomial, such that $f(u \geq 1) = 0$
 and normalised such that
@@ -43,7 +40,6 @@ and normalised such that
 \begin{equation}
   \int f(|\vec{u}|){\rm d}^3u = 1.
 \end{equation}
-
 $H$ is the kernel's support radius and is used as the ``smoothing
 length'' in the Gadget code( {i.e.} $H=h$). This definition is,
 however, not very physical and makes comparison of kernels at a
@@ -56,13 +52,11 @@ standard deviation
   \sigma^2 \equiv \frac{1}{3}\int \vec{u}^2 W(\vec{u},h) {\rm d}^3u.
   \label{eq:sph:sigma}
 \end{equation}
-
 The smoothing length is then:
 \begin{equation}
   h\equiv2\sigma.
     \label{eq:sph:h}
 \end{equation}
-
 Each kernel, {\it i.e.} defintion of $f(u)$, will have a different
 ratio $\gamma = H/h$. So for a \emph{fixed resolution} $h$, one will
 have different kernel support sizes, $H$, and a different number of
@@ -73,18 +67,15 @@ separation:
 \begin{equation}
   h = \eta \langle x \rangle = \eta \left(\frac{m}{\rho}\right)^{1/3},
 \end{equation}
-
 where $\rho$ is the local density of the fluid and $m$ the SPH
 particle mass. 
-
 The (weighted) number of neighbours within the kernel support is a
-useful quantity to use in implementations of SPH. It is defined as (in
-3D)
+useful quantity to use in implementations of SPH. It is defined (in
+3D) as:
 
 \begin{equation}
   N_{\rm ngb} \equiv \frac{4}{3}\pi \left(\frac{H}{h}\eta\right)^3.
 \end{equation}
-
 Once the fixed ratio $\gamma= H/h$ is known (via equations
 \ref{eq:sph:sigma} and \ref{eq:sph:h}) for a given kernel, the number
 of neighbours only depends on the resolution parameter $\eta$.  For
@@ -97,23 +88,22 @@ resolution to $\eta=1.2348$ yields the commonly used value $N_{\rm
 The \swift kernels are split into two categories, the B-splines
 ($M_{4,5,6}$) and the Wendland kernels ($C2$, $C4$ and $C6$). In all
 cases we impose $f(u>1) = 0$.\\
-
 The spline kernels are defined as:
 
 \begin{align}
-  f(u) &= C M_n(u), \\
+  f(u) &= C M_n(u), \nonumber\\
   M_n(u) &\equiv \frac{1}{2\pi}
   \int_{-\infty}^{\infty}
   \left(\frac{\sin\left(k/n\right)}{k/n}\right)^n\cos\left(ku\right){\rm
-  d}k,
+  d}k \nonumber,
 \end{align}
-
-whilst the Wendland kernels read
+whilst the Wendland kernels are constructed from:
 
 \begin{align}
-  f(u) &= C \Psi_{i,j}(u), \\
-  \Psi_{i,j}(u) &\equiv \mathcal{I}^k\left[\max\left(1-u,0\right)\right],\\
-  \mathcal{I}[f](u) &\equiv \int_u^\infty f\left(k\right)k{\rm d}k.
+  f(u) &= C \Psi_{i,j}(u), \nonumber\\
+  \Psi_{i,j}(u) &\equiv
+  \mathcal{I}^k\left[\max\left(1-u,0\right)\right], \nonumber\\
+  \mathcal{I}[f](u) &\equiv \int_u^\infty f\left(k\right)k{\rm d}k. \nonumber
 \end{align}
 
 \subsubsection{Cubic spline ($M_4$) kernel}
@@ -231,11 +221,14 @@ All kernels available in \swift are shown on Fig.~\ref{fig:sph:kernels}.
 The derivatives of the kernel function have relatively simple
 expressions and are shown on Fig.~\ref{fig:sph:kernel_derivatives}:
 
-\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) + 
+\begin{align}
+ \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*}
+\end{align}
 Note that for all the kernels listed above, $f'(0) = 0$. 
 
+\bibliographystyle{mnras}
+\bibliography{./bibliography.bib}
+
 \end{document}

theory/SPH/kernel/run.sh

 #!/bin/bash
 python kernels.py
 pdflatex kernel_definitions.tex
+bibtex kernel_definitions.aux
+pdflatex kernel_definitions.tex
 pdflatex kernel_definitions.tex