本文介绍了在Fortran中使用r2c和c2r FFTW时，前向和后向尺寸是否相同?的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
打击是主文件
PROGRAM SPHEROID
USE nrtype
USE SUB_INFO
INCLUDE "/usr/local/include/fftw3.f"
INTEGER(I8B) :: plan_forward, plan_backward
INTEGER(I4B) :: i, t, int_N
REAL(DP) :: cth_i, sth_i, real_i, perturbation
REAL(DP) :: PolarEffect, dummy, x1, x2, x3
REAL(DP), DIMENSION(4096) :: dummy1, dummy2, gam, th, ph
REAL(DP), DIMENSION(4096) :: k1, k2, k3, k4, l1, l2, l3, l4, f_in
COMPLEX(DPC), DIMENSION(2049) :: output1, output2, f_out
CHARACTER(1024) :: baseOutputFilename
CHARACTER(1024) :: outputFile, format_string
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
int_N = 4096
! File Open Section
format_string = '(I5.5)'
! Write the coodinates at t = 0
do i = 1, N
real_i = real(i)
gam(i) = 2d0*pi/real_N
perturbation = 0.01d0*dsin(2d0*pi*real_i/real_N)
ph(i) = 2d0*pi*real_i/real_N + perturbation
th(i) = pi/3d0 + perturbation
end do
! Initialization Section for FFTW PLANS
call dfftw_plan_dft_r2c_1d(plan_forward, int_N, f_in, f_out, FFTW_ESTIMATE)
call dfftw_plan_dft_c2r_1d(plan_backward, int_N, f_out, f_in, FFTW_ESTIMATE)
! Runge-Kutta 4th Order Method Section
do t = 1, Iter_N
call integration(th, ph, gam, k1, l1)
do i = 1, N
dummy1(i) = th(i) + 0.5d0*dt*k1(i)
end do
do i = 1, N
dummy2(i) = ph(i) + 0.5d0*dt*l1(i)
end do
call integration(dummy1, dummy2, gam, k2, l2)
do i = 1, N
dummy1(i) = th(i) + 0.5d0*dt*k2(i)
end do
do i = 1, N
dummy2(i) = ph(i) + 0.5d0*dt*l2(i)
end do
call integration(dummy1, dummy2, gam, k3, l3)
do i = 1, N
dummy1(i) = th(i) + dt*k3(i)
end do
do i = 1, N
dummy2(i) = ph(i) + dt*l3(i)
end do
call integration(dummy1, dummy2, gam, k4, l4)
do i = 1, N
cth_i = dcos(th(i))
sth_i = dsin(th(i))
PolarEffect = (nv-sv)*dsqrt(1d0+a*sth_i**2) + (nv+sv)*cth_i
PolarEffect = PolarEffect/(sth_i**2)
th(i) = th(i) + dt*(k1(i) + 2d0*k2(i) + 2d0*k3(i) + k4(i))/6d0
ph(i) = ph(i) + dt*(l1(i) + 2d0*l2(i) + 2d0*l3(i) + l4(i))/6d0
ph(i) = ph(i) + dt*0.25d0*PolarEffect/pi
end do
!! Fourier Filtering Section
call dfftw_execute_dft_r2c(plan_forward, th, output1)
do i = 1, N/2+1
dummy = abs(output1(i))
if (dummy.lt.threshhold) then
output1(i) = dcmplx(0.0d0)
end if
end do
call dfftw_execute_dft_c2r(plan_backward, output1, th)
do i = 1, N
th(i) = th(i)/real_N
end do
call dfftw_execute_dft_r2c(plan_forward, ph, output2)
do i = 1, N/2+1
dummy = abs(output2(i))
if (dummy.lt.threshhold) then
output2(i) = dcmplx(0.0d0)
end if
end do
call dfftw_execute_dft_c2r(plan_backward, output2, ph)
do i = 1, N
ph(i) = ph(i)/real_N
end do
!! Data Writing Section
write(baseOutputFilename, format_string) t
outputFile = "xyz" // baseOutputFilename
open(unit=7, file=outputFile)
outputFile = "Fsptrm" // baseOutputFilename
open(unit=8, file=outputFile)
do i = 1, N
x1 = dsin(th(i))*dcos(ph(i))
x2 = dsin(th(i))*dsin(ph(i))
x3 = dsqrt(1d0+a)*dcos(th(i))
write(7,*) x1, x2, x3
end do
do i = 1, N/2+1
write(8,*) abs(output1(i)), abs(output2(i))
end do
close(7)
close(8)
do i = 1, N/2+1
output1(i) = dcmplx(0.0d0)
end do
do i = 1, N/2+1
output2(i) = dcmplx(0.0d0)
end do
end do
! Destroying Process for FFTW PLANS
call dfftw_destroy_plan(plan_forward)
call dfftw_destroy_plan(plan_backward)
END PROGRAM
下面是用于集成的子例程文件
Below is a subroutine file for integration
! We implemented Shelly's spectrally accurate convergence method
SUBROUTINE integration(in1,in2,in3,out1,out2)
USE nrtype
USE SUB_INFO
INTEGER(I4B) :: i, j
REAL(DP) :: th_i, th_j, gi, ph_i, ph_j, gam_j, v1, v2
REAL(DP), DIMENSION(N), INTENT(INOUT) :: in1, in2, in3, out1, out2
REAL(DP) :: ui, uj, part1, part2, gj, cph, sph
REAL(DP) :: denom, numer, temp
do i = 1, N
out1(i) = 0d0
end do
do i = 1, N
out2(i) = 0d0
end do
do i = 1, N
th_i = in1(i)
ph_i = in2(i)
ui = dcos(th_i)
part1 = dsqrt(1d0+a)/(dsqrt(-a)*ui+dsqrt(1d0+a-a*ui*ui))
part1 = part1**(dsqrt(-a))
part2 = (dsqrt(1d0+a-a*ui*ui)+ui)/(dsqrt(1d0+a-a*ui*ui)-ui)
part2 = dsqrt(part2)
gi = dsqrt(1d0-ui*ui)*part1*part2
do j = 1, N
if (mod(i+j,2).eq.1) then
th_j = in1(j)
ph_j = in2(j)
gam_j = in3(j)
uj = dcos(th_j)
part1 = dsqrt(1d0+a)/(dsqrt(-a)*uj+dsqrt(1d0+a-a*uj*uj))
part1 = part1**(dsqrt(-a))
part2 = (dsqrt(1d0+a-a*uj*uj)+uj)/(dsqrt(1d0+a-a*uj*uj)-uj)
part2 = dsqrt(part2)
gj = dsqrt(1d0-ui*ui)*part1*part2
cph = dcos(ph_i-ph_j)
sph = dsin(ph_i-ph_j)
numer = dsqrt(1d0-uj*uj)*sph
denom = (gj/gi*(1d0-ui*ui) + gi/gj*(1d0-uj*uj))*0.5d0
denom = denom - dsqrt((1d0-ui*ui)*(1d0-uj*uj))*cph
denom = denom + krasny_delta
v1 = -0.25d0*gam_j*numer/denom/pi
temp = dsqrt(1d0+(1d0-ui*ui)*a)
numer = -(gj/gi)*(temp+ui)
numer = numer + (gi/gj)*((1d0-uj*uj)/(1d0-ui*ui))*(temp-ui)
numer = numer + 2d0*ui*dsqrt((1d0-uj*uj)/(1d0-ui*ui))*cph
numer = 0.5d0*numer
v2 = -0.25d0*gam_j*numer/denom/pi
out1(i) = out1(i) + 2d0*v1
out2(i) = out2(i) + 2d0*v2
end if
end do
end do
END
下面是一个模块文件
module nrtype
Implicit none
!integer
INTEGER, PARAMETER :: I8B = SELECTED_INT_KIND(20)
INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4)
INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2)
!real
INTEGER, PARAMETER :: SP = KIND(1.0)
INTEGER, PARAMETER :: DP = KIND(1.0D0)
!complex
INTEGER, PARAMETER :: SPC = KIND((1.0,1.0))
INTEGER, PARAMETER :: DPC = KIND((1.0D0,1.0D0))
!defualt logical
INTEGER, PARAMETER :: LGT = KIND(.true.)
!mathematical constants
REAL(DP), PARAMETER :: pi = 3.141592653589793238462643383279502884197_dp
!derived data type s for sparse matrices,single and double precision
!User-Defined Constants
INTEGER(I4B), PARAMETER :: N = 4096, Iter_N = 20000
REAL(DP), PARAMETER :: real_N = 4096d0
REAL(DP), PARAMETER :: a = -0.1d0, dt = 0.001d0, krasny_delta = 0.01d0
REAL(DP), PARAMETER :: nv = 0d0, sv = 0d0, threshhold = 0.00000000001d0
!N : The Number of Point Vortices, Iter_N * dt = Total time, dt : Time Step
!krasny_delta : Smoothing Parameter introduced by R.Krasny
!nv : Northern Vortex Strength, sv : Southern Vortex Strength
!a : The Eccentricity in the direction of z , threshhold : Filtering Threshhold
end module nrtype
下面是一个子例程信息文件
Below is a subroutine info file
MODULE SUB_INFO
INTERFACE
SUBROUTINE integration(in1,in2,in3,out1,out2)
USE nrtype
INTEGER(I4B) :: i, j
REAL(DP) :: th_i, th_j, gi, ph_i, ph_j, gam_j, v1, v2
REAL(DP), DIMENSION(N), INTENT(INOUT) :: in1, in2, in3, out1, out2
REAL(DP) :: ui, uj, part1, part2, gj, cph, sph
REAL(DP) :: denom, numer, temp
END SUBROUTINE
END INTERFACE
END MODULE
我使用以下命令对其进行了编译
I compiled them using the below command
gfortran -o p0 -fbounds-check nrtype.f90 spheroid_sub_info.f90 spheroid_sub_integration.f90 spheroid_main.f90 -lfftw3 -lm -Wall -pedantic -pg
nohup ./p0&
nohup ./p0 &
请注意2049 = 4096/2 +1
Note that 2049 = 4096 / 2 + 1
在进行plan_backward时，由于输出尺寸为2049，我们使用2049而不是4096是否正确?
When making plan_backward, isn't it correct that we use 2049 instead of 4096 since the dimension of output is 2049?
但是当我这样做时，它会炸毁. (爆炸意味着NAN错误)
But when I do that, it blows up. (Blowing up means NAN error)
如果我在制作plan_backward时使用4096，那么一切都很好，除了某些傅立叶系数异常大而不会发生.
If I use 4096 in making plan_backward, Everything is fine except that some Fourier coefficients are abnormally big which should not happen.
请帮助我在Fortran中正确使用FFTW.这个问题使我很久没气了.
Please help me use FFTW in Fortran correctly. This issue has discouraged me for a long time.
推荐答案
一个问题可能是dfftw_execute_dft_c2r可以破坏输入数组的内容，如本.关键摘录是
One issue may be that dfftw_execute_dft_c2r can destroy the content of the input array, as described in this page. The key excerpt is
例如，我们可以通过@VladimirF修改示例代码，以使其在第一次FFT(r2c)调用后立即将data_out保存为data_save，然后在第二次FFT( c2r)致电.因此，对于OP的代码，在进入第二个FFT(c2r)之前将output1和output2保存到不同的数组似乎更安全.
We can verify this, for example, by modifying the sample code by @VladimirF such that it saves data_out to data_save right after the first FFT(r2c) call, and then calculating their difference after the second FFT (c2r) call. So, in the case of OP's code, it seems safer to save output1 and output2 to different arrays before entering the second FFT (c2r).
这篇关于在Fortran中使用r2c和c2r FFTW时，前向和后向尺寸是否相同?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！