It turns out that the decision on when to box-wrap a particle's position is not always consistent with the cell index returned by cell_getid if a particle is very close to the edge of the box. The following code from engine_redistribute.c:
for (int k = 0; k < num_elements; k++) { \
for (int j = 0; j < 3; j++) { \
if (parts[k].x[j] < 0.0) \
parts[k].x[j] += s->dim[j]; \
else if (parts[k].x[j] >= s->dim[j]) \
parts[k].x[j] -= s->dim[j]; \
} \
const int cid = cell_getid(s->cdim, parts[k].x[0] * s->iwidth[0], \
parts[k].x[1] * s->iwidth[1], \
parts[k].x[2] * s->iwidth[2]); \Will sometimes pass on indices to cell_getid that are equal to the corresponding entry in s->cdim. A similar block of code can be found (multiple times) in space_cell_index.c. I think this makes sense: there is no reason why a condition involving x[j] and s->dim[j] should exactly match the result of x[j] * s->iwidth[j].
As a result of this inconsistency, particles can get moved to the wrong node, triggering all kinds of unexpected bugs.
In this MR, I try to address this issue by making the decision on when to box-wrap in integer space, by pre-computing the arguments to cell_getid. However, the current version still triggers some issues, so I am clearly missing something.
An alternative solution would be to make sure that all arguments to cell_getid are limited to the range [0,s->cdim[j]-1]. But that cannot be done with the current macro.
@matthieu any thoughts on this? I will try some more things tomorrow.