adding engine unskip stuff. Still in debug mode, syncing for comparisons

src/cell.c

@@ -3838,6 +3838,19 @@ void cell_activate_subcell_external_grav_tasks(struct cell *ci,
   }
 }
 
+/**
+ * @brief Traverse a sub-cell task and activate the radiative transfer tasks
+ *
+ * @param ci The first #cell we recurse in.
+ * @param cj The second #cell we recurse in.
+ * @param s The task #scheduler.
+ */
+void cell_activate_subcell_rt_tasks(struct cell *ci, struct cell *cj,
+                                    struct scheduler *s) {
+  /* TODO: everything */
+}
+
+
 /**
  * @brief Un-skips all the hydro tasks associated with a given cell and checks
  * if the space needs to be rebuilt.
@@ -4819,6 +4832,60 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
   return rebuild;
 }
 
+/**
+ * @brief Un-skips all the RT tasks associated with a given cell and checks
+ * if the space needs to be rebuilt.
+ *
+ * @param c the #cell.
+ * @param s the #scheduler.
+ *
+ * @return 1 If the space needs rebuilding. 0 otherwise.
+ */
+int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s) {
+  struct engine *e = s->space->e;
+  const int nodeID = e->nodeID;
+
+  int rebuild = 0;
+  int counter = 0;
+
+  for (struct link *l = c->hydro.rt_inject; l != NULL; l = l->next) {
+    counter ++;
+    struct task *t = l->t;
+    struct cell *ci = t->ci;
+    struct cell *cj = t->cj;
+    const int ci_active = cell_is_active_hydro(ci, e);
+    const int cj_active = (cj != NULL) ? cell_is_active_hydro(cj, e) : 0;
+#ifdef WITH_MPI
+    const int ci_nodeID = ci->nodeID;
+    const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
+#else
+    const int ci_nodeID = nodeID;
+    const int cj_nodeID = nodeID;
+#endif
+
+    /* Only activate tasks that involve a local active cell. */
+    if ((ci_active && ci_nodeID == nodeID) ||
+        (cj_active && cj_nodeID == nodeID)) {
+
+      message("ACTIVATING %d", counter);
+      scheduler_activate(s, t);
+
+      if (t->type == task_type_sub_self) {
+        message("trying to access subcell rt tasks activation");
+        cell_activate_subcell_rt_tasks(ci, NULL, s);
+      }
+
+      else if (t->type == task_type_sub_pair) {
+        cell_activate_subcell_rt_tasks(ci, cj, s);
+        message("trying to access subcell rt tasks activation");
+      }
+    }
+  }
+
+  return rebuild;
+}
+
+
 /**
  * @brief Set the super-cell pointers for all cells in a hierarchy.
  *

src/cell.h

@@ -978,6 +978,7 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s);
 int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s,
                             const int with_star_formation);
 int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s);
+int cell_unskip_rt_tasks(struct cell *c, struct scheduler *s);
 int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s);
 int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s);
 void cell_drift_part(struct cell *c, const struct engine *e, int force);
@@ -1007,6 +1008,8 @@ void cell_activate_subcell_black_holes_tasks(struct cell *ci, struct cell *cj,
                                              const int with_timestep_sync);
 void cell_activate_subcell_external_grav_tasks(struct cell *ci,
                                                struct scheduler *s);
+void cell_activate_subcell_rt_tasks(struct cell *ci, struct cell *cj,
+                                    struct scheduler *s);
 void cell_activate_super_spart_drifts(struct cell *c, struct scheduler *s);
 void cell_activate_drift_part(struct cell *c, struct scheduler *s);
 void cell_activate_drift_gpart(struct cell *c, struct scheduler *s);

src/engine_maketasks.c

@@ -1777,7 +1777,8 @@ void engine_link_gravity_tasks(struct engine *e) {
           scheduler_addunlock(sched, t, cj_parent->grav.down_in);
         }
       }
-    }
+
+}
 
     /* Otherwise M-M interaction? */
     else if (t_type == task_type_grav_mm) {

src/engine_unskip.c

@@ -46,6 +46,7 @@ enum task_broad_types {
   task_broad_types_stars,
   task_broad_types_sinks,
   task_broad_types_black_holes,
+  task_broad_types_rt, 
   task_broad_types_count,
 };
 
@@ -238,6 +239,39 @@ static void engine_do_unskip_gravity(struct cell *c, struct engine *e) {
   cell_unskip_gravity_tasks(c, &e->sched);
 }
 
+/**
+ * @brief Unskip any radiative transfer tasks associated with active cells.
+ *
+ * @param c The cell.
+ * @param e The engine.
+ */
+static void engine_do_unskip_rt(struct cell *c, struct engine *e) {
+
+  /* Early abort (are we below the level where tasks are)? */
+  if (!cell_get_flag(c, cell_flag_has_tasks)) return;
+
+  /* Ignore empty cells. */
+  if (c->hydro.count == 0) return;
+
+  /* Skip inactive cells. */
+  if (!cell_is_active_hydro(c, e)) return;
+
+  /* Recurse */
+  if (c->split) {
+    for (int k = 0; k < 8; k++) {
+      if (c->progeny[k] != NULL) {
+        struct cell *cp = c->progeny[k];
+        engine_do_unskip_rt(cp, e);
+      }
+    }
+  }
+
+  /* Unskip any active tasks. */
+  const int forcerebuild = cell_unskip_rt_tasks(c, &e->sched);
+  if (forcerebuild) atomic_inc(&e->forcerebuild);
+}
+
+
 /**
  * @brief Mapper function to unskip active tasks.
  *
@@ -316,6 +350,13 @@ void engine_do_unskip_mapper(void *map_data, int num_elements,
 #endif
         engine_do_unskip_black_holes(c, e);
         break;
+      case task_broad_types_rt:
+#ifdef SWIFT_DEBUG_CHECKS
+        if (!(e->policy & engine_policy_rt))
+          error("Trying to unskip radiative transfer tasks in a non-rt run!");
+#endif
+        engine_do_unskip_rt(c, e);
+        break;
       default:
 #ifdef SWIFT_DEBUG_CHECKS
         error("Invalid broad task type!");
@@ -343,6 +384,7 @@ void engine_unskip(struct engine *e) {
   const int with_sinks = e->policy & engine_policy_sinks;
   const int with_feedback = e->policy & engine_policy_feedback;
   const int with_black_holes = e->policy & engine_policy_black_holes;
+  const int with_rt = e->policy & engine_policy_rt;
 
 #ifdef WITH_PROFILER
   static int count = 0;
@@ -396,6 +438,10 @@ void engine_unskip(struct engine *e) {
     data.task_types[multiplier] = task_broad_types_black_holes;
     multiplier++;
   }
+  if (with_rt) {
+    data.task_types[multiplier] = task_broad_types_rt;
+    multiplier++;
+  }
 
   /* Should we duplicate the list of active cells to better parallelise the
      unskip over the threads ? */

src/rt/debug/rt.h

@@ -33,6 +33,8 @@ __attribute__((always_inline)) INLINE static void rt_first_init_xpart(
   xp->rt_data.iact_stars = 0;
   xp->rt_data.calls_tot = 0;
   xp->rt_data.calls_per_step = 0;
+  xp->rt_data.calls_self = 0;
+  xp->rt_data.calls_pair = 0;
 }
 
 /**
@@ -43,6 +45,8 @@ __attribute__((always_inline)) INLINE static void rt_init_xpart(
 
   xp->rt_data.iact_stars = 0;
   xp->rt_data.calls_per_step = 0;
+  xp->rt_data.calls_self = 0;
+  xp->rt_data.calls_pair = 0;
 }
 
 /**
@@ -54,6 +58,8 @@ __attribute__((always_inline)) INLINE static void rt_first_init_spart(
   sp->rt_data.iact_hydro = 0;
   sp->rt_data.calls_tot = 0;
   sp->rt_data.calls_per_step = 0;
+  sp->rt_data.calls_self = 0;
+  sp->rt_data.calls_pair = 0;
 }
 
 /**
@@ -64,5 +70,7 @@ __attribute__((always_inline)) INLINE static void rt_init_spart(
 
   sp->rt_data.iact_hydro = 0;
   sp->rt_data.calls_per_step = 0;
+  sp->rt_data.calls_self = 0;
+  sp->rt_data.calls_pair = 0;
 }
 #endif /* SWIFT_RT_DEBUG_H */

src/rt/debug/rt_iact.h

@@ -48,7 +48,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_rt_inject(
 
   pd->iact_stars += 1;
   pd->calls_tot += 1;
-  sd->calls_per_step += 1;
+  pd->calls_per_step += 1;
+
+  if (r2 > 0.f){
+    sd->calls_self += 1;
+    pd->calls_self += 1;
+  } else {
+    sd->calls_pair += 1;
+    pd->calls_pair += 1;
+  }
 }
 
 #endif /* SWIFT_RT_IACT_DEBUG_H */

src/rt/debug/rt_io.h

@@ -34,7 +34,7 @@ INLINE static void rt_write_xparts(const struct xpart* xparts,
                                    struct io_props* list, int* num_fields) {
 
   list += *num_fields;
-  *num_fields += 3;
+  *num_fields += 5;
 
   list[0] = io_make_output_field("RT_star_iact", INT, 1, UNIT_CONV_NO_UNITS,
                                  1.f, xparts, rt_data.iact_stars,
@@ -48,6 +48,14 @@ INLINE static void rt_write_xparts(const struct xpart* xparts,
       "RT_calls_this_step", INT, 1, UNIT_CONV_NO_UNITS, 1.f, xparts,
       rt_data.calls_per_step,
       "number of calls to this particle during one time step");
+  list[3] = io_make_output_field(
+      "RT_calls_self", INT, 1, UNIT_CONV_NO_UNITS, 1.f, xparts,
+      rt_data.calls_self,
+      "number of calls to this particle during one time step in self task");
+  list[4] = io_make_output_field(
+      "RT_calls_pair", INT, 1, UNIT_CONV_NO_UNITS, 1.f, xparts,
+      rt_data.calls_pair,
+      "number of calls to this particle during one time step in self task");
 }
 
 /**
@@ -58,7 +66,7 @@ INLINE static void rt_write_stars(const struct spart* sparts,
                                   struct io_props* list, int* num_fields) {
 
   list += *num_fields;
-  *num_fields += 3;
+  *num_fields += 5;
 
   list[0] = io_make_output_field("RT_hydro_iact", INT, 1, UNIT_CONV_NO_UNITS,
                                  1.f, sparts, rt_data.iact_hydro,
@@ -72,5 +80,13 @@ INLINE static void rt_write_stars(const struct spart* sparts,
       "RT_calls_this_step", INT, 1, UNIT_CONV_NO_UNITS, 1.f, sparts,
       rt_data.calls_per_step,
       "number of calls to this particle during one time step");
+  list[3] = io_make_output_field(
+      "RT_calls_self", INT, 1, UNIT_CONV_NO_UNITS, 1.f, sparts,
+      rt_data.calls_self,
+      "number of calls to this particle during one time step in self task");
+  list[4] = io_make_output_field(
+      "RT_calls_pair", INT, 1, UNIT_CONV_NO_UNITS, 1.f, sparts,
+      rt_data.calls_pair,
+      "number of calls to this particle during one time step in self task");
 }
 #endif /* SWIFT_RT_IO_DEBUG_H */

src/rt/debug/rt_struct.h

@@ -28,12 +28,16 @@ struct rt_xpart_data {
   int iact_stars; /* how many stars this particle interacted with */
   int calls_tot;  /* total number of calls to this particle during entire run */
   int calls_per_step; /* calls per time step to this particle */
+  int calls_self;
+  int calls_pair;
 };
 
 struct rt_spart_data {
   int iact_hydro; /* how many hydro particles this particle interacted with */
   int calls_tot;  /* total number of calls to this particle during entire run */
   int calls_per_step; /* calls per time step to this particle */
+  int calls_self;
+  int calls_pair;
 };
 
 #endif /* SWIFT_RT_STRUCT_DEBUG_H */

src/runner_doiact_functions_rt.h

@@ -38,6 +38,8 @@ void DOSELF1_RT(struct runner *r, struct cell *c, int timer) {
 
   TIMER_TIC;
 
+  message("Called self for step %d", r->e->step);
+
   /* Anything to do here? */
   if (c->hydro.count == 0 || c->stars.count == 0) return;
 
@@ -89,6 +91,7 @@ void DOSELF1_RT(struct runner *r, struct cell *c, int timer) {
  * @param cj the second cell, where we take hydro particles from
  */
 void DOPAIR1_NONSYM_RT(struct runner *r, struct cell *ci, struct cell *cj) {
+  message("Called pair for step %d", r->e->step);
 
   const struct engine *e = r->e;
 
@@ -137,7 +140,7 @@ void DOPAIR1_NONSYM_RT(struct runner *r, struct cell *ci, struct cell *cj) {
       float dx[3] = {six[0] - pjx[0], six[1] - pjx[1], six[2] - pjx[2]};
       const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
 
-      if (r2 < hjg2) IACT_RT(r2, dx, hi, hj, si, xpj);
+      if (r2 < hjg2) IACT_RT(-1.f, dx, hi, hj, si, xpj);
 
     } /* loop over the parts in cj. */
   }   /* loop over the parts in ci. */