From 6e406b2a3c0ba481c27eaf6c0ffac705a7bda6db Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Wed, 18 Sep 2019 18:58:06 +0200
Subject: [PATCH] Fixes to the black holes time integration.
---
src/kick.h | 37 +++++++++++++++
src/runner_time_integration.c | 84 +++++++++++++++++++++++++++++++++++
2 files changed, 121 insertions(+)
diff --git a/src/kick.h b/src/kick.h
index 404f7f3117..0e28e271aa 100644
--- a/src/kick.h
+++ b/src/kick.h
@@ -155,4 +155,41 @@ __attribute__((always_inline)) INLINE static void kick_spart(
stars_kick_extra(sp, dt_kick_grav);
}
+/**
+ * @brief Perform the 'kick' operation on a #bpart
+ *
+ * @param bp The #bpart to kick.
+ * @param dt_kick_grav The kick time-step for gravity accelerations.
+ * @param ti_start The starting (integer) time of the kick (for debugging
+ * checks).
+ * @param ti_end The ending (integer) time of the kick (for debugging checks).
+ */
+__attribute__((always_inline)) INLINE static void kick_bpart(
+ struct bpart *restrict bp, double dt_kick_grav, integertime_t ti_start,
+ integertime_t ti_end) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (bp->ti_kick != ti_start)
+ error(
+ "s-particle has not been kicked to the current time bp->ti_kick=%lld, "
+ "ti_start=%lld, ti_end=%lld id=%lld",
+ bp->ti_kick, ti_start, ti_end, bp->id);
+
+ bp->ti_kick = ti_end;
+#endif
+
+ /* Kick particles in momentum space */
+ bp->v[0] += bp->gpart->a_grav[0] * dt_kick_grav;
+ bp->v[1] += bp->gpart->a_grav[1] * dt_kick_grav;
+ bp->v[2] += bp->gpart->a_grav[2] * dt_kick_grav;
+
+ /* Give the gpart friend the same velocity */
+ bp->gpart->v_full[0] = bp->v[0];
+ bp->gpart->v_full[1] = bp->v[1];
+ bp->gpart->v_full[2] = bp->v[2];
+
+ /* Kick extra variables */
+ black_holes_kick_extra(bp, dt_kick_grav);
+}
+
#endif /* SWIFT_KICK_H */
diff --git a/src/runner_time_integration.c b/src/runner_time_integration.c
index e1f5de709d..7fd22424a8 100644
--- a/src/runner_time_integration.c
+++ b/src/runner_time_integration.c
@@ -88,9 +88,11 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
struct xpart *restrict xparts = c->hydro.xparts;
struct gpart *restrict gparts = c->grav.parts;
struct spart *restrict sparts = c->stars.parts;
+ struct bpart *restrict bparts = c->black_holes.parts;
const int count = c->hydro.count;
const int gcount = c->grav.count;
const int scount = c->stars.count;
+ const int bcount = c->black_holes.count;
const integertime_t ti_current = e->ti_current;
const double time_base = e->time_base;
@@ -249,6 +251,44 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
kick_spart(sp, dt_kick_grav, ti_begin, ti_begin + ti_step / 2);
}
}
+
+ /* Loop over the black hole particles in this cell. */
+ for (int k = 0; k < bcount; k++) {
+
+ /* Get a handle on the s-part. */
+ struct bpart *restrict bp = &bparts[k];
+
+ /* If particle needs to be kicked */
+ if (bpart_is_starting(bp, e)) {
+
+ const integertime_t ti_step = get_integer_timestep(bp->time_bin);
+ const integertime_t ti_begin =
+ get_integer_time_begin(ti_current + 1, bp->time_bin);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ const integertime_t ti_end =
+ get_integer_time_end(ti_current + 1, bp->time_bin);
+
+ if (ti_begin != ti_current)
+ error(
+ "Particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, "
+ "ti_step=%lld time_bin=%d ti_current=%lld",
+ ti_end, ti_begin, ti_step, bp->time_bin, ti_current);
+#endif
+
+ /* Time interval for this half-kick */
+ double dt_kick_grav;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_begin,
+ ti_begin + ti_step / 2);
+ } else {
+ dt_kick_grav = (ti_step / 2) * time_base;
+ }
+
+ /* do the kick */
+ kick_bpart(bp, dt_kick_grav, ti_begin, ti_begin + ti_step / 2);
+ }
+ }
}
if (timer) TIMER_TOC(timer_kick1);
@@ -273,10 +313,12 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
const int count = c->hydro.count;
const int gcount = c->grav.count;
const int scount = c->stars.count;
+ const int bcount = c->black_holes.count;
struct part *restrict parts = c->hydro.parts;
struct xpart *restrict xparts = c->hydro.xparts;
struct gpart *restrict gparts = c->grav.parts;
struct spart *restrict sparts = c->stars.parts;
+ struct bpart *restrict bparts = c->black_holes.parts;
const integertime_t ti_current = e->ti_current;
const double time_base = e->time_base;
@@ -453,6 +495,48 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
stars_reset_predicted_values(sp);
}
}
+
+ /* Loop over the particles in this cell. */
+ for (int k = 0; k < bcount; k++) {
+
+ /* Get a handle on the part. */
+ struct bpart *restrict bp = &bparts[k];
+
+ /* If particle needs to be kicked */
+ if (bpart_is_active(bp, e)) {
+
+ const integertime_t ti_step = get_integer_timestep(bp->time_bin);
+ const integertime_t ti_begin =
+ get_integer_time_begin(ti_current, bp->time_bin);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ti_begin + ti_step != ti_current)
+ error("Particle in wrong time-bin");
+#endif
+
+ /* Time interval for this half-kick */
+ double dt_kick_grav;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(
+ cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
+ } else {
+ dt_kick_grav = (ti_step / 2) * time_base;
+ }
+
+ /* Finish the time-step with a second half-kick */
+ kick_bpart(bp, dt_kick_grav, ti_begin + ti_step / 2,
+ ti_begin + ti_step);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that kick and the drift are synchronized */
+ if (bp->ti_drift != bp->ti_kick)
+ error("Error integrating b-part in time.");
+#endif
+
+ /* Prepare the values to be drifted */
+ black_holes_reset_predicted_values(bp);
+ }
+ }
}
if (timer) TIMER_TOC(timer_kick2);
}
--
GitLab