From 17eb2b1628393389c129324cd99dd83e413adb39 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Wed, 6 Jun 2018 18:26:06 +0200
Subject: [PATCH] Correct conversion of peculiar velocities to internal
velocities (and vice-versa) in the i/o modules.
---
src/gravity/Default/gravity_io.h | 6 +++---
src/hydro/Default/hydro_io.h | 6 +++---
src/hydro/Gadget2/hydro_io.h | 6 +++---
src/hydro/GizmoMFM/hydro_io.h | 6 +++---
src/hydro/GizmoMFV/hydro_io.h | 6 +++---
src/hydro/Minimal/hydro_io.h | 6 +++---
src/hydro/MinimalMultiMat/hydro_io.h | 6 +++---
src/hydro/PressureEnergy/hydro_io.h | 6 +++---
src/hydro/PressureEntropy/hydro_io.h | 6 +++---
src/space.c | 6 +++---
theory/Cosmology/coordinates.tex | 28 +++++++++++++++++++---------
11 files changed, 49 insertions(+), 39 deletions(-)
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
index 5016803b1f..7f45317964 100644
--- a/src/gravity/Default/gravity_io.h
+++ b/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;
}
/**
diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h
index d6b6f36ba2..d47c96fbf3 100644
--- a/src/hydro/Default/hydro_io.h
+++ b/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,
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 8f503230e1..3f2af41dc7 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/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,
diff --git a/src/hydro/GizmoMFM/hydro_io.h b/src/hydro/GizmoMFM/hydro_io.h
index 6ce5a23eee..59d579f70c 100644
--- a/src/hydro/GizmoMFM/hydro_io.h
+++ b/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,
diff --git a/src/hydro/GizmoMFV/hydro_io.h b/src/hydro/GizmoMFV/hydro_io.h
index 7b1c651698..92e4378f07 100644
--- a/src/hydro/GizmoMFV/hydro_io.h
+++ b/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,
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index 833103e7ab..879255640f 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/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,
diff --git a/src/hydro/MinimalMultiMat/hydro_io.h b/src/hydro/MinimalMultiMat/hydro_io.h
index 5542b93ae4..61301157e1 100644
--- a/src/hydro/MinimalMultiMat/hydro_io.h
+++ b/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,
diff --git a/src/hydro/PressureEnergy/hydro_io.h b/src/hydro/PressureEnergy/hydro_io.h
index 4f168ed819..78967faec2 100644
--- a/src/hydro/PressureEnergy/hydro_io.h
+++ b/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;
}
/**
diff --git a/src/hydro/PressureEntropy/hydro_io.h b/src/hydro/PressureEntropy/hydro_io.h
index f2247c234a..8c11bf6e33 100644
--- a/src/hydro/PressureEntropy/hydro_io.h
+++ b/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,
diff --git a/src/space.c b/src/space.c
index aaea061e8e..3efb1080c3 100644
--- a/src/space.c
+++ b/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 */
diff --git a/theory/Cosmology/coordinates.tex b/theory/Cosmology/coordinates.tex
index e3a22eaae0..bc59360602 100644
--- a/theory/Cosmology/coordinates.tex
+++ b/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}
--
GitLab