FFTW

Question

私はプラズマ物理学のために2Dでポアソン問題を解決しようとしています

フーリエ変換をdiscretによると、私はその後

を持っています私が得たポアソンシステムを解いた後に

と の4係数です。

私はfftwを使ってEを計算しようとしています。私はgcc poisson.c -lfftw3 -lm -o poissonでコンパイルされたファイル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; 
} 

を持っています。

入力は すべてが正常に動作しますが、入力は であるならば、それは、私は理由を理解していないしないでください。

Answer 1

私はfftw3のドキュメントの画像で誤解されました。ここに訂正があります。一つは、次のように宣言

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)]; 

を持つ必要がありますし、メインループは

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