fftw pour Poisson en physique des plasmas

Question

Je suis en train de résoudre le problème de Poisson en 2D pour la physique des plasmas

Selon le discret Fourrier Transform, je

Ensuite, après avoir résolu le système de Poisson, j'obtiens

avec les coefficients de Fourrier.

J'essaie de calculer E en utilisant fftw. Je le fichier poisson.c

#include <complex.h> 
#include <fftw3.h> 
#include <stdlib.h> 
#include <stdio.h> 
#include <math.h> 

#define PI (double)(3.1415926535897932384626433832795028841971693993751058) //from Wolfram Alpha 

#define SQ(var) ((var)*(var)) 

#define PTS_PER_DIR 32 // power of 2 
#define PTS_2D 1024 // PTS_PER_DIR*PTS_PER_DIR 

int main() 
{ 
    double rho_in[PTS_2D], ex_out[PTS_2D], ey_out[PTS_2D]; 
    fftw_complex rho_out[PTS_2D], ex_in[PTS_2D], ey_in[PTS_2D]; 

    double posX[PTS_PER_DIR], posY[PTS_PER_DIR]; 
    double lx=2.*PI, ly=2.*PI; 
    double dx=lx/(double)(PTS_PER_DIR), dy=ly/(double)(PTS_PER_DIR); 

    fftw_plan forward, backward_ex, backward_ey; 

    forward=fftw_plan_dft_r2c_2d(PTS_PER_DIR, PTS_PER_DIR, rho_in, rho_out, 
           FFTW_EXHAUSTIVE | FFTW_PRESERVE_INPUT); // need to keep rho in full code 
    backward_ex=fftw_plan_dft_c2r_2d(PTS_PER_DIR, PTS_PER_DIR, ex_in, ex_out, FFTW_EXHAUSTIVE); 
    backward_ey=fftw_plan_dft_c2r_2d(PTS_PER_DIR, PTS_PER_DIR, ey_in, ey_out, FFTW_EXHAUSTIVE); 

    int i1, i2, k1, k2; 

    for (i1=0; i1<PTS_PER_DIR; ++i1) { 
    posX[i1]=(double)(i1)*dx; 
    posY[i1]=(double)(i1)*dy; 
    } 

    for (i1=0; i1<PTS_PER_DIR; ++i1) 
    for (i2=0; i2<PTS_PER_DIR; ++i2) 
     rho_in[i1*PTS_PER_DIR+i2]=sin(posY[i2]); 

    fftw_execute(forward); 

    for (i1=0; i1<PTS_PER_DIR; ++i1) { 
    k1=i1; 
    if (i1>PTS_PER_DIR/2) 
     k1-=PTS_PER_DIR; 

    for (i2=0; i2<PTS_PER_DIR; ++i2) { 
     k2=i2; 
     if (i2>PTS_PER_DIR/2) 
     k2-=PTS_PER_DIR; 

     if (i1==0 && i2==0) { 
     ex_in[0]=0; 
     ey_in[0]=0; 

     continue; 
     } 

     ex_in[i1*PTS_PER_DIR+i2]=rho_out[i1*PTS_PER_DIR+i2]*(double)(k1)*I/(2.*PI*lx)/
          (SQ((double)(k1)/lx)+SQ((double)(k2)/ly)); 
     ey_in[i1*PTS_PER_DIR+i2]=rho_out[i1*PTS_PER_DIR+i2]*(double)(k2)*I/(2.*PI*ly)/
          (SQ((double)(k1)/lx)+SQ((double)(k2)/ly)); 
    } // for i2 
    } // for i1 

    fftw_execute(backward_ex); 
    fftw_execute(backward_ey); 

    for (i1=0; i1<PTS_2D; ++i1) { 
    ex_out[i1]/=(double)(PTS_2D); 
    ey_out[i1]/=(double)(PTS_2D); 
    } 

    for (i1=0; i1<PTS_PER_DIR; ++i1) { 
    for (i2=0; i2<PTS_PER_DIR; ++i2) 
     printf("%lg %lg %lg %lg %lg\n", posX[i1], posY[i2], 
      rho_in[i1*PTS_PER_DIR+i2], ex_out[i1*PTS_PER_DIR+i2], ey_out[i1*PTS_PER_DIR+i2]); 
    printf("\n"); 
    } 

    return 0; 
} 

compilé avec gcc poisson.c -lfftw3 -lm -o poisson.

Si l'entrée est tout travail bien, mais si l'entrée est alors il ne fonctionne pas et je ne comprends pas pourquoi.

Answer 1

J'ai été trompé par une image sur la documentation de fftw3. Voici les corrections. On devrait avoir la déclaration suivante

fftw_complex rho_out[PTS_PER_DIR*(PTS_PER_DIR/2+1)], 
      ex_in[PTS_PER_DIR*(PTS_PER_DIR/2+1)], ey_in[PTS_PER_DIR*(PTS_PER_DIR/2+1)]; 

et la boucle principale devient

for (i1=0; i1<PTS_PER_DIR; ++i1) { 
    k1=i1; 
    if (i1>PTS_PER_DIR/2) 
    k1-=PTS_PER_DIR; 

    for (i2=0; i2<=PTS_PER_DIR/2; ++i2) { 
    k2=i2; 

    if (i1==0 && i2==0) { 
     ex_in[0]=0; 
     ey_in[0]=0; 

     continue; 
    } 

    ex_in[i1*(PTS_PER_DIR/2+1)+i2]=rho_out[i1*(PTS_PER_DIR/2+1)+i2]*(double)(k1)*I/(2.*PI*lx)/
          (SQ((double)(k1)/lx)+SQ((double)(k2)/ly)); 
    ey_in[i1*(PTS_PER_DIR/2+1)+i2]=rho_out[i1*(PTS_PER_DIR/2+1)+i2]*(double)(k2)*I/(2.*PI*ly)/
          (SQ((double)(k1)/lx)+SQ((double)(k2)/ly)); 
    } // for i2 
} // for i1