2D fissw poisson solver в различных размерах - PullRequest
Новая
       37

2D fissw poisson solver в различных размерах

1 голос
/ 23 января 2020

Я пытаюсь расширить решатель Пуассона, используя fft, предоставленную в Тестирование путаницы fftw3 - тест 2d уравнения Пуассона , для различных размеров L, так как оригинальный автор и ответ работают только с L = 2pi. Основное изменение происходит при f = g / (kx² + ky²), где kx теперь равно i * 2pi / L или (Ni) * 2pi / L. Это частота Nquist.

Аналогичные изменения сделаны для ky, но мой код работает только для L = 2n * pi для любого целого числа n. Таким образом, ошибка намного больше для других L, таких как L = 3pi или 5 pi.

Я проверил код, запустив fft и ifft снова, не применяя f = g / (kx² + ky²). Это прекрасно для любого значения L, но проблема появляется снова, когда я включаю (kx² + ky²) в код.

Я знаю, что мне не хватает чего-то простого. Любое мнение приветствуется, потому что я из идей .... Спасибо

Вот полный код 

#include <fstream>
#include <iostream> 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>

using namespace std;

int main()
{
    int N = 64;
    int i, j ;
    double pi = 3.14159265359;
    double *X, *Y  ;
    X = (double*) malloc(N*sizeof(double));
    Y = (double*) malloc(N*sizeof(double));
    fftw_complex  *out1, *in2, *out2, *in1;
    fftw_plan     p1, p2;
    double L  = 4.*pi;  //only works for 2n*pi for integer n
    double dx = L/(N);

    in1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
    out1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
    in2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
    out2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );

    p1 = fftw_plan_dft_2d(N, N, in1, out1, FFTW_FORWARD,FFTW_MEASURE );
    p2 = fftw_plan_dft_2d(N, N, in2, out2, FFTW_BACKWARD,FFTW_MEASURE);

    for(i = 0; i < N; i++){
        X[i] = -(L/2) + i*dx ;
        for(j = 0; j < N; j++){
            Y[j] = -(L/2) + j*dx ;
            in1[i*N + j][0] = sin(X[i]) + sin(Y[j]) ; // row major ordering
            in1[i*N + j][1] = 0 ;
        }
    }

    fftw_execute(p1); // FFT forward

    for ( i = 0; i < N; i++){   // f = g / ( kx² + ky² )
        for( j = 0; j < N; j++){

            double fact=0;
            if(2*i<N){fact=pow((double)i*2*pi/L,2);}
            else{fact=pow((double)(N-i)*2*pi/L,2);} 

            if(2*j<N){fact+=pow((double)j*2*pi/L,2);}
            else{fact+=pow((double)(N-j)*2*pi/L,2);}

            if(fact!=0){
                in2[i*N + j][0] = out1[i*N + j][0]/fact;
                in2[i*N + j][1] = out1[i*N + j][1]/fact;
            }else{
                in2[i*N + j][0] = 0;
                in2[i*N + j][1] = 0;
            }
        }
    }

    fftw_execute(p2); //FFT backward

    // checking the results computed
    double erl1 = 0.;
    for ( i = 0; i < N; i++) {
        for( j = 0; j < N; j++){
            erl1 += fabs( in1[i*N + j][0] -  out2[i*N + j][0]/(N*N))*dx*dx; 
            //cout<< i <<" "<< j<<" "<< sin(X[i])+sin(Y[j])<<" "<<  out2[i*N+j][0]/(N*N) <<" "<< endl; // > output
        }
    }
    cout<< "error is "<< erl1 << endl ;  // L1 error

    fftw_destroy_plan(p1);
    fftw_destroy_plan(p2);
    fftw_free(out1);
    fftw_free(out2);
    fftw_free(in1); 
    fftw_free(in2);
    return 0;
}

Аналогичный вопрос был задан в уравнении Пуассона с использованием FFTW с rectanguar домен , но автор использовал fftw_r2r вместо чисто сложного к сложному fftw.

...