/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see .
*
******************************************************************************/
#ifndef SWIFT_COLIBRE_BLACK_HOLES_H
#define SWIFT_COLIBRE_BLACK_HOLES_H
/* Local includes */
#include "black_holes_properties.h"
#include "black_holes_struct.h"
#include "cooling_properties.h"
#include "cosmology.h"
#include "dimension.h"
#include "dust.h"
#include "gravity.h"
#include "kernel_hydro.h"
#include "minmax.h"
#include "physical_constants.h"
#include "random.h"
#include "rays.h"
/* Standard includes */
#include
#include
/**
* @brief Computes the AGN heating temperature in the varialbe AGN dT model.
*
* @param bp Pointer to the s-particle data.
* @param props The properties of the black hole scheme.
*/
__attribute__((always_inline)) INLINE static float black_hole_feedback_delta_T(
const struct bpart* const bp, const struct black_holes_props* props) {
/* The heating temperature reached at `props->AGN_delta_T_Mass' */
double delta_T = props->AGN_delta_T_pivot;
/* Scale pivot heating temperature of AGN in proportion to BH mass */
const double AGN_dT_factor =
min(bp->subgrid_mass / props->AGN_delta_T_Mass, 1.);
delta_T *= AGN_dT_factor;
/* Ensure the new value is greater than the floor value */
delta_T = max(delta_T, props->AGN_delta_T_floor);
return delta_T;
}
/**
* @brief Computes the time-step of a given black hole particle.
*
* @param bp Pointer to the s-particle data.
* @param props The properties of the black hole scheme.
* @param constants The physical constants (in internal units).
*/
__attribute__((always_inline)) INLINE static float black_holes_compute_timestep(
const struct bpart* const bp, const struct black_holes_props* props,
const struct phys_const* constants, const struct cosmology* cosmo) {
/* Do something if and only if the accretion rate is non-zero */
if (bp->accretion_rate > 0.f) {
/* Gather some physical constants (all in internal units) */
const double c = constants->const_speed_light_c;
/* Compute instantaneous energy supply rate to the BH energy reservoir
* which is proportional to the BH mass accretion rate */
const double Energy_rate =
props->epsilon_r * props->epsilon_f * bp->accretion_rate * c * c;
/* Compute average heating energy in AGN feedback */
/* Average particle mass in BH's kernel */
const double mean_ngb_mass = bp->ngb_mass / ((double)bp->num_ngbs);
/* Without multiplying by mean_ngb_mass we'd get energy per unit mass */
/* Compute AGN heating temperature */
const double dT_AGN = black_hole_feedback_delta_T(bp, props);
/* Energy in one AGN heating event */
const double E_heat = dT_AGN * props->temp_to_u_factor * mean_ngb_mass;
/* Compute average time between energy injections for the given accretion
* 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 * props->num_ngbs_to_heat / Energy_rate;
/* The new timestep of the BH cannot be smaller than the miminum allowed
* time-step */
const double bh_timestep = max(dt_heat, props->time_step_min);
return bh_timestep;
}
return FLT_MAX;
}
/**
* @brief Initialises the b-particles for the first time
*
* This function is called only once just after the ICs have been
* read in to do some conversions.
*
* @param bp The particle to act upon
* @param props The properties of the black holes model.
*/
__attribute__((always_inline)) INLINE static void black_holes_first_init_bpart(
struct bpart* bp, const struct black_holes_props* props) {
bp->time_bin = 0;
if (props->use_subgrid_mass_from_ics == 0)
bp->subgrid_mass = bp->mass;
else if (props->with_subgrid_mass_check && bp->subgrid_mass <= 0)
error(
"Black hole %lld has a subgrid mass of %f (internal units).\n"
"If this is because the ICs do not contain a 'SubgridMass' data "
"set, you should set the parameter "
"'EAGLEAGN:use_subgrid_mass_from_ics' to 0 to initialize the "
"black hole subgrid masses to the corresponding dynamical masses.\n"
"If the subgrid mass is intentionally set to this value, you can "
"disable this error by setting 'EAGLEAGN:with_subgrid_mass_check' "
"to 0.",
bp->id, bp->subgrid_mass);
bp->total_accreted_mass = 0.f;
bp->accretion_rate = 0.f;
bp->formation_time = -1.f;
bp->energy_reservoir = 0.f;
bp->cumulative_number_seeds = 1;
bp->number_of_mergers = 0;
bp->number_of_gas_swallows = 0;
bp->number_of_direct_gas_swallows = 0;
bp->cumulative_mass_lost_to_gw = 0.f;
bp->number_of_repositions = 0;
bp->number_of_reposition_attempts = 0;
bp->number_of_time_steps = 0;
bp->last_high_Eddington_fraction_scale_factor = -1.f;
bp->last_minor_merger_time = -1.;
bp->last_major_merger_time = -1.;
bp->last_AGN_event_time = -1.;
bp->swallowed_angular_momentum[0] = 0.f;
bp->swallowed_angular_momentum[1] = 0.f;
bp->swallowed_angular_momentum[2] = 0.f;
bp->accreted_angular_momentum[0] = 0.f;
bp->accreted_angular_momentum[1] = 0.f;
bp->accreted_angular_momentum[2] = 0.f;
bp->dt_heat = 0.f;
bp->AGN_number_of_AGN_events = 0;
bp->AGN_number_of_energy_injections = 0;
black_holes_mark_bpart_as_not_swallowed(&bp->merger_data);
}
/**
* @brief Prepares a b-particle for its interactions
*
* @param bp The particle to act upon
*/
__attribute__((always_inline)) INLINE static void black_holes_init_bpart(
struct bpart* bp) {
#ifdef DEBUG_INTERACTIONS_BLACK_HOLES
for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
bp->ids_ngbs_density[i] = -1;
bp->num_ngb_density = 0;
#endif
bp->density.wcount = 0.f;
bp->density.wcount_dh = 0.f;
bp->rho_gas = 0.f;
bp->sound_speed_gas = 0.f;
bp->internal_energy_gas = 0.f;
bp->velocity_gas[0] = 0.f;
bp->velocity_gas[1] = 0.f;
bp->velocity_gas[2] = 0.f;
bp->spec_angular_momentum_gas[0] = 0.f;
bp->spec_angular_momentum_gas[1] = 0.f;
bp->spec_angular_momentum_gas[2] = 0.f;
bp->curl_v_gas[0] = 0.f;
bp->curl_v_gas[1] = 0.f;
bp->curl_v_gas[2] = 0.f;
bp->velocity_dispersion_gas = 0.f;
bp->ngb_mass = 0.f;
bp->num_ngbs = 0;
bp->reposition.delta_x[0] = -FLT_MAX;
bp->reposition.delta_x[1] = -FLT_MAX;
bp->reposition.delta_x[2] = -FLT_MAX;
bp->reposition.min_potential = FLT_MAX;
bp->reposition.potential = FLT_MAX;
bp->accretion_rate = 0.f; /* Optionally accumulated ngb-by-ngb */
bp->f_visc = FLT_MAX;
bp->mass_at_start_of_step = bp->mass; /* bp->mass may grow in nibbling mode */
bp->min_ngb_time_bin = num_time_bins + 1;
/* Reset the rays carried by this BH */
ray_init(bp->rays, colibre_blackhole_number_of_rays);
}
/**
* @brief Predict additional particle fields forward in time when drifting
*
* The fields do not get predicted but we move the BH to its new position
* if a new one was calculated in the repositioning loop.
*
* @param bp The particle
* @param dt_drift The drift time-step for positions.
*/
__attribute__((always_inline)) INLINE static void black_holes_predict_extra(
struct bpart* restrict bp, float dt_drift) {
/* Are we doing some repositioning? */
if (bp->reposition.min_potential != FLT_MAX) {
#ifdef SWIFT_DEBUG_CHECKS
if (bp->reposition.delta_x[0] == -FLT_MAX ||
bp->reposition.delta_x[1] == -FLT_MAX ||
bp->reposition.delta_x[2] == -FLT_MAX) {
error("Something went wrong with the new repositioning position");
}
const double dx = bp->reposition.delta_x[0];
const double dy = bp->reposition.delta_x[1];
const double dz = bp->reposition.delta_x[2];
const double d = sqrt(dx * dx + dy * dy + dz * dz);
if (d > 1.01 * kernel_gamma * bp->h)
error("Repositioning BH beyond the kernel support!");
#endif
/* Move the black hole */
bp->x[0] += bp->reposition.delta_x[0];
bp->x[1] += bp->reposition.delta_x[1];
bp->x[2] += bp->reposition.delta_x[2];
/* Move its gravity properties as well */
bp->gpart->x[0] += bp->reposition.delta_x[0];
bp->gpart->x[1] += bp->reposition.delta_x[1];
bp->gpart->x[2] += bp->reposition.delta_x[2];
/* Store the delta position */
bp->x_diff[0] -= bp->reposition.delta_x[0];
bp->x_diff[1] -= bp->reposition.delta_x[1];
bp->x_diff[2] -= bp->reposition.delta_x[2];
/* Reset the reposition variables */
bp->reposition.delta_x[0] = -FLT_MAX;
bp->reposition.delta_x[1] = -FLT_MAX;
bp->reposition.delta_x[2] = -FLT_MAX;
bp->reposition.min_potential = FLT_MAX;
/* Count the jump */
bp->number_of_repositions++;
}
}
/**
* @brief Sets the values to be predicted in the drifts to their values at a
* kick time
*
* @param bp The particle.
*/
__attribute__((always_inline)) INLINE static void
black_holes_reset_predicted_values(struct bpart* bp) {}
/**
* @brief Kick the additional variables
*
* @param bp The particle to act upon
* @param dt The time-step for this kick
*/
__attribute__((always_inline)) INLINE static void black_holes_kick_extra(
struct bpart* bp, float dt) {}
/**
* @brief Finishes the calculation of density on black holes
*
* @param bp The particle to act upon
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void black_holes_end_density(
struct bpart* bp, const struct cosmology* cosmo) {
/* Some smoothing length multiples. */
const float h = bp->h;
const float h_inv = 1.0f / h; /* 1/h */
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
/* --- Finish the calculation by inserting the missing h factors --- */
bp->density.wcount *= h_inv_dim;
bp->density.wcount_dh *= h_inv_dim_plus_one;
bp->rho_gas *= h_inv_dim;
const float rho_inv = 1.f / bp->rho_gas;
/* For the following, we also have to undo the mass smoothing
* (N.B.: bp->velocity_gas is in BH frame, in internal units). */
bp->sound_speed_gas *= h_inv_dim * rho_inv;
bp->internal_energy_gas *= h_inv_dim * rho_inv;
bp->velocity_gas[0] *= h_inv_dim * rho_inv;
bp->velocity_gas[1] *= h_inv_dim * rho_inv;
bp->velocity_gas[2] *= h_inv_dim * rho_inv;
bp->velocity_dispersion_gas *= h_inv_dim * rho_inv;
bp->spec_angular_momentum_gas[0] *= h_inv_dim * rho_inv;
bp->spec_angular_momentum_gas[1] *= h_inv_dim * rho_inv;
bp->spec_angular_momentum_gas[2] *= h_inv_dim * rho_inv;
/* ... and for the curl, we also need to divide by an extra h factor */
bp->curl_v_gas[0] *= h_inv_dim_plus_one * rho_inv;
bp->curl_v_gas[1] *= h_inv_dim_plus_one * rho_inv;
bp->curl_v_gas[2] *= h_inv_dim_plus_one * rho_inv;
/* Calculate circular velocity at the smoothing radius from specific
* angular momentum (extra h_inv) */
bp->circular_velocity_gas[0] = bp->spec_angular_momentum_gas[0] * h_inv;
bp->circular_velocity_gas[1] = bp->spec_angular_momentum_gas[1] * h_inv;
bp->circular_velocity_gas[2] = bp->spec_angular_momentum_gas[2] * h_inv;
/* Calculate (actual) gas velocity dispersion. Currently, the variable
* 'velocity_dispersion_gas' holds instead. */
const double speed_gas2 = bp->velocity_gas[0] * bp->velocity_gas[0] +
bp->velocity_gas[1] * bp->velocity_gas[1] +
bp->velocity_gas[2] * bp->velocity_gas[2];
bp->velocity_dispersion_gas -= speed_gas2;
bp->velocity_dispersion_gas = sqrt(fabs(bp->velocity_dispersion_gas));
}
/**
* @brief Sets all particle fields to sensible values when the #spart has 0
* ngbs.
*
* @param bp The particle to act upon
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void
black_holes_bpart_has_no_neighbours(struct bpart* bp,
const struct cosmology* cosmo) {
/* Some smoothing length multiples. */
const float h = bp->h;
const float h_inv = 1.0f / h; /* 1/h */
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
/* Re-set problematic values */
bp->density.wcount = kernel_root * h_inv_dim;
bp->density.wcount_dh = 0.f;
bp->velocity_gas[0] = FLT_MAX;
bp->velocity_gas[1] = FLT_MAX;
bp->velocity_gas[2] = FLT_MAX;
bp->velocity_dispersion_gas = FLT_MAX;
bp->curl_v_gas[0] = FLT_MAX;
bp->curl_v_gas[1] = FLT_MAX;
bp->curl_v_gas[2] = FLT_MAX;
bp->internal_energy_gas = -FLT_MAX;
}
/**
* @brief Return the current instantaneous accretion rate of the BH.
*
* @param bp the #bpart.
*/
__attribute__((always_inline)) INLINE static double
black_holes_get_accretion_rate(const struct bpart* bp) {
return bp->accretion_rate;
}
/**
* @brief Return the total accreted gas mass of this BH.
*
* @param bp the #bpart.
*/
__attribute__((always_inline)) INLINE static double
black_holes_get_accreted_mass(const struct bpart* bp) {
return bp->total_accreted_mass;
}
/**
* @brief Return the subgrid mass of this BH.
*
* @param bp the #bpart.
*/
__attribute__((always_inline)) INLINE static double
black_holes_get_subgrid_mass(const struct bpart* bp) {
return bp->subgrid_mass;
}
/**
* @brief Return the current bolometric luminosity of the BH.
*
* @param bp the #bpart.
*/
__attribute__((always_inline)) INLINE static double
black_holes_get_bolometric_luminosity(const struct bpart* bp,
const struct phys_const* constants) {
return 0.f;
}
/**
* @brief Return the current kinetic jet power of the BH.
*
* No jet in this model so we return 0.
*
* @param bp the #bpart.
*/
__attribute__((always_inline)) INLINE static double black_holes_get_jet_power(
const struct bpart* bp, const struct phys_const* constants) {
return 0.;
}
/**
* @brief Return the cumulative mass lost to GWs of this BH.
*
* @param bp the #bpart.
*/
__attribute__((always_inline)) INLINE static double
black_holes_get_GW_mass_loss(const struct bpart* bp) {
return bp->cumulative_mass_lost_to_gw;
}
/**
* @brief Update the properties of a black hole particles by swallowing
* a gas particle.
*
* @param bp The #bpart to update.
* @param p The #part that is swallowed.
* @param xp The #xpart that is swallowed.
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void black_holes_swallow_part(
struct bpart* bp, const struct part* p, const struct xpart* xp,
const struct cosmology* cosmo) {
/* Get the current dynamical masses */
const float gas_mass = hydro_get_mass(p);
const float BH_mass = bp->mass;
/* Increase the dynamical mass of the BH. */
bp->mass += gas_mass;
bp->gpart->mass += gas_mass;
/* Physical velocity difference between the particles */
const float dv[3] = {(bp->v[0] - p->v[0]) * cosmo->a_inv,
(bp->v[1] - p->v[1]) * cosmo->a_inv,
(bp->v[2] - p->v[2]) * cosmo->a_inv};
/* Physical distance between the particles */
const float dx[3] = {(bp->x[0] - p->x[0]) * cosmo->a,
(bp->x[1] - p->x[1]) * cosmo->a,
(bp->x[2] - p->x[2]) * cosmo->a};
/* Collect the swallowed angular momentum */
bp->swallowed_angular_momentum[0] +=
gas_mass * (dx[1] * dv[2] - dx[2] * dv[1]);
bp->swallowed_angular_momentum[1] +=
gas_mass * (dx[2] * dv[0] - dx[0] * dv[2]);
bp->swallowed_angular_momentum[2] +=
gas_mass * (dx[0] * dv[1] - dx[1] * dv[0]);
/* Update the BH momentum */
const float BH_mom[3] = {BH_mass * bp->v[0] + gas_mass * p->v[0],
BH_mass * bp->v[1] + gas_mass * p->v[1],
BH_mass * bp->v[2] + gas_mass * p->v[2]};
bp->v[0] = BH_mom[0] / bp->mass;
bp->v[1] = BH_mom[1] / bp->mass;
bp->v[2] = BH_mom[2] / bp->mass;
bp->gpart->v_full[0] = bp->v[0];
bp->gpart->v_full[1] = bp->v[1];
bp->gpart->v_full[2] = bp->v[2];
const float dr = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
message(
"BH %lld swallowing gas particle %lld "
"(Delta_v = [%f, %f, %f] U_V, "
"Delta_x = [%f, %f, %f] U_L, "
"Delta_v_rad = %f)",
bp->id, p->id, -dv[0], -dv[1], -dv[2], -dx[0], -dx[1], -dx[2],
(dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]) / dr);
/* Update the BH metal masses */
struct chemistry_bpart_data* bp_chem = &bp->chemistry_data;
const struct chemistry_part_data* p_chem = &p->chemistry_data;
chemistry_add_part_to_bpart(bp_chem, p_chem, gas_mass);
/* Update the BH dust masses */
struct dust_bpart_data* bp_dust = &bp->dust_data;
const struct dust_part_data* p_dust = &p->dust_data;
dust_add_part_to_bpart(bp_dust, p_dust, gas_mass);
/* This BH swallowed a gas particle */
bp->number_of_gas_swallows++;
bp->number_of_direct_gas_swallows++;
/* This BH lost a neighbour */
bp->num_ngbs--;
bp->ngb_mass -= gas_mass;
/* The ray(s) should not point to the no-longer existing particle */
ray_reset_part_id(bp->rays, colibre_blackhole_number_of_rays, p->id);
}
/**
* @brief Update the properties of a black hole particles by swallowing
* a BH particle.
*
* @param bpi The #bpart to update.
* @param bpj The #bpart that is swallowed.
* @param cosmo The current cosmological model.
* @param time Time since the start of the simulation (non-cosmo mode).
* @param with_cosmology Are we running with cosmology?
* @param props The properties of the black hole scheme.
*/
__attribute__((always_inline)) INLINE static void black_holes_swallow_bpart(
struct bpart* bpi, const struct bpart* bpj, const struct cosmology* cosmo,
const double time, const int with_cosmology,
const struct black_holes_props* props, const struct phys_const* constants) {
/* Get the current masses, dynamical and subgrid */
const float bpi_dyn_mass = bpi->mass;
const float bpj_dyn_mass = bpj->mass;
const float bpi_sub_mass = bpi->subgrid_mass;
const float bpj_sub_mass = bpj->subgrid_mass;
/* Is this merger ratio above the threshold for recording? */
const double merger_ratio = bpj->subgrid_mass / bpi->subgrid_mass;
if (merger_ratio > props->major_merger_threshold) {
if (with_cosmology) {
bpi->last_major_merger_scale_factor = cosmo->a;
} else {
bpi->last_major_merger_time = time;
}
} else if (merger_ratio > props->minor_merger_threshold) {
if (with_cosmology) {
bpi->last_minor_merger_scale_factor = cosmo->a;
} else {
bpi->last_minor_merger_time = time;
}
}
/* Increase the masses of the BH. */
bpi->mass += bpj->mass;
bpi->gpart->mass += bpj->mass;
bpi->subgrid_mass += bpj->subgrid_mass;
/* Mass loss into Gravitational waves
*
* Expression from Barausse et al. 2012 ApJ 758,
* eq. 18, 24 & 26 */
const double M = (bpi_sub_mass + bpj_sub_mass);
const double eta = bpi_sub_mass * bpj_sub_mass / (M * M);
const double mass_frac_loss = 0.0572 * eta + 0.54352 * eta * eta;
const float mass_loss_to_gw = mass_frac_loss * M;
message(
"BH %lld (mass %f/%f) swallowing BH %lld (mass %f/%f) - GW mass loss: %f",
bpi->id, bpi->mass, bpi->subgrid_mass, bpj->id, bpj->mass,
bpj->subgrid_mass, mass_loss_to_gw);
bpi->mass -= mass_loss_to_gw;
bpi->gpart->mass -= mass_loss_to_gw;
bpi->subgrid_mass -= mass_loss_to_gw;
bpi->cumulative_mass_lost_to_gw += bpj->cumulative_mass_lost_to_gw;
bpi->cumulative_mass_lost_to_gw += mass_loss_to_gw;
/* Collect the swallowed angular momentum */
bpi->swallowed_angular_momentum[0] += bpj->swallowed_angular_momentum[0];
bpi->swallowed_angular_momentum[1] += bpj->swallowed_angular_momentum[1];
bpi->swallowed_angular_momentum[2] += bpj->swallowed_angular_momentum[2];
/* Update the BH momentum */
const float BH_mom[3] = {bpi_dyn_mass * bpi->v[0] + bpj_dyn_mass * bpj->v[0],
bpi_dyn_mass * bpi->v[1] + bpj_dyn_mass * bpj->v[1],
bpi_dyn_mass * bpi->v[2] + bpj_dyn_mass * bpj->v[2]};
bpi->v[0] = BH_mom[0] / bpi->mass;
bpi->v[1] = BH_mom[1] / bpi->mass;
bpi->v[2] = BH_mom[2] / bpi->mass;
bpi->gpart->v_full[0] = bpi->v[0];
bpi->gpart->v_full[1] = bpi->v[1];
bpi->gpart->v_full[2] = bpi->v[2];
/* Update the BH metal masses */
struct chemistry_bpart_data* bpi_chem = &bpi->chemistry_data;
const struct chemistry_bpart_data* bpj_chem = &bpj->chemistry_data;
chemistry_add_bpart_to_bpart(bpi_chem, bpj_chem);
/* Update the BH dust masses */
struct dust_bpart_data* bpi_dust = &bpi->dust_data;
const struct dust_bpart_data* bpj_dust = &bpj->dust_data;
dust_add_bpart_to_bpart(bpi_dust, bpj_dust);
/* Update the energy reservoir */
bpi->energy_reservoir += bpj->energy_reservoir;
/* Add up all the BH seeds */
bpi->cumulative_number_seeds += bpj->cumulative_number_seeds;
/* Add up all the gas particles we swallowed */
bpi->number_of_gas_swallows += bpj->number_of_gas_swallows;
/* Add the subgrid angular momentum that we swallowed */
bpi->accreted_angular_momentum[0] += bpj->accreted_angular_momentum[0];
bpi->accreted_angular_momentum[1] += bpj->accreted_angular_momentum[1];
bpi->accreted_angular_momentum[2] += bpj->accreted_angular_momentum[2];
/* We had another merger */
bpi->number_of_mergers++;
}
/**
* @brief Compute the accretion rate of the black hole and all the quantites
* required for the feedback loop.
*
* @param bp The black hole particle.
* @param props The properties of the black hole scheme.
* @param constants The physical constants (in internal units).
* @param cosmo The cosmological model.
* @param cooling Properties of the cooling model.
* @param floor_props Properties of the entropy fllor.
* @param time Time since the start of the simulation (non-cosmo mode).
* @param with_cosmology Are we running with cosmology?
* @param dt The time-step size (in physical internal units).
* @param ti_begin Integer time value at the beginning of timestep
*/
__attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
struct bpart* restrict bp, const struct black_holes_props* props,
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 integertime_t ti_begin) {
/* Record that the black hole has another active time step */
bp->number_of_time_steps++;
if (dt == 0. || bp->rho_gas == 0.) return;
/* Gather some physical constants (all in internal units) */
const double G = constants->const_newton_G;
const double c = constants->const_speed_light_c;
const double proton_mass = constants->const_proton_mass;
const double sigma_Thomson = constants->const_thomson_cross_section;
/* (Subgrid) mass of the BH (internal units) */
const double BH_mass = bp->subgrid_mass;
/* Convert the quantities we gathered to physical frame (all internal units).
* Note: for the velocities this means peculiar velocities */
const double gas_c_phys = bp->sound_speed_gas * cosmo->a_factor_sound_speed;
const double gas_c_phys2 = gas_c_phys * gas_c_phys;
const double gas_v_circular[3] = {
bp->circular_velocity_gas[0] * cosmo->a_inv,
bp->circular_velocity_gas[1] * cosmo->a_inv,
bp->circular_velocity_gas[2] * cosmo->a_inv};
/* Norm of the circular velocity of the gas around the BH */
const double tangential_velocity2 = gas_v_circular[0] * gas_v_circular[0] +
gas_v_circular[1] * gas_v_circular[1] +
gas_v_circular[2] * gas_v_circular[2];
const double tangential_velocity = sqrt(tangential_velocity2);
/* We can now compute the accretion rate (internal units) */
double accr_rate;
if (props->use_multi_phase_bondi) {
/* In this case, we are in 'multi-phase-Bondi' mode -- otherwise,
* the accretion_rate is still zero (was initialised to this) */
const float hi_inv = 1.f / bp->h;
const float hi_inv_dim = pow_dimension(hi_inv); /* 1/h^d */
accr_rate = bp->accretion_rate *
(4. * M_PI * G * G * BH_mass * BH_mass * hi_inv_dim);
} else {
/* Standard approach: compute accretion rate for all gas simultaneously.
*
* Convert the quantities we gathered to physical frame (all internal
* units). Note: velocities are already in black hole frame. */
const double gas_rho_phys = bp->rho_gas * cosmo->a3_inv;
const double gas_v_phys[3] = {bp->velocity_gas[0] * cosmo->a_inv,
bp->velocity_gas[1] * cosmo->a_inv,
bp->velocity_gas[2] * cosmo->a_inv};
const double gas_v_norm2 = gas_v_phys[0] * gas_v_phys[0] +
gas_v_phys[1] * gas_v_phys[1] +
gas_v_phys[2] * gas_v_phys[2];
/* In the Bondi-Hoyle-Lyttleton formula, the bulk flow of gas is
* added to the sound speed in quadrature. Treated separately (below)
* in the Krumholz et al. (2006) prescription */
const double denominator2 =
props->use_krumholz ? gas_c_phys2 : gas_v_norm2 + gas_c_phys2;
#ifdef SWIFT_DEBUG_CHECKS
/* Make sure that the denominator is strictly positive */
if (denominator2 <= 0)
error(
"Invalid denominator for BH particle %lld in Bondi rate "
"calculation.",
bp->id);
#endif
const double denominator_inv = 1. / sqrt(denominator2);
accr_rate = 4. * M_PI * G * G * BH_mass * BH_mass * gas_rho_phys *
denominator_inv * denominator_inv * denominator_inv;
if (props->use_krumholz) {
/* Compute the additional correction factors from Krumholz+06,
* accounting for bulk flow and turbulence of ambient gas. */
const double lambda = 1.1;
const double gas_v_dispersion =
bp->velocity_dispersion_gas * cosmo->a_inv;
const double mach_turb = gas_v_dispersion / gas_c_phys;
const double mach_bulk = sqrt(gas_v_norm2) / gas_c_phys;
const double mach2 = mach_turb * mach_turb + mach_bulk * mach_bulk;
const double m1 = 1. + mach2;
const double mach_factor =
sqrt((lambda * lambda + mach2) / (m1 * m1 * m1 * m1));
accr_rate *= mach_factor;
}
if (props->with_krumholz_vorticity) {
/* Change the accretion rate to equation (3) of Krumholz et al. (2006)
* by adding a vorticity-dependent term in inverse quadrature */
/* Convert curl to vorticity in physical units */
const double gas_curlv_phys[3] = {bp->curl_v_gas[0] * cosmo->a2_inv,
bp->curl_v_gas[1] * cosmo->a2_inv,
bp->curl_v_gas[2] * cosmo->a2_inv};
const double gas_vorticity = sqrt(gas_curlv_phys[0] * gas_curlv_phys[0] +
gas_curlv_phys[1] * gas_curlv_phys[1] +
gas_curlv_phys[2] * gas_curlv_phys[2]);
const double Bondi_radius = G * BH_mass / gas_c_phys2;
const double omega_star = gas_vorticity * Bondi_radius / gas_c_phys;
const double f_omega_star = 1.0 / (1.0 + pow(omega_star, 0.9));
const double mdot_omega = 4. * M_PI * gas_rho_phys * G * G * BH_mass *
BH_mass * denominator_inv * denominator_inv *
denominator_inv * 0.34 * f_omega_star;
const double accr_rate_inv = 1. / accr_rate;
const double mdot_omega_inv = 1. / mdot_omega;
accr_rate = 1. / sqrt(accr_rate_inv * accr_rate_inv +
mdot_omega_inv * mdot_omega_inv);
} /* ends calculation of vorticity addition to Krumholz prescription */
} /* ends section without multi-phase accretion */
/* Compute the boost factor from Booth, Schaye (2009) */
if (props->with_boost_factor) {
const double gas_rho_phys = bp->rho_gas * cosmo->a3_inv;
const double n_H = gas_rho_phys * 0.75 / proton_mass;
const double boost_ratio = n_H / props->boost_n_h_star;
double boost_factor;
if (boost_ratio > 1.0) {
boost_factor = props->boost_alpha * pow(boost_ratio, props->boost_beta);
} else {
boost_factor = props->boost_alpha;
}
accr_rate *= boost_factor;
bp->accretion_boost_factor = boost_factor;
} else {
bp->accretion_boost_factor = 1.;
}
/* Compute the reduction factor from Rosas-Guevara et al. (2015) */
if (props->with_angmom_limiter) {
const double Bondi_radius = G * BH_mass / gas_c_phys2;
const double Bondi_time = Bondi_radius / gas_c_phys;
const double r_times_v_tang = Bondi_radius * tangential_velocity;
const double r_times_v_tang_3 =
r_times_v_tang * r_times_v_tang * r_times_v_tang;
const double viscous_time =
2. * M_PI * r_times_v_tang_3 /
(1e-6 * props->alpha_visc * G * G * BH_mass * BH_mass);
const double f_visc = min(Bondi_time / viscous_time, 1.);
bp->f_visc = f_visc;
/* Limit the accretion rate by the Bondi-to-viscous time ratio */
accr_rate *= f_visc;
} else {
bp->f_visc = 1.0;
}
/* Compute the Eddington rate (internal units) */
const double epsilon_r = props->epsilon_r;
const double Eddington_rate =
4. * M_PI * G * BH_mass * proton_mass / (epsilon_r * c * sigma_Thomson);
/* Should we record this time as the most recent high accretion rate? */
if (accr_rate > props->f_Edd_recording * Eddington_rate) {
if (with_cosmology) {
bp->last_high_Eddington_fraction_scale_factor = cosmo->a;
} else {
bp->last_high_Eddington_fraction_time = time;
}
}
/* Limit the accretion rate to a fraction of the Eddington rate */
bp->eddington_fraction = accr_rate / Eddington_rate;
accr_rate = min(accr_rate, props->f_Edd * Eddington_rate);
bp->accretion_rate = accr_rate;
/* Factor in the radiative efficiency */
const double mass_rate = (1. - epsilon_r) * accr_rate;
const double luminosity = epsilon_r * accr_rate * c * c;
/* Integrate forward in time */
bp->subgrid_mass += mass_rate * dt;
bp->total_accreted_mass += mass_rate * dt;
bp->energy_reservoir += luminosity * props->epsilon_f * dt;
if (props->use_nibbling) {
/* There is some loss in this due to radiative losses, so we must decrease
* the particle mass in proprtion to its current accretion rate. We do not
* account for this in the swallowing approach, however. */
bp->mass -= epsilon_r * accr_rate * dt;
if (bp->mass < 0)
error(
"Black hole %lld has reached a negative mass (%f). This is "
"not a great situation, so I am stopping.",
bp->id, bp->mass);
}
/* Increase the subgrid angular momentum according to what we accreted
* (already in physical units, a factors from velocity and radius cancel) */
bp->accreted_angular_momentum[0] +=
bp->spec_angular_momentum_gas[0] * mass_rate * dt;
bp->accreted_angular_momentum[1] +=
bp->spec_angular_momentum_gas[1] * mass_rate * dt;
bp->accreted_angular_momentum[2] +=
bp->spec_angular_momentum_gas[2] * mass_rate * dt;
/* 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. */
/* Mean gas particle mass in the BH's kernel */
const double mean_ngb_mass = bp->ngb_mass / ((double)bp->num_ngbs);
/* Compute AGN heating temperature */
const double dT_AGN = black_hole_feedback_delta_T(bp, props);
/* Energy per unit mass corresponding to the temperature jump delta_T */
double delta_u = dT_AGN * props->temp_to_u_factor;
/* Number of energy injections at this time-step (will be computed below) */
int number_of_energy_injections;
/* Average total energy needed to heat the target number of Ngbs */
const double E_feedback_event =
delta_u * mean_ngb_mass * props->num_ngbs_to_heat;
/* 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) {
/* Probability of heating. Relevant only for the stochastic model. */
double prob = bp->energy_reservoir / (delta_u * bp->ngb_mass);
/* Compute the number of energy injections based on probability if and
* only if we are using the stochastic (i.e. not deterministic) model
* and the probability prob < 1. */
if (prob < 1. && !props->AGN_deterministic) {
/* Initialise counter of energy injections */
number_of_energy_injections = 0;
/* 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_index(
bp->id, i, ti_begin, random_number_BH_feedback);
/* Increase the counter if we are lucky */
if (rand < prob) number_of_energy_injections++;
}
}
/* Determenistic model or prob >= 1. */
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));
}
/* Maximum number of energy injections allowed */
const int N_energy_injections_allowed =
min(colibre_blackhole_number_of_rays, bp->num_ngbs);
/* 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 {
/* Flag that we don't want to heat anyone */
bp->to_distribute.AGN_delta_u = 0.f;
bp->to_distribute.AGN_number_of_energy_injections = 0;
}
}
/**
* @brief Computes the (maximal) repositioning speed for a black hole.
*
* Calculated as upsilon * (m_BH / m_ref) ^ beta_m * (n_H_BH / n_ref) ^ beta_n
* where m_BH = BH subgrid mass, n_H_BH = physical gas density around BH
* and upsilon, m_ref, beta_m, n_ref, and beta_n are parameters.
*
* @param bp The #bpart.
* @param props The properties of the black hole model.
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static double
black_holes_get_repositioning_speed(const struct bpart* restrict bp,
const struct black_holes_props* props,
const struct cosmology* cosmo) {
const double n_gas_phys = bp->rho_gas * cosmo->a3_inv * props->rho_to_n_cgs;
const double v_repos =
props->reposition_coefficient_upsilon *
pow(bp->subgrid_mass / props->reposition_reference_mass,
props->reposition_exponent_mass) *
pow(n_gas_phys / props->reposition_reference_n_gas,
props->reposition_exponent_n_gas);
/* Make sure the repositioning is not back-firing... */
if (v_repos < 0)
error(
"BH %lld wants to reposition at negative speed (%g U_V). Do you "
"think you are being funny? No-one is laughing.",
bp->id, v_repos);
return v_repos;
}
/**
* @brief Finish the calculation of the new BH position.
*
* Here, we check that the BH should indeed be moved in the next drift.
*
* @param bp The black hole particle.
* @param props The properties of the black hole scheme.
* @param constants The physical constants (in internal units).
* @param cosmo The cosmological model.
* @param dt The black hole particle's time step.
* @param ti_begin The time at the start of the step
*/
__attribute__((always_inline)) INLINE static void black_holes_end_reposition(
struct bpart* restrict bp, const struct black_holes_props* props,
const struct phys_const* constants, const struct cosmology* cosmo,
const double dt, const integertime_t ti_begin) {
/* First check: did we find any eligible neighbour particle to jump to? */
if (bp->reposition.min_potential != FLT_MAX) {
/* Record that we have a (possible) repositioning situation */
bp->number_of_reposition_attempts++;
/* Is the potential lower (i.e. the BH is at the bottom already)
* OR is the BH massive enough that we don't reposition? */
const float potential = gravity_get_comoving_potential(bp->gpart);
if (potential < bp->reposition.min_potential ||
bp->subgrid_mass > props->max_reposition_mass) {
/* No need to reposition */
bp->reposition.min_potential = FLT_MAX;
bp->reposition.delta_x[0] = -FLT_MAX;
bp->reposition.delta_x[1] = -FLT_MAX;
bp->reposition.delta_x[2] = -FLT_MAX;
} else if (props->set_reposition_speed) {
/* If we are re-positioning, move the BH a fraction of delta_x, so
* that we have a well-defined re-positioning velocity (repos_vel
* is verified to be non-negative). */
double repos_vel = black_holes_get_repositioning_speed(bp, props, cosmo);
/* Convert target reposition velocity to a fractional reposition
* along reposition.delta_x */
const double dx = bp->reposition.delta_x[0];
const double dy = bp->reposition.delta_x[1];
const double dz = bp->reposition.delta_x[2];
const double d = sqrt(dx * dx + dy * dy + dz * dz);
/* Exclude the pathological case of repositioning by zero distance */
if (d > 0) {
double repos_frac = repos_vel * dt / d;
/* We should never get negative repositioning fractions... */
if (repos_frac < 0)
error("Wanting to reposition by negative fraction (%g)?", repos_frac);
/* ... but fractions > 1 can occur if the target velocity is high.
* We do not want this, because it could lead to overshooting the
* actual potential minimum. */
if (repos_frac > 1) {
repos_frac = 1.;
repos_vel = repos_frac * d / dt;
}
bp->last_repos_vel = (float)repos_vel;
bp->reposition.delta_x[0] *= repos_frac;
bp->reposition.delta_x[1] *= repos_frac;
bp->reposition.delta_x[2] *= repos_frac;
}
/* ends section for fractional repositioning */
} else {
/* We _should_ reposition, but not fractionally. Here, we will
* reposition exactly on top of another gas particle - which
* could cause issues, so we add on a small fractional offset
* of magnitude 0.001 h in the reposition delta. */
/* Generate three random numbers in the interval [-0.5, 0.5[; id,
* id**2, and id**3 are required to give unique random numbers (as
* random_unit_interval is completely reproducible). */
const float offset_dx =
random_unit_interval(bp->id, ti_begin, random_number_BH_reposition) -
0.5f;
const float offset_dy =
random_unit_interval(bp->id * bp->id, ti_begin,
random_number_BH_reposition) -
0.5f;
const float offset_dz =
random_unit_interval(bp->id * bp->id * bp->id, ti_begin,
random_number_BH_reposition) -
0.5f;
const float length_inv =
1.0f / sqrtf(offset_dx * offset_dx + offset_dy * offset_dy +
offset_dz * offset_dz);
const float norm = 0.001f * bp->h * length_inv;
bp->reposition.delta_x[0] += offset_dx * norm;
bp->reposition.delta_x[1] += offset_dy * norm;
bp->reposition.delta_x[2] += offset_dz * norm;
}
} /* ends section if we found eligible repositioning target(s) */
}
/**
* @brief Reset acceleration fields of a particle
*
* This is the equivalent of hydro_reset_acceleration.
* We do not compute the acceleration on black hole, therefore no need to use
* it.
*
* @param bp The particle to act upon
*/
__attribute__((always_inline)) INLINE static void black_holes_reset_feedback(
struct bpart* restrict bp) {
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)
bp->ids_ngbs_force[i] = -1;
bp->num_ngb_force = 0;
#endif
}
/**
* @brief Store the gravitational potential of a black hole by copying it from
* its #gpart friend.
*
* @param bp The black hole particle.
* @param gp The black hole's #gpart.
*/
__attribute__((always_inline)) INLINE static void
black_holes_store_potential_in_bpart(struct bpart* bp, const struct gpart* gp) {
#ifdef SWIFT_DEBUG_CHECKS
if (bp->gpart != gp) error("Copying potential to the wrong black hole!");
#endif
bp->reposition.potential = gp->potential;
}
/**
* @brief Store the gravitational potential of a particle by copying it from
* its #gpart friend.
*
* @param p_data The black hole data of a gas particle.
* @param gp The black hole's #gpart.
*/
__attribute__((always_inline)) INLINE static void
black_holes_store_potential_in_part(struct black_holes_part_data* p_data,
const struct gpart* gp) {
p_data->potential = gp->potential;
}
/**
* @brief Initialise a BH particle that has just been seeded.
*
* @param bp The #bpart to initialise.
* @param props The properties of the black hole scheme.
* @param constants The physical constants in internal units.
* @param cosmo The current cosmological model.
* @param p The #part that became a black hole.
* @param xp The #xpart that became a black hole.
*/
INLINE static void black_holes_create_from_gas(
struct bpart* bp, const struct black_holes_props* props,
const struct phys_const* constants, const struct cosmology* cosmo,
const struct part* p, const struct xpart* xp,
const integertime_t ti_current) {
/* All the non-basic properties of the black hole have been zeroed
* in the FOF code. We update them here.
* (i.e. position, velocity, mass, time-step have been set) */
/* Birth time */
bp->formation_scale_factor = cosmo->a;
/* Initial seed mass */
bp->subgrid_mass = props->subgrid_seed_mass;
/* We haven't accreted anything yet */
bp->total_accreted_mass = 0.f;
bp->cumulative_number_seeds = 1;
bp->number_of_mergers = 0;
bp->number_of_gas_swallows = 0;
bp->number_of_direct_gas_swallows = 0;
bp->cumulative_mass_lost_to_gw = 0.f;
bp->number_of_time_steps = 0;
/* We haven't repositioned yet, nor attempted it */
bp->number_of_repositions = 0;
bp->number_of_reposition_attempts = 0;
/* Copy over the splitting struct */
bp->split_data = xp->split_data;
/* Initial metal masses */
const float gas_mass = hydro_get_mass(p);
struct chemistry_bpart_data* bp_chem = &bp->chemistry_data;
const struct chemistry_part_data* p_chem = &p->chemistry_data;
chemistry_bpart_from_part(bp_chem, p_chem, gas_mass);
/* Initial dust masses */
struct dust_bpart_data* bp_dust = &bp->dust_data;
const struct dust_part_data* p_dust = &p->dust_data;
dust_bpart_from_part(bp_dust, p_dust, gas_mass);
/* No swallowed angular momentum */
bp->swallowed_angular_momentum[0] = 0.f;
bp->swallowed_angular_momentum[1] = 0.f;
bp->swallowed_angular_momentum[2] = 0.f;
/* Last time this BH had a high Eddington fraction */
bp->last_high_Eddington_fraction_scale_factor = -1.f;
/* Last time of mergers */
bp->last_minor_merger_time = -1.;
bp->last_major_merger_time = -1.;
/* First initialisation */
black_holes_init_bpart(bp);
black_holes_mark_bpart_as_not_swallowed(&bp->merger_data);
}
/**
* @brief Store the halo mass in the fof algorithm for the black
* hole particle.
*
* Nothing to do here.
*
* @param p_data The black hole particle data.
* @param halo_mass The halo mass to update.
*/
__attribute__((always_inline)) INLINE static void black_holes_update_halo_mass(
struct bpart* bp, float halo_mass) {}
#endif /* SWIFT_COLIBRE_BLACK_HOLES_H */