diff --git a/src/kick.h b/src/kick.h index 404f7f31170dcec4728b8fbc52717dbf4a35bb0c..0e28e271aad04e9b16a4e988112aa95b36cb4f20 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 e1f5de709da804330953b47a647d0f0ce13de7bb..7fd22424a8c5cb5a5d87d60c5a234feb892906a3 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); }