diff --git a/configure.ac b/configure.ac
index 4b6308e96b81bbfd0a9256bb4914f1356fbfa6f8..386a18c0dfa7cb160046f748865cf19990479892 100644
--- a/configure.ac
+++ b/configure.ac
@@ -16,7 +16,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Init the project.
-AC_INIT([SWIFT],[0.4.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
+AC_INIT([SWIFT],[0.5.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
swift_config_flags="$*"
# Need to define this, instead of using fifth argument of AC_INIT, until 2.64.
diff --git a/examples/ExternalPointMass/energy_plot.py b/examples/ExternalPointMass/energy_plot.py
index a75fcb835d33b3695170aab822092556f12db7d1..25640bcb5af2966dcd57efbe1a814bb18ac4f263 100644
--- a/examples/ExternalPointMass/energy_plot.py
+++ b/examples/ExternalPointMass/energy_plot.py
@@ -91,6 +91,9 @@ for i in range(402):
E_tot_snap[i] = E_kin_snap[i] + E_pot_snap[i]
Lz_snap[i] = np.sum(Lz)
+print "Starting energy:", E_kin_stats[0], E_pot_stats[0], E_tot_stats[0]
+print "Ending energy:", E_kin_stats[-1], E_pot_stats[-1], E_tot_stats[-1]
+
# Plot energy evolution
figure()
plot(time_stats, E_kin_stats, "r-", lw=0.5, label="Kinetic energy")
diff --git a/examples/main.c b/examples/main.c
index 3a8a7a123428abed7227da0cdd17a5f43f43fc09..3c1f4aab261bd7e031a42d408f4275d766948e86 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -334,10 +334,10 @@ int main(int argc, char *argv[]) {
MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
#endif
-/* Prepare the domain decomposition scheme */
+ /* Prepare the domain decomposition scheme */
+ enum repartition_type reparttype = REPART_NONE;
#ifdef WITH_MPI
struct partition initial_partition;
- enum repartition_type reparttype;
partition_init(&initial_partition, &reparttype, params, nr_nodes);
/* Let's report what we did */
@@ -515,8 +515,9 @@ int main(int argc, char *argv[]) {
if (myrank == 0) clocks_gettime(&tic);
struct engine e;
engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff,
- engine_policies, talking, &us, &prog_const, &hydro_properties,
- &gravity_properties, &potential, &cooling_func, &sourceterms);
+ engine_policies, talking, reparttype, &us, &prog_const,
+ &hydro_properties, &gravity_properties, &potential, &cooling_func,
+ &sourceterms);
if (myrank == 0) {
clocks_gettime(&toc);
message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
@@ -590,11 +591,6 @@ int main(int argc, char *argv[]) {
/* Main simulation loop */
for (int j = 0; !engine_is_done(&e) && e.step != nsteps; j++) {
-/* Repartition the space amongst the nodes? */
-#ifdef WITH_MPI
- if (j % 100 == 2) e.forcerepart = reparttype;
-#endif
-
/* Reset timers */
timers_reset(timers_mask_all);
@@ -696,6 +692,7 @@ int main(int argc, char *argv[]) {
#endif
/* Write final output. */
+ engine_drift_all(&e);
engine_dump_snapshot(&e);
#ifdef WITH_MPI
diff --git a/src/active.h b/src/active.h
index 0c22a745fed4fbdf72ef1377fad45b78c86f178f..38020fde8ed53a638231097643476b7ef72d6d49 100644
--- a/src/active.h
+++ b/src/active.h
@@ -50,8 +50,10 @@ __attribute__((always_inline)) INLINE static int cell_is_drifted(
return (c->ti_old == e->ti_current);
}
+/* Are cells / particles active for regular tasks ? */
+
/**
- * @brief Does a cell contain any active particle ?
+ * @brief Does a cell contain any particle finishing their time-step now ?
*
* @param c The #cell.
* @param e The #engine containing information about the current time.
@@ -73,7 +75,7 @@ __attribute__((always_inline)) INLINE static int cell_is_active(
}
/**
- * @brief Are *all* particles in a cell active ?
+ * @brief Are *all* particles in a cell finishing their time-step now ?
*
* @param c The #cell.
* @param e The #engine containing information about the current time.
@@ -94,7 +96,7 @@ __attribute__((always_inline)) INLINE static int cell_is_all_active(
}
/**
- * @brief Is this particle active ?
+ * @brief Is this particle finishing its time-step now ?
*
* @param p The #part.
* @param e The #engine containing information about the current time.
@@ -118,7 +120,7 @@ __attribute__((always_inline)) INLINE static int part_is_active(
}
/**
- * @brief Is this g-particle active ?
+ * @brief Is this g-particle finishing its time-step now ?
*
* @param gp The #gpart.
* @param e The #engine containing information about the current time.
@@ -142,7 +144,7 @@ __attribute__((always_inline)) INLINE static int gpart_is_active(
}
/**
- * @brief Is this s-particle active ?
+ * @brief Is this s-particle finishing its time-step now ?
*
* @param sp The #spart.
* @param e The #engine containing information about the current time.
@@ -157,7 +159,7 @@ __attribute__((always_inline)) INLINE static int spart_is_active(
#ifdef SWIFT_DEBUG_CHECKS
if (ti_end < ti_current)
error(
- "s-particle in an impossible time-zone! gp->ti_end=%lld "
+ "s-particle in an impossible time-zone! sp->ti_end=%lld "
"e->ti_current=%lld",
ti_end, ti_current);
#endif
@@ -165,4 +167,102 @@ __attribute__((always_inline)) INLINE static int spart_is_active(
return (ti_end == ti_current);
}
+/* Are cells / particles active for kick1 tasks ? */
+
+/**
+ * @brief Does a cell contain any particle starting their time-step now ?
+ *
+ * @param c The #cell.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #cell contains at least an active particle, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int cell_is_starting(
+ const struct cell *c, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->ti_beg_max > e->ti_current)
+ error(
+ "cell in an impossible time-zone! c->ti_beg_max=%lld (t=%e) and "
+ "e->ti_current=%lld (t=%e)",
+ c->ti_beg_max, c->ti_beg_max * e->timeBase, e->ti_current,
+ e->ti_current * e->timeBase);
+#endif
+
+ return (c->ti_beg_max == e->ti_current);
+}
+
+/**
+ * @brief Is this particle starting its time-step now ?
+ *
+ * @param p The #part.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #part is active, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int part_is_starting(
+ const struct part *p, const struct engine *e) {
+
+ const integertime_t ti_current = e->ti_current;
+ const integertime_t ti_beg =
+ get_integer_time_begin(ti_current + 1, p->time_bin);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ti_beg > ti_current)
+ error(
+ "particle in an impossible time-zone! p->ti_beg=%lld "
+ "e->ti_current=%lld",
+ ti_beg, ti_current);
+#endif
+
+ return (ti_beg == ti_current);
+}
+
+/**
+ * @brief Is this g-particle starting its time-step now ?
+ *
+ * @param gp The #gpart.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #gpart is active, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int gpart_is_starting(
+ const struct gpart *gp, const struct engine *e) {
+
+ const integertime_t ti_current = e->ti_current;
+ const integertime_t ti_beg =
+ get_integer_time_begin(ti_current + 1, gp->time_bin);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ti_beg > ti_current)
+ error(
+ "g-particle in an impossible time-zone! gp->ti_beg=%lld "
+ "e->ti_current=%lld",
+ ti_beg, ti_current);
+#endif
+
+ return (ti_beg == ti_current);
+}
+
+/**
+ * @brief Is this s-particle starting its time-step now ?
+ *
+ * @param sp The #spart.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #spart is active, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int spart_is_starting(
+ const struct spart *sp, const struct engine *e) {
+
+ const integertime_t ti_current = e->ti_current;
+ const integertime_t ti_beg =
+ get_integer_time_begin(ti_current + 1, sp->time_bin);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ti_beg > ti_current)
+ error(
+ "s-particle in an impossible time-zone! sp->ti_beg=%lld "
+ "e->ti_current=%lld",
+ ti_beg, ti_current);
+#endif
+
+ return (ti_beg == ti_current);
+}
#endif /* SWIFT_ACTIVE_H */
diff --git a/src/cell.c b/src/cell.c
index 941ee7a7623b0006f0291878da58f7c8fa966565..c5f414817e20f5857fd836cb86b99aea923f0b2b 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1026,11 +1026,34 @@ void cell_clean_links(struct cell *c, void *data) {
*/
void cell_check_drift_point(struct cell *c, void *data) {
- integertime_t ti_current = *(integertime_t *)data;
+#ifdef SWIFT_DEBUG_CHECKS
+
+ const integertime_t ti_drift = *(integertime_t *)data;
+
+ /* Only check local cells */
+ if (c->nodeID != engine_rank) return;
+
+ if (c->ti_old != ti_drift)
+ error("Cell in an incorrect time-zone! c->ti_old=%lld ti_drift=%lld",
+ c->ti_old, ti_drift);
+
+ for (int i = 0; i < c->count; ++i)
+ if (c->parts[i].ti_drift != ti_drift)
+ error("part in an incorrect time-zone! p->ti_drift=%lld ti_drift=%lld",
+ c->parts[i].ti_drift, ti_drift);
+
+ for (int i = 0; i < c->gcount; ++i)
+ if (c->gparts[i].ti_drift != ti_drift)
+ error("g-part in an incorrect time-zone! gp->ti_drift=%lld ti_drift=%lld",
+ c->gparts[i].ti_drift, ti_drift);
- if (c->ti_old != ti_current && c->nodeID == engine_rank)
- error("Cell in an incorrect time-zone! c->ti_old=%lld ti_current=%lld",
- c->ti_old, ti_current);
+ for (int i = 0; i < c->scount; ++i)
+ if (c->sparts[i].ti_drift != ti_drift)
+ error("s-part in an incorrect time-zone! sp->ti_drift=%lld ti_drift=%lld",
+ c->sparts[i].ti_drift, ti_drift);
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
}
/**
@@ -1155,6 +1178,8 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
struct engine *e = s->space->e;
#endif
+ int rebuild = 0;
+
/* Un-skip the density tasks involved with this cell. */
for (struct link *l = c->density; l != NULL; l = l->next) {
struct task *t = l->t;
@@ -1180,7 +1205,7 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
(max(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin ||
ci->dx_max > space_maxreldx * ci->h_max ||
cj->dx_max > space_maxreldx * cj->h_max))
- return 1;
+ rebuild = 1;
#ifdef WITH_MPI
/* Activate the send/recv flags. */
@@ -1277,7 +1302,7 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if (c->cooling != NULL) scheduler_activate(s, c->cooling);
if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms);
- return 0;
+ return rebuild;
}
/**
@@ -1424,5 +1449,7 @@ void cell_check_timesteps(struct cell *c) {
if (c->parts[i].time_bin == 0)
error("Particle without assigned time-bin");
}
+#else
+ error("Calling debugging code without debugging flag activated.");
#endif
}
diff --git a/src/cell.h b/src/cell.h
index 0d5a46006bd34cce4e65e9c2f9761292c31f40be..6c90f38c51980fd742480c90864f2f901005c0f1 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -74,7 +74,7 @@ struct pcell {
/* Stats on this cell's particles. */
double h_max;
- integertime_t ti_end_min, ti_end_max, ti_old;
+ integertime_t ti_end_min, ti_end_max, ti_beg_max, ti_old;
/* Number of particles in this cell. */
int count, gcount, scount;
@@ -227,6 +227,9 @@ struct cell {
/*! Maximum end of (integer) time step in this cell. */
integertime_t ti_end_max;
+ /*! Maximum beginning of (integer) time step in this cell. */
+ integertime_t ti_beg_max;
+
/*! Last (integer) time the cell's content was drifted forward in time. */
integertime_t ti_old;
diff --git a/src/drift.h b/src/drift.h
index 687f8d8885a5fedca489f76d65ea8113101626c6..9957d829fd28e033df9db325186195a3b396738c 100644
--- a/src/drift.h
+++ b/src/drift.h
@@ -41,6 +41,18 @@
__attribute__((always_inline)) INLINE static void drift_gpart(
struct gpart *restrict gp, float dt, double timeBase, integertime_t ti_old,
integertime_t ti_current) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (gp->ti_drift != ti_old)
+ error(
+ "g-particle has not been drifted to the current time "
+ "gp->ti_drift=%lld, "
+ "c->ti_old=%lld, ti_current=%lld",
+ gp->ti_drift, ti_old, ti_current);
+
+ gp->ti_drift = ti_current;
+#endif
+
/* Drift... */
gp->x[0] += gp->v_full[0] * dt;
gp->x[1] += gp->v_full[1] * dt;
@@ -69,7 +81,7 @@ __attribute__((always_inline)) INLINE static void drift_part(
#ifdef SWIFT_DEBUG_CHECKS
if (p->ti_drift != ti_old)
error(
- "Particle has not been drifted to the current time p->ti_drift=%lld, "
+ "particle has not been drifted to the current time p->ti_drift=%lld, "
"c->ti_old=%lld, ti_current=%lld",
p->ti_drift, ti_old, ti_current);
@@ -108,6 +120,17 @@ __attribute__((always_inline)) INLINE static void drift_spart(
struct spart *restrict sp, float dt, double timeBase, integertime_t ti_old,
integertime_t ti_current) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (sp->ti_drift != ti_old)
+ error(
+ "s-particle has not been drifted to the current time "
+ "sp->ti_drift=%lld, "
+ "c->ti_old=%lld, ti_current=%lld",
+ sp->ti_drift, ti_old, ti_current);
+
+ sp->ti_drift = ti_current;
+#endif
+
/* Drift... */
sp->x[0] += sp->v[0] * dt;
sp->x[1] += sp->v[1] * dt;
diff --git a/src/engine.c b/src/engine.c
index 531b84c471514341736b0e846efec4e6c734e43e..fbb6da306eb86ee69c47bb86b17e86d768ac24ae 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -157,12 +157,12 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
0, 0, c, NULL, 0);
scheduler_addunlock(s, c->kick2, c->timestep);
+ scheduler_addunlock(s, c->timestep, c->kick1);
/* Add the drift task and its dependencies. */
c->drift = scheduler_addtask(s, task_type_drift, task_subtype_none, 0, 0,
c, NULL, 0);
- scheduler_addunlock(s, c->kick1, c->drift);
scheduler_addunlock(s, c->drift, c->init);
/* Generate the ghost task. */
@@ -822,19 +822,18 @@ void engine_repartition(struct engine *e) {
fflush(stdout);
/* Check that all cells have been drifted to the current time */
- space_check_drift_point(e->s, e->ti_current);
+ space_check_drift_point(e->s, e->ti_old);
#endif
/* Clear the repartition flag. */
- enum repartition_type reparttype = e->forcerepart;
- e->forcerepart = REPART_NONE;
+ e->forcerepart = 0;
/* Nothing to do if only using a single node. Also avoids METIS
* bug that doesn't handle this case well. */
if (e->nr_nodes == 1) return;
/* Do the repartitioning. */
- partition_repartition(reparttype, e->nodeID, e->nr_nodes, e->s,
+ partition_repartition(e->reparttype, e->nodeID, e->nr_nodes, e->s,
e->sched.tasks, e->sched.nr_tasks);
/* Now comes the tricky part: Exchange particles between all nodes.
@@ -2599,8 +2598,15 @@ void engine_rebuild(struct engine *e) {
if (engine_marktasks(e))
error("engine_marktasks failed after space_rebuild.");
- /* Print the status of the system */
- if (e->verbose) engine_print_task_counts(e);
+/* Print the status of the system */
+// if (e->verbose) engine_print_task_counts(e);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that all cells have been drifted to the current time.
+ * That can include cells that have not
+ * previously been active on this rank. */
+ space_check_drift_point(e->s, e->ti_old);
+#endif
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -2611,48 +2617,28 @@ void engine_rebuild(struct engine *e) {
* @brief Prepare the #engine by re-building the cells and tasks.
*
* @param e The #engine to prepare.
- * @param drift_all Whether to drift particles before rebuilding or not. Will
- * not be necessary if all particles have already been
- * drifted (before repartitioning for instance).
- * @param postrepart If we have just repartitioned, if so we need to defer the
- * skip until after the rebuild and not check the if all
- * cells have been drifted.
*/
-void engine_prepare(struct engine *e, int drift_all, int postrepart) {
+void engine_prepare(struct engine *e) {
TIMER_TIC;
- /* Unskip active tasks and check for rebuild */
- if (!postrepart) engine_unskip(e);
-
- /* Run through the tasks and mark as skip or not. */
- int rebuild = e->forcerebuild;
-
-/* Collect the values of rebuild from all nodes. */
-#ifdef WITH_MPI
- int buff = 0;
- if (MPI_Allreduce(&rebuild, &buff, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD) !=
- MPI_SUCCESS)
- error("Failed to aggregate the rebuild flag across nodes.");
- rebuild = buff;
-#endif
-
- /* And rebuild if necessary. */
- if (rebuild) {
-
- /* Drift all particles to the current time if needed. */
- if (drift_all) engine_drift_all(e);
-
#ifdef SWIFT_DEBUG_CHECKS
- /* Check that all cells have been drifted to the current time, unless
- * we have just repartitioned, that can include cells that have not
+ if (e->forcerepart || e->forcerebuild) {
+ /* Check that all cells have been drifted to the current time.
+ * That can include cells that have not
* previously been active on this rank. */
- if (!postrepart) space_check_drift_point(e->s, e->ti_current);
+ space_check_drift_point(e->s, e->ti_old);
+ }
#endif
- engine_rebuild(e);
- }
- if (postrepart) engine_unskip(e);
+ /* Do we need repartitioning ? */
+ if (e->forcerepart) engine_repartition(e);
+
+ /* Do we need rebuilding ? */
+ if (e->forcerebuild) engine_rebuild(e);
+
+ /* Unskip active tasks and check for rebuild */
+ engine_unskip(e);
/* Re-rank the tasks every now and then. */
if (e->tasks_age % engine_tasksreweight == 1) {
@@ -2663,7 +2649,7 @@ void engine_prepare(struct engine *e, int drift_all, int postrepart) {
TIMER_TOC(timer_prepare);
if (e->verbose)
- message("took %.3f %s (including drift all, rebuild and reweight).",
+ message("took %.3f %s (including unskip and reweight).",
clocks_from_ticks(getticks() - tic), clocks_getunit());
}
@@ -2722,7 +2708,7 @@ void engine_collect_kick(struct cell *c) {
/* Counters for the different quantities. */
int updated = 0, g_updated = 0, s_updated = 0;
- integertime_t ti_end_min = max_nr_timesteps;
+ integertime_t ti_end_min = max_nr_timesteps, ti_end_max = 0, ti_beg_max = 0;
/* Collect the values from the progeny. */
for (int k = 0; k < 8; k++) {
@@ -2734,6 +2720,8 @@ void engine_collect_kick(struct cell *c) {
/* And update */
ti_end_min = min(ti_end_min, cp->ti_end_min);
+ ti_end_max = max(ti_end_max, cp->ti_end_max);
+ ti_beg_max = max(ti_beg_max, cp->ti_beg_max);
updated += cp->updated;
g_updated += cp->g_updated;
s_updated += cp->s_updated;
@@ -2747,6 +2735,8 @@ void engine_collect_kick(struct cell *c) {
/* Store the collected values in the cell. */
c->ti_end_min = ti_end_min;
+ c->ti_end_max = ti_end_max;
+ c->ti_beg_max = ti_beg_max;
c->updated = updated;
c->g_updated = g_updated;
c->s_updated = s_updated;
@@ -2762,7 +2752,7 @@ void engine_collect_timestep(struct engine *e) {
const ticks tic = getticks();
int updates = 0, g_updates = 0, s_updates = 0;
- integertime_t ti_end_min = max_nr_timesteps;
+ integertime_t ti_end_min = max_nr_timesteps, ti_end_max = 0, ti_beg_max = 0;
const struct space *s = e->s;
/* Collect the cell data. */
@@ -2775,6 +2765,8 @@ void engine_collect_timestep(struct engine *e) {
/* And aggregate */
ti_end_min = min(ti_end_min, c->ti_end_min);
+ ti_end_max = max(ti_end_max, c->ti_end_max);
+ ti_beg_max = max(ti_beg_max, c->ti_beg_max);
updates += c->updated;
g_updates += c->g_updated;
s_updates += c->s_updated;
@@ -2812,6 +2804,8 @@ void engine_collect_timestep(struct engine *e) {
#endif
e->ti_end_min = ti_end_min;
+ e->ti_end_max = ti_end_max;
+ e->ti_beg_max = ti_beg_max;
e->updates = updates;
e->g_updates = g_updates;
e->s_updates = s_updates;
@@ -2885,7 +2879,7 @@ void engine_skip_force_and_kick(struct engine *e) {
*
* @param e The #engine to act on.
*/
-void engine_skip_drift_and_kick(struct engine *e) {
+void engine_skip_drift(struct engine *e) {
struct task *tasks = e->sched.tasks;
const int nr_tasks = e->sched.nr_tasks;
@@ -2895,7 +2889,7 @@ void engine_skip_drift_and_kick(struct engine *e) {
struct task *t = &tasks[i];
/* Skip everything that updates the particles */
- if (t->type == task_type_drift || t->type == task_type_kick1) t->skip = 1;
+ if (t->type == task_type_drift) t->skip = 1;
}
}
@@ -2956,13 +2950,15 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
if (e->nodeID == 0) message("Computing initial gas densities.");
- engine_prepare(e, 0, 0);
-
- engine_marktasks(e);
+ /* Construct all cells and tasks to start everything */
+ engine_rebuild(e);
/* No time integration. We just want the density and ghosts */
engine_skip_force_and_kick(e);
+ /* Print the number of active tasks ? */
+ if (e->verbose) engine_print_task_counts(e);
+
/* Now, launch the calculation */
TIMER_TIC;
engine_launch(e, e->nr_threads);
@@ -2974,7 +2970,6 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
if (e->nodeID == 0) message("Converting internal energy variable.");
/* Apply the conversion */
- // space_map_cells_pre(s, 0, cell_convert_hydro, NULL);
for (size_t i = 0; i < s->nr_parts; ++i)
hydro_convert_quantities(&s->parts[i], &s->xparts[i]);
@@ -2989,12 +2984,21 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
/* Now time to get ready for the first time-step */
if (e->nodeID == 0) message("Running initial fake time-step.");
+ /* Prepare all the tasks again for a new round */
engine_marktasks(e);
- engine_skip_drift_and_kick(e);
+ /* No drift this time */
+ engine_skip_drift(e);
+ /* Print the number of active tasks ? */
+ if (e->verbose) engine_print_task_counts(e);
+
+ /* Run the 0th time-step */
engine_launch(e, e->nr_threads);
+ /* Recover the (integer) end of the next time-step */
+ engine_collect_timestep(e);
+
clocks_gettime(&time2);
#ifdef SWIFT_DEBUG_CHECKS
@@ -3004,7 +3008,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
#endif
/* Ready to go */
- e->step = 0;
+ e->step = -1;
e->forcerebuild = 1;
e->wallclock_time = (float)clocks_diff(&time1, &time2);
@@ -3018,8 +3022,6 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
*/
void engine_step(struct engine *e) {
- double snapshot_drift_time = 0.;
-
TIMER_TIC2;
struct clocks_time time1, time2;
@@ -3027,36 +3029,13 @@ void engine_step(struct engine *e) {
e->tic_step = getticks();
- /* Recover the (integer) end of the next time-step */
- engine_collect_timestep(e);
-
- /* Check for output */
- while (e->ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0) {
-
- e->ti_old = e->ti_current;
- e->ti_current = e->ti_nextSnapshot;
- e->time = e->ti_current * e->timeBase + e->timeBegin;
- e->timeOld = e->ti_old * e->timeBase + e->timeBegin;
- e->timeStep = (e->ti_current - e->ti_old) * e->timeBase;
- snapshot_drift_time = e->timeStep;
-
- /* Drift everybody to the snapshot position */
- engine_drift_all(e);
-
- /* Dump... */
- engine_dump_snapshot(e);
-
- /* ... and find the next output time */
- engine_compute_next_snapshot_time(e);
- }
-
/* Move forward in time */
e->ti_old = e->ti_current;
e->ti_current = e->ti_end_min;
e->step += 1;
e->time = e->ti_current * e->timeBase + e->timeBegin;
e->timeOld = e->ti_old * e->timeBase + e->timeBegin;
- e->timeStep = (e->ti_current - e->ti_old) * e->timeBase + snapshot_drift_time;
+ e->timeStep = (e->ti_current - e->ti_old) * e->timeBase;
if (e->nodeID == 0) {
@@ -3072,31 +3051,59 @@ void engine_step(struct engine *e) {
fflush(e->file_timesteps);
}
- /* Drift only the necessary particles, that means all particles
- * if we are about to repartition. */
- const int repart = (e->forcerepart != REPART_NONE);
- const int drift_all = (e->policy & engine_policy_drift_all);
- if (repart || drift_all) engine_drift_all(e);
-
- /* Re-distribute the particles amongst the nodes? */
- if (repart) engine_repartition(e);
+ /* Prepare the tasks to be launched, rebuild or repartition if needed. */
+ engine_prepare(e);
- /* Prepare the space. */
- engine_prepare(e, !(drift_all || repart), repart);
+/* Repartition the space amongst the nodes? */
+#if defined(WITH_MPI) && defined(HAVE_METIS)
+ if (e->step % 100 == 2) e->forcerepart = 1;
+#endif
+ /* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
- /* Save some statistics */
+ /* Save some statistics ? */
if (e->time - e->timeLastStatistics >= e->deltaTimeStatistics) {
engine_print_stats(e);
e->timeLastStatistics += e->deltaTimeStatistics;
}
- /* Send off the runners. */
+ /* Start all the tasks. */
TIMER_TIC;
engine_launch(e, e->nr_threads);
TIMER_TOC(timer_runners);
+/* Collect the values of rebuild from all nodes. */
+#ifdef WITH_MPI
+ int buff = 0;
+ if (MPI_Allreduce(&e->forcerebuild, &buff, 1, MPI_INT, MPI_MAX,
+ MPI_COMM_WORLD) != MPI_SUCCESS)
+ error("Failed to aggregate the rebuild flag across nodes.");
+ e->forcerebuild = buff;
+#endif
+
+ /* Do we want a snapshot? */
+ if (e->ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0)
+ e->dump_snapshot = 1;
+
+ /* Drift everybody (i.e. what has not yet been drifted) */
+ /* to the current time */
+ if (e->dump_snapshot || e->forcerebuild || e->forcerepart)
+ engine_drift_all(e);
+
+ /* Write a snapshot ? */
+ if (e->dump_snapshot) {
+
+ /* Dump... */
+ engine_dump_snapshot(e);
+
+ /* ... and find the next output time */
+ engine_compute_next_snapshot_time(e);
+ }
+
+ /* Recover the (integer) end of the next time-step */
+ engine_collect_timestep(e);
+
TIMER_TOC2(timer_step);
clocks_gettime(&time2);
@@ -3353,6 +3360,13 @@ void engine_dump_snapshot(struct engine *e) {
struct clocks_time time1, time2;
clocks_gettime(&time1);
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that all cells have been drifted to the current time.
+ * That can include cells that have not
+ * previously been active on this rank. */
+ space_check_drift_point(e->s, e->ti_current);
+#endif
+
if (e->verbose) message("writing snapshot at t=%e.", e->time);
/* Dump... */
@@ -3371,6 +3385,8 @@ void engine_dump_snapshot(struct engine *e) {
e->snapshotUnits);
#endif
+ e->dump_snapshot = 0;
+
clocks_gettime(&time2);
if (e->verbose)
message("writing particle properties took %.3f %s.",
@@ -3445,6 +3461,7 @@ void engine_unpin() {
* @param with_aff use processor affinity, if supported.
* @param policy The queuing policy to use.
* @param verbose Is this #engine talkative ?
+ * @param reparttype What type of repartition algorithm are we using ?
* @param internal_units The system of units used internally.
* @param physical_constants The #phys_const used for this run.
* @param hydro The #hydro_props used for this run.
@@ -3456,6 +3473,7 @@ void engine_unpin() {
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
int nr_threads, int with_aff, int policy, int verbose,
+ enum repartition_type reparttype,
const struct unit_system *internal_units,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
@@ -3477,7 +3495,9 @@ void engine_init(struct engine *e, struct space *s,
e->proxy_ind = NULL;
e->nr_proxies = 0;
e->forcerebuild = 1;
- e->forcerepart = REPART_NONE;
+ e->forcerepart = 0;
+ e->reparttype = reparttype;
+ e->dump_snapshot = 0;
e->links = NULL;
e->nr_links = 0;
e->timeBegin = parser_get_param_double(params, "TimeIntegration:time_begin");
diff --git a/src/engine.h b/src/engine.h
index 8a719a8952c157e2c567addce2f6fc6491730eca..2fec4b3d34bea58dc8a04cf84b9b09e75f01895b 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -132,6 +132,12 @@ struct engine {
/* Minimal ti_end for the next time-step */
integertime_t ti_end_min;
+ /* Maximal ti_end for the next time-step */
+ integertime_t ti_end_max;
+
+ /* Maximal ti_beg for the next time-step */
+ integertime_t ti_beg_max;
+
/* Number of particles updated */
size_t updates, g_updates, s_updates;
@@ -180,7 +186,13 @@ struct engine {
/* Force the engine to rebuild? */
int forcerebuild;
- enum repartition_type forcerepart;
+
+ /* Force the engine to repartition ? */
+ int forcerepart;
+ enum repartition_type reparttype;
+
+ /* Need to dump a snapshot ? */
+ int dump_snapshot;
/* How many steps have we done with the same set of tasks? */
int tasks_age;
@@ -223,6 +235,7 @@ void engine_dump_snapshot(struct engine *e);
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
int nr_threads, int with_aff, int policy, int verbose,
+ enum repartition_type reparttype,
const struct unit_system *internal_units,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
@@ -231,8 +244,7 @@ void engine_init(struct engine *e, struct space *s,
const struct cooling_function_data *cooling,
struct sourceterms *sourceterms);
void engine_launch(struct engine *e, int nr_runners);
-void engine_prepare(struct engine *e, int drift_all, int postrepart);
-void engine_print(struct engine *e);
+void engine_prepare(struct engine *e);
void engine_init_particles(struct engine *e, int flag_entropy_ICs);
void engine_step(struct engine *e);
void engine_maketasks(struct engine *e);
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index a0bfee05f8b7f93cce65e8b9a3e7e322e166569d..3f97070fbb80c8920a929e6df8d93def6ac6041a 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -85,6 +85,15 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
__attribute__((always_inline)) INLINE static void gravity_kick_extra(
struct gpart* gp, float dt) {}
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param gp The particle.
+ */
+__attribute__((always_inline)) INLINE static void
+gravity_reset_predicted_values(struct gpart* gp) {}
+
/**
* @brief Initialises the g-particles for the first time
*
diff --git a/src/kick.h b/src/kick.h
index d6c85b5eab92a288f78f22fce2f03862bc34604f..7ccea7d26974297cfebc605808c4443633140ec1 100644
--- a/src/kick.h
+++ b/src/kick.h
@@ -43,6 +43,16 @@ __attribute__((always_inline)) INLINE static void kick_gpart(
/* Time interval for this half-kick */
const float dt = (ti_end - ti_start) * timeBase;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (gp->ti_kick != ti_start)
+ error(
+ "g-particle has not been kicked to the current time gp->ti_kick=%lld, "
+ "ti_start=%lld, ti_end=%lld",
+ gp->ti_kick, ti_start, ti_end);
+
+ gp->ti_kick = ti_end;
+#endif
+
/* Kick particles in momentum space */
gp->v_full[0] += gp->a_grav[0] * dt;
gp->v_full[1] += gp->a_grav[1] * dt;
@@ -71,7 +81,7 @@ __attribute__((always_inline)) INLINE static void kick_part(
#ifdef SWIFT_DEBUG_CHECKS
if (p->ti_kick != ti_start)
error(
- "Particle has not been kicked to the current time p->ti_kick=%lld, "
+ "particle has not been kicked to the current time p->ti_kick=%lld, "
"ti_start=%lld, ti_end=%lld",
p->ti_kick, ti_start, ti_end);
@@ -116,6 +126,16 @@ __attribute__((always_inline)) INLINE static void kick_spart(
/* Time interval for this half-kick */
const float dt = (ti_end - ti_start) * timeBase;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (sp->ti_kick != ti_start)
+ error(
+ "s-particle has not been kicked to the current time sp->ti_kick=%lld, "
+ "ti_start=%lld, ti_end=%lld",
+ sp->ti_kick, ti_start, ti_end);
+
+ sp->ti_kick = ti_end;
+#endif
+
/* Acceleration from gravity */
const float a[3] = {sp->gpart->a_grav[0], sp->gpart->a_grav[1],
sp->gpart->a_grav[2]};
diff --git a/src/runner.c b/src/runner.c
index 9910c77f2c31b65014334c3df183854a5b7912eb..6f2f32c1f54a904deb7805b52b1d93bfa82c57b8 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -861,7 +861,7 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
/* Anything to do here? */
- if (!cell_is_active(c, e)) return;
+ if (!cell_is_starting(c, e)) return;
/* Recurse? */
if (c->split) {
@@ -877,17 +877,17 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
struct xpart *restrict xp = &xparts[k];
/* If particle needs to be kicked */
- if (part_is_active(p, e)) {
+ if (part_is_starting(p, e)) {
const integertime_t ti_step = get_integer_timestep(p->time_bin);
const integertime_t ti_begin =
- get_integer_time_begin(ti_current, p->time_bin);
+ get_integer_time_begin(ti_current + 1, p->time_bin);
#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_end =
- get_integer_time_end(ti_current, p->time_bin);
+ get_integer_time_end(ti_current + 1, p->time_bin);
- if (ti_end - ti_begin != ti_step)
+ 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",
@@ -906,17 +906,21 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
struct gpart *restrict gp = &gparts[k];
/* If the g-particle has no counterpart and needs to be kicked */
- if (gp->type == swift_type_dark_matter && gpart_is_active(gp, e)) {
+ if (gp->type == swift_type_dark_matter && gpart_is_starting(gp, e)) {
const integertime_t ti_step = get_integer_timestep(gp->time_bin);
const integertime_t ti_begin =
- get_integer_time_begin(ti_current, gp->time_bin);
+ get_integer_time_begin(ti_current + 1, gp->time_bin);
#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_end =
- get_integer_time_end(ti_current, gp->time_bin);
+ get_integer_time_end(ti_current + 1, gp->time_bin);
- if (ti_end - ti_begin != ti_step) error("Particle in wrong 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, gp->time_bin, ti_current);
#endif
/* do the kick */
@@ -931,17 +935,21 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
struct spart *restrict sp = &sparts[k];
/* If particle needs to be kicked */
- if (spart_is_active(sp, e)) {
+ if (spart_is_starting(sp, e)) {
const integertime_t ti_step = get_integer_timestep(sp->time_bin);
const integertime_t ti_begin =
- get_integer_time_begin(ti_current, sp->time_bin);
+ get_integer_time_begin(ti_current + 1, sp->time_bin);
#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_end =
- get_integer_time_end(ti_current, sp->time_bin);
+ get_integer_time_end(ti_current + 1, sp->time_bin);
- if (ti_end - ti_begin != ti_step) error("Particle in wrong 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, sp->time_bin, ti_current);
#endif
/* do the kick */
@@ -1011,6 +1019,11 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
/* Finish the time-step with a second half-kick */
kick_part(p, xp, ti_begin + ti_step / 2, ti_begin + ti_step, timeBase);
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that kick and the drift are synchronized */
+ if (p->ti_drift != p->ti_kick) error("Error integrating part in time.");
+#endif
+
/* Prepare the values to be drifted */
hydro_reset_predicted_values(p, xp);
}
@@ -1036,6 +1049,15 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
/* Finish the time-step with a second half-kick */
kick_gpart(gp, ti_begin + ti_step / 2, ti_begin + ti_step, timeBase);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that kick and the drift are synchronized */
+ if (gp->ti_drift != gp->ti_kick)
+ error("Error integrating g-part in time.");
+#endif
+
+ /* Prepare the values to be drifted */
+ gravity_reset_predicted_values(gp);
}
}
@@ -1060,6 +1082,12 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
/* Finish the time-step with a second half-kick */
kick_spart(sp, ti_begin + ti_step / 2, ti_begin + ti_step, timeBase);
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that kick and the drift are synchronized */
+ if (sp->ti_drift != sp->ti_kick)
+ error("Error integrating s-part in time.");
+#endif
+
/* Prepare the values to be drifted */
star_reset_predicted_values(sp);
}
@@ -1091,7 +1119,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
int updated = 0, g_updated = 0, s_updated = 0;
- integertime_t ti_end_min = max_nr_timesteps, ti_end_max = 0;
+ integertime_t ti_end_min = max_nr_timesteps, ti_end_max = 0, ti_beg_max = 0;
/* No children? */
if (!c->split) {
@@ -1129,6 +1157,9 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
/* What is the next sync-point ? */
ti_end_min = min(ti_current + ti_new_step, ti_end_min);
ti_end_max = max(ti_current + ti_new_step, ti_end_max);
+
+ /* What is the next starting point for this cell ? */
+ ti_beg_max = max(ti_current, ti_beg_max);
}
else { /* part is inactive */
@@ -1139,6 +1170,12 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
/* What is the next sync-point ? */
ti_end_min = min(ti_end, ti_end_min);
ti_end_max = max(ti_end, ti_end_max);
+
+ const integertime_t ti_beg =
+ get_integer_time_begin(ti_current + 1, p->time_bin);
+
+ /* What is the next starting point for this cell ? */
+ ti_beg_max = max(ti_beg, ti_beg_max);
}
}
@@ -1176,6 +1213,9 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
ti_end_min = min(ti_current + ti_new_step, ti_end_min);
ti_end_max = max(ti_current + ti_new_step, ti_end_max);
+ /* What is the next starting point for this cell ? */
+ ti_beg_max = max(ti_current, ti_beg_max);
+
} else { /* gpart is inactive */
const integertime_t ti_end =
@@ -1184,6 +1224,12 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
/* What is the next sync-point ? */
ti_end_min = min(ti_end, ti_end_min);
ti_end_max = max(ti_end, ti_end_max);
+
+ const integertime_t ti_beg =
+ get_integer_time_begin(ti_current + 1, gp->time_bin);
+
+ /* What is the next starting point for this cell ? */
+ ti_beg_max = max(ti_beg, ti_beg_max);
}
}
}
@@ -1220,6 +1266,9 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
ti_end_min = min(ti_current + ti_new_step, ti_end_min);
ti_end_max = max(ti_current + ti_new_step, ti_end_max);
+ /* What is the next starting point for this cell ? */
+ ti_beg_max = max(ti_current, ti_beg_max);
+
} else { /* star particle is inactive */
const integertime_t ti_end =
@@ -1228,6 +1277,12 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
/* What is the next sync-point ? */
ti_end_min = min(ti_end, ti_end_min);
ti_end_max = max(ti_end, ti_end_max);
+
+ const integertime_t ti_beg =
+ get_integer_time_begin(ti_current + 1, sp->time_bin);
+
+ /* What is the next starting point for this cell ? */
+ ti_beg_max = max(ti_beg, ti_beg_max);
}
}
} else {
@@ -1246,6 +1301,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
s_updated += cp->s_updated;
ti_end_min = min(cp->ti_end_min, ti_end_min);
ti_end_max = max(cp->ti_end_max, ti_end_max);
+ ti_beg_max = max(cp->ti_beg_max, ti_beg_max);
}
}
@@ -1255,6 +1311,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
c->s_updated = s_updated;
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
+ c->ti_beg_max = ti_beg_max;
if (timer) TIMER_TOC(timer_timestep);
}
@@ -1433,6 +1490,11 @@ void runner_do_recv_gpart(struct runner *r, struct cell *c, int timer) {
get_integer_time_end(ti_current, gparts[k].time_bin);
ti_end_min = min(ti_end_min, ti_end);
ti_end_max = max(ti_end_max, ti_end);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (gparts[k].ti_drift != ti_current)
+ error("Received un-drifted g-particle !");
+#endif
}
}
@@ -1584,7 +1646,8 @@ void *runner_main(void *data) {
} else if (cj == NULL) { /* self */
if (!cell_is_active(ci, e) && t->type != task_type_sort &&
- t->type != task_type_send && t->type != task_type_recv)
+ t->type != task_type_send && t->type != task_type_recv &&
+ t->type != task_type_kick1)
error(
"Task (type='%s/%s') should have been skipped ti_current=%lld "
"c->ti_end_min=%lld",
@@ -1600,6 +1663,15 @@ void *runner_main(void *data) {
taskID_names[t->type], subtaskID_names[t->subtype], e->ti_current,
ci->ti_end_min, t->flags);
+ /* Special case for kick1 */
+ if (!cell_is_starting(ci, e) && t->type == task_type_kick1 &&
+ t->flags == 0)
+ error(
+ "Task (type='%s/%s') should have been skipped ti_current=%lld "
+ "c->ti_end_min=%lld t->flags=%d",
+ taskID_names[t->type], subtaskID_names[t->subtype], e->ti_current,
+ ci->ti_end_min, t->flags);
+
} else { /* pair */
if (!cell_is_active(ci, e) && !cell_is_active(cj, e))
diff --git a/src/space.c b/src/space.c
index 68d96397a3c85ed11d5e506519af699dfdd2cecf..7f0d9bcbf6341bd6074ddebf4b9d8c7cb8c04902 100644
--- a/src/space.c
+++ b/src/space.c
@@ -269,7 +269,7 @@ void space_regrid(struct space *s, int verbose) {
const size_t nr_parts = s->nr_parts;
const ticks tic = getticks();
- const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0;
+ const integertime_t ti_old = (s->e != NULL) ? s->e->ti_old : 0;
/* Run through the cells and get the current h_max. */
// tic = getticks();
@@ -437,7 +437,7 @@ void space_regrid(struct space *s, int verbose) {
c->gcount = 0;
c->scount = 0;
c->super = c;
- c->ti_old = ti_current;
+ c->ti_old = ti_old;
if (s->gravity) c->multipole = &s->multipoles_top[cid];
}
@@ -522,7 +522,7 @@ void space_rebuild(struct space *s, int verbose) {
size_t nr_gparts = s->nr_gparts;
size_t nr_sparts = s->nr_sparts;
struct cell *restrict cells_top = s->cells_top;
- const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0;
+ const integertime_t ti_old = (s->e != NULL) ? s->e->ti_old : 0;
/* Run through the particles and get their cell index. Allocates
an index that is larger than the number of particles to avoid
@@ -918,7 +918,7 @@ void space_rebuild(struct space *s, int verbose) {
struct spart *sfinger = s->sparts;
for (int k = 0; k < s->nr_cells; k++) {
struct cell *restrict c = &cells_top[k];
- c->ti_old = ti_current;
+ c->ti_old = ti_old;
c->parts = finger;
c->xparts = xfinger;
c->gparts = gfinger;
@@ -1958,7 +1958,7 @@ void space_split_recursive(struct space *s, struct cell *c,
const int depth = c->depth;
int maxdepth = 0;
float h_max = 0.0f;
- integertime_t ti_end_min = max_nr_timesteps, ti_end_max = 0;
+ integertime_t ti_end_min = max_nr_timesteps, ti_end_max = 0, ti_beg_max = 0;
struct part *parts = c->parts;
struct gpart *gparts = c->gparts;
struct spart *sparts = c->sparts;
@@ -2066,6 +2066,7 @@ void space_split_recursive(struct space *s, struct cell *c,
h_max = max(h_max, c->progeny[k]->h_max);
ti_end_min = min(ti_end_min, c->progeny[k]->ti_end_min);
ti_end_max = max(ti_end_max, c->progeny[k]->ti_end_max);
+ ti_beg_max = max(ti_beg_max, c->progeny[k]->ti_beg_max);
if (c->progeny[k]->maxdepth > maxdepth)
maxdepth = c->progeny[k]->maxdepth;
}
@@ -2121,22 +2122,28 @@ void space_split_recursive(struct space *s, struct cell *c,
const float h = p->h;
const integertime_t ti_end =
get_integer_time_end(e->ti_current, p->time_bin);
+ const integertime_t ti_beg =
+ get_integer_time_begin(e->ti_current + 1, p->time_bin);
xp->x_diff[0] = 0.f;
xp->x_diff[1] = 0.f;
xp->x_diff[2] = 0.f;
if (h > h_max) h_max = h;
if (ti_end < ti_end_min) ti_end_min = ti_end;
if (ti_end > ti_end_max) ti_end_max = ti_end;
+ if (ti_beg > ti_beg_max) ti_beg_max = ti_beg;
}
for (int k = 0; k < gcount; k++) {
struct gpart *gp = &gparts[k];
const integertime_t ti_end =
get_integer_time_end(e->ti_current, gp->time_bin);
+ const integertime_t ti_beg =
+ get_integer_time_begin(e->ti_current + 1, gp->time_bin);
gp->x_diff[0] = 0.f;
gp->x_diff[1] = 0.f;
gp->x_diff[2] = 0.f;
if (ti_end < ti_end_min) ti_end_min = ti_end;
if (ti_end > ti_end_max) ti_end_max = ti_end;
+ if (ti_beg > ti_beg_max) ti_beg_max = ti_beg;
}
for (int k = 0; k < scount; k++) {
struct spart *sp = &sparts[k];
@@ -2154,6 +2161,7 @@ void space_split_recursive(struct space *s, struct cell *c,
c->h_max = h_max;
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
+ c->ti_beg_max = ti_beg_max;
c->maxdepth = maxdepth;
/* Set ownership according to the start of the parts array. */
@@ -2436,6 +2444,11 @@ void space_init_gparts(struct space *s) {
#endif
gravity_first_init_gpart(&gp[i]);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gp->ti_drift = 0;
+ gp->ti_kick = 0;
+#endif
}
}
@@ -2462,6 +2475,11 @@ void space_init_sparts(struct space *s) {
#endif
star_first_init_spart(&sp[i]);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ sp->ti_drift = 0;
+ sp->ti_kick = 0;
+#endif
}
}
@@ -2794,17 +2812,17 @@ void space_link_cleanup(struct space *s) {
}
/**
- * @brief Checks that all cells have been drifted to the current point in time
+ * @brief Checks that all cells have been drifted to a given point in time
*
* Expensive function. Should only be used for debugging purposes.
*
* @param s The #space to check.
- * @param ti_current The (integer) time.
+ * @param ti_drift The (integer) time.
*/
-void space_check_drift_point(struct space *s, integertime_t ti_current) {
+void space_check_drift_point(struct space *s, integertime_t ti_drift) {
/* Recursively check all cells */
- space_map_cells_pre(s, 1, cell_check_drift_point, &ti_current);
+ space_map_cells_pre(s, 1, cell_check_drift_point, &ti_drift);
}
/**
diff --git a/src/space.h b/src/space.h
index cb16d7df56d4fb4c9f0206de7abb8216aacc1de9..4bbb8f32353c795772a461061096a366d0848f7e 100644
--- a/src/space.h
+++ b/src/space.h
@@ -212,7 +212,7 @@ void space_init_parts(struct space *s);
void space_init_gparts(struct space *s);
void space_init_sparts(struct space *s);
void space_link_cleanup(struct space *s);
-void space_check_drift_point(struct space *s, integertime_t ti_current);
+void space_check_drift_point(struct space *s, integertime_t ti_drift);
void space_check_timesteps(struct space *s);
void space_replicate(struct space *s, int replicate, int verbose);
void space_clean(struct space *s);
diff --git a/src/stars/Default/star_part.h b/src/stars/Default/star_part.h
index e958e3d68bc58855a4f57f24d876cfaf73362bd6..68dd4869c257e35b3be7dc21f36e6dcdb725dc17 100644
--- a/src/stars/Default/star_part.h
+++ b/src/stars/Default/star_part.h
@@ -47,6 +47,16 @@ struct spart {
/*! Particle time bin */
timebin_t time_bin;
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Time of the last drift */
+ integertime_t ti_drift;
+
+ /* Time of the last kick */
+ integertime_t ti_kick;
+
+#endif
+
} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_DEFAULT_STAR_PART_H */