Reversed initial conditions leads to large errors in Fourier space

[svstopyra]

I seem to have come across a quite significant problem with the reverse function.

When comparing the Fourier space overdensity fields in both the reversed and non-reversed cases, the fields should add to zero in Fourier space (and real space), but currently it appear that they do not. The errors don't seem to be small, either, but appear only on certain slices of the Fourier space grid (with most slices adding to exactly zero as expected).

The problem seems to occur with the attached parameter file for example (comparing with and without the "reverse" line).

paramfile_unreversed.txt

The attached image shows the sum of the density fields for two different slices (slice 0 and slice 255 if you import them directly into numpy). Most of them are like the top slice (0) - all exactly zero everywhere as expected. Some look more like the bottom slice (255).

There doesn't seem to be anything in the reverse code (multilevelfield.hpp, line 229) that would cause this, as it's simply flipping the sign. I suspect that it could be an error related to performing the Fourier transform (because the code does this back and forth between real and Fourier space in few places).

In real space, the fields don't appear to be reversals of each other at all as the differences lead to large fluctuations, and this seems to affect the resulting particle offsets in the Zeldovich approximation quite significantly.

[apontzen]

Thanks -- very interesting! I wonder when this problem was introduced... could you do some more digging to see why precisely it's happening? Is there anything strange looking about the density field itself in the problematic slices?

[svstopyra]

Will do - I've narrowed down its appearance to some time in April 2017 (before that the sum is zero everywhere, afterwards the problem arises). There seem to be a lot of commits that month so it might take me a while to bisect through them. It's hard to say whether the density fields look strange, unfortunately, because they're in Fourier space which makes them hard to interpret (in real space it's less obvious where the problem is, as the error on a few levels in Fourier space is scattered throughout the grid in real space).

[apontzen]

Great, thanks for working through this. To see whether it is strange, perhaps just compare to the expected amplitude from the power spectrum. (This wouldn't capture every possible problem, but it might capture some!)

[svstopyra]

Ok, I've managed to find and fix the bug. As it turns out, the problem is annoyingly simple. Basically, the reverse assumed at all times that the density vector had N^3 elements in it. This was true in some early versions of the code, but at some point in April 2017, the way Fourier transforms were implemented was changed to more efficiently treat real fields. This means that the Fourier space field actually stores the only the unique real/imaginary parts, rather than N^3 complex numbers, to save space. As it turns out, that is always more than N^3 numbers, so reverse has been flipping the first N^3, but not the rest ever since 2017. The fix is simply to always loop up to the size of the data storage vector, not the size of the grid, when flipping the signs of the density contrast. I'll submit a pull request with the updated code once I've written a new test for this.

[apontzen]

Thanks for tracking down and fixing this