diff --git a/src/timestep.h b/src/timestep.h index ac44ff3d313963b776fb48cfcb22b28a222144c5..6f41d60ed615749a053e458f4a8cf3371513e3f1 100644 --- a/src/timestep.h +++ b/src/timestep.h @@ -175,6 +175,21 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep( new_dt_grav = min(new_dt_self_grav, new_dt_ext_grav); } + /* Compute the next timestep (force condition) */ + float new_dt_force = FLT_MAX; + if (e->policy & engine_policy_external_gravity){ + + float a[3]; + a[0] = p->a_hydro[0] + p->gpart->a_grav[0]; + a[1] = p->a_hydro[1] + p->gpart->a_grav[1]; + a[2] = p->a_hydro[2] + p->gpart->a_grav[2]; + + const float norm_a = sqrtf(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]); + + new_dt_force = norm_a != 0.0f ? 0.25f * sqrtf(p->h / norm_a) : FLT_MAX; + + }; + /* Compute the next timestep (forcing terms condition) */ const float new_dt_forcing = forcing_terms_timestep( e->time, e->forcing_terms, e->physical_constants, p, xp); @@ -196,6 +211,11 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep( new_dt = min(new_dt, dt_h_change); + /* Apply the force criterion */ + new_dt = min(new_dt, new_dt_force); + + if (new_dt == new_dt_force) printf("The force criterion was used ! \n"); + /* Apply the maximal displacement constraint (FLT_MAX if non-cosmological)*/ new_dt = min(new_dt, e->dt_max_RMS_displacement);