Commit d7fd5959 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Also use a time-step condition based on the top-level mesh cell size.

parent 6a93d071
......@@ -59,11 +59,12 @@ Scheduler:
# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 1. # The end time of the simulation (in internal units).
dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 1. # The end time of the simulation (in internal units).
dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
max_dt_RMS_factor: 0.25 # (Optional) Dimensionless factor for the maximal displacement allowed based on the RMS velocities.
# Parameters governing the snapshots
Snapshots:
basename: output # Common part of the name of output files
......
......@@ -5429,7 +5429,7 @@ void engine_init(struct engine *e, struct space *s,
e->dt_max = parser_get_param_double(params, "TimeIntegration:dt_max");
e->dt_max_RMS_displacement = FLT_MAX;
e->max_RMS_displacement_factor = parser_get_opt_param_double(
params, "TimeIntegration:max_dt_RMS_factor", 1.);
params, "TimeIntegration:max_dt_RMS_factor", 0.25);
e->a_first_statistics =
parser_get_opt_param_double(params, "Statistics:scale_factor_first", 0.1);
e->time_first_statistics =
......@@ -6138,6 +6138,7 @@ void engine_compute_next_statistics_time(struct engine *e) {
*/
void engine_recompute_displacement_constraint(struct engine *e) {
/* Get the cosmological information */
const struct cosmology *cosmo = e->cosmology;
const float Om = cosmo->Omega_m;
const float Ob = cosmo->Omega_b;
......@@ -6151,6 +6152,16 @@ void engine_recompute_displacement_constraint(struct engine *e) {
FLT_MAX,
e->s->min_spart_mass,
FLT_MAX};
#ifdef SWIFT_DEBUG_CHECKS
/* Check that the minimal mass collection worked */
float min_part_mass_check = FLT_MAX;
for (size_t i = 0; i < e->s->nr_parts; ++i)
min_part_mass_check =
min(min_part_mass_check, hydro_get_mass(&e->s->parts[i]));
if (min_part_mass_check != min_mass[swift_type_gas])
error("Error collecting minimal mass of gas particles.");
#endif
#ifdef WITH_MPI
MPI_Allreduce(MPI_IN_PLACE, min_mass, swift_type_count, MPI_FLOAT, MPI_MIN,
MPI_COMM_WORLD);
......@@ -6174,14 +6185,6 @@ void engine_recompute_displacement_constraint(struct engine *e) {
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];
......@@ -6190,13 +6193,43 @@ void engine_recompute_displacement_constraint(struct engine *e) {
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;
/* Mesh forces smoothing scale */
const float a_smooth =
e->gravity_properties->a_smooth * e->s->dim[0] / e->s->cdim[0];
float dt_dm = FLT_MAX, dt_b = FLT_MAX;
/* DM case */
if (N_dm > 0.f) {
/* Minimal mass for the DM */
const float min_mass_dm = min_mass[1];
/* Inter-particle sepration for the DM */
const float d_dm = cbrtf(min_mass_dm / ((Om - Ob) * rho_crit));
/* 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);
/* RMS peculiar motion for the DM */
const float rms_vel_dm = vel_norm_dm / N_dm;
/* Time-step based on maximum displacement */
dt_dm = a * a * min(a_smooth, d_dm) / sqrtf(rms_vel_dm);
}
/* Baryon case */
if (N_b > 0.f) {
/* Minimal mass for the baryons */
const float min_mass_b = min(min_mass[0], min_mass[4]);
/* Inter-particle sepration for the baryons */
const float d_b = cbrtf(min_mass_b / (Ob * rho_crit));
/* RMS peculiar motion for the baryons */
const float rms_vel_b = vel_norm_b / N_b;
/* Time-step based on maximum displacement */
dt_b = a * a * min(a_smooth, d_b) / sqrtf(rms_vel_b);
}
/* Use the minimum */
const float dt = min(dt_dm, dt_b);
......
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