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

【问题标题】：Poisson equation solver in 3-D Sine transform of r2r type using FFTW in Fortran在 Fortran 中使用 FFTW 的 r2r 类型的 3-D 正弦变换中的泊松方程求解器
【发布时间】：2018-08-29 21:52:05
【问题描述】：

我在 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 所示进行输出的 3-D 绘图时，幅度现在减小到小于 1。我无法找出代码中的错误在哪里，因为在以下情况下幅度会减小正弦变换。

【问题讨论】：

欢迎，请始终使用标签fortran。请注意，kind=4 和 kind=8 中的幻数 4 和 8 很丑陋且不可移植。我做了很多类似的事情，而且我的特征值大多是错误的。仔细检查它们，尝试不同的波长输入是如何受到影响的。确保哪个波长存储在转换后的数组中的哪个数组索引中。
另外，你对c2r-r2c版本的测试是不完整的，实际上是很不完整的。您必须测试存在多个波数的右侧。不仅仅是一个简单的正弦波。一个正弦波表明你有一个正确的特征值。您必须测试其中的几个。

标签： c fortran fft fftw dft

【解决方案1】：

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

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

现在，让我们考虑将 DST III 用于正向变换的边界条件，即标志 RODFT01。如图there，在documentation of FFTW 中表示：

FFTW_RODFT01 (DST-III)：在 j=-1 左右奇数，在 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): odd around j=-1 and odd around j=n.。
在正确的感应点评估输入。对于 DST-I，如果函数在 0 和 Lx 处为奇数，则使用 n 个点，dx=Lx/(n+1) 和 x_i=dx*(i+1) 对于 i=0,1,..,n-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 的建议无疑是一个很好的建议。事实上，虽然测试单个频率是理解和实现算法的理想选择，但最终测试必须覆盖多个频率以确保程序可靠。

【讨论】：