diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 1228aa69abfe812e5b2e73f066bb13be3292aa20..d541ee777423ef13af00452f550e8cf12b839742 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -160,6 +160,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
return p->mass;
}
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part *restrict p, float m) {
+
+ p->mass = m;
+}
+
/**
* @brief Returns the velocities drifted to the current time of a particle.
*
@@ -500,9 +512,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
* Requires the density to be known
*
* @param p The particle to act upon
+ * @param xp The extended data.
+ * @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part *restrict p, struct xpart *restrict xp) {}
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {}
/**
* @brief Initialises the particles for the first time
@@ -529,4 +544,21 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
hydro_init_part(p, NULL);
}
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part *p, float u_init) {
+
+ p->u = u_init;
+}
+
#endif /* SWIFT_DEFAULT_HYDRO_H */
diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h
index bfff08dcdf061d86986ffdeabee6e7a84cba1c7b..d3b3f08423d2395b006a5b10316abc4a865d72e5 100644
--- a/src/hydro/Default/hydro_io.h
+++ b/src/hydro/Default/hydro_io.h
@@ -54,7 +54,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
}
void convert_part_pos(const struct engine* e, const struct part* p,
- double* ret) {
+ const struct xpart* xp, double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
@@ -67,6 +67,40 @@ void convert_part_pos(const struct engine* e, const struct part* p,
}
}
+void convert_part_vel(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ 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;
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
@@ -74,16 +108,17 @@ void convert_part_pos(const struct engine* e, const struct part* p,
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
-void hydro_write_particles(struct part* parts, struct io_props* list,
- int* num_fields) {
+void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
+ struct io_props* list, int* num_fields) {
*num_fields = 7;
/* List what we want to write */
- list[0] = io_make_output_field_convert_part(
- "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
- list[1] =
- io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
+ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+ UNIT_CONV_LENGTH, parts, xparts,
+ convert_part_pos);
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
list[2] =
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 7df5f35e2317275fb98ec93381ad1b9ccd27ff56..ad60465951f41d96ade35a01d7c1cfe8828f90ec 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -170,6 +170,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
return p->mass;
}
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part *restrict p, float m) {
+
+ p->mass = m;
+}
+
/**
* @brief Returns the velocities drifted to the current time of a particle.
*
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 0faa4bf39ae533ee9942898bdc08118ee4078868..66cfac1be545f2aa26bf08c7b4786c6153c606bf 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -517,7 +517,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part* p, struct xpart* xp) {}
+ struct part* p, struct xpart* xp, const struct cosmology* cosmo) {}
/**
* @brief Extra operations to be done during the drift
@@ -812,6 +812,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
return p->conserved.mass;
}
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part* restrict p, float m) {
+
+ p->conserved.mass = m;
+}
+
/**
* @brief Returns the velocities drifted to the current time of a particle.
*
@@ -917,4 +929,29 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
p->primitives.P = S * pow_gamma(p->primitives.rho);
}
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part* p, float u_init) {
+
+ p->conserved.energy = u_init * p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+ /* add the kinetic energy */
+ p->conserved.energy += 0.5f * p->conserved.mass *
+ (p->conserved.momentum[0] * p->primitives.v[0] +
+ p->conserved.momentum[1] * p->primitives.v[1] +
+ p->conserved.momentum[2] * p->primitives.v[2]);
+#endif
+ p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u_init;
+}
+
#endif /* SWIFT_GIZMO_HYDRO_H */
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index b8e2e30646c3aafb73d387010ab7ea0b0a663a77..4a675b19512fea6e78e69b654c423513bfe3aabc 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -72,7 +72,8 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
* @param p Particle.
* @param ret (return) Internal energy of the particle
*/
-void convert_u(const struct engine* e, const struct part* p, float* ret) {
+void convert_u(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_internal_energy(p);
}
@@ -84,7 +85,8 @@ void convert_u(const struct engine* e, const struct part* p, float* ret) {
* @param p Particle.
* @param ret (return) Entropic function of the particle
*/
-void convert_A(const struct engine* e, const struct part* p, float* ret) {
+void convert_A(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_entropy(p);
}
@@ -95,7 +97,8 @@ void convert_A(const struct engine* e, const struct part* p, float* ret) {
* @param p Particle.
* @return Total energy of the particle
*/
-void convert_Etot(const struct engine* e, const struct part* p, float* ret) {
+void convert_Etot(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
#ifdef GIZMO_TOTAL_ENERGY
ret[0] = p->conserved.energy;
#else
@@ -110,7 +113,7 @@ void convert_Etot(const struct engine* e, const struct part* p, float* ret) {
}
void convert_part_pos(const struct engine* e, const struct part* p,
- double* ret) {
+ const struct xpart* xp, double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
@@ -123,6 +126,40 @@ void convert_part_pos(const struct engine* e, const struct part* p,
}
}
+void convert_part_vel(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ 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;
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
@@ -130,33 +167,35 @@ void convert_part_pos(const struct engine* e, const struct part* p,
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
-void hydro_write_particles(struct part* parts, struct io_props* list,
- int* num_fields) {
+void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
+ struct io_props* list, int* num_fields) {
*num_fields = 10;
/* List what we want to write */
- list[0] = io_make_output_field_convert_part(
- "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
- list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts,
- primitives.v);
+ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+ UNIT_CONV_LENGTH, parts, xparts,
+ convert_part_pos);
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
+
list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts,
conserved.mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
parts, h);
list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
- parts, convert_u);
+ parts, xparts, convert_u);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
list[6] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
primitives.rho);
list[7] = io_make_output_field_convert_part(
- "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, convert_A);
+ "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, xparts, convert_A);
list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
parts, primitives.P);
list[9] = io_make_output_field_convert_part(
- "TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, convert_Etot);
+ "TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, xparts, convert_Etot);
}
/**
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 97cd5ef577dcfd7b07c5183a92349313096a6412..1070df105af19da546ed2814029a759713f4f1d0 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -197,6 +197,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
return p->mass;
}
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part *restrict p, float m) {
+
+ p->mass = m;
+}
+
/**
* @brief Returns the velocities drifted to the current time of a particle.
*
@@ -540,9 +552,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*
* @param p The particle to act upon
* @param xp The extended particle to act upon
+ * @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part *restrict p, struct xpart *restrict xp) {
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
@@ -580,4 +594,21 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
hydro_init_part(p, NULL);
}
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part *p, float u_init) {
+
+ p->u = u_init;
+}
+
#endif /* SWIFT_MINIMAL_HYDRO_H */
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index 2ff70fd2f7e95344c17ba0f3237c3b2ee131752c..e22cbd0f541ba1718c528434e6f7a839baa7d12f 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/src/hydro/Minimal/hydro_io.h
@@ -69,18 +69,20 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, rho);
}
-void convert_S(const struct engine* e, const struct part* p, float* ret) {
+void convert_S(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_entropy(p);
}
-void convert_P(const struct engine* e, const struct part* p, float* ret) {
+void convert_P(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_pressure(p);
}
void convert_part_pos(const struct engine* e, const struct part* p,
- double* ret) {
+ const struct xpart* xp, double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
@@ -93,23 +95,59 @@ void convert_part_pos(const struct engine* e, const struct part* p,
}
}
+void convert_part_vel(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ 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;
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
* @param parts The particle array.
+ * @param xparts The extended particle array.
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
-void hydro_write_particles(struct part* parts, struct io_props* list,
- int* num_fields) {
+void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
+ struct io_props* list, int* num_fields) {
*num_fields = 9;
/* List what we want to write */
- list[0] = io_make_output_field_convert_part(
- "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
- list[1] =
- io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
+ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+ UNIT_CONV_LENGTH, parts, xparts,
+ convert_part_pos);
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
list[2] =
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
@@ -120,10 +158,11 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
UNIT_CONV_NO_UNITS, parts, id);
list[6] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
- list[7] = io_make_output_field_convert_part(
- "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, convert_S);
+ list[7] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
+ UNIT_CONV_ENTROPY_PER_UNIT_MASS,
+ parts, xparts, convert_S);
list[8] = io_make_output_field_convert_part(
- "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
+ "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, xparts, convert_P);
}
/**
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index 37de855b375a5ac11f8a33f3993f7d2101c58522..78cc9c121700e9f92a0c2d1863e88a2174686857 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -170,6 +170,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
return p->mass;
}
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part *restrict p, float m) {
+
+ p->mass = m;
+}
+
/**
* @brief Returns the velocities drifted to the current time of a particle.
*
@@ -556,12 +568,16 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
* Requires the density to be known
*
* @param p The particle to act upon
+ * @param xp The extended data.
+ * @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part *restrict p, struct xpart *restrict xp) {
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
/* We read u in the entropy field. We now get S from u */
- xp->entropy_full = gas_entropy_from_internal_energy(p->rho_bar, p->entropy);
+ xp->entropy_full =
+ gas_entropy_from_internal_energy(p->rho_bar * cosmo->a3_inv, p->entropy);
p->entropy = xp->entropy_full;
p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
@@ -605,4 +621,21 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
hydro_init_part(p, NULL);
}
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part *p, float u_init) {
+
+ p->entropy = u_init;
+}
+
#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_H */
diff --git a/src/hydro/PressureEntropy/hydro_io.h b/src/hydro/PressureEntropy/hydro_io.h
index ed8e0e45d240ea37ebfb7df6339afd454dc0e4c4..6657e50a5bf071380dc4e33909e42521b5ce6037 100644
--- a/src/hydro/PressureEntropy/hydro_io.h
+++ b/src/hydro/PressureEntropy/hydro_io.h
@@ -67,18 +67,20 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, rho);
}
-void convert_u(const struct engine* e, const struct part* p, float* ret) {
+void convert_u(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_internal_energy(p);
}
-void convert_P(const struct engine* e, const struct part* p, float* ret) {
+void convert_P(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_pressure(p);
}
void convert_part_pos(const struct engine* e, const struct part* p,
- double* ret) {
+ const struct xpart* xp, double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
@@ -91,6 +93,40 @@ void convert_part_pos(const struct engine* e, const struct part* p,
}
}
+void convert_part_vel(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ 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;
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
@@ -98,16 +134,17 @@ void convert_part_pos(const struct engine* e, const struct part* p,
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
-void hydro_write_particles(struct part* parts, struct io_props* list,
- int* num_fields) {
+void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
+ struct io_props* list, int* num_fields) {
*num_fields = 10;
/* List what we want to write */
- list[0] = io_make_output_field_convert_part(
- "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
- list[1] =
- io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
+ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+ UNIT_CONV_LENGTH, parts, xparts,
+ convert_part_pos);
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
list[2] =
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
@@ -120,9 +157,9 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
list[7] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
- parts, convert_u);
+ parts, xparts, convert_u);
list[8] = io_make_output_field_convert_part(
- "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
+ "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, xparts, convert_P);
list[9] = io_make_output_field("WeightedDensity", FLOAT, 1, UNIT_CONV_DENSITY,
parts, rho_bar);
}
diff --git a/src/space.c b/src/space.c
index bbf97ea163cddd5fa7638931dad0e7a6e77bece3..ea599a45d3957c2ea224c0e7c3278f6ddfbe163e 100644
--- a/src/space.c
+++ b/src/space.c
@@ -2853,6 +2853,7 @@ void space_convert_quantities(struct space *s, int verbose) {
*
* @param s The #space to initialize.
* @param params The parsed parameter file.
+ * @param cosmo The current cosmological model.
* @param dim Spatial dimensions of the domain.
* @param parts Array of Gas particles.
* @param gparts Array of Gravity particles.
@@ -2863,7 +2864,7 @@ void space_convert_quantities(struct space *s, int verbose) {
* @param periodic flag whether the domain is periodic or not.
* @param replicate How many replications along each direction do we want?
* @param generate_gas_in_ics Are we generating gas particles from the gparts?
- * @param gravity flag whether we are doing gravity or not.
+ * @param self_gravity flag whether we are doing gravity or not?
* @param verbose Print messages to stdout or not.
* @param dry_run If 1, just initialise stuff, don't do anything with the parts.
*
@@ -3240,7 +3241,7 @@ void space_generate_gas(struct space *s, const struct cosmology *cosmo,
/* Set the masses */
gp_dm->mass *= (1. - mass_ratio);
gp_gas->mass *= mass_ratio;
- p->mass = gp_gas->mass;
+ hydro_set_mass(p, gp_gas->mass);
/* Set the new positions */
gp_dm->x[0] += shift_dm;