Added init_grav task to initalize the multipole structure independently of the part/gparts

src/cell.c

@@ -1356,6 +1356,7 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
   if (c->extra_ghost != NULL) scheduler_activate(s, c->extra_ghost);
   if (c->ghost != NULL) scheduler_activate(s, c->ghost);
   if (c->init != NULL) scheduler_activate(s, c->init);
+  if (c->init_grav != NULL) scheduler_activate(s, c->init_grav);
   if (c->drift != NULL) scheduler_activate(s, c->drift);
   if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
   if (c->kick2 != NULL) scheduler_activate(s, c->kick2);

src/cell.h

@@ -148,9 +148,12 @@ struct cell {
   /*! Linked list of the tasks computing this cell's gravity forces. */
   struct link *grav;
 
-  /*! The initialistation task */
+  /*! The particle initialistation task */
   struct task *init;
 
+  /*! The multipole initialistation task */
+  struct task *init_grav;
+
   /*! The ghost task */
   struct task *ghost;
 

src/engine.c

@@ -170,6 +170,10 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
 
       if (is_self_gravity) {
 
+        /* Initialisation of the multipoles */
+        c->init_grav = scheduler_addtask(s, task_type_init_grav,
+                                         task_subtype_none, 0, 0, c, NULL, 0);
+
         /* Gravity non-neighbouring pm calculations */
         c->grav_long_range = scheduler_addtask(
             s, task_type_grav_long_range, task_subtype_none, 0, 0, c, NULL, 0);
@@ -182,8 +186,8 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
         c->grav_down = scheduler_addtask(s, task_type_grav_down,
                                          task_subtype_none, 0, 0, c, NULL, 0);
 
-        scheduler_addunlock(s, c->init, c->grav_long_range);
-        scheduler_addunlock(s, c->init, c->grav_top_level);
+        scheduler_addunlock(s, c->init_grav, c->grav_long_range);
+        scheduler_addunlock(s, c->init_grav, c->grav_top_level);
         scheduler_addunlock(s, c->grav_long_range, c->grav_down);
         scheduler_addunlock(s, c->grav_top_level, c->grav_down);
         scheduler_addunlock(s, c->grav_down, c->kick2);
@@ -1892,6 +1896,7 @@ static inline void engine_make_self_gravity_dependencies(
 
   /* init --> gravity --> grav_down --> kick */
   scheduler_addunlock(sched, c->super->init, gravity);
+  scheduler_addunlock(sched, c->super->init_grav, gravity);
   scheduler_addunlock(sched, gravity, c->super->grav_down);
 }
 
@@ -2544,7 +2549,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
 
     /* Kick/Drift/Init? */
     else if (t->type == task_type_kick1 || t->type == task_type_kick2 ||
-             t->type == task_type_drift || t->type == task_type_init) {
+             t->type == task_type_drift || t->type == task_type_init ||
+             t->type == task_type_init_grav) {
       if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
     }
 
@@ -3252,6 +3258,8 @@ void engine_step(struct engine *e) {
   gravity_exact_force_check(e->s, e, 1e-1);
 #endif
 
+  /* Let's trigger a rebuild every-so-often for good measure */  // MATTHIEU
+                                                                 // improve
   if (!(e->policy & engine_policy_hydro) &&
       (e->policy & engine_policy_self_gravity) && e->step % 20 == 0)
     e->forcerebuild = 1;

src/runner.c

@@ -444,7 +444,44 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
 }
 
 /**
- * @brief Initialize the particles before the density calculation
+ * @brief Initialize the multipoles before the gravity calculation.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer 1 if the time is to be recorded.
+ */
+void runner_do_init_grav(struct runner *r, struct cell *c, int timer) {
+
+  const struct engine *e = r->e;
+
+  TIMER_TIC;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (!(e->policy & engine_policy_self_gravity))
+    error("Grav-init task called outside of self-gravity calculation");
+#endif
+
+  /* Anything to do here? */
+  if (!cell_is_active(c, e)) return;
+
+  /* Drift the multipole */
+  cell_drift_multipole(c, e);
+  
+  /* Reset the gravity acceleration tensors */
+  gravity_field_tensors_init(&c->multipole->pot);
+
+  /* Recurse? */
+  if (c->split) {
+    for (int k = 0; k < 8; k++) {
+      if (c->progeny[k] != NULL) runner_do_init_grav(r, c->progeny[k], 0);
+    }
+  }
+
+  if (timer) TIMER_TOC(timer_init_grav);
+}
+
+/**
+ * @brief Initialize the particles before the density calculation.
  *
  * @param r The runner thread.
  * @param c The cell.
@@ -464,10 +501,6 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
   /* Anything to do here? */
   if (!cell_is_active(c, e)) return;
 
-  /* Reset the gravity acceleration tensors */
-  if (e->policy & engine_policy_self_gravity)
-    gravity_field_tensors_init(&c->multipole->pot);
-
   /* Recurse? */
   if (c->split) {
     for (int k = 0; k < 8; k++)
@@ -1765,6 +1798,9 @@ void *runner_main(void *data) {
         case task_type_sort:
           runner_do_sort(r, ci, t->flags, 1);
           break;
+        case task_type_init_grav:
+          runner_do_init_grav(r, ci, 1);
+          break;
         case task_type_init:
           runner_do_init(r, ci, 1);
           break;

src/space.c

@@ -211,6 +211,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
     c->gcount = 0;
     c->scount = 0;
     c->init = NULL;
+    c->init_grav = NULL;
     c->extra_ghost = NULL;
     c->ghost = NULL;
     c->kick1 = NULL;

src/task.c

@@ -47,27 +47,15 @@
 #include "lock.h"
 
 /* Task type names. */
-const char *taskID_names[task_type_count] = {"none",
-                                             "sort",
-                                             "self",
-                                             "pair",
-                                             "sub_self",
-                                             "sub_pair",
-                                             "init",
-                                             "ghost",
-                                             "extra_ghost",
-                                             "drift",
-                                             "kick1",
-                                             "kick2",
-                                             "timestep",
-                                             "send",
-                                             "recv",
-                                             "grav_top_level",
-                                             "grav_long_range",
-                                             "grav_mm",
-                                             "grav_down",
-                                             "cooling",
-                                             "sourceterms"};
+const char *taskID_names[task_type_count] = {
+    "none",        "sort",           "self",
+    "pair",        "sub_self",       "sub_pair",
+    "init",        "init_grav",      "ghost",
+    "extra_ghost", "drift",          "kick1",
+    "kick2",       "timestep",       "send",
+    "recv",        "grav_top_level", "grav_long_range",
+    "grav_mm",     "grav_down",      "cooling",
+    "sourceterms"};
 
 /* Sub-task type names. */
 const char *subtaskID_names[task_subtype_count] = {
@@ -181,6 +169,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
         error("Task without particles");
       break;
 
+    case task_type_init_grav:
     case task_type_grav_top_level:
     case task_type_grav_long_range:
     case task_type_grav_mm:
@@ -283,6 +272,10 @@ void task_unlock(struct task *t) {
   /* Act based on task type. */
   switch (type) {
 
+    case task_type_init_grav:
+      cell_munlocktree(ci);
+      break;
+
     case task_type_init:
     case task_type_kick1:
     case task_type_kick2:
@@ -370,6 +363,11 @@ int task_lock(struct task *t) {
 #endif
       break;
 
+    case task_type_init_grav:
+      if (ci->mhold) return 0;
+      if (cell_mlocktree(ci) != 0) return 0;
+      break;
+
     case task_type_init:
     case task_type_kick1:
     case task_type_kick2:

src/task.h

@@ -43,6 +43,7 @@ enum task_types {
   task_type_sub_self,
   task_type_sub_pair,
   task_type_init,
+  task_type_init_grav,
   task_type_ghost,
   task_type_extra_ghost,
   task_type_drift,

src/timers.c

@@ -34,6 +34,7 @@ char *timers_names[timer_count] = {
     "none",
     "prepare",
     "init",
+    "init_grav",
     "drift",
     "kick1",
     "kick2",

src/timers.h

@@ -37,6 +37,7 @@ enum {
   timer_none = 0,
   timer_prepare,
   timer_init,
+  timer_init_grav,
   timer_drift,
   timer_kick1,
   timer_kick2,