First implementation

src/cell.c

@@ -1028,7 +1028,8 @@ void cell_clear_drift_flags(struct cell *c, void *data) {
                          cell_flag_do_bh_drift | cell_flag_do_bh_sub_drift |
                          cell_flag_do_stars_drift |
                          cell_flag_do_stars_sub_drift |
-                         cell_flag_do_sink_drift | cell_flag_do_sink_sub_drift);
+                         cell_flag_do_sink_drift | cell_flag_do_sink_sub_drift |
+                  cell_flag_do_recursion_gravity_self | cell_flag_do_recursion_gravity_pair);
 }
 
 /**

src/cell.h

@@ -337,7 +337,9 @@ enum cell_flags {
   cell_flag_do_stars_resort = (1UL << 15),
   cell_flag_has_tasks = (1UL << 16),
   cell_flag_do_hydro_sync = (1UL << 17),
-  cell_flag_do_hydro_sub_sync = (1UL << 18)
+  cell_flag_do_hydro_sub_sync = (1UL << 18),
+  cell_flag_do_recursion_gravity_self = (1UL << 19),
+  cell_flag_do_recursion_gravity_pair = (1UL << 20)
 };
 
 /**
@@ -568,7 +570,7 @@ void cell_activate_sink_formation_tasks(struct cell *c, struct scheduler *s);
 void cell_activate_subcell_hydro_tasks(struct cell *ci, struct cell *cj,
                                        struct scheduler *s,
                                        const int with_timestep_limiter);
-void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
+int cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
                                       struct scheduler *s);
 void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
                                        struct scheduler *s,
@@ -1375,4 +1377,31 @@ __attribute__((always_inline)) INLINE void cell_assign_cell_index(
 #endif
 }
 
+/**
+ * @brief Clear the flags related to the gravity recursion in the current cell and parents.
+ */
+INLINE static void cell_clear_flag_recursion_gravity(
+    struct cell *c) {
+
+  /* /\* Check if the cell is already cleared. *\/ */
+  /* if (!cell_get_flag(c, cell_flag_do_recursion_gravity_self) && */
+  /*     !cell_get_flag(c, cell_flag_do_recursion_gravity_pair)) */
+  /*   return; */
+
+  /* Clear the cell */
+  cell_clear_flag(c, cell_flag_do_recursion_gravity_self |
+                  cell_flag_do_recursion_gravity_pair);
+
+  /* /\* Now do the same with the parents. *\/ */
+  /* if (c->parent != NULL) { */
+  /*   cell_clear_flag_recursion_gravity(c->parent); */
+  /* } */
+
+  /* Now the same for the children */
+  for(int i = 0; i < 8; i++) {
+    if (c->progeny[i] != NULL)
+      cell_clear_flag_recursion_gravity(c->progeny[i]);
+  }
+}
+
 #endif /* SWIFT_CELL_H */

src/cell_drift.c

@@ -284,6 +284,7 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
   if (c->grav.count == 0) {
     /* Clear the drift flags. */
     cell_clear_flag(c, cell_flag_do_grav_drift | cell_flag_do_grav_sub_drift);
+    cell_clear_flag_recursion_gravity(c);
 
     /* Update the time of the last drift */
     c->grav.ti_old_part = ti_current;
@@ -392,6 +393,7 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
 
   /* Clear the drift flags. */
   cell_clear_flag(c, cell_flag_do_grav_drift | cell_flag_do_grav_sub_drift);
+  cell_clear_flag_recursion_gravity(c);
 }
 
 /**

src/cell_unskip.c

@@ -1121,7 +1121,7 @@ void cell_activate_subcell_sinks_tasks(struct cell *ci, struct cell *cj,
  * @param cj The second #cell we recurse in.
  * @param s The task #scheduler.
  */
-void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
+int cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
                                      struct scheduler *s) {
   /* Some constants */
   const struct space *sp = s->space;
@@ -1130,7 +1130,11 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
   /* Self interaction? */
   if (cj == NULL) {
     /* Do anything? */
-    if (ci->grav.count == 0 || !cell_is_active_gravity(ci, e)) return;
+    if (ci->grav.count == 0 || !cell_is_active_gravity(ci, e)) return 1;
+
+    /* /\* Is it already done? *\/ */
+    /* if (cell_get_flag(ci, cell_flag_do_recursion_gravity_self)) return 1; */
+    /* cell_set_flag(ci, cell_flag_do_recursion_gravity_self); */
 
     /* Recurse? */
     if (ci->split) {
@@ -1148,14 +1152,19 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
       /* We have reached the bottom of the tree: activate gpart drift */
       cell_activate_drift_gpart(ci, s);
     }
+    return 1;
   }
 
   /* Pair interaction */
   else {
     /* Anything to do here? */
     if (!cell_is_active_gravity(ci, e) && !cell_is_active_gravity(cj, e))
-      return;
-    if (ci->grav.count == 0 || cj->grav.count == 0) return;
+      return 1;
+    if (ci->grav.count == 0 || cj->grav.count == 0) return 1;
+
+    /* Is it already done? */
+    if (cell_get_flag(ci, cell_flag_do_recursion_gravity_pair) &&
+        cell_get_flag(cj, cell_flag_do_recursion_gravity_pair)) return 1;
 
     /* Atomically drift the multipole in ci */
     lock_lock(&ci->grav.mlock);
@@ -1170,9 +1179,8 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
     /* Can we use multipoles ? */
     if (cell_can_use_pair_mm(ci, cj, e, sp, /*use_rebuild_data=*/0,
                              /*is_tree_walk=*/1)) {
-
       /* Ok, no need to drift anything */
-      return;
+      return 0;
     }
     /* Otherwise, activate the gpart drifts if we are at the bottom. */
     else if (!ci->split && !cj->split) {
@@ -1181,6 +1189,11 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
         if (ci->nodeID == engine_rank) cell_activate_drift_gpart(ci, s);
         if (cj->nodeID == engine_rank) cell_activate_drift_gpart(cj, s);
       }
+      /* Flag the cells as being done. */
+      cell_set_flag(ci, cell_flag_do_recursion_gravity_pair);
+      cell_set_flag(cj, cell_flag_do_recursion_gravity_pair);
+      return 1;
+
     }
     /* Ok, we can still recurse */
     else {
@@ -1190,20 +1203,45 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
       const double ri_max = multi_i->r_max;
       const double rj_max = multi_j->r_max;
 
+      int ci_number_children = 0;
+      if (ci->split) {
+        for(int k = 0; k < 8; k++) {
+          if (ci->progeny[k] != NULL)
+            ci_number_children += 1;
+        }
+      }
+      int cj_number_children = 0;
+      if (cj->split) {
+        for(int k = 0; k < 8; k++) {
+          if (cj->progeny[k] != NULL)
+            cj_number_children += 1;
+        }
+      }
+
       if (ri_max > rj_max) {
         if (ci->split) {
           /* Loop over ci's children */
+          int rec = 0;
           for (int k = 0; k < 8; k++) {
             if (ci->progeny[k] != NULL)
-              cell_activate_subcell_grav_tasks(ci->progeny[k], cj, s);
+              rec += cell_activate_subcell_grav_tasks(ci->progeny[k], cj, s);
           }
+          /* Flag the cells as being done. */
+          if (rec == ci_number_children)
+            cell_set_flag(ci, cell_flag_do_recursion_gravity_pair);
+          return rec == ci_number_children;
 
         } else if (cj->split) {
           /* Loop over cj's children */
+          int rec = 0;
           for (int k = 0; k < 8; k++) {
             if (cj->progeny[k] != NULL)
-              cell_activate_subcell_grav_tasks(ci, cj->progeny[k], s);
+              rec += cell_activate_subcell_grav_tasks(ci, cj->progeny[k], s);
           }
+          /* Flag the cells as being done. */
+          if (rec == cj_number_children)
+            cell_set_flag(cj, cell_flag_do_recursion_gravity_pair);
+          return rec == cj_number_children;
 
         } else {
           error("Fundamental error in the logic");
@@ -1211,17 +1249,27 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
       } else if (rj_max >= ri_max) {
         if (cj->split) {
           /* Loop over cj's children */
+          int rec = 0;
           for (int k = 0; k < 8; k++) {
             if (cj->progeny[k] != NULL)
-              cell_activate_subcell_grav_tasks(ci, cj->progeny[k], s);
+              rec += cell_activate_subcell_grav_tasks(ci, cj->progeny[k], s);
           }
+          /* Flag the cells as being done. */
+          if (rec == cj_number_children)
+            cell_set_flag(cj, cell_flag_do_recursion_gravity_pair);
+          return rec == cj_number_children;
 
         } else if (ci->split) {
           /* Loop over ci's children */
+          int rec = 0;
           for (int k = 0; k < 8; k++) {
             if (ci->progeny[k] != NULL)
-              cell_activate_subcell_grav_tasks(ci->progeny[k], cj, s);
+              rec += cell_activate_subcell_grav_tasks(ci->progeny[k], cj, s);
           }
+          /* Flag the cells as being done. */
+          if (rec == ci_number_children)
+            cell_set_flag(ci, cell_flag_do_recursion_gravity_pair);
+          return rec == ci_number_children;
 
         } else {
           error("Fundamental error in the logic");
@@ -1229,6 +1277,8 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
       }
     }
   }
+  error("Not happening");
+  return -1;
 }
 
 /**

src/engine_unskip.c

@@ -365,6 +365,7 @@ void engine_do_unskip_mapper(void *map_data, int num_elements,
   }
 }
 
+
 /**
  * @brief Unskip all the tasks that act on active cells at this time.
  *
@@ -385,6 +386,10 @@ void engine_unskip(struct engine *e) {
   const int with_black_holes = e->policy & engine_policy_black_holes;
   const int with_rt = e->policy & engine_policy_rt;
 
+  for(int i = 0; i < s->nr_cells; i++) {
+    cell_clear_flag_recursion_gravity(&s->cells_top[i]);
+  }
+
 #ifdef WITH_PROFILER
   static int count = 0;
   char filename[100];

src/runner_main.c

@@ -183,6 +183,14 @@ void *runner_main(void *data) {
 #endif
 
       /* Different types of tasks... */
+      cell_clear_flag(ci, cell_flag_do_recursion_gravity_self |
+                      cell_flag_do_recursion_gravity_pair);
+      if (cj != NULL) {
+        cell_clear_flag(cj, cell_flag_do_recursion_gravity_self |
+                        cell_flag_do_recursion_gravity_pair);
+
+      }
+
       switch (t->type) {
         case task_type_self:
           if (t->subtype == task_subtype_density)