三维周期流伪谱法分解程序的高效实现

软件X 20（2022）101244原始软件出版物一种高效的三维周期流并行谱程序王连平a，王晓波a，王晓波b，王晓波a，王晓波b，王晓波ca广东省湍流研究与应用重点实验室，南方科技大学力学与航天工程系复杂流与软物质研究中心，广东深圳518055b南方科技大学粤港澳数据驱动流体力学及工程应用联合实验室，技术支持：深圳市华企网络科技有限公司哈萨克斯坦共和国国家工程学院，阿拉木图，哈萨克斯坦ar t i cl e i nf o文章历史记录：接收日期：2022年收到修订版，2022年4月20日接受，2022年保留字：伪谱法分解区域分解Taylor-Green涡流快速傅立叶变换a b st ra ct在计算流体力学领域，谱程序的数值结果是许多流体流动问题的事实标准。它们的受欢迎程度是由最高的精度加上体面的计算性能。在这项工作中，我们介绍了基于伪谱求解器的Fortran代码对于三周期流。该代码依赖于快速傅立叶变换（FFT）和FFT-FFTW版本3的最快实现之一。使用二维区域分解和FFTW库执行代码并行化。只要有可能，我们就尽量让代码更明确，即使是在为变量和数组分配新内存而不是重用现有内存导致内存使用量增加的情况使用衰减三维泰勒-格林涡流的结果进行代码验证此外，该代码还采用了CUDA并行化技术。该代码可以使用和扩展的研究社区的各种周期性流动模拟。代码的CPU和GPU并行版本的可用性为选择硬件平台提供了额外的灵活性。©2022作者（S）。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本v0.1指向此代码版本所用代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX-D-22-00023Code Ocean compute capsule不可用法律代码许可证LGPL v3.0使用git的代码版本控制系统使用Fortran、CUDA Fortran、MPI的软件代码语言、工具和服务编译要求、操作环境依赖性Linux发行版、GFortran、OpenMPI、FFTW 3、2decomp FFT、NVIDIA HPC SDK如果可用开发人员文档/手册链接https://github.com/Kairzhan/spectral/blob/main/README.md问题支持电子邮件11930952@mail.sustech.edu.cn1. 动机和意义光谱模拟在CFD界有着悠久的应用历史。高性能计算机的发展和快速傅立叶变换（FFT）算法的使用加速了它们的使用[1]。可用计算能力的不断增加导致了基础研究中的许多新发现通讯作者：广东省湍流研究与应用重点实验室，南方科技大学力学与航天工程系，复杂流动与软物质研究中心，广 东 深圳518055。电子邮件地址：wanglp@sustech.edu.cn（王连平）.https://doi.org/10.1016/j.softx.2022.101244在高雷诺数下的湍流例如，在1969年进行的第一次湍流光谱模拟中，[2]仅使用了大小为323的离散网格，而现在研究人员通常使用许多处理器内核，每个空间方向上的空间分辨率为数千个点[3]。在[4]中使用对流和对流扩散方程的全局谱分析方法[5]对伪谱方法的精度进行了分析本工作的动机是生产准备使用，简单，但快速的代码模拟周期域湍流。另外，我们想把我们以前的2352-7110/©2022作者。 由Elsevier B.V.出版。这是一篇开放获取的文章，使用CC BY许可证（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softxKairzhan Karzhaubayev，Lian-ping Wang，Dauren Zhakebayev软件X 20（2022）1012442] × [：]ijk[：] × [：==∑- − ≤ ≤=∑=×| 2k|NNNρ2频谱代码的最新版本的FFTW库的使用[6]。选择Fortran编程语言的动机是因为Fortran仍然活跃在CFD社区中，并得到主要高性能计算（HPC）供应商的支持。C/C++语言也是CFD社区的热门选择，提供类似的性能。通常，选择C/C++语言的动机是在代码中引入一些抽象级别，目的是生成通用CFD求解器。我们的目标只是为特定类型的流开发代码，而不是试图推广到不同的流问题。另一类编程语言使用解释器，技术而不是直接编译成机器可执行程序，通常为运行时环境提供Python和Julia等编程语言成功地用于复杂的流体流动问题[7，8]。它们的性能与Fortran/C/C++代码的性能然而，要实现这样的性能，需要许多优化步骤，更多的程序员我们相信Fortran仍然是CFD代码开发的一个很好的选择，从长远来看是值得的资源投资。在这里，我们打算以简单，非常明确的方式实现频谱代码，只有一些外部依赖关系，例如MPI [9]，FFTW和2decomp FFT（这里2decomp代表2D域分解）[10]。前两个库在任何HPC系统上都很容易获得，2decomp FFT可以很容易地编译到所需的计算机系统。可选的CUDA版本的代码可以使用免费分发的NVIDIA高性能计算软件开发工具包（HPC SDK）进行编译。谱流解算器的发展吸引了at-现代硬件。这些例子代表了在开发光谱模拟代码和扩展相关应用方面的各种努力。上述伪谱模拟代码的概述表明，研究人员可以访问几个开源代码，并可能在他们的研究中采用这些代码。这些代码的核心是三维场数据的FFT的有效实现在这里，我们发布了我们使用FFTW3.3.8的伪谱代码的实现代码的简单结构应该允许其他人快速修改和扩展我们的代码以满足他们的研究需求。该代码可以在串行、分布式存储器或CUDA支持的并行计算机上使用，这在选择硬件平台方面提供了额外的灵活性。许多研究人员也可能发现使用这里采用的Fortran编程语言更方便。因此，我们的工作丰富了可用的程序库的伪谱流模拟。我们将介绍代码的主要功能，包括并行化和后处理脚本，并说明单个Fortran文件的简单性和代码大小。这对广大的研究界可能非常有用。2. 软件描述应用著名的伪谱方法[16离散傅里叶变换（DFT）的使用使得将立方域的坐标范围设置为0 2π 02π 0 2π。具有以单元为中心的节点分布的物理空间的离散表示被定义为：X=（L（ i−0. 5 ），L （ j−0. 5 ），L （ k− 0. （ 5 ））、（ 1）湍流。在这里，我们回顾了这种类型的几个例子，研究代码。最新的发展之一是spectralDNS代码[7]。spectralDNS是一个简单而强大的代码，用于直接数值模拟湍流，用Python编写spectralDNS的作者还提供了C++版本用于性能比较。第二个代码HIT3D [11]是用Fortran 90编写的频谱代码，目标是湍流流的DNS 和 LES ， 可 在 https://code.google.com/p/hit3d 上 获 得 。HIT3D 能 够 跟 踪 拉 格 朗 日 粒子 和 被 动 标 量 。 教育 系 列 代 码ProtoDNS是阿尔伯塔大学的第三个代码，它是用C++语言编写的，并使用自己版本的FFTW库，具有称为fftw++的附加功能[12]。代码源可在https://github.com/dealias/dns上获得。第四个代码PsDNS（https：//github.com/lanl/PsDNS/），来自洛斯阿拉莫斯国家实验室其中L2π是域大小，N3是网格分辨率，i，j，k 1，. . . ，N.类似地，波数网格的离散表示被定义为kn1n2n3=（（2π/L）n1，（2π/L）n2，（2π/L）n3）（2）其中，对于j，波数为（N/2 1）nj N/2一二三通过向后和向前离散傅立叶变换[19]在u=u∈（k，t）eik·x，（3）K是用Python编写的，为一般的谱Galerkin方法提供了一个框架。正如PsDNS的作者所指出的，1u=N3u（x，t）e−ik·x。（四）XPsDNS的目标是设计透明和易于使用的代码。第五个代码是并行化代码Tarang [13]，其被写为该程序适用于不可压缩、磁流体动力学和热问题。作者做了代码可伸缩性测试，空间分辨率高达40963点。另一代码集（https://github.com/openmichigan/PSNM）从大学-流场由不可压缩的Navier-Stokes方程控制，体积力被吸收到压力项中，tu+ω×u= −<$（p+1u2）+ν<$2u，（5）密歇根大学提供了许多代码的集合，包 括 Navier-Stokes 方 程 的 伪 谱 求 解 器 。 代 码 可 用 Fortran 、Matlab和C语言[14]。三维版本使用MPI并行化，二维示例使用CUDA Fortran并行化。最后，利用Fortran语言在Turbo码中实现了三阶Runge-Kutta时间推进的伪谱方法https://aqua.ulb.ac.be/home/turbo/该au-thors把这个程序应用于磁流体动力学的研究和连续性方程n·u=0，（ 6）其中ωu是涡度矢量，p是压力，ρ是密度，ν是运动粘度。通过应用连续性方程去除压力项后，Navier-Stokes方程在谱空间中的表示湍流[15]。然而，代码Turbo使用FFTW库的过时版本2，这可能会降低其tu|K2|u−kk·T（七）许多从事基础物理研究的研究人员的注意力Kairzhan Karzhaubayev，Lian-ping Wang，Dauren Zhakebayev软件X 20（2022）1012443×=====××==-ˆ||其中T是u的DFTω。 首先在物理空间中计算非线性项，然后使用前向DFT将其变换到离散傅立叶空间。对非线性项采用Adams-Bashforth 方 法 进 行 时 间 积 分 ， 对 线 性 粘 性 项 采 用 隐 式 Crank-Nicolson格式，时间精度均为二阶。通过一个简单的过程来执行去混叠，该过程包括在每个时间步长处移除u 中的离散傅立叶系数（k>N/3）。 可以采用不同的时间积分方法，如多步光谱求解器是用Fortran 2003语言编写的，并在基于Linux内核的操作系统（Ubuntu 20.04和 RHEL 8）上进行了测试。源代码可在https://github.com/Kairzhan/spectral上获得。MPI版本的代码由一个名为main.f90的文件组成CUDA并行版本的代码是main.cuf，其Fortrans-C互操作性接口定义在cufft. cuf文件中。下面 的 所 有 代 码 示 例 都 引 用 文 件 main.f90 。 代 码 编 译 是 使 用Makefile执行的，Makefile负责正确的编译参数和库链接。代码输出和数据以Tecplot ASCII格式存储在文件中。可视化在开放源代码程序Paraview和商业程序Tecplot中都是可能的。在大数组的情况下，用户可以使用二进制输出与适当编写的后处理实用程序。用于二进制输出的后处理Python脚本也作为单独的程序提供。我们限制自己保持代码源代码简单越好.这导致只有一个Fortran文件实现。模拟的所有输入参数都在源文件的开头定义。数值过程的离散化是使用区域分解技术进行在域分解技术中，空间域被划分为一组子域，每个子域由单个CPU处理。由于三维DFT计算由非分解空间方向上的一维DFT组成最近在[21]中提出了关于域分解和存档中的耦合的讨论，具有高性能的更高精度。对于区域分解，在本代码中使用2decomp FFT库[102decomp& FFT库是一个Fortran库，提供1D和2D域分解框架，以及包括FFTW在内的几个FFT库的包装器。这个库用于所谓的模板格式的2D域分解。该代码使用X铅笔方向，这是2decomp FFT的FFT部分 由于在FFT子程序中利用了Hermitian对称性，所以在谱空间中的数据的阵列索引被分配有减少的维度，即，不需要存储具有n 1 0的傅立叶系数<，另外，2decomp& FFT对变换的阵列使用Z-铅笔分布。 通过子例程调用de-comp_2d_fft_get_size（fft_start，fft_end，fft_size）计算这些数组的正确索引，其中fft_start和fft_end是定义全局索引的一维三分量数组.主要工作数组被定义为三维实数 *8和复数 *8 Fortran数组。子例程alloc_x（array，opt_global=. true.） 用于使用全局索引分配分布式数组，即原始（全局）数组的每个分布式子数组具有与全局数组中匹配的访问索引，换句话说，如果数组a是数组a global的子数组，则a（i，j，k）与a_global（i，j，k）是相同的值。子例程output_fstep（nx，ny，nz，fx，fy，fz，is- tep）用于将迭代istep时的三个三维数组保存到相应的ASCII Tecplot文件中。 为了降低磁盘使用率，只能存储不同方向的三个2D切片。为了将整个单个数组以二进制格式保存为字节 流 ， 使 用 2decompFFT 库 提 供 的 子 例 程decomp_2d_write_one（1，array，filename）单GPU版本的代码使用NVIDIA CUDA技术并行化。特别是，使用了NVIDIA HPC SDK版本22.1，其中包括CUDA Fortran编译器和CUFFT FFT库。CUDA版本的代码与CPU版本的代码结构相同。算术密集型循环的消除是使用内核循环指令功能，允许类似于OpenMP的并行化。FFT程序是使用Fortran 2003的标准化C-Fortran互操作性功能，以独立于编译器的方式直接从Fortran代码中调用的。虽然我们使用Linux工作站进行代码开发，并在基于Linux的工作站上进行测试，但代码应在其他系统上运行，前提是这些系统支持代码的线性依赖关系，如MPI、FFTW3、2decomp FFT和可选的NVIDIA HPC SDK。3. 说明性实例首先，衰减的3D3D TGV流是研究由于扰动的增长通过不同物理机制的影响而产生不同尺度的涡流的最简单的系统之一，包括作为不稳定性的附加源的涡流拉伸，不仅说明了湍流的产生，而且还说明了湍流衰减。最近，在[24]中报告了3D TGV流的非线性稳定性分析。在这个问题中，初始3D速度场以以下形式给出：u（x，t=0）=sin（x）cos（y）cos（z）（8）v（x，t=0）= −cos（x）sin（y）cos（z）（9）w（x，t=0）=0（10）周期性边界条件应用于尺寸为2π2π2π的立方区域。流动的控制参数是雷诺数Reu0L0/v，其中U0，L0，v是特征速度，长度和运动粘度。对于该特定基准，使用Re600π，其中U01和L02π。平均动能随时间的变化如图所示。 1表示网格分辨率为64 3.来自短时理论解的结果，称为DUGKS的介观方法[25]，以及我们以前的光谱代码[26]也包括在内以进行比较。由于粘性耗散，平均动能单调减小。在不同的人之间达成了很好的协议，数值方法证实了代码的正确实现。给出了不同程序的t_1~5能谱图二、同样，我们以前的代码，当前的实现和CUDA版本之间的协议非常好。两个图中的模拟结果。1和2是使用网格分辨率为64 3，Courant-Friedrichs-Lewy数为CFL 0。05.作为一个额外的验证案例，将当前代码的结果与Orszag [27]和Sharma等人的结果（案例2）[28]进行了比较，用于Re下的泰勒-格林涡流200.在这如果使用643的空间分辨率，固定时间步长为不超过0。005. 图3均方涡度的演变结果表明，与已发表的结果吻合良好。 图4显示了耗散率演变的非常好的一致性。能量和耗散各向异性测量的演变如图5所示，并再次与Orszag [27]和Sharma等人[28]的结果进行比较，再次显示出非常好的一致性。Kairzhan Karzhaubayev，Lian-ping Wang，Dauren Zhakebayev软件X 20（2022）1012444·····=×=图1.一、衰 减 泰 勒 - 格 林 涡 流 的 动能演化。图二. t=5时的能谱该代码的串行版本也测试了一年-2011年4核台式机级CPU，DDR3 1333内存，结果表明，在2 ×10 - 3、64× 10 - 3和128×10 - 3网格分辨率下，每一时间步长的挂钟时间分别为0.0037、0.0396和0.384 s2处理器网格。 时间比率为110个。598，而理论比值应为1864如果不考虑内存增加的影响。Mortensen等人。[7]报告了0.0092，0.083，0.80 s的挂钟时间，表明我们的实现与他们的C++实现相当。同时，这些结果不应被视为真正的性能基准，而只是为了给出代码的总体性能预期。图中比较了通用GPU与4核和16核桌面CPU的代码性能。6，网格分辨率为643、1283和2563。 图6、标签CPU4对应4核Intel Corei5 2320 3.0 GHz CPU使用DDR3 1333图三. Re 200时均方涡度的演变，并与[27，28]的结果进行比较。见图4。Re 200耗散率的演变，并与[27，28]的结果进行比较。CPU 16对应于16核AMD Ryzen 9 5950X CPU，使用DDR4 2666MHz，GPU对应于配备12 GB GDDR6内存的NVIDIA GeForce3060。一般来说，CUDA版本的代码优于基于CPU的工作站。然而，CUDA GPU结果仅包括图中（CPU 4、CPU 16、GPU）的数值对于网格分辨率643，6为（1，6.77，85），对于分辨率，6Kairzhan Karzhaubayev，Lian-ping Wang，Dauren Zhakebayev软件X 20（2022）1012445=图六、代 码 的 并行性能。图五. Re = 200时动能（E3/E1）和耗散（D3/D1）的 各 向 异 性 比 的 演 变 ， 并与文献[27，28]的结果进行了比较。1283和（1，1.9，106）为网格分辨率2563。随着问题规模的增加，CPU上的并行化效率降低。在N256挂钟的情况下，CPU 16的时间仅比CPU 4快1.9倍。这种行为可以通过增加内存占用来解释，并考虑到DDR4内存的带宽为21333 MB/秒，而DDR3内存的带宽为10667 MB/秒。内存占用的增加使代码受到内存的限制而不是CPU的限制。为了支持这一假设，在CPU 16平台上测量了挂钟时间在这些情况下，与支持内存受限行为假设的16核运行相比，测量的挂钟时间没有显著变化。单GPU实现的代码不受数据换位程序所涉及的MPI版本代码的影响。GDDR 6的峰值带宽为360 GB/秒，并且使用的网格分辨率上的数据适合GPU内存，这使得GPU的性能优于CPU。4. 影响结论本文所介绍的程序是在分布式并行计算机集群上设计的，用于求解周期性立方区域内的三维紊流。我们特别避免了使用许多库进行输入/输出操作或数值积分，并将这些方面的选择留给最终用户，以防必要。具有良好代码可读性的单一文件格式将允许用户在应用程序和扩展功能方面扩展代码。提供了一个后处理Python脚本，用于演示等高线图生成和流场统计计算的基本例程，适用于所有代码（串行、并行和GPU版本）的输出。该代码可用于生成比较数据的发展不可压缩的三维周期性流动的情况下。给出的脱粘性泰勒-格林涡流的计算结果作为程序的验证实例。通过特定的初始条件和附加的后处理例程，代码能够模拟衰减均匀各向同性湍流（DHIT）。在一个具有I/O功能的现代16核工作站上，对网格分辨率为3843点的DHIT问题进行了大约半天的初步研究。该代码的未来扩展包括在分布式内存系统上纳入通用图形处理单元 架 构 （ GPGPU ） 、 混 合 MPI-OpenMP 并 行 化 和 高 级 I/O 。NetCDF、HDF 5和并行HDF 5 [29，30]等库和格式我们未来的发展计划包括将我们自己的铅笔分解技术，而不是什么是由2decompFFT库提供。在我们之前的代码[26]中，这种铅笔分解技术在大规模模拟竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作致谢本工作得到了国家自然科学基金（基金号11988102，91852205，&11961131006 ） 、 广 东 省 湍 流 研 究 与 应 用 重 点 实 验 室（2019B21203001）、粤港澳数据驱动流体力学及工程应用联合实验室 （ 2020B1212030001 ） 、 深 圳 市 科 学 技 术 计 划 （ 批 准 号 ：KQTD20180 411143441009）。D. Zhakebayev得到了哈萨克斯坦共和国教育和科学部的支持，资助号为AP 09260528。计算资源由南方科技大学计算科学与工程中心提供。引用[1]Cooley JW，Tukey JW.复傅里叶级数的机器计算算法。Math Comp1965;19（90）：297-301.Kairzhan Karzhaubayev，Lian-ping Wang，Dauren Zhakebayev软件X 20（2022）1012446[2] Orszag SA.湍流模拟的数值方法。物理流体1969;12（12）：II-250。[3] 作者：Shihara T，Gotoh T，Kaneda Y.高雷诺数各向同性湍流的直接数值模拟研究。流体机械年鉴2009;41：165-80。[4] 张文，张文，等.湍流数值模拟中的拟谱方法.北京：科学出版社，2000，14（1）：100 - 101. 10，2022，p. 9-24，[5] 森古普塔湖高精度计算方法：流体流动和波动现象。剑桥大学出版社;2013.[6] Frigo M，Johnson SG. FFTW：FFT的自适应软件架构。1998年IEEE声学、语音 和 信 号 处 理 国 际 会 议 论 文 集 ， ICASSP'98 。 卷 3.第 三 章 。 IEEE; 1998 ，p.1381-4.[7] Mortensen M，Langtangen HP.高性能Python为二-湍流的直接数值模拟。Comput Phys Comm2016;203：53-65.[8] Pawar S，San O. CFD Julia：构建计算流体动力学入门课程的学习模块。流体2019;4（3）：159。[9] 论坛MPI。MPI：消息传递接口标准。版本3.1。诺克斯维尔，田纳西州：田纳 西 大学; 2015年。[10]Li N，Laizet S.一&个高度可扩展的二维分解库和fft接口。在：克雷用户组2010年会议。2010年，p.1-13号。[11]丘马科夫各向同性均匀湍流中被动标量亚网格尺度通量的先验研究。物理学修订版E2008;78（3）：036313。[12]放 大 图 片 作 者 ： Robert M. 无 填 充 的 高 效 去 锯 齿 卷 积 。 SIAM J SciComput2011;33（1）：386-406.[13]Verma MK， Chatterjee A， Reddy KS， Yadav RK ， Paul S， Chandra M等。湍流模拟伪谱代码Tarang的基准和缩放研究。Pramana2013;81（4）：617-29.[14]陈刚，李宁，李文。并行谱数值方法。密歇根大学;2012年。[15][10]张文辉，张文辉，张文辉.各向异性磁流体动力学湍流中的能量传递。PhysRev E2009;79（4）：046312。[16]放大图片创作者：M. 均匀各向同性湍流中涡量和被动标量场的小尺度特征。物理流体A1991;3（6）：1587-97.[17]放大图片创作者：M.均匀各向同性湍流中小尺度结构的演化。物理流体A1992;4（12）：2747-60.[18]Canuto C，Hussaini MY，Quarteroni A，Thomas Jr A，et al.流体动力学中的光谱方法。Springer Science Business Media; 2012.[19]Briggs WL，Henson V.，DFT：离散傅立叶变换的用户手册。Philadelphia：SIAM;1995.[20]彭毅，廖伟，罗立生，王立平。衰减湍流格子玻尔兹曼方法与伪谱方法的比较：低阶统计。 Comput& Fluids 2010;39（4）：568-91.[21]Sengupta TK，Sundaram P，Suman VK，Bhaumik S. DNS紧致格式的高精度并行求解算法. ACM跨并行计算（TOPC）2020;7（4）：1-32。[22]NVIDIA HPC编译器用户2022年，可在https：//docs.nvidia.com/hpc-sdk/网站。[23]泰勒GI，绿色AE。大涡产生小涡的机理。Proc R Soc Lond Ser A Math Phys Sci1937;158（895）：499[24]Sharma N，Sengupta TK，Brinkerhoff JR. 三维Navier-Stokes方程的非线性不稳定性分析：Taylor-Green涡问题。物理流体2020;32（6）：064102。[25]薄燕，王平，郭志，王立平。三维Taylor-Green涡流和槽道湍流的DUGKS模拟。计算机&流体2017;155：9-21.[26]王立平.基于二维区域分解的三维快速傅立叶变换的并行实现及可扩展性分析.北 京 ： 计 算 机 科 学 出 版 社 ， 2001 ， 17 （ 3 ） ： 117 - 118. 并 行 计 算 2013;39（1）：58-77.[27]Orszag SA.泰勒-格林涡旋的数值模拟。应用科学与工程中的计算方法第2部分。 Springer; 1974年，p. 50比64[28]Sharma N，Sengupta TK.三维泰勒-格林涡问题的涡度动力学。物理流体2019;31（3）：035106。[29]RewR ， Davis G.NetCDF ： 科 学 数 据 访 问 接 口IEEE ComputGraphAppl1990;10（4）：76-82.[30]Rew R，Davis G. HDF5参考手册。IEEE Comput Graph Appl1990;10（4）：76-82.