Skip to content
GitLab
Commit 73d4e5db authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Sort particles consistently between edges.

parent 52861646
Hide whitespace changes
Inline Side-by-side
src/cell.c
View file @ 73d4e5db
......@@ -803,8 +803,8 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
/* Fill the buffer with the indices. */
for (int k = 0; k < count; k++) {
const int bid = (buff[k].x[0] > pivot[0]) * 4 +
(buff[k].x[1] > pivot[1]) * 2 + (buff[k].x[2] > pivot[2]);
const int bid = (buff[k].x[0] >= pivot[0]) * 4 +
(buff[k].x[1] >= pivot[1]) * 2 + (buff[k].x[2] >= pivot[2]);
bucket_count[bid]++;
buff[k].ind = bid;
}
......@@ -876,44 +876,44 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
/* Verify a few sub-cells. */
for (int k = 0; k < c->progeny[0]->count; k++)
if (c->progeny[0]->parts[k].x[0] > pivot[0] ||
c->progeny[0]->parts[k].x[1] > pivot[1] ||
c->progeny[0]->parts[k].x[2] > pivot[2])
if (c->progeny[0]->parts[k].x[0] >= pivot[0] ||
c->progeny[0]->parts[k].x[1] >= pivot[1] ||
c->progeny[0]->parts[k].x[2] >= pivot[2])
error("Sorting failed (progeny=0).");
for (int k = 0; k < c->progeny[1]->count; k++)
if (c->progeny[1]->parts[k].x[0] > pivot[0] ||
c->progeny[1]->parts[k].x[1] > pivot[1] ||
c->progeny[1]->parts[k].x[2] <= pivot[2])
if (c->progeny[1]->parts[k].x[0] >= pivot[0] ||
c->progeny[1]->parts[k].x[1] >= pivot[1] ||
c->progeny[1]->parts[k].x[2] < pivot[2])
error("Sorting failed (progeny=1).");
for (int k = 0; k < c->progeny[2]->count; k++)
if (c->progeny[2]->parts[k].x[0] > pivot[0] ||
c->progeny[2]->parts[k].x[1] <= pivot[1] ||
c->progeny[2]->parts[k].x[2] > pivot[2])
if (c->progeny[2]->parts[k].x[0] >= pivot[0] ||
c->progeny[2]->parts[k].x[1] < pivot[1] ||
c->progeny[2]->parts[k].x[2] >= pivot[2])
error("Sorting failed (progeny=2).");
for (int k = 0; k < c->progeny[3]->count; k++)
if (c->progeny[3]->parts[k].x[0] > pivot[0] ||
c->progeny[3]->parts[k].x[1] <= pivot[1] ||
c->progeny[3]->parts[k].x[2] <= pivot[2])
if (c->progeny[3]->parts[k].x[0] >= pivot[0] ||
c->progeny[3]->parts[k].x[1] < pivot[1] ||
c->progeny[3]->parts[k].x[2] < pivot[2])
error("Sorting failed (progeny=3).");
for (int k = 0; k < c->progeny[4]->count; k++)
if (c->progeny[4]->parts[k].x[0] <= pivot[0] ||
c->progeny[4]->parts[k].x[1] > pivot[1] ||
c->progeny[4]->parts[k].x[2] > pivot[2])
if (c->progeny[4]->parts[k].x[0] < pivot[0] ||
c->progeny[4]->parts[k].x[1] >= pivot[1] ||
c->progeny[4]->parts[k].x[2] >= pivot[2])
error("Sorting failed (progeny=4).");
for (int k = 0; k < c->progeny[5]->count; k++)
if (c->progeny[5]->parts[k].x[0] <= pivot[0] ||
c->progeny[5]->parts[k].x[1] > pivot[1] ||
c->progeny[5]->parts[k].x[2] <= pivot[2])
if (c->progeny[5]->parts[k].x[0] < pivot[0] ||
c->progeny[5]->parts[k].x[1] >= pivot[1] ||
c->progeny[5]->parts[k].x[2] < pivot[2])
error("Sorting failed (progeny=5).");
for (int k = 0; k < c->progeny[6]->count; k++)
if (c->progeny[6]->parts[k].x[0] <= pivot[0] ||
c->progeny[6]->parts[k].x[1] <= pivot[1] ||
c->progeny[6]->parts[k].x[2] > pivot[2])
if (c->progeny[6]->parts[k].x[0] < pivot[0] ||
c->progeny[6]->parts[k].x[1] < pivot[1] ||
c->progeny[6]->parts[k].x[2] >= pivot[2])
error("Sorting failed (progeny=6).");
for (int k = 0; k < c->progeny[7]->count; k++)
if (c->progeny[7]->parts[k].x[0] <= pivot[0] ||
c->progeny[7]->parts[k].x[1] <= pivot[1] ||
c->progeny[7]->parts[k].x[2] <= pivot[2])
if (c->progeny[7]->parts[k].x[0] < pivot[0] ||
c->progeny[7]->parts[k].x[1] < pivot[1] ||
c->progeny[7]->parts[k].x[2] < pivot[2])
error("Sorting failed (progeny=7).");
#endif
......
src/debug.c
View file @ 73d4e5db
......@@ -224,10 +224,11 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
struct part *const p = &parts[k];
struct xpart *const xp = &xparts[k];
if (p->x[0] < loc_min[0] || p->x[0] > loc_max[0] || p->x[1] < loc_min[1] ||
p->x[1] > loc_max[1] || p->x[2] < loc_min[2] || p->x[2] > loc_max[2]) {
if (p->x[0] < loc_min[0] || p->x[0] >= loc_max[0] ||
p->x[1] < loc_min[1] || p->x[1] >= loc_max[1] ||
p->x[2] < loc_min[2] || p->x[2] >= loc_max[2]) {
message(
error(
"Inconsistent part position p->x=[%e %e %e], c->loc=[%e %e %e] "
"c->width=[%e %e %e]",
p->x[0], p->x[1], p->x[2], c->loc[0], c->loc[1], c->loc[2],
......
src/periodic.h
View file @ 73d4e5db
......@@ -36,7 +36,7 @@
const __typeof__(x) _x = (x); \
const __typeof__(a) _a = (a); \
const __typeof__(b) _b = (b); \
_x < _a ? (_x + _b) : ((_x > _b) ? (_x - _b) : _x); \
_x < _a ? (_x + _b) : ((_x >= _b) ? (_x - _b) : _x);\
})
/**
......
src/space.c
View file @ 73d4e5db
......@@ -1081,7 +1081,11 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts,
ind[k] = index;
#ifdef SWIFT_DEBUG_CHECKS
if (pos_x > dim_x || pos_y > dim_y || pos_z > pos_z || pos_x < 0. ||
if(index < 0 || index > cdim[0] * cdim[1] * cdim[2])
error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index,
cdim[0], cdim[1], cdim[2], pos_x, pos_y, pos_z);
if (pos_x >= dim_x || pos_y >= dim_y || pos_z >= dim_z || pos_x < 0. ||
pos_y < 0. || pos_z < 0.)
error("Particle outside of simulation box. p->x=[%e %e %e]", pos_x, pos_y,
pos_z);
......@@ -1138,6 +1142,17 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts,
cell_getid(cdim, pos_x * ih_x, pos_y * ih_y, pos_z * ih_z);
ind[k] = index;
#ifdef SWIFT_DEBUG_CHECKS
if(index < 0 || index > cdim[0] * cdim[1] * cdim[2])
error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index,
cdim[0], cdim[1], cdim[2], pos_x, pos_y, pos_z);
if (pos_x >= dim_x || pos_y >= dim_y || pos_z >= dim_z || pos_x < 0. ||
pos_y < 0. || pos_z < 0.)
error("Particle outside of simulation box. p->x=[%e %e %e]", pos_x, pos_y,
pos_z);
#endif
/* Update the position */
gp->x[0] = pos_x;
gp->x[1] = pos_y;
......@@ -1189,6 +1204,17 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts,
cell_getid(cdim, pos_x * ih_x, pos_y * ih_y, pos_z * ih_z);
ind[k] = index;
#ifdef SWIFT_DEBUG_CHECKS
if(index < 0 || index > cdim[0] * cdim[1] * cdim[2])
error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index,
cdim[0], cdim[1], cdim[2], pos_x, pos_y, pos_z);
if (pos_x >= dim_x || pos_y >= dim_y || pos_z >= dim_z || pos_x < 0. ||
pos_y < 0. || pos_z < 0.)
error("Particle outside of simulation box. p->x=[%e %e %e]", pos_x, pos_y,
pos_z);
#endif
/* Update the position */
sp->x[0] = pos_x;
sp->x[1] = pos_y;
......@@ -2609,11 +2635,15 @@ void space_first_init_parts(struct space *s) {
p[i].v[1] = p[i].v[2] = 0.f;
#endif
double temp = p[i].v[0];
p[i].v[0] = p[i].v[1];
p[i].v[1] = temp;
hydro_first_init_part(&p[i], &xp[i]);
#ifdef SWIFT_DEBUG_CHECKS
p->ti_drift = 0;
p->ti_kick = 0;
p[i].ti_drift = 0;
p[i].ti_kick = 0;
#endif
}
}
......
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