Commit 7da12817 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Backported the new interface to the chemistry from the COLIBRE fork.

parent c009fbb9
......@@ -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 */
......@@ -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 */
......
......@@ -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 */];
}
/**
......
......@@ -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;
......
......@@ -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;
......
......@@ -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? */
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment