← 回到首页 Skip to main content

9.6 并行化

9.6 并行化

如今，几乎所有计算机——台式机、笔记本电脑和平板电脑——的中央处理器（CPU）都由多个核心组成。此外，现代图形处理器（GPU）还提供数百乃至数千个可进行通用编程的处理单元。然而，要利用这种大幅提升的计算能力，就必须合理设计软件。虽然并行化并不直接对应某一种数值方法，但它可能影响数值格式的选择，并且必然影响其实现。

目前，流体动力学程序采用多种编程范式实现并行化。首先是消息传递接口（MPI），它在程序层面工作。异步代理库也通过并行任务之间的消息传递实现并行计算2。不过，与 MPI 相比，它的命令集要受限得多 [136]。另一方面，开放式多处理（OpenMP）、统一计算设备架构（CUDA）、开放计算语言（OpenCL）和 OpenACC [137] 等方法属于更细粒度的方案，作用于函数或循环层面。微软提出、于 2012 年末刚发布的加速大规模并行编程模型（C++ AMP）也属于这一类 [138]。例如，可以将 MPI 与其中一种方法结合，以获得最大加速比。当然，具体效果始终取决于现有硬件的能力。

Image
Image
图 9.11 对称 Joukowsky 翼型（相对厚度 10%）的二维无黏绕流。采用结构网格，$ M_{\infty} = 10^{-2} $，$ \alpha = 3^{\circ} $，中心空间离散及显式多级时间推进格式，无预处理。上图为收敛历程，下图为压力系数与精确势流解的比较。图内 Residual：残差；Lift coefficient：升力系数；Iteration：迭代次数；Euler solver：Euler 求解器；Exact solution：精确解；Log(residual)：残差的对数。
Image
Image
图 9.12 对称 Joukowsky 翼型（相对厚度 10%）的二维无黏绕流。采用结构网格，$ M_{\infty}=10^{-2} $，$ \alpha=3^{\circ} $，二阶 Roe 迎风离散及显式多级时间推进格式，无预处理。上图为收敛历程，下图为压力系数与精确势流解的比较。图内 Residual：残差；Lift coefficient：升力系数；Iteration：迭代次数；Euler solver：Euler 求解器；Exact solution：精确解；Log(residual)：残差的对数。
Image
Image
图 9.13 对称 Joukowsky 翼型（相对厚度 10%）的二维无黏绕流。采用结构网格，$ M_{\infty}=10^{-2} $，$ \alpha=3^{\circ} $，中心空间离散、显式多级时间推进格式及 Weiss–Smith 预处理。上图为收敛历程，下图为压力系数与精确势流解的比较。图内 Residual：残差；Lift coefficient：升力系数；Iteration：迭代次数；Euler solver：Euler 求解器；Exact solution：精确解；Log(residual)：残差的对数。

下面将较详细地介绍其中四种方法。为说明实际用法，我们采用一个简单程序，通过 Jacobi 迭代求解固定边界值的二维 Laplace 方程。为了演示归约操作，还计算残差的 2-范数。示例程序的源代码位于随附软件的 Parallelization 目录，采用 C 和 C++ 编写。其组织方式与编译方法的更多细节见第12章及12.8节。

2 文中涉及的名称，如适用，均为其各自所有者的注册商标。

9.6.1 MPI

MPI 已成为用 C/C++ 或 Fortran 编写的数值程序实现并行化的经典编程模型。其基本思路是：按照处理器（或 CPU 核心）的数量，将计算网格划分成相应数量的部分；每个部分分配给一个处理器，然后并行计算流场解。一个非常简单的 C 程序可以写成：

#include <mpi.h>
...
MPI_Init(&argc,&argv); // 初始化 MPI
MPI_Comm_rank(MPI_COMM_WORLD,&myProc); // 本进程的编号
MPI_Comm_size(MPI_COMM_WORLD,&nProc); // 进程总数
...
err = ReadData(myProc);
...
err = Solve(myProc);
...
err = SaveSolution(myProc);
...
MPI_Finalize(); // 结束所有 MPI 操作

一个网格部分称为分区，其流动通常依赖其他一个或多个分区的流动；同时，每个进程使用自己的独立内存，即分布式内存模型。因此，处理器之间必须交换数据。MPI 为此提供对进程间通信的显式控制。

各分区可发送消息（MPI_Send、MPI_Isend），也可接收消息（MPI_Recv）；消息内容包括流动变量、几何数据和残差等。可汇集所有其他处理器的数据（MPI_Allreduce，例如残差范数、力和力矩），也可将一个分区的数据广播或分发给其他分区（MPI_Bcast 或 MPI_Scatter）。另一项重要功能是进程同步（MPI_Barrier）。编程接口说明、可下载的库和示例见文献 [139]；文献 [140] 提供了一个很值得参考的教程。

虽然 MPI 的概念简明易懂，困难却在实现细节上。首先需要考虑，如何将给定网格合理地拆分为多个分区。要获得良好的并行效率，各处理器必须承担大致相同的工作量，这就是负载均衡。对于结构网格，计算工作量与控制体数量成正比，即与网格单元数量（单元中心格式）或网格节点数量（节点中心格式）成正比。若网格由多个网格块组成，这些块便是分区的基础。宜沿网格块尺寸最大的方向进行切分，以尽量减少相邻分区共有的面或节点，从而减少通信数据量。

对于非结构网格，工作量还取决于分区内的边数（中位对偶格式）或单元面数（单元中心格式）。当空间精度高于二阶时，计算量还与最近邻单元的邻居数量有关。加上节点和单元编号无序，非结构网格的分区可能相当复杂。幸好已有若干库可用于这一目的。在作者看来，最方便的是 Metis $ [141] $。它运行快速、可靠，并提供多种有助于优化负载均衡的选项。

各网格分区至少需要沿共有边界交换数据。根据分区方式及所采用的离散方法，有时还需要在仅沿单元边、甚至仅在单个节点处接触的分区之间通信。结构多块流动求解器即使不采用并行计算，也已经需要相应的处理逻辑（见8.9节）。然而，为适应 MPI，非结构网格求解器需要预先采用合适的程序设计，或对现有软件作大幅重构。这与其他并行化方法有很大区别：后者可以逐步加入现有程序。

分区之间究竟需要传递哪些数据？以5.2.2节的非结构中位对偶格式为例。首先，各分区独立计算的控制体，在进程间边界处是不完整的。因此，必须将共享某个边界节点的所有分区内的局部控制体相加。如前所述，这可能要求

仅共享一个节点的分区之间也进行通信。MPI 不提供搜索相关分区的函数；该功能及其所需的全部逻辑，必须由程序员实现。

此外，高阶空间离散以及 Navier–Stokes 方程的求解，需要计算流动梯度和限制器函数。在这两种情况下，都必须知道与边界节点相连的所有边的端点上的流动变量，因此需要从周围分区传入这些数据（参见随附 Unstruct2D 示例代码中周期边界之间的数据传递位置）。残差计算完成后，也必须像上述控制体一样，在边界节点处进行求和。

显式时间推进格式（见6.1节）的时间更新不需要分区之间通信。但是，隐式残差光顺（IRS，9.3.2节）和所有隐式格式则不同。它们至少需要以近似方式耦合相邻分区，否则收敛性、甚至常常连稳定性都会明显下降。对于6.2.4节介绍的常用 LU-SGS 格式，一种方法是在分区边界沿线增加光顺扫描 $ [142] $。对于6.2.5节的 Newton–Krylov 方法，由于系统矩阵现已分布在各处理器上，必须修改预条件器。最常用的是采用分区重叠的加性 Schwarz 方法 $ [143–146] $。

PETSc 库，即可移植、可扩展的科学计算工具包，可能有助于实现流动求解器的并行化 $ [147] $。PETSc 包含一系列并行线性、非线性方程求解器和时间积分器。该库使用 MPI，可由以 C/C++ 或 Fortran 编写的 MPI 应用程序调用。

9.6.2 OpenMP

OpenMP 面向 C/C++ 和 Fortran 的共享内存并行编程 $ [148] $。它将计算任务划分为多个线程，各线程共享应用程序的一部分内存，并行执行。OpenMP 支持粗粒度和细粒度并行。与 MPI 的消息传递范式不同，它允许程序员逐步将串行程序并行化。

OpenMP 提供一组简单、数量有限的编译器指令，因此只需添加几行代码，就可实现显著的并行化。随附源代码 Parallelization 目录的 openmp 示例程序，其主循环可简化表示为：

#include <omp.h>

#pragma omp parallel default(none) shared(phi,phiOld,Ni,Nj) \
private(i,j,ij)

{
    ...
    #pragma omp for schedule(static) nowait
    for (j=1; j<(Nj-1); j++)
    {
        for (i=1; i<(Ni-1); i++)
        {
            ij = i + j*Ni;
            phi[ij] = 0.25*(phiOld[ij+1]+phiOld[ij-1]+phiOld[ij+Ni]+phiOld[ij-Ni]);
        }
    }
    ...
} // 并行区域

只需两条编译器指令（#pragma omp）即可将 j 循环并行化。第一条指令定义一个并行区域，其中的代码将并行执行；在这个例子中，它还声明各线程共享的变量。第二条指令建立一个并行 for 循环。循环本身与串行代码完全相同。只有启用相应编译器选项时，这些指令才会被处理，例如 GNU 编译器的 -fopenmp；否则，程序按串行方式执行。

OpenMP 的相对简便也伴随着一些不足。例如，它不像 MPI 那样提供并行 I/O，而是由程序员自行实现。不过，如前所述，OpenMP 与 MPI 可以结合使用，以获得最佳并行效率。对于多物理场模拟程序，这种组合尤其值得考虑，近期的含颗粒流动研究即给出了实例 $ [149] $。

9.6.3 CUDA

CUDA 是面向 C/C++ 和 Fortran 的并行编程与计算平台 $ [150] $。自 2007 年公开发布以来，CUDA 使程序开发者能够直接访问 NVIDIA GPU 并行计算单元的虚拟指令集和内存，从而将 GPU 用于通用计算，即 GPGPU。GPU 采用并行吞吐架构，适于以相对较低的单线程速度执行大量并发线程；CPU 则适于高速执行单个线程。因此，若一项计算任务可分解为大量并行线程，就能获得显著加速。

CUDA 工具包可在文献 $ [150] $ 所列地址免费获取，其中包含编译器、数学库（如提供基本线性代数功能的 CUBLAS），以及用于调试和优化应用程序性能的工具。对于 CFD 程序员，基于 CUDA 的模板库 CUSP $ [151] $ 也很有价值，其中包含多种稀疏线性系统求解器，例如 GMRES。它可由文献 $ [152] $ 所列地址免费下载。文献 $ [153] $ 详细介绍了一个采用 CUDA 加速的 CFD 程序。

CUDA 平台通过对标准 C、C++ 和 Fortran 语言的扩展提供访问。以下是随附源代码 Parallelization 目录的 cuda 示例程序片段：

#include <cuda_runtime.h>
#include <device_launch_parameters.h>
// 核函数体
__global__ void Jacobi(int Ni, int Nj, float *phi, float *phiOld)
{
    int i = blockDim.x * blockIdx.x + threadIdx.x; // i、j 是
    int j = blockDim.y * blockIdx.y + threadIdx.y; // phi、phiOld 中的下标

    int ij = i + j * Ni;
    phi[ij] = 0.25f * (phiOld[ij + 1] + phiOld[ij - 1] + phiOld[ij + Ni] + phiOld[ij - Ni]);
}
...
// 调用核函数，在 GPU 上执行
Jacobi<<<blocksPerGrid,threadsPerBlock>>>(Ni,Nj,phiDev,phiDevOld);

CUDA 允许程序员定义称为核函数的函数，由 N 个不同的 CUDA 线程并行执行 N 次。核函数用 __global__ 限定符 声明，正如上例 Jacobi 函数声明前的标记。指定核函数在 GPU 上执行时采用的块数和线程数，称为执行配置；它是置于函数名与括号参数列表之间、形如 <<<...>>> 的表达式。每个执行核函数的线程都有一个唯一索引，可在核函数内部通过内置变量 threadIdx 访问。为便于处理二维、三维网格或矩阵，可使用具有 x、y、z 分量的 dim3 结构，为各空间方向分别指定块数和线程数。上例的物理空间与计算空间均采用二维网格，因此也使用了这一功能。

CUDA 的计算效率、CUBLAS 和 CUSP 等优化库，以及开发和调试工具，都很有吸引力。工具包的另一优势是可混合使用普通代码与 CUDA 代码，两者能够一起编译。不过，CUDA 依赖 NVIDIA 图形硬件，可能构成一个不足。因此，如果 GPU 来自 ATI，或需要编写可移植软件，就需要考虑 C++ AMP 或下面介绍的 OpenCL。

9.6.4 OpenCL

OpenCL 最初由 Apple 公司发起开发。它为包含 CPU 与 GPU 的异构系统提供统一的并行编程环境 $ [154] $。OpenCL 将计算系统视为由多个计算设备组成，这些设备可以是 CPU，也可以是 GPU 等“加速器”。它采用类似 C 的语言 OpenCL C 编写在计算设备上执行的程序，称为核函数。每个计算设备通常由若干处理单元组成，GPU 通常具有数百或数千个处理器；核函数在这些处理单元上并行执行。核函数编程语言以 C99 为基础，但为适应 OpenCL 设备模型作了调整，省去了函数指针、位域和变长数组，并禁止递归。C 标准库也被一组面向数学计算的自定义函数替代。

此外，OpenCL 定义了一个应用程序编程接口（API），使主机上运行的程序能够在计算设备上启动核函数并管理设备内存。设备可能与主机 CPU 共享内存，也可能不共享。API 提供设备内存缓冲区的句柄，以及在主机与计算设备之间双向传输数据的函数。OpenCL 标准定义了 C 和 C++ 接口；另有 Python、Java 等语言的第三方接口。OpenCL C 在运行时编译，从而使应用程序可移植到不同设备的实现上。不过，这也意味着 OpenCL 应用程序必须随附纯文本形式的核函数代码，对专有软件而言是个问题。

Back to top