Re-arranged the star-formation task to split into simpler logical pieces.

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,

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);

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;