FourierFlows / GeophysicalFlows.jl

Geophysical fluid dynamics pseudospectral solvers with Julia and FourierFlows.jl.
https://fourierflows.github.io/GeophysicalFlowsDocumentation/stable/
MIT License
155 stars 32 forks source link

construct a random field with prescribed spectral slope #295

Open liasiegelman opened 2 years ago

liasiegelman commented 2 years ago

I would like to construct a random initial field with a prescribed spectral slope x. If Energy = \int E(k) k dk then E(k) ~ k^{-x}. As suggested by Navid, maybe we can iterate on peakedisotropicspectrum to implement this.

liasiegelman commented 2 years ago

I think this should do

function slopedisotropicspectrum(grid::TwoDGrid{T, A}, slope::Real; mask=ones(size(grid.Krsq)), allones=false) where {T, A}
 if grid.Lx !== grid.Ly
   error("the domain is not square")
 else
   modk = sqrt.(grid.Krsq)
   modψ = A(zeros(T, (grid.nk, grid.nl)))
   modψ = @. modk^((slope-1)/2)
   #modψ[1, 1] = 0.0
   CUDA.@allowscalar modψ[1, 1] = 0.0

   phases = randn(Complex{T}, size(grid.Krsq))
   phases_real, phases_imag = real.(phases), imag.(phases)
   phases = A(phases_real) + im * A(phases_imag)
   ψh = @. phases * modψ
   if allones; ψh = modψ; end
   ψh = ψh .* A(mask)
   q = A(irfft(ψh, grid.nx))
 end

 return q
end
navidcy commented 2 years ago

So what you want is a function that returns a field q for which if energy is $E = \int q(x, y)^2 \mathrm{d}x \mathrm{d}y = \int \mathcal{E}(k) k \mathrm{d}k$ you want $\mathcal{E}(k) \propto k^{-\chi}$. Am I correct?