poisson eq。用譜法

Question

我寫了poisson eq。用譜方法求解。然而，得到的結果與週期性邊界條件差分法的結果並不一致。 我想我錯誤地使用了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

在FFT前進過程中的「2d_c2r」功能可以修改輸入值，然後，如果你使用再到後來的結果將是不正確的，你可以執行該功能之前做數據的副本。

希望幫助

安東尼

Answer 2

在這裏，你有你所求的，考慮到原始數據是在我的情況下，「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 

如果你需要更多的東西只是告訴我！

Antonio