5. PETSc#
PETSc, the Portable, Extensible Toolkit for Scientific Computation, pronounced PET-see (/ˈpɛt-siː/), is for the scalable (parallel) solution of scientific applications modeled by partial differential equations.
Overview
A tour of the PETSc objects used by Firedrake (via petsc4py): vectors and matrices, KSP solvers, viewers, the star forest PetscSF, and the PCPatch preconditioner.
PETSc 是由阿贡国家实验室开发的便携可扩展科学计算工具包, 提供了并行求解大规模方程组的许多算法, 并且可调用外部包对方程组进行求解. 另外, PETSc 也提供了用于数值求解偏微分方程的组件, 包括结构化网格数据结构, 非结构化网格数据结构, 和有限元空间等.
Firedrake 可以看作是 PETSc 的高层次封装. 在 Firedrake 中, 一般无需操作 PETSc 对象, 但有些特殊情况必须直接操纵 PETSc 对象, 并且有时候直接操纵 PETSc 对象会更高效.
PETSc 是 c 语言包, 也提供了 Fortran 接口. petsc4py 是 PETSc 的 python 封装.
学习 PETSc 可以从 Texas Advanced Computing Center (TACC) 发布的 PETSc 入门课程开始
PETSc 网站的手册和入门讲义
C/Fortran API: https://petsc.org/release/manualpages/
petsc4py: https://petsc.org/release/petsc4py/
PETSc 代码库有许多示例可以作为学习材料, 如 PETSc 仓库中 petsc4py 的示例:
Firedrake 已经把 PETSc 包好了, 日常使用时通过 solver_parameters 字典设置求解器就够了. 需要往下一层直接操作 PETSc 对象的, 通常是下面几种情形:
需要矩阵或向量本身: 例如把组装好的矩阵取出来做谱分析、算行列式, 或者按 PETSc 的格式存盘交给别的程序处理;
需要弄清求解器为什么失败: Firedrake 只会告诉你没有收敛, 具体的收敛状态和预条件子的失败原因要从
KSP、PC对象上读;需要操作网格拓扑: 网格分区、粗化、编号以及并行的数据交换都发生在
DMPlex和PetscSF这一层, Firedrake 的Mesh只是它们的包装;需要 Firedrake 尚未包装的功能: 例如
Viewer的各种输出格式, 或者PCPatch这类要求用户自行提供拓扑信息的预条件子.
本章各节就是按这几类需求组织的, 叙述上彼此独立, 可以按需查阅:
Vector和Matrix: 取出底层的Mat、Vec, 以及与numpy/scipy之间的转换;Options: PETSc 的选项数据库, 是所有 PETSc 参数的统一入口;KSP: 线性求解器对象, 重点是求解失败之后如何取得原因;DMPlex、Viewer、PetscSF: 网格拓扑、对象输出与进程间的数据交换;PCPatch: 单元片预条件子的内部数据结构, 是最深的一层.
不过叙述独立不等于代码独立: 后面带 %%px 的单元都跑在 Viewer 一节里启动的 ipyparallel 集群上, 想单独执行它们, 要先运行启动集群的那一格.
初次阅读建议先看 Options 和 KSP 两节, 它们对应日常的调参和排错; DMPlex 之后的内容涉及并行数据结构, 可以等到真正需要时再回头看.
Tip
PETSc 目录中有用的工具, 如 h5dump, petsc_gen_xdmf.py, PetscBinaryIO.py 等.
在 PETSc 环境中, 运行如下命令添加这些工具所在路径到 PATH:
export PATH="$PATH:$PETSC_DIR/lib/petsc/bin"
export PATH="$PATH:$PETSC_DIR/${PETSC_ARCH-default}/bin"
在激活的 Firedrake 环境可以运行如下命令的输出, 添加工具所在路径到环境变量 PATH.
python -c "from firedrake import *; \
import os; \
PETSC_DIR = os.environ['PETSC_DIR']; \
PETSC_ARCH = os.environ['PETSC_ARCH']; \
print('\nRun the follwoing code to add petsc/bin to path:\n'); \
print(f' export PATH=\"\$PATH:{PETSC_DIR}/lib/petsc/bin\"'); \
print(f' export PATH=\"\$PATH:{PETSC_DIR}/{PETSC_ARCH}/bin\"'); \
print('');"
5.1. Vector 和 Matrix#
Firedrake 组装出来的矩阵和右端项都是 PETSc 对象的包装. 本节说明如何取出底层的 Mat 和 Vec, 如何把矩阵转成 numpy/scipy 的数据结构, 以及如何把矩阵存到文件. 下面用一个 4×4 的矩形网格上的 Poisson 问题作为例子.
PETSc 官方的矩阵、向量读写示例见 matvecio.py.
from firedrake import *
from firedrake.petsc import PETSc
test_mesh = RectangleMesh(nx=4, ny=4, Lx=1, Ly=1)
x, y = SpatialCoordinate(test_mesh)
f = sin(pi*x)*sin(pi*y)
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
assemble 作用在双线性形式上得到矩阵, 作用在线性形式上得到右端项. 注意两者的类型不同: 矩阵是 firedrake.matrix.Matrix, 右端项是 firedrake.cofunction.Cofunction (对偶空间中的元素, 不是 Function). 它们都只是包装, 数据存在内部的 PETSc 对象里.
A = assemble(a)
b = assemble(L)
type(A), type(b)
(firedrake.matrix.Matrix, firedrake.cofunction.Cofunction)
5.1.1. Matrix#
type(A.petscmat)
petsc4py.PETSc.Mat
单进程运行且矩阵不大时, 可以把 PETSc 矩阵转换为 numpy 数组
import numpy as np
from scipy.sparse import csr_matrix
m, n = A.petscmat.getSize()
indptr, indices, data = A.petscmat.getValuesCSR()
A_numpy = csr_matrix((data, indices, indptr), shape=(m, n)).toarray()
下面把同一行的两种表示放在一起对照. getRow 直接从 PETSc 矩阵读出一行, 返回 (列号数组, 数值数组); A_numpy[0, :] 则是转换后的稠密行. 对照可以看出, 稀疏表示里只出现属于矩阵稀疏模式的列, 而稀疏模式内的元素即使数值恰好为 0 也会被存下来.
A.petscmat.getRow(0), A_numpy[0, :]
((array([0, 1, 2, 3, 5], dtype=int32), array([ 2. , 0. , -0.5, -1. , -0.5])),
array([ 2. , 0. , -0.5, -1. , 0. , -0.5, 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. ]))
保存矩阵到文件 MatViewFromOptions
https://petsc.org/main/manualpages/Mat/MatViewFromOptions/
在代码中加入如下行
A.petscmat.viewFromOptions('-A_view')
那么在命令行可以通过选项 -A_view binary:A.bin 保存 A 到文件 A.bin.
5.1.2. 把 PETSc 矩阵转换为 csr_matrix 并计算行列式#
下面是一个完整的脚本: 照常组装并求解 Poisson 问题, 再从求解器的 KSP 上取回系数矩阵 (ksp.getOperators() 返回算子矩阵和预条件子矩阵两个对象), 转成 scipy 的稀疏矩阵后做 LU 分解, 用 L、U 对角元的乘积得到行列式的绝对值. 这样取到的矩阵已经施加了 Dirichlet 边界条件, 与直接 assemble(a) 得到的矩阵不同.
这里有一个容易忽略的陷阱: splu 分解的是 Pr*A*Pc = L*U, 因此 L、U 对角元的乘积只给出行列式的绝对值, 符号还取决于 lu.perm_r 和 lu.perm_c 两个置换的奇偶性. 本例中这两个置换的奇偶性恰好相同, 符号互相抵消, 直接相乘才看起来是对的; 换一个矩阵就未必了.
后面的 mpiexec 命令演示了如何用命令行选项观察求解过程: -ksp_view 打印求解器的完整配置, -ksp_monitor 打印每步残差, -ksp_error_if_not_converged 让不收敛时直接报错, -mat_mumps_icntl_33 1 让 MUMPS 顺便算出行列式, -mat_mumps_icntl_4 4 提高 MUMPS 的输出详细程度.
# %load py/poisson_example1.py
from firedrake import *
from firedrake.petsc import PETSc
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import splu
import numpy as np
N = 8
test_mesh = RectangleMesh(nx=N, ny=N, Lx=1, Ly=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 # or f*v*dx
bc = DirichletBC(V, g=g, sub_domain='on_boundary')
u_h = Function(V, name='u_h')
problem = LinearVariationalProblem(a, L, u_h, bcs=bc)
solver = LinearVariationalSolver(problem, options_prefix='')
solver.solve()
ksp = solver.snes.getKSP()
A, P = ksp.getOperators()
m, n = A.getSize()
indptr, indices, data = A.getValuesCSR()
A_numpy = csr_matrix((data, indices, indptr), shape=(m, n)).toarray()
print(A_numpy)
print(data)
lu = splu(A_numpy)
diagL = lu.L.diagonal()
diagU = lu.U.diagonal()
det=diagL.prod()*diagU.prod()
print(det)
mpiexec --bind-to none -n 1 \
python3 a.py \
-ksp_view \
-ksp_error_if_not_converged \
-mat_mumps_icntl_33 1 \
-mat_mumps_icntl_4 4 \
-ksp_monitor 2>&1 | tee "my-$(date +%Y%m%d-%H:%M:%S).log"
5.1.3. Vector#
Function 和 Cofunction 把数据存在 dat 中, 通过上下文管理器借出底层的 PETSc Vec: dat.vec 可读写, dat.vec_ro 只读, dat.vec_wo 只写. 以可写方式退出 with 块时, Firedrake 会把 halo 标记为需要重新交换 (并不是当场同步, 真正的交换推迟到下次用到 halo 时), 因此对 Vec 的操作要写在块内, 不要把它保存到块外使用. 另外这三个属性都带 @mpi.collective, 并行时必须所有进程一起进入 with 块.
with b.dat.vec_ro as vec:
print(type(vec))
<class 'petsc4py.PETSc.Vec'>
5.1.4. ISLocalToGlobalMapping#
局部到全局映射 (local-to-global map) 记录本进程上局部编号与全局编号的对应关系, 并行组装时靠它把本地的贡献累加到全局矩阵、向量的正确位置. 下面手工构造一个: 每个进程持有 rank + 3 个自由度, 先用 Scan 求出本进程在全局编号中的起始偏移, 再据此生成一段连续的全局编号. lgmap.view() 输出的每一行是 “进程号 局部编号 全局编号”. 笔记本里是单进程执行, 因此只有 3 个自由度, 全局编号就是 0, 1, 2.
import firedrake as fd
from firedrake.petsc import PETSc
from pyop2.datatypes import IntType, ScalarType
import numpy as np
rank = COMM_WORLD.rank
owned_sz = np.array(rank+3, dtype=IntType)
offset = np.empty_like(owned_sz)
COMM_WORLD.Scan(owned_sz, offset)
offset -= owned_sz
indices = np.arange(offset, offset + owned_sz, dtype=IntType)
lgmap = PETSc.LGMap()
lgmap.create(indices, bsize=1, comm=COMM_WORLD)
lgmap.view()
ISLocalToGlobalMapping Object: 1 MPI process
type not yet set
[0] 0 0
[0] 1 1
[0] 2 2
5.2. Options#
PETSc.Options() 用于读取命令行 (或脚本内设置) 传入的选项, 在 Firedrake 脚本中很常用. Firedrake 的 solver_parameters 最终也是写进同一个选项数据库的, 因此求解器参数既可以写在脚本里, 也可以在命令行上给出, 调参时不必反复改代码. 但命令行要生效有一个前提: Firedrake 默认会给每个求解器自动生成一个选项前缀 (形如 firedrake_1_), 此时命令行上裸写 -ksp_type cg 匹配不上, 而且不会报错, 求解器仍按默认的 preonly + lu 运行. 想在命令行上调参, 建求解器时要显式写 options_prefix='' (本章后面的两个脚本都是这么做的), 或者自定一个前缀并在命令行上带上它. 下面的例子读取两个自定义选项 -lcs 和 -datapath:
import os
import sys
from firedrake import *
if __name__ == '__main__':
opts = PETSc.Options()
if opts.hasName('lcs'): # -lcs
lcs_str = opts.getString('lcs')
lcs = eval(lcs_str)
else:
lcs = None
if opts.hasName('datapath'): # -datapath
datapath = opts.getString('datapath')
assert os.path.exists(datapath)
else:
datapath = None
print(lcs, datapath)
5.3. KSP#
自定义 KSP 进行线性方程组求解请参考 PETSc 的文档. Firedrake 通过 solver_parameters 配出来的求解器最终也是一个 PETSc KSP 对象, 下面记录两个常用的做法.
求解完成需要检查是否收敛
下面的例子故意构造一个无法完成 LU 分解的 2×2 矩阵, 用来演示失败之后如何取得原因: ksp.getConvergedReason() 给出 KSP 层面的收敛状态, ksp.getPC().getFailedReason() 给出预条件子层面的失败原因. 本例输出的是 DIVERGED_PCSETUP_FAILED 和 FACTOR_NUMERIC_ZEROPIVOT, 即预条件子 (这里是 LU 分解) 在数值分解时遇到了零主元.
代码中的 OptionsManager 负责把参数字典写进 PETSc 的选项数据库, with om.inserted_options() 保证这些选项只在该代码块内生效; 它在新版本中已经移到 petsctools 包, 因此这里用 try/except ImportError 做了兼容导入. 另外, PETSc 自身的错误以 PETSc.Error 抛出, 其 ierr 属性是 PETSc 的错误码, 例子中按错误码分别给出提示.
from firedrake.exceptions import ConvergenceError
from firedrake.petsc import PETSc
from firedrake.solving_utils import KSPReasons
import numpy as np
try:
from petsctools.options import OptionsManager
except ImportError:
from firedrake.petsc import OptionsManager
def _make_reasons(reasons):
return dict([(getattr(reasons, r), r)
for r in dir(reasons) if not r.startswith('_')])
PCFailedReason = _make_reasons(PETSc.PC.FailedReason())
def get_ksp_reason(ksp):
r = ksp.getConvergedReason()
pc = ksp.getPC()
r_pc = pc.getFailedReason()
return KSPReasons[r], PCFailedReason[r_pc]
A = PETSc.Mat()
A.create(PETSc.COMM_WORLD)
A.setSizes([2, 2])
A.setType('aij') # sparse
# A.setPreallocationNNZ(4)
A.setUp()
A.setValue(1, 0, 1)
A.setValue(0, 1, np.inf) # to make the solver failed
A.assemble()
ksp = PETSc.KSP().create()
ksp.setOperators(A) # solve A*x=b by ksp.solve(b,x)
om = OptionsManager(
{
'ksp_type': 'preonly',
'pc_type': 'lu',
# 'ksp_view': None,
'pc_factor_mat_solver_type': 'mumps',
# 'ksp_error_if_not_converged': None,
},
options_prefix='test')
om.set_from_options(ksp)
x, b = A.createVecs()
b.setValue(0, 1)
# ksp.view()
with om.inserted_options():
try:
ksp.solve(b, x)
r = ksp.getConvergedReason()
if r < 0:
raise ConvergenceError(KSPReasons[r])
except ConvergenceError as e:
r, r_pc = get_ksp_reason(ksp)
PETSc.Sys.Print(f"Error: solver did not converged: {r}, PC: {r_pc}")
except PETSc.Error as e:
if e.ierr == 91: # https://petsc.org/release/include/petscerror.h.html
PETSc.Sys.Print(f"Error from PETSc: solver did not converged: {KSPReasons[ksp.getConvergedReason()]}")
elif e.ierr == 76:
PETSc.Sys.Print(f"Error from PETSc:")
PETSc.Sys.Print(f" ksp reason: {KSPReasons[ksp.getConvergedReason()]}")
PETSc.Sys.Print(f" error in library called by PETSc:")
PETSc.Sys.Print(" "*4 + str(e).replace("\n", "\n" + " "*4))
# We should terminate the process when an error occured in petsc
# as suggested by Matt https://lists.mcs.anl.gov/pipermail/petsc-users/2023-March/048146.html
raise
Error: solver did not converged: DIVERGED_PCSETUP_FAILED, PC: FACTOR_NUMERIC_ZEROPIVOT
查看特征值和残差变化, 并保存图片
python test.py -ksp_type gmres -pc_type jacobi -ksp_view_eigenvalues draw -ksp_monitor draw::draw_lg -draw_save .png
5.3.1. 在 Firedrake 中检查 ksp 状态#
在 Firedrake 里, 同一个 KSP 对象可以从求解器上取到: solver.snes.getKSP(). 下面故意使用不带预条件子的共轭梯度法, 并把迭代次数限制为 ksp_max_it: 4, 让求解必定不收敛, 再在循环中捕获异常继续往下跑. 网格规模由 PETSc.Options() 读到的 -N 决定, 默认取 32*size.
输出中有两行来源不同: Linear solve did not converge due to DIVERGED_ITS iterations 4 是选项 ksp_converged_reason 让 PETSc 自己打印的; DIVERGED_MAX_IT 则是用 KSPReasons 按收敛原因码查出来的名字. 二者对应的是同一个原因码 (-3), 只是 PETSc 的 C 端与 petsc4py 端对它的命名不同, 不是两个不同的失败.
from firedrake import *
from firedrake.petsc import PETSc
from firedrake.solving_utils import KSPReasons
import numpy as np
def printf(*args, **kwargs):
PETSc.Sys.Print(*args, **kwargs)
def get_ksp_reason(solver):
r = solver.snes.getKSP().getConvergedReason()
return KSPReasons[r]
rank, size = COMM_WORLD.rank, COMM_WORLD.size
opts = PETSc.Options()
N = opts.getInt('N', 32*size)
test_mesh = RectangleMesh(nx=N, ny=N, Lx=1, Ly=1)
x, y = SpatialCoordinate(test_mesh)
f = sin(pi*x)*sin(pi*y)
V = FunctionSpace(test_mesh, 'CG', degree=1)
u, v = TrialFunction(V), TestFunction(V)
a = inner(grad(u), grad(v))*dx - inner(f, v)*dx
bc = DirichletBC(V, 0, sub_domain='on_boundary')
u_h = Function(V, name='u_h')
problem = LinearVariationalProblem(lhs(a), rhs(a), u_h, bcs=bc)
solver_parameters = {'ksp_type': 'cg',
'ksp_max_it': 4,
'ksp_converged_reason': None,
# 'ksp_error_if_not_converged': None,
'pc_type': 'none'}
solver = LinearVariationalSolver(problem, solver_parameters=solver_parameters, options_prefix='')
for i in range(3):
printf(f"Loop i = {i}")
try:
solver.solve()
except ConvergenceError:
printf(f" Error from Firedrake: solver did not converged: {get_ksp_reason(solver)}")
except PETSc.Error as e:
if e.ierr == 91: # https://petsc.org/release/include/petscerror.h.html
printf(f" Error from PETSc: solver did not converged: {get_ksp_reason(solver)}")
elif e.ierr == 76:
printf(" Error from PETSc:")
printf(f" ksp reason: {get_ksp_reason(solver)}")
printf(" error in library called by PETSc:")
printf(" "*6 + str(e).replace("\n", "\n" + " "*6))
# We should terminate the process when an error occured in petsc
# as suggested by Matt https://lists.mcs.anl.gov/pipermail/petsc-users/2023-March/048146.html
raise
Loop i = 0
Linear solve did not converge due to DIVERGED_ITS iterations 4
Error from Firedrake: solver did not converged: DIVERGED_MAX_IT
Loop i = 1
Linear solve did not converge due to DIVERGED_ITS iterations 4
Error from Firedrake: solver did not converged: DIVERGED_MAX_IT
Loop i = 2
Linear solve did not converge due to DIVERGED_ITS iterations 4
Error from Firedrake: solver did not converged: DIVERGED_MAX_IT
5.4. DMPlex#
DMPlex 是管理非结构化网格的数据结构, 内部使用有向无环图表示网格的拓扑关系 [HKA+21, LMKG16].
DMPlex 中的 Plex 表示 Complex. 并行计算时, 网格会被划分成不同的块, 分配到各个进程.
Firedrake 的常规网格背后都有一个 DMPlex, 可以用 mesh.topology_dm 取到 (VertexOnlyMesh 是例外, 它背后是 DMSwarm). 下面先给出一个直接用 DMPlex 做网格粗化的脚本, 借助 Viewer 观察 DMPlex 的例子放在下一节.
5.4.1. 网格粗化#
下面的脚本用 PETSc.DMPlex 建一张立方体网格, 分发到各个进程, 再调用 coarsen() 得到一张更粗的网格. 三处 viewFromOptions 分别在初始、分发后和粗化后输出网格, 输出目标由命令行选项给出 (格式见下一节), 因此不加选项时它们什么都不做. coarsen() 依赖 PETSc 的网格适应模块, 用哪个库由 -dm_adaptor 决定 (末尾注释中的命令用的是 parmmg), 需要 PETSc 构建时带入相应的库. 这段代码以独立脚本的形式给出 (可以用 mpiexec 并行运行), 没有在笔记本中执行.
import sys
import petsc4py
petsc4py.init(sys.argv)
from petsc4py import PETSc
def output_vtk(dmplex, filename):
viewer = PETSc.Viewer().createVTK(filename, 'w')
viewer.view(dmplex)
opts = PETSc.Options()
N = opts.getInt('N', 4)
dim = opts.getInt('dim', 3)
overlap = opts.getInt('overlap', 1)
faces = [N for _ in range(dim)]
plex = PETSc.DMPlex().createBoxMesh(faces, simplex=True)
plex.setName('Init DM')
plex.viewFromOptions('-init_dm_view')
sf = plex.distribute(overlap=overlap)
plex.setName('Distribue DM')
plex.viewFromOptions('-dist_dm_view')
new_plex = plex.coarsen()
new_plex.setName('Coarsen DM')
new_plex.viewFromOptions('-coarsen_dm_view')
# mpiexec -n 2 python test_coarsen.py -dim 3 -overlap 0 -dm_adaptor parmmg -coarsen_dm_view vtk:data/test.vtu
5.5. Viewer#
Viewer 是 PETSc 输出对象的统一接口. 大多数 PETSc 对象都有两个查看方法: view() 直接打印到屏幕, viewFromOptions('-xxx_view') 则把输出目标和格式交给命令行选项决定. 后者的好处是不改代码就能把同一个对象打印到屏幕、写成二进制或 HDF5 文件, 甚至画成图.
选项字符串的语法 (
PetscObjectViewFromOptions): https://petsc.org/main/manualpages/Sys/PetscObjectViewFromOptions/
Viewer 的选项取值:
If no value is provided ascii:stdout is used
ascii[:[filename][:[format][:append]]] defaults to stdout - format can be one of ascii_info, ascii_info_detail, or ascii_matlab,
for example ascii::ascii_info prints just the information about the object not all details
unless :append is given filename opens in write mode, overwriting what was already there
binary[:[filename][:[format][:append]]] defaults to the file binaryoutput
draw[:drawtype[:filename]] for example, draw:tikz, draw:tikz:figure.tex or draw:x
socket[:port] defaults to the standard output port
saws[:communicatorname] publishes object to the Scientific Application
Webserver (SAWs)
绘图类 viewer 的选项 (
PetscDrawSetFromOptions): https://petsc.org/main/manualpages/Draw/PetscDrawSetFromOptions/
-nox - do not use X graphics (ignore graphics calls, but run program correctly)
-nox_warning - when X Windows support is not installed this prevents the warning message from being printed
-draw_pause <pause amount> -- -1 indicates wait for mouse input,
-2 indicates pause when window is to be destroyed
-draw_marker_type - <x,point>
-draw_save [optional filename] - (X Windows only) saves each image before it is cleared to a file
-draw_save_final_image [optional filename] - (X Windows only) saves the final image displayed in a window
-draw_save_movie - converts image files to a movie at the end of the run. See PetscDrawSetSave()
-draw_save_single_file - saves each new image in the same file, normally each new image is saved in a new file with 'filename/filename_%d.ext'
-draw_save_on_clear - saves an image on each clear, mainly for debugging
-draw_save_on_flush - saves an image on each flush, mainly for debugging
5.5.1. 用 petsc4py 读取并查看网格文件#
下面不经过 Firedrake, 直接用 petsc4py 把一个 gmsh 网格读成 DMPlex, 用 dm.view() 打印它的基本信息, 再用不同的 Viewer 写出几种文件: createHDF5 写成 HDF5 便于以后重新载入; 换成 HDF5_XDMF 格式再写一份, 可以用 PETSc 自带的 petsc_gen_xdmf.py 生成 .xdmf 文件供 ParaView 可视化; createVTK 则直接写成 VTK 文件.
dm.view() 输出中的 Number of 0-cells、1-cells、2-cells 分别是顶点、边和单元的个数, Labels 一节列出网格上的标签, 其中 depth 和 celltype 是 DMPlex 自己维护的拓扑信息.
import os
import sys
import petsc4py
petsc4py.init(sys.argv)
from petsc4py import PETSc
import numpy as np
dm = PETSc.DMPlex().createFromFile('gmsh/Lshape.msh', plexname='test')
dm.view()
# make sure dir "data" exists
os.makedirs("data", exist_ok=True)
# hdf5 for load
viewer = PETSc.Viewer().createHDF5('data/Lshape.h5', mode='w')
viewer(dm)
# hdf5 for visualization:
# You can generate xdmf file from this file by
# `petsc/lib/petsc/bin/petsc_gen_xdmf.py`
# Then load the xdmf file to paraview to visualize the mesh.
viewer = PETSc.Viewer().createHDF5('data/Lshape_xdmf.h5', mode='w')
viewer.pushFormat(viewer.Format.HDF5_XDMF)
viewer(dm)
viewer.popFormat()
# vtk file
viewer = PETSc.Viewer().createVTK('data/Lshape.vtk', mode='w')
viewer(dm)
# draw on X window
# viewer = PETSc.Viewer().createDraw()
# viewer(dm)
DM Object: test 1 MPI process
type: plex
test in 2 dimensions:
Number of 0-cells per rank: 274
Number of 1-cells per rank: 755
Number of 2-cells per rank: 482
Labels:
celltype: 3 strata with value/size (0 (274), 1 (755), 3 (482))
depth: 3 strata with value/size (0 (274), 1 (755), 2 (482))
petsc4py 没有把 PetscDraw 对象暴露出来, 它的 Viewer().createDraw() 也无法指定 tikz 这类输出格式和目标文件名, 因此想用 draw 类型的 viewer 只能走选项数据库: 先用 insertString 把 -dm_view 写进 PETSc.Options(), 再调用 viewFromOptions('-dm_view') 让对象按该选项输出. 这里用 draw:tikz:data/Lshape.tex 把网格画成 TikZ 代码存成 .tex 文件, 可以直接插进 LaTeX 文档. 注意 insertString 写入的是全局的选项数据库, 之后同一进程中所有 viewFromOptions('-dm_view') 都会受它影响.
# since the petsc_draw is not in petsc4py, we use options to save the images
opts = PETSc.Options()
opts_old = opts.getAll()
opts.insertString('-dm_view draw:tikz:data/Lshape.tex')
dm.viewFromOptions('-dm_view')
5.5.2. 用 DMPlex 查看 Firedrake 网格#
Firedrake 网格底层的 DMPlex 用 mesh.topology_dm 取到. 网格分区的信息只有在多进程下才看得出来, 因此下面借助 ipyparallel 启动 2 个 MPI 进程, 之后带 %%px 的单元都在这两个进程上执行.
import ipyparallel as ipp
cluster = ipp.Cluster(engines="mpi", n=2)
client = cluster.start_and_connect_sync()
Starting 2 engines with <class 'ipyparallel.cluster.launcher.MPIEngineSetLauncher'>
下面在 2 个进程上建同一张 8×8 的网格并打印它的 DMPlex. 输出中 Number of 0-cells per rank 这几行后面各有两个数字, 分别是两个进程本地 DMPlex 上的顶点、边、单元个数, 它们之和比串行时的总数大, 因为分区界面附近的实体在相邻进程上各有一份. 要注意这是本地的数量 (含 halo), 不是 “本进程持有” 的数量: 真正持有的数量求和应当正好等于串行时的总数. pyop2_core、pyop2_owned、pyop2_ghost 三个标签是 Firedrake 加上去的, 把网格点分成三类: ghost 是从别的进程借来的一层, owned 是本进程持有并且需要发送给别人的部分, core 是本进程持有且不需要发送的部分, 因此 core 上的计算可以在通信进行的同时先做.
%%px --block
from firedrake import *
mesh = RectangleMesh(8, 8, 1, 1)
mesh.topology_dm.view()
[stdout:0] DM Object: DM_0x24a822c0_0 2 MPI processes
type: plex
DM_0x24a822c0_0 in 2 dimensions:
Number of 0-cells per rank: 61 71
Number of 1-cells per rank: 142 153
Number of 2-cells per rank: 82 83
Labels:
depth: 3 strata with value/size (0 (61), 1 (142), 2 (82))
celltype: 3 strata with value/size (0 (61), 1 (142), 3 (82))
marker: 1 strata with value/size (1 (18))
Face Sets: 3 strata with value/size (2 (5), 3 (4), 4 (9))
exterior_facets: 1 strata with value/size (1 (18))
interior_facets: 1 strata with value/size (1 (199))
pyop2_core: 1 strata with value/size (1 (72))
pyop2_owned: 1 strata with value/size (1 (104))
pyop2_ghost: 1 strata with value/size (1 (109))
5.6. Star Forest 结构体 PetscSF#
PetscSF 是 PETSc 中用于进程间数据交换的数据结构.
它存储了根节点(当前进程上)和叶子节点(任意进程上)的对应关系, 可以把数据从根节点发送到叶子节点, 也可以把数据从叶子节点收集到根节点 [ZBB+22].
下面的例子在 rank 0 上建一张单位正方形的三角形网格 (9 个顶点, 8 个单元), 在顶点上定义一个每点 1 个自由度的 Section, 再把网格分发到 2 个进程: distribute() 返回的 SF 描述网格点的对应关系, 用它调用 distributeSection 可以把 Section 一并分发, 得到叶子进程上的 Section 与偏移, 最后由 createSectionSF 得到自由度层面的 SF. 函数里反复调用 distributeSection, 是为了验证其中几种调用写法 (是否预先创建 leafSection) 给出的结果一致.
%%px --block
from firedrake import *
from firedrake.petsc import PETSc
from petsc4py import PETSc
import numpy as np
# 6--------7--------8
# | | |
# 3--------4--------5
# | | |
# 0--------1--------2
def test_SFDistributeSection():
comm = COMM_WORLD
if comm.rank == 0:
cells = np.asarray(
[[0, 1, 3],
[1, 2, 4],
[1, 4, 3],
[2, 5, 4],
[3, 4, 6],
[4, 5, 7],
[4, 7, 6],
[5, 8, 7]], dtype=np.int32)
coords = np.asarray(
[[0. , 0. ],
[0.5, 0. ],
[1. , 0. ],
[0. , 0.5],
[0.5, 0.5],
[1.0, 0.5],
[0. , 1. ],
[0.5, 1. ],
[1. , 1. ]], dtype=np.double)
else:
cells = np.zeros([0, 3], dtype=np.int32)
coords = np.zeros([0, 2], dtype=np.double)
dim = 2
plex = PETSc.DMPlex().createFromCellList(dim, cells, coords, comm=comm)
rootSection = PETSc.Section().create(comm=comm)
pStart, pEnd = plex.getHeightStratum(2)
rootSection.setChart(*plex.getChart())
for p in range(pStart, pEnd):
rootSection.setDof(p, 1)
rootSection.setUp()
rootSection.viewFromOptions('-section_view')
dplex = plex.clone()
msf = dplex.distribute()
if msf is None:
PETSc.Sys.Print("Warning: plex has not been distributed!")
return
dplex.viewFromOptions('-dm_view')
def isEqualSF(ssf0, ssf1):
nroots0, local0, remote0 = ssf0.getGraph()
nroots1, local1, remote1 = ssf1.getGraph()
return (nroots0 == nroots1) \
and np.array_equal(local0, local1) \
and np.array_equal(remote0, remote1)
remoteOffsets0, leafSection0 = msf.distributeSection(rootSection)
ssf0 = msf.createSectionSF(rootSection, remoteOffsets0, leafSection0)
remoteOffsets1, leafSection1 = msf.distributeSection(rootSection, None)
ssf1 = msf.createSectionSF(rootSection, remoteOffsets1, leafSection1)
leafSection2 = PETSc.Section()
remoteOffsets2, leafSection2 = msf.distributeSection(rootSection, leafSection2)
ssf2 = msf.createSectionSF(rootSection, remoteOffsets2, leafSection2)
leafSection3 = PETSc.Section()
remoteOffsets3, _ = msf.distributeSection(rootSection, leafSection3)
ssf3 = msf.createSectionSF(rootSection, remoteOffsets3, leafSection3)
leafSection4 = PETSc.Section().create(dplex.getComm())
remoteOffsets4, leafSection4 = msf.distributeSection(rootSection, leafSection4)
ssf4 = msf.createSectionSF(rootSection, remoteOffsets4, leafSection4)
leafSection5 = PETSc.Section().create(dplex.getComm())
remoteOffsets5, _ = msf.distributeSection(rootSection, leafSection5)
ssf5 = msf.createSectionSF(rootSection, remoteOffsets5, leafSection5)
assert isEqualSF(ssf0, ssf1)
assert isEqualSF(ssf0, ssf2)
assert isEqualSF(ssf0, ssf3)
assert isEqualSF(ssf0, ssf4)
ssf0.view()
上一个单元只是定义函数, 下面真正调用它. 输出的是分发之后的自由度 SF: rank 0 上有 9 个根节点 (原网格的 9 个顶点各 1 个自由度, 分发前整张网格都在 rank 0 上, 所以 rank 1 的根节点数是 0), 两个进程各有 6 个叶子节点; 0 <- (0,3) 表示本进程的第 0 个叶子对应 rank 0 上的第 3 个根. 两边叶子加起来是 12 而不是 9, 多出来的正是分区界面上被两个进程共享的顶点.
%%px --block
test_SFDistributeSection()
[stdout:0] PetscSF Object: 2 MPI processes
type: basic
[0] Number of roots=9, leaves=6, remote ranks=1
[0] 0 <- (0,1)
[0] 1 <- (0,2)
[0] 2 <- (0,4)
[0] 3 <- (0,5)
[0] 4 <- (0,7)
[0] 5 <- (0,8)
[1] Number of roots=0, leaves=6, remote ranks=1
[1] 0 <- (0,0)
[1] 1 <- (0,1)
[1] 2 <- (0,3)
[1] 3 <- (0,4)
[1] 4 <- (0,6)
[1] 5 <- (0,7)
MultiSF sort=rank-order
5.7. PCPatch#
PCPatch 是在单元片 (patch) 上构造的预条件子: 把整体问题拆成许多个小的局部问题, 分别求解再拼起来, 常用作多重网格的光滑子. 用户可以使用 PETSc 提供的单元片类型 (star, vanka, pardecomp), 也可以自行构造.
它根据用户输入的关于解向量的 dm 和 sf 信息, 构造各个单元片上的自由度和全局自由度的对应关系. 下面记录的是 PCPatch 的内部变量, 只有在需要自行构造单元片, 或者要读懂 PETSc 相关源码时才需要了解.
/* Topology */
PetscInt dim, codim; /* Dimension or codimension of mesh points to loop over; only one of them can be set */
PetscSection cellCounts; /* Maps patch -> # cells in patch */
IS cells; /* [patch][cell in patch]: Cell number */
PetscSection pointCounts; /* Maps patch -> # points with dofs in patch */
IS points; /* [patch][point in patch]: Point number */
PetscSection intFacetCounts;
PetscSection extFacetCounts;
PetscSection cellNumbering; /* Plex: NULL Firedrake: Numbering of cells in DM */
/* Dof layout */
PetscInt nsubspaces; /* Number of fields */
PetscSF sectionSF; /* Combined SF mapping process local to global */
PetscSection *dofSection; /* ?? For each field, patch -> # dofs in patch */
PetscInt *subspaceOffsets; /* Plex: NULL Firedrake: offset of each field in concatenated process local numbering for mixed spaces */
PetscInt **cellNodeMap; /* [field][cell][dof in cell]: global dofs in cell TODO Free this after its use in PCPatchCreateCellPatchDiscretisationInfo() */
IS dofs; /* [patch][cell in patch][dof in cell]: patch local dof */
IS offs; /* [patch][point in patch]: patch local offset (same layout as 'points', used for filling up patchSection) */
PetscSection gtolCounts; /* ?? Indices to extract from local to patch vectors */
IS gtol;
gtol 是从本进程自由度编号到单个 patch 上局部自由度的映射. 本进程自由度相对整个解是局部编号, 但相对单个 patch 则是全局编号.
dofs 是 patch 上局部的自由度映射关系 (可以看作 patch 上的 cell_node_map, 和传入的 map 一致)
References
Vaclav Hapla, Matthew G. Knepley, Michael Afanasiev, Christian Boehm, Martin van Driel, Lion Krischer, and Andreas Fichtner. Fully parallel mesh I/O using PETSc DMPlex with an application to waveform modeling. SIAM Journal on Scientific Computing, 43(2):C127–C153, jan 2021. doi:10.1137/20m1332748.
Michael Lange, Lawrence Mitchell, Matthew G. Knepley, and Gerard J. Gorman. Efficient mesh management in Firedrake using PETSc DMPlex. SIAM Journal on Scientific Computing, 38(5):S143–S155, jan 2016. doi:10.1137/15m1026092.
Junchao Zhang, Jed Brown, Satish Balay, Jacob Faibussowitsch, Matthew Knepley, Oana Marin, Richard Tran Mills, Todd Munson, Barry F. Smith, and Stefano Zampini. The PetscSF scalable communication layer. IEEE Transactions on Parallel and Distributed Systems, 33(4):842–853, apr 2022. doi:10.1109/tpds.2021.3084070.