核数增加64倍，效率却降到16.3% — 以细化网格方式测量weak scaling的陷阱

weak scaling曲线平不平，取决于预处理器的层数，而不是并行化的质量。

效率16.3%，但并行化没有缺陷#

申请超算机时需要附上scaling曲线。固定每核的单元数，把核数从1增到4、16、64，逐个测时间。理想情况下线应该是平的。我拿到的线不平。效率从1.000掉到了0.163。MPI通信量不是元凶，负载不均衡也不是。元凶藏在测量方式本身。本文沿着Hammond、Lichtner与Mills(2014)的PFLOTRAN性能论文所指出的那个位置走一遍，再用约一百行Python重现同一现象，并用一层把它拉回来。

读完剩下两个判断。倾斜的weak scaling指向代码，还是指向预处理。以及"多重网格对weak scaling有利"这句话，只在层数跟着网格一起增加时成立，层数一旦锁死就站不住。

选择细化网格那一侧，问题本身就变了

weak scaling有两种测法。固定网格间距而扩大区域，或者固定区域而细化网格。PFLOTRAN的作者选了后者。理由很实际。要建模的是一个含水层，而含水层不会随核数成比例变大。区域始终是5000 m × 2500 m。多拿到核，只是把这个区域切得更细。

可是这个选择悄悄越了线。扩大区域那一侧只增加自由度，所求解的问题不变。细化网格那一侧则是每个数据点求解的问题都不一样。离散化矩阵的特征值谱会跟着网格间距一起移动。64核的运行与1核的运行并不是同一件事分成64份，而是解了两个不同的线性系统。

在下面的模拟中亲手改变预处理器与网格细化程度。

refinement: 4 processes, 128×64 grid, 2048 cells per process

sweeps to reach 1e-8: 0 and counting

把滑块向右推，核数与网格同时各增4倍，每核单元数固定为2048。用三个按钮切换预处理器，比较sweeps计数器。在point Jacobi下，注入井与抽水井附近的短波动很快消失，但横跨整个区域的平缓污斑一直残留到最后。正是那片污斑决定了迭代次数。

strong与weak回答的是不同的问题#

先把两个指标的定义钉住。strong scaling固定问题规模而增加核数。相对加速比与相对效率这样写。

S_p = \frac{T_{\mathrm{base}}}{T_p}, \qquad E_p = \frac{T_{\mathrm{base}}\, p_{\mathrm{base}}}{T_p\, p} = \frac{S_p}{S_{p,\mathrm{ideal}}}

$T_p$ 是 $p$ 个核上的墙钟时间， $p_{\mathrm{base}}$ 是基准核数。在数千核的规模上，问题装不进一个核，所以基准总是相对的。

weak scaling固定每核的问题规模，并把两者按同一倍率 $\alpha$ 放大。

E_{n,p} = \frac{T_{n_{\mathrm{base}},\, p_{\mathrm{base}}}}{T_{n,p}}, \qquad n = \alpha\, n_{\mathrm{base}}, \quad p = \alpha\, p_{\mathrm{base}}

这里的理想值是 $E_{n,p} = 1$ ，也就是一条平线。负载不均衡另行测量： $\mathrm{LI} = (\mathrm{dof}_{\max} - \mathrm{dof}_{\min}) / \mathrm{dof}_{\min}$ 接近0表示分得均匀。我这里的区域分解完全均匀，所以 $\mathrm{LI} = 0$ 。于是不均衡从嫌疑名单中排除。

迭代次数按 $h^{-1}$ 增长的算术#

剩下的嫌疑对象是线性求解器。看一眼离散拉普拉斯算子的特征值，答案就出来了。边长为 $L$ 的区域以间距 $h$ 划分的五点差分矩阵中

\lambda_{\min} \simeq \frac{\pi^2 h^2}{L^2}, \qquad \lambda_{\max} \simeq 8, \qquad \kappa(A) = \frac{\lambda_{\max}}{\lambda_{\min}} \simeq \frac{8 L^2}{\pi^2 h^2}

$\kappa$ 是条件数，即最大与最小特征值之比。网格加密一半， $\kappa$ 变成4倍。共轭梯度法的迭代次数上界与条件数的平方根成正比。

k_{\mathrm{CG}} \lesssim \tfrac{1}{2}\sqrt{\kappa}\, \ln\frac{2}{\varepsilon} \sim \mathcal{O}(h^{-1}) = \mathcal{O}(N^{1/2})

$N$ 是未知数个数， $\varepsilon$ 是目标相对残差。在二维中， $h$ 减半则单元数变4倍，迭代次数变2倍。正如GMRES与Arnoldi过程中所见，Krylov子空间每做一次矩阵向量乘才长一维。信息穿过网格的速度就是一次迭代一个单元。

点Jacobi预处理改不了这个算术。用对角元去除只是把谱平移和缩放，并不触动两端之比。块Jacobi的情形也类似。它在块内精确求解，但跨越块边界的长波长仍然只能靠迭代传递。核数越多块切得越碎，反而更吃亏。

在同一个doublet上挂了三种预处理#

我把论文的regional doublet缩小后重现：5000 m × 2500 m的均质含水层，自西向东施加的水头差，一口注入井和一口抽水井。稳态饱和流，所以模型就是 $-\nabla\cdot(K\nabla h) = q$ 一行。固定每核2048个单元，把网格分四级加密，数PCG的迭代次数。

import numpy as np
 
def apply_laplace(x):
    """h^2-scaled 5-point operator: Dirichlet on x-faces, no-flow on y-faces."""
    ny, nx = x.shape
    xp = np.empty((ny + 2, nx + 2))
    xp[1:-1, 1:-1] = x
    xp[0, 1:-1] = x[0, :]        # no-flow (ghost = mirror)
    xp[-1, 1:-1] = x[-1, :]
    xp[1:-1, 0] = -x[:, 0]       # Dirichlet 0 (ghost = -x)
    xp[1:-1, -1] = -x[:, -1]
    return 4.0 * x - xp[1:-1, :-2] - xp[1:-1, 2:] - xp[:-2, 1:-1] - xp[2:, 1:-1]
 
def diag_of(ny, nx):
    d = np.full((ny, nx), 4.0)
    d[:, 0] += 1.0
    d[:, -1] += 1.0
    d[0, :] -= 1.0
    d[-1, :] -= 1.0
    return d
 
def make_rhs(ny, nx, head_w, head_e, well_rate):
    """regional gradient + injection/extraction doublet, h^2-scaled."""
    b = np.zeros((ny, nx))
    b[:, 0] += 2.0 * head_w
    b[:, -1] += 2.0 * head_e
    b[ny // 2, nx // 10 * 3] += well_rate
    b[ny // 2, nx // 10 * 7] -= well_rate
    return b
 
def jacobi_smooth(x, b, d, nu, omega=0.8):
    for _ in range(nu):
        x += omega * (b - apply_laplace(x)) / d
    return x
 
def restrict_avg(r):
    return 0.25 * (r[0::2, 0::2] + r[1::2, 0::2] + r[0::2, 1::2] + r[1::2, 1::2])
 
def prolong_const(e):
    return np.repeat(np.repeat(e, 2, axis=0), 2, axis=1)
 
def vcycle(x, b, level, max_level):
    ny, nx = x.shape
    d = diag_of(ny, nx)
    if level == max_level or min(ny, nx) <= 4:
        return jacobi_smooth(x, b, d, 40)
    x = jacobi_smooth(x, b, d, 2)
    rc = 4.0 * restrict_avg(b - apply_laplace(x))
    ec = vcycle(np.zeros_like(rc), rc, level + 1, max_level)
    x += prolong_const(ec)
    return jacobi_smooth(x, b, d, 2)
 
def pcg_solve(b, precond, rtol=1e-8, maxit=20000):
    x = np.zeros_like(b)
    r = b - apply_laplace(x)
    nrm0 = np.linalg.norm(r)
    z = precond(r)
    p, rz = z.copy(), float(np.sum(r * z))
    for k in range(1, maxit + 1):
        Ap = apply_laplace(p)
        alpha = rz / float(np.sum(p * Ap))
        x += alpha * p
        r -= alpha * Ap
        if np.linalg.norm(r) / nrm0 < rtol:
            return x, k
        z = precond(r)
        rz_new = float(np.sum(r * z))
        p = z + (rz_new / rz) * p
        rz = rz_new
    return x, maxit
 
def weak_scaling_table(dof_per_proc=2048):
    print(f"{'procs':>6}{'grid':>12}{'dof':>9} | {'Jacobi':>7}{'E':>7} | "
          f"{'MG-3lv':>7}{'E':>7} | {'MG-grow':>8}{'E':>7}")
    base = {}
    for k in range(4):
        p = 4 ** k
        ny, nx = 32 * 2 ** k, 64 * 2 ** k
        b = make_rhs(ny, nx, 10.0, 0.0, 40.0)
        assert ny * nx == dof_per_proc * p
        row = []
        for name, pre in (
            ("jac", lambda r: r / diag_of(*r.shape)),
            ("mg3", lambda r: vcycle(np.zeros_like(r), r, 0, 2)),
            ("mgg", lambda r: vcycle(np.zeros_like(r), r, 0, 99)),
        ):
            _, it = pcg_solve(b, pre)
            base.setdefault(name, it)
            row.append((it, base[name] / it))
        print(f"{p:>6}{f'{nx}x{ny}':>12}{ny*nx:>9} | " + " | ".join(
            f"{it:>7}{eff:>7.3f}" for it, eff in row))
 
weak_scaling_table()

 procs        grid      dof |  Jacobi      E |  MG-3lv      E |  MG-grow      E
     1       64x32     2048 |     210  1.000 |      10  1.000 |       9  1.000
     4      128x64     8192 |     407  0.516 |      16  0.625 |       9  1.000
    16     256x128    32768 |     754  0.279 |      28  0.357 |      10  0.900
    64     512x256   131072 |    1285  0.163 |      47  0.213 |      10  0.900

层数滑块和网格滑块分开推，右侧被高亮的效率曲线会换一条。

multigrid levels = 3grid = 256×128 on 16 processes

coarsest grid 64×32

层数放到 $1$ ，衰减条只剩一根。网格往右推，最粗网格变大，Jacobi效率往 $0.163$ 倾斜。层数跟着网格一起加，曲线留在 $0.900$ 附近。

点Jacobi那一列的迭代次数是210→407→754→1285。网格加密8倍的过程中变成了6.1倍。逐级之比为1.94、1.85、1.70，与上一节 $\mathcal{O}(h^{-1})$ 的估计一致。每核单元数固定，所以每次迭代的每核开销也固定。于是效率按迭代次数的反比一路掉到0.163。这不是靠调通信或缓存能夺回来的数字。

点Jacobi的 $210$ 、 $407$ 、 $754$ 、 $1285$ 相除，得到 $407/210=1.938$ ， $754/407=1.853$ ， $1285/754=1.704$ 。与表里的 $1.94$ 、 $1.85$ 、 $1.70$ 一致。三者相乘是 $6.10$ ， $1285/210=6.12$ 。效率是这个数的倒数。 $210/1285=0.1634$ ，也就是表里的 $0.163$ 。MG-3lv是 $10$ 、 $16$ 、 $28$ 、 $47$ 。 $47/10=4.7$ ，比Jacobi的 $6.1$ 小，但效率是 $10/47=0.213$ 。层数跟着增长的那一列是 $9$ 、 $9$ 、 $10$ 、 $10$ 。在 $64$ 核上 $10/9=1.11$ ，效率 $0.900$ 就是 $9/10$ 。同一份 $131072$ 个自由度上，Jacobi的 $1285$ 与增长层数的 $10$ 之比是 $128.5$ 倍。每核 $2048$ 个单元四行都一样。 $64\times32=2048$ ， $131072/64=2048$ 。 $h$ 每减半一次，单元数乘 $4$ ， $\kappa$ 乘 $4$ ，迭代上界按 $\sqrt{\kappa}$ 大约乘 $2$ 。三次减半就是加密 $8$ 倍，测到的 $6.1$ 倍是这个 $2$ 稍微缺了一点的乘积，即 $1.94\times1.85\times1.70$ 。

一层把47拉回了10#

第三列是答案。把V-cycle的层数固定为3层，迭代从10长到47。让层数随网格一起增加，使最粗的网格始终保持8 × 4个单元，则是9→10。效率0.900。实际上已经平了。

原因一句话就能写完。Krylov求解器的行为像单层网格上的平滑子，很快擦掉高频误差，把低频误差留下。层数跟到网格那么深，一个循环就能碰到长波。预处理漏掉谱的这一端，迭代就留下，和声速CFL下块预处理把误差每个sweep放大100倍是同一类缺口。把层数锁在3层，512 × 256的网格下最粗的网格是128 × 64。宽度为128个单元的误差波在那层网格上仍然是低频。没有人去擦它。

多重网格并非免费的地方

就此收尾会读成"永远用多重网格"。论文也写下了另一面。

增加层数是有代价的。问题每变大一次就要多叠一层才能维持最优收敛，而那份设置开销会随核数增长。在hypre或GAMG这类代数多重网格中，coarsening本身就是通信密集的。在Kollet等(2010)的结果里，solver setup的效率在16,384核附近几乎掉到0。迭代次数是平的，设置时间却吃掉了整体。也不是所有PDE都适合这套层次。像地下水流这样属于椭圆型、抛物型且初始猜测离解很远时，多重网格才划算。以对流为主导的地球化学输运正相反：时间步很小，矩阵强烈对角占优，点Jacobi和ILU(0)已经够用。没有理由去叠层。把多重网格用在描述多相流的方程组上，仍然是一个研究课题。

预处理的花费也挂在网格形状上。拉伸网格上各方向的耦合强度差别很大，像linelet预处理器那样只把一个方向解精确，更便宜。

PFLOTRAN作者要人先看strong scaling的理由#

论文的立场很清楚。实际工作中给你的是一个固定区域，你想多花核把同一个问题解得更快，回答这个问题的是strong scaling。可多重网格的求解器论文往往只给weak scaling。多重网格在问题变大时才突出，这个选择并不意外。问题在于只给一侧，就无法评估它相对于Krylov求解器的strong scaling表现。反过来，只给strong scaling的论文也有同样的毛病。

PFLOTRAN自身的strong scaling在铜浸出场景下到4096核几乎理想，然后拐弯。拐点处每核持有768个自由度。3.1 M自由度的问题和49 K的问题都在同一比值处拐弯。这个数字在实务上比整条效率曲线更有用。它的意思是：每核自由度不要掉到1000以下。

迭代涨到1285时，MPI拿不回0.163#

64核那一行，点Jacobi是1285次，效率0.163。从210到1285的倍数是6.12，效率是倒数， $210/1285=0.1634$ 。每核单元数每一行都是2048。一次迭代的花费相同，墙钟就跟着迭代次数走。把消息削掉，这6.12倍也不会回来。

分叉在层数。锁在3层，次数从10到47，比值4.7，效率0.213。层数跟着网格一起加，从9到10，比值1.11，效率0.900。同一份131072个自由度上，1285与10之比是128.5倍。曲线一斜，图注先写是扩大了区域还是细化了网格。迭代列涨了6.12倍，要动的是Krylov线性解法。迭代是平的、只有时间在涨，那才是通信这一侧。