Fortran-FFTW 2D FFT 结果不同于 2D 矩阵的 2 1D FFT

问题描述

使用 FFTW 实现 2D FFT/IFFT。我正在使用 Fortran 和 FFTW 开发二维快速泊松求解器。对于我的求解器，我需要它能够沿行运行 1D FFT，然后沿列运行 1D FFT。为了进行基准测试，我使用了 2D FFT 函数，并且能够使用 2D FFTW 函数求解：dfftw_plan_dft_r2c_2d 但是当我尝试沿行使用 dfftw_plan_dft_r2c_1d 和 dfftw_plan_dft_1d（复杂到-complex) 沿着列我没有收到相同的结果。

使用时：dfftw_plan_dft_r2c_2d Correct solution

使用 2 个 fft 时：Incorrect solution

我相信正在发生的是，当我沿着行运行第一个 FFT（实数到复数）时，它的大小从 Nx:128 变为 Nh:(128/2)+1 并运行第二个 FFT (复杂到复杂）沿着 i:1,Nh 的列，我会丢失沿 x 的这些数据，因为 x 现在定义为 i:1,Nh。

我的主要问题是：dfftw_plan_dft_r2c_2d 与沿行运行 dfftw_plan_dft_r2c_1d 和沿列运行 dfftw_plan_dft_1d 时有何不同？

Program poisson2d
implicit none
integer ( kind = 4 ),parameter :: Nx = 128
integer ( kind = 4 ),parameter :: Ny = 128
integer ( kind = 4 ),parameter :: Nh = ( Nx / 2 ) + 1
real ( kind = 8 ),parameter :: pi=3.14159265358979323846d0
include "fftw3.f"
integer ( kind = 4 ) i,j
real ( kind = 8 ) Lx,Ly,dx,dy,kx,ky
real ( kind = 8 ) x(Nx),y(Ny),phi(Nx,Ny),rho(Nx,rho_dum(Nx,phi1(Nx)
complex ( kind = 8 ) rhok(Nh,rhok_dum(Nh,phik(Nh,phik_dum(Nh,Ny)
complex ( kind = 8 ) phiy(Ny),phiy_dumk(Ny),phix_dumk(Nh),rhox(Nh),rhoy(Ny)
complex ( kind = 8 ) rho_dumx(Nx),rho_dumy(Ny)
integer ( kind = 8 ) plan_forward,plan_backward,plan_backwardz
integer ( kind = 8 ) plan_forwardy1,plan_backwardy1

open(unit=10,file='poisson.dat',status='unkNown')
open(unit=11,file='poisson2.dat',status='unkNown')

Lx = 2.0*pi
Ly = 2.0*pi
dx = Lx/float(Nx)
dy = Ly/float(Ny)
do i = 1,Nx
x(i)=0.0d0+real(i-1)*dx
  do j = 1,Ny
  y(j)=0.0+real(j-1)*dy
  rho(i,j) = 2.0d0*sin(x(i))*cos(y(j))
  rho_dum(i,j) = rho(i,j)
  end do
end do

!!FFT along the rows

do j = 1,Ny
    rho_dumx = rho_dum(:,j)
    call dfftw_plan_dft_r2c_1d(plan_forward,Nx,rho_dumx,rhox,FFTW_ESTIMATE)
    call dfftw_execute_ (plan_forward)
    rhok(:,j) = rhox
enddo

!!FFT along the columns

do i = 1,Nh
    rho_dumy = rhok(i,:)
    call dfftw_plan_dft_1d(plan_forwardy1,Ny,rho_dumy,rhoy,FFTW_FORWARD,FFTW_ESTIMATE)
    call dfftw_execute_ (plan_forwardy1)
    rhok(i,:) = rhoy
enddo

do i = 1,Nx
  do j = 1,Ny
  rhok(i,j) = rhok(i,j)
  write(11,*) rhok(i,j)
  end do
end do

!!2D FFT commented out

  !call dfftw_plan_dft_r2c_2d(plan_forward,rho_dum,rhok,FFTW_ESTIMATE)
  !call dfftw_execute_ (plan_forward)

do i = 1,Nx/2+1
  do j = 1,Ny/2
  kx = 2.0d0*pi*float(i-1)/Lx
  ky = 2.0d0*pi*float(j-1)/Ly
    if (i == 1 .and. j == 1) then
    phik(i,j) = (0.0d0,0.0d0)
    else
    phik(i,j)/( kx*kx + ky*ky )
    endif
  phik_dum(i,j) = phik(i,j) 
  enddo
  do j = Ny/2+1,Ny
  kx = 2.0d0*pi*float(i-1)/Lx
  ky = 2.0d0*pi*float((j-1)-Ny)/Ly
  phik(i,j)/( kx*kx + ky*ky )
  phik_dum(i,j) 
  enddo
enddo

do i = 1,Nh
    phiy_dumk = phik_dum(i,:)
    call dfftw_plan_dft_1d(plan_backwardz,phiy_dumk,phiy,FFTW_BACKWARD,FFTW_ESTIMATE)
    call dfftw_execute_ (plan_backwardz)
    phik_dum(i,:) = phiy
enddo

do j = 1,Ny
    phix_dumk = phik_dum(:,j)
    call dfftw_plan_dft_c2r_1d(plan_backward,phix_dumk,phi1,FFTW_ESTIMATE)
    call dfftw_execute_ (plan_backward)
    phi(:,j) = phi1
enddo

  call dfftw_destroy_plan(plan_backwardz)
  call dfftw_destroy_plan(plan_forward)
  call dfftw_destroy_plan(plan_backward)
  call dfftw_destroy_plan(plan_forwardy1)

!!2D FFT commented out

  !call dfftw_plan_dft_c2r_2d_ (plan_backward,phik_dum,phi,FFTW_ESTIMATE)
  !call dfftw_execute_ (plan_backward)
  !call dfftw_destroy_plan_ (plan_forward)
  !call dfftw_destroy_plan_ (plan_backward)
do i = 1,Ny
  phi(i,j) = phi(i,j)/(float(Nx)*float(Ny))
  write(10,*) x(i),y(j),rho(i,j),phi(i,j)
  end do
end do

open(100,file='field4_51.plt')
write(100,*)'variables ="x","y","rho","phi"'
write(100,*)'zone f=point i=',',j=',Ny

do j=1,Ny
do i=1,Nx
write(100,j)
end do
end do
close(100)

end program poisson2d

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

fftw fortran