search

FourierFlows / FourierFlows.jl

Tools for building fast, hackable, pseudospectral partial differential equation solvers on periodic domains
https://bit.ly/FourierFlows
MIT License
202 stars 29 forks source link

Floating Point Error with radialspectrum #367

Open ndefilippis opened 5 months ago

ndefilippis commented 5 months ago

Hi,

It looks like in some cases there can be a floating point error in utils.jl/radialspectrum. Here is a MWE

using FourierFlows
Lx = 3.75e5 * 2π
nx = 256
grid = TwoDGrid(CPU(); Lx, nx)

field = zeros(Complex{Float64}, (grid.nkr, grid.nl))

FourierFlows.radialspectrum(field, grid)

which gives the error

BoundsError: attempt to access 129×256 scale(interpolate(::Matrix{ComplexF64}, BSpline(Linear())), (0.0:2.666666666666667e-6:0.00034133333333333335, -0.00034133333333333335:2.666666666666667e-6:0.00033866666666666664)) with element type ComplexF64 at index [2.073735246556185e-20, 0.0003386666666666667]

Stacktrace:
 [1] throw_boundserror(A::ScaledInterpolation{ComplexF64, 2, Interpolations.BSplineInterpolation{ComplexF64, 2, Matrix{ComplexF64}, BSpline{Linear{Throw{OnGrid}}}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, BSpline{Linear{Throw{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, I::Tuple{Float64, Float64})
   @ Base .\abstractarray.jl:651
 [2] ScaledInterpolation
   @ \Interpolations\nDwIa\src\scaling\scaling.jl:72 [inlined]
 [3] radialspectrum(fh::Matrix{ComplexF64}, grid::TwoDGrid{Float64, Matrix{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, FFTW.cFFTWPlan{ComplexF64, -1, false, 2, UnitRange{Int64}}, FFTW.rFFTWPlan{Float64, -1, false, 2, UnitRange{Int64}}, UnitRange{Int64}, CPU}; n::Nothing, m::Nothing, refinement::Int64)
   @ FourierFlows \FourierFlows\0jzfd\src\utils.jl:269
 [4] radialspectrum(fh::Matrix{ComplexF64}, grid::TwoDGrid{Float64, Matrix{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, FFTW.cFFTWPlan{ComplexF64, -1, false, 2, UnitRange{Int64}}, FFTW.rFFTWPlan{Float64, -1, false, 2, UnitRange{Int64}}, UnitRange{Int64}, CPU})
   @ FourierFlows \FourierFlows\0jzfd\src\utils.jl:241
 [5] top-level scope
   @ In[74]:8

It looks like the error stems from the creation of the interpolation scalings, e.g.

utils.jl:245    lshift = range(-nl/2, stop=nl/2-1, length=nl) * 2π/Ly

not always matching ρmax.

I think a potential fix would be to do the scaling by 2π/Ly inside the range, as in:

utils.jl:245    lshift = range(-nl/2 * 2π/Ly, stop=(nl/2-1) * 2π/Ly, length=nl)

and similarly for kshift

navidcy commented 5 months ago

Thanks for pointing this out!

I'm puzzling to understand why this appears with Lx = 3.75e5 and not at Lx=4e5... In both cases ρmax ≈ maximum(lshift) = true...