diff --git a/examples/main.c b/examples/main.c
index 2f76b7098e37141771cb58daf46df972a1b3a119..9523af49ed30c54d256d287ea2846a854650bc05 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -544,8 +544,8 @@ int main(int argc, char *argv[]) {
/* Legend */
if (myrank == 0)
printf(
- "# Step Time time-step Number of updates CPU Wall-clock time "
- "[%s]\n",
+ "# Step Time time-step Number of updates Number of updates "
+ "CPU Wall-clock time [%s]\n",
clocks_getunit());
/* Let loose a runner on the space. */
diff --git a/src/cell.h b/src/cell.h
index f0e99e8b45b0225764267f9cde5eb2d3a244136f..a471eac44bfd3533c4220ab8c5ff2ddec724e87f 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -141,7 +141,7 @@ struct cell {
double mass, e_pot, e_int, e_kin;
/* Number of particles updated in this cell. */
- int updated;
+ int updated, g_updated;
/* Linking pointer for "memory management". */
struct cell *next;
diff --git a/src/engine.c b/src/engine.c
index 86de519283c94cafa6f7dead757d34e726902f2e..4069875cc4de958ceb0d8689b916acddf6738db6 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1558,6 +1558,7 @@ int engine_marktasks(struct engine *e) {
else if (t->type == task_type_kick) {
t->skip = (t->ci->ti_end_min > ti_end);
t->ci->updated = 0;
+ t->ci->g_updated = 0;
}
/* Drift? */
@@ -1743,7 +1744,7 @@ void engine_collect_kick(struct cell *c) {
if (c->kick != NULL) return;
/* Counters for the different quantities. */
- int updated = 0;
+ int updated = 0, g_updated = 0;
double e_kin = 0.0, e_int = 0.0, e_pot = 0.0;
float mom[3] = {0.0f, 0.0f, 0.0f}, ang[3] = {0.0f, 0.0f, 0.0f};
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
@@ -1766,6 +1767,7 @@ void engine_collect_kick(struct cell *c) {
ti_end_min = min(ti_end_min, cp->ti_end_min);
ti_end_max = max(ti_end_max, cp->ti_end_max);
updated += cp->updated;
+ g_updated += cp->g_updated;
e_kin += cp->e_kin;
e_int += cp->e_int;
e_pot += cp->e_pot;
@@ -1783,6 +1785,7 @@ void engine_collect_kick(struct cell *c) {
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
c->updated = updated;
+ c->g_updated = g_updated;
c->e_kin = e_kin;
c->e_int = e_int;
c->e_pot = e_pot;
@@ -1928,7 +1931,7 @@ void engine_init_particles(struct engine *e) {
*/
void engine_step(struct engine *e) {
- int updates = 0;
+ int updates = 0, g_updates = 0;
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
double e_pot = 0.0, e_int = 0.0, e_kin = 0.0;
float mom[3] = {0.0, 0.0, 0.0};
@@ -1955,6 +1958,7 @@ void engine_step(struct engine *e) {
e_int += c->e_int;
e_pot += c->e_pot;
updates += c->updated;
+ g_updates += c->g_updated;
mom[0] += c->mom[0];
mom[1] += c->mom[1];
mom[2] += c->mom[2];
@@ -1980,18 +1984,20 @@ void engine_step(struct engine *e) {
ti_end_max = in_i[0];
}
{
- double in_d[4], out_d[4];
+ double in_d[5], out_d[5];
out_d[0] = updates;
- out_d[1] = e_kin;
- out_d[2] = e_int;
- out_d[3] = e_pot;
- if (MPI_Allreduce(out_d, in_d, 4, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
+ out_d[1] = g_updates;
+ out_d[2] = e_kin;
+ out_d[3] = e_int;
+ out_d[4] = e_pot;
+ if (MPI_Allreduce(out_d, in_d, 5, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate energies.");
updates = in_d[0];
- e_kin = in_d[1];
- e_int = in_d[2];
- e_pot = in_d[3];
+ g_updates = in_d[1];
+ e_kin = in_d[2];
+ e_int = in_d[3];
+ e_pot = in_d[4];
}
#endif
@@ -2016,8 +2022,8 @@ void engine_step(struct engine *e) {
if (e->nodeID == 0) {
/* Print some information to the screen */
- printf("%d %e %e %d %.3f\n", e->step, e->time, e->timeStep, updates,
- e->wallclock_time);
+ printf("%d %e %e %d %d %.3f\n", e->step, e->time, e->timeStep, updates,
+ g_updates, e->wallclock_time);
fflush(stdout);
/* Write some energy statistics */
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 41d01653fc17583b62c592237948c9d2f778fb84..8945b05f4a7cb7667e4095f9604cc2ed25ccbbed 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -27,14 +27,13 @@
*
*/
-__attribute__((always_inline)) INLINE static float gravity_compute_timestep(
- struct part* p, struct xpart* xp) {
+__attribute__((always_inline))
+ INLINE static float gravity_compute_timestep(struct gpart* gp) {
/* Currently no limit is imposed */
return FLT_MAX;
}
-
/**
* @brief Initialises the g-particles for the first time
*
@@ -44,8 +43,7 @@ __attribute__((always_inline)) INLINE static float gravity_compute_timestep(
* @param gp The particle to act upon
*/
__attribute__((always_inline))
- INLINE static void gravity_first_init_gpart(struct gpart* gp) {
-}
+ INLINE static void gravity_first_init_gpart(struct gpart* gp) {}
/**
* @brief Prepares a g-particle for the gravity calculation
@@ -63,3 +61,23 @@ __attribute__((always_inline))
gp->a_grav[1] = 0.f;
gp->a_grav[2] = 0.f;
}
+
+/**
+ * @brief Finishes the gravity calculation.
+ *
+ * Multiplies the forces and accelerations by the appropiate constants
+ *
+ * @param gp The particle to act upon
+ */
+__attribute__((always_inline))
+ INLINE static void gravity_end_force(struct gpart* gp) {}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @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, float half_dt) {}
diff --git a/src/runner.c b/src/runner.c
index dea26533fd7eab10bd1263a76945cc3cdaadc74f..65a14b1a8bed01992d844b17ffb7017005a43f9f 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -788,27 +788,91 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
const double timeBase = r->e->timeBase;
const double timeBase_inv = 1.0 / r->e->timeBase;
const int count = c->count;
- // const int gcount = c->gcount;
+ const int gcount = c->gcount;
+ struct part *const parts = c->parts;
+ struct xpart *const xparts = c->xparts;
+ struct gpart *const gparts = c->gparts;
const int is_fixdt =
(r->e->policy & engine_policy_fixdt) == engine_policy_fixdt;
- int new_dti;
- int dti_timeline;
-
- int updated = 0;
+ 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};
- struct part *const parts = c->parts;
- struct xpart *const xparts = c->xparts;
- // struct gpart *const gparts = c->gparts;
TIMER_TIC
/* 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 && (is_fixdt || 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) */
+ float new_dt = gravity_compute_timestep(gp);
+
+ /* Limit timestep within the allowed range */
+ new_dt = fminf(new_dt, global_dt_max);
+ new_dt = fmaxf(new_dt, global_dt_min);
+
+ /* Convert to integer time */
+ new_dti = new_dt * timeBase_inv;
+
+ /* Recover the current timestep */
+ const int current_dti = gp->ti_end - gp->ti_begin;
+
+ /* Limit timestep increase */
+ if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
+
+ /* Put this timestep on the time line */
+ int dti_timeline = max_nr_timesteps;
+ while (new_dti < dti_timeline) dti_timeline /= 2;
+
+ /* Now we have a time step, proceed with the kick */
+ new_dti = dti_timeline;
+ }
+
+ /* 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;
+
+ /* 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++;
+ }
+ }
+
+ /* Now do the hydro ones... */
+
/* Loop over the particles and kick the active ones. */
for (int k = 0; k < count; k++) {
@@ -824,8 +888,10 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* And do the same of the extra variable */
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) {
@@ -834,9 +900,13 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
} else {
- /* Compute the next timestep */
+ /* Compute the next timestep (hydro condition) */
const float new_dt_hydro = hydro_compute_timestep(p, xp);
- const float new_dt_grav = gravity_compute_timestep(p, xp);
+
+ /* Compute the next timestep (gravity condition) */
+ float new_dt_grav = FLT_MAX;
+ if (p->gpart != NULL)
+ new_dt_grav = gravity_compute_timestep(p->gpart);
float new_dt = fminf(new_dt_hydro, new_dt_grav);
@@ -861,7 +931,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
/* Put this timestep on the time line */
- dti_timeline = max_nr_timesteps;
+ int dti_timeline = max_nr_timesteps;
while (new_dti < dti_timeline) dti_timeline /= 2;
/* Now we have a time step, proceed with the kick */
@@ -871,28 +941,38 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* 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 float dt = (ti_end - ti_start) * timeBase;
- const float half_dt = (ti_end - p->ti_end) * timeBase;
+ 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;
+ /* 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] += p->a_hydro[0] * dt;
- xp->v_full[1] += p->a_hydro[1] * dt;
- xp->v_full[2] += p->a_hydro[2] * dt;
+ xp->v_full[0] += a_tot[0] * dt;
+ xp->v_full[1] += a_tot[1] * dt;
+ xp->v_full[2] += a_tot[2] * dt;
/* Go back by half-step for the hydro velocity */
- p->v[0] = xp->v_full[0] - half_dt * p->a_hydro[0];
- p->v[1] = xp->v_full[1] - half_dt * p->a_hydro[1];
- p->v[2] = xp->v_full[2] - half_dt * p->a_hydro[2];
+ 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 */
@@ -940,6 +1020,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* 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;
@@ -957,6 +1038,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* 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;