, (1)
北京大学学报(自然科学版) 第62卷 第3期 2026年5月
Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 62, No. 3 (May 2026)
doi: 10.13209/j.0479-8023.2026.005
国家重点研发计划(2022YFF1300803)和天目山实验室自主科研项目(TK-2024-D-004)资助
收稿日期: 2025–03–28;
修回日期:2025–08–21
摘要 对于周期边界条件 Poisson 问题, 采用传统的 SOR 求解算法使 Poisson 求解器难以快速收敛。为解决这一问题, 针对 SOR 算法的缺点, 利用多网格加速技术, 开发一种新的有限差分法多网格迭代泊松求解器Poisson-SOR-MG。此外, 提出一种基于 1D 快速傅里叶变换(FFT)和 3D 并行通信算法耦合的 3D 傅里叶谱变换算法, 摒弃繁复的转置技术。在此基础上, 开发和验证基于标准 2π 周期的 Poisson-Fourier 求解器, 将Poisson 求解器的精度提高至谱精度。然而, 基于标准 2π 周期的 Poisson-Fourier 求解器无法应用到广泛的工程问题, 因此, 基于坐标拉伸变换的理论分析, 进一步开发和验证基于一般周期的 Poisson-Fourier 求解器。
关键词 多重网格技术; 3D傅里叶谱变换; Poisson求解器; SOR迭代法; 坐标拉伸变换
不可压缩的流体力学数值仿真应用非常广泛, 从燃油喷雾[1]、河流泥沙输运[2–3]和港口淤积[4]到海底采矿[5]等领域, 都亟需可靠的工业仿真软件进行优化设计。
一般地, 不可压缩流体力学求解器[6]包含: 非定常对流扩散求解器(式(1))和 Poisson 求解器(式(2))两部分:
, (1)
, (2)
式(1)中, 
表示爱因斯坦张量形式的速度矢量; v 表示流体运动粘度;
表示广义外源力矢量, 包括但不局限于颗粒群反馈力、气液两相表面张力和密度梯度诱导浮力等。式(2)中, p 表示压力, r表示 密度。
最终的不可压缩流体求解器的结果通过压力修正对流扩散求解器的结果获得:
。 (3)
在数值求解器开发过程中, 采用多步求解方法, 将对流方程、扩散方程和压力方程解耦, 将式(2)和(3)化为等效的式(4) [6]:
(4)
。 (5)
式(4)中, f表示压力势。式(5)为 Poisson 方程的一般形式, f表示广义 Poisson 方程源项。
值得注意的是, 除流体力学领域外, Poisson 方程在地壳重力波勘测[7]和可压缩流场 PIV 成像修 正[8]等领域也有重要或潜在的应用价值。因此, 高效且精准地求解 Poisson 方程(式(5))至关重要[9–12]。
求解 Poisson 方程的一类方案是松弛方法, 例如 Gauss-Seidel 迭代法或 Gauss-Seidel 超松弛迭代法(successive over-relaxation, SOR)等。这种方法主要依赖解的初始估计, 然后通过迭代衰减残差 r=Ñ2p– f来连续地改进估计。松弛方法的优点是在程序上易于实现, 占用存储空间少, 适合广泛的椭圆问题, 并且适用于复杂的网格。经典的迭代方法在计算大规模网格时非常耗时, 收敛速度缓慢, 达到给定残差水平所需迭代次数通常随着每个维度网格点的数量增加而增加。为改善这一问题, Fedoren-ko[13]提出多重网格法(multigrid method, MG), 此方法不依赖网格的大小, 可以显著地提高迭代松弛方法的收敛速度[14–17]。多重网格方法使用空间尺寸递增的离散化层次, 在解决初始问题时(在最细的网格上), 通过减去从较粗的网格获得的校正来迭代地改进解。由于多重网格方法的特性, 数值稳定性较好, 算法易于并行。
另一类 Poisson 求解方案是直接方法, 即直接从式(5)求解离散势。具体而言, 在预先选定的目标网格上, 依据卷积定理, 将 Poisson 方程求解过程转换至傅里叶空间进行求解。在傅里叶空间的矩阵求解转化为形式相对简易的代数运算, 从而有效地避开迭代求解方法可能面临的复杂计算和迭代过程。在计算应用中, 此过程通常使用快速傅里叶变换(Fast Fourier Transform, FFT)[18]来完成, 该变换多用于解决周期性边界条件问题[19–23]。
在大规模并行计算环境下, FFT 算法能够充分利用多核处理器以及分布式计算资源, 实现更高效的计算。一直以来, FFT 都依赖开源代码 FFTW, 该技术体系是传统的基于转置技术的 3D 傅里叶变换并行算法, 如图 1 所示, 通过 3 次转置, 分 x, y和 z三个方向进行 1D Fourier 变换, 从而实现最终的 3D傅里叶谱变换, 以 
网格的小型工业问题为例, 对 FFTW 而言, 当 x和 y方向的并行块数为 16 时, 其 Poisson 方程至少需要联合全流场求解十亿量级的变量, 计算量依旧巨大。因此, 本研究摒弃开源代码 FFTW 的技术体系和转置技术, 开发自主可控的 3D 傅里叶谱变换并行技术。
图1 采用开源 FFTW[24]进行傅里叶卷积运算的 3D 并行算法示意图
Fig. 1 Schematic diagram of 3D parallel algorithm using open source FFTW[24] for Fourier convolution operation
尤其值得关注的是, 工业场景中普遍存在的非标准周期性问题难以直接套用现有 FFT 求解框架, 而目前针对一般周期性的 Poisson 求解器的研究较少[25–29], 严重地制约了 Poisson 求解器在复杂工程场景中的实用性。
针对以上挑战, 本研究围绕一般周期性 Poisson方程的高效求解, 开发基于 1D FFT 与 3D 并行通信算法耦合的新型谱变换算法, 突破传统 FFTW 技术体系的转置限制。通过引入坐标拉伸变换理论, 将标准 2π 周期的 Poisson-Fourier 求解器拓展到一般周期的 Poisson-Fourier 求解器。结合多重网格加速技术, 开发适用于一般周期性的并行 Poisson-SOR-MG 求解器。本文核心技术创新为 3D FFT 通信的无转置并行实现; 一般周期问题的傅里叶谱方法扩展与周期适配的多重网格加速技术为经典方法的自主化工程应用, 旨在不依赖外部源代码或库。
本研究的控制方程为 Poisson 方程, 该方程基于笛卡尔直角坐标系的形式如下:
, (6)
式(6)中, p为广义 Poisson 函数解, f为广义 Poisson 源项。
傅里叶谱方法通过将微分方程转换至频域空间, 利用快速傅里叶变换(FFT)实现高效求解, 其核心优势在于将复杂的空间导数运算转化为代数操作, 从而规避传统迭代法的收敛性问题。
针对标准周期域
的 Poisson 方程(式(6)), 对物理量
进行 3D 离散傅里叶变换(Discrete Fourier Transform, DFT), 投影至波数空间 k1, k2和 k3:
(7)
是 p 的傅里叶频谱投影(本文中, 所有带上标^的物理量表示其在傅里叶变换空间中的表示); I为复指数, I2=–1;x1, x2, x3Î[0, 2p]; 波数 k1, k2和 k3分别表示
沿 3 个正交方向 e1, e2和 e3的分量, k1, k2, k3Î
, l = 1, 2, 3; N1表示坐标系
内 x1方向上网格点的总数, 也表示坐标系内的
方向上波数的总数, N2 和 N3 的定义同上。式(7)为 3D DFT 的一般算法, 该算法通过连续执行 3次一维 DFT 实现。简言之, 三维 FFT 在算法结构上与三维 DFT 一致, 其区别在于采用一维FFT 代替一维 DFT 来完成各方向的变换。
通过傅里叶变换, Poisson 方程中的导数项可简化为波数形式:
(8)
然而, 标准傅里叶谱方法受限于 2p周期域, 难以直接应用于工程中常见的任意周期域(如x, y, zÎ [0, L])。为此, 引入坐标拉伸变换技术, 将一般周期域映射至标准 2p域, 进而扩展傅里叶谱方法的适用性。映射转换关系如下:
(9)
其中, 
是 Oxyz 和坐标系之间的变换拉伸因子。特别地, 当
时, 
。
需要明确的是, 本文提出的坐标拉伸变换方法仅适用于长方体区域(边界平行于笛卡尔坐标轴的六面体区域)。该限制源于傅里叶谱变换对空间周期性的严格要求——只有在轴对齐的长方体域中, 各方向周期才能独立地分解, 确保坐标拉伸后波数空间的代数运算仍然保持谱精度。
Poisson 方程求解的核心是求解离散的 Poisson偏微分方程组。对于非复杂边界条件, 有限差分法(finite difference method, FDM)是最常用的离散数值技术之一[30]。
将 Poisson 方程写成稀疏矩阵 Lp=F 形式:
, (10)
其中, 
为稀疏矩阵中心差分算子;
为网格点(i, j, k)处的函数解;
为网格点(i, j, k)处的源项。
在 Poisson 求解器开发过程中, 稀疏线性方程(式(10))的求解占据求解器计算时间的大部分。在航空航天[31]、船舶[32]和防弹爆破[33]等实际工业问题中, 稀疏线性方程求解涉及成百万上千万的网格。受限于计算和复杂度, 迭代方法是实际应用中求解大规模稀疏方程组的主流算法[34–36]。
在笛卡尔坐标系下, 采用二阶中心差分方法在均匀网格上离散, 得到 3D 数值离散相容方程:
(11)
将 Poisson 方程的二阶中心差分数值离散相容方程写为式(10)的形式, 可得如下二阶稀疏矩阵中心差分算子矩阵:
,
,
。 (12)
式(12)的二阶稀疏矩阵中心差分算子为内部网格点的通用离散形式, 针对周期边界条件, 求解器通过虚拟点(ghost point)技术实现边界邻点值的周期映射, 通过 MPI 通信同步至本地虚拟点, 确保差分格式在周期边界下的一致性。
求解大规模稀疏矩阵的方法中, 一类算法为迭代算法, 其中 SOR 迭代方法收敛快, 且具有内在的并行性, 容易编程实现。白中治[37]验证了并行矩阵块 SOR 收敛性较好且具有很强的并行计算功能。Tang 等[38]发现在具有相同的算法和代码结构以及相同的编译优化级别的背景下, SOR 在并行模式下比 Jacobi 迭代快将近两个数量级。基于此, 本研究自研并行 Poisson-SOR 求解器。
Poisson-SOR 求解器采用 SOR 迭代方法对 3D数值离散相容方程(式(11))进行求解, 可得
(13)
式(13)表示第 m次迭代中网格点(i, j, k)处的数值解更新式。
是第 m次迭代网格点(i, j, k)处的解; 
是第 m–1 次迭代网格点(i, j, k)处的解; w是松弛因子, 用于控制解的更新速度; 
分别是x,y和 z方向上的网格步长。
Poisson-SOR 求解器在每个并行分块的边界引入 ghost point, 用于存储从相邻子块接收到的边界值。在求解器并行计算时, 稀疏矩阵被划分为多个子块, 每个子块均由独立的处理器计算。图 2 为子块 2D 网络通信图。可以看出, 在迭代过程中, 相邻子块通过通信将边界数据同步到子块的 ghost point, 从而使后续计算在本处理器中独立完成。这种设计在避免频繁跨子块通信开销的同时, 还可以确保全局解的准确性, 能够显著地减少处理器之间的等待时间和数据传递开销, 为大规模并行计算提供更高的性能支持。
对于非周期性边界这种较为简单的边界, SOR迭代算法[39–40]不失为一种既高效又经济的选择, 其优势在于通过逐次逼近的方式, 能够在相对较少的迭代步数内快速收敛至可接受的精度范围。特别是在处理大规模稀疏矩阵时, 采用并行 SOR 迭代算法能够充分利用稀疏矩阵的局部特性和迭代规则, 逐步细化解的精度, 以一种较为经济且直观的方式获得满足实际需求的结果, 为相关领域的数值计算提供了一种极具性价比的解决方案[41]。
红色虚线框内为子块的 ghost point; nx, ny 和 nz 指子块在x, y 和 z 方向上的内部网格点数
图2 Poisson-SOR 求解器并行通信技术示意图
Fig. 2 Schematic diagram of parallel communication technology Poisson-SOR solver
在笛卡尔坐标系下, 对 Poisson 方程(式(6))采用有限差分法(FDM)离散后得到线性方程组:
(14)
其中, h表示细网格步长。MG 迭代方法采用构建网格层次结构
, 可以提高求解大尺度线性方程的效率[39–40]。如图 3 所示, 本文开发的Poisson-SOR-MG 迭代求解器构建两层 V 型多重网格的层次结构{h, 2h}, 其中 h和 2h用来定义算子和解的网格, h表示细网格尺寸, 2h表示粗网格尺寸。
首先考虑细网格 h上的问题(式(14))。Poisson-SOR-MG 求解器采用并行 SOR 迭代法用于预平滑处理(smoothing), 求解得到细网格数值解 ph。经过预平滑处理, 细网格解
的高频误差被大幅度衰减, 剩余较为平滑的低频误差。细网格 h的残差(主要由解的低频部分构成)如下:
(15)
将残差数值传递到粗网格 2h中, 在粗网格尺寸下对残差进行修正。在粗网格上的残差由式(16)计算得到:
(16)
其中, 
为限制算子(restriction operator), 含义为从细网格向粗网格传递数值的操作。这一步将粗网格函数简单地等同于其对应网格值在细网格中的数值。图 4 为一维问题限制算子示意图, 图 5 为二维均匀网格的限制算子示意图。
图3 多重网格法示意图
Fig. 3 Schematic diagram of multiple network method
特别地, 在细网格 h中的低频误差在粗网格 2h层上表现为高频误差。式(17)为粗网格上的残差方程。类似平滑处理, 通过并行 SOR 迭代, 将残差的高频误差数值耗散衰减, 得到数值解的修正量:
, (17)
其中, Dp2h 为粗网格 2h中的数值解的修正量。式(18)通过延拓(插值)算子
将修正量从粗网格 2h转移到细网格 h:
(18)
图 6 和 7 给出一维和二维延拓(插值)算子示意图。
在程序中, 具体的延拓(插值)运算过程分点、线和面 3 个步骤逐步完成, 具体操作如下:
图4 1D 问题的限制算子示意图
Fig. 4 Restriction operator for one-dimensional problem
图5 2D 均匀网格的限制算子示意图
Fig. 5 Restriction operator in two-dimensions on a uniform mesh
图6 1D 延拓(插值)算子示意图
Fig. 6 Prolongation (interpolation) operator for one-dimensional problem
图7 2D 均匀网格的延拓(插值)算子示意图
Fig. 7 Prolongation (interpolation) operator in two-dimensions on a uniform mesh
(19)
在细网格 h层, 将预平滑阶段数值解
与粗网格的修正解
相加, 即可得到每个网格点更加精确的解 pʹ。
与预平滑阶段类似, 在细网格 h中, 将解 pʹ进行并行的后平滑处理, 消除由延拓带来的高频误差, 即可得到更精确的 Poisson 解。
综上所述, 求解器的计算步骤如下。
步骤 1 预平滑处理: 对细网格解
上采用SOR迭代方法执行松弛迭代, 快速衰减解的高频误差分量, 使残差
呈现低频主导特性。
步骤 2 粗网格修正: 通过限制–求解–延拓操作计算修正量
, 更新解为
。
步骤 3 后平滑处理: 对修正后的解
采用SOR 迭代方法执行松弛迭代, 主要消除延拓插值过程引入的高频伪振荡。
步骤 4 返回步骤 1, 计算新的时间步。
针对一般周期性问题, 本求解器采用 1.1 节中的坐标拉伸变换方法, 将物理空间中的直角坐标系O(x, y, z)通过坐标拉伸, 变换成傅里叶空间中的直角坐标系 O(x1, y2, z3), 将物理空间中非 2π 周期的Poisson 方程转换成傅里叶空间中 2π 周期的形式。图 8 为三维非周期性坐标系的转换示意图。
使用映射技术, 将 Poisson 方程投影到傅里叶空间, 表示为
。 (20)
分别在 x, y和 z方向上采用等距网格, 使用 FFT 方法, 可将全场联立的 Poisson 方程在 x, y和 z方向上解耦, 以便求解。谱变换如下:
图8 坐标拉伸变换示意图
Fig. 8 Grid stretching transformation diagram of coordinate system
其中, 
表示傅里叶变换后函数
的实部(real part)和虚部(imaginary part); FFT(p)表示 FFT 算子, 采用经典的蝶形算法[42]进行计算; 
是坐标转换后傅里叶空间的波数; 表示傅里叶空间的广义函数。
式(21)给出在傅里叶空间中求解 Poisson 方程的一般算法:
(21)
通过 iFFT 傅里叶反变换算法, 将计算结果从傅里叶空间转化到物理空间, 得到一般周期 Poisson方程的解:
(22)
其中, iFFT 表示傅里叶逆变算子。
在多维离散数据处理中, Poisson-Fourier 分解通过频域方法, 将计算复杂性大幅度降低, 避免直接在空间域进行数值计算的高昂代价[24]。如图 9 所示, 本文开发的Poisson-Fourier 求解器采用 3D 傅里叶并行技术, 实现对三维数据的快速处理。
传统的基于转置技术的 3D 傅里叶变换并行算法(FFTW)所需计算复杂度是
, 本研究的 Poisson 方程求解的并行算法所需计算复杂度是
, 计算复杂度整体上降低 1/ npx, 其中 N 是各方向的网格点数, npx, npy和 npz 分别表示 x, y和 z方向的并行块数, 鉴于本研究 x, y 和z三个方向的并行块数相等, 故而 y和方向的并行块数以 x方向的并行块数 npx为代表来分析。当网格点数 N 为 512, npx=16 时, 本求解器算法的计算复杂度降低 10 倍。
本研究的三维傅里叶谱变换技术摒弃传统转置依赖, 通过 MPI_AllReduce 实现 x, y和 z方向的全局数据归约, 直接耦合 1D FFT 完成三维变换, 计算复杂度显著降低, 且自主可控, 不依赖外部开源库。
图 10 示意 Poisson-Fourier 求解傅里叶卷积运算的 3D 并行通信归约技术。在基于 MPI 的并行算法中, 各处理单元初始仅存储局部数据分布。为实现x方向全局数组的获取, 系统采用 MPI_AllRe-duce关键通信操作。该操作首先对 x方向的本地数据进行求和归约计算, 随后将结果分发给所有处理单元。通过此机制, 各处理单元均可获得完整的 x 方向全局数据。该通信流程随后沿 y方向进行重复操作, 最终扩展至 z方向。全局通信策略有效地实现了跨处理单元的高效数据交互与整合。经三维数据重构后, 系统为全局 3D 傅里叶谱变换等计算提供了完整的数据基础。
基于 MPI 的设计使得代码具有较好的并行扩展性, 能够在多核处理器以及集群等不同规模的并行计算平台上方便地运行。随着计算资源的增加, 通过动态地调整 MPI 相关参数, 计算任务能灵活地分配至新增计算单元。该并行化方法可以显著地提升大规模数据的处理效率, 从而加速整个计算过程, 提高对大规模数据的处理效率。
图9 采用自主可控源代码进行傅里叶卷积运算的 3D 并行算法示意图
Fig. 9 Schematic diagram of 3D parallel algorithm for Fourier convolution operation using self-developed code
图10 采用自主可控源代码进行 X方向的傅里叶卷积运算的 3D 并行通信归约技术示意图
Fig. 10 Schematic diagram of 3D parallel communication reduction technology using self-developed code for X-direction Fourier convolution operation
本研究所有数值测试均在配备 AMD RPYC 7542 32-Core Process, 2.9GHz 处理器的 linux 系统上运行。
2.1.1 准确性验证
我们选取较为典型的一般周期的 3D Poisson 方程作为测试基准, 其表达式如式(23)所示。由于该方程形式较为简单, 可通过理论分析得到其理论解(式(24))。将理论解与求解器得到的数值解进行对比, 直接验证本文求解器程序的准确性。采用等距网格划分, 其中网格步长 dx=dy=dz=L/N(L表示计算域在 x, y和 z三个方向上的长度, N表示网格数量)。令网格规模为 323, 643, 1283和 2563。
考虑如下一般周期的 3D Poisson 方程:
(23)
其中, (x, y, z)∈[0, 2], 其理论解为
(24)
使用离散误差 err 来衡量数值解与理论解之间的误差:
(25)
其中,
分别为 x, y和 z方向的网格点数; 
表示离散点
的理论解;
表示离散点的数值解。
显然, 理论解
在 x, y和 z方向均为周期函数, 其周期分别为 Tx=2, Ty=1, Tz=2。在此基础上, 采用数值方法求解该方程, 并在 z方向选取网格点进行分析。图 11 展示网格规模为 1283 时数值解与理论解的对比结果。
通过式(25)定义的误差计算式, 在[0, 2]范围内对误差进行统计, 结果如图 12 所示。通过误差收敛阶, 进一步量化数值方法收敛特性, 定义收敛阶:
其中, h1 和 h2 为相邻两次网格加密后的步长, 
和 errh2 为对应的离散误差。对于 Poisson-SOR-MG求解器, 当网格从 323 加密至 643时, 步长 h由 6.25× 10–2 变为 3.13×10–2, 误差由 1.5×10–3降至 1.1×10–3, 代入收敛阶式, 得到收敛阶 p≈0.45。随着网格进一步加密至 1283 和 2563, 其收敛阶逐步趋近理论二阶, 从 1283 增至 2563 时, p≈1.95, 表明该求解器采用二阶精度的有限差分格式, 具有良好的网格收敛特性。
图11 网格规模为 1283 时的 Poisson-SOR-MG和 Poisson-Fourier 求解器验证
Fig. 11 Verification of Iteration errors of the Poisson-SOR-MG and Poisson-Fourier solver under grid resolutions of 1283
从左到右, 网格规模依次为323, 643, 1283和2563
图12 Poisson-SOR-MG 和 Poisson-Fourier 求解器迭代误差
Fig. 12 Iteration errors of the Poisson-SOR-MG solver and Poisson-Fourier solver
相比之下, 在不同网格规模下, Poisson-Fourier求解器的误差基本上稳定在 10⁻7量级, 并未随网格加密而发生显著变化。这是由于 Poisson-Fourier 求解器基于谱方法, 其精度已达到谱精度上限, 因此误差不再随网格细化而进一步减小。这种特性表明Poisson-Fourier 求解器在周期性问题的求解中具有极高的精度, 且对网格规模的依赖性较弱。
从以上分析可以得出, 在不同的网格规模下, Poisson-SOR-MG 求解器和 Poisson-Fourier 求解器均能有效且准确地求解一般周期的 Poisson 方程。
2.1.2 收敛性对比
为评估本文提出的两种 Poisson 求解器的性能, 将其与基于 SOR 迭代方法的并行求解器(Poisson-SOR 求解器)进行对比测试。该求解器采用二阶中心差分格式离散 Poisson 方程, 并通过并行 SOR 方法求解离散方程组。
SOR 迭代方法的收敛速率与松弛因子参数 ω存在密切的关联: 当松弛因子选择合理时, SOR 算法可以显著地加速迭代收敛; 若松弛因子选取不当, 可能导致收敛速度下降或迭代发散[39]。本文采用较为保守的松弛因子 ω=0.5 来保证方法的稳定性和一定程度的收敛性, 避免因过度松弛导致的震荡或发散现象。MG 算法通过在不同分辨率的网格间有效地转移误差, 显著地加快低频误差的消除速度, 从而提升收敛效率。傅里叶谱变换则将原问题转化为频域中的代数方程, 规避空间离散形成的稀疏矩阵直接求解问题。该算法借助 FFT 算法, 显著地降低计算复杂度, 在保持指数收敛精度的同时, 显著地提升大规模计算的效率[43–44]。
图 13 对比 Poisson 求解器在相同计算条件下的迭代步数和迭代时间。计算参数设置为网格规模N3 =1283, 网格步长 dx=1.56×10–2。横坐标的并行核数对应 3 种网格负载配置: 128grid/core 使用 13个核心, 64grid/core 使用 23个核心, 32grid/core 使用 43个 核心。
在网格规模为 1283, 步长 dx=1.56×10–2的条件下, Poisson-SOR 求解器的迭代步数高达 1630 步, 而Poisson-SOR-MG 求解器仅需 39 步, Poisson-Fourier求解器则仅需单次迭代即可达到收敛。这一结果与以下理论分析情况一致: 传统的 SOR 方法受限于局部误差传播效率, 迭代步数随网格规模线性增长; 多重网格算法通过粗细网格交替迭代, 加速低频误差消除, 迭代效率提升约 40 倍; 傅里叶谱法则通过全局频域变换, 将问题转化为代数方程, 实现一步收敛。同时, Poisson-Fourier 求解器和 Poisson-SOR-MG 求解器表现出较好的收敛特性, 两者的迭代时间均维持在 10–1s 数量级, 迭代步数均维持在较低的水平。多重网格法通过粗细网格交替迭代, 有效地提高收敛速度, 傅里叶变换则在迭代步数上表现出明显的优势, 能够有效地减少通信开销。因此, Poisson-SOR-MG 求解器和 Poisson-Fourier 求解器在迭代步数和迭代时间方面表现出明显的优势。
2.1.3 并行效率
加速比和并行效率是衡量并行计算性能的关键指标, 能够反映求解器在多处理器环境下的加速效果和资源利用率[43,45]。对于衡量相同网格规模问题实例, 加速比 Speedup=T1/TP, 其中 T1表示单核运行时间, TP表示 P核并行运行时间, P表示处理器(核心)的数量。理想情况下, 加速比等于处理器数量, 即线性加速比, 表明并行算法具有良好的可扩展性。并行效率
, 其中 P 为执行并行程序时使用的处理器的数量。理想并行效率为100%, 意味着工作负载平均分配在 P个处理器上。然而, 在实际计算中, 当有多个处理器协作处理任务时, 不可避免地出现通信开销。这种开销通常通过消息传递或共享资源的方式产生, 占用部分处理器的执行时间, 从而使加速比和并行效率难以达到理想状态。本文选取式(23)对提出的 Poisson-Fourier求解器的加速比和并行效率进行测试, 并与Poisson-SOR 和 Poisson-SOR-MG 求解器进行对比, 考察不同求解器在不同网格规模情况下的并行综合效率。图 14 给出 Poisson 求解器在相同网格规模(N=1283)的加速比。
图13 Poisson 求解器迭代步数和迭代时间的对比
Fig. 13 Comparison of iteration steps and iteration time of the Poisson solver
当核数从 13增至 43时, Poisson-Fourier 求解器加速比仅提高 6.11 倍, 而 Poisson-SOR-MG 求解器提高 31.17 倍, 说明 MG 算法具有更好的并行可扩展性。相较于 Poisson-SOR 和 Poisson-SOR-MG 求解器, Poisson-Fourier 求解器并行加速比线性扩展性较弱。这是由于 Poisson 方程求解时需要全流场联立, 在执行三维 FFT 时, 需要 x, y和 z三个方向在物理空间和傅里叶空间进行频繁的数据交换, 这种全局通信模式会随着处理器数量的增加而显著地增加通信开销。此外, 同步等待机制也对求解器的并行效率产生一定的影响。尽管 Poisson-Fourier 求解器的加速比性能表现一般, 但其单核运行时间极短, 仅为 1.86s, 是 Poisson-SOR-MG 求解器运行时间的40%, 是Poisson-SOR 求解器运行时间的 20%, 这表明 Poisson-Fourier 求解器并行加速比的单核运算优势过于突出, 以至多核心协同工作时的效率提升不够明显。此外, Poisson-Fourier 求解器在不同网格规模情况下的并行时间处于较低水平。当并行核数为23和 43时, 并行时间均处于 O(10–1)量级。在相同网格负载情况下, Poisson-Fourier 求解器和 Poisson-SOR-MG 求解器的并行时间量级基本上保持一致, 比 Poisson-SOR 求解器快 10 倍。综上所述可知, 在需要快速响应的场景中, Poisson-Fourier 求解器展现出卓越的性能。
图14 Poisson 求解器在 1283网格规模下的加速比
Fig. 14 Acceleration (speedup) of different Poisson solvers on a 1283 grid
综上所述, 本文提出的 Poisson-Fourier 求解器和Poisson-SOR-MG 求解器在处理一般周期性 Pois-son 问题时展现出良好的性能优势。Poisson-Fou-rier 求解器凭借频域变换特性, 能够实现与网格规模无关的恒定高精度(误差量级稳定在 10–7), 且其单核计算效率显著优于传统迭代方法, 适用于大规模周期性问题的快速求解。然而, 其并行加速比受限于全局数据交换需求, 在强扩展场景下存在效率瓶颈。相比之下, Poisson-SOR-MG 求解器通过多重网格加速策略, 在保证收敛精度的同时, 将迭代步数降低至传统 SOR 方法的约 2.4%, 并展现出更高的并行可扩展性(最高加速比达到 31.17 倍), 适用于复杂边界条件的工程问题。两种方法在计算时间上均比传统 SOR 方法缩短 1～2 个数量级, 体现算法创新对实际工程仿真的显著推动作用。
在三维 Taylor-Green 涡验证算例中, 鉴于湍流高波数脉动的高保真解析需求, 本研究采用基于谱方法的 Poisson-Fourier 求解器。由于 SOR 迭代方法和 MG 迭代算法受限于其二阶空间离散精度, 在湍流能谱的高波数区域易引入显著的数值耗散, 难以准确地捕捉各向同性湍流的高频脉动特性。傅里叶谱方法凭借其指数收敛特性及优异的高波数分辨率, 可有效地抑制数值耗散, 从而实现对该算例的准确模拟, 故选定 Poisson-Fourier 求解器作为本算例的基准求解器。
2.2.1 准确性验证
本研究选取不可压缩、各向同性 3D Taylor-Green 涡问题, 使用 Poisson-Fourier 求解器来模拟不同分辨率下的湍流演变过程。初始条件定义如下:
(26)
其中, (x, y, z)∈[0, 2p]3; u, v 和 w 分别为速度分量。
Taylor-Green 涡是经典的湍流研究问题, 初始条件简单且具有周期性, 伴随明显的能量耗散特性, 因此, 将 Taylor-Green 涡的能量耗散率作为程序验证的基准测试案例。Taylor-Green 涡流求解步骤见附录。流体体积为 Ω 的耗散率 ε 为
(27)
其中, ε 为流体耗散率; Ek 为流体湍动能(turbulent kinetic energy), 表达式为
, (28)
其中, uʹ, vʹ和 wʹ分别为 3 个方向上的脉动值, 
表示时间平均。Re 为无量纲参数雷诺数, 表征流体流动中惯性力与粘性力的比值:
(29)
其中, r表示流体密度, U表示特征速度, L表示特征长度, m表示流体的动力粘度。Brachet 等[46]采用不连续 Galerkin 方法得到不同雷诺数流动的解, 并建立耗散率的高精度基准数据。本研究采用 Brachet等的模拟结果作为参考, 通过直接数值模拟(DNS)来验证 Poisson 求解器。
基于 Poisson-Fourier 求解器得到的 Taylor-Green涡数值验证结果见图 15。CFL(Courant-Friedrichs-Lewy 数)用于控制时间步长, 表征时间步长与网格尺度和局部流速之间的稳定性约束关系:
(30)
其中, Dt 为时间步长, hx为 x方向网格间距, Uc为特征速度。本文中 Taylor-Green 涡算例采用无量纲变量, 因此初始特征速度尺度取 Uc=1, CFL 取值为 0.5, 总模拟时间设置为 20 (tend=20 s)。在 Re=100 和 Re= 400 工况下, 不同网格规模的湍动能耗散率演化曲线与 Brachet 等[46]的高精度基准数据均吻合良好。这表明, 尽管 Poisson 方程求解模块直接影响速度场和压力场的耦合精度, 但通过严格的谱精度控制, Poisson-Fourier 求解器能够有效地降低数值离散误差对物理量演化过程的影响。因此, 本文提出的Poisson-Fourier 求解器在捕捉湍流能量级联和耗散特性方面的精度较高, 且具有较高的稳定性, 并展现出处理实际物理问题的适用性和准确性。
2.2.2 并行效率
图 16 给出 Poisson-Fourier 求解器在不同雷诺数和网格规模下的并行加速性能。实验分为 4 组: Re= 100 时, 网格规模分别为 643 和 1283; Re=400 时, 网格规模分别为 2563 和 5123。可以看出, 在低雷诺数(Re=100)条件下, 小网格(643)的加速比随处理器数量增加而缓慢地增长, 而较大网格(1283)在相同 Re下的加速比的显著改善, 并行加速比更接近线性加速比。相比之下, 高雷诺数(Re=400)条件中, 2563 网格的加速比明显衰减, 5123 网格衰减幅度明显减缓。在低 Re 条件下, 网格规模的增大能够改善并行加速性能, 表现为更接近线性加速。然而, 高 Re湍流模拟中, 由于涡旋多尺度化和非线性耦合效应, 流动复杂度指数级增长, 通信延迟成为并行效率的主要制约因素, 网格扩展对加速比提升的作用受到抑制, 并行加速效果受限, 需针对网格配置和通信策略进一步优化。
验证结果表明, Poisson-Fourier 求解器在三维Taylor-Green 涡的高保真湍流模拟中表现出优异的精度和并行效率。在准确性方面, 求解器成功地复现不同雷诺数下湍动能耗散率的演化过程, 计算结果与高精度基准数据高度吻合, 验证了其在复杂流动多尺度解析中的可靠性。在并行性能方面, 在低雷诺数(Re=100)条件下, 1283 网格的加速比可达21.17 倍(43 核), 改善了并行加速性能; 在高雷诺数(Re=400)条件下, 受限于湍流多尺度耦合与全局通信开销, 加速比降至 14.36 倍(5123 网格), 但仍然展现出一定的可扩展性。研究结果表明, 提升网格规模可以部分地缓解高雷诺数下的并行效率衰减, 但需结合通信优化策略来应对强非线性流动问题。
图15 基于 Poisson-Fourier 求解器的 Taylor-Green 验证
Fig. 15 Taylor-Green verification based on Poisson-Fourier solver
图16 不同雷诺数(Re=100 和 Re=400)下 Poisson-Fourier 求解器加速性能对比
Fig. 16 Comparison of Speedup Ratio of the Poisson-Fourier Solver at Re=100 and Re=400
本文针对不可压缩流体力学中 Poisson 方程的求解问题, 开发基于一般周期的 Poisson-SOR-MG求解器以及 Poisson-Fourier 求解器, 并进行验证。Poisson-Fourier 快速求解器基于 1D FFT 快速傅里叶变换和 3D 并行通信算法耦合的 3D 傅里叶谱变换算法, 将计算复杂度优化为。相较于传统的基于转置技术的 3D 傅里叶变换并行算法(FFTW), 计算复杂度整体上降低 1/nps。在同一网格规模、不同并行块数的情况下, 验证了基于一般周期的 Poisson-SOR-MG 求解器和 Poisson-Fourier求解器的准确性和收敛性; 在不同雷诺数(Re)和网格规模下, 与湍流标模进行湍流能谱验证。此外, 求解器采用集体通信操作 MPI_allReduce 获取 x, y和 z方向的全局数据, 实现多核并行计算, 在一定程度上提高了计算效率, 可为大规模流场模拟提供高效的压力求解方案。
随着当前 GPU 性能的快速提升及其在科学计算中的广泛应用, 未来的研究将重点开发基于 3D-GPU 架构的 Poisson-Fourier 求解器。通过优化 GPU内存管理与任务调度算法, 充分利用其并行计算能力和高带宽优势, 以便支持更大规模网格的高效求解。同时, 结合异步通信与计算重叠技术, 探索混合 MPI-CUDA 框架, 将 GPU 的局部计算与分布式通信相结合, 降低全局通信延迟, 减少通信瓶颈, 提高并行效率, 以期进一步提升三维 FFT 在气候模拟、流体动力学和工程仿真中的计算效率, 为大规模科学与工程问题提供更高效精准的解决方案。
参考文献
[1] 钱丽娟. 雾化射流场中粒子运动和传热特性的研究[D]. 浙江: 浙江大学, 2010
[2] 傅旭东, 王光谦. 低浓度固液两相流的颗粒相动理学模型. 力学学报, 2003, 35(6): 650–659
[3] 卢新华. 波流作用下浮泥及泥沙运动的大涡数值模拟[D]. 武汉: 武汉大学, 2015
[4] 李大鸣, 李孟辉, 张彤宇, 等. 两相流全沙数学模型理论及其在台山电厂港口工程研究中的应用. 中国港湾建设, 2002, 22(5): 6–11
[5] 李炜. 大洋采矿的气力提升特性. 中国有色金属学报, 1993, 3(1): 81–86
[6] Popinet S. Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries. Journal of Computational Physics, 2003, 190(2): 572–600
[7] Scheidegger S U. Gravitational waves from 3D MHD core-collapse supernova simulations with neutrino transport [D]. Basel: University of Basel, 2011
[8] Mani A, Moin P, Wang M. Computational study of op-tical distortions by separated shear layers and turbulent wakes. Journal of Fluid Mechanics, 2009, 625: 273–298
[9] Afzal A, Ansari Z, Faizabadi A R, et al. Parallelization strategies for computational fluid dynamics software: state of the art review. Archives of Computational Methods in Engineering, 2017, 24(2): 337–363
[10] 张江涛, 王正华, 车永刚. CFD显式差分程序的自动并行技术研究. 计算机工程, 2002, 28(7): 102–103
[11] Amritkar A, Deb S, Tafti D. Efficient parallel CFD-DEM simulations using OpenMP. Journal of Computa-tional Physics, 2014, 256501–519
[12] Bykov N Y, Fyodorov S A. Data parallelization algori-thms for the direct simulation Monte Carlo method for rarefied gas flows on the basis of OpenMP technology. Computational Mathematics and Mathematical Phy-sics, 2024, 63(12): 2275–2296
[13] Fedorenko R P. A relaxation method for solving elliptic difference equations. USSR Computational Mathema-tics and Mathematical Physics, 1962, 1(4): 1092–1096
[14] Botto L. A geometric multigrid Poisson solver for do-mains containing solid inclusions. Computer Physics Communications, 2013, 184(3): 1033–1044
[15] Guillet T, Teyssier R. A simple multigrid scheme for solving the Poisson equation with arbitrary domain boundaries. Journal of Computational Physics, 2011, 230(12): 4756–4771
[16] Mehrabani M T, Nobari M R H, Tryggvason G. Acce-lerating Poisson solvers in front tracking method using parallel direct methods. Computers & Fluids, 2015, 118: 101–113
[17] Othman M, Abdullah A R. An efficient multigrid Pois-son solver. International Journal of Computer Mathe-matics, 1999, 71(4): 541–553
[18] Costa P. A FFT-accelerated multi-block finite-differen-ce solver for massively parallel simulations of incom-pressible flows. Computer Physics Communications, 2022, 271: 108194
[19] Ekici K, Djeddi R, Li H, et al. Modeling periodic and non-periodic response of dynamical systems using an efficient Chebyshev-based time-spectral approach. Journal of Computational Physics, 2020, 417: 109560
[20] Sierra-Ausin J, Citro V, Giannetti F, et al. Efficient computation of time-periodic compressible flows with spectral techniques. Computer Methods in Applied Mechanics and Engineering, 2022, 393: 114736
[21] Kozień M S. Using the Fourier methods for cycle counting of bimodal stress histories with variable in time amplitudes of components. Materials, 2022, 16 (1): no. 254
[22] Chai Guo, Wang Tian Jun. Mixed generalized Hermite-Fourier spectral method for Fokker-Planck equation of periodic field. Applied Numerical Mathematics, 2018, 133: 25–40
[23] Fu H, Liu C. A buffered Fourier spectral method for non-periodic PDE. The Global Science Journal, 2012, 9(2): 460–478
[24] Liu H R, Ng S C, Chong K L, et al. An efficient phase-field method for turbulent multiphase flows. Journal of Computational Physics, 2021, 446: 110659
[25] Mckenney A, Greengard L, Mayo A. A fast Poisson solver for complex geometries. Journal of Computa-tional Physics, 1995, 118(2): 348–355
[26] Braverman E, Israeli M, Averbuch A, et al. A fast 3D Poisson solver of arbitrary order accuracy. Journal of Computational Physics, 1998, 144(1): 109–136
[27] Hockney R W. A fast direct solution of Poisson’s equation using Fourier analysis. Journal of the ACM, 1965, 12(1): 95–113
[28] Monchiet V, Bonnet G, Lauriat G. A FFT-based method to compute the permeability induced by a Stokes slip flow through a porous medium. Comptes Rendus-Mé-canique, 2009, 337(4): 192–197
[29] Gorobets A, Trias F, Soria M, et al. A scalable parallel Poisson solver for three-dimensional problems with one periodic direction. Computers and Fluids, 2009, 39 (3): 525–538
[30] Zaman M A. Numerical solution of the Poisson equa-tion using finite difference matrix operators. Electro-nics, 2022, 11(15): 2365
[31] 阎超, 屈峰, 赵雅甜, 等. 航空航天 CFD 物理模型和计算方法的述评与挑战. 空气动力学学报, 2020, 38(5): 829–857
[32] 程宣恺, 周国平, 张雨新, 等. 模型尺度下海工船舶月池对阻力性能影响的数值模拟研究. 船舶力学, 2020, 24(5): 589–598
[33] 贺锋, 梁一鸣, 李季, 等. 异型陶瓷负泊松比复合结构的防弹防爆性能数值模拟研究. 安全与环境学报, 2024(8): 2919–2928
[34] 廖臣, 祝大军, 刘盛纲. 五点差分格式求解泊松方程并行算法的研究. 电子科技大学学报, 2008(1): 81–83
[35] Brandt A. Multi-level adaptive solutions to boundary-value problems. Mathematics of Computation, 1977, 31: 333–390
[36] Miller G H. An iterative boundary potential method for the infinite domain Poisson problem with interior Dirichlet boundaries. Journal of Computational Phy-sics, 2008, 227(16): 7917–7928
[37] 白中治. 并行矩阵多分裂块松弛迭代算法. 计算数学, 1995(3): 238–252
[38] Tang T, LIU W, Mcdonough J M. Parallelization of linear iterative methods for solving the 3-D pressure Poisson equation using various programming langua-ges. Procedia Engineering, 2013, 61: 136–143
[39] 易大义. 数值分析. 第 5 版. 北京: 清华大学出版社, 2008
[40] Chassaing P. Fundamentals of fluid mechanics: for scientists and engineers. Cham: Springer, 2022
[41] Kamowitz D. SOR and MGR[ν] experiments on the crystal multicomputer. Parallel Computing, 1987, 4(2): 117–142
[42] Cooley J W, Tukey J W. An algorithm for the machine calculation of complex Fourier series. Math Comput, 1965, 19(90): 297–301
[43] 明平剑, 张文平, 雷国东, 等. 并行多重网格方法在 CFD 中的应用. 哈尔滨工程大学学报, 2008, 29 (5): 469–473
[44] Caprace D G, Gillis T, Chatelain P. FLUPS: a Fourier-based library of unbounded Poisson solvers. SIAM Journal on Scientific Computing, 2021, 43(1): C31–C60
[45] Budiardja R D, Cardall C Y. Parallel FFT-based Pois-son solver for isolated three-dimensional systems. Computer Physics Communications, 2011, 182(10): 2265–2275
[46] Brachet M E, Meiron D I, Orszag S A, et al. Small-scale structure of the Taylor-Green vortex. Journal of Fluid Mechanics, 1983, 130(1): 411–452
附录
请访问《北京大学学报(自然科学版)》官方网站(https://xbna.pku.edu.cn)查看
Development and Validation of a Poisson Solver Library for Periodic Boundary Problems
Abstract For Poisson problems with periodic boundary conditions, the use of the traditional SOR solution algorithm does not permit the rapid convergence of the Poisson solver. To address this limitation, a novel multigrid iterative Poisson solver (Poisson-SOR-MG) based on finite difference methods was developed, leveraging multigrid acceleration techniques to overcome the drawbacks of SOR. Additionally, a 3D Fourier spectral transform algorithm was designed by coupling 1D Fast Fourier Transform (FFT) with 3D parallel communication algorithms, eliminating cumbersome transpose techniques. Building on this, a standard 2π-periodic Poisson-Fourier solver was developed and validated, achieving spectral accuracy. However, the standard 2π-periodic Poisson-Fourier solver cannot be directly applied to general engineering problems. To resolve this, theoretical analysis based on coordinate stretching transformations was conducted, leading to the further development and validation of a generalized arbitrary-periodic Poisson-Fourier solver.
Key words multigrid methods; 3D Fourier spectral transformation; Poisson solvers; SOR iterative methods; coordinate mapping transformation