fortran 2d-FFTW 与 C 2d-FFTW 结果不一致

Question

我正在学习如何使用 fortran 处理 FFTW 包。为了产生一个易于验证的例子，我计算了一个二维平面的功率谱，我用两个不同的叠加波填充了它。这样，我就可以准确地知道功率谱中预期峰值的位置。

由于 FFTW 文档对 C 的内容更广泛，我首先在 C 中实现了该算法，这给了我相当满意的结果： “lambda1”和“lambda2”是已知的预期波长。蓝线是得到的功率谱。

然后我尝试对 fortran 做同样的事情，结果很奇怪：

我不知道去哪里寻找可能的错误。代码执行顺利。有人可以帮忙吗？

这是 C 代码，编译为：gcc stackexchange.c -o a.out -I/home/myname/.local/include -lfftw3 -lm -g -Wall (gcc 5.4.0)

#include <fftw3.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

int main(void)
{
    // parameters
    int Nx = 200;
    int Ny = 100;            
    int nsamples = 200;
    float pi = 3.1415926;
    float physical_length_x = 20;
    float physical_length_y = 10;
    float lambda1 = 0.5;
    float lambda2 = 0.7;
    float dx = physical_length_x/Nx;
    float dy = physical_length_y/Ny;
    float dkx = 2*pi/physical_length_x;
    float dky = 2*pi/physical_length_y;

    // counters/iterators
    int ind, i, j, ix, iy, ik;

    // power spectra stuff
    float *Pk, *distances_k, d, kmax;
    float *Pk_field;


    // FFTW vars and arrays
    fftw_complex *in, *out;
    fftw_plan my_plan;

    // allocate arrays for input/output
    in  = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * Nx * Ny);
    out = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * Nx * Ny);

    // Create Plan
    int n[2]; n[0] = Nx; n[1] = Ny;
    my_plan = fftw_plan_dft(2, n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

    // fill up array with a wave
    for (i=0; i<Nx; i++){
      for (j=0; j<Ny; j++){
        ind = i*Ny + j;
        in[ind][0] = cos(2.0*pi/lambda1*i*dx) + sin(2.0*pi/lambda2*j*dy);
      }
    }

    // execute fft
    fftw_execute(my_plan); /* repeat as needed */

    // Calculate power spectrum: P(kx, ky) = |F[delta(x,y)]|^2
    Pk_field = malloc(sizeof(float)*Nx*Ny);
    for (i=0; i<Nx; i++){
      for (j=0; j<Ny; j++){
        ind = i*Ny+j;
        Pk_field[ind] = out[ind][0]*out[ind][0] + out[ind][1]*out[ind][1];
      }
    }

    Pk = calloc(nsamples, sizeof(float));
    distances_k = malloc(nsamples*sizeof(float));
    kmax = sqrt(pow((Nx/2+1)*dkx, 2) + pow((Ny/2+1)*dky, 2));
    for(i=0; i<nsamples; i++){
      distances_k[i]= 1.0001*i/nsamples*kmax; // add a little more to make sure kmax will fit
    }

    // histogrammize P(|k|)
    for (i=0; i<Nx; i++){

      if (i<Nx/2+1)
        { ix = i; }
      else
        { ix = -Nx+i; }

      for (j=0; j<Ny; j++){

        if (j<Ny/2+1)
          { iy = j; }
        else
          { iy = -Ny+j; }

        ind = i*Ny + j;
        d = sqrt(pow(ix*dkx,2)+pow(iy*dky,2));

        for(ik=0; ik<nsamples; ik++){
          if (d<=distances_k[ik] || ik==nsamples){
            break;
          }
        }

        Pk[ik] += Pk_field[ind];
      }
    }


    //-----------------------------------
    // write arrays to file.
    // can plot them with plot_fftw.py
    //-----------------------------------
    FILE *filep;
    filep = fopen("./fftw_output_2d_complex.txt", "w");
    for (i=0; i<nsamples; i++){
      fprintf(filep, "%f\t%f\n", distances_k[i], Pk[i]);
    }
    fclose(filep);

    printf("Finished! Written results to ./fftw_output_2d_complex.txt\n");


    //----------------------
    // deallocate arrays
    //----------------------
    fftw_destroy_plan(my_plan);
    fftw_free(in); fftw_free(out);

    return 0;
}

这是 Fortran 代码，编译时使用：gfortran stackexchange.f03 -o a.out -L/home/my_name/.local/lib/ -I/home/my_name/.local/include/ -lfftw3 -g -Wall

program use_fftw

  use,intrinsic :: iso_c_binding
  implicit none
  include 'fftw3.f03'

  integer, parameter     :: dp=kind(1.0d0)
  integer, parameter     :: Nx = 200
  integer, parameter     :: Ny = 100
  integer, parameter     :: nsamples = 200
  real(dp), parameter    :: pi = 3.1415926d0
  real(dp), parameter    :: physical_length_x = 20.d0
  real(dp), parameter    :: physical_length_y = 10.d0
  real(dp), parameter    :: lambda1 = 0.5d0
  real(dp), parameter    :: lambda2 = 0.7d0
  real(dp), parameter    :: dx = physical_length_x/real(Nx,dp)
  real(dp), parameter    :: dy = physical_length_y/real(Ny,dp)
  real(dp), parameter    :: dkx = 2.d0 * pi / physical_length_x
  real(dp), parameter    :: dky = 2.d0 * pi / physical_length_y

  integer :: i, j, ix, iy, ik
  real(dp):: kmax, d

  complex(dp), allocatable, dimension(:,:)    :: arr_in, arr_out, Pk_field
  real(dp), allocatable, dimension(:)         :: Pk, distances_k
  integer*8                                   :: my_plan
  integer                                     :: n(2) = (/Nx, Ny/)


  allocate(arr_in( 1:Nx,1:Ny))
  allocate(arr_out(1:Nx,1:Ny))

  ! call dfftw_plan_dft_2d(my_plan, Nx, Ny, arr_in, arr_out, FFTW_FORWARD, FFTW_ESTIMATE) ! doesn't help either
  call dfftw_plan_dft(my_plan, 2, n, arr_in, arr_out, FFTW_FORWARD, FFTW_ESTIMATE)

  ! fill up wave
  do i = 1, Nx
    do j = 1, Ny
      arr_in(i,j) =cmplx(cos(2.0*pi/lambda1*i*dx)+sin(2.0*pi/lambda2*j*dy) , 0.d0, kind=dp)
    enddo
  enddo

  ! execute fft
  call dfftw_execute_dft(my_plan, arr_in, arr_out)


  allocate(Pk_field(1:Nx, 1:Ny))
  allocate(distances_k(1:nsamples))
  allocate(Pk(1:nsamples))
  Pk=0
  distances_k=0

  ! Get bin distances
  kmax = sqrt(((Nx/2+1)*dkx)**2+((Ny/2+1)*dky)**2)*1.001
  do i = 1, nsamples
    distances_k(i) = i*kmax/nsamples
  enddo

  ! Compute P(k) field, distances field
  do i = 1, Nx
    do j = 1, Ny
      Pk_field(i,j) = arr_out(i,j)*conjg(arr_out(i,j))

      if (i<=Nx/2+1) then
        ix = i
      else
        ix = -Nx+i
      endif
      if (j<=Ny/2+1) then
        iy = j
      else
        iy = -Ny+j
      endif

      d = sqrt((ix*dkx)**2+(iy*dky)**2)

      do ik=1, nsamples
        if (d<=distances_k(ik) .or. ik==nsamples) exit
      enddo

      Pk(ik) = Pk(ik)+real(Pk_field(i,j))
    enddo
  enddo

  ! write file
  open(unit=666,file='./fftw_output_2d_complex.txt', form='formatted')
  do i = 1, nsamples
    write(666, '(2E14.5,x)') distances_k(i), Pk(i)
  enddo
  close(666)


  deallocate(arr_in, arr_out, Pk, Pk_field, distances_k)
  call dfftw_destroy_plan(my_plan)

  write(*,*) "Finished! Written results to fftw_output_2d_complex.txt"

end program use_fftw

为方便起见，这是我用来绘制结果的简短 Python 脚本：

#!/usr/bin/python3

#====================================
# Plots the results of the FFTW
# example programs.
#====================================

import numpy as np
import matplotlib.pyplot as plt
from sys import argv
from time import sleep


errormessage="""
I require an argument: Which output file to plot.
Usage: ./plot_fftw.py <case>
options for case:
    1   fftw_1d_complex.txt
    2   fftw_2d_complex.txt
    3   fftw_1d_real.txt
    4   fftw_2d_real.txt

Please select a case: """



#----------------------
# Hardcoded stuff
#----------------------

file_dict={}
file_dict['1'] = ('fftw_output_1d_complex.txt', '1d complex fftw')
file_dict['2'] = ('fftw_output_2d_complex.txt', '2d complex fftw')
file_dict['3'] = ('fftw_output_1d_real.txt', '1d real fftw')
file_dict['4'] = ('fftw_output_2d_real.txt', '2d real fftw')

lambda1=0.5
lambda2=0.7




#------------------------
# Get case from cmdline
#------------------------

case = ''

def enforce_integer():
    global case
    while True:
        case = input(errormessage)
        try:
            int(case)
            break
        except ValueError:
            print("\n\n!!! Error: Case must be an integer !!!\n\n")
            sleep(2)


if len(argv) != 2:
    enforce_integer()
else:
    try:
        int(argv[1])
        case = argv[1]
    except ValueError:
        enforce_integer()


filename,title=file_dict[case]




#-------------------------------
# Read and plot data
#-------------------------------

k, Pk = np.loadtxt(filename, dtype=float, unpack=True)

fig = plt.figure()

ax = fig.add_subplot(111)
ax.set_title("Power spectrum for "+title)
ax.set_xlabel("k")
ax.set_ylabel("P(k)")
#  ax.plot(k, Pk, label='power spectrum')
ax.semilogx(k[k>0], Pk[k>0], label='power spectrum') # ignore negative k
ax.plot([2*np.pi/lambda1]*2, [Pk.min()-1, Pk.max()+1], label='expected lambda1')
ax.plot([2*np.pi/lambda2]*2, [Pk.min()-1, Pk.max()+1], label='expected lambda2')
ax.legend()

plt.show()

Answer 1

在计算直方图时，循环的第一个索引对应平均值，i。 e. “零频率”。

在 C 中，当 i=0、j=0、ix=0、iy=0 和 d=0 时。那是正确的。对于fortran中数组的同一项，索引i=1，然后ix=1，d不再为空。

这是针对零频率的，但由于 C 和 Fortran 之间的索引差异，整个直方图受到轻微污染。特别是，修改可能会以不同的方式影响正频率和负频率，从而触发图中观察到的双峰。

您能否尝试将 Fortran 中 ix 和 iy 的计算修改为：

  if (i-1<Nx/2+1) then
    ix = i-1
  else
    ix = -Nx+i-1
  endif
  if (j-1<Ny/2+1) then
    iy = j-1
  else
    iy = -Ny+j-1
  endif