Commit 311eef6a authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

In the black hole interactions functions, added a loop over BH pairs that is...

In the black hole interactions functions, added a loop over BH pairs that is only activated in the case of swallowing/merging.
parent a26f42b5
......@@ -32,11 +32,14 @@
* @param cosmo The cosmological model.
* @param ti_current Current integer time value (for random numbers).
*/
__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 struct xpart *restrict xpj, const struct cosmology *cosmo,
const integertime_t ti_current) {
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
const float hi, const float hj,
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;
......@@ -79,11 +82,37 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_density(
* @param cosmo The cosmological model.
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_swallow(
const float r2, const float *dx, const float hi, const float hj,
struct bpart *restrict bi, struct part *restrict pj,
struct xpart *restrict xpj, const struct cosmology *cosmo,
const integertime_t ti_current) {}
__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 *restrict bi,
struct part *restrict pj,
struct xpart *restrict xpj,
const struct cosmology *cosmo,
const integertime_t ti_current) {}
/**
* @brief Swallowing interaction between two BH particles (non-symmetric).
*
* Function used to flag the BH particles that will be swallowed
* by the black hole particle.
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param bi First particle (black hole).
* @param bj Second particle (black hole)
* @param cosmo The cosmological model.
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
const float hi, const float hj,
struct bpart *restrict bi,
struct bpart *restrict bj,
const struct cosmology *cosmo,
const integertime_t ti_current) {}
/**
* @brief Feedback interaction between two particles (non-symmetric).
......@@ -99,12 +128,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_swallow(
* @param ti_current Current integer time value (for random numbers).
*/
__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 struct xpart *restrict xpj,
const struct cosmology *cosmo,
const integertime_t ti_current) {
runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
const float hi, const float hj,
struct bpart *restrict bi,
struct part *restrict pj,
const struct xpart *restrict xpj,
const struct cosmology *cosmo,
const integertime_t ti_current) {
#ifdef DEBUG_INTERACTIONS_BH
/* Update ngb counters */
if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_BH)
......
......@@ -37,11 +37,14 @@
* @param cosmo The cosmological model.
* @param ti_current Current integer time value (for random numbers).
*/
__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 struct xpart *restrict xpj, const struct cosmology *cosmo,
const integertime_t ti_current) {
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
const float hi, const float hj,
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;
......@@ -116,11 +119,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_density(
* @param cosmo The cosmological model.
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_swallow(
const float r2, const float *dx, const float hi, const float hj,
struct bpart *restrict bi, struct part *restrict pj,
struct xpart *restrict xpj, const struct cosmology *cosmo,
const integertime_t ti_current) {
__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 *restrict bi,
struct part *restrict pj,
struct xpart *restrict xpj,
const struct cosmology *cosmo,
const integertime_t ti_current) {
float wi;
......@@ -169,6 +175,29 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_swallow(
}
}
/**
* @brief Swallowing interaction between two BH particles (non-symmetric).
*
* Function used to flag the BH particles that will be swallowed
* by the black hole particle.
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param bi First particle (black hole).
* @param bj Second particle (black hole)
* @param cosmo The cosmological model.
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void
runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
const float hi, const float hj,
struct bpart *restrict bi,
struct bpart *restrict bj,
const struct cosmology *cosmo,
const integertime_t ti_current) {}
/**
* @brief Feedback interaction between two particles (non-symmetric).
*
......@@ -183,12 +212,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_swallow(
* @param ti_current Current integer time value (for random numbers).
*/
__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,
struct xpart *restrict xpj,
const struct cosmology *cosmo,
const integertime_t ti_current) {
runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
const float hi, const float hj,
struct bpart *restrict bi,
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;
......
......@@ -79,8 +79,11 @@
#define _TIMER_DOSUB_PAIR_BH(f) PASTE(timer_dosub_pair_bh, f)
#define TIMER_DOSUB_PAIR_BH _TIMER_DOSUB_PAIR_BH(FUNCTION)
#define _IACT_BH(f) PASTE(runner_iact_nonsym_bh, f)
#define IACT_BH _IACT_BH(FUNCTION)
#define _IACT_BH_GAS(f) PASTE(runner_iact_nonsym_bh_gas, f)
#define IACT_BH_GAS _IACT_BH_GAS(FUNCTION)
#define _IACT_BH_BH(f) PASTE(runner_iact_nonsym_bh_bh, f)
#define IACT_BH_BH _IACT_BH_BH(FUNCTION)
/**
* @brief Calculate the number density of #part around the #bpart
......@@ -146,16 +149,73 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check that particles have been drifted to the current time */
if (bi->ti_drift != e->ti_current)
error("Particle bi not drifted to current time");
if (pj->ti_drift != e->ti_current)
error("Particle pj not drifted to current time");
#endif
if (r2 < hig2) {
IACT_BH(r2, dx, hi, hj, bi, pj, xpj, cosmo, ti_current);
IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, ti_current);
}
} /* loop over the parts in ci. */
} /* loop over the bparts in ci. */
/* When doing BH swallowing, we need a quick loop also over the BH
* neighbours */
#if (FUNCTION_TASK_LOOP == TASK_LOOP_SWALLOW)
/* Loop over the bparts in ci. */
for (int bid = 0; bid < bcount; bid++) {
/* Get a hold of the ith bpart in ci. */
struct bpart *restrict bi = &bparts[bid];
/* Skip inactive particles */
if (!bpart_is_active(bi, e)) continue;
const float hi = bi->h;
const float hig2 = hi * hi * kernel_gamma2;
const float bix[3] = {(float)(bi->x[0] - c->loc[0]),
(float)(bi->x[1] - c->loc[1]),
(float)(bi->x[2] - c->loc[2])};
/* Loop over the parts in cj. */
for (int bjd = 0; bjd < bcount; bjd++) {
/* Skip self interaction */
if (bid == bjd) continue;
/* Get a pointer to the jth particle. */
struct bpart *restrict bj = &bparts[bjd];
const float hj = bj->h;
/* Early abort? */
if (bpart_is_inhibited(bj, e)) continue;
/* Compute the pairwise distance. */
const float bjx[3] = {(float)(bj->x[0] - c->loc[0]),
(float)(bj->x[1] - c->loc[1]),
(float)(bj->x[2] - c->loc[2])};
float dx[3] = {bix[0] - bjx[0], bix[1] - bjx[1], bix[2] - bjx[2]};
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
/* Check that particles have been drifted to the current time */
if (bi->ti_drift != e->ti_current)
error("Particle bi not drifted to current time");
if (bj->ti_drift != e->ti_current)
error("Particle bj not drifted to current time");
#endif
if (r2 < hig2) {
IACT_BH_BH(r2, dx, hi, hj, bi, bj, cosmo, ti_current);
}
} /* loop over the bparts in ci. */
} /* loop over the bparts in ci. */
#endif /* (FUNCTION_TASK_LOOP == TASK_LOOP_SWALLOW) */
TIMER_TOC(TIMER_DOSELF_BH);
}
......@@ -235,15 +295,72 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
#ifdef SWIFT_DEBUG_CHECKS
/* Check that particles have been drifted to the current time */
if (bi->ti_drift != e->ti_current)
error("Particle bi not drifted to current time");
if (pj->ti_drift != e->ti_current)
error("Particle pj not drifted to current time");
#endif
if (r2 < hig2) {
IACT_BH(r2, dx, hi, hj, bi, pj, xpj, cosmo, ti_current);
IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, ti_current);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
} /* loop over the bparts in ci. */
/* When doing BH swallowing, we need a quick loop also over the BH
* neighbours */
#if (FUNCTION_TASK_LOOP == TASK_LOOP_SWALLOW)
const int bcount_j = cj->black_holes.count;
struct bpart *restrict bparts_j = cj->black_holes.parts;
/* Loop over the bparts in ci. */
for (int bid = 0; bid < bcount_i; bid++) {
/* Get a hold of the ith bpart in ci. */
struct bpart *restrict bi = &bparts_i[bid];
/* Skip inactive particles */
if (!bpart_is_active(bi, e)) continue;
const float hi = bi->h;
const float hig2 = hi * hi * kernel_gamma2;
const float bix[3] = {(float)(bi->x[0] - (cj->loc[0] + shift[0])),
(float)(bi->x[1] - (cj->loc[1] + shift[1])),
(float)(bi->x[2] - (cj->loc[2] + shift[2]))};
/* Loop over the bparts in cj. */
for (int bjd = 0; bjd < bcount_j; bjd++) {
/* Get a pointer to the jth particle. */
struct bpart *restrict bj = &bparts_j[bjd];
const float hj = bj->h;
/* Skip inhibited particles. */
if (bpart_is_inhibited(bj, e)) continue;
/* Compute the pairwise distance. */
const float bjx[3] = {(float)(bj->x[0] - cj->loc[0]),
(float)(bj->x[1] - cj->loc[1]),
(float)(bj->x[2] - cj->loc[2])};
float dx[3] = {bix[0] - bjx[0], bix[1] - bjx[1], bix[2] - bjx[2]};
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
/* Check that particles have been drifted to the current time */
if (bi->ti_drift != e->ti_current)
error("Particle bi not drifted to current time");
if (bj->ti_drift != e->ti_current)
error("Particle bj not drifted to current time");
#endif
if (r2 < hig2) {
IACT_BH_BH(r2, dx, hi, hj, bi, bj, cosmo, ti_current);
}
} /* loop over the bparts in cj. */
} /* loop over the bparts in ci. */
#endif /* (FUNCTION_TASK_LOOP == TASK_LOOP_SWALLOW) */
}
void DOPAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
......@@ -350,7 +467,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, bi, pj, xpj, cosmo, ti_current);
IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, ti_current);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
......@@ -425,7 +542,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, bi, pj, xpj, cosmo, ti_current);
IACT_BH_GAS(r2, dx, hi, pj->h, bi, pj, xpj, cosmo, ti_current);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment