From 73d4e5db66fe4f387fb9448dc875b3e98e3baf9c Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Thu, 11 Jan 2018 18:03:35 +0100
Subject: [PATCH] Sort particles consistently between edges.
---
src/cell.c | 52 +++++++++++++++++++++++++-------------------------
src/debug.c | 7 ++++---
src/periodic.h | 2 +-
src/space.c | 36 +++++++++++++++++++++++++++++++---
4 files changed, 64 insertions(+), 33 deletions(-)
diff --git a/src/cell.c b/src/cell.c
index 488637599f..a36b8fd2fd 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -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
diff --git a/src/debug.c b/src/debug.c
index f63ed557d9..067ca156e7 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -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],
diff --git a/src/periodic.h b/src/periodic.h
index 5874b8742e..ba34534a4c 100644
--- a/src/periodic.h
+++ b/src/periodic.h
@@ -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);\
})
/**
diff --git a/src/space.c b/src/space.c
index 02f285edb7..8a07e2e32d 100644
--- a/src/space.c
+++ b/src/space.c
@@ -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
}
}
--
GitLab