FFTW3高级与基础接口2D C2R变换结果不一致问题排查
FFTW3 Fortran接口多平面C2R变换异常问题
我正在维护一套采用伪谱求解器的遗留CFD代码,计划通过替换更高效的FFT算法实现升级,目前正在学习FFTW3的Fortran接口。我有DNS代码输出的多种格式3D速度数组,用于验证FFT变换的正确性。
当前需求是对一个第一维度为物理空间、后两维度为谱空间的3D数组执行2D逆傅里叶变换(C2R),得到全物理空间的3D数组,但代码运行出现问题:
- 最初用
fftw_plan_dft_c2r成功完成单平面2D变换,但遍历所有平面时结果异常;尝试保留输入数据时触发段错误 - 转而使用高级接口
fftw_plan_many_dft_c2r批量处理多平面DFT,先测试单平面场景:预期输出为mz×mx的约150常量数组,基础接口fftw_plan_dft_c2r符合预期,但高级接口输出异常数值 - 尝试调整
istride/ostride、idist/odist、inembed/onembed等参数后问题仍未解决;更新代码明确网格尺寸和谱数组定义后,原本正常的2D变换也失效了
附测试代码:
program fftw_test use, intrinsic :: iso_c_binding implicit none include 'fftw3.f03' ! Grid sizes integer, parameter :: nxh = 128, nz = 128 integer, parameter :: mx = 384, mz = 192 ! FFTW Variables type(C_PTR) :: plan2,plan3 complex(C_DOUBLE_COMPLEX), dimension(nz,nxh) :: u_phys complex(C_DOUBLE_COMPLEX), dimension(nz,nxh) :: u1_2d,u3_2d real(C_DOUBLE), dimension(mz,mx) :: u2_2d,u4_2d integer(C_INT), dimension(2) :: msize_2d = (mx,mz) integer(C_INT) :: rank, howmany, idist, odist, istride, ostride ! ------------------------------------------------------------------------------------ ! ! ------------------------------------------------------------------------------------ ! ! Create FFT plan (before initializing variable) print *,'Create FFT Plan(s)' rank = 2; howmany = 1; ! 1 2D c2r transform on an (mz) x (mx) array ! (msize_2d is flipped because of column-major arrays) idist = 0; odist = 0; ! Unused since howmany = 1 istride = 1; ostride = 1; ! Array is contiguous in memory (I think?) plan2 = fftw_plan_dft_c2r (rank, msize_2d,u1_2d,u2_2d,FFTW_ESTIMATE) plan3 = fftw_plan_many_dft_c2r(rank, msize_2d,howmany, & u3_2d,msize_2d,istride,idist, & u4_2d,msize_2d,ostride,odist, FFTW_ESTIMATE) u_phys = (0.0,0.0) u_phys(1,1) = (150.569261744544,1.136868377216160E-013) u1_2d = u_phys ! temp storage u3_2d = u_phys ! temp storage ! Execute FFTs print *,'Execute FFTs' call fftw_execute_dft_c2r(plan2,u1_2d,u2_2d) write(200,*) u2_2d call fftw_execute_dft_c2r(plan3,u3_2d,u4_2d) write(201,*) u4_2d ! Destroy FFT plans call fftw_destroy_plan(plan2) call fftw_destroy_plan(plan3) end program fftw_test
内容的提问来源于stack exchange,提问作者Ryan Kelly
相关产品推荐
相关产品推荐

