Added a more general function to compute reproducible pseudo-random numbers.

src/Makefile.am

@@ -50,7 +50,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h \
     mesh_gravity.h cbrt.h exp10.h velociraptor_interface.h swift_velociraptor_part.h outputlist.h \
     logger_io.h tracers_io.h tracers.h tracers_struct.h star_formation_io.h \
-    star_formation_struct.h star_formation.h velociraptor_struct.h velociraptor_io.h
+    star_formation_struct.h star_formation.h velociraptor_struct.h velociraptor_io.h \
+    random.h
 
 # source files for EAGLE cooling
 EAGLE_COOLING_SOURCES =

src/random.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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_RANDOM_H
+#define SWIFT_RANDOM_H
+
+/* COde configuration */
+#include "../config.h"
+
+/* Standard header */
+#include <stdlib.h>
+
+/**
+ * @brief The categories of random number generated.
+ */
+enum random_number_type {
+  random_number_star_formation = 0,
+  random_number_stellar_feedback,
+  random_number_stellar_enrichment,
+  random_number_BH_feedback
+};
+
+/**
+ * @brief Returns a pseudo-random number in the range [0, 1[.
+ *
+ * We generate numbers that are always reproducible for a given particle ID and
+ * simulation time. If more than one number per time-step per particle is
+ * needed, additional randomness can be obtained by using the type argument.
+ *
+ * @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.
+ */
+INLINE static double random_unit_interval(const long long int id,
+                                          const integertime_t ti_current,
+                                          const enum random_number_type type) {
+
+  /* Range used for the seeds. Best if prime */
+  static const long long seed_range = 131071LL;
+  static const double RAND_MAX_inv = 1. / ((double)RAND_MAX);
+
+  /* Calculate the seed */
+  unsigned int seed = (id + ti_current + (long long)type) % seed_range;
+
+  /* Generate a random number between 0 and 1. */
+  return rand_r(&seed) / RAND_MAX_inv;
+}
+
+#endif /* SWIFT_RANDOM_H */

src/star_formation/EAGLE/star_formation.h

@@ -29,6 +29,7 @@
 #include "parser.h"
 #include "part.h"
 #include "physical_constants.h"
+#include "random.h"
 #include "stars.h"
 #include "units.h"
 
@@ -296,23 +297,25 @@ INLINE static int star_formation_convert_to_star(
     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) {
+      return 1;
+    }
+
     /* Pressure on the 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 star formation rate */
+
       SFRpergasmass =
           starform->SF_normalization * pow(pressure, starform->SF_power_law);
-    } else if (hydro_get_physical_density(p, cosmo) >
-               starform->max_gas_density * phys_const->const_proton_mass) {
-      /* We give the star formation tracers values of the mass and -1 for sSFR
-       * to be able to trace them, we don't want random variables there*/
-      xp->sf_data.SFR = p->mass;
-      xp->sf_data.sSFR = -1.f;
-      return 1.f;
+
     } else {
+
       SFRpergasmass = starform->SF_high_den_normalization *
                       pow(pressure, starform->SF_high_den_power_law);
     }
@@ -320,18 +323,16 @@ INLINE static int star_formation_convert_to_star(
     /* Store the SFR */
     xp->sf_data.SFR = SFRpergasmass * p->mass;
     xp->sf_data.sSFR = SFRpergasmass;
+
     /* Calculate the propability of forming a star */
     const double prop = SFRpergasmass * dt_star;
 
-    /* Calculate the seed */
-    unsigned int seed = (p->id + e->ti_current) % 8191;
-
-    /* Generate a random number between 0 and 1. */
-    const double randomnumber =
-        rand_r(&seed);  // MATTHIEU: * starform->inv_RAND_MAX;
+    /* 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 if we form a star */
-    return (prop > randomnumber);
+    /* Have we been lucky and need to form a star? */
+    return (prop > random_number);
   }
 
   /* Check if it is the first time steps after star formation */