Fix 2/3 rule filter

[Pibe97]

Currently, the 2/3 filter retains more negative than positive frequencies in the kx direction. For example, working on a 128 by 128 grid (i.e. 128 by 65 in rfft2 space), in the kx direction (i.e. rows) the current filter retains the kx = 0 and 41 positive frequencies, but 43 negative frequencies. This error can be spotted by checking that applying jnp.fft.irfft2 followed by jnp.fft.rfft2 to a trajectory output, which are inverse operations, changes the output. Alternatively, the 0th column (ky=0) of a trajectory snapshot (in rfft2 space), should have complex conjugate symmetry between the positive and negative kx frequencies.

The code has been adjusted to retain the kx = 0 frequency (hence the + 1 in line 132), int(2 / 3 * n) // 2 positive and int(2 / 3 * n) // 2 negative frequencies (fixed by adding parentheses in line 133, to have correct order of operations) in the kx direction. For the 128 by 128 grid, we have int(2 / 3 * n) // 2 = 42.

[gideonite]

Interesting, thanks for the PR. Looks like there might be an off-by-one error here, is that right? Does this work for even and odd number of grid points?

checking that applying jnp.fft.irfft2 followed by jnp.fft.rfft2 to a trajectory output, which are inverse operations, changes the output

I'm not sure I completely follow here. Isn't filtering not invertible?

I'm not clear on this line. This looks like the same computation but with different parenthesis?

I don't want to set the bar too high for this, but am curious, has this affected any experiments or results that you have been working on?

[Pibe97]

Hi Gideon,

Thanks for your reply and questions!

Looks like there might be an off-by-one error here, is that right? Does this work for even and odd number of grid points?

There were actually two different frequencies, at least in the case of n = 128. The 41 first positive frequencies were retained and the 43 first negative ones. The fix should work for both even and odd number of grid points, as what's important is just that the index in line 132 goes one step further compared to line 133, to account for the 0 frequency at the 0th entry, which should always be the case.

I'm not sure I completely follow here. Isn't filtering not invertible?

Apologies for being unclear. Filtering is indeed not invertible. What I meant is applying these operations to the filtered output "trajectory". With the old filtered "trajectory", I found that jnp.fft.rfft2(jnp.fft.irfft2(trajectory)) != trajectory, which should not be the case, since rfft2 and irfft2 are inverses of each other. This happens because the number of positive frequencies and negative frequencies after the old filtered trajectory do not match. With the suggested fix I found that there's equality.

I'm not clear on this line. This looks like the same computation but with different parenthesis?

The parentheses seem to matter as the expressions give different outputs (at least for me). E.g. for n =128, -(int(2 / 3 * n) // 2) gives -42, while -int(2 / 3 * n) // 2 gives -43.

I don't want to set the bar too high for this, but am curious, has this affected any experiments or results that you have been working on?

I haven't re-run all my computations yet, but I doubt there will be a huge difference. This will also creep into temporal derivatives though, where the high-frequency terms matter more due to the Laplacian in Navier-Stokes. I ran a quick test and it seems that the average absolute difference between temporal derivatives with the old and new filters is only $O(10^{-2})$. However this was at a relatively low Reynolds number, maybe for larger Re it would matter more.

Cheers!