Code seems to work when not forming stars. Switching this on next.

src/cell.c

@@ -1990,8 +1990,11 @@ void cell_activate_hydro_sorts(struct cell *c, int sid, struct scheduler *s) {
     }
   }
 
+  //message("c->hydro.sorted=%d", c->hydro.sorted);
+  
   /* Has this cell been sorted at all for the given sid? */
   if (!(c->hydro.sorted & (1 << sid)) || c->nodeID != engine_rank) {
+
     atomic_or(&c->hydro.do_sort, (1 << sid));
     cell_activate_hydro_sorts_up(c, s);
   }
@@ -2004,19 +2007,14 @@ void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
 
   if (c == c->hydro.super) {
 #ifdef SWIFT_DEBUG_CHECKS
-    if ((c->nodeID == engine_rank && c->stars.sorts_local == NULL) ||
-        (c->nodeID != engine_rank && c->stars.sorts_foreign == NULL))
+    if (c->stars.sorts == NULL)
       error("Trying to activate un-existing c->stars.sorts");
 #endif
+    scheduler_activate(s, c->stars.sorts);
     if (c->nodeID == engine_rank) {
-      scheduler_activate(s, c->stars.sorts_local);
       // MATTHIEU: to do: do we actually need both drifts here?
-      cell_activate_drift_part(c, s);
       cell_activate_drift_spart(c, s);
     }
-    if (c->nodeID != engine_rank) {
-      scheduler_activate(s, c->stars.sorts_foreign);
-    }
   } else {
 
     for (struct cell *parent = c->parent;
@@ -2025,20 +2023,11 @@ void cell_activate_stars_sorts_up(struct cell *c, struct scheduler *s) {
       parent->stars.do_sub_sort = 1;
       if (parent == c->hydro.super) {
 #ifdef SWIFT_DEBUG_CHECKS
-        if ((parent->nodeID == engine_rank &&
-             parent->stars.sorts_local == NULL) ||
-            (parent->nodeID != engine_rank &&
-             parent->stars.sorts_foreign == NULL))
+        if (parent->stars.sorts == NULL) 
           error("Trying to activate un-existing parents->stars.sorts");
 #endif
-        if (parent->nodeID == engine_rank) {
-          scheduler_activate(s, parent->stars.sorts_local);
-          cell_activate_drift_part(parent, s);
-          cell_activate_drift_spart(parent, s);
-        }
-        if (parent->nodeID != engine_rank) {
-          scheduler_activate(s, parent->stars.sorts_foreign);
-        }
+	scheduler_activate(s, parent->stars.sorts);
+	if (parent->nodeID == engine_rank) cell_activate_drift_spart(parent, s);
         break;
       }
     }
@@ -2063,8 +2052,13 @@ void cell_activate_stars_sorts(struct cell *c, int sid, struct scheduler *s) {
     }
   }
 
+  //message("c->stars.sorted=%d", c->stars.sorted);
+  
   /* Has this cell been sorted at all for the given sid? */
   if (!(c->stars.sorted & (1 << sid)) || c->nodeID != engine_rank) {
+
+    //message("bbb");
+    
     atomic_or(&c->stars.do_sort, (1 << sid));
     cell_activate_stars_sorts_up(c, s);
   }

src/cell.h

@@ -482,8 +482,7 @@ struct cell {
     struct link *feedback;
 
     /*! The task computing this cell's sorts before the density. */
-    struct task *sorts_local;
-    struct task *sorts_foreign;
+    struct task *sorts;
 
     /*! The drift task for sparts */
     struct task *drift;
@@ -1008,13 +1007,19 @@ cell_can_split_self_gravity_task(const struct cell *c) {
  */
 __attribute__((always_inline)) INLINE static int
 cell_need_rebuild_for_hydro_pair(const struct cell *ci, const struct cell *cj) {
-
+  
   /* Is the cut-off radius plus the max distance the parts in both cells have */
   /* moved larger than the cell size ? */
   /* Note ci->dmin == cj->dmin */
-  return (kernel_gamma * max(ci->hydro.h_max, cj->hydro.h_max) +
-              ci->hydro.dx_max_part + cj->hydro.dx_max_part >
-          cj->dmin);
+  /*return*/
+
+  if(kernel_gamma * max(ci->hydro.h_max, cj->hydro.h_max) +
+     ci->hydro.dx_max_part + cj->hydro.dx_max_part >
+     cj->dmin){
+    //error("Need rebuild hydro!");
+    return 1;
+  }
+return 0;
 }
 /**
  * @brief Have star particles in a pair of cells moved too much and require a
@@ -1029,9 +1034,16 @@ cell_need_rebuild_for_stars_pair(const struct cell *ci, const struct cell *cj) {
   /* Is the cut-off radius plus the max distance the parts in both cells have */
   /* moved larger than the cell size ? */
   /* Note ci->dmin == cj->dmin */
-  return (kernel_gamma * max(ci->stars.h_max, cj->hydro.h_max) +
+  /* return */
+
+  if(kernel_gamma * max(ci->stars.h_max, cj->hydro.h_max) +
               ci->stars.dx_max_part + cj->hydro.dx_max_part >
-          cj->dmin);
+     cj->dmin) {
+    //error("Need rebuild stars!");
+    return 1;
+
+  }
+  return 0;
 }
 
 /**

src/engine.c

@@ -1908,9 +1908,10 @@ void engine_print_task_counts(const struct engine *e) {
   int counts[task_type_count + 1];
   for (int k = 0; k <= task_type_count; k++) counts[k] = 0;
   for (int k = 0; k < nr_tasks; k++) {
-    if (tasks[k].skip)
+    if (tasks[k].skip) {
       counts[task_type_count] += 1;
-    else
+      //if (e->step == 0) message("Skipped %s/%s", taskID_names[tasks[k].type], subtaskID_names[tasks[k].subtype]);
+    }    else
       counts[(int)tasks[k].type] += 1;
   }
   message("Total = %d  (per cell = %d)", nr_tasks,
@@ -2675,6 +2676,7 @@ void engine_skip_force_and_kick(struct engine *e) {
 
     /* Skip everything that updates the particles */
     if (t->type == task_type_drift_part || t->type == task_type_drift_gpart ||
+	t->type == task_type_drift_spart ||
         t->type == task_type_kick1 || t->type == task_type_kick2 ||
         t->type == task_type_timestep ||
         t->type == task_type_timestep_limiter ||
@@ -2710,7 +2712,8 @@ void engine_skip_drift(struct engine *e) {
     struct task *t = &tasks[i];
 
     /* Skip everything that moves the particles */
-    if (t->type == task_type_drift_part || t->type == task_type_drift_gpart)
+    if (t->type == task_type_drift_part || t->type == task_type_drift_gpart
+	|| t->type == task_type_drift_spart)
       t->skip = 1;
   }
 

src/engine_maketasks.c

@@ -479,7 +479,7 @@ void engine_addtasks_recv_stars(struct engine *e, struct cell *c,
                                c->mpi.tag, 0, c, NULL);
 
     /* Need to sort task before feedback loop */
-    scheduler_addunlock(s, t_feed, c->hydro.super->stars.sorts_foreign);
+    scheduler_addunlock(s, t_feed, c->hydro.super->stars.sorts);
   }
 
   c->mpi.hydro.recv_xv = t_xv;
@@ -911,7 +911,7 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
         scheduler_addtask(s, task_type_sort, task_subtype_none, 0, 0, c, NULL);
 
     if (with_feedback) {
-      c->stars.sorts_local = scheduler_addtask(
+      c->stars.sorts = scheduler_addtask(
           s, task_type_stars_sort, task_subtype_none, 0, 0, c, NULL);
     }
 
@@ -1261,8 +1261,8 @@ void engine_count_and_link_tasks_mapper(void *map_data, int num_elements,
     if (t_type == task_type_stars_sort) {
       for (struct cell *finger = t->ci->parent; finger != NULL;
            finger = finger->parent) {
-        if (finger->stars.sorts_local != NULL)
-          scheduler_addunlock(sched, t, finger->stars.sorts_local);
+        if (finger->stars.sorts != NULL)
+          scheduler_addunlock(sched, t, finger->stars.sorts);
       }
     }
 
@@ -1626,7 +1626,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
     struct task *t = &((struct task *)map_data)[ind];
     const enum task_types t_type = t->type;
     const enum task_subtypes t_subtype = t->subtype;
-    const int flags = t->flags;
+    const long long flags = t->flags;
     struct cell *ci = t->ci;
     struct cell *cj = t->cj;
 
@@ -1812,7 +1812,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
           scheduler_addunlock(sched, ci->hydro.super->stars.drift,
                               t_star_density);
-          scheduler_addunlock(sched, ci->hydro.super->stars.sorts_local,
+          scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
                               t_star_density);
           scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
                               t_star_density);
@@ -1845,7 +1845,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
           scheduler_addunlock(sched, cj->hydro.super->stars.drift,
                               t_star_density);
-          scheduler_addunlock(sched, cj->hydro.super->stars.sorts_local,
+          scheduler_addunlock(sched, cj->hydro.super->stars.sorts,
                               t_star_density);
           scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
                               t_star_density);
@@ -1938,7 +1938,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
         scheduler_addunlock(sched, ci->hydro.super->stars.drift,
                             t_star_density);
-        scheduler_addunlock(sched, ci->hydro.super->stars.sorts_local,
+        scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
                             t_star_density);
         scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
                             t_star_density);
@@ -2058,7 +2058,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
           scheduler_addunlock(sched, ci->hydro.super->stars.drift,
                               t_star_density);
-          scheduler_addunlock(sched, ci->hydro.super->stars.sorts_local,
+          scheduler_addunlock(sched, ci->hydro.super->stars.sorts,
                               t_star_density);
           scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
                               t_star_density);
@@ -2091,7 +2091,7 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 
           scheduler_addunlock(sched, cj->hydro.super->stars.drift,
                               t_star_density);
-          scheduler_addunlock(sched, cj->hydro.super->stars.sorts_local,
+          scheduler_addunlock(sched, cj->hydro.super->stars.sorts,
                               t_star_density);
           scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
                               t_star_density);

src/engine_marktasks.c

@@ -157,7 +157,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
       /* Local pointer. */
       struct cell *ci = t->ci;
 
-      if (ci->nodeID != engine_rank) error("Non-local self task found");
+      if (ci->nodeID != nodeID) error("Non-local self task found");
 
       /* Activate the hydro drift */
       if (t_type == task_type_self && t_subtype == task_subtype_density) {
@@ -197,7 +197,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
       }
 
-#ifdef EXTRA_HYDRO_LOOP
       else if (t_type == task_type_self && t_subtype == task_subtype_gradient) {
         if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
       }
@@ -206,7 +205,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
                t_subtype == task_subtype_gradient) {
         if (cell_is_active_hydro(ci, e)) scheduler_activate(s, t);
       }
-#endif
 
       /* Activate the star density */
       else if (t_type == task_type_self &&
@@ -227,7 +225,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         }
       }
 
-      /* Activate the star feedback */
       else if (t_type == task_type_self &&
                t_subtype == task_subtype_stars_feedback) {
         if (cell_is_active_stars(ci, e)) {
@@ -235,13 +232,9 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         }
       }
 
-      /* Store current values of dx_max and h_max. */
       else if (t_type == task_type_sub_self &&
                t_subtype == task_subtype_stars_feedback) {
-        if (cell_is_active_stars(ci, e)) {
-          scheduler_activate(s, t);
-          cell_activate_subcell_stars_tasks(ci, NULL, s);
-        }
+        if (cell_is_active_stars(ci, e))  scheduler_activate(s, t);
       }
 
       /* Activate the gravity drift */
@@ -332,8 +325,9 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
 
       /* Stars density and feedback */
       if ((t_subtype == task_subtype_stars_density ||
-           t_subtype == task_subtype_stars_feedback) &&
-          (ci_active_stars || cj_active_stars)) {
+	  t_subtype == task_subtype_stars_feedback) &&
+          ((ci_active_stars && ci_nodeID == nodeID) ||
+	   (cj_active_stars && cj_nodeID == nodeID))) {  // MATTHIEU: check MPI condition here
 
         scheduler_activate(s, t);
 
@@ -351,12 +345,12 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
 
             /* Activate the hydro drift tasks. */
             if (ci_nodeID == nodeID) cell_activate_drift_spart(ci, s);
-
             if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
 
+	    message("do ci flags=%lld", t->flags);
+	    
             /* Check the sorts and activate them if needed. */
             cell_activate_hydro_sorts(cj, t->flags, s);
-
             cell_activate_stars_sorts(ci, t->flags, s);
           }
 
@@ -371,9 +365,11 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
 
             /* Activate the hydro drift tasks. */
             if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
-
             if (cj_nodeID == nodeID) cell_activate_drift_spart(cj, s);
 
+
+	    message("do cj flags=%lld", t->flags);
+	    
             /* Check the sorts and activate them if needed. */
             cell_activate_hydro_sorts(ci, t->flags, s);
             cell_activate_stars_sorts(cj, t->flags, s);
@@ -381,8 +377,10 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
         }
 
         /* Store current values of dx_max and h_max. */
-        else if (t_type == task_type_sub_pair) {
+        else if (t_type == task_type_sub_pair &&
+		 t_subtype == task_subtype_stars_density) {
           cell_activate_subcell_stars_tasks(t->ci, t->cj, s);
+	  cell_activate_subcell_hydro_tasks(t->ci, t->cj, s);
         }
       }
 
@@ -506,8 +504,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
       if (t->subtype == task_subtype_stars_density) {
 
         /* Too much particle movement? */
-        if (cell_need_rebuild_for_stars_pair(ci, cj)) *rebuild_space = 1;
-        if (cell_need_rebuild_for_stars_pair(cj, ci)) *rebuild_space = 1;
+        //if (cell_need_rebuild_for_stars_pair(ci, cj)) *rebuild_space = 1;
+        //if (cell_need_rebuild_for_stars_pair(cj, ci)) *rebuild_space = 1;
 
 #ifdef WITH_MPI
         engine_activate_stars_mpi(e, s, ci, cj);

src/runner_doiact.h

@@ -1410,7 +1410,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
         /* Compute the pairwise distance. */
         const float dx[3] = {pjx - pix, pjy - piy, pjz - piz};
         const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
-
+	
 #ifdef SWIFT_DEBUG_CHECKS
         /* Check that particles are in the correct frame after the shifts */
         if (pix > shift_threshold_x || pix < -shift_threshold_x)

src/scheduler.h

@@ -120,6 +120,10 @@ struct scheduler {
  */
 __attribute__((always_inline)) INLINE static void scheduler_activate(
     struct scheduler *s, struct task *t) {
+
+  /* if(t->type == task_type_stars_sort) */
+  /*   message("Activating stars sort!!!"); */
+  
   if (atomic_cas(&t->skip, 1, 0)) {
     t->wait = 0;
     int ind = atomic_inc(&s->active_count);

src/space.c

@@ -181,8 +181,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
       space_recycle_list(s, cell_rec_begin, cell_rec_end, multipole_rec_begin,
                          multipole_rec_end);
     c->hydro.sorts = NULL;
-    c->stars.sorts_local = NULL;
-    c->stars.sorts_foreign = NULL;
+    c->stars.sorts = NULL;
     c->nr_tasks = 0;
     c->grav.nr_mm_tasks = 0;
     c->hydro.density = NULL;

src/stars/Default/stars_iact.h

@@ -78,6 +78,8 @@ runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
   /* Get the time derivative for h. */
   si->feedback.h_dt -= mj * dvdr * ri / rhoj * wi_dr;
 
+  //message("Feedback!");
+  
 #ifdef DEBUG_INTERACTIONS_STARS
   /* Update ngb counters */
   if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_STARS)