我在使用 FFTW 求解泊松方程时遇到问题。我试图求解的方程如下:
问题是我得到的解决方案不再是径向对称的,特别是沿着轴失去对称性,正如您从输出图中看到的那样:
代码如下:
double dx = 25. / ( N*1. ), dy = 25. / ( N*1. ), dkx = 2*pi/(N*dx), dky = 2*pi/(N*dy);//define steps in real and Fourier space
double Lx = N*dx, Ly = N*dy, Lkx = N*dkx, Lky = N*dky;//define box sides in real and Fourier space
int i, j, Nh = N/2 + 1;
double *inr, *in2r, *in2d;
fftw_complex *outr;
fftw_complex *outd;
fftw_plan rplan_forward;
fftw_plan dplan_backward;
fftw_plan rplan_backward;
for ( i = 0; i < N; i++ )
{
x[i] = dx*i - Lx/2.;
kx[i] = dkx*i - Lkx/2.;
y[i] = dy*i - Ly/2.;
ky[i] = dky*i - Lky/2.;
}
outd = (fftw_complex*)fftw_malloc ( sizeof ( fftw_complex ) * N * Nh );
inr = ( double * ) malloc ( sizeof ( double ) * N * N );
for ( i = 0; i < N; i++ ) //Initialize ine to rho
{
for ( j = 0; j < N; j++ )
{
inr[i*N+j] = rho[i*N+j] * pow( -1, i + j );
}
}
outr = (fftw_complex*)fftw_malloc ( sizeof ( fftw_complex ) * N * Nh );
rplan_forward = fftw_plan_dft_r2c_2d ( N, N, inr, outr, FFTW_ESTIMATE );
fftw_execute ( rplan_forward ); //fft rho
for ( i = 0; i < N; i++ ) //Evaluate d in Fourier Space according to Poisson eq
{
for ( j = 0; j < Nh; j++ )
{
outd[i*Nh+j][0] = g * outr[i*Nh+j][0] / ( kx[i]*kx[i] + ky[j]*ky[j] );
outd[i*Nh+j][1] = g * outr[i*Nh+j][1] / ( kx[i]*kx[i] + ky[j]*ky[j] );
}
}
outd[Nh*N/2+N/2][0] = 0.; //Avoid singularity at k = 0
outd[Nh*N/2+N/2][1] = 0.;
in2d = ( double * ) malloc ( sizeof ( double ) * N * N );
dplan_backward = fftw_plan_dft_c2r_2d ( N, N, outd, in2d, FFTW_ESTIMATE );
fftw_execute ( dplan_backward ); //Antitrasform d
in2r = ( double * ) malloc ( sizeof ( double ) * N * N );
rplan_backward = fftw_plan_dft_c2r_2d ( N, N, outr, in2r, FFTW_ESTIMATE );
fftw_execute ( rplan_backward ); //Antitrasform rho
for( i = 0; i < N; i++ )
{
for ( j = 0; j < N; j++ )
{
printf( " %le %le %le \n", x[i], y[j], in2d[i*N+j] * pow( -1, i + j ) / ( N*N*1. ) - in2d[0] / ( N*N*1. ) ); //print d in such a way that it is zero at the boundaries
}
}
fftw_destroy_plan ( rplan_forward );
fftw_destroy_plan ( dplan_backward );
fftw_destroy_plan ( rplan_backward );
fftw_free ( outd );
fftw_free ( inr );
fftw_free ( outr );
fftw_free ( in2d );
fftw_free ( in2r);
fftw_cleanup();
请注意,我想要一个在边界处为零的解决方案。
您能告诉我问题出在哪里吗?我怀疑这只是一个数字问题,但我不确定。
这是边界的放大图(特别是 -12.5 < y < -11.5 区域)。
最佳答案
径向对称性的丧失与域的形状有关:它是正方形,而不是圆形。此外,由于应用了离散傅立叶变换,因此研究的问题是使用周期性边界条件的方形域上的泊松方程。因此,边界附近的温度可能不会降至零。
这就是为什么不使用零频率(即源项的平均值)的原因。如果不是这种情况,则无法求解稳态热方程。事实上,由于周期性边界条件,从一侧流出的热通量会从另一侧进入。因为当域完全绝缘时,温度会上升并且永远不会达到稳定的解决方案。
要应用零狄利克雷边界解,必须使用正弦变换之一来代替傅立叶变换。 例如,请参见:
http://www.ds.ewi.tudelft.nl/vakken/in4049TU/sheets/lecture10.pdf
这里是一个答案(我的!),其中包含一段代码,通过 FFTW 的正弦变换求解二维泊松方程、零狄利克雷边界条件。 它是 C++,但将其转换为 C 很简单。
fftw3 for poisson with dirichlet boundary condition for all side of computational domain
要恢复接近径向对称的物体,源的支撑必须远离域的边缘。如果您知道域、边界条件和源项遵循径向对称性,最好的解决方案可能是切换到一维径向方程。
关于c - 在 C 中使用 FFTW 求解泊松方程时失去径向对称性,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/50622162/