From 7da1281703d43f0c5d83d8a7ec98a9ebb155f840 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Tue, 13 Aug 2019 23:13:32 +0200
Subject: [PATCH] Backported the new interface to the chemistry from the
COLIBRE fork.
---
src/chemistry/EAGLE/chemistry.h | 105 +++++++++++++++++++++-
src/cooling/EAGLE/cooling.c | 16 ++--
src/cooling/EAGLE/cooling_rates.h | 70 +++++++--------
src/cooling/EAGLE/newton_cooling.c | 5 +-
src/feedback/EAGLE/feedback.c | 24 ++---
src/star_formation/EAGLE/star_formation.h | 12 ++-
6 files changed, 169 insertions(+), 63 deletions(-)
diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index a04d7f94e7..c06b29f441 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -170,9 +170,10 @@ __attribute__((always_inline)) INLINE static void chemistry_first_init_part(
p->chemistry_data.metal_mass_fraction_total =
data->initial_metal_mass_fraction_total;
- for (int elem = 0; elem < chemistry_element_count; ++elem)
+ for (int elem = 0; elem < chemistry_element_count; ++elem) {
p->chemistry_data.metal_mass_fraction[elem] =
data->initial_metal_mass_fraction[elem];
+ }
}
chemistry_init_part(p, data);
}
@@ -241,4 +242,106 @@ static INLINE void chemistry_print_backend(
chemistry_element_count);
}
+/**
+ * @brief Updates the metal mass fractions after diffusion at the end of the
+ * force loop.
+ *
+ * @param p The particle to act upon.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_end_force(
+ struct part* restrict p, const struct cosmology* cosmo) {}
+
+/**
+ * @brief Computes the chemistry-related time-step constraint.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param cosmo The current cosmological model.
+ * @param us The internal system of units.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float chemistry_timestep(
+ 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 chemistry_global_data* cd, const struct part* restrict p) {
+ return FLT_MAX;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * star particle to be used in feedback/enrichment related routines.
+ *
+ * EAGLE uses smooth abundances for everything.
+ *
+ * @param sp Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_total_metal_mass_fraction_for_feedback(
+ const struct spart* restrict sp) {
+
+ return sp->chemistry_data.smoothed_metal_mass_fraction_total;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * gas particle to be used in cooling related routines.
+ *
+ * EAGLE uses smooth abundances for everything.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_total_metal_mass_fraction_for_cooling(
+ const struct part* restrict p) {
+
+ return p->chemistry_data.smoothed_metal_mass_fraction_total;
+}
+
+/**
+ * @brief Returns the abundance array (metal mass fractions) of the
+ * gas particle to be used in cooling related routines.
+ *
+ * EAGLE uses smooth abundances for everything.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float const*
+chemistry_get_metal_mass_fraction_for_cooling(const struct part* restrict p) {
+
+ return p->chemistry_data.smoothed_metal_mass_fraction;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * gas particle to be used in star formation related routines.
+ *
+ * EAGLE uses smooth abundances for everything.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_total_metal_mass_fraction_for_star_formation(
+ const struct part* restrict p) {
+
+ return p->chemistry_data.smoothed_metal_mass_fraction_total;
+}
+
+/**
+ * @brief Returns the abundance array (metal mass fractions) of the
+ * gas particle to be used in star formation related routines.
+ *
+ * EAGLE uses smooth abundances for everything.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float const*
+chemistry_get_metal_mass_fraction_for_star_formation(
+ const struct part* restrict p) {
+
+ return p->chemistry_data.smoothed_metal_mass_fraction;
+}
+
#endif /* SWIFT_CHEMISTRY_EAGLE_H */
diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
index 674631ae52..c5595ef514 100644
--- a/src/cooling/EAGLE/cooling.c
+++ b/src/cooling/EAGLE/cooling.c
@@ -423,10 +423,10 @@ void cooling_cool_part(const struct phys_const *phys_const,
abundance_ratio_to_solar(p, cooling, abundance_ratio);
/* Get the Hydrogen and Helium mass fractions */
- const float XH =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H];
- const float XHe =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_He];
+ const float *const metal_fraction =
+ chemistry_get_metal_mass_fraction_for_cooling(p);
+ const float XH = metal_fraction[chemistry_element_H];
+ const float XHe = metal_fraction[chemistry_element_He];
/* Get the Helium mass fraction. Note that this is He / (H + He), i.e. a
* metal-free Helium mass fraction as per the Wiersma+08 definition */
@@ -600,10 +600,10 @@ float cooling_get_temperature(
const double u_cgs = u * cooling->internal_energy_to_cgs;
/* Get the Hydrogen and Helium mass fractions */
- const float XH =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H];
- const float XHe =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_He];
+ const float *const metal_fraction =
+ chemistry_get_metal_mass_fraction_for_cooling(p);
+ const float XH = metal_fraction[chemistry_element_H];
+ const float XHe = metal_fraction[chemistry_element_He];
/* Get the Helium mass fraction. Note that this is He / (H + He), i.e. a
* metal-free Helium mass fraction as per the Wiersma+08 definition */
diff --git a/src/cooling/EAGLE/cooling_rates.h b/src/cooling/EAGLE/cooling_rates.h
index bf908c3a7c..39978c9bfa 100644
--- a/src/cooling/EAGLE/cooling_rates.h
+++ b/src/cooling/EAGLE/cooling_rates.h
@@ -43,61 +43,57 @@
* We also re-order the elements such that they match the order of the
* tables. This is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe].
*
+ * The solar abundances table (from the cooling struct) is arranged as
+ * [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe].
+ *
* @param p Pointer to #part struct.
* @param cooling #cooling_function_data struct.
* @param ratio_solar (return) Array of ratios to solar abundances.
*/
-__attribute__((always_inline)) INLINE void abundance_ratio_to_solar(
+__attribute__((always_inline)) INLINE static void abundance_ratio_to_solar(
const struct part *p, const struct cooling_function_data *cooling,
float ratio_solar[eagle_cooling_N_abundances]) {
- ratio_solar[0] =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H] *
- cooling->SolarAbundances_inv[0 /* H */];
+ /* Get the individual metal mass fractions from the particle */
+ const float *const metal_fraction =
+ chemistry_get_metal_mass_fraction_for_cooling(p);
+
+ ratio_solar[0] = metal_fraction[chemistry_element_H] *
+ cooling->SolarAbundances_inv[0 /* H */];
- ratio_solar[1] =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_He] *
- cooling->SolarAbundances_inv[1 /* He */];
+ ratio_solar[1] = metal_fraction[chemistry_element_He] *
+ cooling->SolarAbundances_inv[1 /* He */];
- ratio_solar[2] =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_C] *
- cooling->SolarAbundances_inv[2 /* C */];
+ ratio_solar[2] = metal_fraction[chemistry_element_C] *
+ cooling->SolarAbundances_inv[2 /* C */];
- ratio_solar[3] =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_N] *
- cooling->SolarAbundances_inv[3 /* N */];
+ ratio_solar[3] = metal_fraction[chemistry_element_N] *
+ cooling->SolarAbundances_inv[3 /* N */];
- ratio_solar[4] =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_O] *
- cooling->SolarAbundances_inv[4 /* O */];
+ ratio_solar[4] = metal_fraction[chemistry_element_O] *
+ cooling->SolarAbundances_inv[4 /* O */];
- ratio_solar[5] =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Ne] *
- cooling->SolarAbundances_inv[5 /* Ne */];
+ ratio_solar[5] = metal_fraction[chemistry_element_Ne] *
+ cooling->SolarAbundances_inv[5 /* Ne */];
- ratio_solar[6] =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Mg] *
- cooling->SolarAbundances_inv[6 /* Mg */];
+ ratio_solar[6] = metal_fraction[chemistry_element_Mg] *
+ cooling->SolarAbundances_inv[6 /* Mg */];
- ratio_solar[7] =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Si] *
- cooling->SolarAbundances_inv[7 /* Si */];
+ ratio_solar[7] = metal_fraction[chemistry_element_Si] *
+ cooling->SolarAbundances_inv[7 /* Si */];
/* For S, we use the same ratio as Si */
- ratio_solar[8] =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Si] *
- cooling->SolarAbundances_inv[7 /* Si */] *
- cooling->S_over_Si_ratio_in_solar;
+ ratio_solar[8] = metal_fraction[chemistry_element_Si] *
+ cooling->SolarAbundances_inv[7 /* Si */] *
+ cooling->S_over_Si_ratio_in_solar;
/* For Ca, we use the same ratio as Si */
- ratio_solar[9] =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Si] *
- cooling->SolarAbundances_inv[7 /* Si */] *
- cooling->Ca_over_Si_ratio_in_solar;
-
- ratio_solar[10] =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Fe] *
- cooling->SolarAbundances_inv[10 /* Fe */];
+ ratio_solar[9] = metal_fraction[chemistry_element_Si] *
+ cooling->SolarAbundances_inv[7 /* Si */] *
+ cooling->Ca_over_Si_ratio_in_solar;
+
+ ratio_solar[10] = metal_fraction[chemistry_element_Fe] *
+ cooling->SolarAbundances_inv[10 /* Fe */];
}
/**
diff --git a/src/cooling/EAGLE/newton_cooling.c b/src/cooling/EAGLE/newton_cooling.c
index e86ebb7a2f..50af3294c3 100644
--- a/src/cooling/EAGLE/newton_cooling.c
+++ b/src/cooling/EAGLE/newton_cooling.c
@@ -621,8 +621,9 @@ INLINE static double eagle_metal_cooling_rate(
const float log_table_bound_low = (cooling->Therm[0] + 0.05) / M_LOG10E;
/* convert Hydrogen mass fraction in Hydrogen number density */
- const float XH =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H];
+ float const *metal_fraction =
+ chemistry_get_metal_mass_fraction_for_cooling(p);
+ const float XH = metal_fraction[chemistry_element_H];
const double n_H = hydro_get_physical_density(p, cosmo) * XH /
phys_const->const_proton_mass;
const double n_H_cgs = n_H * cooling->number_density_to_cgs;
diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c
index cf6848760f..69f92fe992 100644
--- a/src/feedback/EAGLE/feedback.c
+++ b/src/feedback/EAGLE/feedback.c
@@ -105,15 +105,15 @@ double eagle_feedback_energy_fraction(const struct spart* sp,
/* Star properties */
- /* Smoothed metallicity (metal mass fraction) at birth time of the star */
- const double Z_smooth = sp->chemistry_data.smoothed_metal_mass_fraction_total;
+ /* Metallicity (metal mass fraction) at birth time of the star */
+ const double Z = chemistry_get_total_metal_mass_fraction_for_feedback(sp);
/* Physical density of the gas at the star's birth time */
const double rho_birth = sp->birth_density;
const double n_birth = rho_birth * props->rho_to_n_cgs;
/* Calculate f_E */
- const double Z_term = pow(max(Z_smooth, 1e-6) / Z_0, n_Z);
+ const double Z_term = pow(max(Z, 1e-6) / Z_0, n_Z);
const double n_term = pow(n_birth / n_0, -n_n);
const double denonimator = 1. + Z_term * n_term;
@@ -333,7 +333,7 @@ INLINE static void evolve_SNIa(const float log10_min_mass,
* and use updated values for the star's age and timestep in this function */
if (log10_max_mass > props->log10_SNIa_max_mass_msun) {
- const float Z = sp->chemistry_data.metal_mass_fraction_total;
+ const float Z = chemistry_get_total_metal_mass_fraction_for_feedback(sp);
const float max_mass = exp10f(props->log10_SNIa_max_mass_msun);
const float lifetime_Gyr = lifetime_in_Gyr(max_mass, Z, props);
@@ -401,6 +401,9 @@ INLINE static void evolve_SNII(float log10_min_mass, float log10_max_mass,
int low_imf_mass_bin_index, high_imf_mass_bin_index, mass_bin_index;
+ /* Metallicity (metal mass fraction) at birth time of the star */
+ const double Z = chemistry_get_total_metal_mass_fraction_for_feedback(sp);
+
/* If mass at beginning of step is less than tabulated lower bound for IMF,
* limit it.*/
if (log10_min_mass < props->log10_SNII_min_mass_msun)
@@ -423,9 +426,7 @@ INLINE static void evolve_SNII(float log10_min_mass, float log10_max_mass,
int iz_low = 0, iz_high = 0, low_index_3d, high_index_3d, low_index_2d,
high_index_2d;
float dz = 0.;
- determine_bin_yield_SNII(&iz_low, &iz_high, &dz,
- log10(sp->chemistry_data.metal_mass_fraction_total),
- props);
+ determine_bin_yield_SNII(&iz_low, &iz_high, &dz, log10(Z), props);
/* compute metals produced */
float metal_mass_released[chemistry_element_count], metal_mass_released_total;
@@ -557,6 +558,9 @@ INLINE static void evolve_AGB(const float log10_min_mass, float log10_max_mass,
int low_imf_mass_bin_index, high_imf_mass_bin_index, mass_bin_index;
+ /* Metallicity (metal mass fraction) at birth time of the star */
+ const double Z = chemistry_get_total_metal_mass_fraction_for_feedback(sp);
+
/* If mass at end of step is greater than tabulated lower bound for IMF, limit
* it.*/
if (log10_max_mass > props->log10_SNII_min_mass_msun)
@@ -574,9 +578,7 @@ INLINE static void evolve_AGB(const float log10_min_mass, float log10_max_mass,
int iz_low = 0, iz_high = 0, low_index_3d, high_index_3d, low_index_2d,
high_index_2d;
float dz = 0.f;
- determine_bin_yield_AGB(&iz_low, &iz_high, &dz,
- log10(sp->chemistry_data.metal_mass_fraction_total),
- props);
+ determine_bin_yield_AGB(&iz_low, &iz_high, &dz, log10(Z), props);
/* compute metals produced */
float metal_mass_released[chemistry_element_count], metal_mass_released_total;
@@ -718,7 +720,7 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props,
const double star_age_Gyr = age * conversion_factor;
/* Get the metallicity */
- const float Z = sp->chemistry_data.metal_mass_fraction_total;
+ const float Z = chemistry_get_total_metal_mass_fraction_for_feedback(sp);
/* Properties collected in the stellar density loop. */
const float ngb_gas_mass = sp->feedback_data.to_collect.ngb_mass;
diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h
index f2c77e0368..3533abcae6 100644
--- a/src/star_formation/EAGLE/star_formation.h
+++ b/src/star_formation/EAGLE/star_formation.h
@@ -232,8 +232,11 @@ INLINE static int star_formation_is_star_forming(
* because we also need to check if the physical density exceeded
* the appropriate limit */
- const double Z = p->chemistry_data.smoothed_metal_mass_fraction_total;
- const double X_H = p->chemistry_data.smoothed_metal_mass_fraction[0];
+ const double Z =
+ chemistry_get_total_metal_mass_fraction_for_star_formation(p);
+ const float* const metal_fraction =
+ chemistry_get_metal_mass_fraction_for_star_formation(p);
+ const double X_H = metal_fraction[chemistry_element_H];
const double n_H = physical_density * X_H;
/* Get the density threshold */
@@ -278,8 +281,9 @@ INLINE static void star_formation_compute_SFR(
}
/* Hydrogen number density of this particle */
- const double physical_density = hydro_get_physical_density(p, cosmo);
- const double X_H = p->chemistry_data.smoothed_metal_mass_fraction[0];
+ const float* const metal_fraction =
+ chemistry_get_metal_mass_fraction_for_star_formation(p);
+ const double X_H = metal_fraction[chemistry_element_H];
const double n_H = physical_density * X_H / phys_const->const_proton_mass;
/* Are we above the threshold for automatic star formation? */
--
GitLab