diff --git a/src/cell.c b/src/cell.c
index 2c4afe81a060103bb3dc239e23233c75c0458ec0..799c3333567653376412f04a6867f2e13a0599a3 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -2028,20 +2028,31 @@ void cell_activate_grav_mm_task(struct cell *ci, struct cell *cj,
struct scheduler *s) {
/* Some constants */
const struct engine *e = s->space->e;
+ const integertime_t ti_current = e->ti_current;
/* Anything to do here? */
if (!cell_is_active_gravity_mm(ci, e) && !cell_is_active_gravity_mm(cj, e))
error("Inactive MM task being activated");
- /* Atomically drift the multipole in ci */
- lock_lock(&ci->mlock);
- if (ci->ti_old_multipole < e->ti_current) cell_drift_multipole(ci, e);
- if (lock_unlock(&ci->mlock) != 0) error("Impossible to unlock m-pole");
+ /* Atomically drift the multipoles in the progenies of ci */
+ for (int i = 0; i < 8; i++) {
+ struct cell *cpi = ci->progeny[i];
+ if (cpi != NULL) {
+ lock_lock(&cpi->mlock);
+ if (cpi->ti_old_multipole < ti_current) cell_drift_multipole(cpi, e);
+ if (lock_unlock(&cpi->mlock) != 0) error("Impossible to unlock m-pole");
+ }
+ }
- /* Atomically drift the multipole in cj */
- lock_lock(&cj->mlock);
- if (cj->ti_old_multipole < e->ti_current) cell_drift_multipole(cj, e);
- if (lock_unlock(&cj->mlock) != 0) error("Impossible to unlock m-pole");
+ /* Atomically drift the multipoles in the progenies of cj */
+ for (int j = 0; j < 8; j++) {
+ struct cell *cpj = cj->progeny[j];
+ if (cpj != NULL) {
+ lock_lock(&cpj->mlock);
+ if (cpj->ti_old_multipole < ti_current) cell_drift_multipole(cpj, e);
+ if (lock_unlock(&cpj->mlock) != 0) error("Impossible to unlock m-pole");
+ }
+ }
}
/**
@@ -2490,7 +2501,7 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
struct cell *ci = t->ci;
struct cell *cj = t->cj;
const int ci_active = cell_is_active_gravity_mm(ci, e);
- const int cj_active = (cj != NULL) ? cell_is_active_gravity_mm(cj, e) : 0;
+ const int cj_active = cell_is_active_gravity_mm(cj, e);
#ifdef WITH_MPI
const int ci_nodeID = ci->nodeID;
const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
diff --git a/src/engine.c b/src/engine.c
index 0b15cc9765b6b59ca82b2bd2c7a8d6c1e9ceebed..148d915cad263fcee1fa5bb12386c61343d1f17d 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2674,9 +2674,29 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
}
#endif
- /* Note that we do not need to link the M-M tasks */
- /* since we already did so when splitting the gravity */
- /* tasks. */
+ /* Multipole-multipole interaction of progenies */
+ } else if (t_type == task_type_grav_mm) {
+
+ atomic_inc(&ci->nr_mm_tasks);
+ atomic_inc(&cj->nr_mm_tasks);
+ engine_addlink(e, &ci->grav_mm, t);
+ engine_addlink(e, &cj->grav_mm, t);
+
+ /* for (int i = 0; i < 8; i++) { */
+ /* struct cell *cpi = ci->progeny[i]; */
+ /* if (cpi != NULL) { */
+ /* atomic_inc(&cpi->nr_mm_tasks); */
+ /* engine_addlink(e, &cpi->grav_mm, t); */
+ /* } */
+ /* } */
+
+ /* for (int j = 0; j < 8; j++) { */
+ /* struct cell *cpj = cj->progeny[j]; */
+ /* if (cpj != NULL) { */
+ /* atomic_inc(&cpj->nr_mm_tasks); */
+ /* engine_addlink(e, &cpj->grav_mm, t); */
+ /* } */
+ /* } */
}
}
}
diff --git a/src/runner.c b/src/runner.c
index 0be6288cbcc0b77fbe9ba3128d820b82a8943480..fc22bd8014edc49612a50fc4ee84d948eaff2205 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -2273,7 +2273,7 @@ void *runner_main(void *data) {
runner_do_grav_long_range(r, t->ci, 1);
break;
case task_type_grav_mm:
- runner_dopair_grav_mm(r, t->ci, t->cj);
+ runner_dopair_grav_mm_progenies(r, t->ci, t->cj);
break;
case task_type_cooling:
runner_do_cooling(r, t->ci, 1);
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index f9d1f833e65f7e6dde6085696dfbcd0719ae2a87..2a9e78a3ede3a16994eb79ddd6bec47ad032e866 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -1300,8 +1300,35 @@ static INLINE void runner_dopair_grav_mm(struct runner *r,
runner_dopair_grav_mm_nonsym(r, ci, cj);
else if (do_j)
runner_dopair_grav_mm_nonsym(r, cj, ci);
- else
- error("Running M-M task with two inactive cells.");
+ /* else */
+ /* error("Running M-M task with two inactive cells."); */
+}
+
+static INLINE void runner_dopair_grav_mm_progenies(struct runner *r,
+ struct cell *restrict ci,
+ struct cell *restrict cj) {
+
+ const struct engine *e = r->e;
+ const struct space *s = e->s;
+
+ /* Loop over all pairs of progenies */
+ for (int i = 0; i < 8; i++) {
+ if (ci->progeny[i] != NULL) {
+ for (int j = 0; j < 8; j++) {
+ if (cj->progeny[j] != NULL) {
+
+ struct cell *cpi = ci->progeny[i];
+ struct cell *cpj = cj->progeny[j];
+
+ /* Did we agree to use an M-M interaction here at the last rebuild? */
+ if (cell_can_use_pair_mm_rebuild(cpi, cpj, e, s)) {
+
+ runner_dopair_grav_mm(r, cpi, cpj);
+ }
+ }
+ }
+ }
+ }
}
static INLINE void runner_dopair_recursive_grav_pm(struct runner *r,
@@ -1631,13 +1658,13 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
if (ci->ti_old_multipole != e->ti_current)
error("Interacting un-drifted multipole");
+ /* Get this cell's multipole information */
+ struct gravity_tensors *const multi_i = ci->multipole;
+
/* Find this cell's top-level (great-)parent */
struct cell *top = ci;
while (top->parent != NULL) top = top->parent;
- /* Flag that contributions will be recieved */
- struct gravity_tensors *const multi_i = ci->multipole;
-
/* Recover the top-level multipole (for distance checks) */
struct gravity_tensors *const multi_top = top->multipole;
const double CoM_rebuild_top[3] = {multi_top->CoM_rebuild[0],
diff --git a/src/scheduler.c b/src/scheduler.c
index abfa6c549818a0e3c2b85a3cdd65e199557cf6e4..a5da777f45f8075c7c577988eaa5dfbd7d1ca5ed 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -888,20 +888,6 @@ static void scheduler_splittask_gravity(struct task *t, struct scheduler *s) {
break;
}
- /* Should we replace it with an M-M task? */
- if (cell_can_use_pair_mm_rebuild(ci, cj, e, sp)) {
-
- t->type = task_type_grav_mm;
- t->subtype = task_subtype_none;
-
- /* Since this task will not be split, we can already link it */
- atomic_inc(&ci->nr_mm_tasks);
- atomic_inc(&cj->nr_mm_tasks);
- engine_addlink(e, &ci->grav_mm, t);
- engine_addlink(e, &cj->grav_mm, t);
- break;
- }
-
/* Should this task be split-up? */
if (cell_can_split_pair_gravity_task(ci) &&
cell_can_split_pair_gravity_task(cj)) {
@@ -910,37 +896,33 @@ static void scheduler_splittask_gravity(struct task *t, struct scheduler *s) {
const long long gcount_j = cj->gcount;
/* Replace by a single sub-task? */
- if (scheduler_dosub && /* Use division to avoid integer overflow. */
+ if (scheduler_dosub &&
gcount_i * gcount_j < ((long long)space_subsize_pair_grav)) {
/* Otherwise, split it. */
} else {
- /* Take a step back (we're going to recycle the current task)... */
- redo = 1;
-
- /* Find the first non-empty childrens of the cells */
- int first_ci_child = 0, first_cj_child = 0;
- while (ci->progeny[first_ci_child] == NULL) first_ci_child++;
- while (cj->progeny[first_cj_child] == NULL) first_cj_child++;
-
- /* Recycle the current pair */
- t->ci = ci->progeny[first_ci_child];
- t->cj = cj->progeny[first_cj_child];
+ /* Turn the task into a M-M task that will take care of all the
+ * progeny pairs */
+ t->type = task_type_grav_mm;
+ t->subtype = task_subtype_none;
/* Make a task for every other pair of progeny */
- for (int i = first_ci_child; i < 8; i++) {
+ for (int i = 0; i < 8; i++) {
if (ci->progeny[i] != NULL) {
- for (int j = first_cj_child; j < 8; j++) {
+ for (int j = 0; j < 8; j++) {
if (cj->progeny[j] != NULL) {
- /* Skip the recycled pair */
- if (i == first_ci_child && j == first_cj_child) continue;
+ /* But only for pairs that cannot be replaced by a M-M
+ * interaction */
+ if (!cell_can_use_pair_mm_rebuild(ci->progeny[i],
+ cj->progeny[j], e, sp)) {
- scheduler_splittask_gravity(
- scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
- ci->progeny[i], cj->progeny[j]),
- s);
+ scheduler_splittask_gravity(
+ scheduler_addtask(s, task_type_pair, task_subtype_grav,
+ 0, 0, ci->progeny[i], cj->progeny[j]),
+ s);
+ }
}
}
}