Merge branch 'debugging_rt' into 'master'

Prevent overflows in time step calculations with RT

See merge request !1819

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.
  *

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 */

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

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);
     }
   }