diff --git a/src/random.h b/src/random.h
index 6a8f45e69bc5c1b1d771a3fd6b5bab494f034008..480827c96df94d47ce84b67cc55217010f6c1954 100644
--- a/src/random.h
+++ b/src/random.h
@@ -44,7 +44,8 @@ enum random_number_type {
*
* @param id The ID of the particle for which to generate a number.
* @param ti_current The time (on the time-line) for which to generate a number.
- * @prarm type The #random_number_type to generate.
+ * @param type The #random_number_type to generate.
+ * @return a random number in the interval [0, 1.[.
*/
INLINE static double random_unit_interval(const long long int id,
const integertime_t ti_current,
diff --git a/src/runner.c b/src/runner.c
index 1087b5a39750ce8b4ca191d88cab2da41e2efa41..df8a52e2cae3f06fb4e59b1cfbce3cf20160f947 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -475,8 +475,8 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
- const struct star_formation *starform = e->star_formation;
- const struct phys_const *constants = e->physical_constants;
+ const struct star_formation *sf_props = e->star_formation;
+ const struct phys_const *phys_const = e->physical_constants;
const int count = c->hydro.count;
struct part *restrict parts = c->hydro.parts;
struct xpart *restrict xparts = c->hydro.xparts;
@@ -505,35 +505,52 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
struct part *restrict p = &parts[k];
struct xpart *restrict xp = &xparts[k];
+ /* Only work on active particles */
if (part_is_active(p, e)) {
- /* Calculate the time step of the current particle for cosmo and no
- * cosmo*/
- double dt_star;
- if (with_cosmology) {
- const integertime_t ti_step = get_integer_timestep(p->time_bin);
- const integertime_t ti_begin =
- get_integer_time_begin(ti_current - 1, p->time_bin);
+ /* Is this particle star forming? */
+ if (star_formation_is_star_forming(p, xp, sf_props, phys_const, cosmo,
+ hydro_props, us, cooling)) {
- dt_star =
- cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
+ /* Time-step size for this particle */
+ double dt_star;
+ if (with_cosmology) {
+ const integertime_t ti_step = get_integer_timestep(p->time_bin);
+ const integertime_t ti_begin =
+ get_integer_time_begin(ti_current - 1, p->time_bin);
- } else {
- dt_star = get_timestep(p->time_bin, time_base);
- }
+ dt_star =
+ cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
+
+ } else {
+ dt_star = get_timestep(p->time_bin, time_base);
+ }
- if (star_formation_should_convert_to_star(
- e, starform, p, xp, constants, cosmo, hydro_props, us, cooling,
- dt_star, with_cosmology)) {
- /* Convert your particle to a star */
- struct spart *sp = cell_convert_part_to_spart(e, c, p, xp);
+ /* Compute the SF rate of the particle */
+ star_formation_compute_SFR(p, xp, sf_props, phys_const, cosmo,
+ dt_star);
- /* Copy the properties of the gas particle to the star particle */
- star_formation_copy_properties(e, p, xp, sp, starform, constants,
- cosmo, with_cosmology);
- }
- }
- }
+ /* Are we forming a star particle from this SF rate? */
+ if (star_formation_should_convert_to_star(p, xp, sf_props, e,
+ dt_star)) {
+
+ /* Convert the gas particle to a star particle */
+ struct spart *sp = cell_convert_part_to_spart(e, c, p, xp);
+
+ /* Copy the properties of the gas particle to the star particle */
+ star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
+ with_cosmology);
+ }
+
+ } else { /* Are we not star-forming? */
+
+ /* Update the particle to flag it as not star-forming */
+ star_formation_update_part_not_SFR(p, xp, e, sf_props,
+ with_cosmology);
+
+ } /* Not Star-forming? */
+ } /* is active? */
+ } /* Loop over particles */
}
if (timer) TIMER_TOC(timer_do_star_formation);
diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h
index 8b32332865d97845901df60e642bed458ab744f6..934c48eb00a5255cb34c3a5da271825bf3bea062 100644
--- a/src/star_formation/EAGLE/star_formation.h
+++ b/src/star_formation/EAGLE/star_formation.h
@@ -213,8 +213,8 @@ INLINE static float EOS_temperature(const float n_H,
*
*/
INLINE static int star_formation_is_star_forming(
- const struct star_formation* starform, const struct part* restrict p,
- const struct xpart* restrict xp, const struct phys_const* const phys_const,
+ const struct part* restrict p, const struct xpart* restrict xp,
+ const struct star_formation* starform, const struct phys_const* phys_const,
const struct cosmology* cosmo,
const struct hydro_props* restrict hydro_props,
const struct unit_system* restrict us,
@@ -263,89 +263,115 @@ INLINE static int star_formation_is_star_forming(
}
/**
- * @brief Calculates if the gas particle gets converted
+ * @brief Compute the star-formation rate of a given particle and store
+ * it into the #xpart.
*
- * @param e the #engine
- * @param starform the star formation law properties to use.
- * @param p the gas particles.
- * @param xp the additional properties of the gas particles.
+ * @param p #part.
+ * @param xp the #xpart.
+ * @param starform the star formation law properties to use
* @param phys_const the physical constants in internal units.
* @param cosmo the cosmological parameters and properties.
- * @param hydro_props The properties of the hydro scheme.
- * @param us The internal system of units.
- * @param cooling The cooling data struct.
* @param dt_star The time-step of this particle.
- * @param with_cosmology Are we running with cosmology on?
*/
-INLINE static int star_formation_should_convert_to_star(
- const struct engine* e, const struct star_formation* starform,
+INLINE static void star_formation_compute_SFR(
const struct part* restrict p, struct xpart* restrict xp,
- const struct phys_const* const phys_const, const struct cosmology* cosmo,
- const struct hydro_props* restrict hydro_props,
- const struct unit_system* restrict us,
- const struct cooling_function_data* restrict cooling, const double dt_star,
- const int with_cosmology) {
+ const struct star_formation* starform, const struct phys_const* phys_const,
+ const struct cosmology* cosmo, const double dt_star) {
/* Abort early if time-step size is 0 */
if (dt_star == 0.f) {
- return 0;
+
+ xp->sf_data.SFR = 0.f;
+ return;
}
- if (star_formation_is_star_forming(starform, p, xp, phys_const, cosmo,
- hydro_props, us, cooling)) {
+ /* Hydrogen number density of this particle */
+ const double physical_density = hydro_get_physical_density(p, cosmo);
+ const double X_H = p->chemistry_data.smoothed_metal_mass_fraction[0];
+ const double n_H = physical_density * X_H / phys_const->const_proton_mass;
- /* Hydrogen number density of this particle */
- const double physical_density = hydro_get_physical_density(p, cosmo);
- const double X_H = p->chemistry_data.smoothed_metal_mass_fraction[0];
- const double n_H = physical_density * X_H / phys_const->const_proton_mass;
+ /* Are we above the threshold for automatic star formation? */
+ if (physical_density >
+ starform->max_gas_density * phys_const->const_proton_mass) {
- /* Are we above the threshold for automatic star formation? */
- if (physical_density >
- starform->max_gas_density * phys_const->const_proton_mass) {
- return 1;
- }
+ xp->sf_data.SFR = p->mass / dt_star;
+ return;
+ }
- /* Pressure on the EOS for this particle */
- const double pressure = EOS_pressure(n_H, starform);
+ /* Pressure on the effective EOS for this particle */
+ const double pressure = EOS_pressure(n_H, starform);
- /* Calculate the specific star formation rate */
- double SFRpergasmass;
- if (hydro_get_physical_density(p, cosmo) <
- starform->KS_high_den_thresh * phys_const->const_proton_mass) {
+ /* Calculate the specific star formation rate */
+ double SFRpergasmass;
+ if (hydro_get_physical_density(p, cosmo) <
+ starform->KS_high_den_thresh * phys_const->const_proton_mass) {
- SFRpergasmass =
- starform->SF_normalization * pow(pressure, starform->SF_power_law);
+ SFRpergasmass =
+ starform->SF_normalization * pow(pressure, starform->SF_power_law);
- } else {
+ } else {
- SFRpergasmass = starform->SF_high_den_normalization *
- pow(pressure, starform->SF_high_den_power_law);
- }
+ SFRpergasmass = starform->SF_high_den_normalization *
+ pow(pressure, starform->SF_high_den_power_law);
+ }
- /* Store the SFR */
- xp->sf_data.SFR = SFRpergasmass * p->mass;
+ /* Store the SFR */
+ xp->sf_data.SFR = SFRpergasmass * p->mass;
+}
- /* Calculate the propability of forming a star */
- const double prop = SFRpergasmass * dt_star;
+/**
+ * @brief Decides whether a particle should be converted into a
+ * star or not.
+ *
+ * Equation 21 of Schaye & Dalla Vecchia 2008.
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param starform The properties of the star formation model.
+ * @param e The #engine (for random numbers).
+ * @param dt_star The time-step of this particle
+ * @return 1 if a conversion should be done, 0 otherwise.
+ */
+INLINE static int star_formation_should_convert_to_star(
+ const struct part* p, const struct xpart* xp,
+ const struct star_formation* starform, const struct engine* e,
+ const double dt_star) {
- /* Get a random number between 0 and 1 for star formation */
- const double random_number = random_unit_interval(
- p->id, e->ti_current, random_number_star_formation);
+ /* Calculate the propability of forming a star */
+ const double prob = xp->sf_data.SFR * dt_star / p->mass;
- /* Have we been lucky and need to form a star? */
- return (prop > random_number);
- }
+ /* Get a unique random number between 0 and 1 for star formation */
+ const double random_number =
+ random_unit_interval(p->id, e->ti_current, random_number_star_formation);
+
+ /* Have we been lucky and need to form a star? */
+ return (prob > random_number);
+}
+
+/**
+ * @brief Update the SF properties of a particle that is not star forming.
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param e The #engine.
+ * @param starform The properties of the star formation model.
+ * @param with_cosmology Are we running with cosmology switched on?
+ */
+INLINE static void star_formation_update_part_not_SFR(
+ struct part* p, struct xpart* xp, const struct engine* e,
+ const struct star_formation* starform, const int with_cosmology) {
/* Check if it is the first time steps after star formation */
if (xp->sf_data.SFR > 0.f) {
+
+ /* Record the current time as an indicator of when this particle was last
+ star-forming. */
if (with_cosmology) {
- xp->sf_data.SFR = -cosmo->a;
+ xp->sf_data.SFR = -e->cosmology->a;
} else {
xp->sf_data.SFR = -e->time;
}
}
-
- return 0;
}
/**
@@ -362,10 +388,9 @@ INLINE static int star_formation_should_convert_to_star(
* @param with_cosmology if we run with cosmology.
*/
INLINE static void star_formation_copy_properties(
- const struct engine* e, const struct part* p, const struct xpart* xp,
- struct spart* sp, const struct star_formation* starform,
- const struct phys_const* const phys_const, const struct cosmology* cosmo,
- const int with_cosmology) {
+ const struct part* p, const struct xpart* xp, struct spart* sp,
+ const struct engine* e, const struct star_formation* starform,
+ const struct cosmology* cosmo, const int with_cosmology) {
/* Store the current mass */
sp->mass = p->mass;