我正在尝试为计算域的四个边求解具有 Dirichlet 边界条件的 Poison 方程。众所周知,我应该使用 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;
最佳答案
我认识到我在 Poisson equation using FFTW with rectanguar domain 中提供给您的技巧以及我在对 Confusion testing fftw3 - poisson equation 2d test 的回答中提供的代码,改编自提问者@Charles_P 的代码!请考虑在以后的问题中添加指向这些 url 的链接!
Confusion testing fftw3 - poisson equation 2d test 的答案致力于周期性边界条件的情况。所以这里有一些修改来解决 Dirichlet 边界条件的情况。
fftw_plan_r2r_2d(N1,N2, in1.data(), out1.data(), FFTW_RODFT00, FFTW_RODFT00,FFTW_EXHAUSTIVE)
对应于类型 I 离散正弦变换 as defined by the FFTW library :
它的含义在https://en.wikipedia.org/wiki/Discrete_sine_transform 中有很好的描述。 .如果 FFTW 数组的大小为 N1=4
且其值为 [a,b,c,d],则包含边界的完整数组为 [0,a,b,c,d,0] .因此空间步长是:
并且类型 I DST 的频率 f_k
是:
I 型 DST 的逆是 I 型 DST,除了比例因子(参见 http://www.fftw.org/doc/1d-Real_002dodd-DFTs-_0028DSTs_0029.html#g_t1d-Real_002dodd-DFTs-_0028DSTs_0029 ),此处 4.(N1+1).(N2+1)
。
最后,测试用例必须适应 Dirichlet 边界条件的情况。实际上,在大小为 L1,L2
的盒子上,函数 不遵守 Diriclet 边界条件。实际上,即使源项相同,满足周期性边界条件的解也可能与满足 Dirichelt 边界条件的解不同。相反,可以测试两个源术语:
源词对应夏令时的单一频率。
源词直接来自解决方案
最后,这是一段使用 FFTW 库的 I 类 DST 求解 2D 泊松方程的代码:
#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 的思想是将解表示为函数的线性组合:
- 在周期性边界条件的情况下,使用 FFT。基本功能是:
在 Dirichlet 边界条件的情况下,使用类型 I DST。在
x=0
和x=L1
处为 0 的基函数是:在 Neumann 边界条件的情况下,使用类型 I DCT。基本功能是:
它们的导数在 x=0
和 x=L1
处为空。
让我们计算 I 型 DST 的分量 f_k
的二阶导数:
因此,解的 DST 的分量 k
对应于源项的 DST 的分量 k
,按一个因子缩放
关于c++ - fftw3 for poisson with dirichlet boundary condition for all side of computational domain,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/35173102/