【发布时间】：2013-05-25 12:02:35
【问题描述】：
我在用 FFTW 处理这个简单的正弦导数。乍一看似乎没问题，但与确切的解决方案相比，有一个很大的错误（5e-6）...... 我确实看到在使用 c2r 之后，复杂的输入都搞砸了，但在我看来，同样的复杂输入是我的问题的原因......我做错了什么？我没有使用任何指针，并试图让一切尽可能简单；我仍然无法弄清楚出了什么问题。 任何帮助表示赞赏！谢谢！！！
program main
! C binding
use, intrinsic :: iso_c_binding
implicit none
double precision, parameter :: pi = 4*ATAN(1.0)
complex, parameter :: ii =(0.0,1.0)
integer(C_INT), parameter :: Nx = 32
integer(C_INT), parameter :: Ny = Nx
double precision, parameter :: Lx = 2*pi, Ly = 2*pi
! Derived paramenter
double precision, parameter :: dx = Lx/Nx, dy = Ly/Ny
real(C_DOUBLE), dimension(Nx,Ny) :: x,y, u0,in,dudx,dudxE, errdU
real(C_DOUBLE), dimension(Nx/2+1,Ny) :: kx, ky
! Fourier space variables
complex(C_DOUBLE_COMPLEX), dimension(Nx/2+1,Ny) :: u_hat_x, out
! indices
integer :: i, j
!---FFTW plans
type(C_PTR) :: pf, pb
! FFTW include
include 'fftw3.f03'
write(*,'(A)') 'Starting...'
! Grid
forall(i=1:Nx,j=1:Ny)
x(i,j) = (i-1)*dx
y(i,j) = (j-1)*dy
end forall
! Compute the wavenumbers
forall(i=1:Nx/2,j=1:Ny) kx(i,j)=2*pi*(i-1)/Lx
kx(Nx/2+1,:) = 0.0
forall(i=1:Nx/2+1,j=1:Ny/2) ky(i,j)=2*pi*(j-1)/Ly
forall(i=1:Nx/2+1,j=Ny/2+1:Ny) ky(i,j)=2*pi*(j-Ny-1)/Ly
! Initial Condition
u0 = sin(2*x)
dudxE = 2*cos(2*x)
! Go to Fourier Space
in = u0
pf = fftw_plan_dft_r2c_2d(Ny, Nx, in,out ,FFTW_ESTIMATE)
call fftw_execute_dft_r2c(pf,in,out)
u_hat_x = out
! Derivative
out = ii*kx*out
! Back to physical space
pb = fftw_plan_dft_c2r_2d(Ny, Nx, out,in,FFTW_ESTIMATE)
call fftw_execute_dft_c2r(pb,out,in)
! rescale
dudx = in/Nx/Ny
! check the derivative
errdU = dudx - dudxE
! Write file
write(*,*) 'Writing to files...'
OPEN(UNIT=1, FILE="out_for.dat", ACTION="write", STATUS="replace", &
FORM="unformatted")
WRITE(1) kx,u0,dudx,errdU,abs(u_hat_x)
CLOSE(UNIT=1)
call fftw_destroy_plan(pf)
call fftw_destroy_plan(pb)
write(*,'(A)') 'Done!'
end program main
【问题讨论】：
-
这个
pi = 4*ATAN(1.0)是有问题的，因为你希望 pi 是双精度的，但它会从单精度转换，你需要使用类似pi = 4*ATAN(1.0_8)的东西，或者更好的是类似pi = 4*ATAN(1.0_C_DOUBLE)或pi = 4*ATAN(1.0d0)。这是因为ATAN函数返回的精度与参数相同，而没有任何_kind后缀的常量1.0只是一个单精度实变量。 -
也许我遗漏了一些东西，但为什么所有的二维数组？您正在对函数进行一维导数 (d/dx)，为什么不只使用一维数组？
-
感谢您的回答，我会尝试更改 pi 的定义，但我不认为这会导致我的问题，因为如果我 FFT 和 IFFT 结果很好。我得到的错误也有一个模式，它看起来非常像 Gibbs 错误，x 方向边界旁边的高频振荡（我正在推导的那个） 1D 工作正常，但我需要 2D ，我从一些更简单的东西开始，在深入研究复杂的功能之前把它弄明白。
标签： fortran fftw derivative