diff --git a/doc/RTD/source/ParameterFiles/parameter_description.rst b/doc/RTD/source/ParameterFiles/parameter_description.rst
index 9f2ed813accc3e4aa91b2b916fb0c27a391544ee..b90c2c26b8fe01ad306106fd3fd2f65612002d8f 100644
--- a/doc/RTD/source/ParameterFiles/parameter_description.rst
+++ b/doc/RTD/source/ParameterFiles/parameter_description.rst
@@ -923,7 +923,7 @@ at some arbitrary steps, and indeed do better than the initial partition
earlier in the run. This can be done using *fixed cost* repartitioning.
Fixed costs are output during each repartition step into the file
-`partition_fixed_costs.h`, this should be created by a test run of your your
+`partition_fixed_costs.h`, this should be created by a test run of your
full simulation (with possibly with a smaller volume, but all the physics
enabled). This file can then be used to replace the same file found in the
`src/` directory and SWIFT should then be recompiled. Once you have that, you
diff --git a/examples/Cooling/ConstantCosmoTempEvolution/const_cosmo_temp_evol.yml b/examples/Cooling/ConstantCosmoTempEvolution/const_cosmo_temp_evol.yml
index 01dcbc96fd329f53b1fd345bc198b05a53861982..50a95809873aa020257c8ef8ebf2444ac0be8be7 100644
--- a/examples/Cooling/ConstantCosmoTempEvolution/const_cosmo_temp_evol.yml
+++ b/examples/Cooling/ConstantCosmoTempEvolution/const_cosmo_temp_evol.yml
@@ -13,7 +13,7 @@ Cosmology:
a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
+ Omega_b: 0.0482519 # Baryon density parameter
# Parameters governing the time integration
TimeIntegration:
diff --git a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
index 97ed3e4ee1299fd25a148c9c2de6f6af6547882c..e94dce66f80b9b7c669e9095b7ce6f76aaf84551 100644
--- a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
@@ -17,7 +17,7 @@ Cosmology:
a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
+ Omega_b: 0.0482519 # Baryon density parameter
# Parameters governing the time integration
TimeIntegration:
diff --git a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
index d83bbd38fbb0d4e572d17e5107f94425d4427b31..4e6f7f6c49040fb2f39f443c7c1c7235570566f3 100644
--- a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
@@ -17,7 +17,7 @@ Cosmology:
a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
+ Omega_b: 0.0482519 # Baryon density parameter
# Parameters governing the time integration
TimeIntegration:
diff --git a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
index a3130edcffe01566d34cf19d39976aff94b65907..9317e9603b106141e9c7a82b9cd2fc36f63698a8 100644
--- a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
@@ -17,7 +17,7 @@ Cosmology:
a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
+ Omega_b: 0.0482519 # Baryon density parameter
# Parameters governing the time integration
TimeIntegration:
diff --git a/examples/EAGLE_low_z/EAGLE_100/eagle_100.yml b/examples/EAGLE_low_z/EAGLE_100/eagle_100.yml
index 5275f709c710d2da1d5396e98f1d4d918e482c6d..746362eb7d573eed75018423b8e33f6835916ccd 100644
--- a/examples/EAGLE_low_z/EAGLE_100/eagle_100.yml
+++ b/examples/EAGLE_low_z/EAGLE_100/eagle_100.yml
@@ -13,7 +13,7 @@ Cosmology:
a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
+ Omega_b: 0.0482519 # Baryon density parameter
# Parameters governing the time integration
TimeIntegration:
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index 685d24d20f14f567175e1d7086a65455017199e0..da6dfc0d094e45146e8bed3e536ca7a0874093fb 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -13,7 +13,7 @@ Cosmology:
a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
+ Omega_b: 0.0482519 # Baryon density parameter
# Parameters governing the time integration
TimeIntegration:
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index b782631d7e264e582455f70b4fb399f7db1f4731..1c2f62ef75209d8e6b3ed787b97d458e6b567dc6 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -21,7 +21,7 @@ Cosmology:
a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
+ Omega_b: 0.0482519 # Baryon density parameter
Scheduler:
links_per_tasks: 500
diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
index 53afda9445b4f8ee7c258b5c42af2a605eba4755..14844f49a5fe56607e0bd95dd57d44c42cb43f9d 100644
--- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
@@ -13,7 +13,7 @@ Cosmology:
a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
+ Omega_b: 0.0482519 # Baryon density parameter
# Parameters governing the time integration
TimeIntegration:
diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index 694c1c111d8f5b2455d4cfcdc567b8854b12ff2c..9aef75d741efa148f2397ba172d1ca6c86dd20d9 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -21,7 +21,7 @@ Cosmology:
a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
+ Omega_b: 0.0482519 # Baryon density parameter
Scheduler:
max_top_level_cells: 8
diff --git a/src/cell.c b/src/cell.c
index bbc9da0ce3f64e66322bb8a278ebefa5a75bd5d3..6af078fcbcfddfc227ddca61eb7dac18f264f410 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -71,6 +71,8 @@
#include "tools.h"
#include "tracers.h"
+extern int engine_star_resort_task_depth;
+
/* Global variables. */
int cell_next_tag = 0;
@@ -2404,6 +2406,77 @@ void cell_clear_limiter_flags(struct cell *c, void *data) {
cell_flag_do_hydro_limiter | cell_flag_do_hydro_sub_limiter);
}
+/**
+ * @brief Recursively clear the stars_resort flag in a cell hierarchy.
+ *
+ * @param c The #cell to act on.
+ */
+void cell_set_star_resort_flag(struct cell *c) {
+
+ cell_set_flag(c, cell_flag_do_stars_resort);
+
+ /* Abort if we reched the level where the resorting task lives */
+ if (c->depth == engine_star_resort_task_depth || c->hydro.super == c) return;
+
+ if (c->split) {
+ for (int k = 0; k < 8; ++k)
+ if (c->progeny[k] != NULL) cell_set_star_resort_flag(c->progeny[k]);
+ }
+}
+
+/**
+ * @brief Recurses in a cell hierarchy down to the level where the
+ * star resort tasks are and activates them.
+ *
+ * The function will fail if called *below* the super-level
+ *
+ * @param c The #cell to recurse into.
+ * @param s The #scheduler.
+ */
+void cell_activate_star_resort_tasks(struct cell *c, struct scheduler *s) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->hydro.super != NULL && c->hydro.super != c)
+ error("Function called below the super level!");
+#endif
+
+ /* The resort tasks are at either the chosen depth or the super level,
+ * whichever comes first. */
+ if (c->depth == engine_star_resort_task_depth || c->hydro.super == c) {
+ scheduler_activate(s, c->hydro.stars_resort);
+ } else {
+ for (int k = 0; k < 8; ++k) {
+ if (c->progeny[k] != NULL) {
+ cell_activate_star_resort_tasks(c->progeny[k], s);
+ }
+ }
+ }
+}
+
+/**
+ * @brief Activate the star formation task as well as the resorting of stars
+ *
+ * Must be called at the top-level in the tree (where the SF task is...)
+ *
+ * @param c The (top-level) #cell.
+ * @param s The #scheduler.
+ */
+void cell_activate_star_formation_tasks(struct cell *c, struct scheduler *s) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->depth != 0) error("Function should be called at the top-level only");
+#endif
+
+ /* Have we already unskipped that task? */
+ if (c->hydro.star_formation->skip == 0) return;
+
+ /* Activate the star formation task */
+ scheduler_activate(s, c->hydro.star_formation);
+
+ /* Activate the star resort tasks at whatever level they are */
+ cell_activate_star_resort_tasks(c, s);
+}
+
/**
* @brief Recurse down in a cell hierarchy until the hydro.super level is
* reached and activate the spart drift at that level.
@@ -3437,9 +3510,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
if (c->logger != NULL) scheduler_activate(s, c->logger);
if (c->top->hydro.star_formation != NULL) {
- scheduler_activate(s, c->top->hydro.star_formation);
- scheduler_activate(s, c->top->hydro.stars_resort);
- cell_activate_drift_spart(c, s);
+ cell_activate_star_formation_tasks(c->top, s);
}
}
@@ -5169,10 +5240,17 @@ void cell_recursively_shift_sparts(struct cell *c,
/* 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)
+ if (main_branch) {
c->stars.count++;
- else
+
+ /* Indicate that the cell is not sorted and cancel the pointer sorting
+ * arrays. */
+ c->stars.sorted = 0;
+ cell_free_stars_sorts(c);
+
+ } else {
c->stars.parts++;
+ }
}
/**
diff --git a/src/cell.h b/src/cell.h
index c6d32c8219b31cadfd20a7110226aa7e15badcf9..3bdae3395452df21455068fa8e7a2707749f067b 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -876,6 +876,8 @@ void cell_drift_multipole(struct cell *c, const struct engine *e);
void cell_drift_all_multipoles(struct cell *c, const struct engine *e);
void cell_check_timesteps(struct cell *c);
void cell_store_pre_drift_values(struct cell *c);
+void cell_set_star_resort_flag(struct cell *c);
+void cell_activate_star_formation_tasks(struct cell *c, struct scheduler *s);
void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
struct scheduler *s);
void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
diff --git a/src/engine.c b/src/engine.c
index 7f6a1639c033ed1e30aefdaced7ae839cb5721e8..145f55c75440ca52a49e6668bc5be8e04ae2c9ef 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3327,7 +3327,7 @@ void engine_skip_force_and_kick(struct engine *e) {
t->type == task_type_drift_gpart_out || t->type == task_type_cooling ||
t->type == task_type_stars_in || t->type == task_type_stars_out ||
t->type == task_type_star_formation ||
- t->type == task_type_extra_ghost ||
+ t->type == task_type_stars_resort || t->type == task_type_extra_ghost ||
t->type == task_type_bh_swallow_ghost1 ||
t->type == task_type_bh_swallow_ghost2 ||
t->type == task_type_bh_swallow_ghost3 ||
diff --git a/src/engine.h b/src/engine.h
index 967a430871648ac4cb91f390f8b85675a314d3a8..dc5c2c9ddd9a843d3e677e9b56c47f7451515b2f 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -113,6 +113,7 @@ enum engine_step_properties {
#define engine_default_timesteps_file_name "timesteps"
#define engine_max_parts_per_ghost_default 1000
#define engine_max_sparts_per_ghost_default 1000
+#define engine_star_resort_task_depth_default 2
#define engine_tasks_per_cell_margin 1.2
/**
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index 036c5aef0adb44a47777a4ef3c661368819096c2..6a8341068cb87d2994a8d658a9a7255462ef294e 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -54,6 +54,7 @@
extern int engine_max_parts_per_ghost;
extern int engine_max_sparts_per_ghost;
+extern int engine_star_resort_task_depth;
/**
* @brief Add send tasks for the gravity pairs to a hierarchy of cells.
@@ -775,9 +776,6 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
if (with_star_formation && c->hydro.count > 0) {
c->hydro.star_formation = scheduler_addtask(
s, task_type_star_formation, task_subtype_none, 0, 0, c, NULL);
- c->hydro.stars_resort = scheduler_addtask(
- s, task_type_stars_resort, task_subtype_none, 0, 0, c, NULL);
- scheduler_addunlock(s, c->hydro.star_formation, c->hydro.stars_resort);
}
}
@@ -979,8 +977,11 @@ void engine_add_ghosts(struct engine *e, struct cell *c, struct task *ghost_in,
*
* @param e The #engine.
* @param c The #cell.
+ * @param star_resort_cell Pointer to the cell where the star_resort task has
+ * been created. NULL above that level or if not running with star formation.
*/
-void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
+void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
+ struct cell *star_resort_cell) {
struct scheduler *s = &e->sched;
const int with_stars = (e->policy & engine_policy_stars);
@@ -989,6 +990,25 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
const int with_star_formation = (e->policy & engine_policy_star_formation);
const int with_black_holes = (e->policy & engine_policy_black_holes);
+ /* Are we are the level where we create the stars' resort tasks?
+ * If the tree is shallow, we need to do this at the super-level if the
+ * super-level is above the level we want */
+ if ((c->nodeID == e->nodeID) && (star_resort_cell == NULL) &&
+ (c->depth == engine_star_resort_task_depth || c->hydro.super == c)) {
+
+ if (with_star_formation && c->hydro.count > 0) {
+
+ /* Record this is the level where we re-sort */
+ star_resort_cell = c;
+
+ c->hydro.stars_resort = scheduler_addtask(
+ s, task_type_stars_resort, task_subtype_none, 0, 0, c, NULL);
+
+ scheduler_addunlock(s, c->top->hydro.star_formation,
+ c->hydro.stars_resort);
+ }
+ }
+
/* Are we in a super-cell ? */
if (c->hydro.super == c) {
@@ -1078,7 +1098,8 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
scheduler_addunlock(s, c->stars.stars_out, c->super->timestep);
if (with_star_formation && c->hydro.count > 0) {
- scheduler_addunlock(s, c->top->hydro.stars_resort, c->stars.stars_in);
+ scheduler_addunlock(s, star_resort_cell->hydro.stars_resort,
+ c->stars.stars_in);
}
}
@@ -1114,7 +1135,8 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
if (c->split)
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL)
- engine_make_hierarchical_tasks_hydro(e, c->progeny[k]);
+ engine_make_hierarchical_tasks_hydro(e, c->progeny[k],
+ star_resort_cell);
}
}
@@ -1131,7 +1153,8 @@ void engine_make_hierarchical_tasks_mapper(void *map_data, int num_elements,
/* Make the common tasks (time integration) */
engine_make_hierarchical_tasks_common(e, c);
/* Add the hydro stuff */
- if (with_hydro) engine_make_hierarchical_tasks_hydro(e, c);
+ if (with_hydro)
+ engine_make_hierarchical_tasks_hydro(e, c, /*star_resort_cell=*/NULL);
/* And the gravity stuff */
if (with_self_gravity || with_ext_gravity)
engine_make_hierarchical_tasks_gravity(e, c);
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index dcf1be29a677ce69d9fa471f7f734583ee6842eb..7f4e2d4c4b4ae4865eb2c2320d021be534cc57c3 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -939,10 +939,9 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
/* Subgrid tasks: star formation */
- else if (t_type == task_type_star_formation ||
- t_type == task_type_stars_resort) {
+ else if (t_type == task_type_star_formation) {
if (cell_is_active_hydro(t->ci, e)) {
- scheduler_activate(s, t);
+ cell_activate_star_formation_tasks(t->ci, s);
cell_activate_super_spart_drifts(t->ci, s);
}
}
diff --git a/src/hydro/AnarchyDU/hydro.h b/src/hydro/AnarchyDU/hydro.h
index a8b84c46aca9a117efe8173c492bbfd678106a5f..ff8eac12510cce6f5e5a578fbf91479aa67b3f5c 100644
--- a/src/hydro/AnarchyDU/hydro.h
+++ b/src/hydro/AnarchyDU/hydro.h
@@ -398,17 +398,23 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
const struct cosmology *cosmo,
const float u) {
+ /* There is no need to use the floor here as this function is called in the
+ * feedback, so the new value of the internal energy should be strictly
+ * higher than the old value. */
+
p->u = u / cosmo->a_factor_internal_energy;
/* Now recompute the extra quantities */
/* Compute the sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
- const float pressure = hydro_get_comoving_pressure(p);
+ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
/* Update variables. */
p->force.soundspeed = soundspeed;
p->force.pressure = pressure;
+
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
@@ -842,6 +848,9 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ /* Update the signal velocity, if we need to. */
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
@@ -903,10 +912,13 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Compute the new sound speed */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ /* Update signal velocity if we need to */
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h
index 806d4aed204c71d3b9265d1bbea1098fdb83ebb1..ca154374a07a99bcd94cee481d1effce33787d2f 100644
--- a/src/hydro/AnarchyPU/hydro.h
+++ b/src/hydro/AnarchyPU/hydro.h
@@ -219,9 +219,7 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) {
/* Compute the sound speed -- see theory section for justification */
/* IDEAL GAS ONLY -- P-U does not work with generic EoS. */
- const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho);
-
- return square_rooted;
+ return gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
}
/**
@@ -408,15 +406,29 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
const struct cosmology *cosmo,
const float u) {
+ /* Store ratio of new internal energy to old internal energy, as we use this
+ * in the drifting of the pressure. */
+ float internal_energy_ratio = 1.f / p->u;
+
+ /* Update the internal energy */
p->u = u / cosmo->a_factor_internal_energy;
+ internal_energy_ratio *= p->u;
+
+ /* Now we can use this to 'update' the value of the smoothed pressure. To
+ * truly update this variable, we would need another loop over neighbours
+ * using the new internal energies of everyone, but that's not feasible. */
+ p->pressure_bar *= internal_energy_ratio;
/* Now recompute the extra quantities */
/* Compute the sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float soundspeed =
+ gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
/* Update variables. */
p->force.soundspeed = soundspeed;
+
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
@@ -467,10 +479,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
(cosmo->a_factor_sound_speed * p->viscosity.v_sig);
- const float dt_u_change =
- (p->u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->u_dt) : FLT_MAX;
-
- return fminf(dt_cfl, dt_u_change);
+ return dt_cfl;
}
/**
@@ -871,8 +880,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props) {
+ /* Store ratio of new internal energy to old internal energy, as we use this
+ * in the drifting of the pressure. */
+ float internal_energy_ratio = 1.f / p->u;
+
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
+ internal_energy_ratio *= p->u;
+
+ /* Now we can use this to 'update' the value of the smoothed pressure. To
+ * truly update this variable, we would need another loop over neighbours
+ * using the new internal energies of everyone, but that's not feasible. */
+ p->pressure_bar *= internal_energy_ratio;
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
@@ -908,9 +927,12 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
}
/* Compute the new sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
-
+ const float soundspeed =
+ gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
p->force.soundspeed = soundspeed;
+
+ /* Update the signal velocity */
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 469fee137ee624ab649718f25c490f61e35cfb59..5a62770d9c4ddc1381ee9b390c06d0e8e607d32b 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -403,12 +403,14 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
/* Now recompute the extra quantities */
/* Compute the sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
- const float pressure = hydro_get_comoving_pressure(p);
+ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
/* Update variables. */
p->force.soundspeed = soundspeed;
p->force.pressure = pressure;
+
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
@@ -782,10 +784,12 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
/* Compute the sound speed */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
@@ -847,10 +851,12 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Compute the new sound speed */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 3fa6f19cc145890feacbd7284368d5378654bf38..00fc59e70c2f421f39a647d104e1369e6c6f626b 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -403,6 +403,8 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
/* Update variables. */
p->force.P_over_rho2 = P_over_rho2;
p->force.soundspeed = soundspeed;
+
+ p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
}
/**
@@ -650,7 +652,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
p->force.h_dt = 0.0f;
/* Reset maximal signal velocity */
- p->force.v_sig = p->force.soundspeed;
+ p->force.v_sig = 2.f * p->force.soundspeed;
}
/**
@@ -757,6 +759,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Update variables */
p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = P_over_rho2;
+
+ p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
}
/**
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 3d7f43579033afd8f5c29e765b31fee145d9c590..cdfbdbf3bcf3790b8cba84ce74084d3885cf75bf 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -390,10 +390,14 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
/* Now recompute the extra quantities */
/* Compute the sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
/* Update variables. */
+ p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
}
/**
@@ -633,7 +637,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
/* Reset the time derivatives. */
p->u_dt = 0.0f;
p->force.h_dt = 0.0f;
- p->force.v_sig = p->force.soundspeed;
+ p->force.v_sig = 2.f * p->force.soundspeed;
}
/**
@@ -665,6 +669,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
/* Update variables */
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
}
/**
@@ -728,6 +734,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
+
+ p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
}
/**
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index e49445f59edd7a3a09860bfffe81a56d6a05abd4..6ce6c5c526ca016884f13d88044ed3a19248de4c 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -226,6 +226,9 @@ __attribute__((always_inline)) INLINE static void hydro_update_soundspeed(
const float comoving_pressure =
pressure_floor_get_comoving_pressure(p, p->pressure_bar, cosmo);
p->force.soundspeed = gas_soundspeed_from_pressure(p->rho, comoving_pressure);
+
+ /* Also update the signal velocity; this could be a huge change! */
+ p->force.v_sig = max(p->force.v_sig, 2.f * p->force.soundspeed);
}
/**
@@ -424,9 +427,18 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
const struct cosmology *cosmo,
const float u) {
+ /* Store ratio of new internal energy to old internal energy, as we use this
+ * in the drifting of the pressure. */
+ float internal_energy_ratio = 1.f / p->u;
+
+ /* Update the internal energy */
p->u = u / cosmo->a_factor_internal_energy;
+ internal_energy_ratio *= p->u;
- /* Now recompute the extra quantities */
+ /* Now we can use this to 'update' the value of the smoothed pressure. To
+ * truly update this variable, we would need another loop over neighbours
+ * using the new internal energies of everyone, but that's not feasible. */
+ p->pressure_bar *= internal_energy_ratio;
/* Update variables. */
hydro_update_soundspeed(p, cosmo);
@@ -476,10 +488,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
(cosmo->a_factor_sound_speed * p->force.v_sig);
- const float dt_u_change =
- (p->u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->u_dt) : FLT_MAX;
-
- return fminf(dt_cfl, dt_u_change);
+ return dt_cfl;
}
/**
@@ -677,7 +686,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
/* Reset the time derivatives. */
p->u_dt = 0.0f;
p->force.h_dt = 0.0f;
- p->force.v_sig = p->force.soundspeed;
+ p->force.v_sig = 2.f * p->force.soundspeed;
}
/**
@@ -727,8 +736,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props) {
+ /* Store ratio of new internal energy to old internal energy, as we use this
+ * in the drifting of the pressure. */
+ float internal_energy_ratio = 1.f / p->u;
+
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
+ internal_energy_ratio *= p->u;
+
+ /* Now we can use this to 'update' the value of the smoothed pressure. To
+ * truly update this variable, we would need another loop over neighbours
+ * using the new internal energies of everyone, but that's not feasible. */
+ p->pressure_bar *= internal_energy_ratio;
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
index e21750bdf24f19b0f784b8652bcb5a18a6a89e8f..ae5da6bfb7d8495cfb2b5298861aebdad31f7b75 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
@@ -164,9 +164,8 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) {
/* Compute the sound speed -- see theory section for justification */
/* IDEAL GAS ONLY -- P-U does not work with generic EoS. */
- const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho);
- return square_rooted;
+ return gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
}
/**
@@ -407,15 +406,28 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
const struct cosmology *cosmo,
const float u) {
+ /* Store ratio of new internal energy to old internal energy, as we use this
+ * in the drifting of the pressure. */
+ float internal_energy_ratio = 1.f / p->u;
+
+ /* Update the internal energy */
p->u = u / cosmo->a_factor_internal_energy;
+ internal_energy_ratio *= p->u;
+
+ /* Now we can use this to 'update' the value of the smoothed pressure. To
+ * truly update this variable, we would need another loop over neighbours
+ * using the new internal energies of everyone, but that's not feasible. */
+ p->pressure_bar *= internal_energy_ratio;
/* Now recompute the extra quantities */
/* Compute the sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float soundspeed =
+ gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
/* Update variables. */
p->force.soundspeed = soundspeed;
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
@@ -463,10 +475,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
(cosmo->a_factor_sound_speed * p->force.v_sig);
- const float dt_u_change =
- (p->u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->u_dt) : FLT_MAX;
-
- return fminf(dt_cfl, dt_u_change);
+ return dt_cfl;
}
/**
@@ -640,14 +649,14 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Artificial viscosity updates */
- /* TODO: Actually work out why this cosmology factor is correct
- * and update the SPH / cosmology theory documents. */
+ /* We perform all of this in physical space. */
+ const float h_inv_physical = cosmo->a_inv * h_inv;
+ const float soundspeed_physical = cosmo->a_factor_sound_speed * soundspeed;
- /* We divide by a^2 here to make this transform under cosmology the
- * same as the velocity (which in SWIFT has an extra 1/a^2 factor.
- * See the cosmology theory documents for more information. */
+ /* Decay rate */
const float inverse_tau =
- (hydro_props->viscosity.length * cosmo->a2_inv) * soundspeed * h_inv;
+ hydro_props->viscosity.length * soundspeed_physical * h_inv_physical;
+ /* Source term (div v) is already in physical co-ordinates for this scheme */
const float source_term = -1.f * min(p->density.div_v, 0.f);
/* Compute da/dt */
@@ -680,7 +689,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
/* Reset the time derivatives. */
p->u_dt = 0.0f;
p->force.h_dt = 0.0f;
- p->force.v_sig = p->force.soundspeed;
+ p->force.v_sig = 2.f * p->force.soundspeed;
}
/**
@@ -732,8 +741,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
const struct hydro_props *hydro_props,
const struct entropy_floor_properties *floor_props) {
+ /* Store ratio of new internal energy to old internal energy, as we use this
+ * in the drifting of the pressure. */
+ float internal_energy_ratio = 1.f / p->u;
+
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
+ internal_energy_ratio *= p->u;
+
+ /* Now we can use this to 'update' the value of the smoothed pressure. To
+ * truly update this variable, we would need another loop over neighbours
+ * using the new internal energies of everyone, but that's not feasible. */
+ p->pressure_bar *= internal_energy_ratio;
/* Check against entropy floor */
const float floor_A = entropy_floor(p, cosmo, floor_props);
@@ -769,9 +788,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
}
/* Compute the new sound speed */
- const float soundspeed = hydro_get_comoving_soundspeed(p);
+ const float soundspeed =
+ gas_soundspeed_from_pressure(p->rho, p->pressure_bar);
p->force.soundspeed = soundspeed;
+ p->viscosity.v_sig = max(p->viscosity.v_sig, 2.f * soundspeed);
}
/**
diff --git a/src/runner.c b/src/runner.c
index ba9ec9b83c9154c8261f4f559a180da43c154dd3..e1c3b5139d6cad9095f46f10b47584daffc86892 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1216,8 +1216,7 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
* re-compute them. */
if (with_feedback && (c == c->top) &&
(current_stars_count != c->stars.count)) {
- cell_set_flag(c, cell_flag_do_stars_resort);
- cell_clear_stars_sort_flags(c, /*clear_unused_flags=*/0);
+ cell_set_star_resort_flag(c);
}
if (timer) TIMER_TOC(timer_do_star_formation);
@@ -1237,7 +1236,6 @@ void runner_do_stars_resort(struct runner *r, struct cell *c, const int timer) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->nodeID != r->e->nodeID) error("Task must be run locally!");
- if (c->depth != 0) error("Task must be run at the top-level");
#endif
TIMER_TIC;
diff --git a/src/space.c b/src/space.c
index d67ab9a0bf3cc425b8939ccabdd61a4ec3274e97..f6b863709d743091afc56e7c4d0ab58d981c34c8 100644
--- a/src/space.c
+++ b/src/space.c
@@ -92,6 +92,9 @@ int space_extra_gparts = space_extra_gparts_default;
int engine_max_parts_per_ghost = engine_max_parts_per_ghost_default;
int engine_max_sparts_per_ghost = engine_max_sparts_per_ghost_default;
+/*! Maximal depth at which the stars resort task can be pushed */
+int engine_star_resort_task_depth = engine_star_resort_task_depth_default;
+
/*! Expected maximal number of strays received at a rebuild */
int space_expected_max_nr_strays = space_expected_max_nr_strays_default;
#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_CELL_GRAPH)
@@ -5424,6 +5427,18 @@ void space_struct_dump(struct space *s, FILE *stream) {
"space_extra_sparts", "space_extra_sparts");
restart_write_blocks(&space_extra_bparts, sizeof(int), 1, stream,
"space_extra_bparts", "space_extra_bparts");
+ restart_write_blocks(&space_expected_max_nr_strays, sizeof(int), 1, stream,
+ "space_expected_max_nr_strays",
+ "space_expected_max_nr_strays");
+ restart_write_blocks(&engine_max_parts_per_ghost, sizeof(int), 1, stream,
+ "engine_max_parts_per_ghost",
+ "engine_max_parts_per_ghost");
+ restart_write_blocks(&engine_max_sparts_per_ghost, sizeof(int), 1, stream,
+ "engine_max_sparts_per_ghost",
+ "engine_max_sparts_per_ghost");
+ restart_write_blocks(&engine_star_resort_task_depth, sizeof(int), 1, stream,
+ "engine_star_resort_task_depth",
+ "engine_star_resort_task_depth");
/* More things to write. */
if (s->nr_parts > 0) {
@@ -5482,6 +5497,14 @@ void space_struct_restore(struct space *s, FILE *stream) {
"space_extra_sparts");
restart_read_blocks(&space_extra_bparts, sizeof(int), 1, stream, NULL,
"space_extra_bparts");
+ restart_read_blocks(&space_expected_max_nr_strays, sizeof(int), 1, stream,
+ NULL, "space_expected_max_nr_strays");
+ restart_read_blocks(&engine_max_parts_per_ghost, sizeof(int), 1, stream, NULL,
+ "engine_max_parts_per_ghost");
+ restart_read_blocks(&engine_max_sparts_per_ghost, sizeof(int), 1, stream,
+ NULL, "engine_max_sparts_per_ghost");
+ restart_read_blocks(&engine_star_resort_task_depth, sizeof(int), 1, stream,
+ NULL, "engine_star_resort_task_depth");
/* Things that should be reconstructed in a rebuild. */
s->cells_top = NULL;
diff --git a/src/timers.c b/src/timers.c
index 8c193f9f2e15746bd551f8a56edafe7c80753334..cc4ebc90d491832178f1e48e96960f489c73f4c1 100644
--- a/src/timers.c
+++ b/src/timers.c
@@ -56,6 +56,9 @@ const char* timers_names[timer_count] = {
"doself_limiter",
"doself_stars_density",
"doself_stars_feedback",
+ "doself_bh_density",
+ "doself_bh_swallow",
+ "doself_bh_feedback",
"doself_grav_pp",
"dopair_density",
"dopair_gradient",
@@ -63,6 +66,9 @@ const char* timers_names[timer_count] = {
"dopair_limiter",
"dopair_stars_density",
"dopair_stars_feedback",
+ "dopair_bh_density",
+ "dopair_bh_swallow",
+ "dopair_bh_feedback",
"dopair_grav_mm",
"dopair_grav_pp",
"dograv_external",
@@ -76,6 +82,9 @@ const char* timers_names[timer_count] = {
"dosub_self_limiter",
"dosub_self_stars_density",
"dosub_self_stars_feedback",
+ "dosub_self_bh_density",
+ "dosub_self_bh_swallow",
+ "dosub_self_bh_feedback",
"dosub_self_grav",
"dosub_pair_density",
"dosub_pair_gradient",
@@ -83,6 +92,9 @@ const char* timers_names[timer_count] = {
"dosub_pair_limiter",
"dosub_pair_stars_density",
"dosub_pair_stars_feedback",
+ "dosub_pair_bh_density",
+ "dosub_pair_bh_swallow",
+ "dosub_pair_bh_feedback",
"dosub_pair_grav",
"doself_subset",
"dopair_subset",
@@ -91,9 +103,11 @@ const char* timers_names[timer_count] = {
"do_ghost",
"do_extra_ghost",
"do_stars_ghost",
+ "do_black_holes_ghost",
"dorecv_part",
"dorecv_gpart",
"dorecv_spart",
+ "dorecv_bpart",
"do_limiter",
"do_cooling",
"do_star_formation",
@@ -106,6 +120,7 @@ const char* timers_names[timer_count] = {
"step",
"logger",
"do_stars_sort",
+ "do_stars_resort",
"fof_self",
"fof_pair",
};
diff --git a/src/timers.h b/src/timers.h
index e36aea2c6269ee07f813fc69764ed14e4d2cc550..96742567f549c156944333fafd4342dba14a4bee 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -57,6 +57,9 @@ enum {
timer_doself_limiter,
timer_doself_stars_density,
timer_doself_stars_feedback,
+ timer_doself_bh_density,
+ timer_doself_bh_swallow,
+ timer_doself_bh_feedback,
timer_doself_grav_pp,
timer_dopair_density,
timer_dopair_gradient,
@@ -64,6 +67,9 @@ enum {
timer_dopair_limiter,
timer_dopair_stars_density,
timer_dopair_stars_feedback,
+ timer_dopair_bh_density,
+ timer_dopair_bh_swallow,
+ timer_dopair_bh_feedback,
timer_dopair_grav_mm,
timer_dopair_grav_pp,
timer_dograv_external,
@@ -77,6 +83,9 @@ enum {
timer_dosub_self_limiter,
timer_dosub_self_stars_density,
timer_dosub_self_stars_feedback,
+ timer_dosub_self_bh_density,
+ timer_dosub_self_bh_swallow,
+ timer_dosub_self_bh_feedback,
timer_dosub_self_grav,
timer_dosub_pair_density,
timer_dosub_pair_gradient,
@@ -84,6 +93,9 @@ enum {
timer_dosub_pair_limiter,
timer_dosub_pair_stars_density,
timer_dosub_pair_stars_feedback,
+ timer_dosub_pair_bh_density,
+ timer_dosub_pair_bh_swallow,
+ timer_dosub_pair_bh_feedback,
timer_dosub_pair_grav,
timer_doself_subset,
timer_dopair_subset,
@@ -92,9 +104,11 @@ enum {
timer_do_ghost,
timer_do_extra_ghost,
timer_do_stars_ghost,
+ timer_do_black_holes_ghost,
timer_dorecv_part,
timer_dorecv_gpart,
timer_dorecv_spart,
+ timer_dorecv_bpart,
timer_do_limiter,
timer_do_cooling,
timer_do_star_formation,