diff --git a/src/black_holes/Default/black_holes.h b/src/black_holes/Default/black_holes.h
index a67d315f0c5deebfddb667aba8c8c08df5089cda..cab7aa70c744881f0bd54b759005a5a824800c37 100644
--- a/src/black_holes/Default/black_holes.h
+++ b/src/black_holes/Default/black_holes.h
@@ -224,7 +224,7 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
const struct phys_const* constants, const struct cosmology* cosmo,
const struct cooling_function_data* cooling,
const struct entropy_floor_properties* floor_props, const double time,
- const int with_cosmology, const double dt) {}
+ const int with_cosmology, const double dt, const integertime_t ti_begin) {}
/**
* @brief Finish the calculation of the new BH position.
diff --git a/src/black_holes/EAGLE/black_holes.h b/src/black_holes/EAGLE/black_holes.h
index 749c9731cdb25ed64ec7a4457dc2966d3a2c99d3..166fda012062ca35e5605635eb3087e26e7852ff 100644
--- a/src/black_holes/EAGLE/black_holes.h
+++ b/src/black_holes/EAGLE/black_holes.h
@@ -29,6 +29,7 @@
#include "kernel_hydro.h"
#include "minmax.h"
#include "physical_constants.h"
+#include "random.h"
/* Standard includes */
#include <float.h>
@@ -72,7 +73,7 @@ __attribute__((always_inline)) INLINE static float black_holes_compute_timestep(
* rate. The time is multiplied by the number of Ngbs to heat because
* if more particles are heated at once then the time between different
* AGN feedback events increases proportionally. */
- const double dt_heat = E_heat * bp->num_ngbs_to_heat / Energy_rate;
+ const double dt_heat = E_heat * props->num_ngbs_to_heat / Energy_rate;
/* The new timestep of the BH cannot be smaller than the miminum allowed
* time-step */
@@ -174,7 +175,6 @@ __attribute__((always_inline)) INLINE static void black_holes_init_bpart(
bp->reposition.potential = FLT_MAX;
bp->accretion_rate = 0.f; /* Optionally accumulated ngb-by-ngb */
bp->f_visc = FLT_MAX;
- bp->accretion_boost_factor = -FLT_MAX;
bp->mass_at_start_of_step = bp->mass; /* bp->mass may grow in nibbling mode */
}
@@ -600,7 +600,7 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
const struct phys_const* constants, const struct cosmology* cosmo,
const struct cooling_function_data* cooling,
const struct entropy_floor_properties* floor_props, const double time,
- const int with_cosmology, const double dt) {
+ const int with_cosmology, const double dt, const integertime_t ti_begin) {
/* Record that the black hole has another active time step */
bp->number_of_time_steps++;
@@ -618,6 +618,7 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
const double f_Edd_recording = props->f_Edd_recording;
const double epsilon_r = props->epsilon_r;
const double epsilon_f = props->epsilon_f;
+ const double num_ngbs_to_heat = props->num_ngbs_to_heat;
const int with_angmom_limiter = props->with_angmom_limiter;
/* (Subgrid) mass of the BH (internal units) */
@@ -841,56 +842,157 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
bp->accreted_angular_momentum[2] +=
bp->circular_velocity_gas[2] * mass_rate * dt / bp->h;
+ /* Below we compute energy required to have a feedback event(s)
+ * Note that we have subtracted the particles we swallowed from the ngb_mass
+ * and num_ngbs accumulators. */
+
/* Now find the temperature increase for a possible feedback event */
const double delta_T = black_hole_feedback_delta_T(bp, props, cosmo);
bp->AGN_delta_T = delta_T;
- const double delta_u = delta_T * props->temp_to_u_factor;
+ double delta_u = delta_T * props->temp_to_u_factor;
const double delta_u_ref =
props->AGN_use_nheat_with_fixed_dT
? props->AGN_delta_T_desired * props->temp_to_u_factor
: delta_u;
- /* Energy required to have a feedback event.
- * Note that we have subtracted particles we may have swallowed from the
- * ngb_mass and num_ngbs accumulators already. */
- const double num_ngbs_to_heat =
- black_hole_energy_reservoir_threshold(bp, props);
+ /* Energy required to have a feedback event
+ * Note that we have subtracted the particles we swallowed from the ngb_mass
+ * and num_ngbs accumulators. */
const double mean_ngb_mass = bp->ngb_mass / ((double)bp->num_ngbs);
const double E_feedback_event =
num_ngbs_to_heat * delta_u_ref * mean_ngb_mass;
+ /* Compute and store BH accretion-limited time-step */
+ if (luminosity > 0.) {
+ const float dt_acc = delta_u * mean_ngb_mass * props->num_ngbs_to_heat /
+ (luminosity * props->epsilon_f);
+ bp->dt_heat = max(dt_acc, props->time_step_min);
+ } else {
+ bp->dt_heat = FLT_MAX;
+ }
+
/* Are we doing some feedback? */
if (bp->energy_reservoir > E_feedback_event) {
- /* Default probability of heating */
- double target_prob = bp->energy_reservoir / (delta_u * bp->ngb_mass);
+ int number_of_energy_injections;
- /* Calculate the change in internal energy of the gas particles that get
- * heated. Adjust the prbability if needed. */
- double gas_delta_u;
- double prob;
- if (target_prob <= 1.) {
+ /* How are we doing feedback? */
+ if (props->AGN_deterministic) {
- /* Normal case */
- prob = target_prob;
- gas_delta_u = delta_u;
+ number_of_energy_injections =
+ (int)(bp->energy_reservoir / (delta_u * mean_ngb_mass));
} else {
- /* Special case: we need to adjust the energy irrespective of the
- * desired deltaT to ensure we inject all the available energy. */
+ /* Probability of heating. */
+ const double prob = bp->energy_reservoir / (delta_u * bp->ngb_mass);
+
+ /* Compute the number of energy injections based on probability */
+ if (prob < 1.) {
+
+ /* Initialise counter of energy injections */
+ number_of_energy_injections = 0;
- prob = 1.;
- gas_delta_u = bp->energy_reservoir / bp->ngb_mass;
+ /* How many AGN energy injections will we get? */
+ for (int i = 0; i < bp->num_ngbs; i++) {
+ const double rand = random_unit_interval_part_ID_and_ray_idx(
+ bp->id, i, ti_begin, random_number_BH_feedback);
+ /* Increase the counter if we are lucky */
+ if (rand < prob) number_of_energy_injections++;
+ }
+
+ } else {
+
+ /* We want to use up all energy avaliable in the reservoir. Therefore,
+ * number_of_energy_injections is > or = props->num_ngbs_to_heat */
+ number_of_energy_injections =
+ (int)(bp->energy_reservoir / (delta_u * mean_ngb_mass));
+ }
}
- /* Store all of this in the black hole for delivery onto the gas. */
- bp->to_distribute.AGN_heating_probability = prob;
- bp->to_distribute.AGN_delta_u = gas_delta_u;
+ /* Maximum number of energy injections allowed */
+ const int N_energy_injections_allowed =
+ min(eagle_blackhole_number_of_rays, bp->num_ngbs);
- /* Decrement the energy in the reservoir by the mean expected energy */
- const double energy_used = bp->energy_reservoir / max(prob, 1.);
- bp->energy_reservoir -= energy_used;
+ /* If there are more energy-injection events than min(the number of Ngbs in
+ * the kernel, maximum number of rays) then lower the number of events &
+ * proportionally increase the energy per event */
+ if (number_of_energy_injections > N_energy_injections_allowed) {
+
+ /* Increase the thermal energy per event */
+ const double alpha_thermal = (double)number_of_energy_injections /
+ (double)N_energy_injections_allowed;
+
+ delta_u *= alpha_thermal;
+
+ /* Lower the maximum number of events to the max allowed value */
+ number_of_energy_injections = N_energy_injections_allowed;
+ }
+
+ /* Compute how much energy will be deposited onto the gas */
+ /* Note that it will in general be different from E_feedback_event if
+ * gas particles are of different mass. */
+ double Energy_deposited = 0.0;
+
+ /* Count the number of unsuccessful energy injections (e.g., if the particle
+ * that the BH wants to heat has been swallowed and thus no longer exists)
+ */
+ int N_unsuccessful_energy_injections = 0;
+
+ for (int i = 0; i < number_of_energy_injections; i++) {
+
+ /* If the gas particle that the BH wants to heat has just been swallowed
+ * by the same BH, increment the counter of unsuccessful injections. If
+ * the particle has not been swallowed by the BH, increase the energy that
+ * will later be subtracted from the BH's energy reservoir. */
+ if (bp->rays[i].id_min_length != -1)
+ Energy_deposited += delta_u * bp->rays[i].mass;
+ else
+ N_unsuccessful_energy_injections++;
+ }
+
+ /* Store all of this in the black hole for delivery onto the gas. */
+ bp->to_distribute.AGN_delta_u = delta_u;
+ bp->to_distribute.AGN_number_of_energy_injections =
+ number_of_energy_injections;
+
+ /* Subtract the deposited energy from the BH energy reservoir. Note
+ * that in the stochastic case, the resulting value might be negative.
+ * This happens when (due to the probabilistic nature of the model) the
+ * BH injects more energy than it actually has in the reservoir. */
+ bp->energy_reservoir -= Energy_deposited;
+
+ /* Total number successful energy injections at this time-step. In each
+ * energy injection, a certain gas particle from the BH's kernel gets
+ * heated. (successful = the particle(s) that is going to get heated by
+ * this BH has not been swallowed by the same BH). */
+ const int N_successful_energy_injections =
+ number_of_energy_injections - N_unsuccessful_energy_injections;
+
+ /* Increase the number of energy injections the black hole has heated so
+ * far. Note that in the isotropic model, a gas particle may receive AGN
+ * energy several times at the same time-step. In this case, the number of
+ * particles heated at this time-step for this BH will be smaller than the
+ * total number of energy injections for this BH. */
+ bp->AGN_number_of_energy_injections += N_successful_energy_injections;
+
+ /* Increase the number of AGN events the black hole has had so far.
+ * If the BH does feedback, the number of AGN events is incremented by one.
+ */
+ bp->AGN_number_of_AGN_events += N_successful_energy_injections > 0;
+
+ /* Update the total (cumulative) energy used for gas heating in AGN feedback
+ * by this BH */
+ bp->AGN_cumulative_energy += Energy_deposited;
+
+ /* Store the time/scale factor when the BH last did AGN feedback */
+ if (N_successful_energy_injections) {
+ if (with_cosmology) {
+ bp->last_AGN_event_scale_factor = cosmo->a;
+ } else {
+ bp->last_AGN_event_time = time;
+ }
+ }
} else {
@@ -1019,8 +1121,8 @@ __attribute__((always_inline)) INLINE static void black_holes_end_reposition(
__attribute__((always_inline)) INLINE static void black_holes_reset_feedback(
struct bpart* restrict bp) {
- bp->to_distribute.AGN_heating_probability = 0.f;
bp->to_distribute.AGN_delta_u = 0.f;
+ bp->to_distribute.AGN_number_of_energy_injections = 0;
#ifdef DEBUG_INTERACTIONS_BLACK_HOLES
for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h
index dd96d2c9a565a24c2b7ded67a13524c1cbfbd85c..21a4da8474c6a54d8136d236ab2e2669bcd141a0 100644
--- a/src/black_holes/EAGLE/black_holes_iact.h
+++ b/src/black_holes/EAGLE/black_holes_iact.h
@@ -24,6 +24,7 @@
#include "gravity.h"
#include "hydro.h"
#include "random.h"
+#include "rays.h"
#include "space.h"
#include "timestep_sync_part.h"
#include "tracers.h"
@@ -149,6 +150,54 @@ runner_iact_nonsym_bh_gas_density(
/* Update ngb counters */
++si->num_ngb_density;
#endif
+
+ /* Gas particle id */
+ const long long gas_id = pj->id;
+
+ /* Choose AGN feedback model */
+ switch (bh_props->feedback_model) {
+ case AGN_isotropic_model: {
+ /* Compute arc lengths in AGN isotropic feedback and collect
+ * relevant data for later use in the feedback_apply loop */
+
+ /* Loop over rays */
+ for (int i = 0; i < eagle_blackhole_number_of_rays; i++) {
+
+ /* We generate two random numbers that we use
+ to randomly select the direction of the ith ray */
+
+ /* Random number in [0, 1[ */
+ const double rand_theta = random_unit_interval_part_ID_and_ray_idx(
+ bi->id, i, ti_current,
+ random_number_isotropic_AGN_feedback_ray_theta);
+
+ /* Random number in [0, 1[ */
+ const double rand_phi = random_unit_interval_part_ID_and_ray_idx(
+ bi->id, i, ti_current,
+ random_number_isotropic_AGN_feedback_ray_phi);
+
+ /* Compute arc length */
+ ray_minimise_arclength(dx, r, bi->rays + i, /*switch=*/-1, gas_id,
+ rand_theta, rand_phi, pj->mass, /*ray_ext=*/NULL,
+ /*v=*/NULL);
+ }
+ break;
+ }
+ case AGN_minimum_distance_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);
+
+ /* 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(r, bi->rays, arr_size, gas_id, pj->mass);
+ break;
+ }
+ }
}
/**
@@ -551,22 +600,34 @@ runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
const integertime_t ti_current,
const double time) {
- /* Get the heating probability */
- const float prob = bi->to_distribute.AGN_heating_probability;
+ /* Number of energy injections per BH per time-step */
+ const int num_energy_injections_per_BH =
+ bi->to_distribute.AGN_number_of_energy_injections;
/* Are we doing some feedback? */
- if (prob > 0.f) {
+ if (num_energy_injections_per_BH > 0) {
- /* Draw a random number (Note mixing both IDs) */
- const float rand = random_unit_interval(bi->id + pj->id, ti_current,
- random_number_BH_feedback);
+ /* Number of energy injections that have reached this gas particle */
+ int num_of_energy_inj_received_by_gas = 0;
- /* Are we lucky? */
- if (rand < prob) {
+ /* Find out how many rays (= energy injections) this gas particle
+ * has received */
+ for (int i = 0; i < num_energy_injections_per_BH; i++) {
+ if (pj->id == bi->rays[i].id_min_length)
+ num_of_energy_inj_received_by_gas++;
+ }
- /* Compute new energy per unit mass of this particle */
+ /* If the number of received rays is non-zero, inject
+ * AGN energy in thermal form */
+ if (num_of_energy_inj_received_by_gas > 0) {
+
+ /* Compute new energy per unit mass of this particle
+ * The energy the particle receives is proportional to the number of rays
+ * (num_of_energy_inj_received_by_gas) to which the particle was found to
+ * be closest. */
const double u_init = hydro_get_physical_internal_energy(pj, xpj, cosmo);
- const float delta_u = bi->to_distribute.AGN_delta_u;
+ const float delta_u = bi->to_distribute.AGN_delta_u *
+ (float)num_of_energy_inj_received_by_gas;
const double u_new = u_init + delta_u;
hydro_set_physical_internal_energy(pj, xpj, cosmo, u_new);
@@ -575,6 +636,9 @@ runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
/* Impose maximal viscosity */
hydro_diffusive_feedback_reset(pj);
+ /* Update cooling properties. */
+ cooling_update_feedback_particle(xpj);
+
/* Store the feedback energy */
const double delta_energy = delta_u * hydro_get_mass(pj);
tracers_after_black_holes_feedback(xpj, with_cosmology, cosmo->a, time,
@@ -583,7 +647,7 @@ runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
/* message( */
/* "We did some AGN heating! id %llu BH id %llu probability " */
/* " %.5e random_num %.5e du %.5e du/ini %.5e", */
- /* pj->id, bi->id, prob, rand, delta_u, delta_u / u_init); */
+ /* pj->id, bi->id, 0.f, 0.f, delta_u, delta_u / u_init); */
/* Synchronize the particle on the timeline */
timestep_sync_part(pj);
diff --git a/src/black_holes/EAGLE/black_holes_io.h b/src/black_holes/EAGLE/black_holes_io.h
index 2ad9232d2756e4f6ff9900e969063510f89fb680..eb6bcd4080560a8824199c0755fb413e7487dbe7 100644
--- a/src/black_holes/EAGLE/black_holes_io.h
+++ b/src/black_holes/EAGLE/black_holes_io.h
@@ -152,7 +152,7 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
int with_cosmology) {
/* Say how much we want to write */
- *num_fields = 38;
+ *num_fields = 42;
/* List what we want to write */
list[0] = io_make_output_field_convert_bpart(
@@ -367,32 +367,59 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
"particles in the most recent feedback event.");
list[33] = io_make_output_field(
- "LastRepositionVelocities", FLOAT, 1, UNIT_CONV_SPEED, 0.f, bparts,
- last_repos_vel,
- "Physical speeds at which the black holes repositioned most recently. "
- "This is 0 for black holes that have never repositioned, or if the "
- "simulation has been run without prescribed repositioning speed.");
+ "NumberOfHeatingEvents", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
+ AGN_number_of_energy_injections,
+ "Integer number of (thermal) energy injections the black hole has had "
+ "so far");
list[34] = io_make_output_field(
+ "NumberOfAGNEvents", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
+ AGN_number_of_AGN_events,
+ "Integer number of AGN events the black hole has had so far"
+ " (the number of times the BH did AGN feedback)");
+
+ if (with_cosmology) {
+ list[35] = io_make_output_field(
+ "LastAGNFeedbackScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
+ bparts, last_AGN_event_scale_factor,
+ "Scale-factors at which the black holes last had an AGN event.");
+ } else {
+ list[35] = io_make_output_field(
+ "LastAGNFeedbackTimes", FLOAT, 1, UNIT_CONV_TIME, 0.f, bparts,
+ last_AGN_event_time,
+ "Times at which the black holes last had an AGN event.");
+ }
+
+ list[36] = io_make_output_field(
+ "AccretionLimitedTimeSteps", FLOAT, 1, UNIT_CONV_TIME, 0.f, bparts,
+ dt_heat, "Accretion-limited time-steps of black holes.");
+
+ list[37] = io_make_output_field(
+ "AGNTotalInjectedEnergies", FLOAT, 1, UNIT_CONV_ENERGY, 0.f, bparts,
+ AGN_cumulative_energy,
+ "Total (cumulative) physical energies injected into gas particles "
+ "in AGN feedback.");
+
+ list[38] = io_make_output_field(
"AccretionBoostFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
accretion_boost_factor,
"Multiplicative factors by which the Bondi-Hoyle-Lyttleton accretion "
"rates have been increased by the density-dependent Booth & Schaye "
"(2009) accretion model.");
- list[35] = io_make_output_field_convert_bpart(
+ list[39] = io_make_output_field_convert_bpart(
"GasTemperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, bparts,
convert_bpart_gas_temperatures,
"Temperature of the gas surrounding the black holes.");
- list[36] = io_make_output_field(
+ list[40] = io_make_output_field(
"EnergyReservoirThresholds", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
num_ngbs_to_heat,
"Minimum energy reservoir required for the black holes to do feedback, "
"expressed in units of the (constant) target heating temperature "
"increase.");
- list[37] = io_make_output_field(
+ list[41] = io_make_output_field(
"EddingtonFractions", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
eddington_fraction,
"Accretion rates of black holes in units of their Eddington rates. "
diff --git a/src/black_holes/EAGLE/black_holes_part.h b/src/black_holes/EAGLE/black_holes_part.h
index 0be3e8fc86f81af77d103a8f88e674e95a89cb89..6fcbdfb2ec3a6ebd0dd2a1109873eeedbf07285a 100644
--- a/src/black_holes/EAGLE/black_holes_part.h
+++ b/src/black_holes/EAGLE/black_holes_part.h
@@ -19,6 +19,9 @@
#ifndef SWIFT_EAGLE_BLACK_HOLE_PART_H
#define SWIFT_EAGLE_BLACK_HOLE_PART_H
+/*! The total number of rays used in AGN feedback */
+#define eagle_blackhole_number_of_rays 50
+
#include "black_holes_struct.h"
#include "chemistry_struct.h"
#include "timeline.h"
@@ -169,12 +172,36 @@ struct bpart {
/*! Instantaneous temperature increase for feedback */
float AGN_delta_T;
- /*! Instantaneous energy reservoir threshold (num-to-heat) */
- float num_ngbs_to_heat;
-
/*! Eddington fractions */
float eddington_fraction;
+ /*! Integer (cumulative) number of energy injections in AGN feedback. At a
+ * given time-step, an AGN-active BH may produce multiple energy injections.
+ * The number of energy injections is equal to or more than the number of
+ * particles heated by the BH during this time-step. */
+ int AGN_number_of_energy_injections;
+
+ /*! Integer (cumulative) number of AGN events. If a BH does feedback at a
+ * given time-step, the number of its AGN events is incremented by 1. Each
+ * AGN event may have multiple energy injections. */
+ int AGN_number_of_AGN_events;
+
+ /* Total energy injected into the gas in AGN feedback by this BH */
+ float AGN_cumulative_energy;
+
+ /*! BH accretion-limited time-step */
+ float dt_heat;
+
+ /*! Union for the last AGN event time and the last AGN event scale factor */
+ union {
+
+ /*! Last AGN event time */
+ float last_AGN_event_time;
+
+ /*! Last AGN event scale-factor */
+ float last_AGN_event_scale_factor;
+ };
+
/*! Union for the last high Eddington ratio point in time */
union {
@@ -208,8 +235,8 @@ struct bpart {
/*! Properties used in the feedback loop to distribute to gas neighbours. */
struct {
- /*! Probability of heating neighbouring gas particles for AGN feedback */
- float AGN_heating_probability;
+ /*! Number of energy injections per time-step */
+ int AGN_number_of_energy_injections;
/*! Change in energy from SNII feedback energy injection */
float AGN_delta_u;
@@ -236,6 +263,9 @@ struct bpart {
/*! Black holes merger information (e.g. merging ID) */
struct black_holes_bpart_data merger_data;
+ /*! Isotropic AGN feedback information */
+ struct ray_data rays[eagle_blackhole_number_of_rays];
+
#ifdef SWIFT_DEBUG_CHECKS
/* Time of the last drift */
diff --git a/src/black_holes/EAGLE/black_holes_properties.h b/src/black_holes/EAGLE/black_holes_properties.h
index afb2ca00bc8e0b1d5ee906ecbd95a9f6a482a934..ba6a07b0bb1d3abf0568dba2ec0fd19739cf7430 100644
--- a/src/black_holes/EAGLE/black_holes_properties.h
+++ b/src/black_holes/EAGLE/black_holes_properties.h
@@ -22,6 +22,13 @@
#include "chemistry.h"
#include "hydro_properties.h"
+#include <string.h>
+
+enum AGN_feedback_models {
+ AGN_isotropic_model, /*< Isotropic model of AGN feedback */
+ AGN_minimum_distance_model /*< Minimum-distance model of AGN feedback */
+};
+
/**
* @brief Properties of black holes and AGN feedback in the EAGEL model.
*/
@@ -106,6 +113,12 @@ struct black_holes_props {
/* ---- Properties of the feedback model ------- */
+ /*! AGN feedback model: isotropic or minimum distance */
+ enum AGN_feedback_models feedback_model;
+
+ /*! Is the AGN feedback model deterministic or stochastic? */
+ int AGN_deterministic;
+
/*! Feedback coupling efficiency of the black holes. */
float epsilon_f;
@@ -348,6 +361,21 @@ 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)
+ 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",
+ temp);
+
+ bp->AGN_deterministic =
+ parser_get_param_int(params, "EAGLEAGN:AGN_use_deterministic_feedback");
+
bp->epsilon_f =
parser_get_param_float(params, "EAGLEAGN:coupling_efficiency");
diff --git a/src/random.h b/src/random.h
index 79c478897fe64d8f418ca47be758ab3a8c5f39b3..1ecb2683f558df319bdbb216511aa4ab14a1d9b5 100644
--- a/src/random.h
+++ b/src/random.h
@@ -59,6 +59,8 @@ enum random_number_type {
random_number_stellar_feedback_3 = 9762399103LL,
random_number_isotropic_SNII_feedback_ray_theta = 3298327511LL,
random_number_isotropic_SNII_feedback_ray_phi = 6311114273LL,
+ random_number_isotropic_AGN_feedback_ray_theta = 8899891613LL,
+ random_number_isotropic_AGN_feedback_ray_phi = 10594523341LL,
random_number_stellar_enrichment = 2936881973LL,
random_number_BH_feedback = 1640531371LL,
random_number_BH_swallow = 4947009007LL,
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index 821725ed5b3183d60fa12dec0f8c17e982840317..8ef7214d80f375cfdd6ace1fcd521c085231e256 100644
--- a/src/runner_ghost.c
+++ b/src/runner_ghost.c
@@ -877,9 +877,10 @@ void runner_do_black_holes_swallow_ghost(struct runner *r, struct cell *c,
e->physical_constants, e->cosmology, dt);
/* Compute variables required for the feedback loop */
- black_holes_prepare_feedback(
- bp, e->black_holes_properties, e->physical_constants, e->cosmology,
- e->cooling_func, e->entropy_floor, e->time, with_cosmology, dt);
+ black_holes_prepare_feedback(bp, e->black_holes_properties,
+ e->physical_constants, e->cosmology,
+ e->cooling_func, e->entropy_floor, e->time,
+ with_cosmology, dt, e->ti_current);
}
}
}