diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/run.sh b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/run.sh
index afe0b60f41b5c2b82a19651a6d972ac05dba9c59..b5f01a8d5a0ef8d6536747189a6279adebcf9ef4 100755
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/run.sh
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/run.sh
@@ -6,7 +6,7 @@ then
./getIC.sh
fi
-if [ ! -e cooling_tables/ ]
+if [ ! -e cooling_tables ]
then
echo "Fetching EAGLE cooling tables for the isolated galaxy example..."
./getEagleCoolingTable.sh
diff --git a/examples/main.c b/examples/main.c
index 5b78a2eed8d7e7d93a396fbdbd71553ac7decbb2..f5622b0615871d72ffadb8dedb50b0de4271fa31 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -96,6 +96,7 @@ int main(int argc, char *argv[]) {
struct gravity_props gravity_properties;
struct hydro_props hydro_properties;
struct stars_props stars_properties;
+ struct feedback_props feedback_properties;
struct entropy_floor_properties entropy_floor;
struct part *parts = NULL;
struct phys_const prog_const;
@@ -728,6 +729,17 @@ int main(int argc, char *argv[]) {
else
bzero(&stars_properties, sizeof(struct stars_props));
+ /* Initialise the feedback properties */
+ if (with_feedback) {
+#ifdef FEEDBACK_NONE
+ error("ERROR: Running with feedback but compiled without it.");
+#endif
+ feedback_props_init(&feedback_properties, &prog_const, &us, params,
+ &hydro_properties, &cosmo);
+ }
+ else
+ bzero(&feedback_properties, sizeof(struct feedback_props));
+
/* Initialise the gravity properties */
if (with_self_gravity)
gravity_props_init(&gravity_properties, params, &cosmo, with_cosmology,
@@ -768,10 +780,6 @@ int main(int argc, char *argv[]) {
chemistry_init(params, &us, &prog_const, &chemistry);
if (myrank == 0) chemistry_print(&chemistry);
- /* Initialise stellar evolution */
- // MATTHIEU
- // if (with_feedback) stars_evolve_init(params, &stars_properties);
-
/* Be verbose about what happens next */
if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
if (myrank == 0 && cleanup_h)
@@ -961,8 +969,8 @@ int main(int argc, char *argv[]) {
engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2],
engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
&hydro_properties, &entropy_floor, &gravity_properties,
- &stars_properties, &mesh, &potential, &cooling_func, &starform,
- &chemistry);
+ &stars_properties, &feedback_properties, &mesh, &potential,
+ &cooling_func, &starform, &chemistry);
engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
talking, restart_file);
diff --git a/src/engine.c b/src/engine.c
index e0747dde19f44f20b370b2d1d806f4a7e757bb3c..a2ad8816e88f94b3b55f6d9d01e691f04b5cfc02 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -4404,6 +4404,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
struct cosmology *cosmo, struct hydro_props *hydro,
const struct entropy_floor_properties *entropy_floor,
struct gravity_props *gravity, const struct stars_props *stars,
+ const struct feedback_props *feedback,
struct pm_mesh *mesh,
const struct external_potential *potential,
struct cooling_function_data *cooling_func,
@@ -4475,6 +4476,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
e->external_potential = potential;
e->cooling_func = cooling_func;
e->star_formation = starform;
+ e->feedback_props = feedback;
e->chemistry = chemistry;
e->parameter_file = params;
#ifdef WITH_MPI
diff --git a/src/engine.h b/src/engine.h
index a4058befcc14097a23928018553882da6da443cc..88f2ae4f9d510c9701ae16b4a769982c58ceb522 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -447,6 +447,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
struct cosmology *cosmo, struct hydro_props *hydro,
const struct entropy_floor_properties *entropy_floor,
struct gravity_props *gravity, const struct stars_props *stars,
+ const struct feedback_props *feedback,
struct pm_mesh *mesh,
const struct external_potential *potential,
struct cooling_function_data *cooling_func,
diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c
index ed57266bb373d9a4c24431151261902ebc46f88f..3cf1313c52d1a68e787727e2126935b26f9a5a04 100644
--- a/src/feedback/EAGLE/feedback.c
+++ b/src/feedback/EAGLE/feedback.c
@@ -25,6 +25,29 @@
#include "interpolate.h"
#include "yield_tables.h"
+/**
+ * @brief Compute number of SN that should go off given the age of the spart
+ *
+ * @param sp spart we're evolving
+ * @param stars_properties stars_props data structure
+ * @param age age of star at the beginning of the step
+ * @param dt length of step
+ */
+float compute_SNe(const struct spart* sp,
+ const struct feedback_props* feedback_props,
+ const float age, const float dt) {
+
+ const float SNII_wind_delay = feedback_props->SNII_wind_delay;
+
+ if (age <= SNII_wind_delay && (age + dt) > SNII_wind_delay) {
+
+ return feedback_props->num_SNII_per_msun * sp->mass_init /
+ feedback_props->const_solar_mass;
+ } else {
+ return 0.f;
+ }
+}
+
/**
* @brief determine which metallicity bin star belongs to for AGB, compute bin
* indices and offsets
@@ -218,6 +241,7 @@ inline static void evolve_SNII(float log10_min_mass, float log10_max_mass,
float* stellar_yields,
const struct feedback_props* restrict feedback_props,
struct spart* restrict sp) {
+
int low_imf_mass_bin_index, high_imf_mass_bin_index, mass_bin_index;
/* If mass at beginning of step is less than tabulated lower bound for IMF,
@@ -376,6 +400,7 @@ inline static void evolve_AGB(float log10_min_mass, float log10_max_mass,
float* stellar_yields,
const struct feedback_props* restrict feedback_props,
struct spart* restrict sp) {
+
int low_imf_mass_bin_index, high_imf_mass_bin_index, mass_bin_index;
/* If mass at end of step is greater than tabulated lower bound for IMF, limit
@@ -517,10 +542,10 @@ inline static void evolve_AGB(float log10_min_mass, float log10_max_mass,
* @param age age of spart at beginning of step
* @param dt length of current timestep
*/
-void compute_stellar_evolution(
- const struct feedback_props* feedback_props,
- struct spart* sp, const struct unit_system* us, float age,
- double dt) {
+void compute_stellar_evolution(const struct feedback_props* feedback_props,
+ const struct cosmology *cosmo, struct spart* sp,
+ const struct unit_system* us, const float age,
+ const float dt) {
/* Allocate temporary array for calculating imf weights */
float* stellar_yields;
@@ -528,149 +553,71 @@ void compute_stellar_evolution(
malloc(feedback_props->n_imf_mass_bins * sizeof(float));
/* Convert dt and stellar age from internal units to Gyr. */
- const double Gyr_in_cgs = 3.154e16;
- double dt_Gyr =
- dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME) / Gyr_in_cgs;
- double star_age_Gyr =
- age * units_cgs_conversion_factor(us, UNIT_CONV_TIME) / Gyr_in_cgs;
-
+ const double Gyr_in_cgs = 1e9 * 365. * 24. * 3600.;
+ const double time_to_cgs = units_cgs_conversion_factor(us, UNIT_CONV_TIME);
+ const float conversion_factor = time_to_cgs / Gyr_in_cgs;
+ const float dt_Gyr = dt * conversion_factor;
+ const float star_age_Gyr = age * conversion_factor;
+
+ /* Get the metallicity */
+ const float Z = sp->chemistry_data.metal_mass_fraction_total;
+
/* calculate mass of stars that has died from the star's birth up to the
* beginning and end of timestep */
- double log10_max_dying_mass_msun = log10(dying_mass_msun(
- star_age_Gyr, sp->chemistry_data.metal_mass_fraction_total,
- feedback_props));
- double log10_min_dying_mass_msun = log10(dying_mass_msun(
- star_age_Gyr + dt_Gyr, sp->chemistry_data.metal_mass_fraction_total,
- feedback_props));
+ const float max_dying_mass_Msun = dying_mass_msun(star_age_Gyr, Z, feedback_props);
+ const float min_dying_mass_Msun = dying_mass_msun(star_age_Gyr + dt_Gyr, Z, feedback_props);
+#ifdef SWIFT_DEBUG_CHECK
/* Sanity check. Worth investigating if necessary as functions for evaluating
* mass of stars dying might be strictly decreasing. */
- if (log10_min_dying_mass_msun > log10_max_dying_mass_msun)
+ if (min_dying_mass_Msun > max_dying_mass_Msun)
error("min dying mass is greater than max dying mass");
-
+#endif
+
/* Integration interval is zero - this can happen if minimum and maximum
- * dying masses are above imf_max_mass_msun. Return without doing any
+ * dying masses are above imf_max_mass_Msun. Return without doing any
* feedback. */
- if (log10_min_dying_mass_msun == log10_max_dying_mass_msun) return;
+ if (min_dying_mass_Msun == max_dying_mass_Msun) return;
- /* Evolve SNIa, SNII, AGB */
- evolve_SNIa(log10_min_dying_mass_msun, log10_max_dying_mass_msun,
+ /* Life is better in log */
+ const float log10_max_dying_mass_Msun = log10f(max_dying_mass_Msun);
+ const float log10_min_dying_mass_Msun = log10f(min_dying_mass_Msun);
+
+ /* Compute elements, energy and momentum to distribute from the
+ * three channels SNIa, SNII, AGB */
+ evolve_SNIa(log10_min_dying_mass_Msun, log10_max_dying_mass_Msun,
feedback_props, sp, star_age_Gyr, dt_Gyr);
- evolve_SNII(log10_min_dying_mass_msun, log10_max_dying_mass_msun,
+ evolve_SNII(log10_min_dying_mass_Msun, log10_max_dying_mass_Msun,
stellar_yields, feedback_props, sp);
- evolve_AGB(log10_min_dying_mass_msun, log10_max_dying_mass_msun,
+ evolve_AGB(log10_min_dying_mass_Msun, log10_max_dying_mass_Msun,
stellar_yields, feedback_props, sp);
- /* Compute the mass to distribute */
+ /* Compute the total mass to distribute (H + He metals) */
sp->feedback_data.to_distribute.mass = sp->feedback_data.to_distribute.total_metal_mass +
sp->feedback_data.to_distribute.metal_mass[chemistry_element_H] +
sp->feedback_data.to_distribute.metal_mass[chemistry_element_He];
- /* Clean up */
- free(stellar_yields);
-}
+ /* Compute the number of type II SNe that went off */
+ sp->feedback_data.to_distribute.num_SNe = compute_SNe(sp, feedback_props, age, dt);
-/**
- * @brief Compute number of SN that should go off given the age of the spart
- *
- * @param sp spart we're evolving
- * @param stars_properties stars_props data structure
- * @param age age of star at the beginning of the step
- * @param dt length of step
- */
-float compute_SNe(struct spart* sp,
- const struct feedback_props* feedback_props,
- float age, double dt) {
-
- if (age <= feedback_props->SNII_wind_delay &&
- age + dt > feedback_props->SNII_wind_delay) {
-
- return feedback_props->num_SNII_per_msun * sp->mass_init /
- feedback_props->const_solar_mass;
- } else {
- return 0;
- }
-}
+ /* Compute energy change due to thermal and kinetic energy of ejecta */
+ sp->feedback_data.to_distribute.d_energy =
+ sp->feedback_data.to_distribute.mass *
+ (feedback_props->ejecta_specific_thermal_energy +
+ 0.5 * (sp->v[0] * sp->v[0] + sp->v[1] * sp->v[1] + sp->v[2] * sp->v[2]) *
+ cosmo->a2_inv);
+ /* Compute probability of heating neighbouring particles */
+ if (dt > 0 && sp->feedback_data.ngb_mass > 0)
+ sp->feedback_data.to_distribute.heating_probability =
+ feedback_props->total_energy_SNe *
+ sp->feedback_data.to_distribute.num_SNe /
+ (feedback_props->temp_to_u_factor *
+ feedback_props->SNe_deltaT_desired * sp->feedback_data.ngb_mass);
-/**
- * @brief Initializes constants related to stellar evolution, initializes imf,
- * reads and processes yield tables
- *
- * @param params swift_params parameters structure
- * @param stars stars_props data structure
- */
-inline static void stars_evolve_init(struct swift_params* params,
- struct feedback_props* feedback_props) {
-
- /* Set number of elements found in yield tables */
- feedback_props->SNIa_n_elements = 42;
- feedback_props->SNII_n_mass = 11;
- feedback_props->SNII_n_elements = 11;
- feedback_props->SNII_n_z = 5;
- feedback_props->AGB_n_mass = 23;
- feedback_props->AGB_n_elements = 11;
- feedback_props->AGB_n_z = 3;
- feedback_props->lifetimes.n_mass = 30;
- feedback_props->lifetimes.n_z = 6;
- feedback_props->element_name_length = 15;
-
- /* Set bounds for imf */
- feedback_props->log10_SNII_min_mass_msun = 0.77815125f; // log10(6).
- feedback_props->log10_SNII_max_mass_msun = 2.f; // log10(100).
- feedback_props->log10_SNIa_max_mass_msun = 0.90308999f; // log10(8).
-
- /* Yield table filepath */
- parser_get_param_string(params, "EAGLEFeedback:filename",
- feedback_props->yield_table_path);
- parser_get_param_string(params, "EAGLEFeedback:imf_model",
- feedback_props->IMF_Model);
-
- /* Initialise IMF */
- init_imf(feedback_props);
-
- /* Allocate yield tables */
- allocate_yield_tables(feedback_props);
-
- /* Set factors for each element adjusting SNII yield */
- feedback_props->typeII_factor[0] = 1.f;
- feedback_props->typeII_factor[1] = 1.f;
- feedback_props->typeII_factor[2] = 0.5f;
- feedback_props->typeII_factor[3] = 1.f;
- feedback_props->typeII_factor[4] = 1.f;
- feedback_props->typeII_factor[5] = 1.f;
- feedback_props->typeII_factor[6] = 2.f;
- feedback_props->typeII_factor[7] = 1.f;
- feedback_props->typeII_factor[8] = 0.5f;
-
- /* Read the tables */
- read_yield_tables(feedback_props);
-
- /* Set yield_mass_bins array */
- const float imf_log10_mass_bin_size =
- (feedback_props->log10_imf_max_mass_msun -
- feedback_props->log10_imf_min_mass_msun) /
- (feedback_props->n_imf_mass_bins - 1);
- for (int i = 0; i < feedback_props->n_imf_mass_bins; i++)
- feedback_props->yield_mass_bins[i] =
- imf_log10_mass_bin_size * i + feedback_props->log10_imf_min_mass_msun;
-
- /* Resample yields from mass bins used in tables to mass bins used in IMF */
- compute_yields(feedback_props);
-
- /* Resample ejecta contribution to enrichment from mass bins used in tables to
- * mass bins used in IMF */
- compute_ejecta(feedback_props);
-
- /* Calculate number of type II SN per solar mass
- * Note: since we are integrating the IMF without weighting it by the yields
- * pass NULL pointer for stellar_yields array */
- feedback_props->num_SNII_per_msun =
- integrate_imf(feedback_props->log10_SNII_min_mass_msun,
- feedback_props->log10_SNII_max_mass_msun, 0,
- /*(stellar_yields=)*/ NULL, feedback_props);
-
- message("initialized stellar feedback");
+
+ /* Clean up */
+ free(stellar_yields);
}
/**
@@ -684,12 +631,12 @@ inline static void stars_evolve_init(struct swift_params* params,
* @param params The parsed parameters.
* @param p The already read-in properties of the hydro scheme.
*/
-void feedbakc_props_init(struct feedback_props *fp,
- const struct phys_const *phys_const,
- const struct unit_system *us,
- struct swift_params *params,
- const struct hydro_props *p,
- const struct cosmology *cosmo) {
+void feedback_props_init(struct feedback_props *fp,
+ const struct phys_const *phys_const,
+ const struct unit_system *us,
+ struct swift_params *params,
+ const struct hydro_props *hydro_props,
+ const struct cosmology *cosmo) {
/* Read SNIa timscale */
fp->SNIa_timescale_Gyr =
@@ -728,7 +675,7 @@ void feedbakc_props_init(struct feedback_props *fp,
/* Calculate temperature to internal energy conversion factor */
fp->temp_to_u_factor =
phys_const->const_boltzmann_k /
- (p->mu_ionised * (hydro_gamma_minus_one)*phys_const->const_proton_mass);
+ (hydro_props->mu_ionised * (hydro_gamma_minus_one)*phys_const->const_proton_mass);
/* Read birth time to set all stars in ICs to (defaults to -1 to indicate star
* present in ICs) */
@@ -737,6 +684,76 @@ void feedbakc_props_init(struct feedback_props *fp,
/* Copy over solar mass */
fp->const_solar_mass = phys_const->const_solar_mass;
+
+
+ /* Set number of elements found in yield tables */
+ fp->SNIa_n_elements = 42;
+ fp->SNII_n_mass = 11;
+ fp->SNII_n_elements = 11;
+ fp->SNII_n_z = 5;
+ fp->AGB_n_mass = 23;
+ fp->AGB_n_elements = 11;
+ fp->AGB_n_z = 3;
+ fp->lifetimes.n_mass = 30;
+ fp->lifetimes.n_z = 6;
+ fp->element_name_length = 15;
+
+ /* Set bounds for imf */
+ fp->log10_SNII_min_mass_msun = 0.77815125f; // log10(6).
+ fp->log10_SNII_max_mass_msun = 2.f; // log10(100).
+ fp->log10_SNIa_max_mass_msun = 0.90308999f; // log10(8).
+
+ /* Yield table filepath */
+ parser_get_param_string(params, "EAGLEFeedback:filename",
+ fp->yield_table_path);
+ parser_get_param_string(params, "EAGLEFeedback:imf_model",
+ fp->IMF_Model);
+
+ /* Initialise IMF */
+ init_imf(fp);
+
+ /* Allocate yield tables */
+ allocate_yield_tables(fp);
+
+ /* Set factors for each element adjusting SNII yield */
+ fp->typeII_factor[0] = 1.f;
+ fp->typeII_factor[1] = 1.f;
+ fp->typeII_factor[2] = 0.5f;
+ fp->typeII_factor[3] = 1.f;
+ fp->typeII_factor[4] = 1.f;
+ fp->typeII_factor[5] = 1.f;
+ fp->typeII_factor[6] = 2.f;
+ fp->typeII_factor[7] = 1.f;
+ fp->typeII_factor[8] = 0.5f;
+
+ /* Read the tables */
+ read_yield_tables(fp);
+
+ /* Set yield_mass_bins array */
+ const float imf_log10_mass_bin_size =
+ (fp->log10_imf_max_mass_msun -
+ fp->log10_imf_min_mass_msun) /
+ (fp->n_imf_mass_bins - 1);
+ for (int i = 0; i < fp->n_imf_mass_bins; i++)
+ fp->yield_mass_bins[i] =
+ imf_log10_mass_bin_size * i + fp->log10_imf_min_mass_msun;
+
+ /* Resample yields from mass bins used in tables to mass bins used in IMF */
+ compute_yields(fp);
+
+ /* Resample ejecta contribution to enrichment from mass bins used in tables to
+ * mass bins used in IMF */
+ compute_ejecta(fp);
+
+ /* Calculate number of type II SN per solar mass
+ * Note: since we are integrating the IMF without weighting it by the yields
+ * pass NULL pointer for stellar_yields array */
+ fp->num_SNII_per_msun =
+ integrate_imf(fp->log10_SNII_min_mass_msun,
+ fp->log10_SNII_max_mass_msun, 0,
+ /*(stellar_yields=)*/ NULL, fp);
+
+ message("initialized stellar feedback");
}
diff --git a/src/feedback/EAGLE/feedback.h b/src/feedback/EAGLE/feedback.h
index d4ea9cb35d100588242343be6738bc87e61c7d9d..ea79dc96ed3581ed89267fd9754d22e0d47c7dc7 100644
--- a/src/feedback/EAGLE/feedback.h
+++ b/src/feedback/EAGLE/feedback.h
@@ -20,20 +20,17 @@
#define SWIFT_FEEDBACK_EAGLE_H
#include "cosmology.h"
+#include "hydro_properties.h"
#include "units.h"
#include "part.h"
#include "feedback_properties.h"
-void compute_stellar_evolution(
- const struct feedback_props* feedback_props,
- struct spart* sp, const struct unit_system* us, float age,
- double dt);
-
-float compute_SNe(struct spart* sp,
- const struct feedback_props* stars_properties,
- float age, double dt);
+void compute_stellar_evolution(const struct feedback_props* feedback_props,
+ const struct cosmology *cosmo,
+ struct spart* sp, const struct unit_system* us,
+ const float age, const float dt);
/**
* @brief Prepares a s-particle for its feedback interactions
@@ -59,8 +56,8 @@ __attribute__((always_inline)) INLINE static void feedback_init_spart(
__attribute__((always_inline)) INLINE static void feedback_first_init_spart(
struct spart* sp, const struct feedback_props* feedback_props) {
- //sp->birth_density = -1.f;
- //sp->birth_time = feedback_props.spart_first_init_birth_time;
+ sp->birth_density = -1.f;
+ sp->birth_time = feedback_props->spart_first_init_birth_time;
feedback_init_spart(sp);
}
@@ -111,28 +108,10 @@ __attribute__((always_inline)) INLINE static void stars_evolve_spart(
/* Compute amount of enrichment and feedback that needs to be done in this
* step */
- compute_stellar_evolution(feedback_props, sp, us, star_age, dt);
+ compute_stellar_evolution(feedback_props, cosmo, sp, us, star_age, dt);
/* Decrease star mass by amount of mass distributed to gas neighbours */
sp->mass -= sp->feedback_data.to_distribute.mass;
-
- /* Compute the number of type II SNe that went off */
- sp->feedback_data.to_distribute.num_SNe = compute_SNe(sp, feedback_props, star_age, dt);
-
- /* Compute energy change due to thermal and kinetic energy of ejecta */
- sp->feedback_data.to_distribute.d_energy =
- sp->feedback_data.to_distribute.mass *
- (feedback_props->ejecta_specific_thermal_energy +
- 0.5 * (sp->v[0] * sp->v[0] + sp->v[1] * sp->v[1] + sp->v[2] * sp->v[2]) *
- cosmo->a2_inv);
-
- /* Compute probability of heating neighbouring particles */
- if (dt > 0 && sp->feedback_data.ngb_mass > 0)
- sp->feedback_data.to_distribute.heating_probability =
- feedback_props->total_energy_SNe *
- sp->feedback_data.to_distribute.num_SNe /
- (feedback_props->temp_to_u_factor *
- feedback_props->SNe_deltaT_desired * sp->feedback_data.ngb_mass);
}
diff --git a/src/feedback/EAGLE/feedback_iact.h b/src/feedback/EAGLE/feedback_iact.h
index ee13c8951e749bc27d4b85846dbd38e1dbcda413..3533e13b462e8eae034ca48dad2758b3ee27549a 100644
--- a/src/feedback/EAGLE/feedback_iact.h
+++ b/src/feedback/EAGLE/feedback_iact.h
@@ -1,8 +1,6 @@
/**
* @brief Density interaction between two particles (non-symmetric).
- * MATTHIEU: check with RGB about concerns with comment (in this and other
- * subgrid schemes)
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
diff --git a/src/feedback/EAGLE/feedback_properties.h b/src/feedback/EAGLE/feedback_properties.h
index 0389edee86ec30994a3c81265f70a49a488222d0..b52bfddc6f3bb9ab0c50186a7a4ab903e5656e6b 100644
--- a/src/feedback/EAGLE/feedback_properties.h
+++ b/src/feedback/EAGLE/feedback_properties.h
@@ -150,3 +150,11 @@ struct feedback_props {
float spart_first_init_birth_time;
};
+
+
+void feedback_props_init(struct feedback_props *fp,
+ const struct phys_const *phys_const,
+ const struct unit_system *us,
+ struct swift_params *params,
+ const struct hydro_props *hydro_props,
+ const struct cosmology *cosmo);
diff --git a/src/feedback/EAGLE/imf.h b/src/feedback/EAGLE/imf.h
index 9914adacae2495337ca1406275a9e906fcd98289..355a97c3de3988def133c6eff215967655aa673a 100644
--- a/src/feedback/EAGLE/imf.h
+++ b/src/feedback/EAGLE/imf.h
@@ -274,194 +274,211 @@ inline static void init_imf(struct feedback_props *restrict feedback_props) {
* @brief Calculate mass (in solar masses) of stars that died from the star
* particle's birth up to its current age (in Gyr).
*
- * Calculation follows Portinari et al. 1998
+ * Calculation uses the tables of Portinari et al. 1998, A&A, 334, 505
*
- * @param age_Gyr age of star in Gyr
- * @param metallicity Star's metallicity
- * @param star_properties the #stars_props data structure
+ * @param age_Gyr age of star in Gyr.
+ * @param Z Star's metallicity (metal mass fraction).
+ * @param star_properties the #stars_props data structure.
+ * @return Mass of stars died up to that age in solar masses.
*/
-inline static double dying_mass_msun(const float age_Gyr,
- const float metallicity,
- const struct feedback_props *star_props) {
+inline static float dying_mass_msun(const float age_Gyr, const float Z,
+ const struct feedback_props *feedback_props) {
- float metallicity_offset, time_offset_low_metallicity = 0,
- time_offset_high_metallicity = 0,
- mass_low_metallicity, mass_high_metallicity;
+ /* Pull out some common terms */
+ const double *lifetime_Z = feedback_props->lifetimes.metallicity;
+ const double *lifetime_m = feedback_props->lifetimes.mass;
+ double **const dying_times = feedback_props->lifetimes.dyingtime;
+ const int n_Z = feedback_props->lifetimes.n_z;
+ const int n_m = feedback_props->lifetimes.n_mass;
+
+ /* Early abort? */
+ if (age_Gyr <= 0.f) {
+ return feedback_props->imf_max_mass_msun;
+ }
- int metallicity_index, time_index_low_metallicity = -1,
- time_index_high_metallicity = -1;
+ const float log10_age_yr = log10f(age_Gyr * 1e9f);
- if (age_Gyr <= 0) {
- return star_props->imf_max_mass_msun;
- }
+ /* Calculate index along the metallicity axis */
+ int Z_index;
+ float Z_offset;
+ if (Z <= lifetime_Z[0]) {
- const float log10_age_yr = log10f(age_Gyr * 1.e9);
+ /* Before start of the table */
+ Z_index = 0;
+ Z_offset = 0.f;
+
+ } else if (Z >= lifetime_Z[n_Z - 1]) {
- if (metallicity <= star_props->lifetimes.metallicity[0]) {
- metallicity_index = 0;
- metallicity_offset = 0.0;
- } else if (metallicity >=
- star_props->lifetimes
- .metallicity[star_props->lifetimes.n_z - 1]) {
- metallicity_index = star_props->lifetimes.n_z - 2;
- metallicity_offset = 1.0;
+ /* After end of the table */
+ Z_index = n_Z - 2;
+ Z_offset = 1.f;
+
} else {
- for (metallicity_index = 0;
- metallicity_index < star_props->lifetimes.n_z - 1;
- metallicity_index++)
- if (star_props->lifetimes.metallicity[metallicity_index + 1] >
- metallicity)
+
+ /* Normal case: Somewhere inside the table */
+ for (Z_index = 0; Z_index < n_Z - 1; Z_index++) {
+ if (lifetime_Z[Z_index + 1] > Z)
break;
+ }
- metallicity_offset =
- (metallicity -
- star_props->lifetimes.metallicity[metallicity_index]) /
- (star_props->lifetimes.metallicity[metallicity_index + 1] -
- star_props->lifetimes.metallicity[metallicity_index]);
+ Z_offset = (Z - lifetime_Z[Z_index]) /
+ (lifetime_Z[Z_index + 1] - lifetime_Z[Z_index]);
}
- if (log10_age_yr >=
- star_props->lifetimes.dyingtime[metallicity_index][0]) {
- time_index_low_metallicity = 0;
- time_offset_low_metallicity = 0.0;
- } else if (log10_age_yr <=
- star_props->lifetimes
- .dyingtime[metallicity_index]
- [star_props->lifetimes.n_mass - 1]) {
- time_index_low_metallicity = star_props->lifetimes.n_mass - 2;
- time_offset_low_metallicity = 1.0;
+ /* Check whether we are not beyond the table */
+ int time_index_lowZ = -1;
+ float time_offset_lowZ = 0.f;
+ if (log10_age_yr >= dying_times[Z_index][0]) {
+
+ /* Before start of the table */
+ time_index_lowZ = 0;
+ time_offset_lowZ = 0.f;
+
+ } else if (log10_age_yr <= dying_times[Z_index][n_m - 1]) {
+
+ /* After end of the table */
+ time_index_lowZ = n_m - 2;
+ time_offset_lowZ = 1.0;
}
- if (log10_age_yr >=
- star_props->lifetimes.dyingtime[metallicity_index + 1][0]) {
- time_index_high_metallicity = 0;
- time_offset_high_metallicity = 0.0;
- } else if (log10_age_yr <=
- star_props->lifetimes
- .dyingtime[metallicity_index + 1]
- [star_props->lifetimes.n_mass - 1]) {
- time_index_high_metallicity = star_props->lifetimes.n_mass - 2;
- time_offset_high_metallicity = 1.0;
+ /* Check whether we are not beyond the table */
+ int time_index_highZ = -1;
+ float time_offset_highZ = 0.f;
+ if (log10_age_yr >= dying_times[Z_index + 1][0]) {
+
+ /* Before start of the table */
+ time_index_highZ = 0;
+ time_offset_highZ = 0.f;
+
+ } else if (log10_age_yr <= dying_times[Z_index + 1][n_m - 1]) {
+
+ /* After end of the table */
+ time_index_highZ = n_m - 2;
+ time_offset_highZ = 1.0;
}
- int i = star_props->lifetimes.n_mass;
- while (i >= 0 && (time_index_low_metallicity == -1 ||
- time_index_high_metallicity == -1)) {
+ /* Search the table starting from the largest times until we reach
+ a solution for the two bounds */
+ int i = n_m;
+ while (i >= 0 && (time_index_lowZ == -1 || time_index_highZ == -1)) {
+
i--;
- if (star_props->lifetimes.dyingtime[metallicity_index][i] >=
- log10_age_yr &&
- time_index_low_metallicity == -1) {
- time_index_low_metallicity = i;
- time_offset_low_metallicity =
- (log10_age_yr -
- star_props->lifetimes
- .dyingtime[metallicity_index][time_index_low_metallicity]) /
- (star_props->lifetimes
- .dyingtime[metallicity_index][time_index_low_metallicity + 1] -
- star_props->lifetimes
- .dyingtime[metallicity_index][time_index_low_metallicity]);
+
+ if (dying_times[Z_index][i] >= log10_age_yr && time_index_lowZ == -1) {
+
+ /* record index */
+ time_index_lowZ = i;
+
+ /* record distance from table element */
+ time_offset_lowZ =
+ (log10_age_yr - dying_times[Z_index][time_index_lowZ]) /
+ (dying_times[Z_index][time_index_lowZ + 1] -
+ dying_times[Z_index][time_index_lowZ]);
}
- if (star_props->lifetimes.dyingtime[metallicity_index + 1][i] >=
- log10_age_yr &&
- time_index_high_metallicity == -1) {
- time_index_high_metallicity = i;
- time_offset_high_metallicity =
- (log10_age_yr -
- star_props->lifetimes
- .dyingtime[metallicity_index + 1][time_index_high_metallicity]) /
- (star_props->lifetimes
- .dyingtime[metallicity_index + 1]
- [time_index_high_metallicity + 1] -
- star_props->lifetimes
- .dyingtime[metallicity_index + 1][time_index_high_metallicity]);
+
+ if (dying_times[Z_index + 1][i] >= log10_age_yr && time_index_highZ == -1) {
+
+ /* record index */
+ time_index_highZ = i;
+
+ /* record distance from table element */
+ time_offset_highZ =
+ (log10_age_yr - dying_times[Z_index + 1][time_index_highZ]) /
+ (dying_times[Z_index + 1][time_index_highZ + 1] -
+ dying_times[Z_index + 1][time_index_highZ]);
}
}
- mass_low_metallicity =
- interpolate_1d(star_props->lifetimes.mass,
- time_index_low_metallicity, time_offset_low_metallicity);
- mass_high_metallicity =
- interpolate_1d(star_props->lifetimes.mass,
- time_index_high_metallicity, time_offset_high_metallicity);
+ /* And now interpolate the solution */
+ const float mass_low_Z = interpolate_1d(lifetime_m, time_index_lowZ, time_offset_lowZ);
+ const float mass_high_Z = interpolate_1d(lifetime_m, time_index_highZ, time_offset_highZ);
- float mass = (1.0 - metallicity_offset) * mass_low_metallicity +
- metallicity_offset * mass_high_metallicity;
+ float mass = (1.f - Z_offset) * mass_low_Z + Z_offset * mass_high_Z;
/* Check that we haven't killed too many stars */
- if (mass > star_props->imf_max_mass_msun)
- mass = star_props->imf_max_mass_msun;
+ mass = min(mass, feedback_props->imf_max_mass_msun);
return mass;
}
/**
- * @brief Calculate lifetime of star poputlation in Gyr. Approach based on
- * Portinari et al. 1998
+ * @brief Calculate lifetime of stellar population in Gyr for a given mass.
*
- * @param mass
- * @param metallicity
- * @param star_properties the #stars_props data struct
+ * Calculation uses the tables of Portinari et al. 1998, A&A, 334, 505
+ *
+ * @param mass in solar masses.
+ * @param Z Metallicity (metal mass fraction).
+ * @param feedback_props the #feedback_props data structure.
+ * @return The life time in Giga-years.
*/
-inline static float lifetime_in_Gyr(float mass, float metallicity,
- const struct feedback_props *restrict
- feedback_props) {
-
- double time_Gyr = 0, mass_offset, metallicity_offset;
-
- int mass_index, metallicity_index;
-
- if (mass <= feedback_props->lifetimes.mass[0]) {
- mass_index = 0;
- mass_offset = 0.0;
- } else if (mass >=
- feedback_props->lifetimes
- .mass[feedback_props->lifetimes.n_mass - 1]) {
- mass_index = feedback_props->lifetimes.n_mass - 2;
- mass_offset = 1.0;
+inline static float lifetime_in_Gyr(const float mass, const float Z,
+ const struct feedback_props *feedback_props) {
+
+ /* Pull out some common terms */
+ const double *lifetime_Z = feedback_props->lifetimes.metallicity;
+ const double *lifetime_m = feedback_props->lifetimes.mass;
+ double **const dying_times = feedback_props->lifetimes.dyingtime;
+ const int n_Z = feedback_props->lifetimes.n_z;
+ const int n_m = feedback_props->lifetimes.n_mass;
+
+ /* Calculate index along the mass axis */
+ int m_index;
+ float m_offset;
+ if (mass <= lifetime_m[0]) {
+
+ /* Before start of the table */
+ m_index = 0;
+ m_offset = 0.f;
+
+ } else if (mass >= lifetime_m[n_m - 1]) {
+
+ /* After end of the table */
+ m_index = n_m - 2;
+ m_offset = 1.f;
+
} else {
- for (mass_index = 0;
- mass_index < feedback_props->lifetimes.n_mass - 1;
- mass_index++)
- if (feedback_props->lifetimes.mass[mass_index + 1] > mass)
- break;
- mass_offset =
- (mass - feedback_props->lifetimes.mass[mass_index]) /
- (feedback_props->lifetimes.mass[mass_index + 1] -
- feedback_props->lifetimes.mass[mass_index]);
+ /* Normal case: Somewhere inside the table */
+ for (m_index = 0; m_index < n_m - 1; m_index++)
+ if (lifetime_m[m_index + 1] > mass)
+ break;
+
+ m_offset = (mass - lifetime_m[m_index]) /
+ (lifetime_m[m_index + 1] - lifetime_m[m_index]);
}
- if (metallicity <= feedback_props->lifetimes.metallicity[0]) {
- metallicity_index = 0;
- metallicity_offset = 0.0;
- } else if (metallicity >=
- feedback_props->lifetimes
- .metallicity[feedback_props->lifetimes.n_z - 1]) {
- metallicity_index = feedback_props->lifetimes.n_z - 2;
- metallicity_offset = 1.0;
+ /* Calculate index along the metallicity axis */
+ int Z_index;
+ float Z_offset;
+ if (Z <= lifetime_Z[0]) {
+
+ /* Before start of the table */
+ Z_index = 0;
+ Z_offset = 0.f;
+
+ } else if (Z >= lifetime_Z[n_Z - 1]) {
+
+ /* After end of the table */
+ Z_index = n_Z - 2;
+ Z_offset = 1.f;
+
} else {
- for (metallicity_index = 0;
- metallicity_index < feedback_props->lifetimes.n_z - 1;
- metallicity_index++)
- if (feedback_props->lifetimes
- .metallicity[metallicity_index + 1] > metallicity)
+
+ for (Z_index = 0; Z_index < n_Z - 1; Z_index++)
+ if (lifetime_Z[Z_index + 1] > Z)
break;
- metallicity_offset =
- (metallicity -
- feedback_props->lifetimes.metallicity[metallicity_index]) /
- (feedback_props->lifetimes
- .metallicity[metallicity_index + 1] -
- feedback_props->lifetimes.metallicity[metallicity_index]);
+ /* Normal case: Somewhere inside the table */
+ Z_offset = (Z - lifetime_Z[Z_index]) /
+ (lifetime_Z[Z_index + 1] - lifetime_Z[Z_index]);
}
- /* time in Gyr */
- time_Gyr =
- exp(M_LN10 * interpolate_2d(feedback_props->lifetimes.dyingtime,
- metallicity_index, mass_index,
- metallicity_offset, mass_offset)) /
- 1.0e9;
+ /* Interpolation of the table to get the time */
+ const float log_time_years = interpolate_2d(dying_times, Z_index, m_index, Z_offset, m_offset);
+ /* Convert to Giga-years */
+ const float time_Gyr = exp10f(log_time_years - 9.f);
+
return time_Gyr;
}
diff --git a/src/feedback/EAGLE/interpolate.h b/src/feedback/EAGLE/interpolate.h
index 5cd6a90df2dd459e1a592cde527eb60e92a55d09..a4b8ff1615a6146ef7084ea0f8c324d595a8c520 100644
--- a/src/feedback/EAGLE/interpolate.h
+++ b/src/feedback/EAGLE/interpolate.h
@@ -33,8 +33,8 @@ __attribute__((always_inline)) static INLINE int row_major_index_2d(const int i,
const int Nx,
const int Ny) {
#ifdef SWIFT_DEBUG_CHECKS
- assert(x < Nx);
- assert(y < Ny);
+ assert(i < Nx);
+ assert(j < Ny);
#endif
return i * Ny + j;
@@ -51,9 +51,9 @@ __attribute__((always_inline)) static INLINE int row_major_index_3d(
const int i, const int j, const int k, const int Nx, const int Ny, const int Nz) {
#ifdef SWIFT_DEBUG_CHECKS
- assert(x < Nx);
- assert(y < Ny);
- assert(z < Nz);
+ assert(i < Nx);
+ assert(j < Ny);
+ assert(k < Nz);
#endif
return i * Ny * Nz + j * Nz + k;
diff --git a/src/feedback_properties.h b/src/feedback_properties.h
new file mode 100644
index 0000000000000000000000000000000000000000..70a5feb41bd69339fb10a235bd915d67d8391cf2
--- /dev/null
+++ b/src/feedback_properties.h
@@ -0,0 +1,34 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_FEEDBACK_PROPERTIES_H
+#define SWIFT_FEEDBACK_PROPERTIES_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Select the correct feedback model */
+#if defined(FEEDBACK_NONE)
+#include "./feedback/Default/feedback_properties.h"
+#elif defined(FEEDBACK_EAGLE)
+#include "./feedback/EAGLE/feedback_properties.h"
+#else
+#error "Invalid choice of feedback model"
+#endif
+
+#endif /* SWIFT_FEEDBACK_PROPERTIES_H */
diff --git a/src/stars/Default/stars_iact.h b/src/stars/Default/stars_iact.h
index a7633d4eae8a84472fce3f76f2e5d5582b25dd8a..38a64d22553caa370dd96abefea0ba772d423330 100644
--- a/src/stars/Default/stars_iact.h
+++ b/src/stars/Default/stars_iact.h
@@ -35,9 +35,7 @@
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_stars_density(
float r2, const float *dx, float hi, float hj, struct spart *restrict si,
- const struct part *restrict pj, const struct cosmology *restrict cosmo,
- const struct stars_props *restrict stars_properties,
- struct xpart *restrict xp, integertime_t ti_current) {
+ const struct part *restrict pj, const float a, const float H) {
float wi, wi_dx;
@@ -78,10 +76,9 @@ runner_iact_nonsym_stars_density(
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_stars_feedback(
- float r2, const float *dx, float hi, float hj, struct spart *restrict si,
- struct part *restrict pj, const struct cosmology *restrict cosmo,
- const struct stars_props *restrict stars_properties,
- struct xpart *restrict xp, integertime_t ti_current) {
+ float r2, const float *dx, float hi, float hj,
+ const struct spart *restrict si, struct part *restrict pj,
+ const float a, const float H) {
const float mj = hydro_get_mass(pj);
const float rhoj = hydro_get_comoving_density(pj);
diff --git a/src/stars/EAGLE/stars_iact.h b/src/stars/EAGLE/stars_iact.h
index 2cf6c8395596e26bf8f5560ba0446e4e392cb6d5..5c82e7777a9103d7c042405476d83c1f84523f62 100644
--- a/src/stars/EAGLE/stars_iact.h
+++ b/src/stars/EAGLE/stars_iact.h
@@ -23,8 +23,6 @@
/**
* @brief Density interaction between two particles (non-symmetric).
- * MATTHIEU: check with RGB about concerns with comment (in this and other
- * subgrid schemes)
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
@@ -34,8 +32,6 @@
* @param pj Second particle (not updated).
* @param a Current scale factor.
* @param H Current Hubble parameter.
- * @param xp Extra particle data
- * @param ti_current Current integer time value
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_stars_density(
@@ -79,13 +75,10 @@ runner_iact_nonsym_stars_density(
* @param dx Comoving vector separating both particles (si - pj).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
- * @param si First (star) particle.
+ * @param si First (star) particle (not updated).
* @param pj Second (gas) particle.
* @param a Current scale factor.
* @param H Current Hubble parameter.
- * @param xp Extra particle data
- * @param ti_current Current integer time used value for seeding random number
- * generator
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_stars_feedback(
diff --git a/src/swift.h b/src/swift.h
index 0b5eda65e73cd6600eb69c25bc6a3f0a2e6e675f..f61fdf4699f743f9f682f64b72c7348dbeea7985 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -40,6 +40,7 @@
#include "engine.h"
#include "entropy_floor.h"
#include "error.h"
+#include "feedback_properties.h"
#include "gravity.h"
#include "gravity_derivatives.h"
#include "gravity_properties.h"