Force an emergency re-sorting of all the star particles below a given super...

Force an emergency re-sorting of all the star particles below a given super cell when star formation happens.

src/cell.c

@@ -1990,8 +1990,8 @@ void cell_activate_hydro_sorts(struct cell *c, int sid, struct scheduler *s) {
     }
   }
 
-  //message("c->hydro.sorted=%d", c->hydro.sorted);
-  
+  // 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) {
 
@@ -2023,11 +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->stars.sorts == NULL) 
+        if (parent->stars.sorts == NULL)
           error("Trying to activate un-existing parents->stars.sorts");
 #endif
-	scheduler_activate(s, parent->stars.sorts);
-	if (parent->nodeID == engine_rank) cell_activate_drift_spart(parent, s);
+        scheduler_activate(s, parent->stars.sorts);
+        if (parent->nodeID == engine_rank) cell_activate_drift_spart(parent, s);
         break;
       }
     }
@@ -2052,13 +2052,13 @@ void cell_activate_stars_sorts(struct cell *c, int sid, struct scheduler *s) {
     }
   }
 
-  //message("c->stars.sorted=%d", c->stars.sorted);
-  
+  // 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");
-    
+    // message("bbb");
+
     atomic_or(&c->stars.do_sort, (1 << sid));
     cell_activate_stars_sorts_up(c, s);
   }

src/cell.h

@@ -1007,19 +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*/
 
-  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!");
+  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;
+  return 0;
 }
 /**
  * @brief Have star particles in a pair of cells moved too much and require a
@@ -1036,12 +1036,11 @@ cell_need_rebuild_for_stars_pair(const struct cell *ci, const struct cell *cj) {
   /* Note ci->dmin == cj->dmin */
   /* 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) {
-    //error("Need rebuild stars!");
+  if (kernel_gamma * max(ci->stars.h_max, cj->hydro.h_max) +
+          ci->stars.dx_max_part + cj->hydro.dx_max_part >
+      cj->dmin) {
+    // error("Need rebuild stars!");
     return 1;
-
   }
   return 0;
 }

src/engine.c

@@ -1910,8 +1910,9 @@ void engine_print_task_counts(const struct engine *e) {
   for (int k = 0; k < nr_tasks; k++) {
     if (tasks[k].skip) {
       counts[task_type_count] += 1;
-      //if (e->step == 0) message("Skipped %s/%s", taskID_names[tasks[k].type], subtaskID_names[tasks[k].subtype]);
-    }    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,
@@ -2676,9 +2677,8 @@ 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_drift_spart || t->type == task_type_kick1 ||
+        t->type == task_type_kick2 || t->type == task_type_timestep ||
         t->type == task_type_timestep_limiter ||
         t->subtype == task_subtype_force ||
         t->subtype == task_subtype_limiter || t->subtype == task_subtype_grav ||
@@ -2712,8 +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
-	|| t->type == task_type_drift_spart)
+    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

@@ -911,8 +911,8 @@ 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 = scheduler_addtask(
-          s, task_type_stars_sort, task_subtype_none, 0, 0, c, NULL);
+      c->stars.sorts = scheduler_addtask(s, task_type_stars_sort,
+                                         task_subtype_none, 0, 0, c, NULL);
     }
 
     /* Local tasks only... */

src/engine_marktasks.c

@@ -234,7 +234,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
 
       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);
+        if (cell_is_active_stars(ci, e)) scheduler_activate(s, t);
       }
 
       /* Activate the gravity drift */
@@ -325,9 +325,10 @@ 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) &&
+           t_subtype == task_subtype_stars_feedback) &&
           ((ci_active_stars && ci_nodeID == nodeID) ||
-	   (cj_active_stars && cj_nodeID == nodeID))) {  // MATTHIEU: check MPI condition here
+           (cj_active_stars &&
+            cj_nodeID == nodeID))) {  // MATTHIEU: check MPI condition here
 
         scheduler_activate(s, t);
 
@@ -347,8 +348,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
             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);
-	    
+            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);
@@ -367,9 +368,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
             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);
 
-	    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);
@@ -378,9 +378,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_pair &&
-		 t_subtype == task_subtype_stars_density) {
+                 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);
+          cell_activate_subcell_hydro_tasks(t->ci, t->cj, s);
         }
       }
 
@@ -503,9 +503,9 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
       /* Only interested in stars_density tasks as of here. */
       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;
+      /* 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;
 
 #ifdef WITH_MPI
         engine_activate_stars_mpi(e, s, ci, cj);

src/runner.c

@@ -119,6 +119,8 @@
 #include "runner_doiact_stars.h"
 #undef FUNCTION
 
+int emergency_sort = 0;
+
 /**
  * @brief Intermediate task after the density to check that the smoothing
  * lengths are correct.
@@ -481,10 +483,31 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
   if (timer) TIMER_TOC(timer_do_cooling);
 }
 
+void runner_clear_stars_sort_flags(struct runner *r, struct cell *c) {
+
+  if (c->split) {
+    for (int k = 0; k < 8; k++)
+      if (c->progeny[k] != NULL)
+        runner_clear_stars_sort_flags(r, c->progeny[k]);
+  }
+
+  c->stars.sorted = 0;
+  for (int i = 0; i < 13; i++) {
+    c->stars.sort[i] = NULL;
+  }
+
+  if (c->hydro.super == c) {
+    for (int i = 0; i < 13; i++) {
+      free(c->stars.sort[i]);
+    }
+  }
+}
+
 /**
  *
  */
-void runner_do_star_formation(struct runner *r, struct cell *c, int timer, int* formed_stars) {
+void runner_do_star_formation(struct runner *r, struct cell *c, int timer,
+                              int *formed_stars) {
 
   struct engine *e = r->e;
   const struct cosmology *cosmo = e->cosmology;
@@ -508,7 +531,8 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer, int*
   /* Recurse? */
   if (c->split) {
     for (int k = 0; k < 8; k++)
-      if (c->progeny[k] != NULL) runner_do_star_formation(r, c->progeny[k], 0, formed_stars);
+      if (c->progeny[k] != NULL)
+        runner_do_star_formation(r, c->progeny[k], 0, formed_stars);
   } else {
 
     /* Loop over the gas particles in this cell. */
@@ -554,8 +578,12 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer, int*
             star_formation_copy_properties(p, xp, sp, e, sf_props, cosmo,
                                            with_cosmology);
 
-	    message("STAR FORMED!!!! super->ID=%d", c->super->cellID);
-	    (*formed_stars)++;
+            message(
+                "STAR FORMED!!!! c->ID=%d super->ID=%d c->depth=%d"
+                "c->maxdepth=%d",
+                c->cellID, c->super->cellID, c->depth, c->maxdepth);
+
+            (*formed_stars)++;
           }
 
         } else { /* Are we not star-forming? */
@@ -570,9 +598,17 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer, int*
   }
 
   if (timer && *formed_stars > 0) {
-    runner_do_stars_sort(r, c, 0x1FFF, 0, 0);
+    message(
+        "Emergency sort! c->ID=%d c->depth=%d c->maxdepth=%d c=%p c->super=%p",
+        c->cellID, c->depth, c->maxdepth, c, c->hydro.super);
+
+    runner_clear_stars_sort_flags(r, c);
+
+    emergency_sort = 0;
+    runner_do_stars_sort(r, c, 0x1FFF, 1, 0);
+    emergency_sort = 0;
   }
-  
+
   if (timer) TIMER_TOC(timer_do_star_formation);
 }
 
@@ -955,6 +991,8 @@ void runner_do_stars_sort(struct runner *r, struct cell *c, int flags,
   if (c->stars.sorted == 0) c->stars.ti_sort = r->e->ti_current;
 #endif
 
+  if (emergency_sort) message("flag=%d", flags);
+
   /* start by allocating the entry arrays in the requested dimensions. */
   for (int j = 0; j < 13; j++) {
     if ((flags & (1 << j)) && c->stars.sort[j] == NULL) {
@@ -985,6 +1023,10 @@ void runner_do_stars_sort(struct runner *r, struct cell *c, int flags,
     c->stars.dx_max_sort = dx_max_sort;
     c->stars.dx_max_sort_old = dx_max_sort_old;
 
+    if (emergency_sort) {
+      message("c->id=%d (split) sorting", c->cellID);
+    }
+
     /* Loop over the 13 different sort arrays. */
     for (int j = 0; j < 13; j++) {
 
@@ -1056,6 +1098,13 @@ void runner_do_stars_sort(struct runner *r, struct cell *c, int flags,
   /* Otherwise, just sort. */
   else {
 
+    if (emergency_sort) {
+      message("c->id=%d (leaf) sorting", c->cellID);
+      /* for (int j = 0; j < 13; j++) */
+      /* 	if (flags & (1 << j)) */
+      /* 	  message("Sorting direction %d", j); */
+    }
+
     /* Reset the sort distance */
     if (c->stars.sorted == 0) {
 
@@ -3216,12 +3265,10 @@ void *runner_main(void *data) {
         case task_type_cooling:
           runner_do_cooling(r, t->ci, 1);
           break;
-        case task_type_star_formation:
-	  {
-	    int formed_stars = 0;
-	    runner_do_star_formation(r, t->ci, 1, &formed_stars);
-	  }
-          break;
+        case task_type_star_formation: {
+          int formed_stars = 0;
+          runner_do_star_formation(r, t->ci, 1, &formed_stars);
+        } break;
         default:
           error("Unknown/invalid task type (%d).", t->type);
       }

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/runner_doiact_stars.h

@@ -1043,9 +1043,12 @@ void DOSELF1_BRANCH_STARS(struct runner *r, struct cell *c) {
             "cj->nodeID=%d "                                                \
             "ci->nodeID=%d d=%e sort_j[pjd].d=%e cj->" #TYPE                \
             ".dx_max_sort=%e "                                              \
-            "cj->" #TYPE ".dx_max_sort_old=%e, %i",                         \
+            "cj->" #TYPE                                                    \
+            ".dx_max_sort_old=%e, cellID=%i super->cellID=%i"               \
+            "cj->depth=%d cj->maxdepth=%d",                                 \
             cj->nodeID, ci->nodeID, d, sort_j[pjd].d, cj->TYPE.dx_max_sort, \
-            cj->TYPE.dx_max_sort_old, cj->cellID);                          \
+            cj->TYPE.dx_max_sort_old, cj->cellID, cj->hydro.super->cellID,  \
+            cj->depth, cj->maxdepth);                                       \
     }                                                                       \
   })
 
@@ -1113,12 +1116,12 @@ void DOPAIR1_BRANCH_STARS(struct runner *r, struct cell *ci, struct cell *cj) {
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (do_ci) {
-    RUNNER_CHECK_SORT(hydro, part, cj, ci, sid);
+    // RUNNER_CHECK_SORT(hydro, part, cj, ci, sid);
     RUNNER_CHECK_SORT(stars, spart, ci, cj, sid);
   }
 
   if (do_cj) {
-    RUNNER_CHECK_SORT(hydro, part, ci, cj, sid);
+    // RUNNER_CHECK_SORT(hydro, part, ci, cj, sid);
     RUNNER_CHECK_SORT(stars, spart, cj, ci, sid);
   }
 #endif /* SWIFT_DEBUG_CHECKS */

src/scheduler.h

@@ -123,7 +123,7 @@ __attribute__((always_inline)) INLINE static void scheduler_activate(
 
   /* 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/stars/Default/stars_iact.h

@@ -78,8 +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!");
-  
+  // message("Feedback!");
+
 #ifdef DEBUG_INTERACTIONS_STARS
   /* Update ngb counters */
   if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_STARS)