引言
单核主频提升停滞之后，数值模拟的算力增长完全依赖并行：从多核 CPU、众核 GPU 到万节点超算集群。一个千万自由度的 CFD 或隐式有限元问题，串行求解需要数天，而合理的并行策略能把它压缩到数小时。但并行不是”免费午餐”——通信开销、负载不均、内存带宽瓶颈都可能吞噬理论加速比。本文梳理并行计算的基本定律与两大主流技术路线（MPI 分布式并行与 GPU 众核加速），给出基于 mpi4py 的二维热传导域分解实现和 CuPy/NumPy 统一的 GPU 模板代码，并讨论性能建模方法。
理论基础
加速比定律：Amdahl 与 Gustafson
设程序串行部分占比为 $s$，在 $N$ 个处理器上理想并行的加速比上限由 Amdahl 定律给出：
$$S(N) = \frac{1}{s + \dfrac{1-s}{N}}$$
即使 $s = 5%$，加速比也被钉死在 20 倍——这解释了为什么大规模并行前必须先消灭串行瓶颈（I/O、全局归约、隐式求解器中的粗网格校正）。Gustafson 定律则指出：实际应用中问题规模随算力同步放大（弱扩展），此时加速比 $S(N) = N – (N-1)s$ 接近线性。对应的两个度量是强扩展（固定规模测加速）与弱扩展（每核固定负载测效率），并行效率 $E(N) = S(N)/N$ 是两者的共同标尺。
MPI 域分解与通信开销
分布式内存并行把计算域剖分为 $N$ 个子域，每个进程持有本地数据加一层”幽灵层”（halo）。以二维五点差分为例，每步需与上下左右邻居交换一条边界数据。设子域边长为 $L_s$，则通信量与计算量之比正比于：
$$\frac{\text{通信}}{\text{计算}} \propto \frac{\text{表面积}}{\text{体积}} = O\left(\frac{1}{L_s}\right)$$
子域越大通信占比越低，这与 Amdahl 定律共同决定了 MPI 程序的扩展极限。点到点阻塞通信（Sendrecv）简单可靠，规模化后应改用非阻塞通信（Isend/Irecv）与通信-计算重叠。
GPU 编程模型与 Roofline
GPU 采用 SIMT（单指令多线程）模型：线程组织为线程块（block），块组成网格（grid），数百个流多处理器同时调度数万个线程。性能要点有三：全局内存访问合并（coalescing）、共享内存复用、以及足够的占用率掩盖访存延迟。判断 kernel 性能上限的 Roofline 模型为：
$$P = \min\left(P_{\text{peak}},; \text{BW} \times \text{AI}\right)$$
其中 AI（算术强度）为每字节访存对应的浮点运算数。差分类 stencil 的 AI 仅约 0.3–0.5 FLOP/byte，远低于 GPU 的临界值（如 A100 约 20），因此 stencil 程序永远是内存带宽受限——优化方向是减少访存（时间分块、寄存器复用）而非增加算力。
| 维度 | MPI | OpenMP | CUDA |
|---|---|---|---|
| 内存模型 | 分布式，显式通信 | 共享内存，线程级 | 设备内存 + 主机拷贝 |
| 典型规模 | 数千节点 | 单节点数十核 | 单卡数万线程 |
| 编程粒度 | 粗（进程间） | 细（循环级） | 细（数据级） |
| 通信开销 | 高（网络） | 低（缓存一致性） | 中（PCIe/NVLink） |
| 适用场景 | 超大规模场求解 | 节点内加速 | 规则密集计算 |
数值方法与实现
第一个例子：mpi4py 实现的二维热方程 Jacobi 迭代，沿 $x$ 方向条带分解、邻居交换 halo：
# heat_mpi.py — 运行方式: mpiexec -n 4 python heat_mpi.py
from mpi4py import MPI
import numpy as np
comm = MPI.COMM_WORLD
rank, size = comm.Get_rank(), comm.Get_size()
up = rank - 1 if rank > 0 else MPI.PROC_NULL
down = rank + 1 if rank < size - 1 else MPI.PROC_NULL
NX, NY, NT = 512, 512, 200
lnx = NX // size # 每个进程的本地行数
u = np.zeros((lnx + 2, NY)) # 上下各一条 halo
u_new = np.zeros_like(u)
gx0 = rank * lnx # 本地首行的全局行号
for i in range(1, lnx + 1):
if abs(gx0 + i - 1 - NX // 2) < 20:
u[i, NY//2 - 20 : NY//2 + 20] = 100.0 # 中心热源
for _ in range(NT):
# halo 交换：上边界送给上方进程，同时接收下方进程的边界
comm.Sendrecv(u[1, :], dest=up, recvbuf=u[lnx + 1, :], source=down)
comm.Sendrecv(u[lnx, :], dest=down, recvbuf=u[0, :], source=up)
u_new[1:-1, 1:-1] = 0.25 * (u[:-2, 1:-1] + u[2:, 1:-1]
+ u[1:-1, :-2] + u[1:-1, 2:])
if rank == 0: u_new[1, :] = 0.0 # 全局上边界
if rank == size - 1: u_new[-2, :] = 0.0 # 全局下边界
u, u_new = u_new, u
total = comm.reduce(np.sum(u[1:-1, :]), op=MPI.SUM, root=0)
if rank == 0:
print(f"{NT} 步后全场热量 = {total:.1f}（零温边界耗散，应单调衰减）")第二个例子：CuPy/NumPy 统一的 stencil 模板——同一份代码在 GPU（装了 CuPy 时）与 CPU 上运行，演示 GPU 移植的最小改动路径：
# heat_gpu.py — 同一代码自动在 GPU/CPU 间切换
try:
import cupy as xp
backend = "CuPy (CUDA GPU)"
except ImportError:
import numpy as xp
backend = "NumPy (CPU)"
import time
def jacobi_2d(n=1024, nsteps=300):
u = xp.zeros((n, n))
u[n//2 - 50 : n//2 + 50, n//2 - 50 : n//2 + 50] = 100.0
for _ in range(nsteps):
u[1:-1, 1:-1] = 0.25 * (u[:-2, 1:-1] + u[2:, 1:-1]
+ u[1:-1, :-2] + u[1:-1, 2:])
u[0, :] = u[-1, :] = u[:, 0] = u[:, -1] = 0.0
return u
t0 = time.perf_counter()
u = jacobi_2d()
if backend.startswith("CuPy"):
_ = xp.asnumpy(u) # 拷回主机并隐式同步
print(f"{backend}: 300 步 Jacobi (1024^2) 耗时 {time.perf_counter()-t0:.2f} s")这种”数组库抽象”路线适合快速原型；追求极限性能则需手写 CUDA kernel，把五点 stencil 的一个时间步融合为单次 kernel 启动，并用共享内存缓存 tile 边界，把全局访存次数降低约 4 倍。
工程案例分析
问题：某电机整机流固耦合传热仿真
某永磁电机整机共轭换热模型约 2×10⁶ 单元，瞬态 2000 步。在集群（每节点 2×16 核 CPU）与 GPU 服务器上测得的性能如下：
| 配置 | 自由度/核（或卡） | 单步耗时 | 相对 1 节点加速比 | 并行效率 |
|---|---|---|---|---|
| 1 节点（32 核） | 62500 | 4.2 s | 1 | 100% |
| 4 节点（128 核） | 15600 | 1.2 s | 3.5 | 88% |
| 16 节点（512 核） | 3900 | 0.42 s | 10 | 63% |
| 1×A100 GPU | — | 0.19 s | ≈22 | — |
| 4×A100（NVLink） | — | 0.07 s | ≈60 | 68%（GPU 侧） |
512 核时效率降至 63%，主因是每核仅约 4000 单元、halo 交换与全局归约占比超过 30%——正好印证表面积-体积律。GPU 单卡以约 22 倍于 CPU 节点的性能胜出，但代价是模型必须装入 40 GB 显存；4 卡场景下卡间 halo 通信经 NVLink 传输，效率仍明显优于跨节点 InfiniBand。工程结论是：百万至千万自由度规模的规则求解，单/多 GPU 是性价比最高的方案；上亿自由度则必须 MPI + GPU 混合编程。
结论与展望
并行计算的决策链条可概括为：先用 Amdahl 定律评估串行瓶颈与扩展上限，再按问题规模选择技术路线——节点内 OpenMP、跨节点 MPI、规则密集 kernel 上 GPU。性能优化永远以测量为依据：MPI 程序先剖析通信占比，GPU 程序先用 Roofline 判定是算力受限还是带宽受限。前沿方向包括：任务图运行时（PaRSEC、StarPU）自动调度异构资源、通信与计算深度重叠的流水线算法、以及用神经网络替代模型压缩粗网格修正环节的多尺度混合求解。对工程团队而言，掌握”域分解 + halo 交换 + 数组库 GPU 移植”这套最小技能组合，已足以覆盖 80% 的加速需求。
发表回复