diff --git a/examples/StellarEvolution/stellar_evolution.yml b/examples/StellarEvolution/stellar_evolution.yml
index f66f7c6cff46e00d631434831d5a9d27c7c73776..e909843047ad958d51289ba51c626711c6bb93a0 100644
--- a/examples/StellarEvolution/stellar_evolution.yml
+++ b/examples/StellarEvolution/stellar_evolution.yml
@@ -47,11 +47,6 @@ Scheduler:
max_top_level_cells: 4
cell_split_size: 50
-Stars:
- feedback_timescale: 1e-4
- continuous_heating: 0 # decide whether using continuous heating in feedback (1) or stochastic (0)
- energy_testing: 1 # Are we testing energy injection? (1: yes, 0: no)
-
Gravity:
mesh_side_length: 32
eta: 0.025
@@ -59,3 +54,10 @@ Gravity:
r_cut_max: 5.
comoving_softening: 0.001
max_physical_softening: 0.001
+
+Stars:
+ feedback_timescale: 1.e-4
+
+EagleStellarEvolution:
+ filename: /cosma5/data/Eagle/BG_Tables/YieldTables/
+ imf_model: Chabrier
diff --git a/examples/main.c b/examples/main.c
index 1d9f3baa0bcc6be64340f8dae25a72d3ed0f8f23..2bfb7ea7600c119e067dc971fc494b74be0d9c44 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -912,6 +912,11 @@ int main(int argc, char *argv[]) {
bzero(&chemistry, sizeof(struct chemistry_global_data));
chemistry_init(params, &us, &prog_const, &chemistry);
if (myrank == 0) chemistry_print(&chemistry);
+
+ /* Initialise the cooling function properties */
+ if (with_feedback)
+ stars_evolve_init(params, &stars_properties);
+ // ALEXEI: add print function!!!
/* Construct the engine policy */
int engine_policies = ENGINE_POLICY | engine_policy_steal;
diff --git a/src/drift.h b/src/drift.h
index aa98e91749a31cf4ecdb9a8688e2e767577c3217..7e874fe0ceabe5b091cc7c5bb53adbef2c9a3efd 100644
--- a/src/drift.h
+++ b/src/drift.h
@@ -92,7 +92,6 @@ __attribute__((always_inline)) INLINE static void drift_part(
p->x[2] += xp->v_full[2] * dt_drift;
/* Predict velocities (for hydro terms) */
- if (p->id == 145267) message("id %llu v_old %.5e %.5e %.5e a_hydro %.5e %.5e %.5e dt_kick_hydro %.5e v_new %.5e %.5e %.5e", p->id, p->v[0], p->v[1], p->v[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2], dt_kick_hydro, p->v[0] + p->a_hydro[0] * dt_kick_hydro, p->v[1] + p->a_hydro[1] * dt_kick_hydro, p->v[2] + p->a_hydro[2] * dt_kick_hydro);
p->v[0] += p->a_hydro[0] * dt_kick_hydro;
p->v[1] += p->a_hydro[1] * dt_kick_hydro;
p->v[2] += p->a_hydro[2] * dt_kick_hydro;
diff --git a/src/runner.c b/src/runner.c
index 65a6705966e5b1a6d4396987334845a077223fd6..a8c431a6d62c10270f31bf27792d40a0931a95d8 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -133,6 +133,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
struct spart *restrict sparts = c->stars.parts;
const struct engine *e = r->e;
+ const struct unit_system *us = e->internal_units;
const int with_cosmology = (e->policy & engine_policy_cosmology);
const struct cosmology *cosmo = e->cosmology;
const float stars_h_max = e->hydro_properties->h_max;
@@ -351,8 +352,11 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
/* We now have a particle whose smoothing length has converged */
stars_reset_feedback(sp);
+ // Get current time
+ float current_time_begin = get_integer_time_begin(e->ti_current - 1, sp->time_bin) * e->time_base;
+
/* Compute the stellar evolution */
- stars_evolve_spart(sp, e->stars_properties, cosmo, dt);
+ stars_evolve_spart(sp, e->stars_properties, cosmo, us, current_time_begin, dt);
}
/* We now need to treat the particles whose smoothing length had not
diff --git a/src/stars/EAGLE/imf.h b/src/stars/EAGLE/imf.h
index f1e257d560d1c637762a42a8b0efee52ed8a6ba5..cc2a6a3377629a47d164aa41346aa3e12d19c11e 100644
--- a/src/stars/EAGLE/imf.h
+++ b/src/stars/EAGLE/imf.h
@@ -213,7 +213,7 @@ inline static void init_imf(struct stars_props *restrict star_properties){
star_properties->imf_mass_bin_log10[i] = log_solar_mass;
}
} else {
- error("Invalid IMF model. Valid models are: PowerLaw and Chabrier\n");
+ error("Invalid IMF model %s. Valid models are: PowerLaw and Chabrier\n", star_properties->IMF_Model);
}
norm = integrate_imf(log_imf_min_solar_mass, log_imf_max_solar_mass,
diff --git a/src/stars/EAGLE/stars.h b/src/stars/EAGLE/stars.h
index 94910f28ea343b4eeae5c9800ad6b1058b1219e0..d99efa1df7e9f11f04a1f7c5181d472f45541bd1 100644
--- a/src/stars/EAGLE/stars.h
+++ b/src/stars/EAGLE/stars.h
@@ -441,8 +441,9 @@ inline static void determine_bin_yield(int *iz_low, int *iz_high, float *dz, flo
inline static void evolve_SNIa(float log10_min_mass, float log10_max_mass,
const struct stars_props *restrict stars,
- struct spart *restrict sp, float dt_Gyr){
- int i;
+ struct spart *restrict sp,
+ const struct unit_system* us,
+ float star_age_Gyr, float dt_Gyr){
/* Check if we're outside the mass range for SNIa */
if (log10_min_mass >= log10_SNIa_max_mass_msun) return;
@@ -451,56 +452,49 @@ inline static void evolve_SNIa(float log10_min_mass, float log10_max_mass,
if (log10_max_mass > log10_SNIa_max_mass_msun) {
log10_max_mass = log10_SNIa_max_mass_msun;
float lifetime_Gyr = lifetime_in_Gyr(exp(M_LN10 * log10_SNIa_max_mass_msun), sp->chemistry_data.metal_mass_fraction_total, stars);
- dt_Gyr = sp->time_since_enrich_Gyr + dt_Gyr - lifetime_Gyr;
- sp->time_since_enrich_Gyr = lifetime_Gyr;
+ dt_Gyr = star_age_Gyr + dt_Gyr - lifetime_Gyr;
+ star_age_Gyr = lifetime_Gyr;
}
- /* compute the fraction of white dwarfs */
- float num_of_SNIa_per_msun;
+ /* compute the number of SNIa per solar mass */
/* Efolding (Forster 2006) */
- num_of_SNIa_per_msun =
+ sp->to_distribute.num_SNIa =
stars->SNIa_efficiency *
- (exp(-sp->time_since_enrich_Gyr / stars->SNIa_timescale) -
- exp(-(sp->time_since_enrich_Gyr + dt_Gyr) / stars->SNIa_timescale));
-
- // switch included in EAGLE but there don't seem to be any alternatives.
- //switch (stars->SNIa_mode) {
- // case 2:
- // /* Efolding (Forster 2006) */
- // num_of_SNIa_per_msun =
- // stars->SNIa_efficiency *
- // (exp(-sp->time_since_enrich_Gyr / stars->SNIa_timescale) -
- // exp(-(sp->time_since_enrich_Gyr + dt_Gyr) / stars->SNIa_timescale));
- // break;
- // default:
- // error("SNIa mode not defined yet %d\n", stars->SNIa_mode);
- // break;
- //}
+ (exp(-star_age_Gyr / stars->SNIa_timescale) -
+ exp(-(star_age_Gyr + dt_Gyr) / stars->SNIa_timescale));
- sp->num_snia = num_of_SNIa_per_msun;
+ /* compute total mass released by SNIa */
+ sp->to_distribute.mass += sp->to_distribute.num_SNIa * stars->yield_SNIa_total_metals_SPH; // / units_cgs_conversion_factor(us, UNIT_CONV_MASS);
+ message("num_snia %.5e mass %.5e", sp->to_distribute.num_SNIa, sp->to_distribute.mass);
- if (stars->SNIa_mass_transfer == 1) {
- for (i = 0; i < chemistry_element_count; i++) {
- sp->metals_released[i] += num_of_SNIa_per_msun * stars->yield_SNIa_SPH[i];
- }
+ /* compute mass fractions of each metal */
+ for (int i = 0; i < chemistry_element_count; i++) {
+ sp->to_distribute.chemistry_data.metal_mass_fraction[i] += sp->to_distribute.num_SNIa * stars->yield_SNIa_SPH[i] / sp->to_distribute.mass;
+ }
- sp->chemistry_data.mass_from_SNIa += num_of_SNIa_per_msun * stars->yield_SNIa_total_metals_SPH;
- sp->chemistry_data.metal_mass_fraction_from_SNIa += num_of_SNIa_per_msun * stars->yield_SNIa_total_metals_SPH;
+ // For diagnostics according to Richard
+ //if (stars->SNIa_mass_transfer == 1) {
+ // for (i = 0; i < chemistry_element_count; i++) {
+ // sp->metals_released[i] += num_of_SNIa_per_msun * stars->yield_SNIa_SPH[i];
+ // }
- sp->metal_mass_released += num_of_SNIa_per_msun * stars->yield_SNIa_total_metals_SPH;
+ // sp->chemistry_data.mass_from_SNIa += num_of_SNIa_per_msun * stars->yield_SNIa_total_metals_SPH;
+ // sp->chemistry_data.metal_mass_fraction_from_SNIa += num_of_SNIa_per_msun * stars->yield_SNIa_total_metals_SPH;
- // Make sure chemistry_element_Fe corresponds to the iron_index used in EAGLE!!!
- sp->chemistry_data.iron_mass_fraction_from_SNIa += num_of_SNIa_per_msun * stars->yield_SNIa_SPH[chemistry_element_Fe];
+ // sp->metal_mass_released += num_of_SNIa_per_msun * stars->yield_SNIa_total_metals_SPH;
- /* metal_mass_released is the yield of ALL metals, not just the
- 11 tabulated in the code. SNIa remnants inject no H or He
- so chemistry_data.mass_from_SNIa == chemistry_data.metal_mass_fraction_from_SNIa */
+ // // Make sure chemistry_element_Fe corresponds to the iron_index used in EAGLE!!!
+ // sp->chemistry_data.iron_mass_fraction_from_SNIa += num_of_SNIa_per_msun * stars->yield_SNIa_SPH[chemistry_element_Fe];
- } else {
- sp->chemistry_data.iron_mass_fraction_from_SNIa = 0;
- sp->chemistry_data.metal_mass_fraction_from_SNIa = 0;
- sp->chemistry_data.mass_from_SNIa = 0;
- }
+ // /* metal_mass_released is the yield of ALL metals, not just the
+ // 11 tabulated in the code. SNIa remnants inject no H or He
+ // so chemistry_data.mass_from_SNIa == chemistry_data.metal_mass_fraction_from_SNIa */
+
+ //} else {
+ // sp->chemistry_data.iron_mass_fraction_from_SNIa = 0;
+ // sp->chemistry_data.metal_mass_fraction_from_SNIa = 0;
+ // sp->chemistry_data.mass_from_SNIa = 0;
+ //}
}
inline static void evolve_SNII(float log10_min_mass, float log10_max_mass,
@@ -743,27 +737,30 @@ inline static void evolve_AGB(float log10_min_mass, float log10_max_mass,
// Analogue of eagle_do_stellar_evolution...
inline static void compute_stellar_evolution(const struct stars_props *restrict star_properties,
- struct spart *restrict sp, float dt) {
+ struct spart *restrict sp,
+ const struct unit_system* us,
+ float age, float dt) {
// Convert dt from internal units to Gyr. (Do correct conversion!)
- float convert_time_to_Gyr = 1.0;
- float dt_Gyr = dt*convert_time_to_Gyr;
- float age_of_star_Gyr = sp->birth_time*convert_time_to_Gyr; // This should be same as age_of_star_in_Gyr_begstep in EAGLE, check. Also, make sure it works with cosmology, might need to use birth_scale_factor.
+ const float Gyr_in_cgs = 3.154e16;
+ float dt_Gyr = dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME) / Gyr_in_cgs;
+ float star_age_Gyr = age * units_cgs_conversion_factor(us, UNIT_CONV_TIME) / Gyr_in_cgs; // This should be same as age_of_star_in_Gyr_begstep in EAGLE, check. Also, make sure it works with cosmology, might need to use birth_scale_factor.
// set max and min mass of dying stars
float log10_max_dying_mass =
- log10(dying_mass_msun(age_of_star_Gyr, sp->chemistry_data.metal_mass_fraction_total, star_properties));
+ log10(dying_mass_msun(star_age_Gyr, sp->chemistry_data.metal_mass_fraction_total, star_properties));
float log10_min_dying_mass = log10(
- dying_mass_msun(age_of_star_Gyr + dt_Gyr, sp->chemistry_data.metal_mass_fraction_total, star_properties));
+ dying_mass_msun(star_age_Gyr + dt_Gyr, sp->chemistry_data.metal_mass_fraction_total, star_properties));
if (log10_min_dying_mass > log10_max_dying_mass) error("min dying mass is greater than max dying mass");
+ message("age %.5e age of star %.5e log10 min max dying mass %.5e %.5e", age, star_age_Gyr, log10_min_dying_mass, log10_max_dying_mass);
/* integration interval is zero - this can happen if minimum and maximum
* dying masses are above imf_max_mass_msun */
if (log10_min_dying_mass == log10_max_dying_mass) return;
/* Evolve SNIa, SNII, AGB */
- evolve_SNIa(log10_min_dying_mass,log10_max_dying_mass,star_properties,sp,dt_Gyr);
+ evolve_SNIa(log10_min_dying_mass,log10_max_dying_mass,star_properties,sp,us,star_age_Gyr,dt_Gyr);
evolve_SNII(log10_min_dying_mass,log10_max_dying_mass,star_properties,sp);
evolve_AGB(log10_min_dying_mass,log10_max_dying_mass,star_properties,sp);
@@ -781,9 +778,14 @@ inline static void compute_stellar_evolution(const struct stars_props *restrict
*/
__attribute__((always_inline)) INLINE static void stars_evolve_spart(
struct spart* restrict sp, const struct stars_props* stars_properties,
- const struct cosmology* cosmo, double dt) {
+ const struct cosmology* cosmo, const struct unit_system* us, float current_time, double dt) {
+
+ // Set birth time for testing purposes
+ sp->birth_time = 0;
+ float star_age = current_time - sp->birth_time;
- sp->num_snia = 0;
+ sp->to_distribute.num_SNIa = 0;
+ sp->to_distribute.mass = 0;
// get mass and initial mass of particle to pass to evolution function
// except we're working in mass fraction so maybe not necessary...
//float mass = hydro_get_mass(p);
@@ -801,9 +803,8 @@ __attribute__((always_inline)) INLINE static void stars_evolve_spart(
sp->chemistry_data.iron_mass_fraction_from_SNIa = 0;
// Evolve the star
- compute_stellar_evolution(stars_properties, sp, dt);
-
- // set_particle_metal_content
+ compute_stellar_evolution(stars_properties, sp, us, star_age, dt);
+
}
inline static void stars_evolve_init(struct swift_params *params, struct stars_props* restrict stars){
@@ -825,11 +826,9 @@ inline static void stars_evolve_init(struct swift_params *params, struct stars_p
stars->AGB_mass_transfer = 1;
stars->SNII_mass_transfer = 1;
- /* Which stellar lifetime model are we using? (To do: read from yml file) */
- stars->stellar_lifetime_flag = 2;
-
/* Yield table filepath */
parser_get_param_string(params, "EagleStellarEvolution:filename", stars->yield_table_path);
+ parser_get_param_string(params, "EagleStellarEvolution:imf_model", stars->IMF_Model);
//stars->yield_SNIa_total_metals_SPH = ;
@@ -848,6 +847,8 @@ inline static void stars_evolve_init(struct swift_params *params, struct stars_p
/* Further calculation on tables to convert them to log10 and compute yields for each element */
compute_yields(stars);
+ message("initialized stellar feedback");
+
}
/**
diff --git a/src/stars/EAGLE/stars_iact.h b/src/stars/EAGLE/stars_iact.h
index 0c009bc9cdca291651e68e6c4f61d442876cec09..bfea2a2c05af669ebd41ad48e3b515b31a1a6ff0 100644
--- a/src/stars/EAGLE/stars_iact.h
+++ b/src/stars/EAGLE/stars_iact.h
@@ -1,3 +1,5 @@
+#include "random.h"
+
/**
* @brief Density interaction between two particles (non-symmetric).
*
@@ -120,8 +122,9 @@ runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
/* Update particle mass */
const float current_mass = hydro_get_mass(pj);
float new_mass = current_mass + si->to_distribute.mass*omega_frac;
- if (stars_properties->const_feedback_energy_testing) new_mass = current_mass;
hydro_set_mass(pj,new_mass);
+
+ message("particle %llu new mass %.5e old mass %.5e star distributed mass %.5e omega_frac %.5e", pj->id, new_mass, current_mass, si->to_distribute.mass, omega_frac);
/* Decrease the mass of star particle */
si->mass -= si->to_distribute.mass*omega_frac;
@@ -185,31 +188,24 @@ runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
d_energy = si->to_distribute.mass * (si->to_distribute.ejecta_specific_thermal_energy
+ 0.5*(si->v[0]*si->v[0] + si->v[1]*si->v[1] + si->v[2]*si->v[2]) * cosmo->a2_inv);
- // If statement temporary for testing, in practice would always be on.
- if (stars_properties->const_feedback_energy_testing) {
- if (stars_properties->continuous_heating) {
- // We're doing ONLY continuous heating
- d_energy += si->to_distribute.num_SNIa * stars_properties->total_energy_SNe * omega_frac * si->mass_init;
- du = d_energy/hydro_get_mass(pj);
- thermal_feedback(du,pj,xp,cosmo);
- } else {
- // We're doing stochastic heating
- heating_probability = stars_properties->SNe_temperature_h * si->to_distribute.num_SNIa *
- stars_properties->SNIa_energy_fraction /
- (stars_properties->SNe_deltaT_desired * si->ngb_mass);
- du = stars_properties->SNe_deltaT_desired * stars_properties->temp_to_u_factor;
- if (heating_probability >= 1) {
- du = stars_properties->SNe_energy_h * si->to_distribute.num_SNIa / si->ngb_mass;
- heating_probability = 1;
- }
+ if (stars_properties->continuous_heating) {
+ // We're doing ONLY continuous heating
+ d_energy += si->to_distribute.num_SNIa * stars_properties->total_energy_SNe * omega_frac * si->mass_init;
+ du = d_energy/hydro_get_mass(pj);
+ thermal_feedback(du,pj,xp,cosmo);
+ } else {
+ // We're doing stochastic heating
+ heating_probability = stars_properties->SNe_temperature_h * si->to_distribute.num_SNIa *
+ stars_properties->SNIa_energy_fraction /
+ (stars_properties->SNe_deltaT_desired * si->ngb_mass);
+ du = stars_properties->SNe_deltaT_desired * stars_properties->temp_to_u_factor;
+ if (heating_probability >= 1) {
+ du = stars_properties->SNe_energy_h * si->to_distribute.num_SNIa / si->ngb_mass;
+ heating_probability = 1;
}
}
- /* pick random number to see if we do stochastic heating */
- // Temporary assignment of random seed. Discuss with Matthieu for better
- // way of generating random numbers
- unsigned int seed = 3*(pj->id + ti_current) % 8191;
- double random_num = rand_r(&seed) * stars_properties->inv_rand_max;
+ double random_num = random_unit_interval(pj->id, ti_current, random_number_stellar_feedback);
if (random_num < heating_probability) {
message("we did some heating! id %llu probability %.5e random_num %.5e du %.5e du/ini %.5e", pj->id, heating_probability, random_num, du, du/hydro_get_physical_internal_energy(pj,xp,cosmo));
thermal_feedback(du, pj, xp, cosmo);
diff --git a/src/stars/EAGLE/stars_io.h b/src/stars/EAGLE/stars_io.h
index b35152d4a1e2bebb5eca620fba0e5a5546ce4b23..0b586e3dbf15b39d4d7284550f7d6a3fc30e95d3 100644
--- a/src/stars/EAGLE/stars_io.h
+++ b/src/stars/EAGLE/stars_io.h
@@ -134,7 +134,12 @@ INLINE static void stars_props_init(struct stars_props *sp,
sp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
sp->stellar_lifetime_flag = parser_get_opt_param_int(
- params, "Stars: lifetime_flag", 0);
+ params, "EAGLEFeedback:lifetime_flag", 0);
+
+ sp->SNIa_timescale = parser_get_opt_param_float(
+ params, "EAGLEFeedback:SNIa_timescale", 2.f);
+ sp->SNIa_efficiency = parser_get_opt_param_float(
+ params, "EAGLEFeedback:SNIa_efficiency", 2.e-3);
}
/**
diff --git a/src/stars/EAGLE/yield_tables.h b/src/stars/EAGLE/yield_tables.h
index 95f63b05aa3904f5b760e2971a893b3fe1118fb8..03c3ece6c4b235701a3b54deb63daef15b3fbb0e 100644
--- a/src/stars/EAGLE/yield_tables.h
+++ b/src/stars/EAGLE/yield_tables.h
@@ -351,7 +351,7 @@ inline static void read_yield_tables(struct stars_props *restrict stars){
}
inline static void allocate_yield_tables(struct stars_props *restrict stars){
-
+
/* Allocate SNIa arrays */
if (posix_memalign((void **)&stars->yields_SNIa, SWIFT_STRUCT_ALIGNMENT, stars->SNIa_n_elements * sizeof(double)) !=0) {
error("Failed to allocate SNIa yields array");
diff --git a/src/stars/const/stars.h b/src/stars/const/stars.h
index 50108b6d4e04c85958587aa406b24610acf42ec8..aec793e0053fa64d768516629f206ef777e175ab 100644
--- a/src/stars/const/stars.h
+++ b/src/stars/const/stars.h
@@ -184,7 +184,7 @@ __attribute__((always_inline)) INLINE static void stars_reset_acceleration(
*/
__attribute__((always_inline)) INLINE static void stars_evolve_spart(
struct spart* restrict sp, const struct stars_props* stars_properties,
- const struct cosmology* cosmo, double dt) {
+ const struct cosmology* cosmo, const struct unit_system* us, float current_time, double dt) {
/* Proportion of quantities to be released each timestep */
float feedback_factor = dt/stars_properties->feedback_timescale;
@@ -221,4 +221,29 @@ __attribute__((always_inline)) INLINE static void stars_evolve_spart(
}
+inline static void stars_evolve_init(struct swift_params *params, struct stars_props* restrict stars){}
+
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * This is the equivalent of hydro_reset_acceleration.
+ * We do not compute the acceleration on star, therefore no need to use it.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void stars_reset_feedback(
+ struct spart* restrict p) {
+
+ /* Reset time derivative */
+ p->feedback.h_dt = 0.f;
+
+#ifdef DEBUG_INTERACTIONS_STARS
+ for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
+ p->ids_ngbs_force[i] = -1;
+ p->num_ngb_force = 0;
+#endif
+}
+
+
#endif /* SWIFT_CONST_STARS_H */