diff --git a/src/rt.h b/src/rt.h index 21482ce3618876c785e2a04926feab2dfb204151..43abd3d0b23d2eec305238ebd88215e3df9b8f8f 100644 --- a/src/rt.h +++ b/src/rt.h @@ -21,7 +21,8 @@ /** * @file src/rt.h - * @brief Branches between the different radiative transfer schemes. + * @brief Branches between the different radiative transfer schemes. Also + * contains some globally valid functions related to time bin data. */ /* Config parameters. */ @@ -44,6 +45,19 @@ #error "Invalid choice of radiation scheme" #endif +/** + * @brief Initialise RT time step data. This struct should be hidden from users, + * so we do it for all schemes here. + * + * @param p The #part. + */ +__attribute__((always_inline)) INLINE static void rt_first_init_timestep_data( + struct part *restrict p) { + + p->rt_time_data.min_ngb_time_bin = num_time_bins + 1; + p->rt_time_data.time_bin = 0; +} + /** * @brief Prepare the rt *time step* quantities for a *hydro force* calculation. * diff --git a/src/runner_time_integration.c b/src/runner_time_integration.c index 152fc8ae342e6b2bdf27b486443e0eb0f07cf152..fd0f0d10a0ab64f877ebef53c7e0e346370802b6 100644 --- a/src/runner_time_integration.c +++ b/src/runner_time_integration.c @@ -740,6 +740,10 @@ void runner_do_timestep(struct runner *r, struct cell *c, const int timer) { * debugging/consistency check. */ rt_debugging_check_timestep(p, &ti_rt_new_step, &ti_new_step, e->max_nr_rt_subcycles, e->time_base); + + if (ti_rt_new_step <= 0LL) + error("Got integer time step <= 0? %lld %lld", + get_part_rt_timestep(p, xp, e), ti_new_step); #endif /* Update particle */ diff --git a/src/space_first_init.c b/src/space_first_init.c index ad98957bfefed0ea209296e7b542fc89ac1b5123..5869ac629318ee9cab825efc2f9337201d796346 100644 --- a/src/space_first_init.c +++ b/src/space_first_init.c @@ -144,6 +144,7 @@ void space_first_init_parts_mapper(void *restrict map_data, int count, particle_splitting_mark_part_as_not_split(&xp[k].split_data, p[k].id); /* And the radiative transfer */ + rt_first_init_timestep_data(&p[k]); rt_first_init_part(&p[k], cosmo, rt_props); #ifdef SWIFT_DEBUG_CHECKS diff --git a/src/timestep.h b/src/timestep.h index a77a2ebdc46857c3fa3d81af63675971489baf88..c09895bac00e2793b621d68697d7d51f33839b7d 100644 --- a/src/timestep.h +++ b/src/timestep.h @@ -218,7 +218,21 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep( /* enforce dt_hydro <= nsubcycles * dt_rt. The rare case where * new_dti_rt > new_dti will be handled in the parent function * that calls this one. */ - const integertime_t max_subcycles = max(e->max_nr_rt_subcycles, 1); + integertime_t max_subcycles = max(e->max_nr_rt_subcycles, 1); + if (max_nr_timesteps / max_subcycles < new_dti_rt) { + /* multiplication new_dti_rt * max_subcycles would overflow. This can + * happen in rare cases, especially if the total physical time the + * simulation should cover is small. So limit max_subcycles to a + * reasonable value. + * First find an integer guess for the maximal permissible number + * of sub-cycles. Then find highest power-of-two below that guess. + * Divide the guess by a factor of 2 to simplify the subsequent while + * loop. The max() is there to prevent bad things happening. */ + const integertime_t max_subcycles_guess = + max(1LL, max_nr_timesteps / (new_dti_rt * 2LL)); + max_subcycles = 1LL; + while (max_subcycles_guess > max_subcycles) max_subcycles *= 2LL; + } new_dti = min(new_dti, new_dti_rt * max_subcycles); } }