0 引言

在民航飞机制造中，翼型设计优化、风洞试验验证、全机气动优化布局、燃油消耗评估等都需要大量的计算流体动力学（computational fluid dynamic，简称CFD）计算，CFD已成为飞机优化设计和复杂流动机制分析的重要工具^[1-6]。大量的CFD应用可以降低开发成本，并为飞机研究提供技术支持，如在前期概念设计阶段，可利用CFD模拟不同机型布局方案的气动性能；在详细设计阶段，可利用CFD优化机翼、发动机整流罩等部件的外形以提高性能和减少风阻；在验证阶段，CFD可与风洞试验相互验证以减少实物试验规模。随着CFD模拟中的流动问题规模的不断扩大，复杂性不断增加，需要的求解精度也越来越高，模拟所需的模型网格大小也在增加，其网格规模可能达到上亿量级^[7-8]。巨大的数据处理需求对于执行CFD计算任务提出了挑战，提升计算效率成为亟待解决的问题。

对于CFD应用所需要的大规模计算资源，可使用多CPU并行化来加速计算。然而，低并行度和功耗效率及较为昂贵的价格可能限制CPU集群的并行性能^[9]。与CPU相比，图形处理单元（graphic processing unit，简称GPU）在数据并行性方面具有强大的浮点运算能力和高内存带宽，同时其功耗和成本也非常低，并被广泛应用于高性能异构计算（HPC）中^[10]。通过在这些异构系统上执行CFD代码，可以获得更快、更好的数值解^[11]。

早期研究为了挖掘异构系统的全部潜力并最大化应用程序性能，会进行CPU-GPU协作并行^[12].相对于GPU并行^[13]，这样可更好地配备强大多核的超级计算机。随着GPU设备的发展与进步，GPU的算力大大提升，其显存也从十年前的3 GB以下提升至现如今的80 GB^[14]，CPU-GPU协同并行计算收益下降。因此本文采用MPI-CUDA混合异构方案，CPU作为主机（Host）管理GPU的计算和通信，GPU作为设备（Device）进行高效计算。

本文将讨论CFD应用中的GPU并行化。其中内存访问和并行优化是性能优化中较为重要的关键步骤。在多GPU计算设备中，GPU通信开销较大，阻塞通信模式实现简单但效率较低，因此发展出了用于GPU并行化的非阻塞通信模式^[15]。同时为了平衡各个GPU的工作负载，可以使用领域分解方法，常用方法有一维域分解^[16]、二维域分解^[17]和三维域分解^[18]。为了减少计算量和简化算法逻辑，本文采用一维域分解方法进行计算。

本文提出了一种基于MPI和CUDA的混合并行算法，实现了CFD应用的多GPU分布式计算。在测试案例中，基于以分割模块为中心的有限体积法对可压缩的Navier-Stokes方程进行离散化，分别采用AUSM + UP迎风方案^[19]三步Runge-Kutta方法^[20]进行空间离散和时间离散，而湍流解则通过剪切应力输运（shear stress transport，简称SST）双方程模型求解^[21]。在GPU分布式并行计算具体实现中，基于异构计算和任务分离的思想，CPU仅管理GPU的执行和通信，而GPU负责数据处理。使用并行执行和内存访问优化来优化基于GPU的CFD代码。同时采用了一种非阻塞通信方法，通过创建两个CUDA流来完全重叠GPU计算，进行跨进程通信和数据传输。使用一维域分解方法来平衡GPU之间的工作以达到GPU负载均衡。此外本文还利用了分布式计算，将单节点的并行程序扩展到两个多GPU节点。所采用的MPI并行化标准可基于点对点通信和集体通信机制实现跨节点的并行协作与通信，以达到多节点集群环境下的大规模高性能计算。

本文评估了该分布式混合并行算法在超音速平板上的可压缩湍流流动计算中的性能。讨论了对于单个GPU相比于CPU获得的加速比，多GPU的强可扩展性，以及采用阻塞通信模式和非阻塞通信模式对并行效率的损失情况。当模型规模处于较大值时，单节点4GPU的加速比稳定在600倍左右，双节点8GPU加速稳定在900倍左右，证明了该并行算法具有较强加速比，良好的可扩展性以及广泛的可移植性，可无缝地在各类多核、多GPU以及多节点系统上扩展运行。这为CFD应用的类似移植和优化工作提供了普遍经验。

1 控制方程和数值模拟

在可压缩湍流的模拟中，可压缩热传导气体的三维Navier-Stockers方程的积分形式可表示为公式：

式中：$\overrightarrow{W}$是保守变量的向量，包含了流体密度、3个速度分量、能量等信息；${\overrightarrow{F}}_c$是对流通量的向量，反映了流体中的压力和速度造成的质量、动量和能量的对流传输；${\overrightarrow{F}}_v$是黏性通量向量，反映了黏性造成的动量和能量的扩散传输；t表示时间；d表示微分元；Ω表示问题的空间域；S表示Ω域内量的任何源。其中：

式中：ρ是流体密度；μ、υ、ω是三个速度分量；Ε是总能量。

式中：n_x、n_y、n_z是单位法向量；P是静压；H是总焓；V是法向量方向的速度。x、y、z分别代表空间直角坐标系中的某一点。

式中：τ_xx、τ_xy、τ_xz、τ_yx、τ_yy、τ_yz、τ_zx、τ_zy、τ_zz表示粘性应力张量的分量；Θ_x、Θ_y、Θ_z表示黏性应力和热通量的功在各个坐标分量上的功。其中：

式中：μ、υ、ω分别代表速度矢量的x方向、y方向、z方向的分量；K代表的是湍流动能；T代表的是总能量。

K-ωSST模型是一种在1994年提出的湍流模拟模型，以K-ω模型为基础，在壁面附近采用K-ω模型，在远离壁面区域采用转换后的K-ω模型，得到其积分形式暨湍流参数的控制方程整体表达式如下：

式中：${\overrightarrow{W}}_T$是湍流参数向量，${\overrightarrow{F}}_{c，T}$是湍流参数的对流通量向量，${\overrightarrow{F}}_{v，T}$是湍流参数的粘性通量向量，${\overrightarrow{Q}}_T$是湍流模型的源项向量，表达式分别如下：

公式（6）和公式（7）共同构成了SST湍流模型的方程组。从这两个公式出发，可以推导出湍流模型离散格式，即公式（8）、（9）、（10）。公式（1）是湍流参数的控制方程整体表达式，使用有限元法对公式（1）进行结构化网格空间离散。采用逆风AUSM + UP格式计算对流通量，该格式在所有速度下都具有较高的分辨率和计算效率。黏性通量由中心格式求解，湍流用K-ωSST双方程模型建模。在解空间上对时间和空间进行了单独离散化，因此可以分别处理时间和空间的积分。公式（1）可写为如下形式：

式中：Ω是控制体体积，${\overrightarrow{W}}_{q，r，s}$和${\overrightarrow{R}}_{q，r，s}$分别代表第（q，r，s）个控制体上的湍流变量向量和残差向量。本文采用三阶精度的Runge-Kutta方法将上面的公式进行时间离散化，该方法将一个时间步长内的解分解成三个子步：

将公式（9）带入公式（8）可得到如下公式：

式中：Δt是Ω的时间步长。

2 算法优化

2.1 MPI和CUDA混合并行策略

本文采用的算法是基于MPI和CUDA的混合GPU并行算法，采用单指令多线程（SIMT）的执行模式。

要在GPU上运行CFD应用程序，必须充分利用CPU和GPU。GPU的计算能力较强，且随着技术进步其显存容量已大大增加。与CPU应用程序相比，GPU数据传输效率较低，因此本文采用了NVLink这种高速、高带宽的GPU连接技术，以实现GPU之间高效的通信和数据传输^[22]。在GPU并行计算程序中，线程是内核执行的最小单位。线程在流式多处理器（SM）中作为warp并发执行，一个warp包含32个线程。因此，每个线程块的最优线程数是32的整数倍。对于当前的设备，一个线程块上的最大线程数是1 024。每个线程块上的线程数量越多，SM可以调用的线程块数量就越少。这种情况将削弱GPU利用多线程来隐藏由于线程同步问题而导致的内存获取和指令执行延迟的优势。而另一方面，线程块中的线程太少将导致线程空闲，从而导致设备利用率不足。在本研究中，使用的线程块大小为256个。线程块的数量由工作负载的规模决定，以确保每个线程都负载一个网格单元的计算。

MPI广泛用于在多GPU平台上运行并行代码的应用程序编程接口，适用于共享和分布式内存，并广泛用于大规模并行计算，在不同类型的硬件上相对容易执行并具有良好的可扩展性^[23]。在这项工作中，采用的是基于MPI和CUDA混合的并行方案，以实现多GPU上的大规模计算。图1给出了基于MPI和CUDA的多GPU集群CFD并行算法。在计算开始时，调用MPI_Init函数进入MPI环境。利用cudaMemcpyHostToDevice函数将数据从CPU传输到GPU，并在GPU上执行边界处理（Boundary_Process）、局部时间步长计算（Time_Step）、基元变量梯度计算（Grad_Primitive_Variables_Compute）、通量计算（Flux）和基元变量更新（Primitive_Variables_Update）等一系列CUDA核函数。在这些CUDA核函数中，梯度和通量的计算消耗的时间最多。对于多GPU并行计算，GPU之间需要进行数据交换，在本程序中使用Primitive_Variables_Exchange函数和Grad_Primitive_Variables_Exchange函数分别交换基元变量及其梯度。迭代结束后，调用cudaMemcpyDeviceToHost函数将数据从GPU传输到CPU进行后处理。最后，调用MPI_Finalize函数以退出MPI环境。

在多GPU并行计算中，数据传输是至关重要的。在本算法中数据的交换是由CPU控制GPU完成的。CPU和GPU之间的数据传输通过PCIe总线进行，GPU之间使用NVLink进行高效传输，过程如图2所示。NVLink在带宽上比PCIe显示出优势，使用NVLink可避免数据回传至CPU，大幅提高数据传输效率^[24]。

2.2 并行执行和内存访问优化

本文通过以下两种方法优化并行程序的性能：并行执行优化以获得最高并行利用率，内存访问优化以获得最大内存吞吐量。

CUDA程序有两种执行模式：同步模式和异步模式。同步模式是指在执行当前CUDA核函数之前，控制不会返回到主机。异步模式意味着一旦CUDA核函数启动，控制权立即返回主机。因此，主机可以启动新的CUDA核函数，同时进行数据交换。流是命令的顺序序列，可以由CUDA程序管理，以控制设备级并行性。一个流中的命令是按顺序执行的，但是来自不同命令的流可以并行执行。因此，对于大规模并行的CFD计算，需要执行一系列的CUDA核函数。这些CUDA核函数的并发执行可以通过流来实现。CUDA通过cudaStreamCreate和cudaStreamDestroy函数创建和销毁流，流同步可使用cudaStreamSynchronize函数。图3给出了CFD异步并行执行算法在GPU上的实现。梯度初始化（Grad_Init）、残差初始化（RHS_Init）、边界条件处理（Boundary_Process）、时间步长计算（Time_Step）通过创建流的方式并行执行。此外，无黏通量（Convective_Flux）的计算也可以与基元变量的梯度计算（Grad_Primitive_Variables_Compute）同时执行。

CPU和GPU之间的数据传输是通过PCIe总线实现，其传输速度明显低于GPU带宽。因此在程序中，内核迭代完全在GPU上进行，主机和GPU之间数据传输只发生在内核迭代的开始和结束，以尽可能地减少主机与设备间的数据交换。中间变量可以在数据内存中创建，并在计算完成后释放。然而，对于多GPU并行计算，GPU之间的数据交换不可避免。GPU提供了不同级别的内存结构：全局内存、纹理内存、共享内存和寄存器。这种存储方式保证了GPU能够减少全局内存与设备之间的数据传输。全局内存位于访问时延较大的显存中。因此，只有合并内存访问才能获得最大带宽。纹理存储器也是显存的一部分。与全局存储器相比，纹理存储器可以在不严格遵守合并存储器访问条件的情况下，利用缓存提高数据访问速度，获得高带宽。共享内存的带宽远高于全局内存和纹理内存。通过将频繁使用的数据存储在共享内存中，可以实现SM中线程之间的数据共享。寄存器是线程的专属存储类型，用于存储内核中声明的变量，以加速数据处理。GPU计算效率可以通过适当使用纹理内存、共享内存和寄存器以及减少对全局内存的访问次数来改进。

2.3 非阻塞通信

在多GPU上进行CFD并行计算时，内核迭代过程需要在边界上交换数据，包括基元变量的梯度。保证基元变量交换$\overrightarrow{U}[\rho ，\mu ，\nu ，\omega ，p，K，W{]}^T$量越小越好。在此过程中，存在CPU跨进程通信和GPU间跨设备数据传输。其中，GPU间数据通信优化对多GPU并行系统的性能影响很大。

传统的方法是阻塞通信模式，其算法逻辑如算法1所示，在阻塞通信模式下，GPU计算、数据传输完全分离。这种通信模式下的通信时间是纯开销，严重降低了并行系统的效率。而非阻塞通信模式可通过计算和通信的重叠，将通信时间与计算时间进行屏蔽，使用多流和调用非阻塞方法，其算法逻辑如算法2所示。

本文选取的为非阻塞通信模式。图4给出了非阻塞通信模式的算法示意图。GPU计算、GPU间数据传输的重叠是通过创建两个CUDA流来实现的。在流0上执行交换基元变量的同时，流1进行边界条件处理和时间步长，以及执行Convective_Flux函数计算。在流0上交换基元变量的梯度时，流1用于计算无黏通量。无黏通量的计算不依赖于GPU之间传递的值所以算法正确性得以保证。因此GPU计算和数据传输可以重叠。

2.4 一维域分解和GPU负载均衡

在多GPU并行计算中，需要对计算网格进行分区。考虑到GPU负载均衡问题，本文采用一维域分解方法，以近似相同数量的计算网格加载每个GPU。1维域分解如图5所示。该方法采用平分法的概念，实现简单，便于负载均衡。同时，由于其边界数据对齐，易于实现合并内存访问，从而有效地提高了边界数据通信的效率。采用非阻塞通信方式时，GPU中的数据分为上、中、下三部分。数据的上、下两部分需要与其他设备交换。上下部分的数据传输可以与中间部分的计算同时发生。这种将通信和计算重叠的非阻塞通信方式，可以显著提升GPU的通信效率，减少因等待通信而造成的计算资源浪费。通过充分利用GPU强大的并行计算能力，非阻塞通信机制使通信的开销隐藏在计算操作的后台，有效地改善了GPU作为协处理器的整体性能。

2.5 分布式策略

本文采用的MPI并行化标准可以实现在多节点集群环境下的大规模高性能计算。MPI作为一种成熟的分布式内存标准，能够充分利用分布式内存系统的并行计算资源，实现程序按节点间的并行执行。借助MPI的点对点通信和集体通信机制，可以轻松实现跨节点的并行协作与通信。因此，MPI程序具有广泛的可移植性，可无缝地在各类多核、多GPU以及多节点系统上扩展运行。

将单节点的并行程序扩展到两个计算节点，采用的分布式架构如图6所示，每个节点包含1个24核Intel® Xeon® Silver 4214 CPU @2.20 GHz（内存128 G）和4个Tesla A100 GPU。为了充分利用多个服务器的并行计算能力，本文采用了一种混合的MPI+CUDA分布式并行化方案。通过此算法在两台服务器之间进行资源分配，每台服务器配备4块Nvidia Tesla A100 GPU，节点内GPU通过NVLink互联。两台服务器之间通过MPI使用网络进行通信和同步。

MPI runtime系统负责将程序进程分配到不同计算节点，进行进程间的信息同步与交换。具体来说，CFD代码采用有限差分方法在结构化网格上求解Navier-Stockers方程。整个计算网格被均匀分区成若干子域，并均匀分配给8块GPU。引入“鬼区”^[25]来存储邻域的边界信息，即在每个计算单元的边界上创建一部分额外区域，用于存储相邻单元所需的边界信息。鬼区的数据通过定期通信操作与相邻处理单元进行同步更新，以保持数据的一致性。这样，在计算过程中计算单元可以直接访问鬼区的数据，而无需进行频繁的通信操作。在每个时间步中，各GPU并发计算自己负责的子域，根据有限差分模板计算通量并更新流场变量。GPU之间使用CUDA stream和CUDA-aware MPI进行点对点的本地边界数据交换。在进入下一个时间步之前，两个服务器节点之间通过MPI进行全局同步。

3 硬件环境和基准测试

本文采用了一个典型的可压缩湍流基准问题——超声速平板流的CFD应用，来测试所提出GPU混合并行算法的效率优化效果。超音速流过平板是CFD应用中压缩湍流流动的众所周知的基准问题，被广泛用于验证CFD代码。

其采用的具体参数为：自由流Ma数为4.5，马赫数6，基于单位长度的雷诺数为6.43×10^-6，静态温度为61.1 K，攻角为0°。在静止平板表面施加无滑移边界条件，并假定其绝热。超音速平板边界层问题在不同网格规模的模型上解决，即模型1（60万）、模型2（120万）、模型3（240万）、模型4（480万）、模型5（960万）、模型6（1 920万）、模型7（3 840万）、模型8（7 680万）。

Tesla A100主要性能参数如下：采用安培架构^[26]，CUDA计算能力为7.0，显存80 GB，SM数量108，有6 912个CUDA core，单精度算力19.5 TFlop/s，双精度算力1 935 TFlop/s，内存带宽1 935 GB/s。对于这类设备，单精度浮点运算远远超过双精度浮点运算。因此，采用单精度数据进行GPU并行化。本文使用CUDA 11.0 版本、 MPICH2 1.4.1用于MPI通信，每节点配置1个2.20 GHz 24核Intel® Xeon® Silver 4214 CPU，4个Telsa A100 GPU。

加速比（speed up，简称SP）和并行效率（parallel efficiency，简称PE）是衡量混合并行算法性能的重要参数，其定义如下^[27]：

式中：t_CPU为CPU串行的一个迭代步骤的运行时间；t_GPU为一个GPU一个迭代步骤的运行时间；W₁和W₂为问题大小；n₁和n₂为GPU数量；t₁和t₂为计算时间。下标1表示单个GPU，下标2表示多个GPU计算集群。对于强扩展率指的是W₁=W₂时（即问题规模保持不变），增加节点数，并行效率的变化情况。

4 性能分析

在不同模型规模上分别执行单CPU和多GPU并行算法，其运算耗时（实现一个迭代步骤所需要的时间，取100次平均）如表1所示。

图7展示了单GPU相比CPU串行代码的加速效果，横轴为模型规模，纵轴为运行时间，次纵轴为SP。可以看出，单个GPU的加速随着模型规模大小的增加而增大。这是因为随着计算规模的增加，计算强度（执行计算的操作数量/访问内存的数据量）增加，能够更充分地利用GPU的并行计算能力。同时，由于SM和CUDA内核数量的限制，GPU计算性能已被充分利用，因此SP渐渐趋于稳定。其中76.8百万网格的加速达到204.7倍，这说明，相较于CPU并行加速算法，GPU具有显著的性能提升。

图8是在单个计算节点上采用不同数量GPU并行所获得的加速，图9是其PE曲线图。可以看出，在模型规模较小时，GPU增加带来的效率收益较小，网格规模为0.6百万时，4个GPU获得的加速是单GPU加速的1.39倍。随着模型规模的增加，4个GPU的加速比也在不断增大并逐步达到稳定状态，当网格达到76.8百万单元时，4个GPU获得的加速是单个GPU加速的3.16倍，且强PE约为86.82%。这说明对于大规模模型，算法能够有效利用增加的GPU资源以获得更高的加速比，能够有效地利用多GPU并行计算资源，从而扩展到大规模网格的CFD计算中。

接下来讨论GPU间的非阻塞通信，以19.2百万网格为例，在一个节点上采用4个GPU进行研究，比较两种通信模式的性能。表2展示了使用阻塞通信和非阻塞通信的计算耗时区别。从表2可知，使用非阻塞通信方式可以提高GPU并行计算CFD程序的性能。随着GPU个数的增加，非阻塞通信相对于阻塞通信的加速效果也越明显。当GPU个数为1时，两种通信方式的运行时间相同，因为此时并没有GPU间数据通信。当GPU个数增加到2时，非阻塞通信的运行时间为63.61 ms，而阻塞通信为71.47 ms，非阻塞比阻塞加速了11.4%。当GPU达到3个时，非阻塞通信比阻塞通信加速了15.1%。最终在4个GPU的情况下，非阻塞通信比阻塞通信的加速效果达到19.3%。

图10展示了纯CPU并行方法、非阻塞通信GPU并行方法和阻塞通信的GPU并行方法的强扩展率曲线。增加非阻塞通信优化后的GPU强扩展率曲线与CPU并行计算强扩展率曲线几乎重合，对于4GPU，并行效率只相对于纯CPU并行计算损失了2.25%，这证明了非阻塞通信的有效性，其可通过计算、通信和重叠来屏蔽通信时间。

当采用2个计算节点（每个节点拥有4个A100 GPU）进行并行计算时，节点间需要通过网络进行通信与数据交换。为了研究跨节点数据传输对性能的影响，测试了不同大小网格在单节点和双节点下的计算性能。双节点效率等于单节点计算时间/（双节点计算时间×2），表3给出了结果。

从表3可知，双节点计算效率优于单节点，且随着网格规模的增加，双节点的计算效率提高，这是因为计算量增加，通信开销占比降低，当模型规模达到38.4百万之后，双节点效率稳定在70%以上。但由于节点间通信影响，双节点效率低于理想的2倍加速。综上所述，网络通信是影响多节点并行效率的一个重要因素。但算法可以通过增大每节点的计算量来减轻通信对性能的影响，在大规模网格场景下仍能获得较好的并行加速比。因此该算法也展现出了良好的节点间可扩展性。

这种混合并行策略最大限度减少了通信开销，同时发挥了并行计算的能力。点对点的数据交换隐藏了通信延迟，计算和通信可以重叠进行。均匀的网格分区天然实现了负载均衡。表3结果表明了本文所提出的方法在多达8个GPU（两个计算节点，每个节点拥有4个GPU）上的效率和可扩展性。这为流体计算代码提供了优秀的可伸缩性，使之能够有效利用高性能计算集群的计算能力，轻松实现由少量节点到上百节点的并行扩展，完成更大规模的高效数值模拟。

5 结论

提出了一种在多GPU设备上实现的CFD应用的MPI和CUDA的混合并行算法，并提高了GPU利用效率。使用并行执行和内存访问优化来优化基于GPU的CFD代码。根据CFD数据格式的特点，使用多CUDA 流的非阻塞通信方法，将计算和数据通信完全重叠。采用一维域分解方法，将问题平均分配到GPU之间，以保证多GPU之间的负载均衡，实现了多GPU的分布式计算。

用平板流应用对所提算法进行了评估，并对并行性能进行了详细分析。对于单个GPU，相比于CPU获得了超过200倍的加速比。同时，与数据排列和通信相比，由于内核的执行比例增加，大规模模型可以达到很好的加速比。采用的非阻塞通信模式可减少14.39%的性能损耗。当模型规模处于较大值时，单节点4个GPU的加速比稳定在640左右，双节点8个GPU加速稳定在900左右。该方法具备较好的可扩展性，CPU作为主机只负责控制GPU，因此无需考虑更换CPU导致的负载不均衡等问题，这为CFD应用的类似移植和优化工作提供了普遍经验。

参考文献