Wrong spectrum calculation?

[parfenyev]

Hello everyone!

In an example the radial energy spectrum is calculated via (line 184):

E = @. 0.5 * (vars.u^2 + vars.v^2)
Eh = rfft(E)
kr, Ehr = FourierFlows.radialspectrum(Eh, grid, refinement=1)

Usually another quantity is called the energy spectrum -- the sum of the squares of the Fourier harmonics of the velocity fields should be taken. It can be also expressed as

@. Eh = 0.5 * grid.invKrsq * abs2(prob.sol)

I suggest correcting the example.

[navidcy]

So what we compute there is what's documented in the FourierFlows.radialspectrum docstring. It's the spectrum of u² which, I'm pretty sure, that with some manipulations it can be shown that is proportional to the spectrum of |û|².

Are you saying that this is not true? Or are you suggesting we could compute another, alternative, spectrum-related quantity?

[parfenyev]

In the natural sciences, people are usually interested in the quantity |fft(u)|^2, which has a clear physical meaning. This quantity is different from |fft(u^2)|, they are not proportional to each other -- you can check this directly on the numerical example.

[glwagner]

I guess we quantify something that we might call the "spectrum of kinetic energy"...

I disagree that the Eh we are plotting doesn't have a clear physical meaning --- the Eh we plot shows the spectrum of kinetic energy, or the amplitudes of the spectral components of kinetic energy.

But I think I see what @parfenyev is saying because the "energy spectrum" is defined so that its integral is equal to the total kinetic energy of the flow. E.g. according to Parseval's theorem,

∫ u² + v² dA = Σ |û|² + |v̂|²

(or something like that), so the integral of Ê = |û|² + |v̂|² in wavenumber space yields the total kinetic energy and maybe lends meaning to spectral slopes, etc. The wavenumber-integral of Eh as its currently defined is probably not equal to the total kinetic energy, right?

So, the motivation for revising the definition of Eh is so that the integral of Eh(r) over the "radial wavenumber" r yields the total kinetic energy.

[parfenyev]

Yes, you are right. Another motivation for revising the definition of Eh is the slope of the radial spectrum Ehr, which is widely discussed in scientific research and for which there are some theoretical predictions -- e.g., see Section 4 and Fig. 9 here.

[parfenyev]

In addition to the initial comment, I found strange behavior of the function FourierFlows.radialspectrum(). As an example, I took some function quant(x,y)=sin(k*x)+sin(k*y) and calculated its Fourier transform quant_h. Then I calculated the radial spectrum via

kr, quant_hr = FourierFlows.radialspectrum(quant_h, grid, refinement=2)

I expected to see a very narrow peak at kr=k, but it turned out to be quite wide. The width does not decrease as the value of the refinement increases. Below I am attaching the corresponding graph and source code.

@navidcy Is this the correct behavior of the function that you would expect?

Source code to reproduce example:

using GeophysicalFlows, Plots
using LinearAlgebra: mul!

dev = CPU()   
n, L  = 512, 2π  

grid = TwoDGrid(dev, n, L)
x, y = grid.x, grid.y

wavenumber = 10.0 * 2π/L 

# expression for some quantity in real (x,y) space:
quant = [sin(wavenumber*i)+sin(wavenumber*j) for i in grid.x, j in grid.y]

# let us move to the Fourier space:
quant_h = Complex.(zeros(Int(n/2)+1, Int(n)))
mul!(quant_h, grid.rfftplan, deepcopy(quant))

#let us calculate the radial spectrum:
quant_hr = Complex.(zeros(2*Int(n),1))
kr, quant_hr = FourierFlows.radialspectrum(quant_h, grid, refinement=2)

plt=plot(kr, abs.(quant_hr),
    linewidth = 2,
        alpha = 0.7,
       xlabel = "k", ylabel = "Radial Spectrum",
        xlims = (5e-1, grid.nx), ylims = (1e-13,1e6),
       xscale = :log10, yscale = :log10,
       legend = :topright, label="k=10.0",
   framestyle = :box)

#Plots.savefig(plt, "./p1.png")

[parfenyev]

In the souce code of function radialspectrum one defines:

lshift = range(-grid.nl/2+1, stop=grid.nl/2, length=grid.nl) * 2π/grid.Ly
kshift = range(-grid.nk/2+1, stop=grid.nk/2, length=grid.nk) * 2π/grid.Lx

However, it seems that the correct definition consistent with the operation fftshift is

lshift = range(-grid.nl/2, stop=grid.nl/2-1, length=grid.nl) * 2π/grid.Ly
kshift = range(-grid.nk/2, stop=grid.nk/2-1, length=grid.nk) * 2π/grid.Lx

Then the result is more like expectations:

[navidcy]

oh I see! you are right.

julia> using FFTW

julia> n = 8
8

julia> fftfreq(n, n)
8-element Frequencies{Float64}:
  0.0
  1.0
  2.0
  3.0
 -4.0
 -3.0
 -2.0
 -1.0

this is from -n/2:n/2-1.

Could you open a PR for this?

[navidcy]

We should construct kshift and lshift using FFTW functionality, e.g.,

julia> using FFTW

julia> fftshift(fftfreq(8, 8))
-4.0:1.0:3.0

[parfenyev]

I'll try to open PR, but I've never done it before. I need some time to understand how to do this and how to work with github. Maybe you know a manual that can help?

[glwagner]

This might help: https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/proposing-changes-to-your-work-with-pull-requests/about-pull-requests

[navidcy]

(Also we are happy to do that if you prefer.)