GEAR: add star formation

examples/parameter_example.yml

@@ -338,6 +338,11 @@ EAGLEChemistry:
 
 # Parameters related to star formation models  -----------------------------------------------
 
+# GEAR star formation model (Revaz and Jablonka 2018)
+GEARStarFormation:
+  star_formation_efficiency: 0.01   # star formation efficiency (c_*)
+  maximal_temperature:  3e4         # Upper limit to the temperature of a star forming particle
+
 # EAGLE star formation model (Schaye and Dalla Vecchia 2008)
 EAGLEStarFormation:
   EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.

src/cell.c

@@ -4413,7 +4413,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
       if (part_is_active(p, e)) {
         hydro_init_part(p, &e->s->hs);
         chemistry_init_part(p, e->chemistry);
-        star_formation_init_part(p, e->star_formation);
+        star_formation_init_part(p, xp, e->star_formation);
         tracers_after_init(p, xp, e->internal_units, e->physical_constants,
                            with_cosmology, e->cosmology, e->hydro_properties,
                            e->cooling_func, e->time);

src/collectgroup.c

@@ -122,7 +122,8 @@ void collectgroup1_apply(const struct collectgroup1 *grp1, struct engine *e) {
   e->total_nr_cells = grp1->total_nr_cells;
   e->total_nr_tasks = grp1->total_nr_tasks;
   e->tasks_per_cell_max = grp1->tasks_per_cell_max;
-  e->sfh = grp1->sfh;
+
+  star_formation_logger_add_to_accumulator(&e->sfh, &grp1->sfh);
 }
 
 /**

src/cooling/grackle/cooling.h

@@ -697,6 +697,17 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
   xp->cooling_data.radiated_energy -= hydro_get_mass(p) * cooling_du_dt * dt;
 }
 
+/**
+ * @brief Compute the temperature of a #part based on the cooling function.
+ *
+ * @param phys_const #phys_const data structure.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param us The internal system of units.
+ * @param cosmo #cosmology data structure.
+ * @param cooling #cooling_function_data struct.
+ * @param p #part data.
+ * @param xp Pointer to the #xpart data.
+ */
 static INLINE float cooling_get_temperature(
     const struct phys_const* restrict phys_const,
     const struct hydro_props* restrict hydro_props,
@@ -704,9 +715,30 @@ static INLINE float cooling_get_temperature(
     const struct cosmology* restrict cosmo,
     const struct cooling_function_data* restrict cooling,
     const struct part* restrict p, const struct xpart* restrict xp) {
+  // TODO use the grackle library
+
+  /* Physical constants */
+  const double m_H = phys_const->const_proton_mass;
+  const double k_B = phys_const->const_boltzmann_k;
 
-  error("This function needs implementing!!");
-  return 0.;
+  /* Gas properties */
+  const double T_transition = hydro_props->hydrogen_ionization_temperature;
+  const double mu_neutral = hydro_props->mu_neutral;
+  const double mu_ionised = hydro_props->mu_ionised;
+
+  /* Particle temperature */
+  const double u = hydro_get_physical_internal_energy(p, xp, cosmo);
+
+  /* Temperature over mean molecular weight */
+  const double T_over_mu = hydro_gamma_minus_one * u * m_H / k_B;
+
+  /* Are we above or below the HII -> HI transition? */
+  if (T_over_mu > (T_transition + 1.) / mu_ionised)
+    return T_over_mu * mu_ionised;
+  else if (T_over_mu < (T_transition - 1.) / mu_neutral)
+    return T_over_mu * mu_neutral;
+  else
+    return T_transition;
 }
 
 /**

src/engine.c

@@ -5062,6 +5062,11 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
         parser_get_opt_param_double(params, "FOF:delta_time", -1.);
   }
 
+  /* Initialize the star formation history structure */
+  if (e->policy & engine_policy_star_formation) {
+    star_formation_logger_accumulator_init(&e->sfh);
+  }
+
   engine_init_output_lists(e, params);
 }
 

src/engine.h

@@ -237,7 +237,7 @@ struct engine {
   long long b_updates_since_rebuild;
 
   /* Star formation logger information */
-  struct star_formation_history sfh;
+  struct star_formation_history_accumulator sfh;
 
   /* Properties of the previous step */
   int step_props;

src/hydro/Gadget2/hydro_part.h

@@ -155,6 +155,9 @@ struct part {
   /*! Black holes information (e.g. swallowing ID) */
   struct black_holes_part_data black_holes_data;
 
+  /* Additional data used by the star formation */
+  struct star_formation_part_data sf_data;
+
   /* Time-step length */
   timebin_t time_bin;
 

src/runner.c

@@ -2272,7 +2272,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
             /* Re-initialise everything */
             hydro_init_part(p, hs);
             chemistry_init_part(p, chemistry);
-            star_formation_init_part(p, star_formation);
+            star_formation_init_part(p, xp, star_formation);
             tracers_after_init(p, xp, e->internal_units, e->physical_constants,
                                with_cosmology, e->cosmology,
                                e->hydro_properties, e->cooling_func, e->time);

src/runner_doiact.h

@@ -884,7 +884,6 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci,
 
   const int count_i = ci->hydro.count;
   struct part *restrict parts_j = ci->hydro.parts;
-
   /* Loop over the parts in ci. */
   for (int pid = 0; pid < count; pid++) {
 

src/space.c

@@ -4365,7 +4365,7 @@ void space_init_parts_mapper(void *restrict map_data, int count,
   for (int k = 0; k < count; k++) {
     hydro_init_part(&parts[k], hs);
     chemistry_init_part(&parts[k], e->chemistry);
-    star_formation_init_part(&parts[k], e->star_formation);
+    star_formation_init_part(&parts[k], &xparts[k], e->star_formation);
     tracers_after_init(&parts[k], &xparts[k], e->internal_units,
                        e->physical_constants, with_cosmology, e->cosmology,
                        e->hydro_properties, e->cooling_func, e->time);

src/star_formation.c

@@ -24,6 +24,7 @@
 #include "part.h"
 #include "restart.h"
 #include "star_formation.h"
+#include "star_formation_io.h"
 #include "units.h"
 
 /**

src/star_formation/EAGLE/star_formation.h

@@ -691,6 +691,7 @@ star_formation_first_init_part(const struct phys_const* restrict phys_const,
  * @param data The global star_formation information.
  */
 __attribute__((always_inline)) INLINE static void star_formation_init_part(
-    struct part* restrict p, const struct star_formation* data) {}
+    struct part* restrict p, struct xpart* restrict xp,
+    const struct star_formation* data) {}
 
 #endif /* SWIFT_EAGLE_STAR_FORMATION_H */

src/star_formation/EAGLE/star_formation_logger.h

@@ -87,6 +87,28 @@ INLINE static void star_formation_logger_add(
   sf_update->SFR_inactive += sf_add->SFR_inactive;
 }
 
+/**
+ * @brief add a star formation history struct to the engine star formation
+ * history accumulator struct
+ *
+ * @param sf_add the star formation accumulator struct which we want to add to
+ * the star formation history
+ * @param sf_update the star formation structure which we want to update
+ */
+INLINE static void star_formation_logger_add_to_accumulator(
+    struct star_formation_history_accumulator *sf_update,
+    const struct star_formation_history *sf_add) {
+
+  /* Update the SFH structure */
+  sf_update->new_stellar_mass += sf_add->new_stellar_mass;
+
+  sf_update->SFR_active += sf_add->SFR_active;
+
+  sf_update->SFRdt_active += sf_add->SFRdt_active;
+
+  sf_update->SFR_inactive += sf_add->SFR_inactive;
+}
+
 /**
  * @brief Initialize the star formation history structure in the #engine
  *
@@ -105,6 +127,24 @@ INLINE static void star_formation_logger_init(
   sfh->SFR_inactive = 0.f;
 }
 
+/**
+ * @brief Initialize the star formation history structure in the #engine
+ *
+ * @param sfh The pointer to the star formation history structure
+ */
+INLINE static void star_formation_logger_accumulator_init(
+    struct star_formation_history_accumulator *sfh) {
+
+  /* Initialize the collecting SFH structure to zero */
+  sfh->new_stellar_mass = 0.f;
+
+  sfh->SFR_active = 0.f;
+
+  sfh->SFRdt_active = 0.f;
+
+  sfh->SFR_inactive = 0.f;
+}
+
 /**
  * @brief Write the final SFH to a file
  *
@@ -117,7 +157,7 @@ INLINE static void star_formation_logger_init(
  */
 INLINE static void star_formation_logger_write_to_log_file(
     FILE *fp, const double time, const double a, const double z,
-    const struct star_formation_history sf, const int step) {
+    const struct star_formation_history_accumulator sf, const int step) {
 
   /* Calculate the total SFR */
   const float totalSFR = sf.SFR_active + sf.SFR_inactive;

src/star_formation/EAGLE/star_formation_logger_struct.h

@@ -34,4 +34,21 @@ struct star_formation_history {
   float SFRdt_active;
 };
 
+/* Starformation history struct for the engine.
+ Allows to integrate in time some values.
+ Nothing to do in EAGLE => copy of star_formation_history */
+struct star_formation_history_accumulator {
+  /*! Total new stellar mass */
+  float new_stellar_mass;
+
+  /*! SFR of all particles */
+  float SFR_inactive;
+
+  /*! SFR of active particles */
+  float SFR_active;
+
+  /*! SFR*dt of active particles */
+  float SFRdt_active;
+};
+
 #endif /* SWIFT_EAGLE_STAR_FORMATION_LOGGER_STRUCT_H */

src/star_formation/EAGLE/star_formation_struct.h

@@ -29,4 +29,6 @@ struct star_formation_xpart_data {
   float SFR;
 };
 
+struct star_formation_part_data {};
+
 #endif /* SWIFT_EAGLE_STAR_FORMATION_STRUCT_H */

src/star_formation/GEAR/star_formation.h

 /*******************************************************************************
  * This file is part of SWIFT.
- * Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
+ * Coypright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+ *               2019 Fabien Jeanquartier (fabien.jeanquartier@epfl.ch)
  *
  * 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
@@ -20,20 +21,24 @@
 #define SWIFT_GEAR_STAR_FORMATION_H
 
 /* Local includes */
+#include "cooling.h"
 #include "cosmology.h"
+#include "engine.h"
 #include "entropy_floor.h"
 #include "error.h"
 #include "hydro_properties.h"
 #include "parser.h"
 #include "part.h"
 #include "physical_constants.h"
+#include "random.h"
+#include "star_formation_struct.h"
 #include "units.h"
 
 /**
  * @brief Calculate if the gas has the potential of becoming
  * a star.
  *
- * No star formation should occur, so return 0.
+ * Use the star formation criterion given by eq. 3 in Revaz & Jablonka 2018.
  *
  * @param starform the star formation law properties to use.
  * @param p the gas particles.
@@ -46,7 +51,7 @@
  *
  */
 INLINE static int star_formation_is_star_forming(
-    const struct part* restrict p, const struct xpart* restrict xp,
+    struct part* restrict p, 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,
@@ -54,14 +59,43 @@ INLINE static int star_formation_is_star_forming(
     const struct cooling_function_data* restrict cooling,
     const struct entropy_floor_properties* restrict entropy_floor) {
 
-  return 0;
+  const float temperature = cooling_get_temperature(phys_const, hydro_props, us,
+                                                    cosmo, cooling, p, xp);
+
+  const float temperature_max = starform->maximal_temperature;
+
+  /* Check the temperature criterion */
+  if (temperature > temperature_max) {
+    return 0;
+  }
+
+  /* Get the required variables */
+  const float sigma2 = p->sf_data.sigma2;
+  const float n_jeans_2_3 = starform->n_jeans_2_3;
+
+  const float h = p->h;
+  const float density = hydro_get_physical_density(p, cosmo);
+
+  // TODO use GRACKLE */
+  const float mu = hydro_props->mu_neutral;
+
+  /* Compute the density criterion */
+  const float coef =
+      M_PI_4 / (phys_const->const_newton_G * n_jeans_2_3 * h * h);
+  const float density_criterion =
+      coef * (hydro_gamma * phys_const->const_boltzmann_k * temperature /
+                  (mu * phys_const->const_proton_mass) +
+              sigma2);
+
+  /* Check the density criterion */
+  return density > density_criterion;
 }
 
 /**
- * @brief Compute the star-formation rate of a given particle and store
- * it into the #xpart.
+ * @brief Compute the star-formation rate of a given particle.
  *
- * Nothing to do here.
+ * Nothing to do here. Everything is done in
+ * #star_formation_should_convert_to_star.
  *
  * @param p #part.
  * @param xp the #xpart.
@@ -71,15 +105,15 @@ INLINE static int star_formation_is_star_forming(
  * @param dt_star The time-step of this particle.
  */
 INLINE static void star_formation_compute_SFR(
-    const struct part* restrict p, struct xpart* restrict xp,
+    struct part* restrict p, struct xpart* restrict xp,
     const struct star_formation* starform, const struct phys_const* phys_const,
     const struct cosmology* cosmo, const double dt_star) {}
 
 /**
  * @brief Decides whether a particle should be converted into a
  * star or not.
- *
- * No SF should occur, so return 0.
+
+ * Compute the star formation rate from eq. 4 in Revaz & Jablonka 2012.
  *
  * @param p The #part.
  * @param xp The #xpart.
@@ -89,18 +123,38 @@ INLINE static void star_formation_compute_SFR(
  * @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) {
+    struct part* p, struct xpart* xp, const struct star_formation* starform,
+    const struct engine* e, const double dt_star) {
+
+  const struct phys_const* phys_const = e->physical_constants;
+  const struct cosmology* cosmo = e->cosmology;
+
+  /* Check that we are running a full time step */
+  if (dt_star == 0.) {
+    return 0;
+  }
+
+  /* Get a few variables */
+  const float G = phys_const->const_newton_G;
+  const float density = hydro_get_physical_density(p, cosmo);
+
+  /* Compute the probability */
+  const float inv_free_fall_time =
+      sqrtf(density * 32.f * G * 0.33333333f * M_1_PI);
+  const float prob = 1.f - exp(-starform->star_formation_efficiency *
+                               inv_free_fall_time * dt_star);
+
+  /* Roll the dice... */
+  const float random_number =
+      random_unit_interval(p->id, e->ti_current, random_number_star_formation);
 
-  return 0;
+  /* Can we form a star? */
+  return random_number < prob;
 }
 
 /**
  * @brief Update the SF properties of a particle that is not star forming.
  *
- * Nothing to do here.
- *
  * @param p The #part.
  * @param xp The #xpart.
  * @param e The #engine.
@@ -115,8 +169,6 @@ INLINE static void star_formation_update_part_not_SFR(
  * @brief Copies the properties of the gas particle over to the
  * star particle.
  *
- * Nothing to do here.
- *
  * @param e The #engine
  * @param p the gas particles.
  * @param xp the additional properties of the gas particles.
@@ -133,21 +185,33 @@ INLINE static void star_formation_copy_properties(
     const struct phys_const* phys_const,
     const struct hydro_props* restrict hydro_props,
     const struct unit_system* restrict us,
-    const struct cooling_function_data* restrict cooling) {}
+    const struct cooling_function_data* restrict cooling) {
 
-/**
- * @brief initialization of the star formation law
- *
- * @param parameter_file The parsed parameter file
- * @param phys_const Physical constants in internal units
- * @param us The current internal system of units
- * @param starform the star formation law properties to initialize
- *
- */
-INLINE static void starformation_init_backend(
-    struct swift_params* parameter_file, const struct phys_const* phys_const,
-    const struct unit_system* us, const struct hydro_props* hydro_props,
-    const struct star_formation* starform) {}
+  /* Store the current mass */
+  sp->mass = hydro_get_mass(p);
+  sp->birth.mass = sp->mass;
+
+  /* Store either the birth_scale_factor or birth_time depending  */
+  if (with_cosmology) {
+    sp->birth_scale_factor = cosmo->a;
+  } else {
+    sp->birth_time = e->time;
+  }
+
+  // TODO copy only metals
+  /* Store the chemistry struct in the star particle */
+  sp->chemistry_data = p->chemistry_data;
+
+  /* Store the tracers data */
+  sp->tracers_data = xp->tracers_data;
+
+  /* Store the birth density in the star particle */
+  sp->birth.density = hydro_get_physical_density(p, cosmo);
+
+  /* Store the birth temperature*/
+  sp->birth.temperature = cooling_get_temperature(phys_const, hydro_props, us,
+                                                  cosmo, cooling, p, xp);
+}
 
 /**
  * @brief Prints the used parameters of the star formation law
@@ -156,7 +220,6 @@ INLINE static void starformation_init_backend(
  */
 INLINE static void starformation_print_backend(
     const struct star_formation* starform) {
-
   message("Star formation law is 'GEAR'");
 }
 
@@ -164,12 +227,22 @@ INLINE static void starformation_print_backend(
  * @brief Finishes the density calculation.
  *
  * @param p The particle to act upon
- * @param cd The global star_formation information.
+ * @param xp The extended particle data to act upon
+ * @param sf The global star_formation information.
  * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void star_formation_end_density(
-    struct part* restrict p, const struct star_formation* cd,
-    const struct cosmology* cosmo) {}
+    struct part* restrict p, const struct star_formation* sf,
+    const struct cosmology* cosmo) {
+
+  // TODO move into pressure floor
+  /* To finish the turbulence estimation we devide by the density */
+  p->sf_data.sigma2 /=
+      pow_dimension(p->h) * hydro_get_physical_density(p, cosmo);
+
+  /* Add the cosmological factor */
+  p->sf_data.sigma2 *= cosmo->a * cosmo->a;
+}
 
 /**
  * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
@@ -183,39 +256,46 @@ __attribute__((always_inline)) INLINE static void
 star_formation_part_has_no_neighbours(struct part* restrict p,
                                       struct xpart* restrict xp,
                                       const struct star_formation* cd,
-                                      const struct cosmology* cosmo) {}
+                                      const struct cosmology* cosmo) {
+
+  // TODO move into pressure floor
+  /* If part has 0 neighbours, the estimation of turbulence is 0 */
+  p->sf_data.sigma2 = 0.f;
+}
 
 /**
  * @brief Sets the star_formation properties of the (x-)particles to a valid
  * start state.
  *
- * Nothing to do here.
- *
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to extended particle data
+ * @param data The global star_formation information.
+ */
+__attribute__((always_inline)) INLINE static void star_formation_init_part(
+    struct part* restrict p, struct xpart* restrict xp,
+    const struct star_formation* data) {
+  p->sf_data.sigma2 = 0.f;
+}
+
+/**
+ * @brief Sets the star_formation properties of the (x-)particles to a valid
+ * start state.
  * @param phys_const The physical constant in internal units.
  * @param us The unit system.
  * @param cosmo The current cosmological model.
  * @param data The global star_formation information used for this run.
  * @param p Pointer to the particle data.
- * @param xp Pointer to the extended particle data.
  */
 __attribute__((always_inline)) INLINE static void
 star_formation_first_init_part(const struct phys_const* restrict phys_const,
                                const struct unit_system* restrict us,
                                const struct cosmology* restrict cosmo,
                                const struct star_formation* data,
-                               const struct part* restrict p,
-                               struct xpart* restrict xp) {}
+                               struct part* restrict p,
+                               struct xpart* restrict xp) {
 
-/**
- * @brief Sets the star_formation properties of the (x-)particles to a valid
- * start state.
- *
- * Nothing to do here.
- *
- * @param p Pointer to the particle data.
- * @param data The global star_formation information.
- */
-__attribute__((always_inline)) INLINE static void star_formation_init_part(
-    struct part* restrict p, const struct star_formation* data) {}
+  /* Nothing special here */
+  star_formation_init_part(p, xp, data);
+}
 
 #endif /* SWIFT_GEAR_STAR_FORMATION_H */

src/star_formation/GEAR/star_formation_iact.h

 /*******************************************************************************
  * This file is part of SWIFT.
- * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+ * Coypright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+ *               2019 Fabien Jeanquartier (fabien.jeanquartier@epfl.ch)
  *
  * 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
@@ -28,6 +29,8 @@
  * @brief do star_formation computation after the runner_iact_density (symmetric
  * version)
  *
+ * Compute the velocity dispersion follow eq. 2 in Revaz & Jablonka 2018.
+ *
  * @param r2 Comoving square distance between the two particles.
  * @param dx Comoving vector separating both particles (pi - pj).
  * @param hi Comoving smoothing-length of particle i.
@@ -39,7 +42,28 @@
  */
 __attribute__((always_inline)) INLINE static void runner_iact_star_formation(
     float r2, const float *dx, float hi, float hj, struct part *restrict pi,
-    struct part *restrict pj, float a, float H) {}
+    struct part *restrict pj, float a, float H) {
+
+  float wi;
+  float wj;
+  /* Evaluation of the SPH kernel */
+  kernel_eval(sqrt(r2) / hi, &wi);
+  kernel_eval(sqrt(r2) / hj, &wj);
+
+  /* Delta v */
+  float dv[3] = {pi->v[0] - pj->v[0], pi->v[1] - pj->v[1], pi->v[2] - pj->v[2]};
+
+  /* Norms at power 2 */
+  const float norm_v2 = dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2];
+  const float norm_x2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+  /* Compute the velocity dispersion */
+  const float sigma2 = norm_v2 + H * norm_x2;
+
+  /* Compute the velocity dispersion */
+  pi->sf_data.sigma2 += sigma2 * wi * hydro_get_mass(pj);
+  pj->sf_data.sigma2 += sigma2 * wj * hydro_get_mass(pi);
+}
 
 /**