From 41f2ac04f5fa9bb1fc2ebdc223894a6b689a8a47 Mon Sep 17 00:00:00 2001
From: Loic Hausammann <loic.hausammann@protonmail.ch>
Date: Mon, 16 Mar 2020 12:24:48 +0000
Subject: [PATCH] Gear star formation split
---
examples/main.c | 5 +
src/cell.c | 250 +++++++++++++++++-
src/cell.h | 7 +
src/engine.h | 2 +
src/runner_others.c | 19 +-
src/space.c | 19 +-
src/space.h | 3 +
src/star_formation/EAGLE/star_formation.h | 34 ++-
src/star_formation/GEAR/star_formation.h | 100 ++++++-
src/star_formation/GEAR/star_formation_io.h | 13 +
.../GEAR/star_formation_struct.h | 15 +-
src/star_formation/QLA/star_formation.h | 34 ++-
src/star_formation/none/star_formation.h | 35 ++-
src/stars/GEAR/stars_io.h | 6 +-
src/stars/GEAR/stars_part.h | 3 +
15 files changed, 521 insertions(+), 24 deletions(-)
diff --git a/examples/main.c b/examples/main.c
index 4784d532f4..2528ce589e 100644
--- a/examples/main.c
+++ b/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] +
diff --git a/src/cell.c b/src/cell.c
index d5b0499d00..65d0c1ad1b 100644
--- a/src/cell.c
+++ b/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.
diff --git a/src/cell.h b/src/cell.h
index b4242bf1f5..f449b92907 100644
--- a/src/cell.h
+++ b/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,
diff --git a/src/engine.h b/src/engine.h
index 27cfd8dd0f..5c3ccc7980 100644
--- a/src/engine.h
+++ b/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,
diff --git a/src/runner_others.c b/src/runner_others.c
index f6d5d6bd5e..6e4f353e4e 100644
--- a/src/runner_others.c
+++ b/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++;
diff --git a/src/space.c b/src/space.c
index b63a261aab..b6d38f88f1 100644
--- a/src/space.c
+++ b/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.");
}
}
diff --git a/src/space.h b/src/space.h
index 885f9ba85e..e98b1b1a0e 100644
--- a/src/space.h
+++ b/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 */
diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h
index 77b3acf3b0..15be7ea6b7 100644
--- a/src/star_formation/EAGLE/star_formation.h
+++ b/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 */
diff --git a/src/star_formation/GEAR/star_formation.h b/src/star_formation/GEAR/star_formation.h
index f2afc77f9d..27980271b3 100644
--- a/src/star_formation/GEAR/star_formation.h
+++ b/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");
}
+/**
+ * @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.
+ *
+ * @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) {
+
+ const struct space* s = e->s;
+ double avg_mass = 0;
+
+ /* Sum the mass over all the particles. */
+ for (size_t i = 0; i < s->nr_parts; i++) {
+ avg_mass += hydro_get_mass(&s->parts[i]);
+ }
+
+#ifdef WITH_MPI
+ /* Compute the contribution from the other nodes. */
+ MPI_Allreduce(MPI_IN_PLACE, &avg_mass, 1, MPI_DOUBLE, MPI_SUM,
+ MPI_COMM_WORLD);
+#endif
+
+ star_form->mass_stars = avg_mass / e->total_nr_parts;
+
+ if (e->nodeID == 0) {
+ message("Average hydro mass: %g", star_form->mass_stars);
+ }
+}
+
#endif /* SWIFT_GEAR_STAR_FORMATION_H */
diff --git a/src/star_formation/GEAR/star_formation_io.h b/src/star_formation/GEAR/star_formation_io.h
index 7e6dcefce6..f99458a931 100644
--- a/src/star_formation/GEAR/star_formation_io.h
+++ b/src/star_formation/GEAR/star_formation_io.h
@@ -64,12 +64,25 @@ INLINE static void starformation_init_backend(
starform->maximal_temperature = parser_get_param_double(
parameter_file, "GEARStarFormation:maximal_temperature");
+ /* Number of stars per particles */
+ starform->n_stars_per_part = parser_get_param_double(
+ parameter_file, "GEARStarFormation:n_stars_per_particle");
+
+ /* Minimal fraction of mass for the last star formed. */
+ starform->min_mass_frac_plus_one = parser_get_param_double(
+ parameter_file, "GEARStarFormation:min_mass_frac");
+ /* Avoid generating gas particle with mass below the fraction => + 1. */
+ starform->min_mass_frac_plus_one += 1.;
+
/* Get the jeans factor */
starform->n_jeans_2_3 = pow(pressure_floor_props.n_jeans, 2. / 3.);
/* Apply unit change */
starform->maximal_temperature *=
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+
+ /* Initialize the mass of the stars to 0 for the stats computation */
+ starform->mass_stars = 0;
}
#endif /* SWIFT_STAR_FORMATION_GEAR_IO_H */
diff --git a/src/star_formation/GEAR/star_formation_struct.h b/src/star_formation/GEAR/star_formation_struct.h
index 9d6d78df5c..c03a8ae0a4 100644
--- a/src/star_formation/GEAR/star_formation_struct.h
+++ b/src/star_formation/GEAR/star_formation_struct.h
@@ -34,14 +34,23 @@ struct star_formation_xpart_data {
struct star_formation {
/*! Number of particle required to resolved the
- * Jeans criterion (at power 2/3) */
+ * Jeans criterion (at power 2/3). */
float n_jeans_2_3;
- /*! Maximal temperature for forming a star */
+ /*! Maximal temperature for forming a star. */
float maximal_temperature;
- /*! Star formation efficiency */
+ /*! Star formation efficiency. */
float star_formation_efficiency;
+
+ /*! Number of possible stars per particle. */
+ int n_stars_per_part;
+
+ /*! Mass of a star. */
+ float mass_stars;
+
+ /*! Minimal fraction of mass_stars for the last star formed by a part. */
+ float min_mass_frac_plus_one;
};
#endif /* SWIFT_GEAR_STAR_FORMATION_STRUCT_H */
diff --git a/src/star_formation/QLA/star_formation.h b/src/star_formation/QLA/star_formation.h
index a0ffa22707..aa94212208 100644
--- a/src/star_formation/QLA/star_formation.h
+++ b/src/star_formation/QLA/star_formation.h
@@ -151,6 +151,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,
@@ -159,7 +160,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);
@@ -297,4 +299,34 @@ star_formation_no_spart_available(const struct engine* e, const struct part* p,
/* Nothing to do since we just turn gas particles into DM. */
}
+/**
+ * @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 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 the quick Lyman-alpha model.
+ *
+ * @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_QLA_STAR_FORMATION_H */
diff --git a/src/star_formation/none/star_formation.h b/src/star_formation/none/star_formation.h
index 589cd0f108..78f8a7c617 100644
--- a/src/star_formation/none/star_formation.h
+++ b/src/star_formation/none/star_formation.h
@@ -101,6 +101,21 @@ INLINE static int star_formation_should_convert_to_star(
return 0;
}
+/**
+ * @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.
*
@@ -130,6 +145,7 @@ 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,
@@ -138,7 +154,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) {}
/**
* @brief initialization of the star formation law
@@ -250,4 +267,20 @@ star_formation_no_spart_available(const struct engine* e, const struct part* p,
const struct xpart* xp) {
/* 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.
+ *
+ * @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_NONE_STAR_FORMATION_H */
diff --git a/src/stars/GEAR/stars_io.h b/src/stars/GEAR/stars_io.h
index a5a6e7368b..83d3f8e418 100644
--- a/src/stars/GEAR/stars_io.h
+++ b/src/stars/GEAR/stars_io.h
@@ -119,7 +119,7 @@ INLINE static void stars_write_particles(const struct spart *sparts,
const int with_cosmology) {
/* Say how much we want to write */
- *num_fields = 9;
+ *num_fields = 10;
/* List what we want to write */
list[0] = io_make_output_field_convert_spart(
@@ -169,6 +169,10 @@ INLINE static void stars_write_particles(const struct spart *sparts,
sparts, birth.mass,
"Masses of the star particles at birth time");
+ list[9] = io_make_output_field("ProgenitorIDs", LONGLONG, 1,
+ UNIT_CONV_NO_UNITS, 0.f, sparts, prog_id,
+ "Unique IDs of the progenitor particle");
+
#ifdef DEBUG_INTERACTIONS_STARS
list += *num_fields;
diff --git a/src/stars/GEAR/stars_part.h b/src/stars/GEAR/stars_part.h
index bee9b5c6a3..0e1e8a983b 100644
--- a/src/stars/GEAR/stars_part.h
+++ b/src/stars/GEAR/stars_part.h
@@ -37,6 +37,9 @@ struct spart {
/*! Particle ID. */
long long id;
+ /*! Progenitor ID */
+ long long prog_id;
+
/*! Pointer to corresponding gravity part. */
struct gpart* gpart;
--
GitLab