Correct conversion of peculiar velocities to internal velocities (and...

Correct conversion of peculiar velocities to internal velocities (and vice-versa) in the i/o modules.

src/gravity/Default/gravity_io.h

@@ -62,9 +62,9 @@ INLINE static void convert_gpart_vel(const struct engine* e,
   ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav;
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 /**

src/hydro/Default/hydro_io.h

@@ -101,9 +101,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,

src/hydro/Gadget2/hydro_io.h

@@ -113,9 +113,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,

src/hydro/GizmoMFM/hydro_io.h

@@ -158,9 +158,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,

src/hydro/GizmoMFV/hydro_io.h

@@ -158,9 +158,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,

src/hydro/Minimal/hydro_io.h

@@ -127,9 +127,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,

src/hydro/MinimalMultiMat/hydro_io.h

@@ -130,9 +130,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,

src/hydro/PressureEnergy/hydro_io.h

@@ -131,9 +131,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 /**

src/hydro/PressureEntropy/hydro_io.h

@@ -125,9 +125,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,

src/space.c

@@ -2354,7 +2354,7 @@ void space_first_init_parts_mapper(void *restrict map_data, int count,
   const struct cosmology *cosmo = s->e->cosmology;
   const struct phys_const *phys_const = s->e->physical_constants;
   const struct unit_system *us = s->e->internal_units;
-  const float a_factor_vel = cosmo->a * cosmo->a;
+  const float a_factor_vel = cosmo->a;
 
   const struct hydro_props *hydro_props = s->e->hydro_properties;
   const float u_init = hydro_props->initial_internal_energy;
@@ -2429,7 +2429,7 @@ void space_first_init_gparts_mapper(void *restrict map_data, int count,
   const struct space *restrict s = (struct space *)extra_data;
 
   const struct cosmology *cosmo = s->e->cosmology;
-  const float a_factor_vel = cosmo->a * cosmo->a;
+  const float a_factor_vel = cosmo->a;
   const struct gravity_props *grav_props = s->e->gravity_properties;
 
   for (int k = 0; k < count; k++) {
@@ -2486,7 +2486,7 @@ void space_first_init_sparts_mapper(void *restrict map_data, int count,
 #endif
 
   const struct cosmology *cosmo = s->e->cosmology;
-  const float a_factor_vel = cosmo->a * cosmo->a;
+  const float a_factor_vel = cosmo->a;
 
   for (int k = 0; k < count; k++) {
     /* Convert velocities to internal units */

theory/Cosmology/coordinates.tex

@@ -40,19 +40,29 @@ $\Psi \equiv \frac{1}{2}a\dot{a}\mathbf{r}_i^2$ and obtain
   -\frac{\phi'}{a},\\
   \phi' &= a\phi + \frac{1}{2}a^2\ddot{a}\mathbf{r}_i'^2.\nonumber
 \end{align}
-Finally, we introduce the velocities $\mathbf{v}' \equiv
-a^2\dot{\mathbf{r}'}$ that are used internally by the code. Note that these
-velocities \emph{do not} have a physical interpretation. We caution that they
-are not the peculiar velocities, nor the Hubble flow, nor the total
-velocities\footnote{One additional inconvenience of our choice of
+Finally, we introduce the velocities
+$\mathbf{v}' \equiv a^2\dot{\mathbf{r}'}$ that are used internally by
+the code. Note that these velocities \emph{do not} have a physical
+interpretation. We caution that they are not the peculiar velocities
+($\mathbf{v}_{\rm p} \equiv a\dot{\mathbf{r}'} =
+\frac{1}{a}\mathbf{v}'$), nor the Hubble flow
+($\mathbf{v}_{\rm H} \equiv \dot{a}\mathbf{r}'$), nor the total
+velocities
+($\mathbf{v}_{\rm tot} \equiv \mathbf{v}_{\rm p} + \mathbf{v}_{\rm H}
+= \dot{a}\mathbf{r}' + \frac{1}{a}\mathbf{v}'$) and also differ from
+the convention used in \gadget snapshots
+($\sqrt{a} \dot{\mathbf{r}'}$) and other related simulation
+codes\footnote{One additional inconvenience of our choice of
   generalised coordinates is that our velocities $\mathbf{v}'$ and
   sound-speed $c'$ do not have the same dependencies on the
   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(
+  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' =
-\sum_jm_jW(\mathbf{r}_{j}'-\mathbf{r}_{i}',h_i') =
+% This choice implies that $\dot{v}' = a \ddot{r}$.
+
+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}