diff --git a/src/black_holes/EAGLE/black_holes.h b/src/black_holes/EAGLE/black_holes.h index 541e0e94c3445cda9fb6ce52d67b14d5080fc17f..e493e89f0e20d82e4304f7c19003ad261a185b93 100644 --- a/src/black_holes/EAGLE/black_holes.h +++ b/src/black_holes/EAGLE/black_holes.h @@ -102,6 +102,7 @@ __attribute__((always_inline)) INLINE static void black_holes_init_bpart( 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; } /** @@ -584,6 +585,36 @@ __attribute__((always_inline)) INLINE static void black_holes_reset_feedback( #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. * diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h index f7aa064c3692f9a1b98a05e882bd66422bc9b3ac..5677ef23dff8cf49b503dc0f00dd91f371a0dbd1 100644 --- a/src/black_holes/EAGLE/black_holes_iact.h +++ b/src/black_holes/EAGLE/black_holes_iact.h @@ -172,7 +172,7 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx, /* Check the velocity criterion */ if (v2_pec < 0.25f * bi->sound_speed_gas * bi->sound_speed_gas) { - const float potential = gravity_get_comoving_potential(pj->gpart); + const float potential = pj->black_holes_data.potential; /* Is the potential lower? */ if (potential < bi->reposition.min_potential) { @@ -269,7 +269,7 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx, /* Check the velocity criterion */ if (v2_pec < 0.25f * bi->sound_speed_gas * bi->sound_speed_gas) { - const float potential = gravity_get_comoving_potential(bj->gpart); + const float potential = bj->reposition.potential; /* Is the potential lower? */ if (potential < bi->reposition.min_potential) { diff --git a/src/black_holes/EAGLE/black_holes_part.h b/src/black_holes/EAGLE/black_holes_part.h index 114baf3163dfe2b261017599c903c45736c02bd9..46dd43718d506537e7515f917ab637b17b334e8c 100644 --- a/src/black_holes/EAGLE/black_holes_part.h +++ b/src/black_holes/EAGLE/black_holes_part.h @@ -157,6 +157,9 @@ struct bpart { struct { + /*! Gravitational potential copied from the #gpart. */ + float potential; + /*! Value of the minimum potential across all neighbours. */ float min_potential; diff --git a/src/black_holes/EAGLE/black_holes_struct.h b/src/black_holes/EAGLE/black_holes_struct.h index 177a54b921bec24f357756e4f777b0ecf1781c30..c987d6608e9893f9502a7aa12cfe4e04cce45428 100644 --- a/src/black_holes/EAGLE/black_holes_struct.h +++ b/src/black_holes/EAGLE/black_holes_struct.h @@ -19,6 +19,10 @@ #ifndef SWIFT_BLACK_HOLES_STRUCT_EAGLE_H #define SWIFT_BLACK_HOLES_STRUCT_EAGLE_H +/* Standard headers */ +#include <float.h> + +/* Local includes */ #include "inline.h" /** @@ -28,6 +32,9 @@ struct black_holes_part_data { /*! ID of the black-hole that will swallow this #part. */ long long swallow_id; + + /*! Gravitational potential of the particle (for repositioning) */ + float potential; }; /** @@ -54,6 +61,11 @@ black_holes_mark_part_as_not_swallowed(struct black_holes_part_data* p_data) { p_data->swallow_id = -1; } +__attribute__((always_inline)) INLINE static void black_holes_init_potential( + struct black_holes_part_data* p_data) { + p_data->potential = FLT_MAX; +} + /** * @brief Update a given #part's BH data field to mark the particle has * having been been swallowed. diff --git a/src/cell.c b/src/cell.c index 96e4e03b1332023a86dc2c966fb25b50478bc1e5..4981ad1baaaf6db0c2a0afb3d29492706d93ecfb 100644 --- a/src/cell.c +++ b/src/cell.c @@ -4620,6 +4620,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) { /* Get ready for a density calculation */ if (part_is_active(p, e)) { hydro_init_part(p, &e->s->hs); + black_holes_init_potential(&p->black_holes_data); chemistry_init_part(p, e->chemistry); pressure_floor_init_part(p, xp); star_formation_init_part(p, e->star_formation); diff --git a/src/runner_others.c b/src/runner_others.c index 60607a1ef07d828d4d70cdc15bc0d722d04b9c9e..60e86bc0390bbdfeb98e01143f8f43e05cc822b8 100644 --- a/src/runner_others.c +++ b/src/runner_others.c @@ -455,6 +455,8 @@ void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) { void runner_do_end_grav_force(struct runner *r, struct cell *c, int timer) { const struct engine *e = r->e; + const int with_self_gravity = (e->policy & engine_policy_self_gravity); + const int with_black_holes = (e->policy & engine_policy_black_holes); TIMER_TIC; @@ -473,7 +475,7 @@ void runner_do_end_grav_force(struct runner *r, struct cell *c, int timer) { /* Potential normalisation in the case of periodic gravity */ float potential_normalisation = 0.; - if (periodic && (e->policy & engine_policy_self_gravity)) { + if (periodic && with_self_gravity) { const double volume = s->dim[0] * s->dim[1] * s->dim[2]; const double r_s = e->mesh->r_s; potential_normalisation = 4. * M_PI * e->total_mass * r_s * r_s / volume; @@ -557,6 +559,17 @@ void runner_do_end_grav_force(struct runner *r, struct cell *c, int timer) { } } #endif + + /* Deal with black holes' need of potentials */ + if (with_black_holes && gp->type == swift_type_black_hole) { + const size_t offset = -gp->id_or_neg_offset; + black_holes_store_potential_in_bpart(&s->bparts[offset], gp); + } + if (with_black_holes && gp->type == swift_type_gas) { + const size_t offset = -gp->id_or_neg_offset; + black_holes_store_potential_in_part( + &s->parts[offset].black_holes_data, gp); + } } } } diff --git a/src/space.c b/src/space.c index 575016fcdb11c85f1099aeeedbd4c13b6d408019..488c7d186727becd279a6a1c819e760f6a25b7b6 100644 --- a/src/space.c +++ b/src/space.c @@ -4554,6 +4554,7 @@ void space_init_parts_mapper(void *restrict map_data, int count, for (int k = 0; k < count; k++) { hydro_init_part(&parts[k], hs); + black_holes_init_potential(&parts[k].black_holes_data); chemistry_init_part(&parts[k], e->chemistry); pressure_floor_init_part(&parts[k], &xparts[k]); star_formation_init_part(&parts[k], e->star_formation);