From 00b54df546fdaba2ab9e9f83a1a605e85533355f Mon Sep 17 00:00:00 2001
From: Yannick Bahe <bahe@strw.leidenuniv.nl>
Date: Fri, 4 Dec 2020 20:53:50 +0100
Subject: [PATCH] Added option to calculate BH sound speed with a limit to gas
particles' internal energy if near the EoS
---
src/black_holes/EAGLE/black_holes_iact.h | 60 ++++++++++++-------
.../EAGLE/black_holes_properties.h | 35 ++++++++---
src/runner_doiact_functions_black_holes.h | 8 +--
3 files changed, 72 insertions(+), 31 deletions(-)
diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h
index 1cf60a80d0..d248dd9156 100644
--- a/src/black_holes/EAGLE/black_holes_iact.h
+++ b/src/black_holes/EAGLE/black_holes_iact.h
@@ -28,6 +28,8 @@
#include "space.h"
#include "timestep_sync_part.h"
#include "tracers.h"
+#include "equation_of_state.h"
+#include "entropy_floor.h"
/**
* @brief Density interaction between two particles (non-symmetric).
@@ -52,7 +54,9 @@ runner_iact_nonsym_bh_gas_density(
struct bpart *bi, const struct part *pj, const struct xpart *xpj,
const int with_cosmology, const struct cosmology *cosmo,
const struct gravity_props *grav_props,
- const struct black_holes_props *bh_props, const integertime_t ti_current,
+ const struct black_holes_props *bh_props,
+ const struct entropy_floor_properties *floor_props,
+ const integertime_t ti_current,
const double time) {
float wi, wi_dx;
@@ -83,7 +87,25 @@ runner_iact_nonsym_bh_gas_density(
bi->ngb_mass += mj;
/* Contribution to the smoothed sound speed */
- const float cj = hydro_get_comoving_soundspeed(pj);
+ float cj = hydro_get_comoving_soundspeed(pj);
+ if (bh_props->with_fixed_T_near_EoS) {
+
+ /* Check whether we are close to the entropy floor. If we are, we
+ * re-calculate the sound speed using the fixed internal energy */
+ const float u_EoS = entropy_floor_temperature(pj, cosmo, floor_props)
+ * bh_props->temp_to_u_factor;
+ if (pj->u < u_EoS * bh_props->fixed_T_above_EoS_factor &&
+ pj->u > bh_props->fixed_u_for_soundspeed) {
+ /*const float cj_old = cj; */
+ cj = gas_soundspeed_from_internal_energy(
+ pj->rho, bh_props->fixed_u_for_soundspeed);
+ /*message("Changing cj=%g to %g for p ID=%lld. u=%g (T=%g); "
+ "rho=%g (nH=%g); u_EoS=%g (T_EoS=%g).",
+ cj_old, cj, pj->id, pj->u, pj->u / bh_props->temp_to_u_factor,
+ pj->rho, pj->rho * cosmo->a3_inv * bh_props->rho_to_n_cgs,
+ u_EoS, u_EoS / bh_props->temp_to_u_factor); */
+ }
+ }
bi->sound_speed_gas += mj * wi * cj;
/* Neighbour internal energy */
@@ -255,15 +277,14 @@ runner_iact_nonsym_bh_gas_density(
* @param time Current physical time in the simulation.
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
- const float hi, const float hj,
- struct bpart *bi, struct part *pj,
- struct xpart *xpj, const int with_cosmology,
- const struct cosmology *cosmo,
- const struct gravity_props *grav_props,
- const struct black_holes_props *bh_props,
- const integertime_t ti_current,
- const double time) {
+runner_iact_nonsym_bh_gas_swallow(
+ const float r2, const float *dx, const float hi, const float hj,
+ struct bpart *bi, struct part *pj, struct xpart *xpj,
+ const int with_cosmology, const struct cosmology *cosmo,
+ const struct gravity_props *grav_props,
+ const struct black_holes_props *bh_props,
+ const struct entropy_floor_properties *floor_props,
+ const integertime_t ti_current, const double time) {
float wi;
@@ -629,15 +650,14 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
* @param time current physical time in the simulation
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
- const float hi, const float hj,
- const struct bpart *bi, struct part *pj,
- struct xpart *xpj, const int with_cosmology,
- const struct cosmology *cosmo,
- const struct gravity_props *grav_props,
- const struct black_holes_props *bh_props,
- const integertime_t ti_current,
- const double time) {
+runner_iact_nonsym_bh_gas_feedback(
+ const float r2, const float *dx, const float hi, const float hj,
+ const struct bpart *bi, struct part *pj, struct xpart *xpj,
+ const int with_cosmology, const struct cosmology *cosmo,
+ const struct gravity_props *grav_props,
+ const struct black_holes_props *bh_props,
+ const struct entropy_floor_properties *floor_props,
+ const integertime_t ti_current, const double time) {
/* Number of energy injections per BH per time-step */
const int num_energy_injections_per_BH =
diff --git a/src/black_holes/EAGLE/black_holes_properties.h b/src/black_holes/EAGLE/black_holes_properties.h
index 5f86bedf34..551df0a96c 100644
--- a/src/black_holes/EAGLE/black_holes_properties.h
+++ b/src/black_holes/EAGLE/black_holes_properties.h
@@ -121,6 +121,16 @@ struct black_holes_props {
/*! Minimum gas particle mass in nibbling mode */
float min_gas_mass_for_nibbling;
+ /*! Switch to calculate the sound speed with a fixed T near the EoS */
+ int with_fixed_T_near_EoS;
+
+ /*! Factor above EoS below which fixed T applies for sound speed */
+ float fixed_T_above_EoS_factor;
+
+ /*! Fixed T (expressed as internal energy) for sound speed near EoS */
+ float fixed_u_for_soundspeed;
+
+
/* ---- Properties of the feedback model ------- */
/*! AGN feedback model: random, isotropic or minimum distance */
@@ -270,6 +280,13 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
const struct hydro_props *hydro_props,
const struct cosmology *cosmo) {
+ /* Calculate temperature to internal energy conversion factor (all internal
+ * units) */
+ const double k_B = phys_const->const_boltzmann_k;
+ const double m_p = phys_const->const_proton_mass;
+ const double mu = hydro_props->mu_ionised;
+ bp->temp_to_u_factor = k_B / (mu * hydro_gamma_minus_one * m_p);
+
/* Read in the basic neighbour search properties or default to the hydro
ones if the user did not provide any different values */
@@ -369,6 +386,17 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
bp->min_gas_mass_for_nibbling *= phys_const->const_solar_mass;
}
+ bp->with_fixed_T_near_EoS =
+ parser_get_param_int(params, "EAGLEAGN:with_fixed_T_near_EoS");
+ if (bp->with_fixed_T_near_EoS) {
+ bp->fixed_T_above_EoS_factor = pow(10.,
+ parser_get_param_float(params, "EAGLEAGN:fixed_T_above_EoS_dex"));
+ bp->fixed_u_for_soundspeed =
+ parser_get_param_float(params, "EAGLEAGN:fixed_T_near_EoS_K") /
+ units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE) *
+ bp->temp_to_u_factor;
+ }
+
/* Feedback parameters ---------------------------------- */
char temp[40];
@@ -530,13 +558,6 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
/* Common conversion factors ----------------------------- */
- /* Calculate temperature to internal energy conversion factor (all internal
- * units) */
- const double k_B = phys_const->const_boltzmann_k;
- const double m_p = phys_const->const_proton_mass;
- const double mu = hydro_props->mu_ionised;
- bp->temp_to_u_factor = k_B / (mu * hydro_gamma_minus_one * m_p);
-
/* Calculate conversion factor from rho to n_H.
* Note this assumes primoridal abundance */
const double X_H = hydro_props->hydrogen_mass_fraction;
diff --git a/src/runner_doiact_functions_black_holes.h b/src/runner_doiact_functions_black_holes.h
index 1cbd4d66d5..cd439300b8 100644
--- a/src/runner_doiact_functions_black_holes.h
+++ b/src/runner_doiact_functions_black_holes.h
@@ -103,7 +103,7 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
if (r2 < hig2) {
IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, with_cosmology, cosmo,
e->gravity_properties, e->black_holes_properties,
- ti_current, e->time);
+ e->entropy_floor, ti_current, e->time);
}
} /* loop over the parts in ci. */
} /* loop over the bparts in ci. */
@@ -257,7 +257,7 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
if (r2 < hig2) {
IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, with_cosmology, cosmo,
e->gravity_properties, e->black_holes_properties,
- ti_current, e->time);
+ e->entropy_floor, ti_current, e->time);
}
} /* loop over the parts in cj. */
} /* loop over the bparts in ci. */
@@ -425,7 +425,7 @@ void DOPAIR1_SUBSET_BH_NAIVE(struct runner *r, struct cell *restrict ci,
if (r2 < hig2) {
IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, with_cosmology, cosmo,
e->gravity_properties, e->black_holes_properties,
- ti_current, e->time);
+ e->entropy_floor, ti_current, e->time);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -503,7 +503,7 @@ void DOSELF1_SUBSET_BH(struct runner *r, struct cell *restrict ci,
if (r2 < hig2) {
IACT_BH_GAS(r2, dx, hi, pj->h, bi, pj, xpj, with_cosmology, cosmo,
e->gravity_properties, e->black_holes_properties,
- ti_current, e->time);
+ e->entropy_floor, ti_current, e->time);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
--
GitLab