Revert "Merge branch 'cosmological_integration' into 'master'"

This reverts commit 28265917, reversing
changes made to 856bcb95.

.gitignore

@@ -113,7 +113,6 @@ tests/testSelectOutput
 tests/testReading
 tests/testSingle
 tests/testTimeIntegration
-tests/testIntegration
 tests/testSPHStep
 tests/testKernel
 tests/testKernelGrav

src/cosmology.h

@@ -163,15 +163,6 @@ struct cosmology {
   /*! Log of final expansion factor */
   double log_a_end;
 
-  /*! Log of starting expansion factor of the interpolation tables */
-  double log_a_table_begin;
-
-  /*! Log of final expansion factor of the interpolation tables */
-  double log_a_table_end;
-
-  /*! Scale-factor interpolation table */
-  double *log_a_interp_table;
-
   /*! Drift factor interpolation table */
   double *drift_fac_interp_table;
 
@@ -184,9 +175,15 @@ struct cosmology {
   /*! Kick factor (hydro correction) interpolation table (GIZMO-MFV only) */
   double *hydro_kick_corr_interp_table;
 
-  /*! Time since Big Bang interpolation table */
+  /*! Time interpolation table */
   double *time_interp_table;
 
+  /*! Scale factor interpolation table */
+  double *scale_factor_interp_table;
+
+  /*! Time between Big Bang and first entry in the table */
+  double time_interp_table_offset;
+
   /*! Time at the present-day (a=1) */
   double universe_age_at_present_day;
 };
@@ -194,8 +191,6 @@ struct cosmology {
 void cosmology_update(struct cosmology *c, const struct phys_const *phys_const,
                       integertime_t ti_current);
 
-double cosmology_get_scale_factor(const struct cosmology *cosmo,
-                                  integertime_t ti);
 double cosmology_get_drift_factor(const struct cosmology *cosmo,
                                   const integertime_t ti_start,
                                   const integertime_t ti_end);
@@ -219,8 +214,7 @@ double cosmology_get_delta_time_from_scale_factors(const struct cosmology *c,
                                                    const double a_start,
                                                    const double a_end);
 
-double cosmology_get_scale_factor_from_time(const struct cosmology *cosmo,
-                                            double t);
+double cosmology_get_scale_factor(const struct cosmology *cosmo, double t);
 
 double cosmology_get_time_since_big_bang(const struct cosmology *c, double a);
 void cosmology_init(struct swift_params *params, const struct unit_system *us,

src/outputlist.c

@@ -106,7 +106,7 @@ void output_list_read_file(struct output_list *outputlist, const char *filename,
     if (type == OUTPUT_LIST_REDSHIFT) *time = 1. / (1. + *time);
 
     if (cosmo && type == OUTPUT_LIST_AGE)
-      *time = cosmology_get_scale_factor_from_time(cosmo, *time);
+      *time = cosmology_get_scale_factor(cosmo, *time);
 
     /* Update size */
     ind += 1;

tests/Makefile.am

@@ -28,7 +28,7 @@ TESTS = testGreetings testMaths testReading.sh testKernel \
         testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
 	testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf \
 	testPotentialPair testEOS testUtilities testSelectOutput.sh \
-	testCbrt testCosmology testIntegration testOutputList testFormat.sh \
+	testCbrt testCosmology testOutputList testFormat.sh \
 	test27cellsStars.sh test27cellsStarsPerturbed.sh testHydroMPIrules
 
 # List of test programs to compile
@@ -40,7 +40,7 @@ check_PROGRAMS = testGreetings testReading testTimeIntegration \
                  testRiemannHLLC testMatrixInversion testDump testLogger \
 		 testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
 		 testGravityDerivatives testPotentialSelf testPotentialPair testEOS testUtilities \
-		 testSelectOutput testCbrt testIntegration testCosmology testOutputList test27cellsStars \
+		 testSelectOutput testCbrt testCosmology testOutputList test27cellsStars \
 		 test27cellsStars_subset testCooling testComovingCooling testFeedback testHashmap testHydroMPIrules
 
 # Rebuild tests when SWIFT is updated.
@@ -59,8 +59,6 @@ testReading_SOURCES = testReading.c
 
 testSelectOutput_SOURCES = testSelectOutput.c
 
-testIntegration_SOURCES = testIntegration.c
-
 testCosmology_SOURCES = testCosmology.c
 
 testOutputList_SOURCES = testOutputList.c

tests/testCosmology.c

@@ -20,15 +20,11 @@
 /* Some standard headers. */
 #include "../config.h"
 
-/* Some standard headers */
-#include <fenv.h>
-#include <math.h>
-
 /* Includes. */
 #include "swift.h"
 
-#define N_CHECK 200
-#define TOLERANCE 1e-6
+#define N_CHECK 20
+#define TOLERANCE 1e-7
 
 void test_params_init(struct swift_params *params) {
   parser_init("", params);
@@ -36,21 +32,12 @@ void test_params_init(struct swift_params *params) {
   parser_set_param(params, "Cosmology:Omega_lambda:0.6910");
   parser_set_param(params, "Cosmology:Omega_b:0.0486");
   parser_set_param(params, "Cosmology:h:0.6774");
-  parser_set_param(params, "Cosmology:a_begin:0.01");
+  parser_set_param(params, "Cosmology:a_begin:0.1");
   parser_set_param(params, "Cosmology:a_end:1.0");
 }
 
 int main(int argc, char *argv[]) {
 
-  /* Initialize CPU frequency, this also starts time. */
-  unsigned long long cpufreq = 0;
-  clocks_set_cpufreq(cpufreq);
-
-/* Choke on FP-exceptions */
-#ifdef HAVE_FE_ENABLE_EXCEPT
-  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
-#endif
-
   message("Initialization...");
 
   /* pseudo initialization of params */
@@ -72,16 +59,15 @@ int main(int argc, char *argv[]) {
   message("Start checking computation...");
 
   for (int i = 0; i < N_CHECK; i++) {
-    double a = 0.01 + 0.99 * i / (N_CHECK - 1.);
-
+    double a = 0.1 + 0.9 * i / (N_CHECK - 1.);
     /* Compute a(t(a)) and check if same results */
-    double time = cosmology_get_time_since_big_bang(&cosmo, a);
-    double my_a = cosmology_get_scale_factor_from_time(&cosmo, time);
+    double tmp = cosmology_get_time_since_big_bang(&cosmo, a);
+    tmp = cosmology_get_scale_factor(&cosmo, tmp);
 
     /* check accuracy */
-    double rel_err = (my_a - a) / a;
-    message("Accuracy of %g at a=%g", rel_err, a);
-    assert(fabs(rel_err) < TOLERANCE);
+    tmp = (tmp - a) / a;
+    message("Accuracy of %g at a=%g", tmp, a);
+    assert(fabs(tmp) < TOLERANCE);
   }
 
   message("Everything seems fine with cosmology.");

tests/testIntegration.c

-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2019 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#include "../config.h"
-
-/* Some standard headers */
-#include <fenv.h>
-#include <math.h>
-
-/* Local headers. */
-#include "swift.h"
-
-#ifdef HAVE_LIBGSL
-#include <gsl/gsl_integration.h>
-#endif
-
-static INLINE double w_tilde(double a, double w0, double wa) {
-  return (a - 1.) * wa - (1. + w0 + wa) * log(a);
-}
-
-static INLINE double E(double Or, double Om, double Ok, double Ol, double w0,
-                       double wa, double a) {
-  const double a_inv = 1. / a;
-  return sqrt(Or * a_inv * a_inv * a_inv * a_inv + Om * a_inv * a_inv * a_inv +
-              Ok * a_inv * a_inv + Ol * exp(3. * w_tilde(a, w0, wa)));
-}
-
-static INLINE double drift_integrand(double a, void *param) {
-
-  const struct cosmology *c = (const struct cosmology *)param;
-  const double Omega_r = c->Omega_r;
-  const double Omega_m = c->Omega_m;
-  const double Omega_k = c->Omega_k;
-  const double Omega_l = c->Omega_lambda;
-  const double w_0 = c->w_0;
-  const double w_a = c->w_a;
-  const double H0 = c->H0;
-
-  const double a_inv = 1. / a;
-  const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w_0, w_a, a);
-  const double H = H0 * E_z;
-  return (1. / H) * a_inv * a_inv * a_inv;
-}
-
-int main(int argc, char *argv[]) {
-
-  /* Initialize CPU frequency, this also starts time. */
-  unsigned long long cpufreq = 0;
-  clocks_set_cpufreq(cpufreq);
-
-/* Choke on FP-exceptions */
-#ifdef HAVE_FE_ENABLE_EXCEPT
-  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
-#endif
-
-  /* Fake parameter file with Planck+13 cosmology
-   * and the usual cosmological system of units */
-  struct swift_params *params =
-      (struct swift_params *)malloc(sizeof(struct swift_params));
-  parser_set_param(params, "InternalUnitSystem:UnitMass_in_cgs:1.98848e43");
-  parser_set_param(params,
-                   "InternalUnitSystem:UnitLength_in_cgs:3.08567758e24");
-  parser_set_param(params, "InternalUnitSystem:UnitVelocity_in_cgs:1e5");
-  parser_set_param(params, "InternalUnitSystem:UnitCurrent_in_cgs:1.");
-  parser_set_param(params, "InternalUnitSystem:UnitTemp_in_cgs:1.");
-  parser_set_param(params, "Cosmology:Omega_m:0.307");
-  parser_set_param(params, "Cosmology:Omega_lambda:0.693");
-  parser_set_param(params, "Cosmology:Omega_b:0.0482519");
-  parser_set_param(params, "Cosmology:h:0.6777");
-  parser_set_param(params, "Cosmology:a_begin:0.0078125");
-  parser_set_param(params, "Cosmology:a_end:1.");
-
-  /* Initialise everything */
-  struct unit_system units;
-  units_init_from_params(&units, params, "InternalUnitSystem");
-
-  struct phys_const phys_const;
-  phys_const_init(&units, params, &phys_const);
-
-  struct cosmology cosmo;
-  cosmology_init(params, &units, &phys_const, &cosmo);
-
-  /* Initalise the GSL workspace */
-  size_t workspace_size = 100000;
-  gsl_integration_workspace *workspace =
-      gsl_integration_workspace_alloc(workspace_size);
-  gsl_function F = {&drift_integrand, &cosmo};
-
-  /* Loop over all the reasonable time-bins */
-  for (int bin = 6; bin < 24; ++bin) {
-
-    const int num_steps = (1LL << bin);
-    const integertime_t time_step_size = max_nr_timesteps / num_steps;
-
-    double min_err = 0.;
-    double max_err = 0.;
-    double sum_err = 0.;
-    integertime_t ti_min = 0;
-    integertime_t ti_max = 0;
-
-    message("Testing %d steps", num_steps);
-
-    /* Cycle through the time-steps */
-    for (integertime_t ti = 0; ti < max_nr_timesteps; ti += time_step_size) {
-
-      const integertime_t ti_beg = ti;
-      const integertime_t ti_end =
-          min(ti + time_step_size, max_nr_timesteps - 1);
-
-      const double a_beg = cosmology_get_scale_factor(&cosmo, ti_beg);
-      const double a_end = cosmology_get_scale_factor(&cosmo, ti_end);
-
-      /* Get the drift factor from SWIFT */
-      const double swift_drift_fac =
-          cosmology_get_drift_factor(&cosmo, ti, ti + time_step_size);
-
-      /* Get the exact drift factor */
-      double exact_drift_fac = 0., abserr;
-      gsl_integration_qag(&F, a_beg, a_end, 0, 1.0e-12, workspace_size,
-                          GSL_INTEG_GAUSS61, workspace, &exact_drift_fac,
-                          &abserr);
-
-      const double rel_err = 0.5 * (swift_drift_fac - exact_drift_fac) /
-                             (swift_drift_fac + exact_drift_fac);
-
-      if (rel_err > max_err) {
-        max_err = rel_err;
-        ti_max = ti;
-      }
-
-      if (rel_err < min_err) {
-        min_err = rel_err;
-        ti_min = ti;
-      }
-
-      sum_err += fabs(rel_err);
-    }
-
-    message("Max  error: %14e at a=[%.9f %.9f] ", max_err,
-            cosmology_get_scale_factor(&cosmo, ti_max),
-            cosmology_get_scale_factor(&cosmo, ti_max + time_step_size));
-    message("Min  error: %14e at a=[%.9f %.9f]", min_err,
-            cosmology_get_scale_factor(&cosmo, ti_min),
-            cosmology_get_scale_factor(&cosmo, ti_min + time_step_size));
-    message("Sum  error: %14e", sum_err);
-    message("Mean error: %14e", sum_err / num_steps);
-
-    if (max_err > 1e-4 || min_err < -1e-4)
-      error("Error too large to be acceptable");
-  }
-
-#ifdef HAVE_LIBGSL
-
-  return 0;
-#else
-  return 0;
-#endif
-}

tests/testOutputList.c

@@ -19,10 +19,6 @@
  ******************************************************************************/
 #include "../config.h"
 
-/* Some standard headers */
-#include <fenv.h>
-#include <math.h>
-
 /* Includes. */
 #include "swift.h"
 
@@ -90,8 +86,7 @@ void test_cosmo(struct engine *e, const char *name, const int with_assert) {
   double output_time = 0;
 
   /* Test Time */
-  e->time_base =
-      log(e->cosmology->a_end / e->cosmology->a_begin) / max_nr_timesteps;
+  e->time_base = log(e->time_end / e->cosmology->a_begin) / max_nr_timesteps;
   e->ti_current = 0;
   e->policy = engine_policy_cosmology;
 
@@ -120,16 +115,6 @@ void test_cosmo(struct engine *e, const char *name, const int with_assert) {
 };
 
 int main(int argc, char *argv[]) {
-
-  /* Initialize CPU frequency, this also starts time. */
-  unsigned long long cpufreq = 0;
-  clocks_set_cpufreq(cpufreq);
-
-/* Choke on FP-exceptions */
-#ifdef HAVE_FE_ENABLE_EXCEPT
-  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
-#endif
-
   /* Create a structure to read file into. */
   struct swift_params params;
 

theory/Cosmology/coordinates.tex

 \subsection{Choice of co-moving coordinates}
 \label{ssec:ccordinates}
 
-Note that, unlike the convention used by many codes, we do not express
-quantities with ``little h'' ($h$) included; for instance units of
-length are expressed in units of $\rm{Mpc}$ and not
-${\rm{Mpc}}/h$. Similarly, the time integration operators (see
-below)include an $h$-factor via the explicit appearance of
+Note that, unlike the \gadget convention, we do not express quantities with
+``little h'' ($h$) included; for instance units of length are expressed in
+units of $\rm{Mpc}$ and not ${\rm{Mpc}}/h$. Similarly, the time integration
+operators (see below) also include an $h$-factor via the explicit appearance of
 the Hubble constant.\\
-
 In physical coordinates, the Lagrangian for a particle $i$ in the
 \cite{Springel2002} flavour of SPH with gravity reads
 \begin{equation}
@@ -63,8 +61,82 @@ codes\footnote{One additional inconvenience of our choice of
   scale-factor. The signal velocity entering the time-step calculation
   will hence read
   $v_{\rm sig} = a\dot{\mathbf{r}'} + c = \frac{1}{a} \left(
-    |\mathbf{v}'| + a^{(5 - 3\gamma)/2}c'\right)$.}.  For the
-derivatives, this choice implies
-$\dot{\mathbf{v}'} = 2H\mathbf{v}' + a^2\ddot{\mathbf{r}'}$.
+    |\mathbf{v}'| + a^{(5 - 3\gamma)/2}c'\right)$.}.
+This choice implies that $\dot{\mathbf{v}'} = 2H\mathbf{v}' +  a^2\ddot{\mathbf{r}'}$.
+
+\subsubsection{SPH equations}
+
+Using the SPH 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
+\begin{alignat}{3}
+  \dot{\mathbf{r}}_i'&= \frac{1}{a^2} \mathbf{v}_i'&  \label{eq:cosmo_eom_r} \\
+  \dot{\mathbf{v}}_i' &= -\sum_j m_j &&\left[\frac{1}{a^{3(\gamma-1)}}f_i'A_i'\rho_i'^{\gamma-2}\mathbf{\nabla}_i'W_{ij}'(h_i)\right. \nonumber\\
+  &   && + \left. \frac{1}{a^{3(\gamma-1)}}f_j'A_j'\rho_j'^{\gamma-2}\mathbf{\nabla}_i'W_{ij}'(h_j)\right. \nonumber\\
+  &   && + \left. \frac{1}{a}\mathbf{\nabla}_i'\phi'\right] \label{eq:cosmo_eom_v}
+\end{alignat}
+with
+\begin{equation}
+    f_i' = \left[1 + \frac{h_i'}{3\rho_i'}\frac{\partial
+      \rho_i'}{\partial h_i'}\right]^{-1}, \qquad \mathbf{\nabla}_i'
+  \equiv \frac{\partial}{\partial \mathbf{r}_{i}'}. \nonumber
+\end{equation}
+These correspond to the equations of motion for density-entropy SPH
+\citep[e.g. eq. 14 of][]{Hopkins2013} with cosmological and
+gravitational terms. SPH flavours that evolve the internal energy $u$ instead of the
+entropy require the additional equation of motion describing the evolution of
+$u'$:
+\begin{equation}
+  \dot{u}_i' = \frac{1}{a^2}\frac{P_i'}{\rho_i'^2} f_i'\sum_jm_j\left(\mathbf{v}_i' -
+    \mathbf{v}_j'\right)\cdot\mathbf{\nabla}_i'W_{ij}'(h_i).
+  \label{eq:cosmo_eom_u}
+\end{equation}
 
+In all these cases, the scale-factors appearing in the equations are
+later absorbed in the time-integration operators
+(Sec.~\ref{ssec:operators}) such that the RHS of the equations of
+motions is identical for the primed quantities to the ones obtained in
+the non-cosmological case for the physical quantities.
+
+Additional terms in the SPH equations of motion (e.g. viscosity
+switches) often rely on the velocity divergence and curl. Their SPH
+estimators $\langle\cdot\rangle$ in physical coordinates can be
+related to their estimators based on our primed-coordinates using:
+\begin{align}
+  \left\langle \mathbf{\nabla}\cdot\dot{\mathbf{r}}_i \right\rangle &=
+  \frac{1}{a^2} \left\langle
+  \mathbf{\nabla}'\cdot\mathbf{v}_i'\right\rangle =
+  \frac{1}{a^2\rho_i'}\sum_j m_j\left(\mathbf{v}_j' -
+  \mathbf{v}_i'\right) \cdot \mathbf{\nabla}_i'W_{ij}'(h_i), \nonumber \\
+  \left\langle \mathbf{\nabla}\times\dot{\mathbf{r}}_i \right\rangle &=
+  \frac{1}{a^2} \left\langle
+  \mathbf{\nabla}'\times\mathbf{v}_i'\right\rangle =
+  \frac{1}{a^2\rho_i'}\sum_j m_j\left(\mathbf{v}_j' -
+  \mathbf{v}_i'\right) \times \mathbf{\nabla}_i'W_{ij}'(h_i). \nonumber
+\end{align}
+We finally give the expressions for the \cite{Monaghan1997} viscosity
+term, that enters most SPH flavours, in our system of coordinates. The
+viscosity ``tensor'' $\Pi_{ij} \equiv \frac{1}{2}\alpha_{\rm visc} v_{\rm
+    sig}'\mu_{ij}/\left(\rho_i + \rho_j\right)$ can be computed
+  in the primmed coordinates from the following quantities
+\begin{align}
+  \omega_{ij}' &= \left(\mathbf{v}_i' - \mathbf{v}_j'\right) \cdot
+  \left(\mathbf{r}_i' - \mathbf{r}_j'\right), \\
+  \mu_{ij}' &=
+  a^{(3\gamma-5)/2} \left(\omega_{ij}' + a^2H\left|\mathbf{r}_i' -
+  \mathbf{r}_j'\right|^2 \right) / |\mathbf{r}_i' - \mathbf{r}_j' |,
+  \\
+  v_{\rm sig}' &= c_i' + c_j' - 3 \mu_{ij}'.
+\end{align}
+which leads to $\Pi_{ij}'=a^{3\gamma}\Pi_{ij}$. Note that he last quantity is
+also used to compute the time-step size. The evolution of entropy is
+then
+\begin{equation}
+  \dot{A}_i = \dot{A}_i' = \frac{1}{a^2}\frac{1}{2}\frac{\gamma-1}{\rho_i'^{\gamma-1}} \sum_j
+  m_j \Pi_{ij}' \left(\mathbf{v}_i' -
+  \mathbf{v}_j'\right)\cdot\mathbf{\nabla}_i'W_{ij}'(h_i),
+\end{equation}
+indicating that the entropy evolves with the same scale-factor
+dependence as the comoving positions (eq.~\ref{eq:cosmo_eom_r}).
 

theory/Cosmology/cosmology_standalone.tex

@@ -26,7 +26,7 @@
 \label{firstpage}
 
 \begin{abstract}
-Notes on the cosmology equations evolved in \swift
+Making cosmology great again.
 \end{abstract}
 
 \begin{keywords}
@@ -39,15 +39,13 @@ Notes on the cosmology equations evolved in \swift
 
 \input{coordinates}
 
-\input{sph}
-
 \input{operators}
 
 \input{timesteps}
 
-%\input{artificialvisc}
+\input{artificialvisc}
 
-%\input{gizmo}
+\input{gizmo}
 
 \bibliographystyle{mnras}
 \bibliography{./bibliography.bib}

theory/Cosmology/flrw.tex

@@ -48,7 +48,7 @@ of $a_{\rm age}$ equally spaced between $\log(a_{\rm start})$ and $\log(a_{\rm
 end})$. The values are obtained via adaptive quadrature using the 61-points
 Gauss-Konrod rule implemented in the {\sc gsl} library \citep{GSL} with a
 relative error limit of $\epsilon=10^{-10}$. The value for a specific $a$ (over
-the course of a simulation run) is then obtained by cubic interpolation of the
+the course of a simulation run) is then obtained by linear interpolation of the
 table.
 
 \subsubsection{Additional quantities}

theory/Cosmology/operators.tex

@@ -47,12 +47,11 @@ and the corresponding change in redshift:
   \Delta z &= \frac{1}{a_n} - \frac{1}{a_{n+1}} \approx -\frac{H}{a} \Delta t_{\rm cosmic}.
 \end{align}
 Following the same method as for the age of the Universe
-(sec. \ref{ssec:flrw}), these three non-trivial integrals are evaluated
-numerically at the start of the simulation for a series $10^4$ values of
-$a$ placed at regular intervals between $\log a_{\rm begin}$ and
-$\log a_{\rm end}$. The values for a specific pair of scale-factors $a_n$
-and $a_{n+1}$ are then obtained by interpolating that table with a cubic
-polynomial. An accuracy better than one part in $10^8$ is achieved for all
-reasonable cosmologies and time-step lengths.
+(sec. \ref{ssec:flrw}), these three non-trivial integrals are
+evaluated numerically at the start of the simulation for a series
+$10^4$ values of $a$ placed at regular intervals between $\log a_{\rm
+  begin}$ and $\log a_{\rm end}$. The values for a specific pair of
+scale-factors $a_n$ and $a_{n+1}$ are then obtained by interpolating
+that table linearly.
 
 

theory/Cosmology/sph.tex

-\subsubsection{SPH equations}
-
-Using the SPH 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
-\begin{alignat}{3}
-  \dot{\mathbf{r}}_i'&= \frac{1}{a^2} \mathbf{v}_i'&  \label{eq:cosmo_eom_r} \\
-  \dot{\mathbf{v}}_i' &= -\sum_j m_j &&\left[\frac{1}{a^{3(\gamma-1)}}f_i'A_i'\rho_i'^{\gamma-2}\mathbf{\nabla}_i'W_{ij}'(h_i)\right. \nonumber\\
-  &   && + \left. \frac{1}{a^{3(\gamma-1)}}f_j'A_j'\rho_j'^{\gamma-2}\mathbf{\nabla}_i'W_{ij}'(h_j)\right. \nonumber\\
-  &   && + \left. \frac{1}{a}\mathbf{\nabla}_i'\phi'\right] \label{eq:cosmo_eom_v}
-\end{alignat}
-with
-\begin{equation}
-    f_i' = \left[1 + \frac{h_i'}{3\rho_i'}\frac{\partial
-      \rho_i'}{\partial h_i'}\right]^{-1}, \qquad \mathbf{\nabla}_i'
-  \equiv \frac{\partial}{\partial \mathbf{r}_{i}'}. \nonumber
-\end{equation}
-These correspond to the equations of motion for density-entropy SPH
-\citep[e.g. eq. 14 of][]{Hopkins2013} with cosmological and
-gravitational terms. SPH flavours that evolve the internal energy $u$ instead of the
-entropy require the additional equation of motion describing the evolution of
-$u'$:
-\begin{equation}
-  \dot{u}_i' = \frac{1}{a^2}\frac{P_i'}{\rho_i'^2} f_i'\sum_jm_j\left(\mathbf{v}_i' -
-    \mathbf{v}_j'\right)\cdot\mathbf{\nabla}_i'W_{ij}'(h_i).
-  \label{eq:cosmo_eom_u}
-\end{equation}
-In all these cases, the scale-factors appearing as the first term on
-the RHS of the equations are later absorbed in the time-integration
-operators (Sec.~\ref{ssec:operators}) such that the RHS of the
-equations of motions is identical for the primed quantities to the
-ones obtained in the non-cosmological case for the physical
-quantities.
-
-Additional terms in the SPH equations of motion (e.g. viscosity
-switches) often rely on the velocity divergence and curl. Their SPH
-estimators $\langle\cdot\rangle$ in physical coordinates can be
-related to their estimators based on our primed-coordinates using:
-\begin{align}
-  \left\langle \mathbf{\nabla}\cdot\dot{\mathbf{r}}_i \right\rangle &=
-  \frac{1}{a^2} \left\langle
-  \mathbf{\nabla}'\cdot\mathbf{v}_i'\right\rangle =
-  \frac{1}{a^2\rho_i'}\sum_j m_j\left(\mathbf{v}_j' -
-  \mathbf{v}_i'\right) \cdot \mathbf{\nabla}_i'W_{ij}'(h_i), \nonumber \\
-  \left\langle \mathbf{\nabla}\times\dot{\mathbf{r}}_i \right\rangle &=
-  \frac{1}{a^2} \left\langle
-  \mathbf{\nabla}'\times\mathbf{v}_i'\right\rangle =
-  \frac{1}{a^2\rho_i'}\sum_j m_j\left(\mathbf{v}_j' -
-  \mathbf{v}_i'\right) \times \mathbf{\nabla}_i'W_{ij}'(h_i). \nonumber
-\end{align}
-We finally give the expressions for the \cite{Monaghan1997} viscosity
-term, that enters most SPH flavours, in our system of coordinates. The
-viscosity ``tensor'' $\Pi_{ij} \equiv \frac{1}{2}\alpha_{\rm visc} v_{\rm
-    sig}'\mu_{ij}/\left(\rho_i + \rho_j\right)$ can be computed
-  in the primmed coordinates from the following quantities
-\begin{align}
-  \omega_{ij}' &= \left(\mathbf{v}_i' - \mathbf{v}_j'\right) \cdot
-  \left(\mathbf{r}_i' - \mathbf{r}_j'\right), \\
-  \mu_{ij}' &=
-  a^{(3\gamma-5)/2} \left(\omega_{ij}' + a^2H\left|\mathbf{r}_i' -
-  \mathbf{r}_j'\right|^2 \right) / |\mathbf{r}_i' - \mathbf{r}_j' |,
-  \\
-  v_{\rm sig}' &= c_i' + c_j' - 3 \mu_{ij}'.
-\end{align}
-which leads to $\Pi_{ij}'=a^{3\gamma}\Pi_{ij}$. Note that he last quantity is
-also used to compute the time-step size. The evolution of entropy is
-then
-\begin{equation}
-  \dot{A}_i = \dot{A}_i' = \frac{1}{a^2}\frac{1}{2}\frac{\gamma-1}{\rho_i'^{\gamma-1}} \sum_j
-  m_j \Pi_{ij}' \left(\mathbf{v}_i' -
-  \mathbf{v}_j'\right)\cdot\mathbf{\nabla}_i'W_{ij}'(h_i),
-\end{equation}
-indicating that the entropy evolves with the same scale-factor
-dependence as the comoving positions (eq.~\ref{eq:cosmo_eom_r}).