Issue with the way `cell_hydro.dx_max_sort` is computed
I want to preface this with saying that I'm pretty sure, this is a real issue, but that I'm also aware that this really is part of the core of SWIFT and has been looked at and tested extensively, so I might be wrong after all...
Description of the issue:
Because of the way cell_hydro.dx_max_sort
is currently computed (see lines 323-326 and lines 365, 371 in cell_drift.c
and 211-213 in drift.h
), it is not monotonically increasing.
When some sorting directions are computed later than others, cell_hydro.dx_max_sort
is no longer guaranteed to be an upper bound of the individual particles' particle shift differences for the sorted sid
's (as checked in e.g. the hydro dopair functions).
This is because cell_hydro.dx_max_sort
is computed from hydro_part.x_diff_sort
which stores the vector offset from the position of the particle when the first sorts where computed (since the last cleanup of the sorts) of the cell containing the particle. Suppose that for some sid
the sorts are updated at a later timestep then the last cleanup (where hydro_part.x_diff_sort
has been reset) and subsequently, the particles move back to a position closer to their positions at cleanup. Then cell_hydro.dx_max_sort
can decrease and actually become smaller than some of the particles' shift differences for the sid
that was sorted later (when the length of hydro_part.x_diff_sort
was larger). This violates one of the core assumptions when using the sorting arrays in a neighbor loop.
Expected behavior:
At any time, for all the particles, the particle's shift difference for any currently sorted 'sid' (indicated by cell_hydro.sorted
) should be smaller than cell_hydro.dx_max_sort
of the cell(s) containing the particle.
This is implicitly assumed when checking whether resorting is needed and explicitly checked in e.g. the the hydro dopair functions.
Observed behavior:
For some specific task structures, some directions will not be sorted immediately after a cleanup of the sorts. For some particular movement of the particles, this can cause cell_hydro.dx_max_sort
to decrease after those directions are sorted until it becomes smaller than some particles' shift differences for those directions (mechanism described above). I caught this behavior in the moving mesh branch with exaggerated steering causing some oscillatory motions in some of the particles.
Why this might have gone unnoticed:
I think this might have gone unnoticed for several reasons:
- It only occurs for specific task structures and specific particle motions
- It is only enforced when debug checks are activated
- The result might not even be wrong if
cell_hydro.dx_max_sort
is too small. This just means that we might miss some particle interactions for some SID's - Even if we do miss a few particle interactions, that will still be hard to see from the results alone
Proposed fix:
Making this change to the computation of cell_hydro.dx_max_sort
fixes the issue for me, but does increase the number of times particles are sorted a bit.
I'm curious if anybody has other suggestions.