Fixes and improvements to h_max_active

src/cell_drift.c

@@ -588,6 +588,7 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force,
   const int periodic = e->s->periodic;
   const double dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
   const int with_cosmology = (e->policy & engine_policy_cosmology);
+  const int with_rt = (e->policy & engine_policy_rt);
   const float stars_h_max = e->hydro_properties->h_max;
   const float stars_h_min = e->hydro_properties->h_min;
   const integertime_t ti_old_spart = c->stars.ti_old_part;
@@ -753,7 +754,8 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force,
         rt_init_spart(sp);
 
         /* Update the maximal active smoothing length in the cell */
-        cell_h_max_active = max(cell_h_max_active, sp->h);
+        if (feedback_is_active(sp, e) || with_rt)
+          cell_h_max_active = max(cell_h_max_active, sp->h);
       }
     }
 

src/runner_ghost.c

@@ -99,7 +99,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
 
   /* Running value of the maximal smoothing length */
   float h_max = c->stars.h_max;
-  float h_max_active = c->stars.h_max_active;
+  float h_max_active = 0.f;
 
   TIMER_TIC;
 
@@ -297,6 +297,13 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
                                                rt_props, phys_const, us);
             }
 
+            if (feedback_is_active(sp, e) || with_rt) {
+
+              /* Check if h_max has increased */
+              h_max = max(h_max, sp->h);
+              h_max_active = max(h_max_active, sp->h);
+            }
+
             /* Ok, we are done with this particle */
             continue;
           }
@@ -561,7 +568,9 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
 
     if (h > c->stars.h_max)
       error("Particle has h larger than h_max (id=%lld)", sp->id);
-    if (spart_is_active(sp, e) && h > c->stars.h_max_active)
+
+    if (spart_is_active(sp, e) && (feedback_is_active(sp, e) || with_rt) &&
+        (h > c->stars.h_max_active))
       error("Active particle has h larger than h_max_active (id=%lld)", sp->id);
   }
 #endif
@@ -602,7 +611,7 @@ void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
 
   /* Running value of the maximal smoothing length */
   float h_max = c->black_holes.h_max;
-  float h_max_active = c->black_holes.h_max_active;
+  float h_max_active = 0.f;
 
   TIMER_TIC;
 
@@ -721,6 +730,10 @@ void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
 
             black_holes_reset_feedback(bp);
 
+            /* Check if h_max has increased */
+            h_max = max(h_max, bp->h);
+            h_max_active = max(h_max_active, bp->h);
+
             /* Ok, we are done with this particle */
             continue;
           }
@@ -1125,7 +1138,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 
   /* Running value of the maximal smoothing length */
   float h_max = c->hydro.h_max;
-  float h_max_active = c->hydro.h_max_active;
+  float h_max_active = 0.f;
 
   TIMER_TIC;
 
@@ -1327,6 +1340,10 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
               rt_reset_part(p, cosmo);
             }
 
+            /* Check if h_max has increased */
+            h_max = max(h_max, p->h);
+            h_max_active = max(h_max_active, p->h);
+
             /* Ok, we are done with this particle */
             continue;
           }
@@ -1749,7 +1766,7 @@ void runner_do_sinks_density_ghost(struct runner *r, struct cell *c,
 
   /* Running value of the maximal smoothing length */
   float h_max = c->sinks.h_max;
-  float h_max_active = c->sinks.h_max_active;
+  float h_max_active = 0.f;
 
   TIMER_TIC;
 
@@ -1887,6 +1904,10 @@ void runner_do_sinks_density_ghost(struct runner *r, struct cell *c,
             if (((sp->h >= sinks_h_max) && (f < 0.f)) ||
                 ((sp->h <= sinks_h_min) && (f > 0.f))) {
 
+              /* Check if h_max has increased */
+              h_max = max(h_max, sp->h);
+              h_max_active = max(h_max_active, sp->h);
+
               /* Ok, we are done with this particle */
               continue;
             }

src/runner_recv.c

@@ -31,7 +31,9 @@
 #include "runner.h"
 
 /* Local headers. */
+#include "active.h"
 #include "engine.h"
+#include "feedback.h"
 #include "timers.h"
 
 /**
@@ -225,10 +227,11 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
 
 #ifdef WITH_MPI
 
+  const struct engine *e = r->e;
+  const int with_rt = (e->policy & engine_policy_rt);
+  const integertime_t ti_current = e->ti_current;
   struct spart *restrict sparts = c->stars.parts;
   const size_t nr_sparts = c->stars.count;
-  const integertime_t ti_current = r->e->ti_current;
-  const timebin_t max_active_bin = r->e->max_active_bin;
 
   TIMER_TIC;
 
@@ -258,7 +261,8 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
       time_bin_max = max(time_bin_max, sparts[k].time_bin);
       h_max = max(h_max, sparts[k].h);
       sparts[k].gpart = NULL;
-      if (sparts[k].time_bin <= max_active_bin)
+      if (spart_is_active(&sparts[k], e) &&
+          (feedback_is_active(&sparts[k], e) || with_rt))
         h_max_active = max(h_max_active, sparts[k].h);
     }