Skip to content

Fix inconsistent box wrapping for particles on edge of box

Bert Vandenbroucke requested to merge safe-box-wrapping into master

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.

Edited by Bert Vandenbroucke

Merge request reports