Gear star formation split

examples/main.c

@@ -1204,6 +1204,11 @@ int main(int argc, char *argv[]) {
       fflush(stdout);
     }
 
+    /* Compute some stats for the star formation */
+    if (with_star_formation) {
+      star_formation_first_init_stats(&starform, &e);
+    }
+
     /* Get some info to the user. */
     if (myrank == 0) {
       const long long N_DM = N_total[swift_type_dark_matter] +

src/cell.c

@@ -5480,10 +5480,44 @@ void cell_recursively_shift_sparts(struct cell *c,
   }
 }
 
+/**
+ * @brief Recursively update the pointer and counter for #gpart after the
+ * addition of a new particle.
+ *
+ * @param c The cell we are working on.
+ * @param progeny_list The list of the progeny index at each level for the
+ * leaf-cell where the particle was added.
+ * @param main_branch Are we in a cell directly above the leaf where the new
+ * particle was added?
+ */
+void cell_recursively_shift_gparts(struct cell *c,
+                                   const int progeny_list[space_cell_maxdepth],
+                                   const int main_branch) {
+  if (c->split) {
+    /* No need to recurse in progenies located before the insestion point */
+    const int first_progeny = main_branch ? progeny_list[(int)c->depth] : 0;
+
+    for (int k = first_progeny; k < 8; ++k) {
+      if (c->progeny[k] != NULL)
+        cell_recursively_shift_gparts(c->progeny[k], progeny_list,
+                                      main_branch && (k == first_progeny));
+    }
+  }
+
+  /* When directly above the leaf with the new particle: increase the particle
+   * count */
+  /* When after the leaf with the new particle: shift by one position */
+  if (main_branch) {
+    c->grav.count++;
+  } else {
+    c->grav.parts++;
+  }
+}
+
 /**
  * @brief "Add" a #spart in a given #cell.
  *
- * This function will a a #spart at the start of the current cell's array by
+ * This function will add a #spart at the start of the current cell's array by
  * shifting all the #spart in the top-level cell by one position. All the
  * pointers and cell counts are updated accordingly.
  *
@@ -5557,7 +5591,7 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) {
     memmove(&c->stars.parts[1], &c->stars.parts[0],
             n_copy * sizeof(struct spart));
 
-    /* Update the gpart->spart links (shift by 1) */
+    /* Update the spart->gpart links (shift by 1) */
     for (size_t i = 0; i < n_copy; ++i) {
 #ifdef SWIFT_DEBUG_CHECKS
       if (c->stars.parts[i + 1].gpart == NULL) {
@@ -5610,6 +5644,140 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) {
   return sp;
 }
 
+/**
+ * @brief "Add" a #gpart in a given #cell.
+ *
+ * This function will add a #gpart at the start of the current cell's array by
+ * shifting all the #gpart in the top-level cell by one position. All the
+ * pointers and cell counts are updated accordingly.
+ *
+ * @param e The #engine.
+ * @param c The leaf-cell in which to add the #gpart.
+ *
+ * @return A pointer to the newly added #gpart. The gpart has a been zeroed
+ * and given a position within the cell as well as set to the minimal active
+ * time bin.
+ */
+struct gpart *cell_add_gpart(struct engine *e, struct cell *c) {
+  /* Perform some basic consitency checks */
+  if (c->nodeID != engine_rank) error("Adding gpart on a foreign node");
+  if (c->grav.ti_old_part != e->ti_current) error("Undrifted cell!");
+  if (c->split) error("Addition of gpart performed above the leaf level");
+
+  struct space *s = e->s;
+
+  /* Progeny number at each level */
+  int progeny[space_cell_maxdepth];
+#ifdef SWIFT_DEBUG_CHECKS
+  for (int i = 0; i < space_cell_maxdepth; ++i) progeny[i] = -1;
+#endif
+
+  /* Get the top-level this leaf cell is in and compute the progeny indices at
+     each level */
+  struct cell *top = c;
+  while (top->parent != NULL) {
+    /* What is the progeny index of the cell? */
+    for (int k = 0; k < 8; ++k) {
+      if (top->parent->progeny[k] == top) {
+        progeny[(int)top->parent->depth] = k;
+      }
+    }
+
+      /* Check that the cell was indeed drifted to this point to avoid future
+       * issues */
+#ifdef SWIFT_DEBUG_CHECKS
+    if (top->grav.super != NULL && top->grav.count > 0 &&
+        top->grav.ti_old_part != e->ti_current) {
+      error("Cell had not been correctly drifted before adding a gpart");
+    }
+#endif
+
+    /* Climb up */
+    top = top->parent;
+  }
+
+  /* Lock the top-level cell as we are going to operate on it */
+  lock_lock(&top->grav.star_formation_lock);
+
+  /* Are there any extra particles left? */
+  if (top->grav.count == top->grav.count_total - 1) {
+    /* Release the local lock before exiting. */
+    if (lock_unlock(&top->grav.star_formation_lock) != 0)
+      error("Failed to unlock the top-level cell.");
+    message("We ran out of gravity particles!");
+    atomic_inc(&e->forcerebuild);
+    return NULL;
+  }
+
+  /* Number of particles to shift in order to get a free space. */
+  const size_t n_copy = &top->grav.parts[top->grav.count] - c->grav.parts;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (c->grav.parts + n_copy > top->grav.parts + top->grav.count)
+    error("Copying beyond the allowed range");
+#endif
+
+  if (n_copy > 0) {
+    // MATTHIEU: This can be improved. We don't need to copy everything, just
+    // need to swap a few particles.
+    memmove(&c->grav.parts[1], &c->grav.parts[0],
+            n_copy * sizeof(struct gpart));
+
+    /* Update the gpart->spart links (shift by 1) */
+    struct gpart *gparts = c->grav.parts;
+    for (size_t i = 0; i < n_copy; ++i) {
+      if (gparts[i + 1].type == swift_type_gas) {
+        s->parts[-gparts[i + 1].id_or_neg_offset].gpart++;
+      } else if (gparts[i + 1].type == swift_type_stars) {
+        s->sparts[-gparts[i + 1].id_or_neg_offset].gpart++;
+      } else if (gparts[i + 1].type == swift_type_black_hole) {
+        s->bparts[-gparts[i + 1].id_or_neg_offset].gpart++;
+      }
+    }
+  }
+
+  /* Recursively shift all the gpart to get a free spot at the start of the
+   * current cell*/
+  cell_recursively_shift_gparts(top, progeny, /* main_branch=*/1);
+
+  /* Make sure the gravity will be recomputed for this particle in the next
+   * step
+   */
+  struct cell *top2 = c;
+  while (top2->parent != NULL) {
+    top2->grav.ti_old_part = e->ti_current;
+    top2 = top2->parent;
+  }
+  top2->grav.ti_old_part = e->ti_current;
+
+  /* Release the lock */
+  if (lock_unlock(&top->grav.star_formation_lock) != 0)
+    error("Failed to unlock the top-level cell.");
+
+  /* We now have an empty gpart as the first particle in that cell */
+  struct gpart *gp = &c->grav.parts[0];
+  bzero(gp, sizeof(struct gpart));
+
+  /* Give it a decent position */
+  gp->x[0] = c->loc[0] + 0.5 * c->width[0];
+  gp->x[1] = c->loc[1] + 0.5 * c->width[1];
+  gp->x[2] = c->loc[2] + 0.5 * c->width[2];
+
+  /* Set it to the current time-bin */
+  gp->time_bin = e->min_active_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Specify it was drifted to this point */
+  gp->ti_drift = e->ti_current;
+#endif
+
+  /* Register that we used one of the free slots. */
+  const size_t one = 1;
+  atomic_sub(&e->s->nr_extra_gparts, one);
+
+  return gp;
+}
+
 /**
  * @brief "Remove" a gas particle from the calculation.
  *
@@ -5935,6 +6103,84 @@ struct spart *cell_convert_part_to_spart(struct engine *e, struct cell *c,
   return sp;
 }
 
+/**
+ * @brief Add a new #spart based on a #part and link it to a new #gpart.
+ * The part and xpart are not changed.
+ *
+ * @param e The #engine.
+ * @param c The #cell from which to remove the #part.
+ * @param p The #part to remove (must be inside c).
+ * @param xp The extended data of the #part.
+ *
+ * @return A fresh #spart with the same ID, position, velocity and
+ * time-bin as the original #part.
+ */
+struct spart *cell_spawn_new_spart_from_part(struct engine *e, struct cell *c,
+                                             const struct part *p,
+                                             const struct xpart *xp) {
+  /* Quick cross-check */
+  if (c->nodeID != e->nodeID)
+    error("Can't spawn a particle in a foreign cell.");
+
+  if (p->gpart == NULL)
+    error("Trying to create a new spart from a part without gpart friend!");
+
+  /* Create a fresh (empty) spart */
+  struct spart *sp = cell_add_spart(e, c);
+
+  /* Did we run out of free spart slots? */
+  if (sp == NULL) return NULL;
+
+  /* Copy over the distance since rebuild */
+  sp->x_diff[0] = xp->x_diff[0];
+  sp->x_diff[1] = xp->x_diff[1];
+  sp->x_diff[2] = xp->x_diff[2];
+
+  /* Create a new gpart */
+  struct gpart *gp = cell_add_gpart(e, c);
+
+  /* Did we run out of free gpart slots? */
+  if (gp == NULL) {
+    /* Remove the particle created */
+    cell_remove_spart(e, c, sp);
+    return NULL;
+  }
+
+  /* Copy the gpart */
+  *gp = *p->gpart;
+
+  /* Assign the ID back */
+  sp->id = p->id;
+  gp->type = swift_type_stars;
+
+  /* Re-link things */
+  sp->gpart = gp;
+  gp->id_or_neg_offset = -(sp - e->s->sparts);
+
+  /* Synchronize clocks */
+  gp->time_bin = sp->time_bin;
+
+  /* Synchronize masses, positions and velocities */
+  sp->mass = p->mass;
+  sp->x[0] = p->x[0];
+  sp->x[1] = p->x[1];
+  sp->x[2] = p->x[2];
+  sp->v[0] = xp->v_full[0];
+  sp->v[1] = xp->v_full[1];
+  sp->v[2] = xp->v_full[2];
+
+#ifdef SWIFT_DEBUG_CHECKS
+  sp->ti_kick = p->ti_kick;
+  sp->ti_drift = p->ti_drift;
+#endif
+
+  /* Set a smoothing length */
+  sp->h = p->h;
+
+  /* Here comes the Sun! */
+  return sp;
+}
+
 /**
  * @brief Re-arrange the #part in a top-level cell such that all the extra
  * ones for on-the-fly creation are located at the end of the array.

src/cell.h

@@ -518,6 +518,9 @@ struct cell {
     /*! Spin lock for various uses (#multipole case). */
     swift_lock_type mlock;
 
+    /*! Spin lock for star formation use. */
+    swift_lock_type star_formation_lock;
+
     /*! Nr of #gpart in this cell. */
     int count;
 
@@ -932,6 +935,10 @@ void cell_remove_spart(const struct engine *e, struct cell *c,
 void cell_remove_bpart(const struct engine *e, struct cell *c,
                        struct bpart *bp);
 struct spart *cell_add_spart(struct engine *e, struct cell *c);
+struct gpart *cell_add_gpart(struct engine *e, struct cell *c);
+struct spart *cell_spawn_new_spart_from_part(struct engine *e, struct cell *c,
+                                             const struct part *p,
+                                             const struct xpart *xp);
 struct gpart *cell_convert_part_to_gpart(const struct engine *e, struct cell *c,
                                          struct part *p, struct xpart *xp);
 struct gpart *cell_convert_spart_to_gpart(const struct engine *e,

src/engine.h

@@ -547,6 +547,8 @@ void engine_launch(struct engine *e, const char *call);
 void engine_prepare(struct engine *e);
 void engine_init_particles(struct engine *e, int flag_entropy_ICs,
                            int clean_h_values);
+void engine_compute_star_formation_stats(struct engine *e,
+                                         struct star_formation *starform);
 void engine_step(struct engine *e);
 void engine_split(struct engine *e, struct partition *initial_partition);
 void engine_exchange_strays(struct engine *e, const size_t offset_parts,

src/runner_others.c

@@ -330,9 +330,17 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
 
             /* Convert the gas particle to a star particle */
             struct spart *sp = NULL;
+            const int spawn_spart =
+                star_formation_should_spawn_spart(p, xp, sf_props);
 
             if (swift_star_formation_model_creates_stars) {
-              sp = cell_convert_part_to_spart(e, c, p, xp);
+              /* Check if we should create a new particle or transform one */
+              if (spawn_spart) {
+                sp = cell_spawn_new_spart_from_part(e, c, p, xp);
+              } else {
+                /* Convert the gas particle to a star particle */
+                sp = cell_convert_part_to_spart(e, c, p, xp);
+              }
             } else {
               cell_convert_part_to_gpart(e, c, p, xp);
             }
@@ -344,9 +352,9 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
                * c->cellID); */
 
               /* Copy the properties of the gas particle to the star particle */
-              star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
-                                             with_cosmology, phys_const,
-                                             hydro_props, us, cooling);
+              star_formation_copy_properties(
+                  p, xp, sp, e, sf_props, cosmo, with_cosmology, phys_const,
+                  hydro_props, us, cooling, spawn_spart);
 
               /* Update the Star formation history */
               star_formation_logger_log_new_spart(sp, &c->stars.sfh);
@@ -550,7 +558,8 @@ void runner_do_end_grav_force(struct runner *r, struct cell *c, int timer) {
 
 #ifdef SWIFT_DEBUG_CHECKS
         if ((e->policy & engine_policy_self_gravity) &&
-            !(e->policy & engine_policy_black_holes)) {
+            !(e->policy & engine_policy_black_holes) &&
+            !(e->policy & engine_policy_star_formation)) {
 
           /* Let's add a self interaction to simplify the count */
           gp->num_interacted++;

src/space.c

@@ -577,12 +577,14 @@ void space_regrid(struct space *s, int verbose) {
         error("Failed to init spinlock for gravity.");
       if (lock_init(&s->cells_top[k].grav.mlock) != 0)
         error("Failed to init spinlock for multipoles.");
+      if (lock_init(&s->cells_top[k].grav.star_formation_lock) != 0)
+        error("Failed to init spinlock for star formation (gpart).");
       if (lock_init(&s->cells_top[k].stars.lock) != 0)
         error("Failed to init spinlock for stars.");
       if (lock_init(&s->cells_top[k].black_holes.lock) != 0)
         error("Failed to init spinlock for black holes.");
       if (lock_init(&s->cells_top[k].stars.star_formation_lock) != 0)
-        error("Failed to init spinlock for star formation.");
+        error("Failed to init spinlock for star formation (spart).");
     }
 
     /* Set the cell location and sizes. */
@@ -781,7 +783,7 @@ void space_allocate_extras(struct space *s, int verbose) {
 
   /* Do we have enough space for the extra gparts (i.e. we haven't used up any)
    * ? */
-  if (nr_gparts + expected_num_extra_gparts > size_gparts) {
+  if (nr_actual_gparts + expected_num_extra_gparts > nr_gparts) {
 
     /* Ok... need to put some more in the game */
 
@@ -817,6 +819,10 @@ void space_allocate_extras(struct space *s, int verbose) {
         if (s->sparts[i].time_bin != time_bin_not_created)
           s->sparts[i].gpart += delta;
       }
+      for (size_t i = 0; i < nr_bparts; ++i) {
+        if (s->bparts[i].time_bin != time_bin_not_created)
+          s->bparts[i].gpart += delta;
+      }
     }
 
     /* Turn some of the allocated spares into particles we can use */
@@ -3869,7 +3875,8 @@ void space_recycle(struct space *s, struct cell *c) {
   if (lock_destroy(&c->hydro.lock) != 0 || lock_destroy(&c->grav.plock) != 0 ||
       lock_destroy(&c->grav.mlock) != 0 || lock_destroy(&c->stars.lock) != 0 ||
       lock_destroy(&c->black_holes.lock) != 0 ||
-      lock_destroy(&c->stars.star_formation_lock))
+      lock_destroy(&c->grav.star_formation_lock) != 0 ||
+      lock_destroy(&c->stars.star_formation_lock) != 0)
     error("Failed to destroy spinlocks.");
 
   /* Lock the space. */
@@ -3922,7 +3929,8 @@ void space_recycle_list(struct space *s, struct cell *cell_list_begin,
         lock_destroy(&c->grav.mlock) != 0 ||
         lock_destroy(&c->stars.lock) != 0 ||
         lock_destroy(&c->black_holes.lock) != 0 ||
-        lock_destroy(&c->stars.star_formation_lock))
+        lock_destroy(&c->stars.star_formation_lock) != 0 ||
+        lock_destroy(&c->grav.star_formation_lock) != 0)
       error("Failed to destroy spinlocks.");
 
     /* Count this cell. */
@@ -4023,7 +4031,8 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
         lock_init(&cells[j]->grav.mlock) != 0 ||
         lock_init(&cells[j]->stars.lock) != 0 ||
         lock_init(&cells[j]->black_holes.lock) != 0 ||
-        lock_init(&cells[j]->stars.star_formation_lock) != 0)
+        lock_init(&cells[j]->stars.star_formation_lock) != 0 ||
+        lock_init(&cells[j]->grav.star_formation_lock) != 0)
       error("Failed to initialize cell spinlocks.");
   }
 }

src/space.h

@@ -40,6 +40,7 @@
 struct cell;
 struct cosmology;
 struct gravity_props;
+struct star_formation;
 
 /* Some constants. */
 #define space_cellallocchunk 1000
@@ -384,5 +385,7 @@ void space_free_foreign_parts(struct space *s, const int clear_cell_pointers);
 void space_struct_dump(struct space *s, FILE *stream);
 void space_struct_restore(struct space *s, FILE *stream);
 void space_write_cell_hierarchy(const struct space *s, int j);
+void space_compute_star_formation_stats(const struct space *s,
+                                        struct star_formation *star_form);
 
 #endif /* SWIFT_SPACE_H */

src/star_formation/EAGLE/star_formation.h

@@ -368,6 +368,21 @@ INLINE static int star_formation_should_convert_to_star(
   return (prob > random_number);
 }
 
+/**
+ * @brief Decides whether a new particle should be created or if the hydro
+ * particle needs to be transformed.
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param starform The properties of the star formation model.
+ *
+ * @return 1 if a new spart needs to be created.
+ */
+INLINE static int star_formation_should_spawn_spart(
+    struct part* p, struct xpart* xp, const struct star_formation* starform) {
+  return 0;
+}
+
 /**
  * @brief Update the SF properties of a particle that is not star forming.
  *
@@ -409,6 +424,7 @@ INLINE static void star_formation_update_part_not_SFR(
  * @param hydro_props The properties of the hydro scheme.
  * @param us The internal system of units.
  * @param cooling The cooling data struct.
+ * @param convert_part Did we convert a part (or spawned one)?
  */
 INLINE static void star_formation_copy_properties(
     const struct part* p, const struct xpart* xp, struct spart* sp,
@@ -417,7 +433,8 @@ 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,
+    const int convert_part) {
 
   /* Store the current mass */
   sp->mass = hydro_get_mass(p);
@@ -752,4 +769,19 @@ star_formation_no_spart_available(const struct engine* e, const struct part* p,
   /* Nothing to do, we just skip it and deal with it next step */
 }
 
+/**
+ * @brief Compute some information for the star formation model based
+ * on all the particles that were read in.
+ *
+ * This is called once on start-up of the code.
+ *
+ * Nothing to do here for EAGLE.
+ *
+ * @param star_form The #star_formation structure.
+ * @param e The #engine.
+ */
+__attribute__((always_inline)) INLINE static void
+star_formation_first_init_stats(struct star_formation* star_form,
+                                const struct engine* e) {}
+
 #endif /* SWIFT_EAGLE_STAR_FORMATION_H */

src/star_formation/GEAR/star_formation.h

@@ -145,11 +145,24 @@ INLINE static int star_formation_should_convert_to_star(
   const float G = phys_const->const_newton_G;
   const float density = hydro_get_physical_density(p, cosmo);
 
+  /* Get the mass of the future possible star */
+  const float mass_gas = hydro_get_mass(p);
+
   /* 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);
+  float prob = 1.f - expf(-starform->star_formation_efficiency *
+                          inv_free_fall_time * dt_star);
+
+  /* Add the mass factor */
+  if (starform->n_stars_per_part != 1) {
+    const float min_mass =
+        starform->mass_stars * starform->min_mass_frac_plus_one;
+    const float mass_star =
+        mass_gas > min_mass ? starform->mass_stars : mass_gas;
+
+    prob *= mass_gas / mass_star;
+  }
 
   /* Roll the dice... */
   const float random_number =
@@ -159,6 +172,29 @@ INLINE static int star_formation_should_convert_to_star(
   return random_number < prob;
 }
 
+/**
+ * @brief Decides whether a new particle should be created or if the hydro
+ * particle needs to be transformed.
+ *
+ * @param p The #part.
+ * @param xp The #xpart.
+ * @param starform The properties of the star formation model.
+ *
+ * @return 1 if a new spart needs to be created.
+ */
+INLINE static int star_formation_should_spawn_spart(
+    struct part* p, struct xpart* xp, const struct star_formation* starform) {
+
+  /* Check if we are splitting the particles or not */
+  if (starform->n_stars_per_part == 1) {
+    return 0;
+  }
+
+  const float mass_min =
+      starform->min_mass_frac_plus_one * starform->mass_stars;
+  return hydro_get_mass(p) > mass_min;
+}
+
 /**
  * @brief Update the SF properties of a particle that is not star forming.
  *
@@ -184,18 +220,35 @@ INLINE static void star_formation_update_part_not_SFR(
  * @param phys_const the physical constants in internal units.
  * @param cosmo the cosmological parameters and properties.
  * @param with_cosmology if we run with cosmology.
+ * @param convert_part Did we convert a part (or spawned one)?
  */
 INLINE static void star_formation_copy_properties(
-    const struct part* p, const struct xpart* xp, struct spart* sp,
+    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,
     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,
+    const int convert_part) {
 
   /* Store the current mass */
-  sp->mass = hydro_get_mass(p);
+  const float mass_gas = hydro_get_mass(p);
+  if (!convert_part) {
+    /* Update the spart */
+    const float min_mass =
+        starform->mass_stars * starform->min_mass_frac_plus_one;
+    const float mass_star =
+        mass_gas > min_mass ? starform->mass_stars : mass_gas;
+    sp->mass = mass_star;
+    sp->gpart->mass = mass_star;
+
+    /* Update the part */
+    hydro_set_mass(p, mass_gas - mass_star);
+    p->gpart->mass = mass_gas - mass_star;
+  } else {
+    sp->mass = mass_gas;
+  }
   sp->birth.mass = sp->mass;
 
   /* Store either the birth_scale_factor or birth_time depending  */
@@ -217,6 +270,9 @@ INLINE static void star_formation_copy_properties(
 
   /* Copy the chemistry properties */
   chemistry_copy_star_formation_properties(p, xp, sp);
+
+  /* Copy the progenitor id */
+  sp->prog_id = p->id;
 }
 
 /**
@@ -325,4 +381,38 @@ star_formation_no_spart_available(const struct engine* e, const struct part* p,
       "or Scheduler:cell_extra_gparts");
 }
 
+/**