Merge branch 'rt-injection-task' into 'master'

RT injection task

See merge request !1155

examples/main.c

@@ -579,7 +579,27 @@ int main(int argc, char *argv[]) {
         "Error: Cannot use radiative transfer without stars, --stars must be "
         "chosen\n");
   }
+  /* Temporary, this dependency will be removed later */
+  if (with_rt && !with_feedback) {
+    error(
+        "Error: Cannot use radiative transfer without feedback for now,"
+        " --feedback must be chosen\n");
+  }
+
+#ifndef GADGET2_SPH
+  /* Temporary, this dependency will be removed later */
+  error("Error: Cannot use radiative transfer without gadget2-sph for now\n");
+#endif
+#ifndef STARS_GEAR
+  /* Temporary, this dependency will be removed later */
+  error(
+      "Error: Cannot use radiative transfer without GEAR star model for now\n");
+#endif
+#ifdef WITH_MPI
+  /* Temporary, this will be removed in due time */
+  error("Error: Cannot use radiative transfer with MPI\n");
 #endif
+#endif /* idfef RT_NONE */
 
 /* Let's pin the main thread, now we know if affinity will be used. */
 #if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)

src/Makefile.am

@@ -102,7 +102,8 @@ endif
 AM_SOURCES = space.c runner_main.c runner_doiact_hydro.c runner_doiact_limiter.c \
     runner_doiact_stars.c runner_doiact_black_holes.c runner_ghost.c runner_recv.c \
     runner_sort.c runner_drift.c runner_black_holes.c runner_time_integration.c \
-    runner_doiact_hydro_vec.c runner_others.c\
+    runner_doiact_hydro_vec.c runner_others.c \
+	runner_doiact_rt.c\
     queue.c task.c cell.c engine.c engine_maketasks.c engine_split_particles.c \
     engine_marktasks.c engine_drift.c engine_unskip.c engine_collect_end_of_step.c \
     engine_redistribute.c engine_fof.c serial_io.c timers.c debug.c scheduler.c \
@@ -127,6 +128,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
 	         runner_doiact_nosort.h runner_doiact_hydro.h runner_doiact_stars.h runner_doiact_black_holes.h runner_doiact_grav.h \
                  runner_doiact_functions_hydro.h runner_doiact_functions_stars.h runner_doiact_functions_black_holes.h \
 		 runner_doiact_functions_limiter.h runner_doiact_limiter.h units.h intrinsics.h minmax.h \
+		 	runner_doiact_rt.h runner_doiact_functions_rt.h \
                  kick.h timestep.h drift.h adiabatic_index.h io_properties.h dimension.h part_type.h periodic.h memswap.h \
                  timestep_limiter.h timestep_limiter_iact.h timestep_sync.h timestep_sync_part.h timestep_limiter_struct.h \
      dump.h logger.h sign.h logger_io.h hashmap.h gravity.h gravity_io.h gravity_logger.h  gravity_cache.h output_options.h\
@@ -208,8 +210,21 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
 	         riemann/riemann_hllc.h riemann/riemann_trrs.h \
 		 riemann/riemann_exact.h riemann/riemann_vacuum.h \
                  riemann/riemann_checks.h \
-		 rt.h rt/none/rt.h rt/debug/rt.h rt/M1closure/rt.h \
-	 	 stars.h stars_io.h stars_logger.h \
+		 rt.h  \
+		 rt_io.h \
+		 rt_struct.h \
+			 rt/none/rt.h  \
+			 rt/none/rt_iact.h \
+			 rt/none/rt_io.h \
+			 rt/none/rt_struct.h\
+			 rt/debug/rt.h  \
+			 rt/debug/rt_iact.h \
+			 rt/debug/rt_io.h \
+			 rt/debug/rt_struct.h\
+			 rt/M1closure/rt.h  \
+			 rt/M1closure/rt_iact.h \
+			 rt/M1closure/rt_io.h \
+			 rt/M1closure/rt_struct.h\
 	   stars.h stars_io.h stars_logger.h \
 		 stars/Default/stars.h stars/Default/stars_iact.h stars/Default/stars_io.h \
      stars/Default/stars_debug.h stars/Default/stars_part.h stars/Default/stars_logger.h  \

src/cell.c

@@ -63,6 +63,7 @@
 #include "minmax.h"
 #include "multipole.h"
 #include "pressure_floor.h"
+#include "rt.h"
 #include "scheduler.h"
 #include "space.h"
 #include "space_getsid.h"
@@ -2214,6 +2215,7 @@ void cell_clean_links(struct cell *c, void *data) {
   c->hydro.gradient = NULL;
   c->hydro.force = NULL;
   c->hydro.limiter = NULL;
+  c->hydro.rt_inject = NULL;
   c->grav.grav = NULL;
   c->grav.mm = NULL;
   c->stars.density = NULL;
@@ -3838,6 +3840,99 @@ 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) {
+  const struct engine *e = s->space->e;
+
+  /* Self interaction? */
+  if (cj == NULL) {
+    /* Do anything? */
+    if (ci->hydro.count == 0 || !cell_is_active_hydro(ci, e)) return;
+
+    /* Recurse? */
+    if (cell_can_recurse_in_self_hydro_task(ci)) {
+      /* Loop over all progenies and pairs of progenies */
+      for (int j = 0; j < 8; j++) {
+        if (ci->progeny[j] != NULL) {
+          cell_activate_subcell_rt_tasks(ci->progeny[j], NULL, s);
+          for (int k = j + 1; k < 8; k++)
+            if (ci->progeny[k] != NULL)
+              cell_activate_subcell_rt_tasks(ci->progeny[j], ci->progeny[k], s);
+        }
+      }
+    } else {
+      /* We have reached the bottom of the tree: activate tasks */
+      for (struct link *l = ci->hydro.rt_inject; l != NULL; l = l->next) {
+        struct task *t = l->t;
+        const int ci_active = cell_is_active_hydro(ci, e);
+#ifdef WITH_MPI
+        const int ci_nodeID = ci->nodeID;
+#else
+        const int ci_nodeID = e->nodeID;
+#endif
+        /* Only activate tasks that involve a local active cell. */
+        if (ci_active && ci_nodeID == e->nodeID) {
+          scheduler_activate(s, t);
+        }
+      }
+    }
+  }
+
+  /* Otherwise, pair interaction */
+  else {
+    /* Should we even bother? */
+    if (!cell_is_active_hydro(ci, e) && !cell_is_active_hydro(cj, e)) return;
+    if (ci->hydro.count == 0 || cj->hydro.count == 0) return;
+
+    /* Get the orientation of the pair. */
+    double shift[3];
+    const int sid = space_getsid(s->space, &ci, &cj, shift);
+
+    /* recurse? */
+    if (cell_can_recurse_in_pair_hydro_task(ci) &&
+        cell_can_recurse_in_pair_hydro_task(cj)) {
+      const struct cell_split_pair *csp = &cell_split_pairs[sid];
+      for (int k = 0; k < csp->count; k++) {
+        const int pid = csp->pairs[k].pid;
+        const int pjd = csp->pairs[k].pjd;
+        if (ci->progeny[pid] != NULL && cj->progeny[pjd] != NULL)
+          cell_activate_subcell_rt_tasks(ci->progeny[pid], cj->progeny[pjd], s);
+      }
+    }
+
+    /* Otherwise, activate the RT tasks. */
+    else if (cell_is_active_hydro(ci, e) || cell_is_active_hydro(cj, e)) {
+
+      /* Activate the drifts if the cells are local. */
+      for (struct link *l = ci->hydro.rt_inject; l != NULL; l = l->next) {
+        struct task *t = l->t;
+        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 = e->nodeID;
+        const int cj_nodeID = e->nodeID;
+#endif
+
+        /* Only activate tasks that involve a local active cell. */
+        if ((ci_active && ci_nodeID == e->nodeID) ||
+            (cj_active && cj_nodeID == e->nodeID)) {
+          scheduler_activate(s, t);
+        }
+      }
+    }
+  }
+}
+
 /**
  * @brief Un-skips all the hydro tasks associated with a given cell and checks
  * if the space needs to be rebuilt.
@@ -4819,6 +4914,56 @@ 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)) {
+
+      scheduler_activate(s, t);
+
+      if (t->type == task_type_sub_self) {
+        cell_activate_subcell_rt_tasks(ci, NULL, s);
+      }
+
+      else if (t->type == task_type_sub_pair) {
+        cell_activate_subcell_rt_tasks(ci, cj, s);
+      }
+    }
+  }
+
+  return rebuild;
+}
+
 /**
  * @brief Set the super-cell pointers for all cells in a hierarchy.
  *
@@ -5137,6 +5282,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
         tracers_after_init(p, xp, e->internal_units, e->physical_constants,
                            with_cosmology, e->cosmology, e->hydro_properties,
                            e->cooling_func, e->time);
+        rt_init_part(p);
       }
     }
 
@@ -5458,6 +5604,7 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force) {
       if (spart_is_active(sp, e)) {
         stars_init_spart(sp);
         feedback_init_spart(sp);
+        rt_init_spart(sp);
       }
     }
 

src/cell.h

@@ -407,6 +407,9 @@ struct cell {
     /*! Task for sorting the stars again after a SF event */
     struct task *stars_resort;
 
+    /*! Task for self/pair injection step of radiative transfer */
+    struct link *rt_inject;
+
     /*! Last (integer) time the cell's part were drifted forward in time. */
     integertime_t ti_old_part;
 
@@ -975,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);
@@ -1004,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/common_io.c

@@ -42,6 +42,7 @@
 #include "kernel_hydro.h"
 #include "part.h"
 #include "part_type.h"
+#include "rt_io.h"
 #include "sink_io.h"
 #include "star_formation_io.h"
 #include "stars_io.h"
@@ -2867,6 +2868,7 @@ int get_ptype_fields(const int ptype, struct io_props* list,
         num_fields += fof_write_parts(NULL, NULL, list + num_fields);
       if (with_stf)
         num_fields += velociraptor_write_parts(NULL, NULL, list + num_fields);
+      num_fields += rt_write_particles(NULL, list + num_fields);
       break;
 
     case swift_type_dark_matter:
@@ -2892,6 +2894,7 @@ int get_ptype_fields(const int ptype, struct io_props* list,
       if (with_fof) num_fields += fof_write_sparts(NULL, list + num_fields);
       if (with_stf)
         num_fields += velociraptor_write_sparts(NULL, list + num_fields);
+      num_fields += rt_write_stars(NULL, list + num_fields);
       break;
 
     case swift_type_sink:

src/distributed_io.c

@@ -54,6 +54,7 @@
 #include "output_options.h"
 #include "part.h"
 #include "part_type.h"
+#include "rt_io.h"
 #include "sink_io.h"
 #include "star_formation_io.h"
 #include "stars_io.h"
@@ -268,6 +269,7 @@ void write_output_distributed(struct engine* e,
   const int with_temperature = e->policy & engine_policy_temperature;
   const int with_fof = e->policy & engine_policy_fof;
   const int with_DM_background = e->s->with_DM_background;
+  const int with_rt = e->policy & engine_policy_rt;
 #ifdef HAVE_VELOCIRAPTOR
   const int with_stf = (e->policy & engine_policy_structure_finding) &&
                        (e->s->gpart_group_data != NULL);
@@ -515,6 +517,9 @@ void write_output_distributed(struct engine* e,
               parts, xparts, list + num_fields, with_cosmology);
           num_fields +=
               star_formation_write_particles(parts, xparts, list + num_fields);
+          if (with_rt) {
+            num_fields += rt_write_particles(parts, list + num_fields);
+          }
 
         } else {
 
@@ -557,6 +562,9 @@ void write_output_distributed(struct engine* e,
               parts_written, xparts_written, list + num_fields, with_cosmology);
           num_fields += star_formation_write_particles(
               parts_written, xparts_written, list + num_fields);
+          if (with_rt) {
+            num_fields += rt_write_particles(parts_written, list + num_fields);
+          }
         }
       } break;
 
@@ -691,6 +699,9 @@ void write_output_distributed(struct engine* e,
           if (with_stf) {
             num_fields += velociraptor_write_sparts(sparts, list + num_fields);
           }
+          if (with_rt) {
+            num_fields += rt_write_stars(sparts, list + num_fields);
+          }
         } else {
 
           /* Ok, we need to fish out the particles we want */
@@ -720,6 +731,9 @@ void write_output_distributed(struct engine* e,
             num_fields +=
                 velociraptor_write_sparts(sparts_written, list + num_fields);
           }
+          if (with_rt) {
+            num_fields += rt_write_stars(sparts_written, list + num_fields);
+          }
         }
       } break;
 

src/engine.c

@@ -1631,6 +1631,10 @@ int engine_estimate_nr_tasks(const struct engine *e) {
     n1 += 1;
   }
 #endif
+  if (e->policy & engine_policy_rt) {
+    /* 1 self (inject), (3^3-1)/2 = 26/2 = 13 inject pairs */
+    n1 += 14;
+  }
 
 #ifdef WITH_MPI
 
@@ -2102,7 +2106,9 @@ void engine_skip_force_and_kick(struct engine *e) {
         t->subtype == task_subtype_tend_spart ||
         t->subtype == task_subtype_tend_sink ||
         t->subtype == task_subtype_tend_bpart ||
-        t->subtype == task_subtype_rho || t->subtype == task_subtype_sf_counts)
+        t->subtype == task_subtype_rho ||
+        t->subtype == task_subtype_sf_counts ||
+        t->subtype == task_subtype_rt_inject)
       t->skip = 1;
   }
 

src/engine_maketasks.c

@@ -1777,6 +1777,7 @@ void engine_link_gravity_tasks(struct engine *e) {
           scheduler_addunlock(sched, t, cj_parent->grav.down_in);
         }
       }
+
     }
 
     /* Otherwise M-M interaction? */
@@ -1873,6 +1874,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
   const int with_timestep_sync = (e->policy & engine_policy_timestep_sync);
   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 EXTRA_HYDRO_LOOP
   struct task *t_gradient = NULL;
 #endif
@@ -1885,6 +1887,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
   struct task *t_do_gas_swallow = NULL;
   struct task *t_do_bh_swallow = NULL;
   struct task *t_bh_feedback = NULL;
+  struct task *t_rt_inject = NULL;
 
   for (int ind = 0; ind < num_elements; ind++) {
 
@@ -1955,6 +1958,11 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                               flags, 0, ci, NULL);
       }
 
+      if (with_rt) {
+        t_rt_inject = scheduler_addtask(
+            sched, task_type_self, task_subtype_rt_inject, flags, 0, ci, NULL);
+      }
+
       /* Link the tasks to the cells */
       engine_addlink(e, &ci->hydro.force, t_force);
       if (with_timestep_limiter) {
@@ -1971,6 +1979,9 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         engine_addlink(e, &ci->black_holes.do_bh_swallow, t_do_bh_swallow);
         engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
       }
+      if (with_rt) {
+        engine_addlink(e, &ci->hydro.rt_inject, t_rt_inject);
+      }
 
 #ifdef EXTRA_HYDRO_LOOP
 
@@ -2058,6 +2069,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
       if (with_timestep_sync && with_black_holes && bcount_i > 0) {
         scheduler_addunlock(sched, t_bh_feedback, ci->super->timestep_sync);
       }
+
+      if (with_rt) {
+        scheduler_addunlock(sched, ci->hydro.super->stars.stars_out,
+                            t_rt_inject);
+        scheduler_addunlock(sched, t_rt_inject, ci->super->timestep);
+      }
     }
 
     /* Otherwise, pair interaction? */
@@ -2116,6 +2133,11 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
             sched, task_type_pair, task_subtype_bh_feedback, flags, 0, ci, cj);
       }
 
+      if (with_rt) {
+        t_rt_inject = scheduler_addtask(
+            sched, task_type_pair, task_subtype_rt_inject, flags, 0, ci, cj);
+      }
+
       engine_addlink(e, &ci->hydro.force, t_force);
       engine_addlink(e, &cj->hydro.force, t_force);
       if (with_timestep_limiter) {
@@ -2140,6 +2162,10 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
         engine_addlink(e, &cj->black_holes.feedback, t_bh_feedback);
       }
+      if (with_rt) {
+        engine_addlink(e, &ci->hydro.rt_inject, t_rt_inject);
+        engine_addlink(e, &cj->hydro.rt_inject, t_rt_inject);
+      }
 
 #ifdef EXTRA_HYDRO_LOOP
 
@@ -2259,6 +2285,13 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
             (bcount_i > 0 || bcount_j > 0)) {
           scheduler_addunlock(sched, t_bh_feedback, ci->super->timestep_sync);
         }
+
+        if (with_rt) {
+          scheduler_addunlock(sched, ci->hydro.super->stars.stars_out,
+                              t_rt_inject);
+          scheduler_addunlock(sched, t_rt_inject, ci->super->timestep);
+        }
+
       } else /*(ci->nodeID != nodeID) */ {
         if (with_feedback) {
           scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
@@ -2332,6 +2365,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                                 cj->hydro.super->black_holes.black_holes_out);
           }
 
+          if (with_rt) {
+            scheduler_addunlock(sched, cj->hydro.super->stars.stars_out,
+                                t_rt_inject);
+            scheduler_addunlock(sched, t_rt_inject, cj->super->timestep);
+          }
+
           if (with_timestep_limiter) {
             scheduler_addunlock(sched, cj->hydro.super->hydro.drift, t_limiter);
           }
@@ -2421,6 +2460,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                               task_subtype_bh_feedback, flags, 0, ci, NULL);
       }
 
+      if (with_rt) {
+        t_rt_inject =
+            scheduler_addtask(sched, task_type_sub_self, task_subtype_rt_inject,
+                              flags, 0, ci, NULL);
+      }
+
       /* Add the link between the new loop and the cell */
       engine_addlink(e, &ci->hydro.force, t_force);
       if (with_timestep_limiter) {
@@ -2437,6 +2482,9 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         engine_addlink(e, &ci->black_holes.do_bh_swallow, t_do_bh_swallow);
         engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
       }
+      if (with_rt) {
+        engine_addlink(e, &ci->hydro.rt_inject, t_rt_inject);
+      }
 
 #ifdef EXTRA_HYDRO_LOOP
 
@@ -2530,6 +2578,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
       if (with_timestep_sync && with_black_holes && bcount_i > 0) {
         scheduler_addunlock(sched, t_bh_feedback, ci->super->timestep_sync);
       }
+
+      if (with_rt) {
+        scheduler_addunlock(sched, ci->hydro.super->stars.stars_out,
+                            t_rt_inject);
+        scheduler_addunlock(sched, t_rt_inject, ci->super->timestep);
+      }
     }
 
     /* Otherwise, sub-pair interaction? */
@@ -2592,6 +2646,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                               task_subtype_bh_feedback, flags, 0, ci, cj);
       }
 
+      if (with_rt) {
+        t_rt_inject =
+            scheduler_addtask(sched, task_type_sub_pair, task_subtype_rt_inject,
+                              flags, 0, ci, cj);
+      }
+
       engine_addlink(e, &ci->hydro.force, t_force);
       engine_addlink(e, &cj->hydro.force, t_force);
       if (with_timestep_limiter) {
@@ -2616,6 +2676,10 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
         engine_addlink(e, &cj->black_holes.feedback, t_bh_feedback);
       }
+      if (with_rt) {
+        engine_addlink(e, &ci->hydro.rt_inject, t_rt_inject);
+        engine_addlink(e, &cj->hydro.rt_inject, t_rt_inject);
+      }
 
 #ifdef EXTRA_HYDRO_LOOP
 
@@ -2734,6 +2798,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
             (bcount_i > 0 || bcount_j > 0)) {
           scheduler_addunlock(sched, t_bh_feedback, ci->super->timestep_sync);
         }
+
+        if (with_rt) {
+          scheduler_addunlock(sched, ci->hydro.super->stars.stars_out,
+                              t_rt_inject);
+          scheduler_addunlock(sched, t_rt_inject, ci->super->timestep);
+        }
       } else /* ci->nodeID != nodeID */ {
 
         if (with_feedback) {
@@ -2807,6 +2877,11 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
             scheduler_addunlock(sched, t_bh_feedback,
                                 cj->hydro.super->black_holes.black_holes_out);
           }
+          if (with_rt) {
+            scheduler_addunlock(sched, cj->hydro.super->stars.stars_out,
+                                t_rt_inject);
+            scheduler_addunlock(sched, t_rt_inject, cj->super->timestep);
+          }
 
           if (with_timestep_limiter) {
             scheduler_addunlock(sched, cj->hydro.super->hydro.drift, t_limiter);

src/engine_marktasks.c

@@ -261,6 +261,13 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         }
       }
 
+      /* Activate RT injection */
+      else if (t_subtype == task_subtype_rt_inject) {
+        if (ci_active_hydro) {
+          scheduler_activate(s, t);
+        }
+      }
+
 #ifdef SWIFT_DEBUG_CHECKS
       else {
         error("Invalid task type / sub-type encountered");
@@ -471,6 +478,24 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
 #endif
       }
 
+      /* RT injection tasks */
+      else if (t_subtype == task_subtype_rt_inject) {
+        /* We only want to activate the task if the cell is active and is
+           going to update some gas on the *local* node */
+        if ((ci_nodeID == nodeID && cj_nodeID == nodeID) &&
+            (ci_active_hydro || cj_active_hydro)) {
+          scheduler_activate(s, t);
+
+        } else if ((ci_nodeID == nodeID && cj_nodeID != nodeID) &&
+                   (cj_active_hydro)) {
+          scheduler_activate(s, t);
+
+        } else if ((ci_nodeID != nodeID && cj_nodeID == nodeID) &&
+                   (ci_active_hydro)) {
+          scheduler_activate(s, t);
+        }
+      }
+
       /* Only interested in density tasks as of here. */
       if (t_subtype == task_subtype_density) {
 
@@ -858,6 +883,11 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
             scheduler_activate_send(s, ci->mpi.send, task_subtype_tend_gpart,
                                     cj_nodeID);
         }
+#endif
+      } /* Only interested in RT tasks as of here. */
+      else if (t_subtype == task_subtype_rt_inject) {
+#ifdef WITH_MPI
+        error("RT doesn't work with MPI yet.");
 #endif
       }
     }

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,38 @@ 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 +349,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 +383,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 +437,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/hydro/Gadget2/hydro_part.h

@@ -37,6 +37,7 @@
 #include "feedback_struct.h"
 #include "logger.h"
 #include "pressure_floor_struct.h"
+#include "rt_struct.h"
 #include "star_formation_struct.h"
 #include "timestep_limiter_struct.h"
 #include "tracers_struct.h"
@@ -167,6 +168,9 @@ struct part {
   /*! Additional data used by the pressure floor */
   struct pressure_floor_part_data pressure_floor_data;
 
+  /* Additional Radiative Transfer Data */
+  struct rt_part_data rt_data;
+
   /*! Time-step length */
   timebin_t time_bin;
 

src/line_of_sight.c

@@ -36,6 +36,7 @@
 #include "kernel_hydro.h"
 #include "line_of_sight.h"
 #include "periodic.h"
+#include "rt_io.h"
 #include "star_formation_io.h"
 #include "tracers_io.h"
 #include "velociraptor_io.h"
@@ -432,6 +433,7 @@ void write_los_hdf5_datasets(hid_t grp, const int j, const size_t N,
 #else
   const int with_stf = 0;
 #endif
+  const int with_rt = e->policy & engine_policy_rt;