From f0728bfb5ea6b4d0b679e98d4de50f991e03715d Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Tue, 14 Jun 2016 18:40:40 +0100
Subject: [PATCH] Modified the mm task to use all the multipoles interacting
 with one given cell

---
 src/const.h              |  2 +-
 src/engine.c             | 35 +++++++++++++++++++++++++++++------
 src/runner.c             | 15 +++++----------
 src/runner_doiact_grav.h | 39 +++++++++++++++++++--------------------
 src/scheduler.c          |  3 +--
 src/task.c               | 15 +++++----------
 src/task.h               |  4 ++--
 7 files changed, 62 insertions(+), 51 deletions(-)

diff --git a/src/const.h b/src/const.h
index ba71e3aeb1..14a4ea6f68 100644
--- a/src/const.h
+++ b/src/const.h
@@ -61,6 +61,6 @@
 /* External gravity properties */
 #define EXTERNAL_POTENTIAL_POINTMASS
 
-//#define SANITY_CHECKS
+#define SANITY_CHECKS
 
 #endif /* SWIFT_CONST_H */
diff --git a/src/engine.c b/src/engine.c
index 438ae251ae..f0878af71b 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1079,7 +1079,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
 }
 
 /**
- * @brief Constructs the top-level pair tasks for the short-range gravity
+ * @brief Constructs the top-level tasks for the short-range gravity
  * interactions.
  *
  * All top-cells get a self task.
@@ -1110,6 +1110,10 @@ void engine_make_gravity_tasks(struct engine *e) {
     scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci, NULL,
                       0);
 
+    /* Let's also build a task for all the non-neighbouring pm calculations */
+    scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, 0, 0, ci,
+                      NULL, 0);
+
     for (int cjd = cid + 1; cjd < nr_cells; ++cjd) {
 
       struct cell *cj = &cells[cjd];
@@ -1123,9 +1127,6 @@ void engine_make_gravity_tasks(struct engine *e) {
       if (cell_are_neighbours(ci, cj))
         scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, ci,
                           cj, 1);
-      else
-        scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, 0, 0, ci,
-                          cj, 1);
     }
   }
 }
@@ -1289,6 +1290,14 @@ void engine_link_gravity_tasks(struct engine *e) {
   const int nodeID = e->nodeID;
   const int nr_tasks = sched->nr_tasks;
 
+  /* Add one task gathering all the multipoles */
+  struct task *gather = scheduler_addtask(
+      sched, task_type_grav_gather_m, task_subtype_none, 0, 0, NULL, NULL, 0);
+
+  /* And one task performing the FFT */
+  struct task *fft = scheduler_addtask(sched, task_type_grav_fft,
+                                       task_subtype_none, 0, 0, NULL, NULL, 0);
+
   for (int k = 0; k < nr_tasks; k++) {
 
     /* Get a pointer to the task. */
@@ -1297,11 +1306,21 @@ void engine_link_gravity_tasks(struct engine *e) {
     /* Skip? */
     if (t->skip) continue;
 
+    /* Multipole construction */
+    if (t->type == task_type_grav_up) {
+      scheduler_addunlock(sched, t, gather);
+      scheduler_addunlock(sched, t, fft);
+    }
+
     /* Long-range interaction */
     if (t->type == task_type_grav_mm) {
 
-      engine_make_gravity_dependencies(sched, t, t->ci);
-      engine_make_gravity_dependencies(sched, t, t->cj);
+      /* Gather the multipoles --> mm interaction --> kick */
+      scheduler_addunlock(sched, gather, t);
+      scheduler_addunlock(sched, t, t->ci->super->kick);
+
+      /* init --> mm interaction */
+      scheduler_addunlock(sched, t->ci->super->init, t);
     }
 
     /* Self-interaction? */
@@ -2205,6 +2224,8 @@ void engine_init_particles(struct engine *e) {
 
     mask |= 1 << task_type_grav_up;
     mask |= 1 << task_type_grav_mm;
+    mask |= 1 << task_type_grav_gather_m;
+    mask |= 1 << task_type_grav_fft;
     mask |= 1 << task_type_self;
     mask |= 1 << task_type_pair;
     mask |= 1 << task_type_sub;
@@ -2342,6 +2363,8 @@ void engine_step(struct engine *e) {
 
     mask |= 1 << task_type_grav_up;
     mask |= 1 << task_type_grav_mm;
+    mask |= 1 << task_type_grav_gather_m;
+    mask |= 1 << task_type_grav_fft;
     mask |= 1 << task_type_self;
     mask |= 1 << task_type_pair;
     mask |= 1 << task_type_sub;
diff --git a/src/runner.c b/src/runner.c
index 772e0478cc..bc4fab5b60 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1109,21 +1109,16 @@ void *runner_main(void *data) {
         case task_type_recv:
           runner_do_recv_cell(r, ci, 1);
           break;
-        /* case task_type_grav_pp: */
-        /*   if (t->cj == NULL) */
-        /*     runner_doself_grav(r, t->ci); */
-        /*   else */
-        /*     runner_dopair_grav(r, t->ci, t->cj); */
-        /*   break; */
         case task_type_grav_mm:
-          runner_do_grav_mm(r, t->ci, t->cj);
+          runner_do_grav_mm(r, t->ci, 1);
           break;
         case task_type_grav_up:
           runner_do_grav_up(r, t->ci);
           break;
-        /* case task_type_grav_down: */
-        /*   runner_dograv_down(r, t->ci); */
-        /*   break; */
+        case task_type_grav_gather_m:
+          break;
+        case task_type_grav_fft:
+          break;
         case task_type_grav_external:
           runner_do_grav_external(r, t->ci, 1);
           break;
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 2a9d126114..b4d06479d5 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -444,14 +444,7 @@ static void runner_dosub_grav(struct runner *r, struct cell *ci,
   }
 }
 
-static void runner_do_grav_mm(struct runner *r, struct cell *ci,
-                              struct cell *cj) {
-
-#ifdef SANITY_CHECKS
-  if (cell_are_neighbours(ci, cj)) {
-    error("Non-neighbouring cells in mm task !");
-  }
-#endif
+static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
 
 #if ICHECK > 0
   for (int pid = 0; pid < ci->gcount; pid++) {
@@ -459,24 +452,30 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci,
     /* Get a hold of the ith part in ci. */
     struct gpart *restrict gp = &ci->gparts[pid];
 
-    if (gp->id == -ICHECK)
-      message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id, cj->loc[0],
-              cj->loc[1], cj->loc[2], cj->h[0], cj->gcount);
-  }
-
-  for (int pid = 0; pid < cj->gcount; pid++) {
-
-    /* Get a hold of the ith part in ci. */
-    struct gpart *restrict gp = &cj->gparts[pid];
-
     if (gp->id == -ICHECK)
       message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id, ci->loc[0],
               ci->loc[1], ci->loc[2], ci->h[0], ci->gcount);
   }
 #endif
 
-  runner_dopair_grav_pm(r, ci, cj);
-  runner_dopair_grav_pm(r, cj, ci);
+  /* Recover the list of top-level cells */
+  const struct engine *e = r->e;
+  struct cell *cells = e->s->cells;
+  const int nr_cells = e->s->nr_cells;
+  const int ti_current = e->ti_current;
+
+  /* Anything to do here? */
+  if (ci->ti_end_min > ti_current) return;
+
+  /* Loop over all the cells and go for a p-m interaction if far enough */
+  for (int i = 0; i < nr_cells; ++i) {
+
+    struct cell *cj = &cells[i];
+
+    if (ci == cj) continue;
+
+    if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_pm(r, ci, cj);
+  }
 }
 
 #endif /* SWIFT_RUNNER_DOIACT_GRAV_H */
diff --git a/src/scheduler.c b/src/scheduler.c
index 6dac16b40c..3549e107cb 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -520,10 +520,9 @@ void scheduler_splittasks(struct scheduler *s) {
 
       /* Get a handle on the cells involved. */
       struct cell *ci = t->ci;
-      struct cell *cj = t->cj;
 
       /* Safety thing */
-      if (ci->gcount == 0 || cj->gcount == 0) t->type = task_type_none;
+      if (ci->gcount == 0) t->type = task_type_none;
 
     } /* gravity interaction? */
 
diff --git a/src/task.c b/src/task.c
index 8b0f0de8ac..68e3f4ee91 100644
--- a/src/task.c
+++ b/src/task.c
@@ -47,10 +47,11 @@
 
 /* Task type names. */
 const char *taskID_names[task_type_count] = {
-    "none",      "sort",       "self",       "pair",    "sub",
-    "init",      "ghost",      "drift",      "kick",    "kick_fixdt",
-    "send",      "recv",       "grav_mm",    "grav_up", "grav_external",
-    "part_sort", "gpart_sort", "split_cell", "rewait"};
+    "none",    "sort",          "self",          "pair",       "sub",
+    "init",    "ghost",         "drift",         "kick",       "kick_fixdt",
+    "send",    "recv",          "grav_gather_m", "grav_fft",   "grav_mm",
+    "grav_up", "grav_external", "part_sort",     "gpart_sort", "split_cell",
+    "rewait"};
 
 const char *subtaskID_names[task_type_count] = {"none", "density", "force",
                                                 "grav"};
@@ -151,7 +152,6 @@ void task_unlock(struct task *t) {
 
     case task_type_grav_mm:
       cell_gunlocktree(ci);
-      cell_gunlocktree(cj);
       break;
     default:
       break;
@@ -252,12 +252,7 @@ int task_lock(struct task *t) {
       break;
 
     case task_type_grav_mm:
-      if (ci->ghold || cj->ghold) return 0;
       if (cell_glocktree(ci) != 0) return 0;
-      if (cell_glocktree(cj) != 0) {
-        cell_gunlocktree(ci);
-        return 0;
-      }
 
     default:
       break;
diff --git a/src/task.h b/src/task.h
index bd0209d912..9f3429518d 100644
--- a/src/task.h
+++ b/src/task.h
@@ -45,10 +45,10 @@ enum task_types {
   task_type_kick_fixdt,
   task_type_send,
   task_type_recv,
-  /* task_type_grav_pp, */
+  task_type_grav_gather_m,
+  task_type_grav_fft,
   task_type_grav_mm,
   task_type_grav_up,
-  /* task_type_grav_down, */
   task_type_grav_external,
   task_type_part_sort,
   task_type_gpart_sort,
-- 
GitLab