diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index e533055dd6fdfbe531eac8d00621783e464a716f..bf8362a573519ecfacb54ae6001aedb74e429db3 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -96,7 +96,6 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
*
* @param gp The particle to act upon
* @param dt The time-step for this kick
- * @param half_dt The half time-step for this kick
*/
__attribute__((always_inline)) INLINE static void gravity_kick_extra(
struct gpart* gp, float dt) {}
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 3b30c1b34c87b32cc28130e3ebb011705fdc2884..eba1a03663b007ca31ca16cac1b7f649cf3bb04b 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -354,6 +354,13 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
p->force.v_sig = 0.0f;
}
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param p The particle.
+ * @param xp The extended data of this particle.
+ */
__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
struct part *restrict p, const struct xpart *restrict xp) {
@@ -433,7 +440,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param p The particle to act upon
* @param xp The particle extended data to act upon
* @param dt The time-step for this kick
- * @param half_dt The half time-step for this kick
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt) {
diff --git a/src/kick.h b/src/kick.h
index 08dd0307ee829fd6079b91f5961aaec4de280407..d532b2a9abb0cc234cd7d55d6636e8dae1b059ab 100644
--- a/src/kick.h
+++ b/src/kick.h
@@ -27,42 +27,29 @@
#include "debug.h"
#include "timeline.h"
-#if 0
-
/**
* @brief Perform the 'kick' operation on a #gpart
*
* @param gp The #gpart to kick.
- * @param new_dti The (integer) time-step for this kick.
+ * @param ti_start The starting (integer) time of the kick
+ * @param ti_end The ending (integer) time of the kick
* @param ti_current The current (integer) time.
* @param timeBase The minimal allowed time-step size.
*/
__attribute__((always_inline)) INLINE static void kick_gpart(
- struct gpart *restrict gp, integertime_t new_dti, integertime_t ti_current,
- double timeBase) {
+ struct gpart *restrict gp, integertime_t ti_start, integertime_t ti_end,
+ integertime_t ti_current, double timeBase) {
- /* Compute the time step for this kick */
- const integertime_t old_ti_begin =
- get_integer_time_begin(ti_current, gp->time_bin);
- const integertime_t old_ti_end =
- get_integer_time_end(ti_current, gp->time_bin);
- const integertime_t ti_start = (old_ti_begin + old_ti_end) / 2;
- const integertime_t ti_end = old_ti_end + new_dti / 2;
+ /* Time interval for this half-kick */
const float dt = (ti_end - ti_start) * timeBase;
- const float half_dt = (ti_end - old_ti_end) * timeBase;
-
- /* Move particle forward in time */
- // gp->ti_begin = gp->ti_end;
- // gp->ti_end = gp->ti_begin + new_dti;
- gp->time_bin = get_time_bin(new_dti);
/* Kick particles in momentum space */
gp->v_full[0] += gp->a_grav[0] * dt;
gp->v_full[1] += gp->a_grav[1] * dt;
gp->v_full[2] += gp->a_grav[2] * dt;
- /* Extra kick work */
- gravity_kick_extra(gp, dt, half_dt);
+ /* Kick extra variables */
+ gravity_kick_extra(gp, dt);
}
/**
@@ -70,66 +57,11 @@ __attribute__((always_inline)) INLINE static void kick_gpart(
*
* @param p The #part to kick.
* @param xp The #xpart of the particle.
- * @param new_dti The (integer) time-step for this kick.
+ * @param ti_start The starting (integer) time of the kick
+ * @param ti_end The ending (integer) time of the kick
* @param ti_current The current (integer) time.
* @param timeBase The minimal allowed time-step size.
*/
-__attribute__((always_inline)) INLINE static void kick_part(
- struct part *restrict p, struct xpart *restrict xp, integertime_t new_dti,
- integertime_t ti_current, double timeBase) {
-
- /* Compute the time step for this kick */
- const integertime_t old_ti_begin =
- get_integer_time_begin(ti_current, p->time_bin);
- const integertime_t old_ti_end =
- get_integer_time_end(ti_current, p->time_bin);
- const integertime_t ti_start = (old_ti_begin + old_ti_end) / 2;
- const integertime_t ti_end = old_ti_end + new_dti / 2;
- const float dt = (ti_end - ti_start) * timeBase;
- const float half_dt = (ti_end - old_ti_end) * timeBase;
-
- /* Move particle forward in time */
- // p->ti_begin = p->ti_end;
- // p->ti_end = p->ti_begin + new_dti;
- p->time_bin = get_time_bin(new_dti);
- // if (p->id == 116650) message("Time bin=%d new_dti=%d", p->time_bin,
- // new_dti);
- if (p->gpart != NULL) {
- // p->gpart->ti_begin = p->ti_begin;
- // p->gpart->ti_end = p->ti_end;
- p->gpart->time_bin = get_time_bin(new_dti);
- }
-
- /* Get the acceleration */
- float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]};
- if (p->gpart != NULL) {
- a_tot[0] += p->gpart->a_grav[0];
- a_tot[1] += p->gpart->a_grav[1];
- a_tot[2] += p->gpart->a_grav[2];
- }
-
- /* Kick particles in momentum space */
- xp->v_full[0] += a_tot[0] * dt;
- xp->v_full[1] += a_tot[1] * dt;
- xp->v_full[2] += a_tot[2] * dt;
- if (p->gpart != NULL) {
- p->gpart->v_full[0] = xp->v_full[0];
- p->gpart->v_full[1] = xp->v_full[1];
- p->gpart->v_full[2] = xp->v_full[2];
- }
-
- /* Go back by half-step for the hydro velocity */
- p->v[0] = xp->v_full[0] - half_dt * a_tot[0];
- p->v[1] = xp->v_full[1] - half_dt * a_tot[1];
- p->v[2] = xp->v_full[2] - half_dt * a_tot[2];
-
- /* Extra kick work */
- hydro_kick_extra(p, xp, dt, half_dt);
- if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt, half_dt);
-}
-
-#endif
-
__attribute__((always_inline)) INLINE static void kick_part(
struct part *restrict p, struct xpart *restrict xp, integertime_t ti_start,
integertime_t ti_end, integertime_t ti_current, double timeBase) {
@@ -160,20 +92,4 @@ __attribute__((always_inline)) INLINE static void kick_part(
if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt);
}
-__attribute__((always_inline)) INLINE static void kick_gpart(
- struct gpart *restrict gp, integertime_t ti_start, integertime_t ti_end,
- integertime_t ti_current, double timeBase) {
-
- /* Time interval for this half-kick */
- const float dt = (ti_end - ti_start) * timeBase;
-
- /* Kick particles in momentum space */
- gp->v_full[0] += gp->a_grav[0] * dt;
- gp->v_full[1] += gp->a_grav[1] * dt;
- gp->v_full[2] += gp->a_grav[2] * dt;
-
- /* Kick extra variables */
- gravity_kick_extra(gp, dt);
-}
-
#endif /* SWIFT_KICK_H */
diff --git a/src/runner.c b/src/runner.c
index 964e5018f4007efbcb8c4f4b40c3b732874d51de..b7c7129f25ed70abbd1c35d577c906b2ff08c6de 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -823,6 +823,13 @@ void runner_do_drift_mapper(void *map_data, int num_elements,
}
}
+/**
+ * @brief Perform the first half-kick on all the active particles in a cell.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
@@ -865,7 +872,7 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
}
}
- /* Loop over the parts in this cell. */
+ /* Loop over the gparts in this cell. */
for (int k = 0; k < gcount; k++) {
/* Get a handle on the part. */
@@ -886,6 +893,16 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_kick1);
}
+/**
+ * @brief Perform the second half-kick on all the active particles in a cell.
+ *
+ * Also computes the next time-step of all active particles, prepare them to be
+ * drifted and update the cell's statistics.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
@@ -957,7 +974,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
}
}
- /* Loop over the particles in this cell. */
+ /* Loop over the g-particles in this cell. */
for (int k = 0; k < gcount; k++) {
/* Get a handle on the part. */
@@ -1023,155 +1040,6 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_kick2);
}
-#ifdef OLD_KICK
-
-/**
- * @brief Kick particles in momentum space and collect statistics.
- *
- * @param r The runner thread.
- * @param c The cell.
- * @param timer Are we timing this ?
- */
-void runner_do_kick(struct runner *r, struct cell *c, int timer) {
-
- const struct engine *e = r->e;
- const double timeBase = e->timeBase;
- const integertime_t ti_current = e->ti_current;
- const int count = c->count;
- const int gcount = c->gcount;
- struct part *restrict parts = c->parts;
- struct xpart *restrict xparts = c->xparts;
- struct gpart *restrict gparts = c->gparts;
- const double const_G = e->physical_constants->const_newton_G;
-
- TIMER_TIC;
-
- /* Anything to do here? */
- if (!cell_is_active(c, e)) {
- c->updated = 0;
- c->g_updated = 0;
- return;
- }
-
- int updated = 0, g_updated = 0;
- integertime_t ti_end_min = max_nr_timesteps, ti_end_max = 0;
-
- /* No children? */
- if (!c->split) {
-
- /* Loop over the g-particles and kick the active ones. */
- for (int k = 0; k < gcount; k++) {
-
- /* Get a handle on the part. */
- struct gpart *restrict gp = &gparts[k];
-
- /* If the g-particle has no counterpart and needs to be kicked */
- if (gp->id_or_neg_offset > 0) {
-
- if (gpart_is_active(gp, e)) {
-
- /* First, finish the force calculation */
- gravity_end_force(gp, const_G);
-
- /* Compute the next timestep */
- const int new_dti = get_gpart_timestep(gp, e);
-
- /* Now we have a time step, proceed with the kick */
- kick_gpart(gp, new_dti, e->ti_current, timeBase);
-
- /* Number of updated g-particles */
- g_updated++;
- }
-
- /* Minimal time for next end of time-step */
- const integertime_t ti_end =
- get_integer_time_end(ti_current, gp->time_bin);
- ti_end_min = min(ti_end, ti_end_min);
- ti_end_max = max(ti_end, ti_end_max);
- }
- }
-
- /* Now do the hydro ones... */
-
- /* Loop over the particles and kick the active ones. */
- for (int k = 0; k < count; k++) {
-
- /* Get a handle on the part. */
- struct part *restrict p = &parts[k];
- struct xpart *restrict xp = &xparts[k];
-
- /* Minimal time for next end of time-step */
- integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
-
- if (p->id == ICHECK)
- message("Particle in kick ti_end=%lld ti_current=%lld", ti_end,
- e->ti_current);
-
- /* If particle needs to be kicked */
- if (part_is_active(p, e)) {
-
- if (p->id == ICHECK) message("Particle active in kick");
-
- /* First, finish the force loop */
- hydro_end_force(p);
- if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
-
- /* Compute the next timestep (hydro condition) */
- const integertime_t new_dti = get_part_timestep(p, xp, e);
-
- if (p->id == ICHECK)
- message("time_step=%lld (%e)", new_dti, new_dti * e->timeBase);
-
- /* Now we have a time step, proceed with the kick */
- kick_part(p, xp, new_dti, e->ti_current, timeBase);
-
- // if (p->id == ICHECK) printParticle_single(p, xp);
-
- /* Number of updated particles */
- updated++;
- if (p->gpart != NULL) g_updated++;
-
- ti_end += get_integer_timestep(p->time_bin);
- }
-
- if (p->id == ICHECK)
- message("ti_current = %lld dti=%lld ti_end=%lld (%f)", ti_current,
- get_integer_timestep(p->time_bin), ti_end,
- ti_end * e->timeBase);
- ti_end_min = min(ti_end, ti_end_min);
- ti_end_max = max(ti_end, ti_end_max);
- }
- }
-
- /* Otherwise, aggregate data from children. */
- else {
-
- /* Loop over the progeny. */
- for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL) {
- struct cell *restrict cp = c->progeny[k];
-
- /* Recurse */
- runner_do_kick(r, cp, 0);
-
- /* And aggregate */
- updated += cp->updated;
- g_updated += cp->g_updated;
- ti_end_min = min(cp->ti_end_min, ti_end_min);
- ti_end_max = max(cp->ti_end_max, ti_end_max);
- }
- }
-
- /* Store the values. */
- c->updated = updated;
- c->g_updated = g_updated;
- c->ti_end_min = ti_end_min;
- c->ti_end_max = ti_end_max;
-
- if (timer) TIMER_TOC(timer_kick);
-}
-#endif
-
/**
* @brief Construct the cell properties from the received particles
*
diff --git a/src/timestep.h b/src/timestep.h
index 10e450987e7b02f84243a9001e903e6930d07b2c..63a21bf3e2c82344f9a0b1ac5f4ac769af10c7d1 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -37,22 +37,15 @@
*/
__attribute__((always_inline)) INLINE static integertime_t
make_integer_timestep(float new_dt, timebin_t old_bin, integertime_t ti_current,
- double timeBase_inv, int verbose) {
+ double timeBase_inv) {
/* Convert to integer time */
integertime_t new_dti = (integertime_t)(new_dt * timeBase_inv);
- /* if (verbose) message("new_dti=%lld", new_dti); */
-
/* Current time-step */
integertime_t current_dti = get_integer_timestep(old_bin);
integertime_t ti_end = get_integer_time_end(ti_current, old_bin);
- /* if (verbose) */
- /* message("current_dti=%lld old_bin=%d ti_end=%lld", current_dti, old_bin,
- */
- /* ti_end); */
-
/* Limit timestep increase */
if (old_bin > 0) new_dti = min(new_dti, 2 * current_dti);
@@ -61,15 +54,10 @@ make_integer_timestep(float new_dt, timebin_t old_bin, integertime_t ti_current,
while (new_dti < dti_timeline) dti_timeline /= 2LL;
new_dti = dti_timeline;
- /* if (verbose) message("new_dti=%lld", new_dti); */
-
/* Make sure we are allowed to increase the timestep size */
if (new_dti > current_dti) {
if ((max_nr_timesteps - ti_end) % new_dti > 0) new_dti = current_dti;
}
-
- /* if (verbose) message("new_dti=%lld", new_dti); */
-
return new_dti;
}
@@ -97,7 +85,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_gpart_timestep(
/* Convert to integer time */
const integertime_t new_dti = make_integer_timestep(
- new_dt, gp->time_bin, e->ti_current, e->timeBase_inv, 0);
+ new_dt, gp->time_bin, e->ti_current, e->timeBase_inv);
return new_dti;
}
@@ -151,13 +139,9 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
new_dt = min(new_dt, e->dt_max);
new_dt = max(new_dt, e->dt_min);
- /* if (p->id == ICHECK) message("new_dt=%e", new_dt); */
-
/* Convert to integer time */
const integertime_t new_dti = make_integer_timestep(
- new_dt, p->time_bin, e->ti_current, e->timeBase_inv, p->id == ICHECK);
-
- /* if (p->id == ICHECK) message("new_dti=%lld", new_dti); */
+ new_dt, p->time_bin, e->ti_current, e->timeBase_inv);
return new_dti;
}