Revert "Merge branch 'gear_star_formation' into 'master'"

This reverts merge request !839

examples/parameter_example.yml

@@ -338,11 +338,6 @@ 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, xp, e->star_formation);
+        star_formation_init_part(p, 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,8 +122,7 @@ 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;
-
-  star_formation_logger_add_to_accumulator(&e->sfh, &grp1->sfh);
+  e->sfh = grp1->sfh;
 }
 
 /**

src/cooling/grackle/cooling.h

@@ -697,17 +697,6 @@ __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,
@@ -715,30 +704,9 @@ 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;
 
-  /* 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;
+  error("This function needs implementing!!");
+  return 0.;
 }
 
 /**

src/engine.c

@@ -5062,11 +5062,6 @@ 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_accumulator sfh;
+  struct star_formation_history sfh;
 
   /* Properties of the previous step */
   int step_props;

src/hydro/Gadget2/hydro_part.h

@@ -155,9 +155,6 @@ 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, xp, star_formation);
+            star_formation_init_part(p, 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,6 +884,7 @@ 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], &xparts[k], e->star_formation);
+    star_formation_init_part(&parts[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,7 +24,6 @@
 #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,7 +691,6 @@ 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, struct xpart* restrict xp,
-    const struct star_formation* data) {}
+    struct part* restrict p, const struct star_formation* data) {}
 
 #endif /* SWIFT_EAGLE_STAR_FORMATION_H */

src/star_formation/EAGLE/star_formation_logger.h

@@ -87,28 +87,6 @@ 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
  *
@@ -127,24 +105,6 @@ 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
  *
@@ -157,7 +117,7 @@ INLINE static void star_formation_logger_accumulator_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_accumulator sf, const int step) {
+    const struct star_formation_history 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,21 +34,4 @@ 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,6 +29,4 @@ 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.
- * Coypright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
- *               2019 Fabien Jeanquartier (fabien.jeanquartier@epfl.ch)
+ * Copyright (c) 2018 Folkert Nobels (nobels@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
@@ -21,24 +20,20 @@
 #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.
  *
- * Use the star formation criterion given by eq. 3 in Revaz & Jablonka 2018.
+ * No star formation should occur, so return 0.
  *
  * @param starform the star formation law properties to use.
  * @param p the gas particles.
@@ -51,7 +46,7 @@
  *
  */
 INLINE static int star_formation_is_star_forming(
-    struct part* restrict p, struct xpart* restrict xp,
+    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,
@@ -59,43 +54,14 @@ INLINE static int star_formation_is_star_forming(
     const struct cooling_function_data* restrict cooling,
     const struct entropy_floor_properties* restrict entropy_floor) {
 
-  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;
+  return 0;
 }
 
 /**
- * @brief Compute the star-formation rate of a given particle.
+ * @brief Compute the star-formation rate of a given particle and store
+ * it into the #xpart.
  *
- * Nothing to do here. Everything is done in
- * #star_formation_should_convert_to_star.
+ * Nothing to do here.
  *
  * @param p #part.
  * @param xp the #xpart.
@@ -105,15 +71,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(
-    struct part* restrict p, struct xpart* restrict xp,
+    const 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.
-
- * Compute the star formation rate from eq. 4 in Revaz & Jablonka 2012.
+ *
+ * No SF should occur, so return 0.
  *
  * @param p The #part.
  * @param xp The #xpart.
@@ -123,38 +89,18 @@ 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(
-    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);
+    const struct part* p, const struct xpart* xp,
+    const struct star_formation* starform, const struct engine* e,
+    const double dt_star) {
 
-  /* Can we form a star? */
-  return random_number < prob;
+  return 0;
 }
 
 /**
  * @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.
@@ -169,6 +115,8 @@ 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.
@@ -185,33 +133,21 @@ 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) {
-
-  /* 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);
+    const struct cooling_function_data* restrict cooling) {}
 
-  /* Store the birth temperature*/
-  sp->birth.temperature = cooling_get_temperature(phys_const, hydro_props, us,
-                                                  cosmo, cooling, p, xp);
-}
+/**
+ * @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) {}
 
 /**
  * @brief Prints the used parameters of the star formation law
@@ -220,6 +156,7 @@ INLINE static void star_formation_copy_properties(
  */
 INLINE static void starformation_print_backend(
     const struct star_formation* starform) {
+
   message("Star formation law is 'GEAR'");
 }
 
@@ -227,22 +164,12 @@ INLINE static void starformation_print_backend(
  * @brief Finishes the density calculation.
  *
  * @param p The particle to act upon
- * @param xp The extended particle data to act upon
- * @param sf The global star_formation information.
+ * @param cd 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* 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;
-}
+    struct part* restrict p, const struct star_formation* cd,
+    const struct cosmology* cosmo) {}
 
 /**
  * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
@@ -256,46 +183,39 @@ __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) {
-
-  // TODO move into pressure floor
-  /* If part has 0 neighbours, the estimation of turbulence is 0 */
-  p->sf_data.sigma2 = 0.f;
-}
+                                      const struct cosmology* cosmo) {}
 
 /**
  * @brief Sets the star_formation properties of the (x-)particles to a valid
  * start state.
  *
- * @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.
+ * Nothing to do here.
+ *
  * @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,
-                               struct part* restrict p,
-                               struct xpart* restrict xp) {
+                               const struct part* restrict p,
+                               struct xpart* restrict xp) {}
 
-  /* Nothing special here */
-  star_formation_init_part(p, xp, data);
-}
+/**
+ * @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) {}
 
 #endif /* SWIFT_GEAR_STAR_FORMATION_H */

src/star_formation/GEAR/star_formation_iact.h

 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
- *               2019 Fabien Jeanquartier (fabien.jeanquartier@epfl.ch)
+ * Copyright (c) 2019 Loic Hausammann (loic.hausammann@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
@@ -29,8 +28,6 @@
  * @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.
@@ -42,28 +39,7 @@
  */
 __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) {
-
-  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);
-}
+    struct part *restrict pj, float a, float H) {}
 
 /**