From 6a93d071b6ecebe75bc9af021406b3524fced0b2 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 9 May 2018 19:02:13 +0200 Subject: [PATCH] Use the values gathered to compute the time-step based on the RMS peculiar velocity in cosmological runs. --- examples/main.c | 4 +-- src/engine.c | 77 ++++++++++++++++++++++++++++++++++++++++++++++--- src/engine.h | 4 +-- 3 files changed, 77 insertions(+), 8 deletions(-) diff --git a/examples/main.c b/examples/main.c index 690cbba113..37d9ff3346 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 bb9a478b49..a3e643fea6 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 a6707e0426..db3db31587 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, -- GitLab