Skip to content
GitLab
Commit 1aee343e authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Sampled stellar evolution

parent bde56afa
Hide whitespace changes
Inline Side-by-side
.gitignore
View file @ 1aee343e
......@@ -32,6 +32,7 @@ examples/*/*/*.dat
examples/*/*/*.png
examples/*/*/*.pdf
examples/*/*/*.mp4
examples/*/*/*.txt
examples/*/*/*.rst
examples/*/*/*.hdf5
examples/*/*/*.csv
......
doc/RTD/source/SubgridModels/EAGLE/index.rst
View file @ 1aee343e
......@@ -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.
......
examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
View file @ 1aee343e
......@@ -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.
......
examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
View file @ 1aee343e
......@@ -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.
......
examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
View file @ 1aee343e
......@@ -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.
......
examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
View file @ 1aee343e
......@@ -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.
......
examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
View file @ 1aee343e
......@@ -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.
......
examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
View file @ 1aee343e
......@@ -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.
......
examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
View file @ 1aee343e
......@@ -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.
......
examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
View file @ 1aee343e
......@@ -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.
......
examples/SubgridTests/StellarEvolution/plot_box_evolution.py
View file @ 1aee343e
......@@ -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")
......
examples/SubgridTests/StellarEvolution/stellar_evolution.yml
View file @ 1aee343e
......@@ -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.
......
examples/parameter_example.yml
View file @ 1aee343e
......@@ -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.
......
src/feedback/EAGLE/feedback.c
View file @ 1aee343e
......@@ -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 */
......
src/feedback/EAGLE/feedback.h
View file @ 1aee343e
......@@ -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);
......
src/feedback/EAGLE/feedback_iact.h
View file @ 1aee343e
......@@ -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 */
......
src/feedback/EAGLE/feedback_properties.h
View file @ 1aee343e
......@@ -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,
......
src/feedback/none/feedback.h
View file @ 1aee343e
......@@ -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
......
src/runner_ghost.c
View file @ 1aee343e
......@@ -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 {