From 02ed3c26dfa96e75488d079d021911039106faae Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Thu, 15 Oct 2020 16:48:22 +0200
Subject: [PATCH] Add a random AGN feedback model to reproduce the old
 behaviour of EAGLE where particles were selected randomly in the kernel

---
 src/black_holes/EAGLE/black_holes_iact.h      | 20 +++++++++++++++++++
 .../EAGLE/black_holes_properties.h            | 11 ++++++----
 2 files changed, 27 insertions(+), 4 deletions(-)

diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h
index 21a4da8474..ae80ac82bc 100644
--- a/src/black_holes/EAGLE/black_holes_iact.h
+++ b/src/black_holes/EAGLE/black_holes_iact.h
@@ -197,6 +197,26 @@ runner_iact_nonsym_bh_gas_density(
       ray_minimise_distance(r, bi->rays, arr_size, gas_id, pj->mass);
       break;
     }
+    case AGN_random_ngb_model: {
+      /* Compute the size of the array that we want to sort. If the current
+       * function is called for the first time (at this time-step for this BH),
+       * then bi->num_ngbs = 1 and there is nothing to sort. Note that the
+       * maximum size of the sorted array cannot be larger then the maximum
+       * number of rays. */
+      const int arr_size = min(bi->num_ngbs, eagle_blackhole_number_of_rays);
+
+      /* To mimic a random draw among all the particles in the kernel, we
+       * draw random distances in [0,1) and then pick the particle(s) with
+       * the smallest of these 'fake' distances */
+      const float dist = random_unit_interval_two_IDs(
+          bi->id, pj->id, ti_current, random_number_BH_feedback);
+
+      /* Minimise separation between the gas particles and the BH. The rays
+       * structs with smaller ids in the ray array will refer to the particles
+       * with smaller distances to the BH. */
+      ray_minimise_distance(dist, bi->rays, arr_size, gas_id, pj->mass);
+      break;
+    }
   }
 }
 
diff --git a/src/black_holes/EAGLE/black_holes_properties.h b/src/black_holes/EAGLE/black_holes_properties.h
index ba6a07b0bb..b00ec57663 100644
--- a/src/black_holes/EAGLE/black_holes_properties.h
+++ b/src/black_holes/EAGLE/black_holes_properties.h
@@ -25,6 +25,7 @@
 #include <string.h>
 
 enum AGN_feedback_models {
+  AGN_random_ngb_model,      /*< Random neighbour model for AGN feedback */
   AGN_isotropic_model,       /*< Isotropic model of AGN feedback */
   AGN_minimum_distance_model /*< Minimum-distance model of AGN feedback */
 };
@@ -362,15 +363,17 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
   /* Feedback parameters ---------------------------------- */
 
   char temp[40];
-  parser_get_param_string(params, "EAGLEAGN:AGN_feedback_model", temp);
-  if (strcmp(temp, "Isotropic") == 0)
+  parser_get_param_string(params, "COLIBREAGN:AGN_feedback_model", temp);
+  if (strcmp(temp, "Random") == 0)
+    bp->feedback_model = AGN_random_ngb_model;
+  else if (strcmp(temp, "Isotropic") == 0)
     bp->feedback_model = AGN_isotropic_model;
   else if (strcmp(temp, "MinimumDistance") == 0)
     bp->feedback_model = AGN_minimum_distance_model;
   else
     error(
-        "The AGN feedback model must be either MinimumDistance or Isotropic, "
-        "not %s",
+        "The AGN feedback model must be either 'Random', 'MinimumDistance' or "
+        "'Isotropic', not %s",
         temp);
 
   bp->AGN_deterministic =
-- 
GitLab