使用带有矩形域的FFTW的泊松方程 [英] Poisson equation using FFTW with rectanguar domain

问题描述

我试图通过使用具有矩形域(-4 ＜＝ x ＜＝ 4和-2 ＜＝ y ＜＝ 2)的FFTW库来求解泊松方程.如果域是正方形，我有正确的结果，如果域是矩形，那是错误的.请给我一些建议.非常感谢.这是我的代码.

I am trying to solve Poisson equation by using FFTW library with rectangular domain(-4<=x<=4 and -2<=y<=2). I have correct result if domain is square and it is wrong if domain is rectangular. Please give me some suggestion. Thank you so much. Here is my code.

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

using namespace std;
int main() {

    int N1=64;
    int N2=32;

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

    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_REDFT00, FFTW_REDFT00, FFTW_EXHAUSTIVE);
    q = fftw_plan_r2r_2d(N1,N2, in2.data(), out2.data(), FFTW_REDFT00, FFTW_REDFT00, FFTW_EXHAUSTIVE);

    int l=-1;
    for(i = 0;i <N1;i++){
        X[i] =-4.0+(double)i*dx ;           
        for(j = 0;j<N2;j++){
            l=l+1;
            Y[j] =-2.0+ (double)j*dy ;
            in1[l]= cos(pi*X[i]) + cos(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);
            }else{
                fact=((double)(N1-i)*(N1-i));
            }
            if(2*j<N2){
                fact+=((double)j*j);
            }else{
                fact+=((double)(N2-j)*(N2-j));
            }
            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 += fabs( in1[l]-  8.*out2[l]/((double)(N1-1))/((double)(N2-1))); 
            printf("%3d %10.5f %10.5f\n", l, in1[l], 8.*out2[l]/((double)(N1-1))/((double)(N2-1)));

        }
    }

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



    return 0;
}

推荐答案

在傅立叶空间中，频率k_x和k_y必须取决于域的大小.由于域是矩形的，因此变得很重要.

In the Fourier space, the frequencies k_x and k_y must depend on the size of the domain. As the domain is rectangular, it becomes a matter of importance.

尝试:

double invL1s=1.0/(L1*L1);
double invL2s=1.0/(L2*L2);
...
        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;
        }

输出应接近预期的输出(比例因子除外).

The output should be closer to the expected one (except for a scale factor).

这篇关于使用带有矩形域的FFTW的泊松方程的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！