c - 在 Fortran 中使用 FFTW 的 r2r 类型的 3-D 正弦变换中的泊松方程求解器

Question

我在 Fortran 中编写了以下代码，使用 r2r(实数到实数)类型的 FFTW sin 变换求解泊松方程。在这段代码中，首先我完成了正常的 FFTW，即 r2c(实数到复数)类型的数学函数 27*sin(3x)*sin(3y)*sin(3z) 然后将它除以27 (3*3+3*3+3*3) 计算输入函数的二阶导数。输入函数振幅的 3-D 图显示了它的正确振幅。 amplitude along z-axis against x-,y- co-ordinate .
然后 c2r(复数到实数)类型的逆 FFTW 重新生成输入函数，但振幅现在减小到 1，如 3-D 图所示 2 .这表示我的泊松方程求解器代码在正常 FFTW 的情况下工作正常。

Program Derivative

! To run this code: gcc dst_3d.c -c -std=c99 && gfortran derivative.f95 dst_3d.o -lfftw3 && ./a.out

implicit none

include "fftw3.f"

integer ( kind = 4 ), parameter :: Nx = 64
integer ( kind = 4 ), parameter :: Ny = 64
integer ( kind = 4 ), parameter :: Nz = 64

real ( kind = 8 ), parameter :: pi=3.14159265358979323846d0

integer ( kind = 4 ) i,j,k
real ( kind = 8 ) Lx,Ly,Lz,dx,dy,dz,kx,ky,kz
real ( kind = 8 ) x(Nx),y(Ny),z(Nz)

real ( kind = 8 ) in_dst(Nx,Ny,Nz),out_dst(Nx,Ny,Nz) ! DST
real ( kind = 8 ) in_k_dst(Nx,Ny,Nz),out_k_dst(Nx,Ny,Nz) ! DST

real ( kind = 8 ) in_dft(Nx,Ny,Nz),out_dft(Nx,Ny,Nz) ! DFT
complex ( kind = 8 ) in_k_dft(Nx/2+1,Ny,Nz),out_k_dft(Nx/2+1,Ny,Nz) ! DFT
integer ( kind = 8 ) plan_forward,plan_backward ! DFT

! System Size.
Lx = 2.0d0*pi; Ly = 2.0d0*pi; Lz = 2.0d0*pi 

! Grid Resolution.
dx = Lx/dfloat(Nx); dy = Ly/dfloat(Ny); dz = Lz/dfloat(Nz)

! =================================== INPUT ===========================================

! Initial Profile Details.
do i = 1, Nx
  x(i) = dfloat(i-1)*dx
  do j = 1, Ny
    y(j) = dfloat(j-1)*dy
    do k = 1, Nz
      z(k) = dfloat(k-1)*dz
      in_dst(i,j,k) = 27.0d0*sin(3.0d0*x(i))*sin(3.0d0*y(j))*sin(3.0d0*z(k))
      in_dft(i,j,k) = in_dst(i,j,k)
      write(10,*) x(i), y(j), z(k), in_dft(i,j,k)
    enddo
  enddo
enddo

! =================================== 3D DFT ===========================================

  call dfftw_plan_dft_r2c_3d_ (plan_forward, Nx, Ny, Nz, in_dft, in_k_dft, FFTW_ESTIMATE)
  call dfftw_execute_ (plan_forward)
  call dfftw_destroy_plan_ (plan_forward)

do i = 1, Nx/2+1
  do j = 1, Ny/2
    do k = 1, Nz/2
      kx = 2.0d0*pi*dfloat(i-1)/Lx
      ky = 2.0d0*pi*dfloat(j-1)/Ly
      kz = 2.0d0*pi*dfloat(k-1)/Lz
      out_k_dft(i,j,k) = in_k_dft(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo 
    do k = Nz/2+1, Nz
      kx = 2.0d0*pi*dfloat(i-1)/Lx
      ky = 2.0d0*pi*dfloat(j-1)/Ly
      kz = 2.0d0*pi*dfloat((k-1)-Nz)/Lz
      out_k_dft(i,j,k) = in_k_dft(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo
  enddo
  do j = Ny/2+1, Ny
    do k = 1, Nz/2
      kx = 2.0d0*pi*dfloat(i-1)/Lx
      ky = 2.0d0*pi*dfloat((j-1)-Ny)/Ly
      kz = 2.0d0*pi*dfloat(k-1)/Lz
      out_k_dft(i,j,k) = in_k_dft(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo
    do k = Nz/2+1, Nz
      kx = 2.0d0*pi*dfloat(i-1)/Lx
      ky = 2.0d0*pi*dfloat((j-1)-Ny)/Ly
      kz = 2.0d0*pi*dfloat((k-1)-Nz)/Lz
      out_k_dft(i,j,k) = in_k_dft(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo
  enddo   
enddo

out_k_dft(1,1,1) = in_k_dft(1,1,1)

  call dfftw_plan_dft_c2r_3d_ (plan_backward, Nx, Ny, Nz, out_k_dft, out_dft, FFTW_ESTIMATE)
  call dfftw_execute_ (plan_backward)
  call dfftw_destroy_plan_ (plan_backward)

do i = 1, Nx
  do j = 1, Ny
    do k = 1, Nz
    out_dft(i,j,k) = out_dft(i,j,k)/dfloat(Nx*Ny*Nz)
    write(20,*) x(i), y(j), z(k), out_dft(i,j,k)    
    enddo
  enddo   
enddo

! =================================== 3D DST ===========================================

  call Forward_FFT(Nx, Ny, Nz, in_dst, in_k_dst)

do k = 1, Nz
  do j = 1, Ny/2
    do i = 1, Nx/2
      kz = 2.0d0*pi*dfloat((k-1))/Lz
      ky = 2.0d0*pi*dfloat((j-1))/Ly
      kx = 2.0d0*pi*dfloat((i-1))/Lx
      out_k_dst(i,j,k) = in_k_dst(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo 
    do i = Nx/2+1, Nx
      kz = 2.0d0*pi*dfloat((k-1))/Lz
      ky = 2.0d0*pi*dfloat((j-1))/Ly
      kx = 2.0d0*pi*dfloat(Nx-(i-1))/Lx
      out_k_dst(i,j,k) = in_k_dst(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo
  enddo  
  do j = Ny/2+1, Ny
    do i = 1, Nx/2
      kz = 2.0d0*pi*dfloat((k-1))/Lz
      ky = 2.0d0*pi*dfloat(Ny-(j-1))/Ly
      kx = 2.0d0*pi*dfloat((i-1))/Lx
      out_k_dst(i,j,k) = in_k_dst(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo
    do i = Nx/2+1, Nx
      kz = 2.0d0*pi*dfloat((k-1))/Lz
      ky = 2.0d0*pi*dfloat(Ny-(j-1))/Ly
      kx = 2.0d0*pi*dfloat(Nx-(i-1))/Lx
      out_k_dst(i,j,k) = in_k_dst(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo
  enddo   
enddo

out_k_dst(1,1,1) = in_k_dst(1,1,1)

  call Backward_FFT(Nx, Ny, Nz, out_k_dst, out_dst)  

do i = 1, Nx
  do j = 1, Ny
    do k = 1, Nz
    out_dst(i,j,k) = out_dst(i,j,k)/dfloat(8*Nx*Ny*Nz)
    write(30,*) x(i), y(j), z(k), out_dst(i,j,k)
    enddo
  enddo   
enddo

end program Derivative

! ================================== FFTW SUBROUTINES 
====================================================

! ================================================================= !
! Wrapper Subroutine to call forward fftw c functions for 3d arrays !
! ================================================================= !

subroutine Forward_FFT(Nx, Ny, Nz, in, out)
implicit none
integer ( kind = 4 ) Nx,Ny,Nz,i,j,k
real ( kind = 8 ) in(Nx, Ny, Nz),out(Nx, Ny, Nz),dum(Nx*Ny*Nz)

do i = 1, Nx
  do j = 1, Ny
    do k = 1, Nz
    dum(((i-1)*Ny+(j-1))*Nz+k) = in(i,j,k)
    enddo
  enddo
enddo

  call dst3f(Nx, Ny, Nz, dum)

do i = 1, Nx
  do j = 1, Ny
    do k = 1, Nz
    out(i,j,k) = dum(((i-1)*Ny+(j-1))*Nz+k)
    enddo
  enddo
enddo

end subroutine

! ================================================================== !
! Wrapper Subroutine to call backward fftw c functions for 3d arrays !
! ================================================================== !

subroutine Backward_FFT(Nx, Ny, Nz, in, out)
implicit none
integer ( kind = 4 ) Nx,Ny,Nz,i,j,k
real ( kind = 8 ) in(Nx, Ny, Nz),out(Nx, Ny, Nz),dum(Nx*Ny*Nz)

do i = 1, Nx
  do j = 1, Ny
    do k = 1, Nz
    dum(((i-1)*Ny+(j-1))*Nz+k) = in(i,j,k)
    enddo
  enddo
enddo

  call dst3b(Nx, Ny, Nz, dum)

do i = 1, Nx
  do j = 1, Ny
    do k = 1, Nz
    out(i,j,k) = dum(((i-1)*Ny+(j-1))*Nz+k)
    enddo
  enddo
enddo

end subroutine
! ==================================================================

这段代码使用下面的 C-wrapper 来计算正向 3D FFTW 正弦变换和反向 3D FFTW 正弦变换，

#include <fftw3.h>

int dst3f_(int *n0, int *n1, int *n2, double *in3cs)
{
    double *out3cs;
    out3cs = (double*) fftw_malloc(sizeof(double) * (*n0) * (*n1) * (*n2));
    fftw_plan p3cs;
    p3cs = fftw_plan_r2r_3d(*n0, *n1, *n2, in3cs, out3cs, FFTW_RODFT10, FFTW_RODFT10, FFTW_RODFT10, FFTW_ESTIMATE);

    fftw_execute(p3cs);
    fftw_destroy_plan(p3cs);

    for(int i0=0;i0<*n0;i0++) {
        for(int i1=0;i1<*n1;i1++) {
            for(int i2=0;i2<*n2;i2++) {
                in3cs[(i0*(*n1)+i1)*(*n2)+i2] = out3cs[(i0*(*n1)+i1)*(*n2)+i2];
            }
        }
    }

    return 0;
}

int dst3b_(int *n0, int *n1, int *n2, double *in3cs)
{
    double *out3cs;
    out3cs = (double*) fftw_malloc(sizeof(double) * (*n0) * (*n1) * (*n2));
    fftw_plan p3cs;
    p3cs = fftw_plan_r2r_3d(*n0, *n1, *n2, in3cs, out3cs, FFTW_RODFT01, FFTW_RODFT01, FFTW_RODFT01, FFTW_ESTIMATE);

    fftw_execute(p3cs);
    fftw_destroy_plan(p3cs);

    for(int i0=0;i0<*n0;i0++) {
        for(int i1=0;i1<*n1;i1++) {
            for(int i2=0;i2<*n2;i2++) {
                in3cs[(i0*(*n1)+i1)*(*n2)+i2] = out3cs[(i0*(*n1)+i1)*(*n2)+i2];
            }
        }
    }

    return 0;
}

然后我尝试使用正弦变换 FFTW 求解相同的泊松方程，即 r2r(实数到实数)类型。当我绘制输出的 3-D 图时，如图所示 3 ，振幅现在减小到小于 1。我无法找出我在代码中的错误在哪里，因为在正弦变换的情况下这个振幅被减小了。

Answer 1

使用实到实变换来求解泊松方程非常有吸引力，因为它允许使用各种边界条件。然而，传感点与周期性边界条件所考虑的点并不完全对应。

对于周期性边界条件，传感点位于规则网格上，例如 0,1..,n-1。如果晶胞的大小是Lx，那么点之间的间距就是Lx/n。故事结束。

现在，让我们考虑使用 DST III 进行正向变换的边界条件，即标记 RODFT01。如图there , 它在 documentation of FFTW 中发出信号那:

FFTW_RODFT01 (DST-III): odd around j=-1 and even around j=n-1.

感应点还是0,1..,n-1。如果长度为Lx，则间距仍为dx=Lx/(n)。 但是DST III的输入函数和基函数在j=-1附近是奇数，在j=n-1附近是偶数。这解释了幅度上的差异:

in_dst(i,j,k) = 27.0d0*sin(3.0d0*x(i))*sin(3.0d0*y(j))*sin(3.0d0*z(k))

事实上，这个输入在 i=-1 附近不是奇数，在 i=n-1 附近是偶数。 i=0 和 i=n 都是奇怪的。因此，以下操作可能会有所帮助:

• 确保所考虑的正向变换是正确的。对于奇-奇边界条件，它可以是 FFTW_RODFT00 (DST-I):奇数在 j=-1 附近和奇数在 j=n 附近。。
• 评估正确传感点的输入。对于 DST-I，如果函数在 0 和 Lx 处为奇数，则使用 n 个点，对于 i=0,1,..,n-1，dx=Lx/(n+1) 和 x_i=dx*(i+1) .
• 计算变换空间中的频率。一路the inverse transform writes有助于找到这些频率。由于DST I的倒数是DST I，第一个频率(k=0)对应的周期为2Lx；第二个(k=1)对应Lx，第k-st个频率对应2Lx/(k+1)的周期。

最后，可能需要对输出进行缩放，因为 FFTW transforms are not normalized.对于 DST-I，FFTW_RODFT00，它是 N=2(n+1)。 VladimirF 的建议无疑是一个很好的建议。事实上，虽然测试单个频率是理解和实现算法的理想选择，但最终测试必须涵盖多个频率以确保程序可靠。