diff --git a/examples/main.c b/examples/main.c
index 690cbba113e0f90037b57d80bc4f3245559efa03..37d9ff3346d94f8a2588f5bdfabf7406be557d26 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -798,8 +798,8 @@ int main(int argc, char *argv[]) {
/* Initialize the engine with the space and policies. */
if (myrank == 0) clocks_gettime(&tic);
- engine_init(&e, &s, params, N_total[0], N_total[1], engine_policies,
- talking, &reparttype, &us, &prog_const, &cosmo,
+ engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2],
+ engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
&hydro_properties, &gravity_properties, &potential,
&cooling_func, &chemistry, &sourceterms);
engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
diff --git a/src/engine.c b/src/engine.c
index bb9a478b4933106c3753a31388bdf78117366147..a3e643fea695ccc33c6b44159569a9995e4fd5ef 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -5361,7 +5361,8 @@ void engine_unpin() {
* @param s The #space in which this #runner will run.
* @param params The parsed parameter file.
* @param Ngas total number of gas particles in the simulation.
- * @param Ndm total number of gravity particles in the simulation.
+ * @param Ngparts total number of gravity particles in the simulation.
+ * @param Nstars total number of star particles in the simulation.
* @param policy The queuing policy to use.
* @param verbose Is this #engine talkative ?
* @param reparttype What type of repartition algorithm are we using ?
@@ -5377,7 +5378,7 @@ void engine_unpin() {
*/
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, long long Ngas,
- long long Ndm, int policy, int verbose,
+ long long Ngparts, long long Nstars, int policy, int verbose,
struct repartition *reparttype,
const struct unit_system *internal_units,
const struct phys_const *physical_constants,
@@ -5396,7 +5397,8 @@ void engine_init(struct engine *e, struct space *s,
e->policy = policy;
e->step = 0;
e->total_nr_parts = Ngas;
- e->total_nr_gparts = Ndm;
+ e->total_nr_gparts = Ngparts;
+ e->total_nr_sparts = Nstars;
e->proxy_ind = NULL;
e->nr_proxies = 0;
e->reparttype = reparttype;
@@ -6136,7 +6138,74 @@ void engine_compute_next_statistics_time(struct engine *e) {
*/
void engine_recompute_displacement_constraint(struct engine *e) {
- message("max_dt_RMS_displacement = %e", e->dt_max_RMS_displacement);
+ const struct cosmology *cosmo = e->cosmology;
+ const float Om = cosmo->Omega_m;
+ const float Ob = cosmo->Omega_b;
+ const float rho_crit = cosmo->critical_density;
+ const float a = cosmo->a;
+
+ /* Start by reducing the minimal mass of each particle type */
+ float min_mass[swift_type_count] = {e->s->min_part_mass,
+ e->s->min_gpart_mass,
+ FLT_MAX,
+ FLT_MAX,
+ e->s->min_spart_mass,
+ FLT_MAX};
+#ifdef WITH_MPI
+ MPI_Allreduce(MPI_IN_PLACE, min_mass, swift_type_count, MPI_FLOAT, MPI_MIN,
+ MPI_COMM_WORLD);
+#endif
+
+ /* Do the same for the velocity norm sum */
+ float vel_norm[swift_type_count] = {e->s->sum_part_vel_norm,
+ e->s->sum_gpart_vel_norm,
+ 0.f,
+ 0.f,
+ e->s->sum_spart_vel_norm,
+ 0.f};
+#ifdef WITH_MPI
+ MPI_Allreduce(MPI_IN_PLACE, vel_norm, swift_type_count, MPI_FLOAT, MPI_SUM,
+ MPI_COMM_WORLD);
+#endif
+
+ /* Get the counts of each particle types */
+ const long long total_nr_dm_gparts =
+ e->total_nr_gparts - e->total_nr_parts - e->total_nr_sparts;
+ float count_parts[swift_type_count] = {
+ e->total_nr_parts, total_nr_dm_gparts, 0.f, 0.f, e->total_nr_sparts, 0.f};
+
+ /* Minimal mass for the two species */
+ const float min_mass_dm = min_mass[1];
+ const float min_mass_b = min(min_mass[0], min_mass[4]);
+
+ /* Inter-particle sepration for the two species */
+ const float d_dm = cbrtf(min_mass_dm / ((Om - Ob) * rho_crit));
+ const float d_b = cbrtf(min_mass_b / (Ob * rho_crit));
+
+ /* Count of particles for the two species */
+ const float N_dm = count_parts[1];
+ const float N_b = count_parts[0] + count_parts[4];
+
+ /* Peculiar motion norm for the two species */
+ const float vel_norm_dm = vel_norm[1];
+ const float vel_norm_b = vel_norm[0] + vel_norm[4];
+
+ /* RMS peculiar motion for the two species */
+ const float rms_vel_dm = (N_dm > 0.f) ? (vel_norm_dm / N_dm) : 0.f;
+ const float rms_vel_b = (N_b > 0.f) ? (vel_norm_b / N_b) : 0.f;
+
+ /* Time-step based on maximum displacement */
+ const float dt_dm = a * a * d_dm / (sqrtf(rms_vel_dm) + FLT_MIN);
+ const float dt_b = a * a * d_b / (sqrtf(rms_vel_b) + FLT_MIN);
+
+ /* Use the minimum */
+ const float dt = min(dt_dm, dt_b);
+
+ /* Apply the dimensionless factor */
+ e->dt_max_RMS_displacement = dt * e->max_RMS_displacement_factor;
+
+ if (e->verbose)
+ message("max_dt_RMS_displacement = %e", e->dt_max_RMS_displacement);
}
/**
diff --git a/src/engine.h b/src/engine.h
index a6707e0426d485a61464c440915c94d6e74a580b..db3db3158748414fb9b048feabcc2695ae0873a5 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -192,7 +192,7 @@ struct engine {
int step_props;
/* Total numbers of particles in the system. */
- long long total_nr_parts, total_nr_gparts;
+ long long total_nr_parts, total_nr_gparts, total_nr_sparts;
/* The internal system of units */
const struct unit_system *internal_units;
@@ -350,7 +350,7 @@ void engine_print_stats(struct engine *e);
void engine_dump_snapshot(struct engine *e);
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, long long Ngas,
- long long Ndm, int policy, int verbose,
+ long long Ngparts, long long Nstars, int policy, int verbose,
struct repartition *reparttype,
const struct unit_system *internal_units,
const struct phys_const *physical_constants,