fftw3 для пуассона с граничным условием Дирихле для всех сторон расчетной области

Question

fftw3 для пуассона с граничным условием Дирихле для всех сторон расчетной области

Я пытаюсь решить уравнение Яда с граничным условием Дирихле для четырех сторон вычислительной области. Как известно, я должен использовать FFTW_RODFT00, чтобы удовлетворить условию. Тем не менее, результат не является правильным. Не могли бы вы помочь мне?

#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;

}

1

c++ fftw poisson dirichlet

Источник

user4235046 03 фев '16 в 09:30

1 ответ

Решение

Другие вопросы по тегам c++ fftw poisson dirichlet

user2999346 03 фев '16 в 22:26 2016-02-03 22:26 · Accepted Answer · 2016-02-03 22:26

Я узнал трюк, который я вам предоставил, в уравнении Пуассона, используя FFTW с прямоугольной областью, и код, который я предоставил в своем ответе на тест Confusion fftw3 - тест 2d уравнения Пуассона, который был адаптирован из кода asker @Charles_P! Пожалуйста, рассмотрите возможность добавления ссылок на эти URL в будущих вопросах!

Ответ на проверку путаницы 2d тест уравнения fftw3 - Пуассона был посвящен случаю периодических граничных условий. Итак, вот несколько модификаций для решения случая граничных условий Дирихле.

fftw_plan_r2r_2d(N1,N2, in1.data(), out1.data(), FFTW_RODFT00, FFTW_RODFT00,FFTW_EXHAUSTIVE) соответствует дискретному синусоидальному преобразованию типа I, как определено библиотекой FFTW:

$Y_k = 2 \ sum_ {J = 0} ^ {N-1} X_j. \ Грех \ слева (\ пи (J + 1) \ гидроразрыва {к + 1} {п + 1} \ справа)$

Это значение хорошо описано в https://en.wikipedia.org/wiki/Discrete_sine_transform. Если размер массива FFTW N1=4 и его значения [a,b,c,d], полный массив, включая границы, равен [0,a,b,c,d,0]. Отсюда пространственный шаг:

$дй = \ гидроразрыв {L_1} {Н_1 +-}$

И частоты f_k типа I DST являются:

$\ Xi_k = \ р (к + 1) \ frac1 {L_1}$

Инверсией DST типа I является DST типа I, за исключением коэффициента масштабирования (см. http://www.fftw.org/doc/1d-Real_002dodd-DFTs-_0028DSTs_0029.html), Вот 4.(N1+1).(N2+1),

Наконец, контрольный пример должен быть адаптирован к случаю граничных условий Дирихле. Действительно, на коробке размера L1,L2 функция $Р (х, у) = \ Sin (\ р \ гидроразрыва {х} {L_1}) + \ Sin (\ р \ гидроразрыва {у} {} L ^ 2)$ не соблюдает граничные условия Дирихле. Действительно, даже если исходный член один и тот же, решение, удовлетворяющее периодическим граничным условиям, может отличаться от решения, удовлетворяющего граничным условиям Дирихелта. Вместо этого можно проверить два исходных термина:

Исходный термин $Р (х, у) = \ Sin (\ р \ гидроразрыва {х} {L_1}). \ Sin (\ р \ гидроразрыва {у} {} L ^ 2)$ соответствует одной частоте ДСТ.
Исходный термин $f(x,y)=2x(L_1-x)+2y(L_2-y)$ напрямую выводится из решения $f(x,y)=x(L_1-x)y(L_2-y)$

Наконец, вот код, решающий двумерное уравнение Пуассона с использованием DST типа I библиотеки 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;
}

Составлено g++ main.cpp -o main -lfftw3 -Wall,

РЕДАКТИРОВАТЬ: Как рассчитать ответ на каждую частоту?

Идея на основе FFT состоит в том, чтобы представить решение в виде линейной комбинации функций:

В случае периодических граничных условий используется БПФ. Основные функции:

$f_k(x)=e^{\sqrt{-1}\frac{2\pi k}{L_1}x} \hspace{0.5cm}\textrm{if}\hspace{0.5cm} 2k<n$

$f_k(x)=e^{-\sqrt{-1}\frac{2\pi (n-k)}{L_1}x} \hspace{0.5cm}\textrm{if}\hspace{0.5cm} 2k>n$

В случае граничных условий Дирихле используется DST типа I. Базовые функции, равные 0 при x=0 а также x=L1, являются:
$f_k (x) = \ sin \ left (\ frac {\ pi (k + 1)} {L_1} x \ right)$
В случае граничных условий Неймана используется DCT типа I. Основные функции:
$f_k (x) = \ cos \ left (\ frac {\ pi k} {L_1} x \ right)$

Их производные нулевые в x=0 а также x=L1,

Давайте вычислим вторую производную компонента f_k типа I DST:

$\frac{d^2f_k}{dx^2}(x)=-\pi^2\frac{(k+1)^2}{L_1^2}sin\left(\pi \frac{k+1}{L_1}x\right)$

Следовательно, компонент k ДСТ раствора соответствует компоненту k DST исходного слагаемого, масштабированный с коэффициентом

$\frac{L_1^2}{\pi^2(k+1)^2}$