From 8d045eef73af072f2f15e4b08e1f6a2bdcb297f8 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Mon, 28 Apr 2025 07:38:09 +0000
Subject: [PATCH] Fixes and improvements to h_max_active

---
 src/cell_drift.c   |  4 +++-
 src/runner_ghost.c | 31 ++++++++++++++++++++++++++-----
 src/runner_recv.c  | 10 +++++++---
 3 files changed, 36 insertions(+), 9 deletions(-)

diff --git a/src/cell_drift.c b/src/cell_drift.c
index 6454dc0392..ea427d140a 100644
--- a/src/cell_drift.c
+++ b/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);
       }
     }
 
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index 77d584a37c..830298e2ca 100644
--- a/src/runner_ghost.c
+++ b/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;
             }
diff --git a/src/runner_recv.c b/src/runner_recv.c
index c36f6d2349..5b019f2d0d 100644
--- a/src/runner_recv.c
+++ b/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);
     }
 
-- 
GitLab