FFTW, fortran의 간단한 파생어

나는 FFTW를 사용한이 사인의 간단한 파생어로 고심하고 있습니다. 처음에는 괜찮아 보이지만 정확한 해결책과 비교할 때 아주 큰 오류 (5e-6)가 있습니다 ... 나는 c2r을 취한 후에 복잡한 입력이 엉망이되는 것을 보았습니다. 같은 복잡한 입력이 내 문제의 원인입니다 ... 내가 뭘 잘못하고 있니? 나는 포인터를 사용하지 않고 모든 것을 가능한 한 단순하게 유지하려고 노력했다. 아직도 나는 틀린 것이 무엇인지 알 수 없다. 도움을 주시면 감사하겠습니다. 감사!!!FFTW, fortran의 간단한 파생어

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

출처

2013-05-25 Lorenzo

이 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'이 단 정밀도 실수 변수이기 때문입니다. – steabert

아마도 내가 놓친 것이지만 모든 2D 배열은 왜 그런가요? 함수의 1D 파생 함수 (d/dx)를 수행하고 있는데 왜 1D 배열을 사용하지 않습니까? –

답변 주셔서 감사합니다, 나는 파이의 정의를 변경하려고하지만, 내가 FFT 및 IFFT 결과가 좋기 때문에, 내 문제가 발생할 생각하지 않습니다. 또한 내가 얻은 오류는 패턴을 가지고 있는데, 깁스 오류와 매우 비슷하게 보입니다. x 방향의 경계선 옆에 고주파 진동이 있습니다. (파생 된 것과 관련하여) 1D는 잘 동작하지만 2D가 필요하고, 복잡한 기능으로 다이빙하기 전에 좀 더 간단하게 이해하기 시작했습니다. – Lorenzo

steabert가 맞았습니다! 문제는 단순히 단 정밀도 인 pi 정의를 사용하는 것입니다. 정말 고맙습니다!

출처

2013-05-26 00:38:09 Lorenzo

FFTW, fortran의 간단한 파생어

답변

관련 문제