diff --git a/src/engine.c b/src/engine.c
index 827778100fcdea0a89afb210677bf26fc9acd8b8..14dd0f068518264fe03d348aae5b019941e67353 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1639,8 +1639,10 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
/**
* @brief Constructs the top-level tasks for the short-range gravity
- * interactions.
+ * and long-range gravity interactions.
*
+ * - One FTT task per MPI rank.
+ * - Multiple gravity ghosts for dependencies.
* - All top-cells get a self task.
* - All pairs within range according to the multipole acceptance
* criterion get a pair task.
@@ -1653,6 +1655,7 @@ void engine_make_self_gravity_tasks(struct engine *e) {
struct scheduler *sched = &e->sched;
const int nodeID = e->nodeID;
const int periodic = s->periodic;
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
const int cdim_ghost[3] = {s->cdim[0] / 4 + 1, s->cdim[1] / 4 + 1,
s->cdim[2] / 4 + 1};
@@ -1704,7 +1707,7 @@ void engine_make_self_gravity_tasks(struct engine *e) {
scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci,
NULL);
- /* Deal with periodicity dependencies */
+ /* Deal with dependencies of the FFT mesh calculation */
const int ghost_id = cell_getid(cdim_ghost, i / 4, j / 4, k / 4);
if (ghost_id > n_ghosts) error("Invalid ghost_id");
if (periodic) {
@@ -1731,8 +1734,9 @@ void engine_make_self_gravity_tasks(struct engine *e) {
if (cj->nodeID != nodeID) continue; // MATTHIEU
/* Are the cells to close for a MM interaction ? */
- if (!gravity_multipole_accept(ci->multipole, cj->multipole,
- theta_crit_inv, 1)) {
+ if (!gravity_multipole_accept_rebuild(
+ ci->multipole, cj->multipole, theta_crit_inv, periodic,
+ dim)) {
scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0,
0, ci, cj);
@@ -1746,6 +1750,11 @@ void engine_make_self_gravity_tasks(struct engine *e) {
if (periodic) free(ghosts);
}
+/**
+ * @brief Constructs the top-level tasks for the external gravity.
+ *
+ * @param e The #engine.
+ */
void engine_make_external_gravity_tasks(struct engine *e) {
struct space *s = e->s;
diff --git a/src/multipole.h b/src/multipole.h
index 23f5194a30b7316aac15073cba36dc404efa21c1..8e8ecae54e5dedde1a789d2a5469b6e39d853b18 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -2636,29 +2636,75 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
/**
* @brief Checks whether a cell-cell interaction can be appromixated by a M-M
- * interaction.
+ * interaction using the CoM and cell radius at rebuild.
+ *
+ * We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179,
+ * Issue 1, pp.27-42, equation 10.
+ *
+ * @param ma The #multipole of the first #cell.
+ * @param mb The #multipole of the second #cell.
+ * @param theta_crit_inv The inverse of the critical opening angle.
+ * @param periodic Are we using periodic boundary conditions ?
+ * @param dim The dimensions of the box.
+ */
+__attribute__((always_inline)) INLINE static int
+gravity_multipole_accept_rebuild(const struct gravity_tensors *ma,
+ const struct gravity_tensors *mb,
+ double theta_crit_inv, int periodic,
+ const double dim[3]) {
+
+ const double r_crit_a = ma->r_max_rebuild * theta_crit_inv;
+ const double r_crit_b = mb->r_max_rebuild * theta_crit_inv;
+
+ double dx = ma->CoM_rebuild[0] - mb->CoM_rebuild[0];
+ double dy = ma->CoM_rebuild[1] - mb->CoM_rebuild[1];
+ double dz = ma->CoM_rebuild[2] - mb->CoM_rebuild[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;
+
+ // MATTHIEU: Make this mass-dependent ?
+
+ /* Multipole acceptance criterion (Dehnen 2002, eq.10) */
+ return (r2 > (r_crit_a + r_crit_b) * (r_crit_a + r_crit_b));
+}
+
+/**
+ * @brief Checks whether a cell-cell interaction can be appromixated by a M-M
+ * interaction using the CoM and cell radius at the current time.
+ *
+ * We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179,
+ * Issue 1, pp.27-42, equation 10.
*
* @param ma The #multipole of the first #cell.
* @param mb The #multipole of the second #cell.
* @param theta_crit_inv The inverse of the critical opening angle.
- * @param rebuild Are we using the current value of CoM or the ones from
- * the last rebuild ?
+ * @param periodic Are we using periodic boundary conditions ?
+ * @param dim The dimensions of the box.
*/
__attribute__((always_inline)) INLINE static int gravity_multipole_accept(
const struct gravity_tensors *ma, const struct gravity_tensors *mb,
- double theta_crit_inv, int rebuild) {
-
- const double r_crit_a =
- (rebuild ? ma->r_max_rebuild : ma->r_max) * theta_crit_inv;
- const double r_crit_b =
- (rebuild ? mb->r_max_rebuild : mb->r_max) * theta_crit_inv;
-
- const double dx = rebuild ? ma->CoM_rebuild[0] - mb->CoM_rebuild[0]
- : ma->CoM[0] - mb->CoM[0];
- const double dy = rebuild ? ma->CoM_rebuild[1] - mb->CoM_rebuild[1]
- : ma->CoM[1] - mb->CoM[1];
- const double dz = rebuild ? ma->CoM_rebuild[2] - mb->CoM_rebuild[2]
- : ma->CoM[2] - mb->CoM[2];
+ double theta_crit_inv, int periodic, const double dim[3]) {
+
+ const double r_crit_a = ma->r_max * theta_crit_inv;
+ const double r_crit_b = mb->r_max * theta_crit_inv;
+
+ double dx = ma->CoM[0] - mb->CoM[0];
+ double dy = ma->CoM[1] - mb->CoM[1];
+ double dz = ma->CoM[2] - mb->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;
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index a66cc5e0c9ed241aba3bb1b4329016b8e505e280..0a5915cd10a170171378d708f3f0cda54c34e9b0 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -415,6 +415,9 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
/* Some constants */
const struct engine *e = r->e;
+ const struct space *s = e->s;
+ const int periodic = s->periodic;
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const struct gravity_props *props = e->gravity_properties;
const double theta_crit_inv = props->theta_crit_inv;
@@ -464,7 +467,8 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
/* Can we use M-M interactions ? */
if (gravity_multipole_accept(ci->multipole, cj->multipole, theta_crit_inv,
- 0)) {
+ periodic, dim)) {
+
/* MATTHIEU: make a symmetric M-M interaction function ! */
runner_dopair_grav_mm(r, ci, cj);
runner_dopair_grav_mm(r, cj, ci);
@@ -617,6 +621,9 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
/* Some constants */
const struct engine *e = r->e;
+ const struct space *s = e->s;
+ const int periodic = s->periodic;
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const struct gravity_props *props = e->gravity_properties;
const double theta_crit_inv = props->theta_crit_inv;
@@ -627,7 +634,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
const int nr_cells = e->s->nr_cells;
/* Anything to do here? */
- if (!cell_is_active(ci, e)) return; // MATTHIEU (should never happen)
+ if (!cell_is_active(ci, e)) return;
/* Check multipole has been drifted */
if (ci->ti_old_multipole != e->ti_current)
@@ -645,14 +652,14 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
/* Check the multipole acceptance criterion */
if (gravity_multipole_accept(ci->multipole, cj->multipole, theta_crit_inv,
- 0)) {
+ periodic, dim)) {
/* Go for a (non-symmetric) M-M calculation */
runner_dopair_grav_mm(r, ci, cj);
}
/* Is the criterion violated now but was OK at the last rebuild ? */
- else if (gravity_multipole_accept(ci->multipole, cj->multipole,
- theta_crit_inv, 1)) {
+ else if (gravity_multipole_accept_rebuild(ci->multipole, cj->multipole,
+ theta_crit_inv, periodic, dim)) {
/* Alright, we have to take charge of that pair in a different way. */
// MATTHIEU: We should actually open the tree-node here and recurse.
diff --git a/src/scheduler.c b/src/scheduler.c
index b07c403e4ecd960b22b51f24372ca0a3420a453f..76db1f2a279828800862629340598857c99e306f 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -682,7 +682,7 @@ static void scheduler_splittask_gravity(struct task *t, struct scheduler *s) {
/* convert to a self-subtask. */
t->type = task_type_sub_self;
- /* Make sure we have a drift task (MATTHIEU temp. fix) */
+ /* Make sure we have a drift task */
lock_lock(&ci->lock);
if (ci->drift_gpart == NULL)
ci->drift_gpart = scheduler_addtask(