diff --git a/.gitignore b/.gitignore
index d8badf0e0bf1e7fa785065aa9480e87dabcad0b3..163d09fd929ab0a4e414915592c68e0293c11eb8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -32,6 +32,7 @@ examples/*/*/*.dat
examples/*/*/*.png
examples/*/*/*.pdf
examples/*/*/*.mp4
+examples/*/*/*.txt
examples/*/*/*.rst
examples/*/*/*.hdf5
examples/*/*/*.csv
diff --git a/doc/RTD/source/SubgridModels/EAGLE/index.rst b/doc/RTD/source/SubgridModels/EAGLE/index.rst
index d9eb19eb7fd41b016443d2cef1edbdbbf8d314b7..3cbfbaccd9cea682bde42f7b615894d023be1075 100644
--- a/doc/RTD/source/SubgridModels/EAGLE/index.rst
+++ b/doc/RTD/source/SubgridModels/EAGLE/index.rst
@@ -666,6 +666,8 @@ Supernova feedback: Dalla Vecchia+2012 & Schaye+2015
SNIa_DTD_exp_timescale_Gyr: 2.0 # Time-scale of the SNIa delay time distribution in the exponential DTD case (in Gyr).
SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ stellar_evolution_age_cut_Gyr: 0.1 # Stellar age in Gyr above which the enrichment is down-sampled.
+ stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
diff --git a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
index 487672ba8ff84168450acd9dbddef4c42f93accc..452c1f84a4cefb4bb907a4dc651cee12f3ed7359 100644
--- a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
@@ -162,6 +162,8 @@ EAGLEFeedback:
SNIa_DTD_power_law_norm_p_Msun: 0.0012 # Normalization of the SNIa delay time distribution (in Msun^-1).
SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ stellar_evolution_age_cut_Gyr: 0.1 # Stellar age in Gyr above which the enrichment is down-sampled.
+ stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
diff --git a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
index 5595fbe78b5bceae983b6d38a0ee30b3f631a359..43034a08a094098320ebe79e482d9a2d1fcdb929 100644
--- a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
@@ -163,6 +163,8 @@ EAGLEFeedback:
SNIa_DTD_power_law_norm_p_Msun: 0.0012 # Normalization of the SNIa delay time distribution (in Msun^-1).
SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ stellar_evolution_age_cut_Gyr: 0.1 # Stellar age in Gyr above which the enrichment is down-sampled.
+ stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
diff --git a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
index a3b518cbac8f076e416d4d659cd3466b89af9b8b..abc092e262cc168c4f74c73a66ae2183e4128722 100644
--- a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
@@ -163,6 +163,8 @@ EAGLEFeedback:
SNIa_DTD_power_law_norm_p_Msun: 0.0012 # Normalization of the SNIa delay time distribution (in Msun^-1).
SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ stellar_evolution_age_cut_Gyr: 0.1 # Stellar age in Gyr above which the enrichment is down-sampled.
+ stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index 40d4aeaec0fc621b40081e810ce22696fcbed8f8..3caab833bebb0f4d7baf306f2acb58a94c9d2fd0 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -157,6 +157,8 @@ EAGLEFeedback:
SNIa_DTD_power_law_norm_p_Msun: 0.0012 # Normalization of the SNIa delay time distribution (in Msun^-1).
SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ stellar_evolution_age_cut_Gyr: 0.1 # Stellar age in Gyr above which the enrichment is down-sampled.
+ stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index 6969a993816263b47c744bdcc84b4dba09e8de71..2c955661478297c741dc0779190252d8b81c7929 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -165,6 +165,8 @@ EAGLEFeedback:
SNIa_DTD_power_law_norm_p_Msun: 0.0012 # Normalization of the SNIa delay time distribution (in Msun^-1).
SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ stellar_evolution_age_cut_Gyr: 0.1 # Stellar age in Gyr above which the enrichment is down-sampled.
+ stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
index 3c709c3875168a015358b17fa8772c3296b968bf..49c631a3d931eaeeaf3bbb63593c69afd8522f62 100644
--- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
@@ -156,6 +156,8 @@ EAGLEFeedback:
SNIa_DTD_power_law_norm_p_Msun: 0.0012 # Normalization of the SNIa delay time distribution (in Msun^-1).
SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ stellar_evolution_age_cut_Gyr: 0.1 # Stellar age in Gyr above which the enrichment is down-sampled.
+ stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index cae4b4313f4d34726115e1060e9f20ff2ba303fa..7b9608b6021bee3687b3585a3a9b02c8773a14cc 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -166,6 +166,8 @@ EAGLEFeedback:
SNIa_DTD_power_law_norm_p_Msun: 0.0012 # Normalization of the SNIa delay time distribution (in Msun^-1).
SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ stellar_evolution_age_cut_Gyr: 0.1 # Stellar age in Gyr above which the enrichment is down-sampled.
+ stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
index 94a1b5a8fb29a46971bf1f025b545e47e2faf9f4..0642e9312f0926cc0d005f9ec4b6d3f5215ea7d4 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
@@ -141,6 +141,8 @@ EAGLEFeedback:
SNIa_DTD_power_law_norm_p_Msun: 0.0012 # Normalization of the SNIa delay time distribution (in Msun^-1).
SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ stellar_evolution_age_cut_Gyr: 0.1 # Age in Gyr above which the enrichment is downsampled.
+ stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
diff --git a/examples/SubgridTests/StellarEvolution/plot_box_evolution.py b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
index a4f601875beb40c8f264cda385815f2d6dd81a41..aa1f07f98ecede241423bc72a3e0cf7d76963979 100644
--- a/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
+++ b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
@@ -107,12 +107,12 @@ t = zeros(n_snapshots)
# Read data from snapshots
for i in range(n_snapshots):
- #print("reading snapshot "+str(i))
- sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
- t[i] = sim["/Header"].attrs["Time"][0]
+ #print("reading snapshot "+str(i))
+ sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
+ t[i] = sim["/Header"].attrs["Time"][0]
- masses = sim["/PartType0/Masses"][:]
- swift_box_gas_mass[i] = np.sum(masses)
+ masses = sim["/PartType0/Masses"][:]
+ swift_box_gas_mass[i] = np.sum(masses)
AGB_mass = sim["/PartType0/MassesFromAGB"][:]
SNII_mass = sim["/PartType0/MassesFromSNII"][:]
@@ -121,13 +121,13 @@ for i in range(n_snapshots):
swift_box_gas_mass_AGB[i] = np.sum(AGB_mass)
swift_box_gas_mass_SNII[i] = np.sum(SNII_mass)
swift_box_gas_mass_SNIa[i] = np.sum(SNIa_mass)
-
+
Z_star = sim["/PartType4/MetalMassFractions"][0]
- star_masses = sim["/PartType4/Masses"][:]
- swift_box_star_mass[i] = np.sum(star_masses)
+ star_masses = sim["/PartType4/Masses"][:]
+ swift_box_star_mass[i] = np.sum(star_masses)
- metallicities = sim["/PartType0/MetalMassFractions"][:]
- swift_box_gas_metal_mass[i] = np.sum(metallicities * masses)
+ metallicities = sim["/PartType0/MetalMassFractions"][:]
+ swift_box_gas_metal_mass[i] = np.sum(metallicities * masses)
AGB_Z_frac = sim["/PartType0/MetalMassFractionsFromAGB"][:]
SNII_Z_frac = sim["/PartType0/MetalMassFractionsFromSNII"][:]
@@ -136,10 +136,10 @@ for i in range(n_snapshots):
swift_box_gas_metal_mass_AGB[i] = np.sum(AGB_Z_frac * masses)
swift_box_gas_metal_mass_SNII[i] = np.sum(SNII_Z_frac * masses)
swift_box_gas_metal_mass_SNIa[i] = np.sum(SNIa_Z_frac * masses)
-
- element_abundances = sim["/PartType0/ElementMassFractions"][:][:]
- for j in range(n_elements):
- swift_element_mass[i,j] = np.sum(element_abundances[:,j] * masses)
+
+ element_abundances = sim["/PartType0/ElementMassFractions"][:][:]
+ for j in range(n_elements):
+ swift_element_mass[i,j] = np.sum(element_abundances[:,j] * masses)
v = sim["/PartType0/Velocities"][:,:]
v2 = v[:,0]**2 + v[:,1]**2 + v[:,2]**2
@@ -218,8 +218,8 @@ colours = ['k','r','g','b','c','y','m','skyblue','plum']
element_names = ['H','He','C','N','O','Ne','Mg','Si','Fe']
subplot(223)
for j in range(n_elements):
- plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (swift_element_mass[1:,j] - swift_element_mass[0,j]) * unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color=colours[j], ms=0.5, label=element_names[j])
- plot(eagle_time_Gyr[1:],eagle_total_element_mass[:-1,j],linewidth=1,color=colours[j],linestyle='--')
+ plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (swift_element_mass[1:,j] - swift_element_mass[0,j]) * unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color=colours[j], ms=0.5, label=element_names[j])
+ plot(eagle_time_Gyr[1:],eagle_total_element_mass[:-1,j],linewidth=1,color=colours[j],linestyle='--')
xlabel("${\\rm Time~[Gyr]}$", labelpad=0)
ylabel("Change in element mass of gas particles ${[\\rm M_\\odot]}$", labelpad=2)
xscale("log")
diff --git a/examples/SubgridTests/StellarEvolution/stellar_evolution.yml b/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
index 52b212de9e72e58a95f7f72a4796cc673a9dd4b5..fb0b385bc9345920e478846f2bfccc07416a2932 100644
--- a/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
+++ b/examples/SubgridTests/StellarEvolution/stellar_evolution.yml
@@ -16,8 +16,8 @@ TimeIntegration:
# Parameters governing the snapshots
Snapshots:
basename: stellar_evolution # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
- delta_time: 3.e-5 # Time difference between consecutive outputs without cosmology (internal units)
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 3.e-5 # Time difference between consecutive outputs (internal units)
compression: 4
# Parameters governing the conserved quantities statistics
@@ -106,6 +106,8 @@ EAGLEFeedback:
SNIa_DTD_exp_norm_p_Msun: 0.002 # Normalisation of the SNIa rates in inverse solar masses.
SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ stellar_evolution_age_cut_Gyr: 0.1 # Age in Gyr above which the enrichment is downsampled.
+ stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index a605a804681770530bf77bfefffdf02655f6ebe8..7e753c4a1cf06f4dc28f9e80534df63be1565cb4 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -412,6 +412,8 @@ EAGLEFeedback:
SNIa_DTD_exp_timescale_Gyr: 2.0 # Time-scale of the SNIa delay time distribution in the exponential DTD case (in Gyr).
SNIa_energy_erg: 1.0e51 # Energy of one SNIa explosion in ergs.
AGB_ejecta_velocity_km_p_s: 10.0 # Velocity of the AGB ejectas in km/s.
+ stellar_evolution_age_cut_Gyr: 0.1 # Stellar age in Gyr above which the enrichment is down-sampled.
+ stellar_evolution_sampling_rate: 10 # Number of time-steps in-between two enrichment events for a star above the age threshold.
SNII_yield_factor_Hydrogen: 1.0 # (Optional) Correction factor to apply to the Hydrogen yield from the SNII channel.
SNII_yield_factor_Helium: 1.0 # (Optional) Correction factor to apply to the Helium yield from the SNII channel.
SNII_yield_factor_Carbon: 0.5 # (Optional) Correction factor to apply to the Carbon yield from the SNII channel.
diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c
index 59612ca2f494f346548e4ec2543cd87e6d9c3941..01c96b83fb4cd78592569d6b38000cac21038ec7 100644
--- a/src/feedback/EAGLE/feedback.c
+++ b/src/feedback/EAGLE/feedback.c
@@ -840,6 +840,9 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props,
#ifdef SWIFT_DEBUG_CHECKS
if (age < 0.f) error("Negative age for a star.");
+
+ if (sp->count_since_last_enrichment != 0)
+ error("Computing feedback on a star that should not");
#endif
/* Convert dt and stellar age from internal units to Gyr. */
@@ -981,6 +984,8 @@ void feedback_props_init(struct feedback_props* fp,
const struct hydro_props* hydro_props,
const struct cosmology* cosmo) {
+ const double Gyr_in_cgs = 1.0e9 * 365.25 * 24. * 3600.;
+
/* Main operation modes ------------------------------------------------- */
fp->with_SNII_feedback =
@@ -1027,7 +1032,6 @@ void feedback_props_init(struct feedback_props* fp,
if (!fp->SNII_sampled_delay) {
/* Set the delay time before SNII occur */
- const double Gyr_in_cgs = 1.0e9 * 365.25 * 24. * 3600.;
fp->SNII_wind_delay =
parser_get_param_double(params, "EAGLEFeedback:SNII_wind_delay_Gyr") *
Gyr_in_cgs / units_cgs_conversion_factor(us, UNIT_CONV_TIME);
@@ -1148,6 +1152,28 @@ void feedback_props_init(struct feedback_props* fp,
fp->AGB_ejecta_specific_kinetic_energy =
0.5f * ejecta_velocity * ejecta_velocity;
+ /* Properties of the enrichment down-sampling ----------------------------- */
+
+ fp->stellar_evolution_age_cut =
+ parser_get_param_double(params,
+ "EAGLEFeedback:stellar_evolution_age_cut_Gyr") *
+ Gyr_in_cgs / units_cgs_conversion_factor(us, UNIT_CONV_TIME);
+
+ fp->stellar_evolution_sampling_rate = parser_get_param_double(
+ params, "EAGLEFeedback:stellar_evolution_sampling_rate");
+
+ if (fp->stellar_evolution_sampling_rate < 1 ||
+ fp->stellar_evolution_sampling_rate >= (1 << (8 * sizeof(char) - 1)))
+ error("Stellar evolution sampling rate too large. Must be >0 and <%d",
+ (1 << (8 * sizeof(char) - 1)));
+
+ /* Check that we are not downsampling before reaching the SNII delay */
+ if (!fp->SNII_sampled_delay &&
+ fp->stellar_evolution_age_cut < fp->SNII_wind_delay)
+ error(
+ "Time at which the enrichment downsampling stars is lower than the "
+ "SNII wind delay!");
+
/* Gather common conversion factors --------------------------------------- */
/* Calculate internal mass to solar mass conversion factor */
diff --git a/src/feedback/EAGLE/feedback.h b/src/feedback/EAGLE/feedback.h
index 2c89a76371b27cbcf0c5232f5bcf9ba6d7462379..c798eb386c564b8911593c08e4d820e73d5551e8 100644
--- a/src/feedback/EAGLE/feedback.h
+++ b/src/feedback/EAGLE/feedback.h
@@ -33,17 +33,6 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props,
const struct unit_system* us, const double age,
const double dt);
-/**
- * @brief Should we do feedback for this star?
- *
- * @param sp The star to consider.
- */
-__attribute__((always_inline)) INLINE static int feedback_do_feedback(
- const struct spart* sp) {
-
- return (sp->birth_time != -1.);
-}
-
/**
* @brief Should this particle be doing any feedback-related operation?
*
@@ -56,13 +45,7 @@ __attribute__((always_inline)) INLINE static int feedback_is_active(
const struct spart* sp, const float time, const struct cosmology* cosmo,
const int with_cosmology) {
- if (sp->birth_time == -1.) return 0;
-
- if (with_cosmology) {
- return ((float)cosmo->a) > sp->birth_scale_factor;
- } else {
- return time > sp->birth_time;
- }
+ return (sp->birth_time != -1.) && (sp->count_since_last_enrichment == 0);
}
/**
@@ -77,6 +60,30 @@ __attribute__((always_inline)) INLINE static void feedback_init_spart(
sp->feedback_data.to_collect.ngb_mass = 0.f;
}
+/**
+ * @brief Returns the length of time since the particle last did
+ * enrichment/feedback.
+ *
+ * @param sp The #spart.
+ * @param with_cosmology Are we running with cosmological time integration on?
+ * @param cosmo The cosmological model.
+ * @param time The current time (since the Big Bang / start of the run) in
+ * internal units.
+ * @param dt_star the length of this particle's time-step in internal units.
+ * @return The length of the enrichment step in internal units.
+ */
+INLINE static double feedback_get_enrichment_timestep(
+ const struct spart* sp, const int with_cosmology,
+ const struct cosmology* cosmo, const double time, const double dt_star) {
+
+ if (with_cosmology) {
+ return cosmology_get_delta_time_from_scale_factors(
+ cosmo, (double)sp->last_enrichment_time, cosmo->a);
+ } else {
+ return time - sp->last_enrichment_time;
+ }
+}
+
/**
* @brief Prepares a star's feedback field before computing what
* needs to be distributed.
@@ -153,11 +160,14 @@ __attribute__((always_inline)) INLINE static void feedback_prepare_spart(
* @param star_age_beg_step The age of the star at the star of the time-step in
* internal units.
* @param dt The time-step size of this star in internal units.
+ * @param time The physical time in internal units.
+ * @param with_cosmology Are we running with cosmology on?
*/
__attribute__((always_inline)) INLINE static void feedback_evolve_spart(
struct spart* restrict sp, const struct feedback_props* feedback_props,
const struct cosmology* cosmo, const struct unit_system* us,
- const double star_age_beg_step, const double dt) {
+ const double star_age_beg_step, const double dt, const double time,
+ const int with_cosmology) {
#ifdef SWIFT_DEBUG_CHECKS
if (sp->birth_time == -1.) error("Evolving a star particle that should not!");
@@ -170,6 +180,64 @@ __attribute__((always_inline)) INLINE static void feedback_evolve_spart(
/* Decrease star mass by amount of mass distributed to gas neighbours */
sp->mass -= sp->feedback_data.to_distribute.mass;
+
+ /* Mark this is the last time we did enrichment */
+ if (with_cosmology)
+ sp->last_enrichment_time = cosmo->a;
+ else
+ sp->last_enrichment_time = time;
+}
+
+/**
+ * @brief Will this star particle want to do feedback during the next time-step?
+ *
+ * @param sp The star of interest.
+ * @param feedback_props The properties of the feedback model.
+ * @param age_of_star Age of the star in internal units.
+ */
+__attribute__((always_inline)) INLINE static int feedback_will_do_feedback(
+ struct spart* restrict sp, const struct feedback_props* feedback_props,
+ const int with_cosmology, const struct cosmology* cosmo,
+ const double time) {
+
+ /* Calculate age of the star at current time */
+ double age_of_star;
+ if (with_cosmology) {
+ age_of_star = cosmology_get_delta_time_from_scale_factors(
+ cosmo, (double)sp->birth_scale_factor, cosmo->a);
+ } else {
+ age_of_star = time - (double)sp->birth_time;
+ }
+
+ /* Is the star still young? */
+ if (age_of_star < feedback_props->stellar_evolution_age_cut) {
+
+ /* Set the counter to "let's do enrichment" */
+ sp->count_since_last_enrichment = 0;
+
+ /* Say we want to do feedback */
+ return 1;
+
+ } else {
+
+ /* Increment counter */
+ sp->count_since_last_enrichment++;
+
+ if ((sp->count_since_last_enrichment %
+ feedback_props->stellar_evolution_sampling_rate) == 0) {
+
+ /* Reset counter */
+ sp->count_since_last_enrichment = 0;
+
+ /* Say we want to do feedback */
+ return 1;
+
+ } else {
+
+ /* Say we don't want to do feedback */
+ return 0;
+ }
+ }
}
void feedback_clean(struct feedback_props* feedback_props);
diff --git a/src/feedback/EAGLE/feedback_iact.h b/src/feedback/EAGLE/feedback_iact.h
index 244f9718d5c62c6b19258fd0b6a078f97732195d..832972aa5adeeda70b5f3015778d00c9cf6b6fe3 100644
--- a/src/feedback/EAGLE/feedback_iact.h
+++ b/src/feedback/EAGLE/feedback_iact.h
@@ -91,6 +91,11 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
const struct cosmology *cosmo,
const integertime_t ti_current) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (si->count_since_last_enrichment != 0)
+ error("Computing feedback from a star that should not");
+#endif
+
/* Get r. */
const float r = sqrtf(r2);
@@ -113,8 +118,10 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
#ifdef SWIFT_DEBUG_CHECKS
if (Omega_frac < 0. || Omega_frac > 1.01)
- error("Invalid fraction of material to distribute. Omega_frac=%e",
- Omega_frac);
+ error(
+ "Invalid fraction of material to distribute for star ID=%lld "
+ "Omega_frac=%e count since last enrich=%d",
+ si->id, Omega_frac, si->count_since_last_enrichment);
#endif
/* Update particle mass */
diff --git a/src/feedback/EAGLE/feedback_properties.h b/src/feedback/EAGLE/feedback_properties.h
index a0107ccdd7f499e750beba34bce6e195312a8d5e..55822742587e9e7f40e25f4d8c2d3bd9c91f0d33 100644
--- a/src/feedback/EAGLE/feedback_properties.h
+++ b/src/feedback/EAGLE/feedback_properties.h
@@ -264,6 +264,15 @@ struct feedback_props {
/*! Slope of the metallicity dependance of the feedback energy fraction model
*/
double n_Z;
+
+ /* ------------ Enrichment sampling properties ------------ */
+
+ /*! Star age above which the enrichment will be downsampled (in internal
+ * units) */
+ double stellar_evolution_age_cut;
+
+ /*! Number of time-steps in-between two enrichment events */
+ int stellar_evolution_sampling_rate;
};
void feedback_props_init(struct feedback_props *fp,
diff --git a/src/feedback/none/feedback.h b/src/feedback/none/feedback.h
index 98cf0aa708320046e2586f203b683069db5f8482..85a8cf9566af25b179490df66d7e9fbeef0ee461 100644
--- a/src/feedback/none/feedback.h
+++ b/src/feedback/none/feedback.h
@@ -64,6 +64,29 @@ __attribute__((always_inline)) INLINE static int feedback_is_active(
return 1;
}
+/**
+ * @brief Returns the length of time since the particle last did
+ * enrichment/feedback.
+ *
+ * We just return the normal time-step here since particles do something every
+ * regular time-step.
+ *
+ * @param sp The #spart.
+ * @param with_cosmology Are we running with cosmological time integration on?
+ * @param cosmo The cosmological model.
+ * @param time The current time (since the Big Bang / start of the run) in
+ * internal units.
+ * @param dt_star the length of this particle's time-step in internal units.
+ * @return The length of the enrichment step in internal units.
+ */
+INLINE static double feedback_get_enrichment_timestep(
+ const struct spart* sp, const int with_cosmology,
+ const struct cosmology* cosmo, const double time, const double dt_star) {
+
+ /* Just return the regular step length */
+ return dt_star;
+}
+
/**
* @brief Prepares a star's feedback field before computing what
* needs to be distributed.
@@ -108,11 +131,30 @@ __attribute__((always_inline)) INLINE static void feedback_prepare_spart(
* @param star_age_beg_step The age of the star at the star of the time-step in
* internal units.
* @param dt The time-step size of this star in internal units.
+ * @param time The physical time in internal units.
+ * @param with_cosmology Are we running with cosmology on?
*/
__attribute__((always_inline)) INLINE static void feedback_evolve_spart(
struct spart* restrict sp, const struct feedback_props* feedback_props,
const struct cosmology* cosmo, const struct unit_system* us,
- const double star_age_beg_step, const double dt) {}
+ const double star_age_beg_step, const double dt, const double time,
+ const int with_cosmology) {}
+
+/**
+ * @brief Will this star particle want to do feedback during the next time-step?
+ *
+ * @param sp The star of interest.
+ * @param feedback_props The properties of the feedback model.
+ * @param with_cosmology Are we running with cosmological time integration?
+ * @param cosmo The #cosmology object.
+ * @param time The current time (since the Big Bang).
+ */
+__attribute__((always_inline)) INLINE static int feedback_will_do_feedback(
+ struct spart* restrict sp, const struct feedback_props* feedback_props,
+ const int with_cosmology, const struct cosmology* cosmo,
+ const double time) {
+ return 1;
+}
/**
* @brief Clean-up the memory allocated for the feedback routines
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index 50bd95e05faab51ce3241a0b299e99faef242d58..a4113b5b897e53d1906ff869e74aca64b00863e3 100644
--- a/src/runner_ghost.c
+++ b/src/runner_ghost.c
@@ -204,19 +204,19 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
stars_reset_feedback(sp);
/* Only do feedback if stars have a reasonable birth time */
- if (feedback_do_feedback(sp)) {
+ if (feedback_is_active(sp, e->time, cosmo, with_cosmology)) {
const integertime_t ti_step = get_integer_timestep(sp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(e->ti_current - 1, sp->time_bin);
/* Get particle time-step */
- double dt;
+ double dt_star;
if (with_cosmology) {
- dt = cosmology_get_delta_time(e->cosmology, ti_begin,
- ti_begin + ti_step);
+ dt_star = cosmology_get_delta_time(e->cosmology, ti_begin,
+ ti_begin + ti_step);
} else {
- dt = get_timestep(sp->time_bin, e->time_base);
+ dt_star = get_timestep(sp->time_bin, e->time_base);
}
/* Calculate age of the star at current time */
@@ -226,19 +226,22 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
cosmology_get_delta_time_from_scale_factors(
cosmo, (double)sp->birth_scale_factor, cosmo->a);
} else {
- star_age_end_of_step = (float)e->time - sp->birth_time;
+ star_age_end_of_step = e->time - (double)sp->birth_time;
}
/* Has this star been around for a while ? */
if (star_age_end_of_step > 0.) {
- /* Age of the star at the start of the step */
+ /* Get the length of the enrichment time-step */
+ const double dt_enrichment = feedback_get_enrichment_timestep(
+ sp, with_cosmology, cosmo, e->time, dt_star);
const double star_age_beg_of_step =
- max(star_age_end_of_step - dt, 0.);
+ star_age_end_of_step - dt_enrichment;
/* Compute the stellar evolution */
feedback_evolve_spart(sp, feedback_props, cosmo, us,
- star_age_beg_of_step, dt);
+ star_age_beg_of_step, dt_enrichment,
+ e->time, with_cosmology);
} else {
/* Reset the feedback fields of the star particle */
@@ -345,40 +348,43 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
stars_reset_feedback(sp);
/* Only do feedback if stars have a reasonable birth time */
- if (feedback_do_feedback(sp)) {
+ if (feedback_is_active(sp, e->time, cosmo, with_cosmology)) {
const integertime_t ti_step = get_integer_timestep(sp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(e->ti_current - 1, sp->time_bin);
/* Get particle time-step */
- double dt;
+ double dt_star;
if (with_cosmology) {
- dt = cosmology_get_delta_time(e->cosmology, ti_begin,
- ti_begin + ti_step);
+ dt_star = cosmology_get_delta_time(e->cosmology, ti_begin,
+ ti_begin + ti_step);
} else {
- dt = get_timestep(sp->time_bin, e->time_base);
+ dt_star = get_timestep(sp->time_bin, e->time_base);
}
/* Calculate age of the star at current time */
double star_age_end_of_step;
if (with_cosmology) {
star_age_end_of_step = cosmology_get_delta_time_from_scale_factors(
- cosmo, sp->birth_scale_factor, (float)cosmo->a);
+ cosmo, (double)sp->birth_scale_factor, cosmo->a);
} else {
- star_age_end_of_step = (float)e->time - sp->birth_time;
+ star_age_end_of_step = e->time - (double)sp->birth_time;
}
/* Has this star been around for a while ? */
if (star_age_end_of_step > 0.) {
- /* Age of the star at the start of the step */
+ /* Get the length of the enrichment time-step */
+ const double dt_enrichment = feedback_get_enrichment_timestep(
+ sp, with_cosmology, cosmo, e->time, dt_star);
const double star_age_beg_of_step =
- max(star_age_end_of_step - dt, 0.);
+ star_age_end_of_step - dt_enrichment;
/* Compute the stellar evolution */
feedback_evolve_spart(sp, feedback_props, cosmo, us,
- star_age_beg_of_step, dt);
+ star_age_beg_of_step, dt_enrichment, e->time,
+ with_cosmology);
} else {
/* Reset the feedback fields of the star particle */
diff --git a/src/runner_time_integration.c b/src/runner_time_integration.c
index 09d0544f1e18e15bfdaf0083cd257b8c9ab3e91c..2a4b138b6c6da7a93a92b8097d2a3dc59e582370 100644
--- a/src/runner_time_integration.c
+++ b/src/runner_time_integration.c
@@ -30,6 +30,7 @@
#include "black_holes.h"
#include "cell.h"
#include "engine.h"
+#include "feedback.h"
#include "kick.h"
#include "timers.h"
#include "timestep.h"
@@ -746,6 +747,9 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
sp->time_bin = get_time_bin(ti_new_step);
sp->gpart->time_bin = get_time_bin(ti_new_step);
+ feedback_will_do_feedback(sp, e->feedback_props, with_cosmology,
+ e->cosmology, e->time);
+
/* Number of updated s-particles */
s_updated++;
g_updated++;
diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h
index 45b3a2ce5436395a6f1576c7e0d89efcf8bf15d3..767fa2715ab35f21ae41c5162787c3a544b87da7 100644
--- a/src/star_formation/EAGLE/star_formation.h
+++ b/src/star_formation/EAGLE/star_formation.h
@@ -447,6 +447,8 @@ INLINE static void star_formation_copy_properties(
/* Flag that this particle has not done feedback yet */
sp->f_E = -1.f;
+ sp->last_enrichment_time = sp->birth_time;
+ sp->count_since_last_enrichment = -1;
}
/**
diff --git a/src/stars/EAGLE/stars.h b/src/stars/EAGLE/stars.h
index 82ad2e25f474e6c66e6cf5e09f39188faa09b084..e9ea6ed0b8d2398f2a75b37430a295f62b3b22e8 100644
--- a/src/stars/EAGLE/stars.h
+++ b/src/stars/EAGLE/stars.h
@@ -68,6 +68,9 @@ __attribute__((always_inline)) INLINE static void stars_first_init_spart(
if (stars_properties->overwrite_birth_time)
sp->birth_time = stars_properties->spart_first_init_birth_time;
+ sp->last_enrichment_time = sp->birth_time;
+ sp->count_since_last_enrichment = -1;
+
stars_init_spart(sp);
}
diff --git a/src/stars/EAGLE/stars_part.h b/src/stars/EAGLE/stars_part.h
index 94b9315ba77cca700a4bfe5109db2873e678829f..13f0783674534514213f6bac406f4b41b9dfaaf1 100644
--- a/src/stars/EAGLE/stars_part.h
+++ b/src/stars/EAGLE/stars_part.h
@@ -56,9 +56,6 @@ struct spart {
/*! Star mass */
float mass;
- /*! Initial star mass */
- float mass_init;
-
/*! Particle smoothing length. */
float h;
@@ -82,6 +79,12 @@ struct spart {
float birth_scale_factor;
};
+ /*! Scale-factor / time at which this particle last did enrichment */
+ float last_enrichment_time;
+
+ /*! Initial star mass */
+ float mass_init;
+
/*! Birth density */
float birth_density;
@@ -103,6 +106,9 @@ struct spart {
/*! Particle time bin */
timebin_t time_bin;
+ /*! Number of time-steps since the last enrichment step */
+ char count_since_last_enrichment;
+
#ifdef WITH_LOGGER
/* Additional data for the particle logger */
struct logger_part_data logger_data;
@@ -149,7 +155,7 @@ struct stars_props {
/*! Smoothing length tolerance */
float h_tolerance;
- /*! Tolerance on neighbour number (for info only) */
+ /*! Tolerance on neighbour number (for info only)*/
float delta_neighbours;
/*! Maximal number of iterations to converge h */