diff --git a/src/runner.c b/src/runner.c
index 7b9872bad3e4fcb9df67983164edb95e808523aa..341370d8cb05b4f5b2721650bf59d575856c84aa 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -906,7 +906,224 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
}
/**
- * @brief Kick particles in momentum space and collect statistics
+ * @brief Kick particles in momentum space and collect statistics (fixed
+ * time-step case)
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
+
+ const double global_dt = r->e->dt_max;
+ const double timeBase = r->e->timeBase;
+ const int count = c->count;
+ const int gcount = c->gcount;
+ struct part *const parts = c->parts;
+ struct xpart *const xparts = c->xparts;
+ struct gpart *const gparts = c->gparts;
+
+ int updated = 0, g_updated = 0;
+ int ti_end_min = max_nr_timesteps, ti_end_max = 0;
+ double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
+ float mom[3] = {0.0f, 0.0f, 0.0f};
+ float ang[3] = {0.0f, 0.0f, 0.0f};
+
+ TIMER_TIC
+
+#ifdef TASK_VERBOSE
+ OUT;
+#endif
+
+ /* The new time-step */
+ const int new_dti = global_dt / timeBase;
+
+ /* 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 *const gp = &gparts[k];
+
+ /* If the g-particle has no counterpart and needs to be kicked */
+ if (gp->id < 0) {
+
+ /* First, finish the force calculation */
+ gravity_end_force(gp);
+
+ /* Compute the time step for this kick */
+ const int ti_start = (gp->ti_begin + gp->ti_end) / 2;
+ const int ti_end = gp->ti_end + new_dti / 2;
+ const double dt = (ti_end - ti_start) * timeBase;
+ const double half_dt = (ti_end - gp->ti_end) * timeBase;
+
+ /* Move particle forward in time */
+ gp->ti_begin = gp->ti_end;
+ gp->ti_end = gp->ti_begin + 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);
+
+ /* Number of updated g-particles */
+ g_updated++;
+
+ /* Minimal time for next end of time-step */
+ ti_end_min = min(gp->ti_end, ti_end_min);
+ ti_end_max = max(gp->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 *const p = &parts[k];
+ struct xpart *const xp = &xparts[k];
+
+ /* First, finish the force loop */
+ p->h_dt *= p->h * 0.333333333f;
+
+ /* And do the same of the extra variable */
+ hydro_end_force(p);
+ if (p->gpart != NULL) gravity_end_force(p->gpart);
+
+ /* Compute the time step for this kick */
+ const int ti_start = (p->ti_begin + p->ti_end) / 2;
+ const int ti_end = p->ti_end + new_dti / 2;
+ const double dt = (ti_end - ti_start) * timeBase;
+ const double half_dt = (ti_end - p->ti_end) * timeBase;
+
+ /* Move particle forward in time */
+ p->ti_begin = p->ti_end;
+ p->ti_end = p->ti_begin + new_dti;
+ if (p->gpart != NULL) {
+ p->gpart->ti_begin = p->ti_begin;
+ p->gpart->ti_end = p->ti_end;
+ }
+
+ /* 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[1] += 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);
+
+ /* Number of updated particles */
+ updated++;
+ if (p->gpart != NULL) g_updated++;
+
+ /* Now collect quantities for statistics */
+
+ const double x[3] = {p->x[0], p->x[1], p->x[2]};
+ const float v_full[3] = {xp->v_full[0], xp->v_full[1], xp->v_full[2]};
+ const float m = p->mass;
+
+ /* Collect mass */
+ mass += m;
+
+ /* Collect momentum */
+ mom[0] += m * v_full[0];
+ mom[1] += m * v_full[1];
+ mom[2] += m * v_full[2];
+
+ /* Collect angular momentum */
+ ang[0] += m * (x[1] * v_full[2] - x[2] * v_full[1]);
+ ang[1] += m * (x[2] * v_full[0] - x[0] * v_full[2]);
+ ang[2] += m * (x[0] * v_full[1] - x[1] * v_full[0]);
+
+ /* Collect total energy. */
+ e_kin += 0.5 * m * (v_full[0] * v_full[0] + v_full[1] * v_full[1] +
+ v_full[2] * v_full[2]);
+ e_pot += 0.f; /* No gravitational potential thus far */
+ e_int += hydro_get_internal_energy(p);
+
+ /* Minimal time for next end of time-step */
+ ti_end_min = min(p->ti_end, ti_end_min);
+ ti_end_max = max(p->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 *const cp = c->progeny[k];
+
+ /* Recurse */
+ runner_do_kick_fixdt(r, cp, 0);
+
+ /* And aggregate */
+ updated += cp->updated;
+ g_updated += cp->g_updated;
+ e_kin += cp->e_kin;
+ e_int += cp->e_int;
+ e_pot += cp->e_pot;
+ mass += cp->mass;
+ mom[0] += cp->mom[0];
+ mom[1] += cp->mom[1];
+ mom[2] += cp->mom[2];
+ ang[0] += cp->ang[0];
+ ang[1] += cp->ang[1];
+ ang[2] += cp->ang[2];
+ 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->e_kin = e_kin;
+ c->e_int = e_int;
+ c->e_pot = e_pot;
+ c->mass = mass;
+ c->mom[0] = mom[0];
+ c->mom[1] = mom[1];
+ c->mom[2] = mom[2];
+ c->ang[0] = ang[0];
+ c->ang[1] = ang[1];
+ c->ang[2] = ang[2];
+ c->ti_end_min = ti_end_min;
+ c->ti_end_max = ti_end_max;
+
+ if (timer) TIMER_TOC(timer_kick);
+}
+
+/**
+ * @brief Kick particles in momentum space and collect statistics (floating
+ * time-step case)
*
* @param r The runner thread.
* @param c The cell.
@@ -914,8 +1131,8 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
*/
void runner_do_kick(struct runner *r, struct cell *c, int timer) {
- const float global_dt_min = r->e->dt_min;
- const float global_dt_max = r->e->dt_max;
+ const double global_dt_min = r->e->dt_min;
+ const double global_dt_max = r->e->dt_max;
const int ti_current = r->e->ti_current;
const double timeBase = r->e->timeBase;
const double timeBase_inv = 1.0 / r->e->timeBase;
@@ -928,8 +1145,6 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
const struct hydro_props *hydro_properties = r->e->hydro_properties;
const struct phys_const *constants = r->e->physical_constants;
const float ln_max_h_change = hydro_properties->log_max_h_change;
- const int is_fixdt =
- (r->e->policy & engine_policy_fixdt) == engine_policy_fixdt;
int updated = 0, g_updated = 0;
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
@@ -955,57 +1170,46 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
/* If the g-particle has no counterpart and needs to be kicked */
if (gp->id < 0) {
- if (is_fixdt || gp->ti_end <= ti_current) {
+ if (gp->ti_end <= ti_current) {
/* First, finish the force calculation */
gravity_end_force(gp);
- /* Now we are ready to compute the next time-step size */
- int new_dti;
-
- if (is_fixdt) {
-
- /* Now we have a time step, proceed with the kick */
- new_dti = global_dt_max * timeBase_inv;
-
- } else {
-
- /* Compute the next timestep (gravity condition) */
- const float new_dt_external =
- gravity_compute_timestep_external(potential, constants, gp);
- const float new_dt_self =
- gravity_compute_timestep_self(constants, gp);
-
- float new_dt = fminf(new_dt_external, new_dt_self);
+ /* Compute the next timestep (gravity condition) */
+ const float new_dt_external =
+ gravity_compute_timestep_external(potential, constants, gp);
+ const float new_dt_self =
+ gravity_compute_timestep_self(constants, gp);
- /* Limit timestep within the allowed range */
- new_dt = fminf(new_dt, global_dt_max);
- new_dt = fmaxf(new_dt, global_dt_min);
+ float new_dt = fminf(new_dt_external, new_dt_self);
- /* Convert to integer time */
- new_dti = new_dt * timeBase_inv;
+ /* Limit timestep within the allowed range */
+ new_dt = fminf(new_dt, global_dt_max);
+ new_dt = fmaxf(new_dt, global_dt_min);
- /* Recover the current timestep */
- const int current_dti = gp->ti_end - gp->ti_begin;
+ /* Convert to integer time */
+ int new_dti = new_dt * timeBase_inv;
- /* Limit timestep increase */
- if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
+ /* Recover the current timestep */
+ const int current_dti = gp->ti_end - gp->ti_begin;
- /* Put this timestep on the time line */
- int dti_timeline = max_nr_timesteps;
- while (new_dti < dti_timeline) dti_timeline /= 2;
+ /* Limit timestep increase */
+ if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
- new_dti = dti_timeline;
+ /* Put this timestep on the time line */
+ int dti_timeline = max_nr_timesteps;
+ while (new_dti < dti_timeline) dti_timeline /= 2;
- /* Make sure we are allowed to increase the timestep size */
- if (new_dti > current_dti) {
- if ((max_nr_timesteps - gp->ti_end) % new_dti > 0)
- new_dti = current_dti;
- }
+ new_dti = dti_timeline;
- /* Now we have a time step, proceed with the kick */
+ /* Make sure we are allowed to increase the timestep size */
+ if (new_dti > current_dti) {
+ if ((max_nr_timesteps - gp->ti_end) % new_dti > 0)
+ new_dti = current_dti;
}
+ /* Now we have a time step, proceed with the kick */
+
/* Compute the time step for this kick */
const int ti_start = (gp->ti_begin + gp->ti_end) / 2;
const int ti_end = gp->ti_end + new_dti / 2;
@@ -1044,7 +1248,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
struct xpart *const xp = &xparts[k];
/* If particle needs to be kicked */
- if (is_fixdt || p->ti_end <= ti_current) {
+ if (p->ti_end <= ti_current) {
/* First, finish the force loop */
p->h_dt *= p->h * 0.333333333f;
@@ -1053,70 +1257,59 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
hydro_end_force(p);
if (p->gpart != NULL) gravity_end_force(p->gpart);
- /* Now we are ready to compute the next time-step size */
- int new_dti;
-
- if (is_fixdt) {
-
- /* Now we have a time step, proceed with the kick */
- new_dti = global_dt_max * timeBase_inv;
+ /* Compute the next timestep (hydro condition) */
+ const float new_dt_hydro =
+ hydro_compute_timestep(p, xp, hydro_properties);
- } else {
+ /* Compute the next timestep (gravity condition) */
+ float new_dt_grav = FLT_MAX;
+ if (p->gpart != NULL) {
- /* Compute the next timestep (hydro condition) */
- const float new_dt_hydro =
- hydro_compute_timestep(p, xp, hydro_properties);
+ const float new_dt_external =
+ gravity_compute_timestep_external(potential, constants, p->gpart);
+ const float new_dt_self =
+ gravity_compute_timestep_self(constants, p->gpart);
- /* Compute the next timestep (gravity condition) */
- float new_dt_grav = FLT_MAX;
- if (p->gpart != NULL) {
+ new_dt_grav = fminf(new_dt_external, new_dt_self);
+ }
- const float new_dt_external = gravity_compute_timestep_external(
- potential, constants, p->gpart);
- const float new_dt_self =
- gravity_compute_timestep_self(constants, p->gpart);
+ /* Final time-step is minimum of hydro and gravity */
+ float new_dt = fminf(new_dt_hydro, new_dt_grav);
- new_dt_grav = fminf(new_dt_external, new_dt_self);
- }
+ /* Limit change in h */
+ const float dt_h_change = (p->h_dt != 0.0f)
+ ? fabsf(ln_max_h_change * p->h / p->h_dt)
+ : FLT_MAX;
- /* Final time-step is minimum of hydro and gravity */
- float new_dt = fminf(new_dt_hydro, new_dt_grav);
+ new_dt = fminf(new_dt, dt_h_change);
- /* Limit change in h */
- const float dt_h_change =
- (p->h_dt != 0.0f) ? fabsf(ln_max_h_change * p->h / p->h_dt)
- : FLT_MAX;
+ /* Limit timestep within the allowed range */
+ new_dt = fminf(new_dt, global_dt_max);
+ new_dt = fmaxf(new_dt, global_dt_min);
- new_dt = fminf(new_dt, dt_h_change);
+ /* Convert to integer time */
+ int new_dti = new_dt * timeBase_inv;
- /* Limit timestep within the allowed range */
- new_dt = fminf(new_dt, global_dt_max);
- new_dt = fmaxf(new_dt, global_dt_min);
+ /* Recover the current timestep */
+ const int current_dti = p->ti_end - p->ti_begin;
- /* Convert to integer time */
- new_dti = new_dt * timeBase_inv;
+ /* Limit timestep increase */
+ if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
- /* Recover the current timestep */
- const int current_dti = p->ti_end - p->ti_begin;
+ /* Put this timestep on the time line */
+ int dti_timeline = max_nr_timesteps;
+ while (new_dti < dti_timeline) dti_timeline /= 2;
- /* Limit timestep increase */
- if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
+ new_dti = dti_timeline;
- /* Put this timestep on the time line */
- int dti_timeline = max_nr_timesteps;
- while (new_dti < dti_timeline) dti_timeline /= 2;
-
- new_dti = dti_timeline;
-
- /* Make sure we are allowed to increase the timestep size */
- if (new_dti > current_dti) {
- if ((max_nr_timesteps - p->ti_end) % new_dti > 0)
- new_dti = current_dti;
- }
-
- /* Now we have a time step, proceed with the kick */
+ /* Make sure we are allowed to increase the timestep size */
+ if (new_dti > current_dti) {
+ if ((max_nr_timesteps - p->ti_end) % new_dti > 0)
+ new_dti = current_dti;
}
+ /* Now we have a time step, proceed with the kick */
+
/* Compute the time step for this kick */
const int ti_start = (p->ti_begin + p->ti_end) / 2;
const int ti_end = p->ti_end + new_dti / 2;
@@ -1372,7 +1565,7 @@ void *runner_main(void *data) {
runner_do_kick(r, ci, 1);
break;
case task_type_kick_fixdt:
- runner_do_kick(r, ci, 1);
+ runner_do_kick_fixdt(r, ci, 1);
break;
case task_type_send:
break;
diff --git a/src/runner.h b/src/runner.h
index f53822cdbf9608a473fc72e6a2d049cfd6d3c5a2..35e5f56a7b94145eec656b02a1c5568ec39352b6 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -51,6 +51,7 @@ void runner_do_ghost(struct runner *r, struct cell *c);
void runner_do_sort(struct runner *r, struct cell *c, int flag, int clock);
void runner_do_gsort(struct runner *r, struct cell *c, int flag, int clock);
void runner_do_kick(struct runner *r, struct cell *c, int timer);
+void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer);
void runner_do_drift(struct runner *r, struct cell *c, int timer);
void runner_do_init(struct runner *r, struct cell *c, int timer);
void *runner_main(void *data);