Commit 1775f39e authored by Pedro Gonnet's avatar Pedro Gonnet
Browse files

track movement since sorting and rebuild separately in cells.

parent 473d5f9a
......@@ -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;
......
......@@ -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;
......
......@@ -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 */
......@@ -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];
......
......@@ -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);
......
......@@ -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. */
......
......@@ -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;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment