Added spectrum computation to 3D RBC

[brownbaerchen]

In turbulent flow simulations people often show a sort of energy spectrum to demonstrate that the spatial resolution is adequate. The main feature of this PR is a method for computing precisely this spectrum. @tlunet, please check if you find this reasonable.
Also, I fixed a bug in the Nusselt number computation. The formula employed there previously was valid for the special case of Ra=1, which is a case nobody cares about.
Now I can generate plots like this for verification:

Some background on the spectrum

People are often worried about aliasing. As an example, consider the grid $x= 0, \pi/2, \pi, 3\pi / 2$. Evaluating the wave $\exp(i 2 x)$ on this grid, gives $1, -1, 1, -1$. Now consider $\exp(i6 \pi)$. Evaluating this on the grid gives the same values $1, -1, 1, -1$ because of the periodicity of the grid. The wave with wavenumber 6 is aliased to the wave with wavenumber 2 on this grid with only four modes.
The fastest wavenumber supported on this grid is 3. This is easy to see because wavenumber 4 gets aliased to wavenumber 0 as both waves are always evaluated as constant 1.

Aliasing is not much of a problem as long as the wavenumbers greater than supported by the resolution never appear. However, they do tend to appear in nonlinear problems. In the above example, the fastest supported wave is wavenumber 3 and we just saw that it's square, which has wavenumber 6 is aliased to wavenumber 2.
When we are dealing with convection, we have terms such as $u\partial_x u$, which introduce waves of twice the speed, which may get aliased.

Now spectral methods are famous for exponential decay of coefficients with wavenumber. This is only true for special cases, but as a rule of thumb, coefficients decay quickly if the problem is well resolved.
If indeed coefficients decay quickly enough, the aliasing will not matter so much because the spurious contributions in low wavenumbers due to aliasing are small compared to the magnitude of the non-spurious slow waves.
This is why turbulence people show these spectrum plots: A lot of energy in the fastest supported wave, often called "2h wave" is associated with instability due to a lot of aliasing, but if there is little energy in the 2h wave, it's fine.

Actually, we get rid of the instability caused by aliasing by doing dealiasing.
The idea is simple: Get a larger grid that supports the waves so they don't get aliased. But larger grid is also more expensive and you may run into a vicious circle of increasing the resolution. So you only increase the size of the grid in physical space.
Before the inverse Fourier transform from frequency space to grid space, you padd the solution with zeros in wave numbers that would be aliased. Then you compute the square of the waves on the grid, compute the Fourier transform to go back to frequency space and throw away all the waves of wave numbers higher than the original resolution.
Orzag figured out that you need to padd by 50% to get rid of aliasing.

So aliasing is not really a problem for us, but we should still show these plots so reviewers are happy and if the coefficients are large that is still a sign of poor resolution.
Don't mind all the above text. Thanks for rubber ducking!

[tlunet]

There is also a physical motivation behind the computation of spectrum, that does not consider aliasing 😅 ... but thanks for the quick novel on it, didn't realize the aliasing problem before 😉

[tlunet]

Could you add just a bit of documentation here ? One can understand from the function how it works, but it would be better to allow that directly from the docstring

[tlunet]

Seems fine, although I think that the spectrum computation should be put in a separate helper module. I do understand why you need to link it to the problem to have the transform available, but then a separate function should be able to compute it from the "pseudo-spectral" 3D field. Would allow to re-use the code for other problems or even other codes (like Dedalus with SDC, etc ...)

[brownbaerchen]

Seems fine, although I think that the spectrum computation should be put in a separate helper module. I do understand why you need to link it to the problem to have the transform available, but then a separate function should be able to compute it from the "pseudo-spectral" 3D field. Would allow to re-use the code for other problems or even other codes (like Dedalus with SDC, etc ...)

No, I think this should stay here because it is computing the spectrum in the horizontal plane for every point in z. This is specific to this RBC configuration.
Also, the distributed computation is a main part of this which requires local wave numbers and communicators, something that the problem class already has. And if you want to use this function with Dedalus data, best load the Dedalus data with this pySDC RBC class and then do the computation.

I'll add some documentation later.

[tlunet]

Looks good, thanks !