并行计算 Parallel Computing

1. 并行计算 Parallel Computing#

本章内容

介绍 Firedrake 如何通过 MPI 进行并行计算: 用 mpiexec 运行脚本, 在 Jupyter 中通过 ipyparallel 执行并行代码并对照串行与并行的行为差异, 以及并行编程的常见陷阱. 部分单元需要额外安装 ipyparallel (pip install ipyparallel).

Overview

How Firedrake computes in parallel with MPI: running scripts with mpiexec, running parallel code inside Jupyter via ipyparallel while comparing serial and parallel execution, and common pitfalls of parallel programming. Some cells need the ipyparallel package (pip install ipyparallel).

1.1. Firedrake 中并行计算简介#

Firedrake 使用 MPI 实现分布式内存并行, 并采用单程序多数据 (SPMD, single program, multiple data) 模式: 所有进程运行同一份程序, 但分别处理全局问题的不同部分.

在并行计算中, Firedrake 借助 PETSc 的 DMPlex 表示和管理网格. 网格会被划分为若干子区域并分配给不同进程; 每个进程主要存储和计算本地网格、自由度以及矩阵和向量中由自己负责的部分. 为完成分区边界附近的计算, 进程还会保存必要的 halo (或 ghost) 数据. 不同进程共享的数据通过 MPI 通信进行更新和同步.

这些并行细节大多由 Firedrake 和 PETSc 自动处理. 因此, 同一份 Firedrake 程序通常既可以串行运行, 也可以通过 MPI 启动多个进程并行运行, 而无须修改有限元问题的主体代码.

需要特别注意的是, 每个进程都会执行整份 Python 程序. 因此, 普通的 print 会被每个进程分别执行, 可能产生重复输出; 直接使用常规绘图库或自行读写文件时, 也可能出现重复操作或文件冲突. 这类与输入输出有关的操作通常需要显式指定由某个进程执行, 或者使用 Firedrake、PETSc 提供的并行安全接口.

1.1.1. 并行启动程序#

Firedrake 启动并行计算只需在终端中运行以下命令:

mpiexec -n <number-of-processes> python3 /path/to/your/script.py

下面以求解 Poisson 方程的脚本 py/poisson.py 为例, 其完整内容如下:

"""并行求解 Poisson 方程的最小示例.

在 firedrake 目录下运行:

    mpiexec -n 2 python3 py/poisson.py
"""
from firedrake import *
from firedrake.petsc import PETSc

N = 4
test_mesh = RectangleMesh(N, N, 1, 1)
x, y = SpatialCoordinate(test_mesh)
f = sin(pi*x)*sin(pi*y)
g = Constant(0)

V = FunctionSpace(test_mesh, 'CG', degree=1)

u, v = TrialFunction(V), TestFunction(V)
a = inner(grad(u), grad(v))*dx
L = inner(f, v)*dx

bc = DirichletBC(V, g=g, sub_domain='on_boundary')

u_h = Function(V, name='u_h')
solve(a == L, u_h, bcs=bc)

# assemble 是集合操作, 所有进程一起调用; Print 只由 0 号进程输出
PETSc.Sys.Print(f'[{COMM_WORLD.size} process(es)] integral of u_h: '
                f'{assemble(u_h*dx):.10e}')

output = VTKFile('data/result.pvd')
output.write(u_h)

若使用 2 个进程进行计算, 在激活了 Firedrake 环境的终端中, 切换到本笔记本所在的 firedrake 目录后运行:

mpiexec -n 2 python3 py/poisson.py

脚本会先打印 u_h 的积分值 (串行与并行运行得到的值在打印的精度内相同, 更细致的比较见如何验证并行结果是对的), 再把解保存到 data/result.pvd. 在 data 目录下会生成一个 result 文件夹, 计算结果保存在该文件夹中, 每个进程对应一个结果文件 (2 个进程时为 result_0_0.vtu 和 result_0_1.vtu). 索引是两级的: result.pvd 指向 result/result_0.pvtu, 后者才列出各进程的 .vtu 文件; 用 ParaView 打开 result.pvd 即可, 不必关心这两级的细节.

如本节开头所述, 每个进程只负责整体网格的一个子区域, 因此每个进程的结果文件仅包含该进程对应区域的结果.

1.1.2. 输出与进程控制#

1.1.2.1. 并行时的输出#

参考 py/intro_utils.py

第一个进程输出 (其他进程的调用会被忽略)

PETSc.Sys.Print("This is the message will only show once!")

多个进程同步输出

rank, size = COMM_WORLD.rank, COMM_WORLD.size
PETSc.Sys.syncPrint(f"[{rank}/{size}] This is the message from rank {rank}!")
if rank == 0:
     PETSc.Sys.syncPrint(f"[{rank}/{size}] Message from rank {rank}!")
PETSc.Sys.syncFlush()

1.1.2.2. 仅在指定进程上运算#

有些事情只需要做一次, 例如打印汇总信息、写日志文件. 可以用条件语句让它们只在某个进程 (通常是 0 号进程) 上执行:

u_int = assemble(u_h*dx)                  # 集合操作: 所有进程一起调用
if COMM_WORLD.rank == 0:
    print(f'integral of u_h: {u_int}')    # 纯本地操作: 可以只在 0 号进程执行

要特别小心的是, 只有纯本地操作 (普通 print、写文本文件、对已经算好的数值做后处理等) 才能放进这样的条件分支; 凡是涉及 MPI 集合通信的操作 (如 assemble、solve、norm, 甚至 triplot 绘图) 都必须所有进程一起调用, 只让部分进程执行会使程序永久挂起, 详见并行常见陷阱.

1.2. 在 Jupyter 中执行并行代码#

在 Jupyter 交互环境中, 我们可以方便地对串行代码进行测试, 从而进行快速验证. 为了在 Jupyter 环境下运行并行程序, ipyparallel 包提供了强大的支持, 它能够帮助我们在 Jupyter 中验证并行代码, 同时加深我们对 Firedrake 并行机制的理解.

1.2.1. `ipyparallel` 简介#

ipyparallel 的工作方式是: 在 notebook 之外另行启动若干独立的 Python 进程 (称为引擎, engine), notebook 通过一个客户端 (client) 把代码发送到这些引擎上执行, 再把结果取回显示. 引擎才是真正做计算的进程, notebook 自身的内核仍然是单进程的.

下面创建并连接一个包含 2 个引擎的集群: engines="mpi" 表示用 MPI 启动引擎, 使这 2 个引擎构成一个 MPI 通信组 (效果相当于 mpiexec -n 2), Firedrake 在其中即可按前面介绍的方式并行计算; start_and_connect_sync() 启动集群并等待所有引擎就绪, 返回与之通信的客户端.

Tip

这里必须写 engines="mpi". 若改用 ipp.Cluster(profile="mpi"), 在本机缺少相应 profile 时会静默退化为若干互不通信的单进程引擎 (每个引擎中 COMM_WORLD 的大小都是 1), 看起来能运行, 但并没有真正并行.

import ipyparallel as ipp

cluster = ipp.Cluster(engines="mpi", n=2)
client = cluster.start_and_connect_sync()

集群就绪后, 用 ipyparallel 提供的魔法指令 %%px 可以把整个单元的代码发送到所有引擎上同时执行. --block 表示同步执行: 等所有引擎都运行完毕后再继续执行后面的单元, 这样输出更便于对照.

1.2.2. 执行并行代码#

下面的单元在 2 个引擎上分别导入 Firedrake 并生成网格, 然后打印各自的进程编号和总进程数: mesh.comm 是该网格所在的 MPI 通信子, 由它取得的 size 为 2, 说明两个引擎确实位于同一个 MPI 通信组中. 输出仍按前面介绍的方式控制: syncPrint 配合 syncFlush 按进程编号顺序输出每个进程的信息, Print 只输出一次.

%%px --block
from firedrake import *
from firedrake.petsc import PETSc

N = 4
mesh = RectangleMesh(N, N, 1, 1)

rank, size = mesh.comm.rank, mesh.comm.size
PETSc.Sys.syncPrint(f"[{rank}/{size}] This is the message from rank {rank}!")
PETSc.Sys.syncFlush()

PETSc.Sys.Print(f"[{rank}/{size}] This message only from rank 0!")

[stdout:0] [0/2] This is the message from rank 0!
[1/2] This is the message from rank 1!
[0/2] This message only from rank 0!

提示: 输出中的 [stdout:0] 前缀

%%px 显示多个引擎的输出时, 会为每个引擎加上 [stdout:<引擎编号>] 前缀. 该前缀由 ipyparallel 在内部硬编码 (见 AsyncResult.display_outputs), 没有提供修改文本的配置项; %%px 只能用 --group-outputs 改变分组顺序, 或用 -t/--targets 限制参与执行的引擎. 另外, 连续来自同一引擎的多行输出只在开头显示一次前缀.

上面的例子中所有信息都由 0 号引擎打印, 因此只出现一个 [stdout:0]: PETSc.Sys.syncPrint 把各进程的消息汇总到 0 号进程, 再由 syncFlush 一次输出, 各行开头的 [rank/size] 是我们自己加的进程标识. 需要区分进程时, 这种写法比依赖 ipyparallel 的前缀更可控.

若想完全自定义显示方式, 可以不用 %%px 的自动显示, 而是直接处理执行结果:

ar = client[:].execute(
    "from mpi4py import MPI; print(f'rank {MPI.COMM_WORLD.rank}')", block=True
)
for eid, out in zip(ar.engine_id, ar.stdout):
    for line in out.rstrip('\n').split('\n'):
        print(f'<engine {eid}> {line}')

输出为:

<engine 0> rank 0
<engine 1> rank 1

1.2.3. `%%px` 的更多用法#

前面用的一直是最简单的 %%px --block, 它把整个单元发给所有引擎并等待全部完成. 这里补充几个常用的参数.

-t/--targets 用来只让部分引擎执行, 例如 --targets 0 表示只在 0 号引擎上运行. 下面用它打印引擎的工作目录:

%%px --block --targets 0
import os
print(f'cwd basename            : {os.path.basename(os.getcwd())}')
print(f'py/intro_utils.py exists: {os.path.exists("py/intro_utils.py")}')

[stdout:0] cwd basename            : firedrake
py/intro_utils.py exists: True

两行输出合起来说明: 引擎当前目录的名字就是 firedrake, 而且在这个目录下用相对路径 py/intro_utils.py 确实找得到文件. 也就是说, 引擎的工作目录是启动本 notebook 的那个目录, 而不是引擎自己的某个临时目录. 这一点很容易被忽略, 却是本章能跑通的前提: 后面并行单元里的 from py.intro_utils import triplot 之所以在引擎里也能成功, 靠的正是这一点; 换个目录启动 notebook, 同样的代码就会在引擎里抛出 ModuleNotFoundError. 引擎里用相对路径读写文件同理, 例如并行下的 I/O 一节用 CheckpointFile 写出的 data/parallel_checkpoint.h5, 也落在这个目录下.

要提醒的是, --targets 只适合放纯本地的操作. 单元里一旦出现集合操作 (assemble、solve、norm 等), 只让部分引擎执行就会挂起, 原因见并行常见陷阱.

--noblock 则相反: 把代码发出去就立刻返回, 不等引擎算完, 单元的值是一个 AsyncResult 对象. 这样可以在引擎慢慢算的同时, 继续在本地内核里做别的事. 稍后用 %pxresult 取回并显示上一次异步执行的输出:

%%px --noblock
import time
time.sleep(1)
slow_value = 'the slow cell has finished'
PETSc.Sys.Print(slow_value)

<AsyncResult(%px): pending>

%pxresult

[stdout:0] 
the slow cell has finished

--noblock 那个单元立刻返回了一个 <AsyncResult(%px): pending>, 说明它没有等待; 直到执行 %pxresult, 才阻塞着等引擎算完并显示输出. 若想在代码里显式等待, 可以把 AsyncResult 存到变量里再等它: %%px --noblock --out ar, 然后在后面的单元中调用 ar.wait().

另外还有 %autopx: 打开之后每一个单元都会自动发到引擎上执行, 直到再次执行 %autopx 关闭. 它省去了反复写 %%px, 但也很容易忘记关闭, 使本该在本地运行的单元跑到引擎上去, 因此本章没有使用它.

最后再强调一次: 引擎和本地内核是不同的进程, 变量互不相通. 要把本地的值送到引擎上, 用 client[:].push({'name': value}); 要把引擎上的结果取回本地, 用 client[0].pull('name', block=True) —— 后者在后面比较串行与并行结果时会实际用到.

1.3. 使用 `ipyparallel` 观察串行和并行过程#

下面每个小节都用几乎相同的代码运行两遍: 不加 %%px 的单元在 notebook 自身的本地内核 (单进程) 中运行, 加了 %%px --block 的单元则发送到前面启动的集群引擎上运行. 两边各自维护自己的一套变量, 互不干扰: 后续不带 %%px 的单元使用的是本地内核里最近一次定义的 mesh/V1/V2/W 等对象, 带 %%px 的单元使用的则是集群引擎里最近一次通过 %%px 定义的同名对象. 通过对照两侧的输出, 可以直观看到同一段 Firedrake 代码在串行和并行下的行为差异.

为便于区分, 下文每个代码单元前都有一行粗体标注: 串行 (本地内核) 表示该单元在本地内核中运行, 并行 (2 个引擎) 表示该单元以 %%px --block 开头、在集群引擎上运行; 当并行单元的代码与其串行版本完全相同 (仅首行多出 %%px --block) 时, 默认折叠并行单元的代码, 需要时点击展开查看.

1.3.1. 生成网格并画出网格#

串行 (本地内核): 先在本地内核里生成 4×4 网格并画出, 作为对照基准; mesh.topology_dm.view() 输出该网格对应 DMPlex 的拓扑信息.

from firedrake import *
from py.intro_utils import triplot
from firedrake.petsc import PETSc

import matplotlib.pyplot as plt

N = 4
mesh = RectangleMesh(N, N, 1, 1)
mesh.topology_dm.view()
fig, axes = plt.subplots(figsize=[4, 3])
triplot(mesh, axes=axes)
axes.set_xlim([-0.1, 1.1])
axes.set_ylim([-0.1, 1.1])
fig.tight_layout()

DM Object: DM_0x31e231d0_0 1 MPI process
  type: plex
DM_0x31e231d0_0 in 2 dimensions:
  Number of 0-cells per rank: 25
  Number of 1-cells per rank: 56
  Number of 2-cells per rank: 32
Labels:
  celltype: 3 strata with value/size (0 (25), 1 (56), 3 (32))
  depth: 3 strata with value/size (0 (25), 1 (56), 2 (32))
  marker: 1 strata with value/size (1 (32))
  Face Sets: 4 strata with value/size (1 (9), 2 (9), 3 (9), 4 (9))
  exterior_facets: 1 strata with value/size (1 (32))
  interior_facets: 1 strata with value/size (1 (63))
  pyop2_core: 1 strata with value/size (1 (113))
  pyop2_owned: 0 strata with value/size ()
  pyop2_ghost: 0 strata with value/size ()

../_images/219b2e05715375019baccc64002b8af5b59a33cca2c812fa343725b78ae7a470.png

并行 (2 个引擎): 在集群的 2 个引擎上重复相同代码; mesh.comm 此时是各引擎自己的 MPI 通信子, 生成的网格会被 DMPlex 自动划分到 2 个进程.

[stdout:0] DM Object: DM_0x2773a8e0_1 2 MPI processes
  type: plex
DM_0x2773a8e0_1 in 2 dimensions:
  Number of 0-cells per rank: 19 19
  Number of 1-cells per rank: 38 38
  Number of 2-cells per rank: 20 20
Labels:
  depth: 3 strata with value/size (0 (19), 1 (38), 2 (20))
  celltype: 3 strata with value/size (0 (19), 1 (38), 3 (20))
  marker: 1 strata with value/size (1 (19))
  Face Sets: 3 strata with value/size (2 (9), 3 (5), 4 (7))
  exterior_facets: 1 strata with value/size (1 (19))
  interior_facets: 1 strata with value/size (1 (47))
  pyop2_core: 1 strata with value/size (1 (26))
  pyop2_owned: 1 strata with value/size (1 (26))
  pyop2_ghost: 1 strata with value/size (1 (25))

[output:0]

../_images/c27bb12b52d820e9d9c70ada929df14c7634a55b25292ad82d887d4d715d7646.png

[output:1]

../_images/342a93706dbd19563dfe48f7461ebca6eb1b5c42b5a71fc595991935c47102ff.png

可以看到上面并行中两个网格图分别只是整体网格的一部分, 且边界处有重叠 (ghost/halo 单元): 这些重叠单元使每个进程在组装局部矩阵时也能访问到相邻单元的信息.

两段 mesh.topology_dm.view() 输出中的数字也值得对照:

串行时网格共有 25 个顶点 (0-cells)、56 条边 (1-cells)、32 个单元 (2-cells); 并行时每个进程各有 19 个顶点、38 条边、20 个单元. 注意 20 + 20 = 40 > 32: 多出的 8 个单元正是分区边界两侧互相重叠的 ghost 单元, 顶点和边同理.
pyop2_core、pyop2_owned、pyop2_ghost 把 DMPlex 的点 (顶点、边、单元的统称) 分为三类 (含义见后文): 串行时全部 113 = 25 + 56 + 32 个点都是 core; 并行时 0 号进程为 26、26、25, 三者之和 77 = 19 + 38 + 20.
Labels 一段只显示 0 号进程的信息, 并行输出中 Face Sets 只列出 2、3、4 三个值: 编号为 1 的那条边界完全落在 1 号进程上, 0 号进程没有它的任何一段. 这正是本章从 py/intro_utils.py 导入 triplot 包装函数的原因: 并行绘图时, 拿不到某条边界的进程会收到 Firedrake 的 “Subdomain is empty” 警告, 包装函数把这条无害的警告过滤掉了.

1.3.2. 分区重叠 (overlap) 的显式控制#

上面看到的 ghost 单元并不是固定不变的, 它由建网格时的 distribution_parameters 控制. 其中 overlap_type 是一个二元组 (类型, 层数): Firedrake 的默认值是 (DistributedMeshOverlapType.FACET, 1), 表示每个进程除了自己拥有的单元, 再多保留一层”与本进程单元共享一条边 (面) 的邻居单元”; 改成 (DistributedMeshOverlapType.NONE, 0) 就完全关闭重叠, 改成 (DistributedMeshOverlapType.FACET, 2) 则保留两层.

下面用同一个 4×4 网格把三种设置各建一次, 只看单元的计数: cell_set.size 是本进程拥有的单元数, cell_set.total_size 还包含 ghost 单元, 两者之差就是 ghost 单元数 (这个三段划分的含义见后文). 注意这里用的是新变量名 mesh_ovlp, 后面各小节继续使用的仍是前面生成的 mesh. 单元开头之所以还要重新取一次 rank 和 size, 原因见 from firedrake import * 会覆盖 rank 等变量.

并行 (2 个引擎)

%%px --block
rank, size = mesh.comm.rank, mesh.comm.size   # 重新取一次: firedrake 也导出了名为 rank 的函数
for label, overlap in (('(NONE, 0) ', (DistributedMeshOverlapType.NONE, 0)),
                       ('(FACET, 1)', (DistributedMeshOverlapType.FACET, 1)),
                       ('(FACET, 2)', (DistributedMeshOverlapType.FACET, 2))):
    mesh_ovlp = RectangleMesh(N, N, 1, 1,
                              distribution_parameters={'overlap_type': overlap})
    cells = mesh_ovlp.cell_set
    PETSc.Sys.syncPrint(f'[{rank}/{size}] overlap_type={label} '
                        f'owned cells: {cells.size:>3}, '
                        f'total cells: {cells.total_size:>3}, '
                        f'ghost cells: {cells.total_size - cells.size:>3}')
    PETSc.Sys.syncFlush()

[stdout:0] [0/2] overlap_type=(NONE, 0)  owned cells:  16, total cells:  16, ghost cells:   0
[1/2] overlap_type=(NONE, 0)  owned cells:  16, total cells:  16, ghost cells:   0
[0/2] overlap_type=(FACET, 1) owned cells:  16, total cells:  20, ghost cells:   4
[1/2] overlap_type=(FACET, 1) owned cells:  16, total cells:  20, ghost cells:   4
[0/2] overlap_type=(FACET, 2) owned cells:  16, total cells:  24, ghost cells:   8
[1/2] overlap_type=(FACET, 2) owned cells:  16, total cells:  24, ghost cells:   8

对照三行输出可以看到: owned cells 一列始终是 16, 变化的是 total cells —— 关掉重叠后它从 20 变回 16, ghost 单元为 0. 这正是 overlap 的含义: 它只增加每个进程本地存储的单元数, 不改变每个进程真正拥有的单元数. 而 16 + 16 = 32 恰是整张网格的单元总数, 上一段”20 + 20 = 40 > 32, 多出的 8 个是 ghost 单元”这句话, 到这里就变成了可以直接读出来的数字. 把层数加到 2, ghost 单元又从 4 增加到 8: 重叠层越厚, 每个进程需要额外保存和同步的数据就越多.

有两点需要说明:

这个对照只在并行下才有意义. 串行运行时整张网格都在一个进程上, 没有分区也就没有重叠, 三种设置得到的计数完全相同 (都是 32 个单元、0 个 ghost).
不要因为”关掉重叠更省内存”就随手改这个参数. 像本章这样只含单元内部积分 (dx) 的 CG 问题, 关掉重叠后仍能算出正确结果 (实测串行与并行的积分值一致), 因为 CG 的自由度交换只依赖分区边界上的节点; 但只要计算需要用到相邻单元的数据 —— 例如间断 Galerkin 方法的内部面积分 dS (见 Poisson 方程 II 的 SIPG 一节) —— 缺少重叠层就会出错, 实测甚至会让程序直接崩溃, 而不是给出清晰的报错. 默认的 (FACET, 1) 正是为了让这些常见用法开箱即用.

1.3.3. 定义变分问题#

定义一个简单的混合变分问题 (向量分量 \(u_1\) 用 \([P_1]^2\), 标量分量 \(u_2\) 用 \(P_2\)), 后续几个小节都基于这里定义的 mesh、V1、V2、W 展开.

串行 (本地内核)

V1 = VectorFunctionSpace(mesh, 'CG', 1)
V2 = FunctionSpace(mesh, 'CG', 2)
W = MixedFunctionSpace([V1, V2])  # W = V1*V2
u1, u2 = TrialFunctions(W)
v1, v2 = TestFunctions(W)
a = dot(u1, v1)*dx + u2*v2*dx

x, y = SpatialCoordinate(mesh)
f = dot(as_vector((sin(x), cos(y))), v1)*dx + cos(y)*v2*dx

bc = DirichletBC(W.sub(0), 0, 1)
uh = Function(W)
problem = LinearVariationalProblem(a, f, uh, bcs=bc)

并行 (2 个引擎)

两段代码完全相同, 区别只在于后者由 %%px 分发到两个引擎上执行: W 在每个进程上只保存该进程负责的那部分自由度.

1.3.4. 函数空间维度#

函数空间的自由度个数可以从几个角度查看, 下面以 V1 为例. V1.dim() 返回全局自由度总数, 所有进程上一致; node_count 和 dof_count 则是含 halo 的本地计数: 除了本进程拥有的节点/自由度外, 还包括从相邻进程同步来的 ghost 节点/自由度.

串行 (本地内核)

rank, size = mesh.comm.rank, mesh.comm.size
PETSc.Sys.Print(f'Number of dofs of V1: {V1.dim()}')
PETSc.Sys.syncPrint(f'[{rank}/{size}] V1 Node count: {V1.node_count}; V1 Dof count: {V1.dof_count}')
PETSc.Sys.syncFlush()

Number of dofs of V1: 50
[0/1] V1 Node count: 25; V1 Dof count: 50

并行 (2 个引擎)

[stdout:0] Number of dofs of V1: 50
[0/2] V1 Node count: 19; V1 Dof count: 38
[1/2] V1 Node count: 19; V1 Dof count: 38

串行时 node_count 就是全局节点数 25, dof_count 是它的 2 倍 (每个节点上有 2 个分量). 并行时两个进程的 node_count 恰好都是 19, 但组成并不相同: 0 号进程是 10 个自身拥有的节点加 9 个 ghost 节点, 1 号进程是 15 个加 4 个. 拥有的节点数之和 10 + 15 = 25 恰为全局节点数; 而 19 + 19 = 38 > 25, 多出的正是分区边界两侧互为 ghost 的节点. 各进程拥有的节点数可以从后文的 node_set.size 读出.

为进一步观察本地和 ghost 自由度的区别, 下面创建一个 Function, 分别用 PETSc Vec 的只读视图 (dat.vec_ro) 和 dat 的只读数组视图 (data_ro、data_ro_with_halos) 查看它的本地大小. 顺带一提, 只读访问数据请用 data_ro 系列属性; 可写的 data 会把 halo 标记为失效, 使下一次访问带 halo 的数据时触发同步通信.

串行 (本地内核)

f_v1 = Function(V1)
with f_v1.dat.vec_ro as vec:
    PETSc.Sys.syncPrint(f'[{rank}/{size}] vec sizes: {vec.getSizes()}')
    PETSc.Sys.syncFlush()
PETSc.Sys.syncPrint(f'[{rank}/{size}] data: {f_v1.dat.data_ro.shape}, '
                    f'data_with_halos: {f_v1.dat.data_ro_with_halos.shape}')
PETSc.Sys.syncFlush()

[0/1] vec sizes: (50, 50)
[0/1] data: (25, 2), data_with_halos: (25, 2)

并行 (2 个引擎)

[stdout:0] [0/2] vec sizes: (20, 50)
[1/2] vec sizes: (30, 50)
[0/2] data: (10, 2), data_with_halos: (19, 2)
[1/2] data: (15, 2), data_with_halos: (19, 2)

vec.getSizes() 返回的二元组是 (本地长度, 全局长度): 串行时为 (50, 50); 并行时全局长度不变, 本地长度 20 和 30 分别是两个进程拥有的自由度数 (10 个和 15 个节点, 每个节点 2 个分量). data_ro 的形状与之对应: (10, 2) 和 (15, 2) 只含本进程拥有的自由度, 10 × 2 = 20 正是 Vec 的本地长度; data_ro_with_halos 的形状 (19, 2) 则与前面的 node_count 一致, 额外包含了 ghost 自由度. 串行时没有 ghost 区域, 两种视图大小相同.

1.3.5. 节点集和自由度数据集#

PyOP2 用 node_set 和 dof_dset 描述一个函数空间在各进程上如何切分: node_set 对应网格节点 (自由度承载点), dof_dset 在此基础上进一步考虑每个节点上的分量数 (如向量场每个节点有多个自由度). 下面打印它们的 str 和 repr: str 中的 size 是本进程拥有的节点数, repr 中的三元组是 (core_size, size, total_size), 含义见下一小节.

串行 (本地内核)

PETSc.Sys.syncPrint(f'[{rank}/{size}] V1: {str(V1.node_set)}')
PETSc.Sys.syncPrint(f'[{rank}/{size}]     {repr(V1.node_set)}')
PETSc.Sys.syncPrint(f'[{rank}/{size}] V1: {str(V1.dof_dset)}')
PETSc.Sys.syncPrint(f'[{rank}/{size}]     {repr(V1.dof_dset)}')
PETSc.Sys.syncFlush()

[0/1] V1: OP2 Set: set_#x7fe9418f3800 with size 25
[0/1]     Set((np.int64(25), np.int64(25), np.int64(25)), 'set_#x7fe9418f3800')
[0/1] V1: OP2 DataSet: None_nodes_dset on set OP2 Set: set_#x7fe9418f3800 with size 25, with dim (2,), False
[0/1]     DataSet(Set((np.int64(25), np.int64(25), np.int64(25)), 'set_#x7fe9418f3800'), (2,), 'None_nodes_dset', False)

并行 (2 个引擎)

[stdout:0] [0/2] V1: OP2 Set: set_#x7fe08473dc40 with size 10
[0/2]     Set((np.int64(5), np.int64(10), np.int64(19)), 'set_#x7fe08473dc40')
[0/2] V1: OP2 DataSet: None_nodes_dset on set OP2 Set: set_#x7fe08473dc40 with size 10, with dim (2,), False
[0/2]     DataSet(Set((np.int64(5), np.int64(10), np.int64(19)), 'set_#x7fe08473dc40'), (2,), 'None_nodes_dset', False)
[1/2] V1: OP2 Set: set_#x7fbda1a61af0 with size 15
[1/2]     Set((np.int64(10), np.int64(15), np.int64(19)), 'set_#x7fbda1a61af0')
[1/2] V1: OP2 DataSet: None_nodes_dset on set OP2 Set: set_#x7fbda1a61af0 with size 15, with dim (2,), False
[1/2]     DataSet(Set((np.int64(10), np.int64(15), np.int64(19)), 'set_#x7fbda1a61af0'), (2,), 'None_nodes_dset', False)

1.3.5.1. Set 与 DataSet 的各种 size#

参考: firedrakeproject/firedrake

Set 中的元素按下标分成三段, 段的边界由 core_size、size、total_size 三个属性给出:

[0, core_size): 本进程拥有、且计算时不依赖任何 halo 数据的元素;
[core_size, size): 本进程拥有的其余元素, 它们与其他进程相邻, 计算时可能用到 halo 数据;
[size, total_size): ghost 元素, 由其他进程拥有, 本进程只保存副本.

属性 sizes 即三元组 (core_size, size, total_size). 对照下面的输出: 0 号进程的 node_set 三元组为 (5, 10, 19), 即它拥有 10 个节点 (其中 5 个完全不与 ghost 节点相邻), 加上 9 个 ghost 节点后本地共存储 19 个; 串行时没有 ghost, 三个数都等于全局节点数 25.

串行 (本地内核)

node_set = V1.node_set
msg = f'core size: {node_set.core_size}, size: {node_set.size}, total size: {node_set.total_size}'
PETSc.Sys.syncPrint(f'[{rank}/{size}] {msg}')
PETSc.Sys.syncFlush()

[0/1] core size: 25, size: 25, total size: 25

并行 (2 个引擎)

[stdout:0] [0/2] core size: 5, size: 10, total size: 19
[1/2] core size: 10, size: 15, total size: 19

串行 (本地内核)

dof_dset = V1.dof_dset
size_msg = f'core size: {dof_dset.core_size}, size: {dof_dset.size}, total size: {dof_dset.total_size}'
# dim: shape tuple of the values for each element, cdim: product of dim tuple
dim_msg = f'dim: {dof_dset.dim}, cdim: {dof_dset.cdim}'
PETSc.Sys.syncPrint(f'[{rank}/{size}] {size_msg}, {dim_msg}')
PETSc.Sys.syncFlush()

[0/1] core size: 25, size: 25, total size: 25, dim: (2,), cdim: 2

并行 (2 个引擎)

[stdout:0] [0/2] core size: 5, size: 10, total size: 19, dim: (2,), cdim: 2
[1/2] core size: 10, size: 15, total size: 19, dim: (2,), cdim: 2

dof_dset 的三段划分与 node_set 完全一致 (它计数的仍是节点), 额外的 dim=(2,) 与 cdim=2 表示每个节点上有 2 个自由度分量: 自由度数等于节点数乘以 cdim, 例如 0 号进程拥有 10 × 2 = 20 个自由度, 全局共 25 × 2 = 50 个.

1.3.5.2. DataSet 的索引集 (IS)#

IS (Index Set) 是 PETSc 用来描述一组整数下标的对象, 这里用来表示混合空间 W 中各子空间对应的自由度下标.

field_ises: firedrakeproject/firedrake A list of PETSc ISes defining the global indices for each set in the DataSet. Used when extracting blocks from matrices for solvers.

local_ises: A list of PETSc ISes defining the local indices for each set in the DataSet. Used when extracting blocks from matrices for assembly.

串行 (本地内核)

local_ises_msg = f'{[_.getIndices() for _ in W.dof_dset.local_ises]}'
PETSc.Sys.syncPrint(f'[{rank}/{size}] {local_ises_msg}')
PETSc.Sys.syncFlush()

[0/1] [array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
       34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49],
      dtype=int32), array([ 50,  51,  52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,
        63,  64,  65,  66,  67,  68,  69,  70,  71,  72,  73,  74,  75,
        76,  77,  78,  79,  80,  81,  82,  83,  84,  85,  86,  87,  88,
        89,  90,  91,  92,  93,  94,  95,  96,  97,  98,  99, 100, 101,
       102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114,
       115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127,
       128, 129, 130], dtype=int32)]

并行 (2 个引擎)

[stdout:0] [0/2] [array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
       34, 35, 36, 37], dtype=int32), array([38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,
       55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71,
       72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88,
       89, 90, 91, 92, 93, 94], dtype=int32)]
[1/2] [array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
       34, 35, 36, 37], dtype=int32), array([38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,
       55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71,
       72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88,
       89, 90, 91, 92, 93, 94], dtype=int32)]

串行 (本地内核)

field_ises_msg = f'{[_.getIndices() for _ in W.dof_dset.field_ises]}'
PETSc.Sys.syncPrint(f'[{rank}/{size}] {field_ises_msg}')
PETSc.Sys.syncFlush()

[0/1] [array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
       34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49],
      dtype=int32), array([ 50,  51,  52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,
        63,  64,  65,  66,  67,  68,  69,  70,  71,  72,  73,  74,  75,
        76,  77,  78,  79,  80,  81,  82,  83,  84,  85,  86,  87,  88,
        89,  90,  91,  92,  93,  94,  95,  96,  97,  98,  99, 100, 101,
       102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114,
       115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127,
       128, 129, 130], dtype=int32)]

并行 (2 个引擎)

[stdout:0] [0/2] [array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19], dtype=int32), array([20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
       37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
       54, 55], dtype=int32)]
[1/2] [array([56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72,
       73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85], dtype=int32), array([ 86,  87,  88,  89,  90,  91,  92,  93,  94,  95,  96,  97,  98,
        99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111,
       112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124,
       125, 126, 127, 128, 129, 130], dtype=int32)]

local_ises 和 field_ises 覆盖的自由度集合并不相同, 它们不是同一组下标的局部/全局两种写法:

local_ises 用含 halo 的本地编号给出各子空间的下标区间: 本例中两个进程都是 V1 占 [0, 38)、V2 占 [38, 95), 共 95 = 38 + 57 个下标 (对应含 halo 的本地自由度); 本地组装时用它从矩阵中提取分块.
field_ises 只包含本进程拥有的自由度, 给出的是全局编号: 0 号进程为 V1 占 [0, 20)、V2 占 [20, 56), 1 号进程为 V1 占 [56, 86)、V2 占 [86, 131); 求解器从全局矩阵中提取分块时用它.

field_ises 的输出还揭示了全局自由度的排列规律: 先按进程、再按场. 0 号进程拥有的全部自由度 (V1 和 V2 的都在内) 占据全局编号区间 [0, 56), 1 号进程的占据 [56, 131). 正因如此, 想按场提取子矩阵/子向量时, 不能简单地把全局向量切成前后两段, 而必须借助这样的索引集; Navier–Stokes 方程一章中基于 fieldsplit 的预条件子正是建立在这套机制上.

1.3.5.3. 局部到全局的映射#

并行时每个进程只保存网格和自由度的一部分, 需要一套局部编号到全局编号的映射, 才能正确组装全局矩阵/向量, 或在 DMPlex 层面同步不同进程间共享的实体.

下面直接查看 Halo 维护的局部到全局编号映射数组: 串行运行时该映射是恒等映射 (每个局部编号就等于全局编号), 并行运行时才能看到两者的差异.

串行 (本地内核)

halo = V1.dof_dset.halo
PETSc.Sys.syncPrint(f'[{rank}/{size}] {halo.local_to_global_numbering}')
PETSc.Sys.syncFlush()
# halo.dm.getGlobalSection().view()

[0/1] [ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
 24]

并行 (2 个引擎)

[stdout:0] [0/2] [ 0  1  2  3  4  5  6  7  8  9 24 19 23 17 22 15 21 20 13]
[1/2] [10 11 12 13 14 15 16 17 18 19 20 21 22 23 24  9  8  7  6]

以 0 号进程的输出为例: 前 10 个是它拥有的节点, 全局编号恰为连续的 0 到 9 (全局编号先按进程排列); 后 9 个是 ghost 节点, 它们的全局编号 (24、19、23 等) 都落在 1 号进程拥有的区间 [10, 25) 内. 1 号进程正好相反: 拥有的 15 个节点编号为 10 到 24, 最后 4 个 ghost 节点的编号来自 0 号进程.

dof_dset.lgmap 是 PETSc 的 LGMap (local-to-global map) 对象, 负责把进程局部下标映射为全局下标. 下面用它的 apply 方法转换 local_ises 中的局部下标:

串行 (本地内核)

applied = [V1.dof_dset.lgmap.apply(_) for _ in V1.dof_dset.local_ises]
PETSc.Sys.syncPrint(f'[{rank}/{size}] {applied}')
PETSc.Sys.syncFlush()

[0/1] [array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
       34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49],
      dtype=int32)]

并行 (2 个引擎)

[stdout:0] [0/2] [array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 48, 49, 38, 39, 46, 47, 34, 35, 44, 45, 30, 31, 42, 43,
       40, 41, 26, 27], dtype=int32)]
[1/2] [array([20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
       37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 18, 19, 16, 17,
       14, 15, 12, 13], dtype=int32)]

lgmap.apply(...) 把给定的局部下标数组转换为全局下标; 若想直接拿到完整的局部到全局映射表, 可以用下面的 .indices.

串行 (本地内核)

PETSc.Sys.syncPrint(f'[{rank}/{size}] {V1.dof_dset.lgmap.indices}')
PETSc.Sys.syncFlush()

[0/1] [ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
 48 49]

并行 (2 个引擎)

[stdout:0] [0/2] [ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 48 49 38 39
 46 47 34 35 44 45 30 31 42 43 40 41 26 27]
[1/2] [20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
 44 45 46 47 48 49 18 19 16 17 14 15 12 13]

注意 lgmap 按自由度而非节点编号: 它以块大小 cdim = 2 建立, 节点全局编号乘以 2 即为该节点第一个自由度的全局编号. 与上一小节 local_to_global_numbering 的输出逐位对照: 0 号进程的第一个 ghost 节点全局编号为 24, 对应这里的自由度编号 48、49; 下一个 ghost 节点 19 对应 38、39, 依此类推.

1.3.5.4. DataSet 的布局向量#

layout_vec 是与 dof_dset 并行分布一致的 PETSc Vec, 可用于需要直接操作 PETSc 向量 (而非 Firedrake Function) 的场合; 这里查看它的本地和全局长度, 本地长度即本进程拥有的自由度数, 与前面 vec.getSizes() 的输出一致.

串行 (本地内核)

vec_msg = f'Local Size: {dof_dset.layout_vec.getLocalSize()}, Size: {dof_dset.layout_vec.getSize()}'
PETSc.Sys.syncPrint(f'[{rank}/{size}] {vec_msg}')
PETSc.Sys.syncFlush()

[0/1] Local Size: 50, Size: 50

并行 (2 个引擎)

[stdout:0] [0/2] Local Size: 20, Size: 50
[1/2] Local Size: 30, Size: 50

1.3.6. 矩阵组装#

自由度按进程分布之后, 组装得到的矩阵/向量在存储上也是分布式的; 下面这些函数负责根据 dof_dset 的分布信息分配矩阵稀疏模式 (sparsity) 和存储.

Firedrake method: ExplicitMatrixAssembler.allocate
PyOp2 class: Sparsity
PyOp2 function: build_sparsity
PyOp2 class: Mat

更多组装细节请看矩阵组装内核

1.4. 通信子: `COMM_WORLD` 与 `mesh.comm`#

本章出现过两种写法: 并行时的输出一节用的是 COMM_WORLD, 而生成网格之后的单元用的都是 mesh.comm. 两者的关系是:

COMM_WORLD 是 MPI 的全局通信子, 包含由 mpiexec (或本章的引擎) 启动的全部进程;
mesh.comm 是这张网格所在的通信子. 不特别指定时, Firedrake 就把网格建在 COMM_WORLD 上, 此时两者是同一个对象;
建网格时可以用 comm= 参数指定别的通信子, 这时 mesh.comm 才与 COMM_WORLD 不同.

因此, 与网格有关的代码优先用 mesh.comm: 默认情形下它和 COMM_WORLD 完全一样, 而一旦网格建在子通信子上, 只有 mesh.comm 才是对的.

下面先验证默认情形, 再用 COMM_WORLD.Split 按进程号把 2 个进程拆成 2 个各含 1 个进程的子通信子, 在每个子通信子上各建一张网格、各解一个右端项不同的问题.

并行 (2 个引擎)

%%px --block
PETSc.Sys.syncPrint(f'[{rank}/{size}] mesh.comm is COMM_WORLD: {mesh.comm is COMM_WORLD}')
PETSc.Sys.syncFlush()

subcomm = COMM_WORLD.Split(color=rank, key=0)     # 每个进程各成一组
mesh_sub = UnitSquareMesh(8, 8, comm=subcomm)
V_sub = FunctionSpace(mesh_sub, 'CG', 1)
u_sub, v_sub = TrialFunction(V_sub), TestFunction(V_sub)
x_sub, y_sub = SpatialCoordinate(mesh_sub)
uh_sub = Function(V_sub)
solve(inner(grad(u_sub), grad(v_sub))*dx
      == inner(Constant(rank + 1)*sin(pi*x_sub)*sin(pi*y_sub), v_sub)*dx,
      uh_sub, bcs=DirichletBC(V_sub, 0, 'on_boundary'))
PETSc.Sys.syncPrint(f'[{rank}/{size}] mesh_sub.comm size: {mesh_sub.comm.size}, '
                    f'local cells: {mesh_sub.num_cells()}, '
                    f'integral: {assemble(uh_sub*dx):.10e}')
PETSc.Sys.syncFlush()

[stdout:0] [0/2] mesh.comm is COMM_WORLD: True
[1/2] mesh.comm is COMM_WORLD: True
[0/2] mesh_sub.comm size: 1, local cells: 128, integral: 1.9751170317e-02
[1/2] mesh_sub.comm size: 1, local cells: 128, integral: 3.9502340634e-02

第一行确认了默认网格的 mesh.comm 就是 COMM_WORLD 本身 (is 判断为真). 随后 COMM_WORLD.Split 按进程号把两个进程分到不同的组, 每组只有 1 个进程: 于是 mesh_sub 在每个进程上都是一张完整的 8×8 网格 (128 个单元, 不再分区), 两个进程各自独立地求解自己那个问题. 两个积分值不同, 且第二个恰好是第一个的两倍 —— 因为方程是线性的, 而右端项差了 rank + 1 这个倍数.

这正是”空间并行 × 参数并行”两级分解的原型: 把全部进程分成若干组, 组内并行求解一个空间问题, 组间并行处理不同的参数、不同的时间片或不同的样本. Firedrake 为此提供了 Ensemble 类, 它替你完成上面这种 Split, 并管理两套通信子: ensemble.comm 用于组内的空间并行 (要传给 Mesh()), ensemble.ensemble_comm 用于组间通信, 还提供了可直接作用于 Function 的 allreduce、send/recv、bcast 等方法. 具体用法见官方文档 Parallelism in Firedrake 和 Ensemble parallelism.

1.5. 并行常见陷阱#

本节汇总几个初学者最容易踩的并行陷阱. 下面的演示单元继续使用前面启动的 2 引擎集群 (以及其中已定义的 mesh 等变量).

1.5.1. 集合操作必须所有进程一起调用#

这是并行程序中最典型也最隐蔽的错误: 涉及 MPI 集合通信的操作, 一旦只让部分进程执行, 参与的进程就会永远等待缺席的进程, 程序挂起且没有任何报错. 例如, 想”只画一次网格”而写出下面的代码:

if COMM_WORLD.rank == 0:
    fig, axes = plt.subplots()
    triplot(mesh, axes=axes)   # 危险: 程序会永久挂起!

用 mpiexec -n 2 运行时, 0 号进程进入 triplot 后需要读取带 halo 的坐标数据, 于是等待与其他进程同步; 而 1 号进程早已跳过该分支继续前进, 双方永远等不到对方. 即使不挂起, 这个写法在语义上也是错的: 0 号进程只保存自己那部分子网格, 画出来的不是整张网格. 想看整张网格, 应当像本章前面那样让所有进程一起调用 triplot (各画各的部分), 或者用 VTKFile 输出后在 ParaView 中查看.

下面这些操作都属于集合操作 (均在 2 进程下实测过, 只在一个进程上调用就会挂起):

assemble、solve、norm、errornorm 等涉及全局归约或求解的计算;
triplot 等绘图函数, VTKFile(...).write(...) 等并行输出;
点求值 (见下文) 以及 comm.allreduce 等显式通信调用.

更难排查的是只在特定状态下才通信的属性访问: 例如通过可写的 dat.data 修改数据后, halo 被标记为失效, 此后第一次访问 dat.data_ro_with_halos 会触发同步通信, 只在单个进程上执行同样会挂起 (halo 有效时访问则不通信、不会挂起). 因此稳妥的原则是: 除非确定操作是纯本地的 (普通 print、写文本文件、处理已经算好的数值), 否则让所有进程一起执行. 正确的写法如前文所示: 集合操作在条件分支外由所有进程一起完成, 分支内只做纯本地的输出.

1.5.2. 局部数据上的归约不是全局结果#

dat.data_ro 只包含本进程拥有的那部分自由度, 直接在上面调用 max、min、sum 得到的只是局部结果, 各进程互不相同. 要得到全局结果, 可以进入 PETSc Vec 视图 (vec.max()、vec.norm() 等都是集合操作, 返回全局值), 或者用 mesh.comm.allreduce 手动归约; assemble、norm、errornorm 这类函数返回的本来就是全局值.

并行 (2 个引擎)

%%px --block
x, y = SpatialCoordinate(mesh)
g = Function(FunctionSpace(mesh, 'CG', 1)).interpolate(x + y)
local_max = g.dat.data_ro.max()
with g.dat.vec_ro as vec:
    _, global_max = vec.max()
PETSc.Sys.syncPrint(f'[{rank}/{size}] local max: {local_max:.4f}, global max: {global_max:.4f}')
PETSc.Sys.syncFlush()

[stdout:0] [0/2] local max: 2.0000, global max: 2.0000
[1/2] local max: 1.5000, global max: 2.0000

函数 g = x + y 的全局最大值在角点 (1, 1) 处取得, 该点只属于某一个进程: 另一个进程的局部最大值达不到 2, 而 vec.max() 在所有进程上都返回同一个全局最大值.

1.5.3. 各进程的随机数各自独立#

每个引擎/进程都是独立的 Python 进程, numpy 的随机数生成器互不同步: 不固定种子时, 各进程得到的随机数完全不同.

并行 (2 个引擎)

%%px --block
import numpy as np
PETSc.Sys.syncPrint(f'[{rank}/{size}] {np.random.default_rng().random(3)}')
PETSc.Sys.syncFlush()

[stdout:0] [0/2] [0.61369732 0.29765761 0.71182686]
[1/2] [0.61218728 0.01528603 0.85504656]

如果程序逻辑依赖”所有进程看到同一个随机数” (例如随机选择扰动幅度), 必须固定种子或由 0 号进程生成后广播, 否则各进程的状态会悄悄发散. 反过来, 即便固定了种子, 用随机数给 Function 的自由度赋值时, 每个进程也只是用相同的随机数流填充自己拥有的那部分自由度, 得到的全局场与串行运行并不相同, 且依赖于进程数和分区方式 (Cahn–Hilliard 方程一章中的随机初值算例即属此类).

1.5.4. 点求值是集合操作#

查询函数在某个点的值时, 该点通常只落在某一个进程的子网格中, 但点求值 (PointEvaluator) 是集合操作: 必须所有进程一起调用, 求得的值会同步到所有进程; 只在单个进程上调用会像前面一样挂起.

并行 (2 个引擎)

%%px --block
point_eval = PointEvaluator(mesh, [[0.5, 0.5]])
PETSc.Sys.syncPrint(f'[{rank}/{size}] g(0.5, 0.5) = {point_eval.evaluate(g)}')
PETSc.Sys.syncFlush()

[stdout:0] [0/2] g(0.5, 0.5) = [1.]
[1/2] g(0.5, 0.5) = [1.]

1.5.5. `from firedrake import *` 会覆盖 `rank` 等变量#

这一条与 MPI 无关, 却几乎只在并行代码里发作: from firedrake import * 导出的名字中, 恰好有几个是并行程序最爱用的变量名. 下面看其中三个 (外加一个没有被导出的):

并行 (2 个引擎)

%%px --block
import firedrake
for name in ('mesh', 'halo', 'rank', 'size'):
    obj = getattr(firedrake, name, None)
    kind = 'not exported' if obj is None else type(obj).__name__
    PETSc.Sys.Print(f'firedrake.{name:<5} -> {kind}')

[stdout:0] firedrake.mesh  -> module
firedrake.halo  -> module
firedrake.rank  -> function
firedrake.size  -> not exported

mesh 和 halo 是 Firedrake 的两个子模块, rank 是 UFL 的一个函数 (返回张量表达式的阶), 而 size 根本没有被导出. 于是, 在写下 rank, size = mesh.comm.rank, mesh.comm.size 之后, 只要再执行一次 from firedrake import *, rank 就会被换成那个函数, size 却原样保留; 此后 f'[{rank}/{size}]' 打印出来的会是 <function rank at 0x...>/2 这种一半正常、一半错乱的东西, 而且全程没有任何报错. 本章”分区重叠 (overlap) 的显式控制”一节的单元, 正是因为它前面的 %%px 单元二次执行了 from firedrake import *, 才在开头重新取了一次 rank 和 size.

变量名 mesh 同理: 若某个单元先 from firedrake import *, 再使用上一个单元定义的 mesh, 拿到的会是 firedrake.mesh 模块而不是网格. 稳妥的做法是把导入集中在最开始的单元里, 不要在后续单元中反复 import *; 若确实需要重复导入, 就在导入之后重新给这些变量赋值.

1.5.6. 迭代法的迭代次数随进程数变化#

不少常用的预条件子 (如 block Jacobi、加性 Schwarz) 是按进程分块构造的: 进程数变了, 预条件子本身就变了. 因此同一个问题用不同进程数求解, Krylov 迭代次数一般不同, 这并不是程序出错. 下面用 CG 加 block Jacobi (子块用 ILU) 求解同一个 Poisson 问题:

串行 (本地内核)

mesh_iter = UnitSquareMesh(32, 32)
V_iter = FunctionSpace(mesh_iter, 'CG', 1)
u, v = TrialFunction(V_iter), TestFunction(V_iter)
x, y = SpatialCoordinate(mesh_iter)
a_iter = inner(grad(u), grad(v))*dx
L_iter = inner(sin(pi*x)*sin(pi*y), v)*dx
uh_iter = Function(V_iter)
solver_iter = LinearVariationalSolver(
    LinearVariationalProblem(a_iter, L_iter, uh_iter, bcs=DirichletBC(V_iter, 0, 'on_boundary')),
    solver_parameters={'ksp_type': 'cg', 'pc_type': 'bjacobi', 'sub_pc_type': 'ilu'})
solver_iter.solve()
PETSc.Sys.Print(f'iterations: {solver_iter.snes.ksp.getIterationNumber()}, '
                f'integral of u_h: {assemble(uh_iter*dx):.10e}')

iterations: 15, integral of u_h: 2.0482533688e-02

并行 (2 个引擎)

[stdout:0] iterations: 24, integral of u_h: 2.0482533691e-02

两次求解的迭代次数不同, 但解在求解容差内是一致的: 两个积分值的前若干位数字完全相同. 使用直接法时 (Firedrake solve 未指定求解器参数时的默认行为) 没有迭代次数的问题, 串行与并行结果的差异在机器精度量级. 需要检查并行结果是否正确时, 就可以像这样比较 assemble、norm、errornorm 等集合运算的返回值: 它们在所有进程上返回同一个全局值, 可以直接与串行结果对照.

1.6. 如何验证并行结果是对的#

前面几节一直在观察网格和自由度是怎么被切分的, 但真正需要回答的问题是: 同一份程序改用多个进程运行之后, 算出来的结果还对不对?

判断的办法是比较集合量. assemble、norm、errornorm 这类函数内部已经做过全局归约, 在所有进程上返回同一个数, 因此可以用 PETSc.Sys.Print 打印一次, 直接与串行运行的结果逐位对照. 反过来, 像 dat.data_ro.max() 这样的局部量在各进程上互不相同, 不能拿来做这种对比 (见并行常见陷阱).

下面在单位正方形上求解 Poisson 方程

\[ -\Delta u = 2\pi^2 \sin(\pi x)\sin(\pi y), \qquad u|_{\partial\Omega} = 0, \]

其精确解为 \(u = \sin(\pi x)\sin(\pi y)\). 求解器显式指定为 preonly + lu, 即由 MUMPS 完成 LU 分解的直接法; 串行单元和并行单元运行的是完全相同的代码.

串行 (本地内核)

mesh_cmp = UnitSquareMesh(16, 16, name='poisson_mesh')
V_cmp = FunctionSpace(mesh_cmp, 'CG', 1)
u, v = TrialFunction(V_cmp), TestFunction(V_cmp)
x, y = SpatialCoordinate(mesh_cmp)
u_exact = sin(pi*x)*sin(pi*y)
uh_cmp = Function(V_cmp, name='u_h')
solve(inner(grad(u), grad(v))*dx == inner(2*pi**2*u_exact, v)*dx, uh_cmp,
      bcs=DirichletBC(V_cmp, 0, 'on_boundary'),
      solver_parameters={'ksp_type': 'preonly', 'pc_type': 'lu',
                         'pc_factor_mat_solver_type': 'mumps'})
int_cmp = assemble(uh_cmp*dx)
norm_cmp = norm(uh_cmp)
err_cmp = errornorm(u_exact, uh_cmp)
PETSc.Sys.Print(f'processes : {mesh_cmp.comm.size}')
PETSc.Sys.Print(f'integral  : {int_cmp:.16e}')
PETSc.Sys.Print(f'L2 norm   : {norm_cmp:.16e}')
PETSc.Sys.Print(f'errornorm : {err_cmp:.16e}')

processes : 1
integral  : 4.0139184844570802e-01
L2 norm   : 4.9521125786396453e-01
errornorm : 5.3774350099476901e-03

tsfc:WARNING Estimated quadrature degree 12 more than tenfold greater than any argument/coefficient degree (max 1)

并行 (2 个引擎)

[stdout:0] processes : 2
integral  : 4.0139184844570841e-01
L2 norm   : 4.9521125786396508e-01
errornorm : 5.3774350099472130e-03

两组输出已经可以逐位比较了. 为了把差异算成具体数字, 下面在本地内核里用 client[0].pull 把 0 号引擎中的三个变量取回来: 引擎和本地内核各有自己的一套变量, 必须显式取回; 而这三个量都是集合量, 从哪个引擎取回都一样.

int_par, norm_par, err_par = client[0].pull(
    ['int_cmp', 'norm_cmp', 'err_cmp'], block=True)
for name, serial_value, parallel_value in (('integral', int_cmp, int_par),
                                           ('L2 norm', norm_cmp, norm_par),
                                           ('errornorm', err_cmp, err_par)):
    diff = abs(serial_value - parallel_value)
    PETSc.Sys.Print(f'{name:>9}: absolute difference {diff:.2e}, '
                    f'relative difference {diff/abs(serial_value):.2e}')

 integral: absolute difference 3.89e-16, relative difference 9.68e-16
  L2 norm: absolute difference 5.55e-16, relative difference 1.12e-15
errornorm: absolute difference 4.77e-16, relative difference 8.87e-14

integral 和 L2 norm 的相对差在 \(10^{-16}\)–\(10^{-15}\) 量级, 也就是双精度机器精度 (约 \(2.2\times 10^{-16}\)) 的几倍到十几倍, 可以认为串行和并行求得的是同一个解. 这点差异不是算法不同造成的, 而是因为浮点加法不满足结合律: 并行时矩阵组装的求和次序、MUMPS 分解中的运算次序都与串行不同, 末尾一两位对不上是正常现象; 换用别的进程数 (例如 3 个进程), 末位数字还会再变一次.

errornorm 的相对差要大得多 (\(10^{-13}\) 量级), 但它的绝对差与前两个量属于同一量级 (三者都在 \(10^{-16}\)–\(10^{-15}\) 之间). 原因是 errornorm 本身只有 \(5\times 10^{-3}\) 左右: 由三角不等式, 两个 errornorm 之差不会超过两个解之差的范数, 所以绝对差仍停留在机器精度量级; 只是 errornorm (约 \(5\times 10^{-3}\)) 比解本身的范数 (约 \(0.5\)) 小了两个数量级, 同样大小的绝对差除以它, 相对差自然会被放大约一百倍. 因此判断串并是否一致, 应当看绝对差异并结合被比较量自身的尺度, 而不是只盯着相对差的位数.

以上是直接法的结论. 换成迭代法, 串行与并行只在求解容差内一致, 而且迭代次数本身就会随进程数改变 (预条件子按进程分块构造), 这在并行常见陷阱的”迭代法的迭代次数随进程数变化”一节中已经演示过. 注意 Firedrake 的 solve 在不指定求解器参数时用的就是直接法, 想观察迭代法的行为必须显式设置 ksp_type 和 pc_type; 用迭代法比较串并结果时, 一致性也只能要求到 ksp_rtol 的量级, 不能要求机器精度.

1.7. 并行下的 I/O#

并行程序写出的文件有两类, 用途完全不同, 不能混用.

第一类是可视化输出. 本章开头用 mpiexec 运行 poisson.py 时已经见过: VTKFile('data/result.pvd') 在并行下会按进程拆成多个文件, 每个文件只含该进程负责的那块子区域, .pvd 只是这些文件的索引. VTKFile 写出的是专为可视化准备的数据, Firedrake 也没有提供把 .pvd/.vtu 读回成 Function 的接口 (VTKFile 只有 write 一个方法). 因此它只能用来在 ParaView 等工具里查看结果, 不能用来续算.

第二类是用于续算的检查点, 应当使用 CheckpointFile. 它基于 HDF5, 保存的是网格拓扑和自由度数据本身, 并且写出的文件与进程数无关: 用 N 个进程写出的文件可以用 M 个进程读回, M 不必等于 N. 下面把上一节并行求得的解存下来, 再用本地内核 (单进程) 读回.

并行 (2 个引擎)

%%px --block
import os
os.makedirs('data', exist_ok=True)
with CheckpointFile('data/parallel_checkpoint.h5', 'w') as chk:
    chk.save_mesh(mesh_cmp)
    chk.save_function(uh_cmp)
PETSc.Sys.Print(f'saved with {mesh_cmp.comm.size} processes')

[stdout:0] saved with 2 processes

串行 (本地内核)

with CheckpointFile('data/parallel_checkpoint.h5', 'r') as chk:
    mesh_load = chk.load_mesh('poisson_mesh')
    uh_load = chk.load_function(mesh_load, 'u_h')

x, y = SpatialCoordinate(mesh_load)
norm_load = norm(uh_load)
PETSc.Sys.Print(f'loaded with {mesh_load.comm.size} process(es), '
                f'cells: {mesh_load.num_cells()}')
PETSc.Sys.Print(f'integral  : {assemble(uh_load*dx):.16e}')
PETSc.Sys.Print(f'L2 norm   : {norm_load:.16e}')
PETSc.Sys.Print(f'errornorm : {errornorm(sin(pi*x)*sin(pi*y), uh_load):.16e}')
PETSc.Sys.Print(f'difference from the parallel L2 norm: '
                f'{abs(norm_load - norm_par):.2e}')

loaded with 1 process(es), cells: 512
integral  : 4.0139184844570847e-01
L2 norm   : 4.9521125786396508e-01
errornorm : 5.3774350099472165e-03
difference from the parallel L2 norm: 0.00e+00

文件由 2 个进程写出, 却在单进程的本地内核里读回成功: 读回后整张网格都落在一个进程上 (cells 是全部单元数, 不再是某个子区域), 分区方式与写出时完全不同, 而三个集合量与上面并行单元打印的数值一致到机器精度. 也就是说, 存进去、取出来的确实是并行算出的那个解; 它与串行求解结果之间的差异, 仍是上一节讨论的机器精度量级.

使用 CheckpointFile 时有几点需要注意:

上面显式调用了 save_mesh 并给网格起了名字 (poisson_mesh), 这是推荐做法而不是硬性要求: 实测只调 save_function 也能读回 (它会把网格一并存进去), 网格不起名时默认叫 firedrake_default, load_mesh() 不传参数也取得到. 但一个文件里存多张网格、多个函数时就只能靠名字定位了, 因此建议养成显式保存并命名的习惯. 函数同样靠 Function 的 name (这里是 u_h) 定位.
读回时先 load_mesh 得到网格, 再把它传给 load_function. 读回的网格是一个新的网格对象, 与串行单元里的 mesh_cmp 不是同一张网格, 因此不能直接对两者的函数做 errornorm (会因为”多个积分区域”而报错), 只能像上面那样比较各自算出的标量.
CheckpointFile 的读写都是集合操作, 必须所有进程一起进入 with 块, 不能包在 if rank == 0 里.
检查点文件通常不小, 应写到不纳入版本库的目录. 上面写在 data/ 下 (本仓库已将该目录排除在版本库之外), 并用 os.makedirs(..., exist_ok=True) 保证目录存在; 打开模式 'w' 会覆盖同名文件, 因此这两个单元可以反复执行.

1.7.1. 只在 0 号进程写日志#

除了解场之外, 并行程序还常常要记录收敛历史、能量、时间步长这类标量日志. 这类文件是普通文本, 若每个进程都写一遍, 轻则内容重复, 重则互相覆盖, 因此应当只让 0 号进程写. 关键是把集合操作和写文件分成两步:

energy = assemble(0.5*inner(grad(u_h), grad(u_h))*dx)  # 集合操作: 所有进程一起调用
if mesh.comm.rank == 0:                                # 纯 Python I/O: 只在 0 号进程执行
    with open('data/history.log', 'a') as log:
        log.write(f'{step} {energy}\n')

只要条件分支里剩下的是纯 Python 的文件操作, 这样写就是安全的; 一旦把 assemble、norm 之类的集合操作挪进分支, 程序就会永久挂起. 详见前文和并行常见陷阱.

1.8. 并行性能与加速比#

最后回到最实际的问题: 多用几个进程, 到底能快多少?

下面固定问题类型 (单位正方形上的 Poisson 方程, \(P_1\) 单元), 在几个不同的网格规模上各求解一次, 计时区间是一次 solver.solve() 调用. 它到底包含什么值得说清楚: LinearVariationalProblem 默认 constant_jacobian=False, 因此每次 solve() 都会重新组装矩阵和右端项, 这部分开销是计在时间里的; 网格和函数空间的建立不在计时区间内, 第一次调用时的即时编译也已被预热排除 —— 每个规模先空跑一次, 再重复若干次取最短的一次 (取最小值可以削弱机器上其他任务造成的干扰).

求解器选用 CG 加点 Jacobi 预条件子, 这个选择是刻意的: 点 Jacobi 只用到矩阵的对角线, 与网格如何分区无关, 因此串行和并行的迭代次数完全相同, 时间才具有可比性. 若换成前面用过的 block Jacobi, 迭代次数本身就随进程数变化, 时间差里就混进了算法的因素.

串行 (本地内核)

import time

def time_solve(n, repeat=10):
    mesh_t = UnitSquareMesh(n, n)
    V_t = FunctionSpace(mesh_t, 'CG', 1)
    u, v = TrialFunction(V_t), TestFunction(V_t)
    x, y = SpatialCoordinate(mesh_t)
    uh_t = Function(V_t)
    solver = LinearVariationalSolver(
        LinearVariationalProblem(inner(grad(u), grad(v))*dx,
                                 inner(sin(pi*x)*sin(pi*y), v)*dx, uh_t,
                                 bcs=DirichletBC(V_t, 0, 'on_boundary')),
        solver_parameters={'ksp_type': 'cg', 'pc_type': 'jacobi',
                           'ksp_rtol': 1e-8})
    solver.solve()                        # 预热: 先触发即时编译
    best = float('inf')
    for _ in range(repeat):
        uh_t.assign(0)
        mesh_t.comm.Barrier()             # 让各进程同时开始计时
        t0 = time.perf_counter()
        solver.solve()
        mesh_t.comm.Barrier()             # 等最慢的进程算完
        best = min(best, time.perf_counter() - t0)
    return V_t.dim(), solver.snes.ksp.getIterationNumber(), best

mesh_sizes = [4, 16, 64, 256]
timings = {n: time_solve(n) for n in mesh_sizes}   # 先全部算完, 最后统一输出

PETSc.Sys.Print(f'{"dofs":>8} {"iterations":>11} {"solve time (s)":>15}')
for n in mesh_sizes:
    PETSc.Sys.Print(f'{timings[n][0]:>8} {timings[n][1]:>11} '
                    f'{timings[n][2]:>15.6f}')

    dofs  iterations  solve time (s)
         2        0.000391
        21        0.000558
        84        0.004998
       315        0.176780

并行 (2 个引擎): 这个单元运行时间较长, 因此加上 --progress-after -1 关闭 %%px 默认在 2 秒后显示的进度条.

[stdout:0]     dofs  iterations  solve time (s)
         2        0.000497
        21        0.000661
        84        0.003440
       315        0.094416

两边的迭代次数逐行相同, 说明求解的是同一个问题、走的是同一条迭代路径. 下面把并行的计时取回本地内核, 算出加速比 (串行时间 ÷ 并行时间):

timings_par = client[0].pull('timings', block=True)
n_engines = len(client)
label = f'{n_engines} processes (s)'
PETSc.Sys.Print(f'{"dofs":>8} {"dofs/process":>13} {"serial (s)":>12} '
                f'{label:>17} {"speedup":>8}')
for n in mesh_sizes:
    dofs, _, t_serial = timings[n]
    t_parallel = timings_par[n][2]
    PETSc.Sys.Print(f'{dofs:>8} {dofs//n_engines:>13} {t_serial:>12.6f} '
                    f'{t_parallel:>17.6f} {t_serial/t_parallel:>8.2f}')

    dofs  dofs/process   serial (s)   2 processes (s)  speedup
          12     0.000391          0.000497     0.79
         144     0.000558          0.000661     0.84
        2112     0.004998          0.003440     1.45
       33024     0.176780          0.094416     1.87

表里的趋势比其中任何一个数字都重要:

规模最小的两行, 并行不但没有收益, 通常还略慢. 几十到几百个自由度分给 2 个进程, 每个进程只算很少的数, 而每步 CG 迭代都要多做全局内积 (MPI 归约) 和 halo 交换 —— 这些是净增的开销, 并没有换来相应的计算量减少, 所以变慢是系统性的. 至于具体慢多少、某一次运行是否恰好越过 1, 则取决于测量抖动, 加速比在 1 上下浮动都属正常. 本章前面一直在用的 4×4 网格正属于这一档 —— 它是用来讲清楚分区和自由度分布的教学网格, 不是用来加速的.
随着规模增大, 加速比明显上升并越过 1. 每个进程分到的计算量足够大之后, 通信开销才被摊薄. 本例中的转折点大致落在每进程几百到几千个自由度之间.
由此可以记住一条经验: 每进程的自由度太少就不值得再往下分. 这里要区分两件事: “加速比越过 1”只要求并行不亏本, 每进程几百到几千个自由度就可能做到; 而”接近线性的并行效率” (2 个进程接近 2 倍、4 个进程接近 4 倍) 门槛要高得多, 常见的经验值是每进程至少要有一两万个自由度. 一味增加进程数, 只会让通信开销吃掉全部收益.

还要提醒的是, 这类计时依赖于运行它的机器: 核数、内存带宽, 以及当时机器上还跑着什么, 都会影响结果. 你在网页上看到的这张表是构建这本笔记时算出来的, 具体数值随构建机器而变, 可能高于也可能低于你在自己机器上得到的结果. 请把注意力放在趋势上, 不要记住某一次的具体秒数. 另外, 这里的 2 个引擎和 notebook 内核同在一台机器上、彼此争抢资源, 严格的性能测量应当在空闲的机器上用 mpiexec 运行独立脚本.

最后区分两个容易混用的概念:

强可扩展性 (strong scaling): 固定问题总规模, 增加进程数, 看时间能缩短多少. 上表每一行的串并对比就是它的最小版本; 进程数继续增加时, 每进程的自由度会越来越少, 迟早退回到”通信占主导”的那一档.
弱可扩展性 (weak scaling): 让每进程的自由度保持不变, 进程数和问题总规模同步放大, 看时间能否基本持平.

想进一步弄清时间到底花在哪里 (组装、求解、通信各占多少), 可以用 PETSc 的 -log_view, 详见性能分析一章.

1.9. 关闭集群#

引擎是 notebook 之外的独立进程, 关掉 notebook 并不会自动把它们收走. 用完之后应当显式停止集群, 否则这些进程会一直占着内存和 CPU; 在自己机器上反复重跑本章时, 还会不断累积出新的引擎.

下面的单元停掉控制器和全部引擎. 它必须是全章最后执行的单元 —— 在它之后再出现 %%px 或用到 client 的代码都会失败.

cluster.stop_cluster_sync()

Show code cell output

Hide code cell output

Stopping controller
Controller stopped: {'exit_code': 0, 'pid': 2700, 'identifier': 'ipcontroller-1785587441-4oi7-2681'}
Stopping engine(s): 1785587442
Output for ipengine-1785587441-4oi7-1785587442-2681:
2026-08-01 13:30:43.580 [KernelNanny.1] Nanny watching parent pid 2736.
2026-08-01 13:30:43.644 [IPEngine.0.0] Loading IPython extension: storemagic
2026-08-01 13:30:43.645 [IPEngine.0.0] Running code in user namespace: 
from mpi4py import MPI
mpi_rank = MPI.COMM_WORLD.Get_rank()
mpi_size = MPI.COMM_WORLD.Get_size()

2026-08-01 13:30:43.646 [IPEngine.0.0] WARNING | debugpy_stream undefined, debugging will not be enabled
2026-08-01 13:30:43.649 [IPEngine.0.0] Starting to monitor the heartbeat signal from the hub every 3500 ms.
2026-08-01 13:30:43.649 [IPEngine.0.0] Completed registration with id 0
2026-08-01 13:30:43.650 [IPEngine.1.1] Loading IPython extension: storemagic
2026-08-01 13:30:43.651 [IPEngine.1.1] Running code in user namespace: 
from mpi4py import MPI
mpi_rank = MPI.COMM_WORLD.Get_rank()
mpi_size = MPI.COMM_WORLD.Get_size()

2026-08-01 13:30:43.653 [IPEngine.1.1] WARNING | debugpy_stream undefined, debugging will not be enabled
2026-08-01 13:30:43.654 [IPEngine.1.1] Starting to monitor the heartbeat signal from the hub every 3500 ms.
2026-08-01 13:30:43.655 [IPEngine.1.1] Completed registration with id 1
2026-08-01 13:30:48.144 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_2
2026-08-01 13:30:48.145 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_1
2026-08-01 13:30:49.744 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_3
2026-08-01 13:30:49.755 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_4
2026-08-01 13:30:49.757 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_5
2026-08-01 13:30:52.811 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_6
2026-08-01 13:30:52.812 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_7
2026-08-01 13:30:52.930 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_9
2026-08-01 13:30:52.930 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_8
2026-08-01 13:30:52.978 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_10
2026-08-01 13:30:52.978 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_11
2026-08-01 13:30:53.134 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_12
2026-08-01 13:30:53.135 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_13
2026-08-01 13:30:53.151 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_14
2026-08-01 13:30:53.151 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_15
2026-08-01 13:30:53.169 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_16
2026-08-01 13:30:53.169 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_17
2026-08-01 13:30:53.184 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_18
2026-08-01 13:30:53.185 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_19
2026-08-01 13:30:53.200 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_20
2026-08-01 13:30:53.201 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_21
2026-08-01 13:30:53.217 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_22
2026-08-01 13:30:53.218 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_23
2026-08-01 13:30:53.234 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_24
2026-08-01 13:30:53.235 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_25
2026-08-01 13:30:53.251 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_26
2026-08-01 13:30:53.253 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_27
2026-08-01 13:30:53.270 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_28
2026-08-01 13:30:53.270 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_29
2026-08-01 13:30:53.285 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_30
2026-08-01 13:30:53.288 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_31
2026-08-01 13:30:53.305 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_32
2026-08-01 13:30:53.305 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_33
2026-08-01 13:30:53.317 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_35
2026-08-01 13:30:53.317 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_34
2026-08-01 13:30:54.396 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_36
2026-08-01 13:30:54.397 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_37
2026-08-01 13:30:54.707 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_38
2026-08-01 13:30:54.708 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_39
2026-08-01 13:30:54.719 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_40
2026-08-01 13:30:54.719 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_41
2026-08-01 13:30:55.223 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_42
2026-08-01 13:30:55.224 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_43
2026-08-01 13:30:55.868 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_45
2026-08-01 13:30:55.868 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_44
2026-08-01 13:30:57.238 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_46
2026-08-01 13:30:57.239 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_47
2026-08-01 13:30:57.334 [IPEngine.0.0] Handling apply_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_48
2026-08-01 13:30:57.343 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_49
2026-08-01 13:30:57.346 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_50
2026-08-01 13:31:03.684 [IPEngine.0.0] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_51
2026-08-01 13:31:03.685 [IPEngine.1.1] Handling execute_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_52
2026-08-01 13:31:07.747 [IPEngine.0.0] Handling apply_request: 87a18dd0-4cc0546de09b1c3cc10f71d4_2681_53
2026-08-01 13:31:08.150 [IPEngine.0.0] WARNING | No heartbeat in the last 3500 ms (1 time(s) in a row).
2026-08-01 13:31:08.869 [IPEngine.0.0] CRITICAL | received signal 15, stopping
Traceback (most recent call last):
  File "<frozen runpy>", line 198, in _run_module_as_main
  File "<frozen runpy>", line 88, in _run_code
  File "/usr/local/lib/python3.12/dist-packages/ipyparallel/engine/__main__.py", line 4, in <module>
    main()
  File "/usr/local/lib/python3.12/dist-packages/traitlets/config/application.py", line 1082, in launch_instance
    app.start()
  File "/usr/local/lib/python3.12/dist-packages/ipyparallel/engine/app.py", line 995, in start
    self.loop.run_sync(self._start)
  File "/usr/local/lib/python3.12/dist-packages/tornado/ioloop.py", line 546, in run_sync
    raise RuntimeError("Event loop stopped before Future completed.")
RuntimeError: Event loop stopped before Future completed.
2026-08-01 13:31:08.156 [IPEngine.1.1] WARNING | No heartbeat in the last 3500 ms (1 time(s) in a row).
2026-08-01 13:31:08.869 [IPEngine.1.1] CRITICAL | received signal 15, stopping
Traceback (most recent call last):
  File "<frozen runpy>", line 198, in _run_module_as_main
  File "<frozen runpy>", line 88, in _run_code
  File "/usr/local/lib/python3.12/dist-packages/ipyparallel/engine/__main__.py", line 4, in <module>
    main()
  File "/usr/local/lib/python3.12/dist-packages/traitlets/config/application.py", line 1082, in launch_instance
    app.start()
  File "/usr/local/lib/python3.12/dist-packages/ipyparallel/engine/app.py", line 995, in start
    self.loop.run_sync(self._start)
  File "/usr/local/lib/python3.12/dist-packages/tornado/ioloop.py", line 546, in run_sync
    raise RuntimeError("Event loop stopped before Future completed.")
RuntimeError: Event loop stopped before Future completed.
engine set stopped 1785587442: {'exit_code': 1, 'pid': 2728, 'identifier': 'ipengine-1785587441-4oi7-1785587442-2681'}

输出是控制器和各引擎的停止日志 (含进程号等信息, 默认已折叠). 如果集群是在终端里用 ipcluster start 启动的, 则用 ipcluster stop 停止.

并行计算 Parallel Computing

Contents

1. 并行计算 Parallel Computing#

1.1. Firedrake 中并行计算简介#

1.1.1. 并行启动程序#

1.1.2. 输出与进程控制#

1.1.2.1. 并行时的输出#

1.1.2.2. 仅在指定进程上运算#

1.2. 在 Jupyter 中执行并行代码#

1.2.1. ipyparallel 简介#

1.2.2. 执行并行代码#

1.2.3. %%px 的更多用法#

1.3. 使用 ipyparallel 观察串行和并行过程#

1.3.1. 生成网格并画出网格#

1.3.2. 分区重叠 (overlap) 的显式控制#

1.3.3. 定义变分问题#

1.3.4. 函数空间维度#

1.3.5. 节点集和自由度数据集#

1.3.5.1. Set 与 DataSet 的各种 size#

1.3.5.2. DataSet 的索引集 (IS)#

1.3.5.3. 局部到全局的映射#

1.3.5.4. DataSet 的布局向量#

1.3.6. 矩阵组装#

1.4. 通信子: COMM_WORLD 与 mesh.comm#

1.5. 并行常见陷阱#

1.5.1. 集合操作必须所有进程一起调用#

1.5.2. 局部数据上的归约不是全局结果#

1.5.3. 各进程的随机数各自独立#

1.5.4. 点求值是集合操作#

1.5.5. from firedrake import * 会覆盖 rank 等变量#

1.5.6. 迭代法的迭代次数随进程数变化#

1.6. 如何验证并行结果是对的#

1.7. 并行下的 I/O#

1.7.1. 只在 0 号进程写日志#

1.8. 并行性能与加速比#

1.9. 关闭集群#

1.2.1. `ipyparallel` 简介#

1.2.3. `%%px` 的更多用法#

1.3. 使用 `ipyparallel` 观察串行和并行过程#

1.4. 通信子: `COMM_WORLD` 与 `mesh.comm`#

1.5.5. `from firedrake import *` 会覆盖 `rank` 等变量#