poisson éq. avec méthode spectrale

Question

J'ai écrit poisson eq. solveur avec la méthode spectrale. Cependant, le résultat obtenu ne coïncide pas avec le résultat de la méthode des différences avec condition aux limites périodiques. Je pense que je me trompe dans l'utilisation de FFTW. Pourriez-vous me dire quelle partie du code suivant contient des erreurs? Merci.

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

Dans la FFT processus avant la fonction « 2d_c2r » peut modifier les valeurs d'entrée, puis, si vous utilisez ensuite plus tard, les résultats ne seront pas corrects, vous pouvez faire une copie des données avant d'exécuter cette fonction.

espoir qui aide

Antonio

Answer 2

Ici vous avez un code qui faites ce que vous demandez, prendre en compte que les données d'origine était « f1 » et « f2 » dans mon cas, l'importante les commentaires sont en anglais, en espagnol un autre, si vous avez des problèmes à comprendre juste me dire :)

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

Si vous avez besoin quelque chose de plus me dire!

Antonio