diff --git a/src/cell.c b/src/cell.c index 4ba22beffc592ddd2eaf189409aabcbb763a5192..8f3fcf8ed7ce10fa087b92f7bc41f5a53b23aa90 100644 --- a/src/cell.c +++ b/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); } diff --git a/src/cell.h b/src/cell.h index f9182ae3921e077a1634179ae8128eea524b3e7f..3fa3181933c9ff118e198518ce0b20a0430060e0 100644 --- a/src/cell.h +++ b/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; } diff --git a/src/engine.c b/src/engine.c index 84b0aff8dcffd6b5b4a9a397809fa0d1c3bb0352..07a448c2afe6348059914d685e06ddfe3b8c2311 100644 --- a/src/engine.c +++ b/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; } diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c index 3da1aa3bb800bc8b155d745c2d21bf58dd68299b..55451bcb35d90899dc8e8569471db802254cba4d 100644 --- a/src/engine_maketasks.c +++ b/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... */ diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c index eddc4b2dcae6779c775afe3d8121327c98d2887b..b7a930979930088f2258a94fb1ab6c0e3a6eca7a 100644 --- a/src/engine_marktasks.c +++ b/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); diff --git a/src/runner.c b/src/runner.c index 8d8a2e6b6a648ba1b3c08376a9d84582cef6fa8d..ceae7616cf8bc4d90af0c5b74d1f6f3109c6dece 100644 --- a/src/runner.c +++ b/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); } diff --git a/src/runner_doiact.h b/src/runner_doiact.h index 0bb6fd18348fcd9d7808ec20f370f0ccf86037af..861798b70b8ba90b9267375253bd8570baec3e9a 100644 --- a/src/runner_doiact.h +++ b/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) diff --git a/src/runner_doiact_stars.h b/src/runner_doiact_stars.h index c23f3f6f0aaf8cfff6a285e14d064191db837b1a..a8cd8a94b1bb409647a73904b32804a3ddd0a73c 100644 --- a/src/runner_doiact_stars.h +++ b/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 */ diff --git a/src/scheduler.h b/src/scheduler.h index 3cd8b3e1fc84dd702a670c9fc5103a754141ffe9..bc87f27f73c3760f4292f0f5bf6534e995b4abee 100644 --- a/src/scheduler.h +++ b/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); diff --git a/src/stars/Default/stars_iact.h b/src/stars/Default/stars_iact.h index 35628e30643768f1d29392ed932076dae898351a..56cd570d6c1338c340ac809ae5302c66de723c85 100644 --- a/src/stars/Default/stars_iact.h +++ b/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)