Added definition of the max RMS displacement in cosmology runs to the theory documents.

src/timestep.h

@@ -84,7 +84,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_gpart_timestep(
   /* Take the minimum of all */
   float new_dt = min(new_dt_self, new_dt_ext);
 
-  /* Apply cosmology correction */
+  /* Apply cosmology correction (This is 1 if non-cosmological) */
   new_dt *= e->cosmology->time_step_factor;
 
   /* Limit timestep within the allowed range */
@@ -148,7 +148,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
 
   new_dt = min(new_dt, dt_h_change);
 
-  /* Apply cosmology correction (H==1 if non-cosmological) */
+  /* Apply cosmology correction (This is 1 if non-cosmological) */
   new_dt *= e->cosmology->time_step_factor;
 
   /* Limit timestep within the allowed range */
@@ -191,7 +191,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
   /* Take the minimum of all */
   float new_dt = min3(new_dt_star, new_dt_self, new_dt_ext);
 
-  /* Apply cosmology correction (H==1 if non-cosmological) */
+  /* Apply cosmology correction (This is 1 if non-cosmological) */
   new_dt *= e->cosmology->time_step_factor;
 
   /* Limit timestep within the allowed range */

theory/Cosmology/coordinates.tex

@@ -51,7 +51,7 @@ velocities\footnote{One additional inconvenience of our choice of
   will hence read $v_{\rm sig} = a\dot{\mathbf{r}'} + c = \frac{1}{a} \left(
     |\mathbf{v}'| + a^{(5 - 3\gamma)/2}c'\right)$.}.
 This choice implies that $\dot{v}' = a \ddot{r}$. Using the SPH
-definition of density, $\rho_i =
+definition of density, $\rho_i' =
 \sum_jm_jW(\mathbf{r}_{j}'-\mathbf{r}_{i}',h_i') =
 \sum_jm_jW_{ij}'(h_i')$, we can follow \cite{Price2012} and apply the
 Euler-Lagrange equations to write

theory/Cosmology/cosmology_standalone.tex

@@ -40,6 +40,8 @@ Making cosmology great again.
 
 \input{operators}
 
+\input{timesteps}
+
 \bibliographystyle{mnras}
 \bibliography{./bibliography.bib}
 

theory/Cosmology/flrw.tex

@@ -51,6 +51,20 @@ relative error limit of $\epsilon=10^{-10}$. The value for a specific $a$ (over
 the course of a simulation run) is then obtained by linear interpolation of the
 table.
 
+\subsubsection{Additional quantities}
+
+\swift computes additional quantities that enter common expressions appearing in
+cosmological simulations. For completeness we give their expressions here. The
+look-back time is defined as
+\begin{equation}
+  t_{\rm lookback} = \int_{a_{\rm lookback}}^{a_0} \frac{da}{a E(a)},
+\end{equation}
+and the critical density for closure at a given redshift as
+\begin{equation}
+  \rho_{\rm crit} = \frac{3H(a)^2}{8\pi G_{\rm{N}}}.
+\end{equation}
+These quantities are computed every time-step.
+
 \subsubsection{Typical Values of the Cosmological Parameters}
 
 Typical values for the constants are: $\Omega_m = 0.3, \Omega_\Lambda=0.7, 0 <

theory/Cosmology/timesteps.tex

+\subsection{Choice of time-step size}
+\label{ssec:timesteps}
+
+When running \swift with cosmological time-integration switched on, the
+time-stepping algorithm gets modified in two ways. An additional criterion is
+used to limit the maximal distance a particle can move and the integer time-line
+used for the time-steps changes meaning and represents jumps in scale-factor $a$, 
+hence requiring an additional conversion.
+
+\subsubsection{Maximal displacement}
+
+to prevent particles from moving on trajectories that do not include the effects
+of the expansion of the Universe, we compute a maximal time-step for the
+particles based on their RMS peculiar motion:
+\begin{equation}
+  \Delta t_{\rm cosmo} \equiv \mathcal{C}_{\rm RMS} \frac{a^2 d_{\rm p}}{\sqrt{\frac{1}{N_{\rm p}}\sum_i | \mathbf{v}_i' |^2}},
+  \label{eq:dt_RMS}
+\end{equation}
+where the sum runs over all particles of a species $p$, $\mathcal{C}_{\rm RMS}$
+is a free parameter, $N_{\rm p}$ is the number of baryonic or non-baryonic
+particles, and $d_{\rm p}$ is the mean inter-particle separation for the
+particle with the lowest mass $m_i$ of a given species:
+\begin{equation}
+  d_{\rm baryons} \equiv \left(\frac{m_i}{\Omega_{\rm b} \rho_{\rm crit}}\right)^{1/3}, \quad d_{\rm DM} \equiv \left(\frac{m_i}{\left(\Omega_{\rm m} - \Omega_{\rm b}\right) \rho_{\rm crit}}\right)^{1/3}.
+  \nonumber
+\end{equation}
+We typically use $\mathcal{C}_{\rm RMS} = 0.25$ and given the slow evolution of
+this maximal time-step size, we only re-compute it every time the tree is
+reconstructed.
+
+We also apply an additional criterion based on the smoothing scale of the forces
+computed from the top-level mesh.  In eq.~\ref{eq:dt_RMS}, we replace
+$d_{\rm p}$ by $a_{\rm smooth} \frac{L_{\rm box}}{N_{\rm mesh}}$, where we used
+the definition of the mesh parameters introduced earlier. Given the rather
+coarse mesh usually used in \swift, this time-step condition rarely dominates
+the overall time-step size calculation.
+
+\subsubsection{Conversion from time to integer time-line} 
+