Aliasing/Polluted Result from rfftplan for Float32 2^n 3D array

[doraemonho]

This issue was originally from the discussion in FourierFlows but I think I should better move it here.

Environment: WSL2 Ubuntu 18.04/20.04
Julia version: Julia 1.7.3/1.8.2
CUDA version: 11.8

The brief summary are,

1. For a 3D Array with size $2^n$ (for example, $128^3$)
2. If we define a function F(k) in 3-dimensional k-space as a gaussian function like $e^{-(k-k_0)^2}e^{-i\phi}$ where $k_0$ is some number and $\phi$ is a random number array $\in [0, 2\pi]$.
3. Construct a rfftplan with type Float32
4. The result of complex to real transformation contains the pollution of the aliasing effect

The Figure below illustrates the issue

---

However, if we change the size from 2^n to other numbers such as 10*n (for example, 130) or use Float64, this aliasing would disappear

To better shows the problem, I construct a test code for that, and here is the graphical comparison

using CUDA
using FFTW
using Random
using LinearAlgebra: mul!, ldiv!
using Test


function getk⁻²_and_rfftplan(dev::Dev; 
                             nx = 64, ny = nx, nz = nx, Lx = 2π, Ly = Lx, Lz = Lx,
                             T = Float32, nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE) where Dev
    ArrayType = dev == "CPU" ? Array : CuArray;
    
    # Define the rfft parameters
    nk = nx
    nl = ny
    nm = nz
    nkr = Int(nx/2 + 1)  
    # Wavenubmer
    k = ArrayType{T}(reshape( fftfreq(nx, 2π/Lx*nx), (nk, 1, 1)))
    l = ArrayType{T}(reshape( fftfreq(ny, 2π/Ly*ny), (1, nl, 1)))
    m = ArrayType{T}(reshape( fftfreq(nz, 2π/Lz*nz), (1, 1, nm)))
    kr = ArrayType{T}(reshape(rfftfreq(nx, 2π/Lx*nx), (nkr, 1, 1)))

       Ksq = @. k^2 + l^2 + m^2
    invKsq = @. 1 / Ksq
    CUDA.@allowscalar  invKsq[1, 1, 1] = 0;

       Krsq = @. kr^2 + l^2 + m^2
    invKrsq = @. 1 / Krsq
    CUDA.@allowscalar invKrsq[1, 1, 1] = 0;

    FFTW.set_num_threads(nthreads);
    rfftplan = plan_rfft(ArrayType{T, 3}(undef, nx, ny, nz))
    
    return invKrsq,rfftplan
end


for nx in [130,128]
    println("---------------------------------------------------")
    println("Teting for nx = $nx")
    println("---------------------------------------------------")
    ny, nz = nx, nx, nx;       # define size of the array
    Random.seed!(1234);        # Set up the random seed
    θ = rand(div(nx,2)+1, ny, nz)*2π; # define \theta \in [0,2\pi]
    
    for T in [Float64,Float32]
          k⁻²,  rfftplan = getk⁻²_and_rfftplan("CPU"; nx = nx, Lx = 2π, ny = ny, nz = ny, T =  T)             
        k⁻²_c,rfftplan_c = getk⁻²_and_rfftplan("GPU"; nx = nx, Lx = 2π, ny = ny, nz = ny, T =  T)

        θT  = convert(Array{T,3},θ)
        eⁱᶿk⁻² = @. exp.(im*θT)*k⁻²;     #define the spectral filtered function in CPU
        eⁱᶿk⁻²_cu = CuArray(eⁱᶿk⁻²);     #define the spectral filtered function in GPU

        println("before the rfft")
        @show typeof(eⁱᶿk⁻²)
        @show typeof(eⁱᶿk⁻²_cu)

        c = zeros(T, nx, ny, nz)
        c_cu = CuArray(c);

        ldiv!(c, rfftplan, deepcopy(eⁱᶿk⁻²))    
        ldiv!(c_cu, rfftplan_c, deepcopy(eⁱᶿk⁻²_cu))
        c_cu = Array(c_cu); #move the result from GPU-> CPU


        println("before the test")
        @show typeof(c)
        @show typeof(c_cu)

        test_title = "eⁱᶿk⁻² CPU & GPU test for " * string(T)

        @testset "$test_title" begin
            #Compare the result using ≈
            @test c_cu  ≈  c
        end 
    end
end

[maleadt]

Do you have any idea what causes this? Why do you think CUDA.jl's wrappers are to blame, and not the underlying CUFFT libraries (in which case it would be better to file a bug with them)?

[doraemonho]

Thanks for checking it out. I am not meaning to blame CUDA.jl.
The reason I am posting in here mainly because I perform the same calculation before in Julia long time ago and it works normally.
I just rerun the same calculation in v1.5.3 Julia to confirm that and the aliasing just disappeared.

[maleadt]

What version of CUDA and CUDA.jl were you using with Julia 1.5.3? If the CUDA.jl version is different, a different toolkit could have been loaded. The version of AbstractFFTs.jl could also matter

To rule out a bug in the CUDA toolkit, try the following: With the latest released version of CUDA.jl, run your code under various settings of the JULIA_CUDA_VERSION environment variable (e.g. 10.2, 11.0, 11.1, etc). Verify with CUDA.versioninfo() that a different toolkit has been loaded.

[doraemonho]

What version of CUDA and CUDA.jl were you using with Julia 1.5.3? If the CUDA.jl version is different, a different toolkit could have been loaded. The version of AbstractFFTs.jl could also matter

To rule out a bug in the CUDA toolkit, try the following: With the latest released version of CUDA.jl, run your code under various settings of the JULIA_CUDA_VERSION environment variable (e.g. 10.2, 11.0, 11.1, etc). Verify with CUDA.versioninfo() that a different toolkit has been loaded.

Thank you for your suggestion. After some trial and error, I located the bug would happen when the CuFFT library version > 10.3.0.
Let me close the issue and report it to NVIDIA.