From d60c59bad12093bb4e65a467c581e2872b1d2f05 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Fri, 8 Nov 2019 14:02:36 +0100 Subject: [PATCH] In the black hole reposition code, calculate a delta position rather than an absolute position. This correctly handles the periodic boundary conditions. Also update the x_diff of the BH particles. --- src/black_holes/EAGLE/black_holes.h | 49 ++++++++++++++++-------- src/black_holes/EAGLE/black_holes_iact.h | 12 +++--- src/black_holes/EAGLE/black_holes_part.h | 4 +- 3 files changed, 40 insertions(+), 25 deletions(-) diff --git a/src/black_holes/EAGLE/black_holes.h b/src/black_holes/EAGLE/black_holes.h index f600f6eeb8..ae012f70aa 100644 --- a/src/black_holes/EAGLE/black_holes.h +++ b/src/black_holes/EAGLE/black_holes.h @@ -90,9 +90,9 @@ __attribute__((always_inline)) INLINE static void black_holes_init_bpart( bp->circular_velocity_gas[2] = 0.f; bp->ngb_mass = 0.f; bp->num_ngbs = 0; - bp->reposition.x[0] = -FLT_MAX; - bp->reposition.x[1] = -FLT_MAX; - bp->reposition.x[2] = -FLT_MAX; + 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; } @@ -112,24 +112,39 @@ __attribute__((always_inline)) INLINE static void black_holes_predict_extra( if (bp->reposition.min_potential != FLT_MAX) { #ifdef SWIFT_DEBUG_CHECKS - if (bp->reposition.x[0] == -FLT_MAX || bp->reposition.x[1] == -FLT_MAX || - bp->reposition.x[2] == -FLT_MAX) { + 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 - bp->x[0] = bp->reposition.x[0]; - bp->x[1] = bp->reposition.x[1]; - bp->x[2] = bp->reposition.x[2]; + /* 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]; - bp->gpart->x[0] = bp->reposition.x[0]; - bp->gpart->x[1] = bp->reposition.x[1]; - bp->gpart->x[2] = bp->reposition.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.x[0] = -FLT_MAX; - bp->reposition.x[1] = -FLT_MAX; - bp->reposition.x[2] = -FLT_MAX; + 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; } } @@ -487,9 +502,9 @@ __attribute__((always_inline)) INLINE static void black_holes_end_reposition( /* No need to reposition */ bp->reposition.min_potential = FLT_MAX; - bp->reposition.x[0] = -FLT_MAX; - bp->reposition.x[1] = -FLT_MAX; - bp->reposition.x[2] = -FLT_MAX; + bp->reposition.delta_x[0] = -FLT_MAX; + bp->reposition.delta_x[1] = -FLT_MAX; + bp->reposition.delta_x[2] = -FLT_MAX; } } diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h index 3e2238dc2a..eee2d648aa 100644 --- a/src/black_holes/EAGLE/black_holes_iact.h +++ b/src/black_holes/EAGLE/black_holes_iact.h @@ -170,9 +170,9 @@ runner_iact_nonsym_bh_gas_swallow( /* Store this as our new best */ bi->reposition.min_potential = potential; - bi->reposition.x[0] = pj->x[0]; - bi->reposition.x[1] = pj->x[1]; - bi->reposition.x[2] = pj->x[2]; + bi->reposition.delta_x[0] = -dx[0]; + bi->reposition.delta_x[1] = -dx[1]; + bi->reposition.delta_x[2] = -dx[2]; } } } @@ -266,9 +266,9 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx, /* Store this as our new best */ bi->reposition.min_potential = potential; - bi->reposition.x[0] = bj->x[0]; - bi->reposition.x[1] = bj->x[1]; - bi->reposition.x[2] = bj->x[2]; + bi->reposition.delta_x[0] = -dx[0]; + bi->reposition.delta_x[1] = -dx[1]; + bi->reposition.delta_x[2] = -dx[2]; } } } diff --git a/src/black_holes/EAGLE/black_holes_part.h b/src/black_holes/EAGLE/black_holes_part.h index 4d8a9c8e17..364621f242 100644 --- a/src/black_holes/EAGLE/black_holes_part.h +++ b/src/black_holes/EAGLE/black_holes_part.h @@ -127,8 +127,8 @@ struct bpart { /*! Value of the minimum potential across all neighbours. */ float min_potential; - /*! New position of the BH following the reposition procedure */ - double x[3]; + /*! Delta position to apply after the reposition procedure */ + double delta_x[3]; } reposition; -- GitLab