在Fortran中使用r2c和c2r FFTW时，前向和后向尺寸是否相同? [英] When using r2c and c2r FFTW in Fortran, are the forward and backward dimensions same?

查看：205
本文介绍了在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

FFTW_PRESERVE_INPUT指定非原位转换一定不能更改其输入数组.通常，这是默认设置，除了c2r和hc2r(即复杂到实数)转换(FFTW_DESTROY_INPUT是默认设置...

FFTW_PRESERVE_INPUT specifies that an out-of-place transform must not change its input array. This is ordinarily the default, except for c2r and hc2r (i.e. complex-to-real) transforms for which FFTW_DESTROY_INPUTis the default...

例如，我们可以通过@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时，前向和后向尺寸是否相同?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！

查看全文
相关文章
其他开发最新文章
热门教程
热门工具
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆