diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h
index 9f241bb7abda99fb6f6e6c56f891b6c819664ab6..267daa91190a70eefe18051bf60ef49c1c24b8ca 100644
--- a/src/black_holes/EAGLE/black_holes_iact.h
+++ b/src/black_holes/EAGLE/black_holes_iact.h
@@ -19,7 +19,9 @@
#ifndef SWIFT_EAGLE_BH_IACT_H
#define SWIFT_EAGLE_BH_IACT_H
+/* Local includes */
#include "hydro.h"
+#include "random.h"
/**
* @brief Density interaction between two particles (non-symmetric).
@@ -35,8 +37,9 @@
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_density(
const float r2, const float *dx, const float hi, const float hj,
- struct bpart *restrict bi, const struct part *restrict pj, const float a,
- const float H) {
+ struct bpart *restrict bi, const struct part *restrict pj,
+ const struct xpart *restrict xpj, const struct cosmology *cosmo,
+ const integertime_t ti_current) {
float wi, wi_dx;
@@ -104,8 +107,34 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_density(
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_feedback(const float r2, const float *dx, const float hi,
const float hj, struct bpart *restrict bi,
- struct part *restrict pj, const float a,
- const float H) {
+ struct part *restrict pj,
+ struct xpart *restrict xpj,
+ const struct cosmology *cosmo,
+ const integertime_t ti_current) {
+
+ /* Get the heating probability */
+ const float prob = bi->to_distribute.AGN_heating_probability;
+
+ /* Are we doing some feedback? */
+ if (prob > 0.f) {
+
+ /* Draw a random number (Note mixing both IDs) */
+ const float rand = random_unit_interval(bi->id + pj->id, ti_current,
+ random_number_BH_feedback);
+
+ /* Are we lucky? */
+ if (rand < prob) {
+
+ /* Compute new energy of this particle */
+ const double u_init = hydro_get_physical_internal_energy(pj, xpj, cosmo);
+ const float delta_u = bi->to_distribute.AGN_delta_u;
+ const double u_new = u_init + delta_u;
+
+ hydro_set_physical_internal_energy(pj, xpj, cosmo, u_new);
+ hydro_set_drifted_physical_internal_energy(pj, cosmo, u_new);
+ }
+ }
+
#ifdef DEBUG_INTERACTIONS_BH
/* Update ngb counters */
if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_BH)
diff --git a/src/runner_doiact_black_holes.h b/src/runner_doiact_black_holes.h
index 47604f2a1ea21d7e89f72a0ef4d441c45254b96c..323ecb254d7d1af9750ee3d7a8d362267bb6e0d6 100644
--- a/src/runner_doiact_black_holes.h
+++ b/src/runner_doiact_black_holes.h
@@ -98,23 +98,18 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
const struct engine *e = r->e;
- // const int with_cosmology = e->policy & engine_policy_cosmology;
- // const integertime_t ti_current = e->ti_current;
+ const integertime_t ti_current = e->ti_current;
const struct cosmology *cosmo = e->cosmology;
/* Anything to do here? */
if (c->hydro.count == 0 || c->black_holes.count == 0) return;
if (!cell_is_active_black_holes(c, e)) return;
- /* Cosmological terms */
- const float a = cosmo->a;
- const float H = cosmo->H;
-
const int bcount = c->black_holes.count;
const int count = c->hydro.count;
struct bpart *restrict bparts = c->black_holes.parts;
struct part *restrict parts = c->hydro.parts;
- // struct xpart *restrict xparts = c->hydro.xparts;
+ struct xpart *restrict xparts = c->hydro.xparts;
/* Loop over the bparts in ci. */
for (int sid = 0; sid < bcount; sid++) {
@@ -136,7 +131,7 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
/* Get a pointer to the jth particle. */
struct part *restrict pj = &parts[pjd];
- // struct xpart *restrict xpj = &xparts[pjd];
+ struct xpart *restrict xpj = &xparts[pjd];
const float hj = pj->h;
/* Early abort? */
@@ -156,7 +151,7 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
#endif
if (r2 < hig2) {
- IACT_BH(r2, dx, hi, hj, si, pj, a, H);
+ IACT_BH(r2, dx, hi, hj, si, pj, xpj, cosmo, ti_current);
}
} /* loop over the parts in ci. */
} /* loop over the bparts in ci. */
@@ -183,23 +178,18 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
#endif
const struct engine *e = r->e;
- // const int with_cosmology = e->policy & engine_policy_cosmology;
- // const integertime_t ti_current = e->ti_current;
+ const integertime_t ti_current = e->ti_current;
const struct cosmology *cosmo = e->cosmology;
/* Anything to do here? */
if (cj->hydro.count == 0 || ci->black_holes.count == 0) return;
if (!cell_is_active_black_holes(ci, e)) return;
- /* Cosmological terms */
- const float a = cosmo->a;
- const float H = cosmo->H;
-
const int bcount_i = ci->black_holes.count;
const int count_j = cj->hydro.count;
struct bpart *restrict bparts_i = ci->black_holes.parts;
struct part *restrict parts_j = cj->hydro.parts;
- // struct xpart *restrict xparts_j = cj->hydro.xparts;
+ struct xpart *restrict xparts_j = cj->hydro.xparts;
/* Get the relative distance between the pairs, wrapping. */
double shift[3] = {0.0, 0.0, 0.0};
@@ -230,7 +220,7 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
/* Get a pointer to the jth particle. */
struct part *restrict pj = &parts_j[pjd];
- // struct xpart *restrict xpj = &xparts_j[pjd];
+ struct xpart *restrict xpj = &xparts_j[pjd];
const float hj = pj->h;
/* Skip inhibited particles. */
@@ -250,7 +240,7 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
#endif
if (r2 < hig2) {
- IACT_BH(r2, dx, hi, hj, si, pj, a, H);
+ IACT_BH(r2, dx, hi, hj, si, pj, xpj, cosmo, ti_current);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -301,16 +291,12 @@ void DOPAIR1_SUBSET_BH_NAIVE(struct runner *r, struct cell *restrict ci,
#endif
const struct engine *e = r->e;
- // const integertime_t ti_current = e->ti_current;
+ const integertime_t ti_current = e->ti_current;
const struct cosmology *cosmo = e->cosmology;
- /* Cosmological terms */
- const float a = cosmo->a;
- const float H = cosmo->H;
-
const int count_j = cj->hydro.count;
struct part *restrict parts_j = cj->hydro.parts;
- // struct xpart *restrict xparts_j = cj->hydro.xparts;
+ struct xpart *restrict xparts_j = cj->hydro.xparts;
/* Early abort? */
if (count_j == 0) return;
@@ -337,7 +323,7 @@ void DOPAIR1_SUBSET_BH_NAIVE(struct runner *r, struct cell *restrict ci,
/* Get a pointer to the jth particle. */
struct part *restrict pj = &parts_j[pjd];
- // struct xpart *restrict xpj = &xparts_j[pjd];
+ struct xpart *restrict xpj = &xparts_j[pjd];
/* Skip inhibited particles */
if (part_is_inhibited(pj, e)) continue;
@@ -359,7 +345,7 @@ void DOPAIR1_SUBSET_BH_NAIVE(struct runner *r, struct cell *restrict ci,
#endif
/* Hit or miss? */
if (r2 < hig2) {
- IACT_BH(r2, dx, hi, hj, bpi, pj, a, H);
+ IACT_BH(r2, dx, hi, hj, bpi, pj, xpj, cosmo, ti_current);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -384,16 +370,12 @@ void DOSELF1_SUBSET_BH(struct runner *r, struct cell *restrict ci,
#endif
const struct engine *e = r->e;
- // const integertime_t ti_current = e->ti_current;
+ const integertime_t ti_current = e->ti_current;
const struct cosmology *cosmo = e->cosmology;
- /* Cosmological terms */
- const float a = cosmo->a;
- const float H = cosmo->H;
-
const int count_i = ci->hydro.count;
struct part *restrict parts_j = ci->hydro.parts;
- // struct xpart *restrict xparts_j = ci->hydro.xparts;
+ struct xpart *restrict xparts_j = ci->hydro.xparts;
/* Early abort? */
if (count_i == 0) return;
@@ -419,7 +401,7 @@ void DOSELF1_SUBSET_BH(struct runner *r, struct cell *restrict ci,
/* Get a pointer to the jth particle. */
struct part *restrict pj = &parts_j[pjd];
- // struct xpart *restrict xpj = &xparts_j[pjd];
+ struct xpart *restrict xpj = &xparts_j[pjd];
/* Early abort? */
if (part_is_inhibited(pj, e)) continue;
@@ -439,7 +421,7 @@ void DOSELF1_SUBSET_BH(struct runner *r, struct cell *restrict ci,
/* Hit or miss? */
if (r2 < hig2) {
- IACT_BH(r2, dx, hi, pj->h, bpi, pj, a, H);
+ IACT_BH(r2, dx, hi, pj->h, bpi, pj, xpj, cosmo, ti_current);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */