[BUG] Periodic boundary conditions bug: Fields with RSDs have kinks in the edges when using `periodic=True`

[nikos-triantafyllou]

Bug description:
At https://github.com/21cmfast/21cmFAST/blob/main/src/py21cmfast/rsds.py#L106

Fields with RSDs have kinks in the edges when using apply_rsds with periodic=True.

Proof:
This is demonstrated in the following images.

The orange line is the mean of the Tb field of a coeval cube for different redshifts (shown here in pixels). Blue line is the same after using the function apply_rsds with periodic=True. There is an apparent unnatural kink close to the 0th pixel.

For reference, here I plot the same thing, but now with green it's the result using the same function apply_rsds, with periodic=False but manually adding 2 copies of the original box at higher and lower z (i.e,. concatenating the same box at the edges to apply periodic boundary conditions).

To reproduce:

Compare:

lc_with_rsds = apply_rsds(
    field = lc,
    los_velocity = lc_velocity_z,
    redshifts= lc_redshifts,
    inputs_dict =  inputs_dict,
    periodic =  True,
    n_rsd_subcells = 1,
)

with

pad=lc.shape[0]
inputs_dict['H_at_all_z'] = Planck18.H(padded_lc_redshifts)

padded_lc = np.concatenate([lc[:, :, -pad:], 
                              lc, 
                              lc[:, :, :pad]], axis=2)

padded_lc_velocity_z = np.concatenate([lc_velocity_z[:, :, -pad:], 
                              lc_velocity_z, 
                              lc_velocity_z[:, :, :pad]], axis=2)

lc_with_rsds_padded = RSDs.apply_rsds(
    field = padded_lc,
    los_velocity = padded_lc_velocity_z,
    redshifts= padded_lc_redshifts,
    inputs_dict =  inputs_dict,
    periodic =  False,
    n_rsd_subcells = 1,
)
lc_with_rsds_correct = lc_with_rsds_padded[:,:,pad:-pad]

[jordanflitter]

Thanks @nikos-triantafyllou for reporting this and providing useful information! I will try to look at this asap.
For reproducibility, could you share your script/notebook? Can you reproduce this kink with a small box? (small HII_DIM)

I'm confused a bit for the inputs_dict argument in your code. Why is this a dictionary and not an object of type InputParameters?

It seems to me that if there's a problem it lies in this block. It seems like we only put the first (last) slice of the velocity array at the end (beginning) when periodic=True, I honestly can't remember why we do it, the wrapping of the masses should be done in cosmotile (this line here). Any thoughts on this @steven-murray?

[jordanflitter]

I think we do what we do in that block is in order not to screw up the interpolation of the velocity field at the subcells at the edges. We do linear interpolation so these subcells must have the information from the edges of the box

[nikos-triantafyllou]

I think I've rewritten it to work with dictionaries so that there are as few dependencies as possible. Unfortunately, my code is specifically applied to the database, so it would be way more confusing if I sent my script. Lmk if you still want that.
But it looks like you found the issue.

[jordanflitter]

No, I haven't found the issue yet. It is going to be difficult to debug this without the inputs for apply_rsds. Can you share a pickle file containing those inputs? (with inputs of type InputParameters)

[steven-murray]

Oof, yeah I remember the vague reasons we did this, but not the specific ones that would help with this. Not sure I'm going to have time to look at this very soon.