Closed
Description
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,
- For a 3D Array with size
$2^n$ (for example,$128^3$ ) - 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]$ . - Construct a rfftplan with type Float32
- 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