I have run the code and looked at some results. One thing I found puzzling was the spectra. For example if I set Nz to 48 I would expect 24 non zero values in the spectra. Instead I find 16 which tells me this is probability an issue on aliasing.
I have run a channel flow code written at the Naval Research Laboratory in the past that is based on the original channel flow code described in the 1986 paper by Kim Moin and Moser. This code handled the aliasing in the x-z plane one plane at a time. For example if the resolution was Nx = 64, Ny = 129, and Nz = 48, it would loop through each of the 129 planes and create a fine grid of 96 x 72 (Multiply Nx and Nz by 3/2) and calculate the non linear terms in this plane, then transform back to spectral space and only store 32 x 48 complex terms in the standard array.
From the numbers I am seeing, I get the feeling that you do not do this but keep all terms in the fine grid which is more costly in memory but memory is cheap enough. The concern I have is understanding exactly the grid I am using. Kim Moin and Moser did their 1986 run on a grid of 192 x 129 x 160 in X, Y and Z respectively. One can look at the spectra on page 137 of their paper to see this relates to the number of non-zero values in the spectra. If I am making the right assumptions about ChannelFlow 2.0, I would need to run on a grid of 288 x 129 x 240 to do the exact same simulation. Is that correct or am I misunderstanding something ?
This would also change grid spacing needed to fully resolve the flow. They (KMM) suggest keeping delta x at 12 viscous units and delta z at 7 viscous units. With your code, I would need to keep them at 8 and 4.66 to fully resolve the flow. I could also operate at a higher CFL number and I see you are happy operating at 0.5 where I have always used the rule of thumb of 0.3 . If my understanding is correct 0.3 in the KMM grid would correspond to 0.45 on your grid.
Any insight people could offer would be greatly appreciated. Any way you do it is fine. Either way would take the same compute time but their way would save on memory but I am not stressed for memory. I just want to be able to understand what I am doing and the fields I have.