Skip to content
Snippets Groups Projects
Commit 6e406b2a authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Fixes to the black holes time integration.

parent bc5a6188
No related branches found
No related tags found
No related merge requests found
......@@ -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 */
......@@ -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);
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment