diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c index 72752df439e27448974265ce6d4585d1ee7b60bf..c5f24a4760efae065ba27850e002e1d1e762fe3e 100644 --- a/src/cooling/EAGLE/cooling.c +++ b/src/cooling/EAGLE/cooling.c @@ -790,7 +790,15 @@ void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling, /* Loop through particles and set new heat */ for (size_t i = 0; i < s->nr_parts; i++) { - hydro_reion_heating(&parts[i], &xparts[i], cosmo, extra_heat); + + struct part *p = &parts[i]; + struct xpart *xp = &xparts[i]; + + const float old_u = hydro_get_physical_internal_energy(p, xp, cosmo); + const float new_u = old_u + extra_heat; + + hydro_set_physical_internal_energy(p, xp, cosmo, new_u); + hydro_set_drifted_physical_internal_energy(p, cosmo, new_u); } } diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index 95e0c7c9df800193ffef977092b2567c068d2192..bc505392e96c48d282b3543c8d47a87202d182d9 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -347,6 +347,57 @@ __attribute__((always_inline)) INLINE static void hydro_set_physical_entropy( xp->entropy_full = entropy; } +/** + * @brief Sets the physical internal energy of a particle + * + * @param p The particle of interest. + * @param xp The extended particle data. + * @param cosmo Cosmology data structure + * @param u The physical entropy + */ +__attribute__((always_inline)) INLINE static void +hydro_set_physical_internal_energy(struct part *p, struct xpart *xp, + const struct cosmology *cosmo, + const float u) { + + xp->entropy_full = + gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, u); +} + +/** + * @brief Sets the drifted physical internal energy of a particle + * + * @param p The particle of interest. + * @param xp The extended particle data. + * @param cosmo Cosmology data structure + * @param u The physical entropy + */ +__attribute__((always_inline)) INLINE static void +hydro_set_drifted_physical_internal_energy(struct part *p, + const struct cosmology *cosmo, + const float u) { + + p->entropy = gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, u); + + /* Now recompute the extra quantities */ + + /* Inverse of the co-moving density */ + const float rho_inv = 1.f / p->rho; + + /* Compute the pressure */ + const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); + + /* Compute the sound speed */ + const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + + /* Divide the pressure by the density squared to get the SPH term */ + const float P_over_rho2 = pressure * rho_inv * rho_inv; + + /* Update variables. */ + p->force.P_over_rho2 = P_over_rho2; + p->force.soundspeed = soundspeed; +} + /** * @brief Computes the hydro time-step of a given particle * @@ -806,25 +857,4 @@ hydro_set_init_internal_energy(struct part *p, float u_init) { __attribute__((always_inline)) INLINE static void hydro_remove_part( const struct part *p, const struct xpart *xp) {} -/** - * @brief Inputs extra heat to a particle at redshift of reionization - * - * We assume a constant density. - * - * @param p The particle of interest. - * @param xp The extended particle data - * @param cosmo Cosmology data structure - * @param extra_heat The extra internal energy given to the particle. - */ -__attribute__((always_inline)) INLINE static void hydro_reion_heating( - struct part *p, struct xpart *xp, const struct cosmology *cosmo, - float extra_heat) { - - const float old_u = gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, - xp->entropy_full); - const float new_u = old_u + extra_heat; - xp->entropy_full = - gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, new_u); -} - #endif /* SWIFT_GADGET2_HYDRO_H */