泊松方程光谱法

Question

我写了泊松方程。谱法求解器。然而，所得结果与周期性边界条件差分法的结果并不相符。 我认为我使用 FFTW 是错误的。 你能告诉我以下代码的哪一部分有错误吗？ 谢谢。

program main
  implicit none
  include 'fftw3.f'
  integer(8) :: plan
  integer, parameter :: j_max = 100, k_max = 100, m_max = j_max/2 + 1, n_max = k_max
  integer :: j, k, m, n, mm, nn
  real(8) :: v(1:j_max, 1:k_max), f(1:j_max, 1:k_max)
  real(8) :: x_max, y_max, dx, dy, x, y, t_max, pi
  complex(8), parameter :: im = (0.d0, 1.d0)
  complex(8) :: vk(1:m_max, 1:n_max), fk(1:m_max, 1:n_max)

  pi = 4.d0*atan(1.d0)
  x_max = 2.d0*pi
  y_max = 2.d0*pi
  dx = x_max/j_max
  dy = y_max/k_max

!*-- Initial Condition ---
  do j = 1, j_max
    x = dx*j
    do k = 1, k_max
      y = dy*k
      f(j, k) = dexp(-(x - x_max/2)**2 -(y - y_max/2)**2)
    enddo
  enddo

!*-- FFT forward ---
  call dfftw_plan_dft_r2c_2d(plan, j_max, k_max, v, vk, FFTW_ESTIMATE)
  call dfftw_execute(plan)
  call dfftw_plan_dft_r2c_2d(plan, j_max, k_max, f, fk, FFTW_ESTIMATE)
  call dfftw_execute(plan)

  do m = 1, m_max
    do n = 1, n_max
      if(m <= m_max/2 + 1) then
        mm = m - 1
      else
        mm = m - 1 - m_max
      endif
      if(n <= n_max/2 + 1) then
        nn = n - 1
      else
        nn = n - 1 - n_max
      endif

      if(mm == 0 .and. nn == 0) then
      else 
        vk(m, n) = fk(m, n)/(mm**2 + nn**2)
      endif
    enddo
      enddo

!*-- FFT backward ---
  call dfftw_plan_dft_c2r_2d(plan, j_max, k_max, vk, v, FFTW_ESTIMATE)
  call dfftw_execute(plan)

!*-- normalization ---
  v = v/j_max/k_max

  call dfftw_destroy_plan(plan)

end program main

Answer 1

这里你有一个代码可以满足你的要求，考虑到原始数据在我的例子中是“f1”和“f2”，重要的注释是英语的，其他一些是西班牙语的，如果你有理解上的问题告诉我吧:)

    // FFT CALCULATION
    // Inicialización de elementos necesarios para el cálculo de la FFT
    fftw_plan p1; // variable para almacenar la planificación de la FFT
    fftw_plan p2; // variable para almacenar la planificación de la FFT
    int N_fft= ancho*alto; //number of points of the image
    fftw_complex *U1 =(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*alto*((ancho/2)+1)); //puntero que apuntará al resultado de la FFT
    fftw_complex *U2 =(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*alto*((ancho/2)+1));
    p1 = fftw_plan_dft_r2c_2d(alto,ancho, f1, U1, FFTW_ESTIMATE); // FFT planning
    p2 = fftw_plan_dft_r2c_2d(alto,ancho, f2, U2, FFTW_ESTIMATE); // FFT planning
    fftw_execute(p1); // FFT calculation
    fftw_execute(p2); // FFT calculation
    fftw_destroy_plan(p1);// Eliminación de la planificación de la FFT
    fftw_destroy_plan(p2);// Eliminación de la planificación de la FFT

    // Security saving of U1 and U2 in auxiliar variables because the ifft modifies the input data
    for (int y = 0 ; y < alto ; y++){
        for (int x = 0 ; x < (ancho/2)+1 ; x++){
            U1_input_save[((ancho/2)+1)*y+x][0] = U1[((ancho/2)+1)*y+x][0];
            U1_input_save[((ancho/2)+1)*y+x][1] = U1[((ancho/2)+1)*y+x][1];
            U2_input_save[((ancho/2)+1)*y+x][0] = U2[((ancho/2)+1)*y+x][0];
            U2_input_save[((ancho/2)+1)*y+x][1] = U2[((ancho/2)+1)*y+x][1];
        }
    }

    // IFFT ( U1,U2 --> u1,u2)
    //----IFFT-----
    double *u1 = (double*) malloc(sizeof(double)*N_fft);//puntero que apuntará al resultado de la IFFT 
    double *u2 = (double*) malloc(sizeof(double)*N_fft);
    fftw_plan p3;// variable para almacenar la planificación de la IFFT
    fftw_plan p4;// variable para almacenar la planificación de la IFFT

    p3 = fftw_plan_dft_c2r_2d(alto, ancho, U1, u1, FFTW_ESTIMATE);//planificación de la fft inversa
    p4 = fftw_plan_dft_c2r_2d(alto, ancho, U2, u2, FFTW_ESTIMATE);//planificación de la fft inversa
    fftw_execute(p3); // Calculo de la fft inversa
    fftw_execute(p4); // Calculo de la fft inversa 
    fftw_destroy_plan(p3); // Eliminación de la planificación de la IFFT
    fftw_destroy_plan(p4); // Eliminación de la planificación de la IFFT


    // Normalización after IFFT important!
    u1 = fftw_normalization(ancho,alto,N_fft,u1);
    u2 = fftw_normalization(ancho,alto,N_fft,u2);

    // Correction of U1 and U2, restoring the original data 
    for (int y = 0 ; y < alto ; y++){
        for (int x = 0 ; x < (ancho/2)+1 ; x++){
            U1[((ancho/2)+1)*y+x][0] = U1_input_save[((ancho/2)+1)*y+x][0];
            U1[((ancho/2)+1)*y+x][1] = U1_input_save[((ancho/2)+1)*y+x][1];
            U2[((ancho/2)+1)*y+x][0] = U2_input_save[((ancho/2)+1)*y+x][0];
            U2[((ancho/2)+1)*y+x][1] = U2_input_save[((ancho/2)+1)*y+x][1];
        }
    }

    // FIN CALCULATION FFT

Answer 2

在FFT正向过程中，'2d_c2r'函数可以修改输入值，然后，如果您稍后使用then，结果将不正确，您可以在执行该函数之前复制数据。