From 3b84d63f998ab42518dddb2dffdcb29b43af1d7a Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Sun, 7 Oct 2018 16:27:14 +0200
Subject: [PATCH] Modify the particle getters to distinguish between kicked and
drifted thermal quantities.
---
doc/Doxyfile.in | 2 +
.../plotTempEvolution.py | 2 +-
src/cooling/const_lambda/cooling.h | 10 ++-
src/cooling/const_lambda/cooling_io.h | 3 +-
src/hydro/Gadget2/hydro.h | 89 +++++++++++++++----
src/hydro/Gadget2/hydro_io.h | 3 +-
src/statistics.c | 4 +-
src/timestep.h | 2 +-
8 files changed, 88 insertions(+), 27 deletions(-)
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index cba52250cc..dada62acab 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -761,11 +761,13 @@ WARN_LOGFILE =
INPUT = @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_srcdir@/examples
INPUT += @top_srcdir@/src/hydro/Minimal
+INPUT += @top_srcdir@/src/hydro/Gadget2
INPUT += @top_srcdir@/src/gravity/Default
INPUT += @top_srcdir@/src/stars/Default
INPUT += @top_srcdir@/src/riemann
INPUT += @top_srcdir@/src/potential/point_mass
INPUT += @top_srcdir@/src/equation_of_state/ideal_gas
+INPUT += @top_srcdir@/src/cooling/const_lambda
INPUT += @top_srcdir@/src/cooling/EAGLE
INPUT += @top_srcdir@/src/chemistry/EAGLE
diff --git a/examples/SmallCosmoVolume_cooling/plotTempEvolution.py b/examples/SmallCosmoVolume_cooling/plotTempEvolution.py
index 526a417dc6..6f7f3bbe76 100644
--- a/examples/SmallCosmoVolume_cooling/plotTempEvolution.py
+++ b/examples/SmallCosmoVolume_cooling/plotTempEvolution.py
@@ -24,7 +24,7 @@ k_in_J_K = 1.38064852e-23
mH_in_kg = 1.6737236e-27
# Number of snapshots generated
-n_snapshots = 160
+n_snapshots = 200
import matplotlib
matplotlib.use("Agg")
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index c1f5d802e4..e38b6a602e 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/src/cooling/const_lambda/cooling.h
@@ -98,7 +98,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
const float u_floor = hydro_props->minimal_internal_energy;
/* Current energy */
- const float u_old = hydro_get_physical_internal_energy2(p, xp, cosmo);
+ const float u_old = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Current du_dt in physical coordinates */
const float hydro_du_dt = hydro_get_physical_internal_energy_dt(p, cosmo);
@@ -151,13 +151,15 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
* @param hydro_props The properties of the hydro scheme.
* @param us The internal system of units.
* @param p Pointer to the particle data.
+ * @param xp Pointer to the extended data of the particle.
*/
__attribute__((always_inline)) INLINE static float cooling_timestep(
const struct cooling_function_data* restrict cooling,
const struct phys_const* restrict phys_const,
const struct cosmology* restrict cosmo,
const struct unit_system* restrict us,
- const struct hydro_props* hydro_props, const struct part* restrict p) {
+ const struct hydro_props* hydro_props, const struct part* restrict p,
+ const struct xpart* restrict xp) {
/* Start with the case where there is no limit */
if (cooling->cooling_tstep_mult == FLT_MAX) {
@@ -165,7 +167,7 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
}
/* Get current internal energy and cooling rate */
- const float u = hydro_get_physical_internal_energy(p, cosmo);
+ const float u = hydro_get_physical_internal_energy(p, xp, cosmo);
const double cooling_du_dt_cgs =
cooling_rate_cgs(cosmo, hydro_props, cooling, p);
@@ -174,7 +176,7 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
cooling_du_dt_cgs * cooling->conv_factor_energy_rate_from_cgs;
/* If we are close to (or below) the limit, we ignore the condition */
- if (u < 1.01f * hydro_props->minimal_internal_energy)
+ if (u < 1.01f * hydro_props->minimal_internal_energy || cooling_du_dt == 0.f)
return FLT_MAX;
else
return cooling->cooling_tstep_mult * u / fabsf(cooling_du_dt);
diff --git a/src/cooling/const_lambda/cooling_io.h b/src/cooling/const_lambda/cooling_io.h
index e26949e62e..bcc675a3a6 100644
--- a/src/cooling/const_lambda/cooling_io.h
+++ b/src/cooling/const_lambda/cooling_io.h
@@ -32,7 +32,8 @@
/**
* @brief Writes the current model of SPH to the file
- * @param h_grpsph The HDF5 group in which to write
+ * @param h_grp The HDF5 group in which to write
+ * @param cooling the parameters of the cooling function.
*/
__attribute__((always_inline)) INLINE static void cooling_write_flavour(
hid_t h_grp, const struct cooling_function_data* cooling) {
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 1b4b1d27e7..36db326cf9 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -42,42 +42,60 @@
#include "minmax.h"
/**
- * @brief Returns the comoving internal energy of a particle
+ * @brief Returns the comoving internal energy of a particle at the last
+ * time the particle was kicked.
*
* @param p The particle of interest
+ * @param xp The extended data of the particle of interest.
*/
__attribute__((always_inline)) INLINE static float
-hydro_get_comoving_internal_energy(const struct part *restrict p) {
+hydro_get_comoving_internal_energy(const struct part *restrict p,
+ const struct xpart *restrict xp) {
- return gas_internal_energy_from_entropy(p->rho, p->entropy);
+ return gas_internal_energy_from_entropy(p->rho, xp->entropy_full);
}
/**
- * @brief Returns the physical internal energy of a particle
+ * @brief Returns the physical internal energy of a particle at the last
+ * time the particle was kicked.
*
* @param p The particle of interest.
+ * @param xp The extended data of the particle of interest.
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float
hydro_get_physical_internal_energy(const struct part *restrict p,
+ const struct xpart *restrict xp,
const struct cosmology *cosmo) {
- return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
+ return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv,
+ xp->entropy_full);
}
/**
- * @brief Returns the physical internal energy of a particle
+ * @brief Returns the comoving internal energy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_comoving_internal_energy(const struct part *restrict p) {
+
+ return gas_internal_energy_from_entropy(p->rho, p->entropy);
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle drifted to the
+ * current time.
*
* @param p The particle of interest.
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float
-hydro_get_physical_internal_energy2(const struct part *restrict p,
- const struct xpart *restrict xp,
- const struct cosmology *cosmo) {
+hydro_get_drifted_physical_internal_energy(const struct part *restrict p,
+ const struct cosmology *cosmo) {
- return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv,
- xp->entropy_full);
+ return gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
}
/**
@@ -94,7 +112,8 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
/**
* @brief Returns the physical pressure of a particle
*
- * @param p The particle of interest
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
const struct part *restrict p, const struct cosmology *cosmo) {
@@ -103,24 +122,57 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
}
/**
- * @brief Returns the comoving entropy of a particle
+ * @brief Returns the comoving entropy of a particle at the last
+ * time the particle was kicked.
*
* @param p The particle of interest.
+ * @param xp The extended data of the particle of interest.
*/
__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
- const struct part *restrict p) {
+ const struct part *restrict p, const struct xpart *restrict xp) {
- return p->entropy;
+ return xp->entropy_full;
}
/**
- * @brief Returns the physical entropy of a particle
+ * @brief Returns the physical entropy of a particl at the last
+ * time the particle was kicked.
*
* @param p The particle of interest.
+ * @param xp The extended data of the particle of interest.
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
- const struct part *restrict p, const struct cosmology *cosmo) {
+ const struct part *restrict p, const struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
+
+ /* Note: no cosmological conversion required here with our choice of
+ * coordinates. */
+ return xp->entropy_full;
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_comoving_entropy(const struct part *restrict p) {
+
+ return p->entropy;
+}
+
+/**
+ * @brief Returns the physical entropy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_physical_entropy(const struct part *restrict p,
+ const struct cosmology *cosmo) {
/* Note: no cosmological conversion required here with our choice of
* coordinates. */
@@ -581,6 +633,9 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param p The particle to act upon
* @param xp The particle extended data to act upon
* @param dt_therm The time-step for this kick (for thermodynamic quantities)
+ * @param dt_grav The time-step for this kick (for gravity forces)
+ * @param dt_hydro The time-step for this kick (for hydro forces)
+ * @param dt_kick_corr The time-step for this kick (for correction of the kick)
* @param cosmo The cosmological model.
* @param hydro_props The constants used in the scheme
*/
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 3f2af41dc7..27aab4b897 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -59,7 +59,7 @@ INLINE static void hydro_read_particles(struct part* parts,
INLINE static void convert_part_u(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
- ret[0] = hydro_get_comoving_internal_energy(p);
+ ret[0] = hydro_get_comoving_internal_energy(p, xp);
}
INLINE static void convert_part_P(const struct engine* e, const struct part* p,
@@ -132,6 +132,7 @@ INLINE static void convert_part_potential(const struct engine* e,
* @brief Specifies which particle fields to write to a dataset
*
* @param parts The particle array.
+ * @param xparts The extended particle data array.
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
diff --git a/src/statistics.c b/src/statistics.c
index 22ddc2e971..47365dd9e3 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -166,8 +166,8 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, v);
const double x[3] = {p->x[0], p->x[1], p->x[2]};
const float m = hydro_get_mass(p);
- const float entropy = hydro_get_physical_entropy(p, cosmo);
- const float u_inter = hydro_get_physical_internal_energy(p, cosmo);
+ const float entropy = hydro_get_drifted_physical_entropy(p, cosmo);
+ const float u_inter = hydro_get_drifted_physical_internal_energy(p, cosmo);
/* Collect mass */
stats.mass += m;
diff --git a/src/timestep.h b/src/timestep.h
index ae90213c1e..e9943a41a0 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -128,7 +128,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
if (e->policy & engine_policy_cooling)
new_dt_cooling =
cooling_timestep(e->cooling_func, e->physical_constants, e->cosmology,
- e->internal_units, e->hydro_properties, p);
+ e->internal_units, e->hydro_properties, p, xp);
/* Compute the next timestep (gravity condition) */
float new_dt_grav = FLT_MAX, new_dt_self_grav = FLT_MAX,
--
GitLab