Reflective boundary conditions show undesirable behavior like discontinuities and asymmetry

[antoniusf]

Describe the bug While extending @LeanderSchlegel's work on tricubic interpolation to support reflective boundary conditions, we discovered some problems with the current implementation, such as:

• the interpolation "going backwards" in mirrored sections of the grid, leading to an incorrectly interpolated grid, along with discontinuities,
• the grid not being centered inside the used boundaries, meaning that points on opposite faces of the original grid being treated differently; some points occur three times in succession, others only once,
• the reduction function returning out-of-bounds indices, leading to memory access violations and garbled results.

For more details, including the formal analysis which resulted in the discovery of these bugs, please see the original issue in Leander's repository. We're opening this issue here to formally report the bug and get a fix merged upstream.

To Reproduce

from crpropa import Grid1f, Vector3d
import numpy as np
import matplotlib.pyplot as plt

grid = Grid1f(Vector3d(0), 3, 1.0)

grid.setValue(0, 0, 0, 1)
grid.setValue(0, 1, 0, 2)
grid.setValue(0, 2, 0, 0.5)

grid.setValue(1, 0, 0, 3)
grid.setValue(1, 1, 0, 4)
grid.setValue(1, 2, 0, 0)

grid.setReflective(True)

# this uses N = 901 to exacerbate the cases where we sample exactly
# on-grid, highlighting the memory access issues
N = 901
sampled_field = np.zeros((N, N))

for (ix, x) in enumerate(np.linspace(0.5, 9.5, N)):
    for (iy, y) in enumerate(np.linspace(0.5, 9.5, N)):
        sampled_field[ix, iy] = grid.interpolate(Vector3d(x, y, 0))

# uncomment vmin and vmax to observe artifacts resulting from oob access
# (note that, naturally, this is not deterministic and might not work reliably)
plt.imshow(sampled_field.T, origin="lower", vmin=0, vmax=2)
plt.colorbar()
plt.show()

Expected behavior Definitely not this:

System (please complete the following information): All – this is an algorithmic issue in Grid.h.

Additional context See our comments in the original issue.

I'm tagging in @LeanderSchlegel, @JulienDoerner, and @reichherzerp since they've been working on the problem over there.

[lukasmerten]

Should be closed with #363.