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

Question

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

使用时：dfftw_plan_dft_r2c_2d 正确的解决方案

使用 2 ffts 时：不正确的解决方案

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

我的主要问题是：与沿行和沿列dfftw_plan_dft_r2c_2d运行时有何不同？dfftw_plan_dft_r2c_1ddfftw_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,Ny),rho_dum(Nx,Ny),phi1(Nx)
complex ( kind = 8 ) rhok(Nh,Ny),rhok_dum(Nh,Ny),phik(Nh,Ny),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, Nx, Ny, 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) = rhok(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) = rhok(i,j)/( kx*kx + ky*ky )
  phik_dum(i,j) = phik(i,j) 
  enddo
enddo

do i = 1,Nh
    phiy_dumk = phik_dum(i,:)
    call dfftw_plan_dft_1d(plan_backwardz, Ny, 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, Nx, phix_dumk, phi1,FFTW_BACKWARD,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, Nx, Ny, 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, Nx
  do j = 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=',Nx,',j=',Ny

do j=1,Ny
do i=1,Nx
write(100,*) x(i),y(j),rho(i,j),phi(i,j)
end do
end do
close(100)

end program poisson2d

标签： fortranfftw