Commit 6a93d071 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Use the values gathered to compute the time-step based on the RMS peculiar...

Use the values gathered to compute the time-step based on the RMS peculiar velocity in cosmological runs.
parent 3379d0e8
......@@ -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,
......
......@@ -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);
}
/**
......
......@@ -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,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment