diff --git a/examples/main.c b/examples/main.c
index ee1253062409ec2e787e064a5fb50da2c830d35d..672683f13c3c53e1b456f6eb65330dcf38233f2e 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -504,6 +504,8 @@ int main(int argc, char *argv[]) {
fflush(stdout);
}
+ periodic = 0;
+
#ifdef SWIFT_DEBUG_CHECKS
/* Check once and for all that we don't have unwanted links */
if (!with_stars) {
diff --git a/src/cell.c b/src/cell.c
index f85e4dbd2a71033b10852bc65a97b8eec6bf75f9..e619ddf246d46d4258dd66e6f9301eaf2dca7b94 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1290,41 +1290,14 @@ void cell_clean(struct cell *c) {
if (c->progeny[k]) cell_clean(c->progeny[k]);
}
-/**
- * @brief Checks whether a given cell needs drifting or not.
- *
- * @param c the #cell.
- * @param e The #engine (holding current time information).
- *
- * @return 1 If the cell needs drifting, 0 otherwise.
- */
-int cell_is_drift_needed(struct cell *c, const struct engine *e) {
-
- /* Do we have at least one active particle in the cell ?*/
- if (cell_is_active(c, e)) return 1;
-
- /* Loop over the pair tasks that involve this cell */
- for (struct link *l = c->density; l != NULL; l = l->next) {
-
- if (l->t->type != task_type_pair && l->t->type != task_type_sub_pair)
- continue;
-
- /* Is the other cell in the pair active ? */
- if ((l->t->ci == c && cell_is_active(l->t->cj, e)) ||
- (l->t->cj == c && cell_is_active(l->t->ci, e)))
- return 1;
- }
-
- /* No neighbouring cell has active particles. Drift not necessary */
- return 0;
-}
-
/**
* @brief Clear the drift flags on the given cell.
*/
void cell_clear_drift_flags(struct cell *c, void *data) {
c->do_drift = 0;
c->do_sub_drift = 0;
+ c->do_grav_drift = 0;
+ c->do_grav_sub_drift = 0;
}
/**
@@ -1359,13 +1332,32 @@ void cell_activate_drift_part(struct cell *c, struct scheduler *s) {
*/
void cell_activate_drift_gpart(struct cell *c, struct scheduler *s) {
- scheduler_activate(s, c->super->drift_gpart);
+ /* If this cell is already marked for drift, quit early. */
+ if (c->do_grav_drift) return;
+
+ /* Mark this cell for drifting. */
+ c->do_grav_drift = 1;
+
+ /* Set the do_grav_sub_drifts all the way up and activate the super drift
+ if this has not yet been done. */
+ if (c == c->super) {
+ scheduler_activate(s, c->drift_gpart);
+ } else {
+ for (struct cell *parent = c->parent;
+ parent != NULL && !parent->do_grav_sub_drift;
+ parent = parent->parent) {
+ parent->do_grav_sub_drift = 1;
+ if (parent == c->super) {
+ scheduler_activate(s, parent->drift_gpart);
+ break;
+ }
+ }
+ }
}
/**
* @brief Activate the sorts up a cell hierarchy.
*/
-
void cell_activate_sorts_up(struct cell *c, struct scheduler *s) {
if (c == c->super) {
scheduler_activate(s, c->sorts);
@@ -1409,7 +1401,8 @@ void cell_activate_sorts(struct cell *c, int sid, struct scheduler *s) {
}
/**
- * @brief Traverse a sub-cell task and activate the sort tasks along the way.
+ * @brief Traverse a sub-cell task and activate the drift and sort tasks along
+ * the way.
*/
void cell_activate_subcell_tasks(struct cell *ci, struct cell *cj,
struct scheduler *s) {
@@ -1676,6 +1669,132 @@ void cell_activate_subcell_tasks(struct cell *ci, struct cell *cj,
}
}
+/**
+ * @brief Traverse a sub-cell task and activate the gravity drift tasks
+ */
+void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
+ struct scheduler *s) {
+ /* Some constants */
+ const struct space *sp = s->space;
+ const struct engine *e = sp->e;
+ const int periodic = sp->periodic;
+ const double dim[3] = {sp->dim[0], sp->dim[1], sp->dim[2]};
+ const double theta_crit = e->gravity_properties->theta_crit;
+
+ /* Self interaction? */
+ if (cj == NULL) {
+
+ /* Do anything? */
+ if (!cell_is_active(ci, e)) return;
+
+ /* Recurse? */
+ if (ci->split) {
+
+ /* Loop over all progenies and pairs of progenies */
+ for (int j = 0; j < 8; j++) {
+ if (ci->progeny[j] != NULL) {
+ cell_activate_subcell_grav_tasks(ci->progeny[j], NULL, s);
+ for (int k = j + 1; k < 8; k++)
+ if (ci->progeny[k] != NULL)
+ cell_activate_subcell_grav_tasks(ci->progeny[j], ci->progeny[k],
+ s);
+ }
+ }
+ } else {
+
+ /* We have reached the bottom of the tree: activate gpart drift */
+ cell_activate_drift_gpart(ci, s);
+ }
+ }
+
+ /* Pair interaction */
+ else {
+
+ /* Anything to do here? */
+ if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
+
+ /* Recover the multipole information */
+ struct gravity_tensors *const multi_i = ci->multipole;
+ struct gravity_tensors *const multi_j = cj->multipole;
+ const double ri_max = multi_i->r_max;
+ const double rj_max = multi_j->r_max;
+
+ /* Get the distance between the CoMs */
+ double dx = multi_i->CoM[0] - multi_j->CoM[0];
+ double dy = multi_i->CoM[1] - multi_j->CoM[1];
+ double dz = multi_i->CoM[2] - multi_j->CoM[2];
+
+ /* Apply BC */
+ if (periodic) {
+ dx = nearest(dx, dim[0]);
+ dy = nearest(dy, dim[1]);
+ dz = nearest(dz, dim[2]);
+ }
+ const double r2 = dx * dx + dy * dy + dz * dz;
+
+ /* Can we use multipoles ? */
+ if (gravity_multipole_accept(multi_i, multi_j, theta_crit, r2)) {
+
+ /* Ok, no need to drift anything */
+ return;
+ }
+ /* Otherwise, activate the gpart drifts if we are at the bottom. */
+ else if (!ci->split && !cj->split) {
+
+ /* Activate the drifts if the cells are local. */
+ if (cell_is_active(ci, e) || cell_is_active(cj, e)) {
+ if (ci->nodeID == engine_rank) cell_activate_drift_gpart(ci, s);
+ if (cj->nodeID == engine_rank) cell_activate_drift_gpart(cj, s);
+ }
+ }
+ /* Ok, we can still recurse */
+ else {
+
+ if (ri_max > rj_max) {
+ if (ci->split) {
+
+ /* Loop over ci's children */
+ for (int k = 0; k < 8; k++) {
+ if (ci->progeny[k] != NULL)
+ cell_activate_subcell_grav_tasks(ci->progeny[k], cj, s);
+ }
+
+ } else if (cj->split) {
+
+ /* Loop over cj's children */
+ for (int k = 0; k < 8; k++) {
+ if (cj->progeny[k] != NULL)
+ cell_activate_subcell_grav_tasks(ci, cj->progeny[k], s);
+ }
+
+ } else {
+ error("Fundamental error in the logic");
+ }
+ } else if (rj_max >= ri_max) {
+ if (cj->split) {
+
+ /* Loop over cj's children */
+ for (int k = 0; k < 8; k++) {
+ if (cj->progeny[k] != NULL)
+ cell_activate_subcell_grav_tasks(ci, cj->progeny[k], s);
+ }
+
+ } else if (ci->split) {
+
+ /* Loop over ci's children */
+ for (int k = 0; k < 8; k++) {
+ if (ci->progeny[k] != NULL)
+ cell_activate_subcell_grav_tasks(ci->progeny[k], cj, s);
+ }
+
+ } else {
+ error("Fundamental error in the logic");
+ }
+ }
+ }
+ }
+}
+
/**
* @brief Un-skips all the tasks associated with a given cell and checks
* if the space needs to be rebuilt.
@@ -1701,8 +1820,13 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
(cj != NULL && cell_is_active(cj, e) && cj->nodeID == engine_rank)) {
scheduler_activate(s, t);
- /* Set the correct sorting flags */
- if (t->type == task_type_pair) {
+ /* Activate hydro drift */
+ if (t->type == task_type_self) {
+ if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s);
+ }
+
+ /* Set the correct sorting flags and activate hydro drifts */
+ else if (t->type == task_type_pair) {
/* Store some values. */
atomic_or(&ci->requires_sorts, 1 << t->flags);
atomic_or(&cj->requires_sorts, 1 << t->flags);
@@ -1863,9 +1987,10 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate(s, t);
/* Set the drifting flags */
- if (t->type == task_type_pair) {
- cell_activate_drift_gpart(ci, s);
- cell_activate_drift_gpart(cj, s);
+ if (t->type == task_type_self) {
+ cell_activate_subcell_grav_tasks(t->ci, NULL, s);
+ } else if (t->type == task_type_pair) {
+ cell_activate_subcell_grav_tasks(t->ci, t->cj, s);
}
}
}
@@ -1955,7 +2080,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
/* Check that we are actually going to move forward. */
if (ti_current < ti_old_part) error("Attempt to drift to the past");
-#endif // SWIFT_DEBUG_CHECKS
+#endif
/* Are we not in a leaf ? */
if (c->split && (force || c->do_sub_drift)) {
@@ -2040,8 +2165,9 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
*
* @param c The #cell.
* @param e The #engine (to get ti_current).
+ * @param force Drift the particles irrespective of the #cell flags.
*/
-void cell_drift_gpart(struct cell *c, const struct engine *e) {
+void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
const double timeBase = e->timeBase;
const integertime_t ti_old_gpart = c->ti_old_gpart;
@@ -2053,11 +2179,19 @@ void cell_drift_gpart(struct cell *c, const struct engine *e) {
const double dt = (ti_current - ti_old_gpart) * timeBase;
float dx_max = 0.f, dx2_max = 0.f;
+ /* Drift irrespective of cell flags? */
+ force |= c->do_grav_drift;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that we only drift local cells. */
+ if (c->nodeID != engine_rank) error("Drifting a foreign cell is nope.");
+
/* Check that we are actually going to move forward. */
if (ti_current < ti_old_gpart) error("Attempt to drift to the past");
+#endif
/* Are we not in a leaf ? */
- if (c->split) {
+ if (c->split && (force || c->do_grav_sub_drift)) {
/* Loop over the progeny and collect their data. */
for (int k = 0; k < 8; k++)
@@ -2065,13 +2199,19 @@ void cell_drift_gpart(struct cell *c, const struct engine *e) {
struct cell *cp = c->progeny[k];
/* Recurse */
- cell_drift_gpart(cp, e);
+ cell_drift_gpart(cp, e, force);
/* Update */
dx_max = max(dx_max, cp->dx_max_gpart);
}
- } else if (ti_current > ti_old_gpart) {
+ /* Store the values */
+ c->dx_max_gpart = dx_max;
+
+ /* Update the time of the last drift */
+ c->ti_old_gpart = ti_current;
+
+ } else if (!c->split && force && ti_current > ti_old_gpart) {
/* Loop over all the g-particles in the cell */
const size_t nr_gparts = c->gcount;
@@ -2111,16 +2251,16 @@ void cell_drift_gpart(struct cell *c, const struct engine *e) {
/* Now, get the maximal particle motion from its square */
dx_max = sqrtf(dx2_max);
- } else {
+ /* Store the values */
+ c->dx_max_gpart = dx_max;
- dx_max = c->dx_max_gpart;
+ /* Update the time of the last drift */
+ c->ti_old_gpart = ti_current;
}
- /* Store the values */
- c->dx_max_gpart = dx_max;
-
- /* Update the time of the last drift */
- c->ti_old_gpart = ti_current;
+ /* Clear the drift flags. */
+ c->do_grav_drift = 0;
+ c->do_grav_sub_drift = 0;
}
/**
diff --git a/src/cell.h b/src/cell.h
index 82cda87444e60fde9eb4b9a6914c2ed1ba3a470d..657ff737ed8e1e53c2addb88a7e80961a11dbd8e 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -331,12 +331,18 @@ struct cell {
/* Bit mask of sort directions that will be needed in the next timestep. */
unsigned int requires_sorts;
- /*! Does this cell need to be drifted? */
+ /*! Does this cell need to be drifted (hydro)? */
char do_drift;
- /*! Do any of this cell's sub-cells need to be drifted? */
+ /*! Do any of this cell's sub-cells need to be drifted (hydro)? */
char do_sub_drift;
+ /*! Does this cell need to be drifted (gravity)? */
+ char do_grav_drift;
+
+ /*! Do any of this cell's sub-cells need to be drifted (gravity)? */
+ char do_grav_sub_drift;
+
/*! Bit mask of sorts that need to be computed for this cell. */
unsigned int do_sort;
@@ -390,17 +396,18 @@ void cell_check_part_drift_point(struct cell *c, void *data);
void cell_check_gpart_drift_point(struct cell *c, void *data);
void cell_check_multipole_drift_point(struct cell *c, void *data);
void cell_reset_task_counters(struct cell *c);
-int cell_is_drift_needed(struct cell *c, const struct engine *e);
int cell_unskip_tasks(struct cell *c, struct scheduler *s);
void cell_set_super(struct cell *c, struct cell *super);
void cell_drift_part(struct cell *c, const struct engine *e, int force);
-void cell_drift_gpart(struct cell *c, const struct engine *e);
+void cell_drift_gpart(struct cell *c, const struct engine *e, int force);
void cell_drift_multipole(struct cell *c, const struct engine *e);
void cell_drift_all_multipoles(struct cell *c, const struct engine *e);
void cell_check_timesteps(struct cell *c);
void cell_store_pre_drift_values(struct cell *c);
void cell_activate_subcell_tasks(struct cell *ci, struct cell *cj,
struct scheduler *s);
+void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
+ struct scheduler *s);
void cell_activate_drift_part(struct cell *c, struct scheduler *s);
void cell_activate_drift_gpart(struct cell *c, struct scheduler *s);
void cell_activate_sorts(struct cell *c, int sid, struct scheduler *s);
diff --git a/src/engine.c b/src/engine.c
index 41b9ad4fa87b438f3179a83a1523ae5c70bd493b..8a66b56007811e2a1327f9b64fddfb478b9e650f 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2653,16 +2653,32 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
struct task *t = &tasks[ind];
/* Single-cell task? */
- if (t->type == task_type_self || t->type == task_type_ghost ||
- t->type == task_type_extra_ghost || t->type == task_type_cooling ||
- t->type == task_type_sourceterms || t->type == task_type_sub_self) {
+ if (t->type == task_type_self || t->type == task_type_sub_self) {
+
+ /* Local pointer. */
+ struct cell *ci = t->ci;
+
+ if (ci->nodeID != engine_rank) error("Non-local self task found");
/* Set this task's skip. */
- if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
+ if (cell_is_active(ci, e)) scheduler_activate(s, t);
+ /* Activate the hydro drift */
+ if (t->type == task_type_self && t->subtype == task_subtype_density) {
+ cell_activate_drift_part(ci, s);
+ }
+ /* Activate the gravity drift */
+ else if (t->type == task_type_self && t->subtype == task_subtype_grav) {
+ cell_activate_subcell_grav_tasks(t->ci, NULL, s);
+ }
/* Store current values of dx_max and h_max. */
- if (t->type == task_type_sub_self && t->subtype == task_subtype_density) {
- cell_activate_subcell_tasks(t->ci, NULL, s);
+ else if (t->type == task_type_sub_self &&
+ t->subtype == task_subtype_density) {
+ cell_activate_subcell_tasks(ci, NULL, s);
+
+ } else if (t->type == task_type_sub_self &&
+ t->subtype == task_subtype_grav) {
+ error("Invalid task sub-type encountered");
}
}
@@ -2673,34 +2689,42 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
struct cell *ci = t->ci;
struct cell *cj = t->cj;
- /* If this task does not involve any active cells, skip it. */
- if (!cell_is_active(t->ci, e) && !cell_is_active(t->cj, e)) continue;
-
/* Only activate tasks that involve a local active cell. */
if ((cell_is_active(ci, e) && ci->nodeID == engine_rank) ||
- (cj != NULL && cell_is_active(cj, e) && cj->nodeID == engine_rank)) {
+ (cell_is_active(cj, e) && cj->nodeID == engine_rank)) {
scheduler_activate(s, t);
/* Set the correct sorting flags */
if (t->type == task_type_pair && t->subtype == task_subtype_density) {
+
/* Store some values. */
atomic_or(&ci->requires_sorts, 1 << t->flags);
atomic_or(&cj->requires_sorts, 1 << t->flags);
ci->dx_max_sort_old = ci->dx_max_sort;
cj->dx_max_sort_old = cj->dx_max_sort;
- /* Activate the drift tasks. */
+ /* Activate the hydro drift tasks. */
if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s);
if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s);
/* Check the sorts and activate them if needed. */
cell_activate_sorts(ci, t->flags, s);
cell_activate_sorts(cj, t->flags, s);
+
+ } else if (t->type == task_type_pair &&
+ t->subtype == task_subtype_grav) {
+ /* Activate the gravity drift */
+ cell_activate_subcell_grav_tasks(t->ci, t->cj, s);
}
+
/* Store current values of dx_max and h_max. */
else if (t->type == task_type_sub_pair &&
t->subtype == task_subtype_density) {
cell_activate_subcell_tasks(t->ci, t->cj, s);
+
+ } else if (t->type == task_type_sub_pair &&
+ t->subtype == task_subtype_grav) {
+ error("Invalid task sub-type encountered");
}
}
@@ -2833,20 +2857,24 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
}
- /* Kick/Drift/init ? */
- if (t->type == task_type_kick1 || t->type == task_type_kick2 ||
- t->type == task_type_drift_part || t->type == task_type_drift_gpart ||
- t->type == task_type_init_grav) {
+ /* Kick/init ? */
+ else if (t->type == task_type_kick1 || t->type == task_type_kick2 ||
+ t->type == task_type_init_grav) {
if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
}
- /* Gravity ? */
+ /* Hydro ghost tasks ? */
+ else if (t->type == task_type_ghost || t->type == task_type_extra_ghost) {
+ if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
+ }
+
+ /* Gravity stuff ? */
else if (t->type == task_type_grav_down ||
t->type == task_type_grav_long_range) {
if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
}
- /* Periodic gravity ? */
+ /* Periodic gravity stuff (Note this is not linked to a cell) ? */
else if (t->type == task_type_grav_top_level ||
t->type == task_type_grav_ghost) {
scheduler_activate(s, t);
@@ -2859,6 +2887,11 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
t->ci->s_updated = 0;
if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
}
+
+ /* Subgrid tasks */
+ else if (t->type == task_type_cooling || t->type == task_type_sourceterms) {
+ if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
+ }
}
}
@@ -3723,7 +3756,7 @@ void engine_do_drift_all_mapper(void *map_data, int num_elements,
cell_drift_part(c, e, 1);
/* Drift all the g-particles */
- cell_drift_gpart(c, e);
+ cell_drift_gpart(c, e, 1);
/* Drift the multipoles */
if (e->policy & engine_policy_self_gravity)
diff --git a/src/multipole.h b/src/multipole.h
index 255be3e538dcf50ced1a62a5fa9fef2a93f917da..03b6d134c865ed217f1c2fe82453de9f4a3406a5 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -213,18 +213,18 @@ INLINE static void gravity_reset(struct gravity_tensors *m) {
*/
INLINE static void gravity_drift(struct gravity_tensors *m, double dt) {
- const double dx = m->m_pole.vel[0] * dt;
- const double dy = m->m_pole.vel[1] * dt;
- const double dz = m->m_pole.vel[2] * dt;
+ /* const double dx = m->m_pole.vel[0] * dt; */
+ /* const double dy = m->m_pole.vel[1] * dt; */
+ /* const double dz = m->m_pole.vel[2] * dt; */
/* Move the whole thing according to bulk motion */
- m->CoM[0] += dx;
- m->CoM[1] += dy;
- m->CoM[2] += dz;
+ /* m->CoM[0] += dx; */
+ /* m->CoM[1] += dy; */
+ /* m->CoM[2] += dz; */
/* Conservative change in maximal radius containing all gpart */
/* MATTHIEU: Use gpart->x_diff here ? */
- m->r_max += sqrt(dx * dx + dy * dy + dz * dz);
+ /* m->r_max += sqrt(dx * dx + dy * dy + dz * dz); */
}
/**
diff --git a/src/runner.c b/src/runner.c
index bb5c5cb91c7ca0ab0040f1ff637e124b1c673aeb..eb968f58440e678df4ac4df1e83ee734100b65f5 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -903,7 +903,7 @@ void runner_do_drift_gpart(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
- cell_drift_gpart(c, r->e);
+ cell_drift_gpart(c, r->e, 0);
if (timer) TIMER_TOC(timer_drift_gpart);
}