diff --git a/src/cell.c b/src/cell.c
index 87f3f447283e49cc18020efffd7c44f470a2dacd..b566584807f6928d906207eb1c4f78181aa423ea 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -127,6 +127,7 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
temp->depth = c->depth + 1;
temp->split = 0;
temp->dx_max = 0.f;
+ temp->dx_max_sort = 0.f;
temp->nodeID = c->nodeID;
temp->parent = c;
temp->ti_old = c->ti_old;
@@ -908,8 +909,7 @@ int cell_is_drift_needed(struct cell *c, const struct engine *e) {
int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
/* Reset the sort flags if the particles have moved too much. */
- if (c->dx_max > space_maxreldx * c->h_max)
- c->sorted = 0;
+ if (c->dx_max_sort > space_maxreldx * c->h_max) c->sorted = 0;
/* Un-skip the density tasks involved with this cell. */
for (struct link *l = c->density; l != NULL; l = l->next) {
@@ -1066,7 +1066,9 @@ void cell_drift(struct cell *c, const struct engine *e) {
/* Drift from the last time the cell was drifted to the current time */
const double dt = (ti_current - ti_old) * timeBase;
- float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
+ float dx_max = 0.f, dx2_max = 0.f;
+ float dx_max_sort = 0.f, dx2_max_sort = 0.f;
+ float h_max = 0.f;
/* Check that we are actually going to move forward. */
if (ti_current < ti_old) error("Attempt to drift to the past");
@@ -1080,6 +1082,7 @@ void cell_drift(struct cell *c, const struct engine *e) {
struct cell *cp = c->progeny[k];
cell_drift(cp, e);
dx_max = max(dx_max, cp->dx_max);
+ dx_max_sort = max(dx_max_sort, cp->dx_max_sort);
h_max = max(h_max, cp->h_max);
}
@@ -1118,6 +1121,10 @@ void cell_drift(struct cell *c, const struct engine *e) {
xp->x_diff[1] * xp->x_diff[1] +
xp->x_diff[2] * xp->x_diff[2];
dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
+ const float dx2_sort = xp->x_diff_sort[0] * xp->x_diff_sort[0] +
+ xp->x_diff_sort[1] * xp->x_diff_sort[1] +
+ xp->x_diff_sort[2] * xp->x_diff_sort[2];
+ dx2_max_sort = (dx2_max_sort > dx2_sort) ? dx2_max_sort : dx2_sort;
/* Maximal smoothing length */
h_max = (h_max > p->h) ? h_max : p->h;
@@ -1125,16 +1132,19 @@ void cell_drift(struct cell *c, const struct engine *e) {
/* Now, get the maximal particle motion from its square */
dx_max = sqrtf(dx2_max);
+ dx_max_sort = sqrt(dx2_max_sort);
} else {
h_max = c->h_max;
dx_max = c->dx_max;
+ dx_max_sort = c->dx_max_sort;
}
/* Store the values */
c->h_max = h_max;
c->dx_max = dx_max;
+ c->dx_max_sort = dx_max_sort;
/* Update the time of the last drift */
c->ti_old = ti_current;
diff --git a/src/cell.h b/src/cell.h
index d99019f58ea20ad83c835178ba71569192be85a1..a89353a2b9105387fdc5f22b130ac1e108498303 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -232,6 +232,9 @@ struct cell {
/*! Maximum particle movement in this cell since last construction. */
float dx_max;
+ /*! Maximum particle movement in this cell since the last sort. */
+ float dx_max_sort;
+
/*! Nr of #part in this cell. */
int count;
diff --git a/src/drift.h b/src/drift.h
index 5ad0a8ad26d431f03dc5332d79df2ff6c17dead8..c87d6b1c0ab3a982203af5e175fed4ee19a147cf 100644
--- a/src/drift.h
+++ b/src/drift.h
@@ -90,9 +90,11 @@ __attribute__((always_inline)) INLINE static void drift_part(
hydro_predict_extra(p, xp, dt);
/* Compute offset since last cell construction */
- xp->x_diff[0] -= xp->v_full[0] * dt;
- xp->x_diff[1] -= xp->v_full[1] * dt;
- xp->x_diff[2] -= xp->v_full[2] * dt;
+ for (int k = 0; k < 3; k++) {
+ const float dx = xp->v_full[k] * dt;
+ xp->x_diff[k] -= dx;
+ xp->x_diff_sort[k] -= dx;
+ }
}
#endif /* SWIFT_DRIFT_H */
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 69ae79666e1db4e4f405c653cfc533606989a73a..571aaf39ed2509d95e9f68c34ab8e6c09ded3fbb 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -39,6 +39,9 @@ struct xpart {
/* Offset between current position and position at last tree rebuild. */
float x_diff[3];
+ /* Offset between the current position and position at the last sort. */
+ float x_diff_sort[3];
+
/* Velocity at the last full step. */
float v_full[3];
diff --git a/src/runner.c b/src/runner.c
index 5753973722f8ecb5c9299e5e2c9eb25cb8531926..167ef3969e531c80075537b89bf69998214c6beb 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -334,11 +334,11 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
struct entry *finger;
struct entry *fingers[8];
struct part *parts = c->parts;
+ struct xpart *xparts = c->xparts;
struct entry *sort;
- int j, k, count = c->count;
- int i, ind, off[8], inds[8], temp_i, missing;
+ const int count = c->count;
+ int off[8], inds[8], temp_i, missing;
float buff[8];
- double px[3];
TIMER_TIC
@@ -360,27 +360,28 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
if (c->split) {
/* Fill in the gaps within the progeny. */
- for (k = 0; k < 8; k++) {
+ for (int k = 0; k < 8; k++) {
if (c->progeny[k] == NULL) continue;
missing = flags & ~c->progeny[k]->sorted;
if (missing) runner_do_sort(r, c->progeny[k], missing, 0);
}
/* Loop over the 13 different sort arrays. */
- for (j = 0; j < 13; j++) {
+ for (int j = 0; j < 13; j++) {
/* Has this sort array been flagged? */
if (!(flags & (1 << j))) continue;
/* Init the particle index offsets. */
- for (off[0] = 0, k = 1; k < 8; k++)
+ off[0] = 0;
+ for (int k = 1; k < 8; k++)
if (c->progeny[k - 1] != NULL)
off[k] = off[k - 1] + c->progeny[k - 1]->count;
else
off[k] = off[k - 1];
/* Init the entries and indices. */
- for (k = 0; k < 8; k++) {
+ for (int k = 0; k < 8; k++) {
inds[k] = k;
if (c->progeny[k] != NULL && c->progeny[k]->count > 0) {
fingers[k] = &c->progeny[k]->sort[j * (c->progeny[k]->count + 1)];
@@ -391,8 +392,8 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
}
/* Sort the buffer. */
- for (i = 0; i < 7; i++)
- for (k = i + 1; k < 8; k++)
+ for (int i = 0; i < 7; i++)
+ for (int k = i + 1; k < 8; k++)
if (buff[inds[k]] < buff[inds[i]]) {
temp_i = inds[i];
inds[i] = inds[k];
@@ -401,7 +402,7 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
/* For each entry in the new sort list. */
finger = &sort[j * (count + 1)];
- for (ind = 0; ind < count; ind++) {
+ for (int ind = 0; ind < count; ind++) {
/* Copy the minimum into the new sort array. */
finger[ind].d = buff[inds[0]];
@@ -412,7 +413,7 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
buff[inds[0]] = fingers[inds[0]]->d;
/* Find the smallest entry. */
- for (k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
+ for (int k = 1; k < 8 && buff[inds[k]] < buff[inds[k - 1]]; k++) {
temp_i = inds[k - 1];
inds[k - 1] = inds[k];
inds[k] = temp_i;
@@ -435,11 +436,11 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
else {
/* Fill the sort array. */
- for (k = 0; k < count; k++) {
- px[0] = parts[k].x[0];
- px[1] = parts[k].x[1];
- px[2] = parts[k].x[2];
- for (j = 0; j < 13; j++)
+ for (int k = 0; k < count; k++) {
+ xparts[k].x_diff_sort[0] = xparts[k].x_diff_sort[1] =
+ xparts[k].x_diff_sort[2] = 0.0f;
+ const double px[3] = {parts[k].x[0], parts[k].x[1], parts[k].x[2]};
+ for (int j = 0; j < 13; j++)
if (flags & (1 << j)) {
sort[j * (count + 1) + k].i = k;
sort[j * (count + 1) + k].d = px[0] * runner_shift[j][0] +
@@ -449,7 +450,7 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
}
/* Add the sentinel and sort. */
- for (j = 0; j < 13; j++)
+ for (int j = 0; j < 13; j++)
if (flags & (1 << j)) {
sort[j * (count + 1) + count].d = FLT_MAX;
sort[j * (count + 1) + count].i = 0;
@@ -460,10 +461,10 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
#ifdef SWIFT_DEBUG_CHECKS
/* Verify the sorting. */
- for (j = 0; j < 13; j++) {
+ for (int j = 0; j < 13; j++) {
if (!(flags & (1 << j))) continue;
finger = &sort[j * (count + 1)];
- for (k = 1; k < count; k++) {
+ for (int k = 1; k < count; k++) {
if (finger[k].d < finger[k - 1].d)
error("Sorting failed, ascending array.");
if (finger[k].i >= count) error("Sorting failed, indices borked.");
@@ -1056,9 +1057,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
/* What is the next sync-point ? */
ti_end_min = min(ti_current + ti_new_step, ti_end_min);
ti_end_max = max(ti_current + ti_new_step, ti_end_max);
- }
-
- else { /* part is inactive */
+ } else { /* part is inactive */
const integertime_t ti_end =
get_integer_time_end(ti_current, p->time_bin);
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index e78d87e4b803b441f86fdd4db4fdcb0b9f61e091..8ab099f268c85d05db70b2af11086aac02bdaf7a 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -461,7 +461,7 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
/* Pick-out the sorted lists. */
const struct entry *restrict sort_j = &cj->sort[sid * (cj->count + 1)];
- const float dxj = cj->dx_max;
+ const float dxj = cj->dx_max_sort;
/* Parts are on the left? */
if (!flipped) {
@@ -749,7 +749,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
struct part *restrict parts_j = cj->parts;
const double di_max = sort_i[count_i - 1].d - rshift;
const double dj_min = sort_j[0].d;
- const float dx_max = (ci->dx_max + cj->dx_max);
+ const float dx_max = (ci->dx_max_sort + cj->dx_max_sort);
/* Loop over the parts in ci. */
for (int pid = count_i - 1;
@@ -960,7 +960,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
struct part *restrict parts_j = cj->parts;
const double di_max = sort_i[count_i - 1].d - rshift;
const double dj_min = sort_j[0].d;
- const double dx_max = (ci->dx_max + cj->dx_max);
+ const double dx_max = (ci->dx_max_sort + cj->dx_max_sort);
/* Collect the number of parts left and right below dt. */
int countdt_i = 0, countdt_j = 0;
@@ -1832,7 +1832,8 @@ void DOSUB_PAIR1(struct runner *r, struct cell *ci, struct cell *cj, int sid,
/* Recurse? */
if (ci->split && cj->split &&
- max(ci->h_max, cj->h_max) * kernel_gamma + ci->dx_max + cj->dx_max <
+ max(ci->h_max, cj->h_max) * kernel_gamma + ci->dx_max_sort +
+ cj->dx_max_sort <
h / 2) {
/* Different types of flags. */
@@ -2118,7 +2119,8 @@ void DOSUB_PAIR2(struct runner *r, struct cell *ci, struct cell *cj, int sid,
/* Recurse? */
if (ci->split && cj->split &&
- max(ci->h_max, cj->h_max) * kernel_gamma + ci->dx_max + cj->dx_max <
+ max(ci->h_max, cj->h_max) * kernel_gamma + ci->dx_max_sort +
+ cj->dx_max_sort <
h / 2) {
/* Different types of flags. */
@@ -2424,7 +2426,8 @@ void DOSUB_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
/* Recurse? */
if (ci->split && cj->split &&
- max(ci->h_max, cj->h_max) * kernel_gamma + ci->dx_max + cj->dx_max <
+ max(ci->h_max, cj->h_max) * kernel_gamma + ci->dx_max_sort +
+ cj->dx_max_sort <
h / 2) {
/* Get the type of pair if not specified explicitly. */
diff --git a/src/space.c b/src/space.c
index 24c2c413198d32e0942a0d26e5a3efbd5571c68e..5aff6bc17265a2c058366804fe8a1adfa0d898df 100644
--- a/src/space.c
+++ b/src/space.c
@@ -207,6 +207,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->force = NULL;
c->grav = NULL;
c->dx_max = 0.0f;
+ c->dx_max_sort = 0.0f;
c->sorted = 0;
c->count = 0;
c->gcount = 0;
@@ -1551,6 +1552,7 @@ void space_split_recursive(struct space *s, struct cell *c,
cp->split = 0;
cp->h_max = 0.0;
cp->dx_max = 0.f;
+ cp->dx_max_sort = 0.f;
cp->nodeID = c->nodeID;
cp->parent = c;
cp->super = NULL;