diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
index 286b7591dd110616ff2f21d711469daf5317641a..0d93406eebc903c880a53ffa15f874e21902675b 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
@@ -51,6 +51,8 @@ Statistics:
Scheduler:
max_top_level_cells: 8
cell_split_size: 50
+ cell_extra_gparts: 100 # (Optional) Number of spare gparts per top-level allocated at rebuild time for on-the-fly creation.
+ cell_extra_sparts: 100 # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation.
# Parameters related to the initial conditions
InitialConditions:
@@ -96,22 +98,40 @@ EAGLEEntropyFloor:
Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
-
+
+GEARPressureFloor:
+ jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
+
+
# Cooling with Grackle 3.0
GrackleCooling:
- cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
- with_UV_background: 1 # Enable or not the UV background
- redshift: -1 # Redshift to use (-1 means time based redshift)
- with_metal_cooling: 1 # Enable or not the metal cooling
- provide_volumetric_heating_rates: 0 # (optional) User provide volumetric heating rates
- provide_specific_heating_rates: 0 # (optional) User provide specific heating rates
- self_shielding_method: 0 # (optional) Grackle (<= 3) or Gear self shielding method
- max_steps: 10000 # (optional) Max number of step when computing the initial composition
- convergence_limit: 1e-2 # (optional) Convergence threshold (relative) for initial composition
- thermal_time_myr: 5
+ cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+ with_UV_background: 1 # Enable or not the UV background
+ redshift: 0 # Redshift to use (-1 means time based redshift)
+ with_metal_cooling: 1 # Enable or not the metal cooling
+ provide_volumetric_heating_rates: 0 # (optional) User provide volumetric heating rates
+ provide_specific_heating_rates: 0 # (optional) User provide specific heating rates
+ max_steps: 10000 # (optional) Max number of step when computing the initial composition
+ convergence_limit: 1e-2 # (optional) Convergence threshold (relative) for initial composition
+ thermal_time_myr: 5 # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
+ self_shielding_method: -1 # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
+ self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
+
+# GEAR chemistry model (Revaz and Jablonka 2018)
GEARChemistry:
- initial_metallicity: 0.01295
+ initial_metallicity: 1 # Initial metallicity of the gas (mass fraction)
+ scale_initial_metallicity: 1 # Should we scale the initial metallicity with the solar one?
-GEARPressureFloor:
- jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
+ # 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
+ n_stars_per_particle: 4
+ min_mass_frac: 0.5
+
+# GEAR feedback model
+GEARFeedback:
+ supernovae_energy_erg: 0.1e51 # Energy released by a single supernovae.
+ yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 # Table containing the yields.
+ discrete_yields: 0 # Should we use discrete yields or the IMF integrated one?
diff --git a/src/cell.c b/src/cell.c
index fceaed8ac83a7754114179d35c7c8f91ad6cf262..802ce7dbf6e00893fafd02aedc65ca84da3b058c 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -233,6 +233,7 @@ int cell_link_gparts(struct cell *c, struct gpart *gparts) {
#endif
c->grav.parts = gparts;
+ c->grav.parts_rebuild = gparts;
/* Fill the progeny recursively, depth-first. */
if (c->split) {
@@ -1154,7 +1155,12 @@ int cell_pack_sf_counts(struct cell *restrict c,
pcells[0].stars.count = c->stars.count;
pcells[0].stars.dx_max_part = c->stars.dx_max_part;
+ /* Pack this cell's data. */
+ pcells[0].grav.delta_from_rebuild = c->grav.parts - c->grav.parts_rebuild;
+ pcells[0].grav.count = c->grav.count;
+
#ifdef SWIFT_DEBUG_CHECKS
+ /* Stars */
if (c->stars.parts_rebuild == NULL)
error("Star particles array at rebuild is NULL! c->depth=%d", c->depth);
@@ -1163,6 +1169,16 @@ int cell_pack_sf_counts(struct cell *restrict c,
if (pcells[0].stars.delta_from_rebuild > 0 && c->depth == 0)
error("Shifting the top-level pointer is not allowed!");
+
+ /* Grav */
+ if (c->grav.parts_rebuild == NULL)
+ error("Grav. particles array at rebuild is NULL! c->depth=%d", c->depth);
+
+ if (pcells[0].grav.delta_from_rebuild < 0)
+ error("Grav part pointer moved in the wrong direction!");
+
+ if (pcells[0].grav.delta_from_rebuild > 0 && c->depth == 0)
+ error("Shifting the top-level pointer is not allowed!");
#endif
/* Fill in the progeny, depth-first recursion. */
@@ -1198,6 +1214,8 @@ int cell_unpack_sf_counts(struct cell *restrict c,
#ifdef SWIFT_DEBUG_CHECKS
if (c->stars.parts_rebuild == NULL)
error("Star particles array at rebuild is NULL!");
+ if (c->grav.parts_rebuild == NULL)
+ error("Grav particles array at rebuild is NULL!");
#endif
/* Unpack this cell's data. */
@@ -1205,6 +1223,9 @@ int cell_unpack_sf_counts(struct cell *restrict c,
c->stars.parts = c->stars.parts_rebuild + pcells[0].stars.delta_from_rebuild;
c->stars.dx_max_part = pcells[0].stars.dx_max_part;
+ c->grav.count = pcells[0].grav.count;
+ c->grav.parts = c->grav.parts_rebuild + pcells[0].grav.delta_from_rebuild;
+
/* Fill in the progeny, depth-first recursion. */
int count = 1;
for (int k = 0; k < 8; k++)
@@ -1967,6 +1988,7 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
c->progeny[k]->grav.count = bucket_count[k];
c->progeny[k]->grav.count_total = c->progeny[k]->grav.count;
c->progeny[k]->grav.parts = &c->grav.parts[bucket_offset[k]];
+ c->progeny[k]->grav.parts_rebuild = c->progeny[k]->grav.parts;
}
}
@@ -5593,11 +5615,12 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) {
lock_lock(&top->stars.star_formation_lock);
/* Are there any extra particles left? */
- if (top->stars.count == top->stars.count_total - 1) {
+ if (top->stars.count >= top->stars.count_total - 1) {
/* Release the local lock before exiting. */
if (lock_unlock(&top->stars.star_formation_lock) != 0)
error("Failed to unlock the top-level cell.");
- message("We ran out of star particles!");
+ if (top->stars.count == top->stars.count_total - 1)
+ message("We ran out of star particles!");
atomic_inc(&e->forcerebuild);
return NULL;
}
@@ -5725,11 +5748,12 @@ struct gpart *cell_add_gpart(struct engine *e, struct cell *c) {
lock_lock(&top->grav.star_formation_lock);
/* Are there any extra particles left? */
- if (top->grav.count == top->grav.count_total - 1) {
+ 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!");
+ if (top->grav.count == top->grav.count_total - 1)
+ message("We ran out of gravity particles!");
atomic_inc(&e->forcerebuild);
return NULL;
}
diff --git a/src/cell.h b/src/cell.h
index a17b2a8e5794d0b092ca47acd9e6b87877820fa5..e84e00bcd77972be5e2a89d89b7d180a96bd637a 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -273,6 +273,17 @@ struct pcell_sf {
float dx_max_part;
} stars;
+
+ /*! Grav. variables */
+ struct {
+
+ /* Distance by which the gpart pointer has moved since the last rebuild */
+ ptrdiff_t delta_from_rebuild;
+
+ /* Number of particles in the cell */
+ int count;
+
+ } grav;
};
/**
@@ -461,6 +472,9 @@ struct cell {
/*! Pointer to the #gpart data. */
struct gpart *parts;
+ /*! Pointer to the #spart data at rebuild time. */
+ struct gpart *parts_rebuild;
+
/*! This cell's multipole. */
struct gravity_tensors *multipole;
diff --git a/src/runner_others.c b/src/runner_others.c
index 0a496a6443c0993327be410a127f3a117e801517..5db756ac2b859d473fa9f07538e53d2ad6b31332 100644
--- a/src/runner_others.c
+++ b/src/runner_others.c
@@ -267,6 +267,14 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
/* Update current cell using child cells */
star_formation_logger_add(&c->stars.sfh, &cp->stars.sfh);
+
+ /* Update the dx_max */
+ if (star_formation_need_update_dx_max) {
+ c->hydro.dx_max_part =
+ max(cp->hydro.dx_max_part, c->hydro.dx_max_part);
+ c->hydro.dx_max_sort =
+ max(cp->hydro.dx_max_sort, c->hydro.dx_max_sort);
+ }
}
} else {
@@ -359,6 +367,25 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
/* Update the Star formation history */
star_formation_logger_log_new_spart(sp, &c->stars.sfh);
+ /* Update the displacement information */
+ if (star_formation_need_update_dx_max) {
+ const float dx2_part = xp->x_diff[0] * xp->x_diff[0] +
+ xp->x_diff[1] * xp->x_diff[1] +
+ xp->x_diff[2] * xp->x_diff[2];
+ const float dx2_sort = xp->x_diff_sort[0] * xp->x_diff_sort[0] +
+ xp->x_diff_sort[1] * xp->x_diff_sort[1] +
+ xp->x_diff_sort[2] * xp->x_diff_sort[2];
+
+ const float dx_part = sqrtf(dx2_part);
+ const float dx_sort = sqrtf(dx2_sort);
+
+ /* No need to climb up the tree as the star formation starts
+ higher in the hierarchy than the hydro and goes to the bottom
+ */
+ c->hydro.dx_max_part = max(c->hydro.dx_max_part, dx_part);
+ c->hydro.dx_max_sort = max(c->hydro.dx_max_sort, dx_sort);
+ }
+
#ifdef WITH_LOGGER
/* Copy the properties back to the stellar particle */
sp->logger_data = xp->logger_data;
diff --git a/src/space.c b/src/space.c
index 765d1d69decb1f59411642cd9a52cb29af32c678..a625b48346670b2983cded7cd9db35a364ec1d94 100644
--- a/src/space.c
+++ b/src/space.c
@@ -265,6 +265,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->hydro.parts = NULL;
c->hydro.xparts = NULL;
c->grav.parts = NULL;
+ c->grav.parts_rebuild = NULL;
c->stars.parts = NULL;
c->stars.parts_rebuild = NULL;
c->black_holes.parts = NULL;
@@ -1954,6 +1955,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* Store the state at rebuild time */
c->stars.parts_rebuild = c->stars.parts;
+ c->grav.parts_rebuild = c->grav.parts;
c->hydro.count_total = c->hydro.count + space_extra_parts;
c->grav.count_total = c->grav.count + space_extra_gparts;
diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h
index f8c42436f3f5b40941718f33e0e22a2b9eef9556..9c9c20b2616ea4c04e19ad6cf4ce07ac9ab389e9 100644
--- a/src/star_formation/EAGLE/star_formation.h
+++ b/src/star_formation/EAGLE/star_formation.h
@@ -41,6 +41,8 @@
* @brief Star formation model used in the EAGLE model
*/
+#define star_formation_need_update_dx_max 0
+
/**
* @brief Properties of the EAGLE star formation model.
*/
diff --git a/src/star_formation/GEAR/star_formation.h b/src/star_formation/GEAR/star_formation.h
index 82a43c9e1ff2ce58f766a1c5e38f837093176987..3ef57147567d6b9d7ef150ffba5893ef9abd0929 100644
--- a/src/star_formation/GEAR/star_formation.h
+++ b/src/star_formation/GEAR/star_formation.h
@@ -34,6 +34,8 @@
#include "star_formation_struct.h"
#include "units.h"
+#define star_formation_need_update_dx_max 1
+
/**
* @brief Calculate if the gas has the potential of becoming
* a star.
@@ -207,6 +209,77 @@ INLINE static void star_formation_update_part_not_SFR(
struct part* p, struct xpart* xp, const struct engine* e,
const struct star_formation* starform, const int with_cosmology) {}
+/**
+ * @brief Separate the #spart and #part by randomly moving both of them.
+ *
+ * @param e The #engine.
+ * @param p The #part generating a star.
+ * @param xp The #xpart generating a star.
+ * @param sp The new #spart.
+ */
+INLINE static void star_formation_separate_particles(const struct engine* e,
+ struct part* p,
+ struct xpart* xp,
+ struct spart* sp) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (p->x[0] != sp->x[0] || p->x[1] != sp->x[1] || p->x[2] != sp->x[2]) {
+ error(
+ "Moving particles that are not at the same location."
+ " (%g, %g, %g) - (%g, %g, %g)",
+ p->x[0], p->x[1], p->x[2], sp->x[0], sp->x[1], sp->x[2]);
+ }
+#endif
+
+ /* Move a bit the particle in order to avoid
+ division by 0.
+ */
+ const float max_displacement = 0.2;
+ const double delta_x =
+ 2.f * random_unit_interval(p->id, e->ti_current,
+ (enum random_number_type)0) -
+ 1.f;
+ const double delta_y =
+ 2.f * random_unit_interval(p->id, e->ti_current,
+ (enum random_number_type)1) -
+ 1.f;
+ const double delta_z =
+ 2.f * random_unit_interval(p->id, e->ti_current,
+ (enum random_number_type)2) -
+ 1.f;
+
+ sp->x[0] += delta_x * max_displacement * p->h;
+ sp->x[1] += delta_y * max_displacement * p->h;
+ sp->x[2] += delta_z * max_displacement * p->h;
+
+ /* Copy the position to the gpart */
+ sp->gpart->x[0] = sp->x[0];
+ sp->gpart->x[1] = sp->x[1];
+ sp->gpart->x[2] = sp->x[2];
+
+ /* Do the gas particle. */
+ const double mass_ratio = sp->mass / hydro_get_mass(p);
+ const double dx[3] = {mass_ratio * delta_x * max_displacement * p->h,
+ mass_ratio * delta_y * max_displacement * p->h,
+ mass_ratio * delta_z * max_displacement * p->h};
+
+ p->x[0] -= dx[0];
+ p->x[1] -= dx[1];
+ p->x[2] -= dx[2];
+
+ /* Compute offsets since last cell construction */
+ xp->x_diff[0] += dx[0];
+ xp->x_diff[1] += dx[1];
+ xp->x_diff[1] += dx[2];
+ xp->x_diff_sort[0] += dx[0];
+ xp->x_diff_sort[1] += dx[1];
+ xp->x_diff_sort[2] += dx[2];
+
+ /* Copy the position to the gpart */
+ p->gpart->x[0] = p->x[0];
+ p->gpart->x[1] = p->x[1];
+ p->gpart->x[2] = p->x[2];
+}
+
/**
* @brief Copies the properties of the gas particle over to the
* star particle.
@@ -225,10 +298,9 @@ INLINE static void star_formation_update_part_not_SFR(
* @param convert_part Did we convert a part (or spawned one)?
*/
INLINE static void star_formation_copy_properties(
- 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,
+ struct part* p, 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,
@@ -247,6 +319,8 @@ INLINE static void star_formation_copy_properties(
/* Update the part */
hydro_set_mass(p, new_mass_gas);
p->gpart->mass = new_mass_gas;
+
+ star_formation_separate_particles(e, p, xp, sp);
} else {
sp->mass = mass_gas;
}
diff --git a/src/star_formation/QLA/star_formation.h b/src/star_formation/QLA/star_formation.h
index aa9421220830df556bf0b56f50e251cbab8c8a7d..5f8a36f2103256464d7a3f469da09fc43a4b5b2a 100644
--- a/src/star_formation/QLA/star_formation.h
+++ b/src/star_formation/QLA/star_formation.h
@@ -34,6 +34,8 @@
* @brief Star formation model used in the EAGLE model
*/
+#define star_formation_need_update_dx_max 0
+
/**
* @brief Properties of the EAGLE star formation model.
*/
diff --git a/src/star_formation/none/star_formation.h b/src/star_formation/none/star_formation.h
index 78f8a7c61789ce0075bf62a2c099e2cc534a5914..c40a2f93b6e2e6aa6cd873a3220d4c4c48c90919 100644
--- a/src/star_formation/none/star_formation.h
+++ b/src/star_formation/none/star_formation.h
@@ -32,6 +32,11 @@
/* Starformation struct */
struct star_formation {};
+/* Does the star formation model move the hydro particles?
+ This will update the c->hydro.dx_max_part and
+ c->hydro.dx_max_sort after forming a star. */
+#define star_formation_need_update_dx_max 0
+
/**
* @brief Calculate if the gas has the potential of becoming
* a star.