fftw3 pour poisson avec condition de limite de dirichlet pour tout le côté du domaine de calcul

Question

J'essaye de résoudre l'équation de poison avec la condition de frontière de Dirichlet pour quatre côtés du domaine de calcul. Comme je le sais que je devrais utiliser FFTW_RODFT00 pour satisfaire la condition. Cependant, le résultat n'est pas correct. Pourriez-vous m'aider s'il vous plaît?

#include <stdio.h> 
#include <math.h> 
#include <cmath> 
#include <fftw3.h> 
#include <iostream> 
#include <vector> 

using namespace std; 
int main() { 

int N1=100; 
int N2=100; 

double pi = 3.141592653589793; 
double L1 = 2.0; 
double dx = L1/(double)(N1-1); 
double L2= 2.0; 
double dy=L2/(double)(N2-1); 

double invL1s=1.0/(L1*L1); 
double invL2s=1.0/(L2*L2); 

std::vector<double> in1(N1*N2,0.0); 
std::vector<double> in2(N1*N2,0.0); 
std::vector<double> out1(N1*N2,0.0); 
std::vector<double> out2(N1*N2,0.0); 
std::vector<double> X(N1,0.0); 
std::vector<double> Y(N2,0.0); 


fftw_plan p, q; 
int i,j; 
p = fftw_plan_r2r_2d(N1,N2, in1.data(), out1.data(), FFTW_RODFT00, FFTW_RODFT00, FFTW_EXHAUSTIVE); 
q = fftw_plan_r2r_2d(N1,N2, in2.data(), out2.data(), FFTW_RODFT00, FFTW_RODFT00, FFTW_EXHAUSTIVE); 

int l=-1; 
for(i = 0;i <N1;i++){ 
    X[i] =-1.0+(double)i*dx ;   
    for(j = 0;j<N2;j++){ 
     l=l+1; 
     Y[j] =-1.0+ (double)j*dy ; 
     in1[l]= sin(pi*X[i]) + sin(pi*Y[j]) ; // row major ordering 
    } 
} 

fftw_execute(p); 

l=-1; 
for (i = 0; i < N1; i++){ // f = g/(kx² + ky²) 
    for(j = 0; j < N2; j++){ 
      l=l+1; 
     double fact=0; 
     in2[l]=0; 

     if(2*i<N1){ 
      fact=((double)i*i)*invL1s;; 
     }else{ 
      fact=((double)(N1-i)*(N1-i))*invL1s; 
     } 
     if(2*j<N2){ 
      fact+=((double)j*j)*invL2s; 
     }else{ 
      fact+=((double)(N2-j)*(N2-j))*invL2s; 
     } 
     if(fact!=0){ 
      in2[l] = out1[l]/fact; 
     }else{ 
      in2[l] = 0.0; 
     } 
    } 
} 

fftw_execute(q); 
l=-1; 
double erl1 = 0.; 
for (i = 0; i < N1; i++) { 
    for(j = 0; j < N2; j++){ 
     l=l+1; 

     erl1 +=1.0/pi/pi*fabs(in1[l]- 0.25*out2[l]/((double)(N1-1))/((double)(N2-1))); 
     printf("%3d %10.5f %10.5f\n", l, in1[l], 0.25*out2[l]/((double)(N1-1))/((double)(N2-1))); 

    } 
} 

cout<<"error=" <<erl1 <<endl ; 
fftw_destroy_plan(p); fftw_destroy_plan(q); fftw_cleanup(); 

return 0; 

}

Answer 1

Je reconnais un truc je vous ai fourni dans Poisson equation using FFTW with rectanguar domain et le code et dans ma réponse à Confusion testing fftw3 - poisson equation 2d test, qui a été adapté à partir du code de la @Charles_P ASKER! S'il vous plaît envisager d'ajouter des liens à ces URL dans les prochaines questions!

La réponse à Confusion testing fftw3 - poisson equation 2d test a été consacrée au cas des conditions aux limites périodiques. Voici donc quelques modifications pour résoudre le cas des conditions aux limites de Dirichlet.

Le fftw_plan_r2r_2d(N1,N2, in1.data(), out1.data(), FFTW_RODFT00, FFTW_RODFT00,FFTW_EXHAUSTIVE) correspond à un type I sine transformée discrète de as defined by the FFTW library:

Il est sens bien décrit dans https://en.wikipedia.org/wiki/Discrete_sine_transform. Si la taille du tableau FFTW est N1=4 et ses valeurs [a, b, c, d], le tableau complet incluant les limites est [0, a, b, c, d, 0]. D'où l'étape spatiale est:

et les fréquences f_k du type I DST sont:

L'inverse du type I DST est le type I DST, à l'exception de un facteur d'échelle (voir http://www.fftw.org/doc/1d-Real_002dodd-DFTs-_0028DSTs_0029.html#g_t1d-Real_002dodd-DFTs-_0028DSTs_0029), ici 4.(N1+1).(N2+1). Enfin, le cas d'essai doit être adapté au cas des conditions aux limites de Dirichlet. En effet, sur la case de taille L1,L2 la fonction ne respecte pas les conditions aux limites de Diriclet. En effet, même si le terme source est le même, la solution satisfaisant aux conditions aux limites périodiques peut être différente de la solution qui satisfait les conditions aux limites de Dirichelt. Au lieu de cela, deux termes sources peuvent être testées:

• Le terme source correspond à une seule fréquence de la DST.

• Le terme source est directement dérivé de la solution

Enfin, voici un code résolution de l'équation de Poisson 2D à l'aide du type I DST de la bibliothèque FFTW:

#include <stdio.h> 
#include <math.h> 
#include <cmath> 
#include <fftw3.h> 
#include <iostream> 
#include <vector> 

using namespace std; 
int main() { 

    int N1=100; 
    int N2=200; 

    double pi = 3.141592653589793; 
    double L1 = 1.0; 
    double dx = L1/(double)(N1+1);//+ instead of -1 
    double L2= 5.0; 
    double dy=L2/(double)(N2+1); 

    double invL1s=1.0/(L1*L1); 
    double invL2s=1.0/(L2*L2); 

    std::vector<double> in1(N1*N2,0.0); 
    std::vector<double> in2(N1*N2,0.0); 
    std::vector<double> out1(N1*N2,0.0); 
    std::vector<double> out2(N1*N2,0.0); 
    std::vector<double> X(N1,0.0); 
    std::vector<double> Y(N2,0.0); 


    fftw_plan p, q; 
    int i,j; 
    p = fftw_plan_r2r_2d(N1,N2, in1.data(), out1.data(), FFTW_RODFT00, FFTW_RODFT00, FFTW_EXHAUSTIVE); 
    q = fftw_plan_r2r_2d(N1,N2, in2.data(), out2.data(), FFTW_RODFT00, FFTW_RODFT00, FFTW_EXHAUSTIVE); 

    int l=0; 

    for(i = 0;i <N1;i++){ 
     for(j = 0;j<N2;j++){ 
      X[i] =dx+(double)i*dx ;  
      Y[j] =dy+ (double)j*dy ; 
      //in1[l]= sin(pi*X[i])*sin(pi*Y[j]) ; // row major ordering 
      in1[l]=2*Y[j]*(L2-Y[j])+2*X[i]*(L1-X[i]); 
      l=l+1; 
     } 
    } 

    fftw_execute(p); 

    l=-1; 
    for (i = 0; i < N1; i++){ // f = g/(kx² + ky²) 
     for(j = 0; j < N2; j++){ 

      l=l+1; 
      double fact=0; 

      fact=pi*pi*((double)(i+1)*(i+1))*invL1s; 

      fact+=pi*pi*((double)(j+1)*(j+1))*invL2s; 

      in2[l] = out1[l]/fact; 

     } 
    } 

    fftw_execute(q); 
    l=-1; 
    double erl1 = 0.; 
    for (i = 0; i < N1; i++) { 
     for(j = 0; j < N2; j++){ 
      l=l+1; 
      X[i] =dx+(double)i*dx ;  
      Y[j] =dy+ (double)j*dy ; 
      //double res=0.5/pi/pi*in1[l]; 
      double res=X[i]*(L1-X[i])*Y[j]*(L2-Y[j]); 
      erl1 +=pow(fabs(res- 0.25*out2[l]/((double)(N1+1))/((double)(N2+1))),2); 
      printf("%3d %10.5g %10.5g\n", l, res, 0.25*out2[l]/((double)(N1+1))/((double)(N2+1))); 

     } 
    } 
    erl1=erl1/((double)N1*N2); 
    cout<<"error=" <<sqrt(erl1) <<endl ; 
    fftw_destroy_plan(p); fftw_destroy_plan(q); fftw_cleanup(); 

    return 0; 
} 

Il est compilé par g++ main.cpp -o main -lfftw3 -Wall.

EDIT: Comment calculer la réponse à chaque fréquence?

L'idée de base de FFT est de représenter la solution sous forme d'une combinaison linéaire de fonctions:

• Dans le cas de conditions limites périodiques, la FFT est utilisée. Les fonctions de base sont:

• Dans le cas des conditions de Dirichlet limites, est utilisé le type I DST. Les fonctions de base, étant de 0 à x=0 et x=L1, sont les suivants:

  

• Dans le cas des conditions aux limites de Neumann, on utilise la DCT de type-I. Les fonctions de base sont les suivants:

  

leurs dérivés sont nulles à x=0 et x=L1.

Calculons la dérivée seconde de la f_k componant de type-I DST:

Par conséquent, la componant k de l'heure d'été de la solution correspond à la componant k du DST du terme source, mis à l'échelle d'un facteur