poisson eq. with spectral method

Question

I wrote poisson eq. solver with spectral method. However, the obtained result does not coincide with the result of difference method with periodic boundary condition. I think I am mistaken in the use of FFTW. Could you tell me which part of the following code contains errors? Thank you.

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

Here you have a code that do what you ask for, take into account that the original data was in 'f1' and 'f2' in my case, the important comments are in english, some other in spanish, if you have problems to understand just tell me :)

    // 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

In the FFT Forward process the '2d_c2r' function can modify the input values, then, if you use then later the results will not be correct, you can do a copy of the data before executing that function.