diff --git a/src/black_holes/Default/black_holes.h b/src/black_holes/Default/black_holes.h
index 38df7ad7b7126bcc4747ebf40ee623cec1d2f5d7..a65651d6b054efe980af8c6ace4970a3a8de73be 100644
--- a/src/black_holes/Default/black_holes.h
+++ b/src/black_holes/Default/black_holes.h
@@ -247,12 +247,13 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
* @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 time-step size (in physical internal units).
+ * @param dt The black hole particle's time step.
+ * @param ti_begin The time at the start of the temp
*/
__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 double dt, const integertime_t ti_begin) {}
/**
* @brief Reset acceleration fields of a particle
diff --git a/src/black_holes/Default/black_holes_iact.h b/src/black_holes/Default/black_holes_iact.h
index 012ba050edda07b70fdb4d1a97628ad51ab078e9..f58825ebdbffc6b492a1b1514da8e5dee74f8920 100644
--- a/src/black_holes/Default/black_holes_iact.h
+++ b/src/black_holes/Default/black_holes_iact.h
@@ -48,8 +48,7 @@ runner_iact_nonsym_bh_gas_density(
float wi, wi_dx;
/* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ const float r = sqrtf(r2);
/* Compute the kernel function */
const float hi_inv = 1.0f / hi;
diff --git a/src/black_holes/EAGLE/black_holes.h b/src/black_holes/EAGLE/black_holes.h
index e26f4d726b99025e7c104bc17cb3ef61892c417e..0f9f1f969e7b0b6853a3adf3948b5677daaff9ee 100644
--- a/src/black_holes/EAGLE/black_holes.h
+++ b/src/black_holes/EAGLE/black_holes.h
@@ -603,7 +603,7 @@ black_hole_energy_reservoir_threshold(struct bpart* bp,
}
/**
- * @brief Compute the accretion rate of the black hole and all the quantites
+ * @brief Compute the accretion rate of the black hole and all the quantities
* required for the feedback loop.
*
* @param bp The black hole particle.
@@ -615,6 +615,7 @@ black_hole_energy_reservoir_threshold(struct bpart* bp,
* @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 The time at which the step begun (ti_current).
*/
__attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
struct bpart* restrict bp, const struct black_holes_props* props,
@@ -1083,11 +1084,12 @@ black_holes_get_repositioning_speed(const struct bpart* restrict bp,
* @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 temp
*/
__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 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) {
@@ -1142,8 +1144,38 @@ __attribute__((always_inline)) INLINE static void black_holes_end_reposition(
bp->reposition.delta_x[1] *= repos_frac;
bp->reposition.delta_x[2] *= repos_frac;
}
- } /* ends section for fractional repositioning */
- } /* ends section if we found eligible repositioning target(s) */
+ } /* 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_y[0] += offset_dy * norm;
+ bp->reposition.delta_z[0] += offset_dz * norm;
+ }
+ } /* ends section if we found eligible repositioning target(s) */
}
/**
diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h
index 7007f3ca08c7b50674de47f17bea0013c065ca29..4e147b008a156f2315688027b9c8cfda15f455d0 100644
--- a/src/black_holes/EAGLE/black_holes_iact.h
+++ b/src/black_holes/EAGLE/black_holes_iact.h
@@ -60,11 +60,9 @@ runner_iact_nonsym_bh_gas_density(
float wi, wi_dx;
- /* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
-
- /* Compute the kernel function */
+ /* Compute the kernel function; note that r cannot be optimised
+ * to r2 / sqrtf(r2) because of 1 / 0 behaviour. */
+ const float r = sqrtf(r2);
const float hi_inv = 1.0f / hi;
const float ui = r * hi_inv;
kernel_deval(ui, &wi, &wi_dx);
@@ -283,11 +281,9 @@ runner_iact_nonsym_bh_gas_swallow(
float wi;
- /* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
-
- /* Compute the kernel function */
+ /* Compute the kernel function; note that r cannot be optimised
+ * to r2 / sqrtf(r2) because of 1 / 0 behaviour. */
+ const float r = sqrtf(r2);
const float hi_inv = 1.0f / hi;
const float hi_inv_dim = pow_dimension(hi_inv);
const float ui = r * hi_inv;
diff --git a/src/feedback/GEAR/feedback_iact.h b/src/feedback/GEAR/feedback_iact.h
index be07e5d3508f52523f8ffcbe00598432e152ff29..e906f3abd13a177c6d6bd76ade76485d80430aa7 100644
--- a/src/feedback/GEAR/feedback_iact.h
+++ b/src/feedback/GEAR/feedback_iact.h
@@ -50,9 +50,8 @@ runner_iact_nonsym_feedback_density(const float r2, const float *dx,
/* Get the gas mass. */
const float mj = hydro_get_mass(pj);
- /* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r */
+ const float r = sqrtf(r2);
/* Compute the kernel function */
const float hi_inv = 1.0f / hi;
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index cccc485ec319c5e592bc6895fb29f9b69e578f72..e6e7e8ba5870fa2fd2635f5a3cbba183a18b2231 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -68,8 +68,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
const float mj = pj->mass;
/* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Compute the kernel function for pi */
const float hi_inv = 1.f / hi;
@@ -164,8 +164,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
const float mj = pj->mass;
/* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Compute the kernel function */
const float hi_inv = 1.0f / hi;
@@ -463,8 +463,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float a2_Hubble = a * a * H;
/* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Get some values in local variables. */
const float mi = pi->mass;
@@ -592,8 +592,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float a2_Hubble = a * a * H;
/* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Get some values in local variables. */
const float mj = pj->mass;
diff --git a/src/hydro/Gizmo/hydro_gradients_gizmo.h b/src/hydro/Gizmo/hydro_gradients_gizmo.h
index c8914129bc58e32104ce400a59f1ac94b5a83b36..99ae8e6b7cdda29059fb826545d167846eec78c0 100644
--- a/src/hydro/Gizmo/hydro_gradients_gizmo.h
+++ b/src/hydro/Gizmo/hydro_gradients_gizmo.h
@@ -50,8 +50,9 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj) {
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
float wi, wj, wi_dx, wj_dx;
float Bi[3][3];
@@ -178,8 +179,9 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
struct part *restrict pi,
struct part *restrict pj) {
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
float Bi[3][3];
float Wi[5], Wj[5];
diff --git a/src/hydro/Gizmo/hydro_gradients_sph.h b/src/hydro/Gizmo/hydro_gradients_sph.h
index 9ece14c32d3c1ca7cd91a668c3142bac5f9b6010..cbe33514bc7b3f8046faf201aee59515f2d3b269 100644
--- a/src/hydro/Gizmo/hydro_gradients_sph.h
+++ b/src/hydro/Gizmo/hydro_gradients_sph.h
@@ -50,8 +50,9 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj) {
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
float wi, wi_dx;
const float hi_inv = 1.0f / hi;
@@ -142,8 +143,9 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
struct part *restrict pi,
struct part *restrict pj) {
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
float wi, wi_dx;
const float hi_inv = 1.0f / hi;
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index bc34bae489fb9eeec4d682c417b4996447424ee5..3760503e19d7d29223a8055191eaafd6af727b0a 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -213,8 +213,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj, int mode, float a, float H) {
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Initialize local variables */
float Bi[3][3];
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 768b0358ab3457813fe503f77d69bb7828381af4..20bc84c33ad21e7699910c0aafcedbbf89030a89 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -61,9 +61,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
error("Inhibited pj in interaction function!");
#endif
- /* Get r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Get the masses. */
const float mi = pi->mass;
@@ -145,9 +145,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* Get the masses. */
const float mj = pj->mass;
- /* Get r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
const float h_inv = 1.f / hi;
const float ui = r * h_inv;
@@ -207,9 +207,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
- /* Get r and r inverse. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Recover some data */
const float mi = pi->mass;
@@ -339,9 +339,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
- /* Get r and r inverse. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Recover some data */
const float mi = pi->mass;
diff --git a/src/hydro/Planetary/hydro_iact.h b/src/hydro/Planetary/hydro_iact.h
index b9a22c8382da6c32608cdf4194ceb015fc0bc7e6..fae1571acffdec1609ceac0b9b4b47ab522f8b4c 100644
--- a/src/hydro/Planetary/hydro_iact.h
+++ b/src/hydro/Planetary/hydro_iact.h
@@ -63,9 +63,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
error("Inhibited pj in interaction function!");
#endif
- /* Get r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Get the masses. */
const float mi = pi->mass;
@@ -147,9 +147,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* Get the masses. */
const float mj = pj->mass;
- /* Get r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
const float h_inv = 1.f / hi;
const float ui = r * h_inv;
@@ -209,9 +209,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
- /* Get r and r inverse. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
+ ;
/* Recover some data */
const float mi = pi->mass;
@@ -338,9 +339,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
- /* Get r and r inverse. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Recover some data */
const float mi = pi->mass;
diff --git a/src/hydro/PressureEntropy/hydro_iact.h b/src/hydro/PressureEntropy/hydro_iact.h
index ddf5297254929f6e0f45cbc568b8d8ca5949caec..6a426ad2c752d0b164f879ad5938e094d5725df1 100644
--- a/src/hydro/PressureEntropy/hydro_iact.h
+++ b/src/hydro/PressureEntropy/hydro_iact.h
@@ -57,8 +57,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
const float mj = pj->mass;
/* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Compute the kernel function for pi */
const float hi_inv = 1.f / hi;
@@ -145,8 +145,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
const float mj = pj->mass;
/* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Compute the kernel function */
const float h_inv = 1.0f / hi;
@@ -207,9 +207,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
- /* Get r and 1/r . */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ /* Get r and 1/r. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Get some values in local variables. */
const float mi = pi->mass;
@@ -321,8 +321,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float a2_Hubble = a * a * H;
/* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
+ ;
/* Get some values in local variables. */
// const float mi = pi->mass;
diff --git a/src/random.h b/src/random.h
index 1ecb2683f558df319bdbb216511aa4ab14a1d9b5..9c4e35f641a508df1c94549f9635017b8b620eb0 100644
--- a/src/random.h
+++ b/src/random.h
@@ -42,7 +42,6 @@
* generator.
* In case new numbers need to be added other possible
* numbers could be:
- * 59969537
* 65610001
* 126247697
* 193877777
@@ -64,6 +63,7 @@ enum random_number_type {
random_number_stellar_enrichment = 2936881973LL,
random_number_BH_feedback = 1640531371LL,
random_number_BH_swallow = 4947009007LL,
+ random_number_BH_reposition = 59969537LL,
};
#ifndef __APPLE__
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index c1096bf023a33469cd0192679b27e54900d4993a..4c69714b6cdd2ac7eadae71e496dff5d67565da0 100644
--- a/src/runner_ghost.c
+++ b/src/runner_ghost.c
@@ -897,7 +897,8 @@ void runner_do_black_holes_swallow_ghost(struct runner *r, struct cell *c,
/* Compute the final operations for repositioning of this BH */
black_holes_end_reposition(bp, e->black_holes_properties,
- e->physical_constants, e->cosmology, dt);
+ e->physical_constants, e->cosmology, dt,
+ e->ti_current);
/* Compute variables required for the feedback loop */
black_holes_prepare_feedback(bp, e->black_holes_properties,
diff --git a/src/stars/Basic/stars_iact.h b/src/stars/Basic/stars_iact.h
index 3efd75e9efbe3d77ed7dff0bedc8dc38c85b94b6..a573837663d1bdb82742c2c79309d1703421fdb3 100644
--- a/src/stars/Basic/stars_iact.h
+++ b/src/stars/Basic/stars_iact.h
@@ -42,8 +42,8 @@ runner_iact_nonsym_stars_density(const float r2, const float *dx,
float wi, wi_dx;
/* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Compute the kernel function */
const float hi_inv = 1.0f / hi;
diff --git a/src/stars/GEAR/stars_iact.h b/src/stars/GEAR/stars_iact.h
index 024ac7a9b9fff9ed707298d60ec2e562a5edf6bd..81cd48110dfc8f35f510a0dc39faca8ee1fa8cd2 100644
--- a/src/stars/GEAR/stars_iact.h
+++ b/src/stars/GEAR/stars_iact.h
@@ -41,8 +41,8 @@ runner_iact_nonsym_stars_density(const float r2, const float *dx,
float wi, wi_dx;
/* Get r and 1/r. */
- const float r_inv = 1.0f / sqrtf(r2);
- const float r = r2 * r_inv;
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
/* Compute the kernel function */
const float hi_inv = 1.0f / hi;