PETSc

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.

概要

Firedrake 底层用到的 PETSc 对象一览 (通过 petsc4py): 向量与矩阵、KSP 求解器、查看器、star forest PetscSF, 以及 PCPatch 预条件子.

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 入门课程开始

  1. 视频: https://youtu.be/4Y8g-DcTreY

  2. 讲义: https://web.corral.tacc.utexas.edu/CompEdu/pdf/pcse/petsc_p_course.pdf

PETSc 网站的手册和入门讲义

  1. 手册: https://petsc.org/release/manual/

  2. 入门讲义: https://petsc.org/release/tutorials/

  3. C/Fortran API: https://petsc.org/release/manualpages/

  4. petsc4py: https://petsc.org/release/petsc4py/

PETSc 代码库有许多示例可以作为学习材料, 如 PETSc 仓库中 petsc4py 的示例:

  1. petsc/petsc/-/tree/main/src/binding/petsc4py/demo

PETSc is a portable, extensible toolkit for scientific computation developed at Argonne National Laboratory. It provides many algorithms for solving large systems of equations in parallel, and it can also call external packages to solve them. PETSc further provides components for the numerical solution of partial differential equations, including structured mesh data structures, unstructured mesh data structures and finite element spaces, and so on.

Firedrake can be regarded as a high-level wrapper around PETSc. In Firedrake one usually does not need to manipulate PETSc objects, but in some special cases PETSc objects have to be manipulated directly, and sometimes doing so is more efficient.

PETSc is a C library that also provides a Fortran interface. petsc4py is the python wrapper of PETSc.

A good starting point for learning PETSc is the introductory PETSc course published by the Texas Advanced Computing Center (TACC)

  1. Video: https://youtu.be/4Y8g-DcTreY

  2. Slides: https://web.corral.tacc.utexas.edu/CompEdu/pdf/pcse/petsc_p_course.pdf

The manual and the introductory tutorials on the PETSc website

  1. Manual: https://petsc.org/release/manual/

  2. Tutorials: https://petsc.org/release/tutorials/

  3. C/Fortran API: https://petsc.org/release/manualpages/

  4. petsc4py: https://petsc.org/release/petsc4py/

The PETSc code base contains many examples that can serve as learning material, such as the petsc4py demos in the PETSc repository:

  1. petsc/petsc/-/tree/main/src/binding/petsc4py/demo

Firedrake 已经把 PETSc 包好了, 日常使用时通过 solver_parameters 字典设置求解器就够了. 需要往下一层直接操作 PETSc 对象的, 通常是下面几种情形:

  1. 矩阵或向量本身: 例如把组装好的矩阵取出来做谱分析、算行列式, 或者按 PETSc 的格式存盘交给别的程序处理;

  2. 求解器为什么失败: Firedrake 只会告诉你没有收敛, 具体的收敛状态和预条件子的失败原因要从 KSP、PC 对象上读;

  3. 网格拓扑: 网格分区、粗化、编号以及并行的数据交换都发生在 DMPlex 和 PetscSF 这一层, Firedrake 的 Mesh 只是它们的包装;

  4. Firedrake 尚未包装的功能: 例如 Viewer 的各种输出格式, 或者 PCPatch 这类要求用户自行提供拓扑信息的预条件子.

各节彼此独立, 可以按需查阅:

  • Vector 和 Matrix: 取出底层的 Mat、Vec, 以及与 numpy/scipy 之间的转换;

  • Options: PETSc 的选项数据库, 是所有 PETSc 参数的统一入口;

  • KSP: 线性求解器对象, 重点是求解失败之后如何取得原因;

  • DMPlex、Viewer、PetscSF: 网格拓扑、对象输出与进程间的数据交换;

  • PCPatch: 单元片预条件子的内部数据结构, 是最深的一层.

不过代码上并不独立: 后面带 %%px 的单元都跑在 Viewer 一节里启动的 ipyparallel 集群上, 想单独执行它们, 要先运行启动集群的那一格.

Firedrake already wraps PETSc, and for everyday use it is enough to configure the solver through the solver_parameters dictionary. Going one level down and manipulating PETSc objects directly is usually needed in these situations:

  1. The matrix or the vector itself: for example taking the assembled matrix out for a spectral analysis or for computing a determinant, or saving it in PETSc format for another program to process;

  2. Why a solver failed: Firedrake only tells you that it did not converge; the actual convergence status and the reason for the failure of the preconditioner have to be read from the KSP and PC objects;

  3. The mesh topology: mesh partitioning, coarsening, numbering and the parallel data exchange all happen at the level of DMPlex and PetscSF, and Firedrake’s Mesh is only a wrapper around them;

  4. Functionality Firedrake has not wrapped yet: for example the various output formats of Viewer, or preconditioners such as PCPatch that require the user to supply the topological information.

The sections are independent of each other and can be consulted as needed:

  • Vector and Matrix: taking out the underlying Mat and Vec, and converting to and from numpy/scipy;

  • Options: the PETSc options database, the single entry point for all PETSc parameters;

  • KSP: the linear solver object, with the emphasis on how to obtain the reason after a failed solve;

  • DMPlex, Viewer, PetscSF: mesh topology, object output and data exchange between processes;

  • PCPatch: the internal data structures of the patch preconditioner, the deepest level.

The code, however, is not independent: all the cells below that carry %%px run on the ipyparallel cluster started in the Viewer section, so in order to execute them on their own the cell that starts the cluster has to be run first.

Tip

PETSc 目录中有用的工具, 如 h5dump, petsc_gen_xdmf.py, PetscBinaryIO.py 等.Useful tools live in the PETSc directory, such as h5dump, petsc_gen_xdmf.py and PetscBinaryIO.py.

在 PETSc 环境中, 运行如下命令添加这些工具所在路径到 PATH:In a PETSc environment, running the following commands adds the directories of these tools to PATH:

export PATH="$PATH:$PETSC_DIR/lib/petsc/bin"
export PATH="$PATH:$PETSC_DIR/${PETSC_ARCH-default}/bin"

在激活的 Firedrake 环境可以运行如下命令的输出, 添加工具所在路径到环境变量 PATH.In an activated Firedrake environment, running the output of the following command adds the directories of these tools to the environment variable 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 Vector and Matrix#

Firedrake 组装出来的矩阵和右端项都是 PETSc 对象的包装. 下面用一个 4×4 的矩形网格上的 Poisson 问题作为例子.

PETSc 官方的矩阵、向量读写示例见 matvecio.py.

The matrices and right-hand sides assembled by Firedrake are all wrappers around PETSc objects. A Poisson problem on a 4×4 rectangular mesh is used as the example below.

The official PETSc example for reading and writing matrices and vectors is 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 对象里.

assemble applied to a bilinear form gives a matrix, and applied to a linear form gives a right-hand side. Note that the two have different types: the matrix is a firedrake.matrix.Matrix, while the right-hand side is a firedrake.cofunction.Cofunction (an element of the dual space, not a Function). Both are only wrappers; the data lives in the PETSc objects inside them.

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 数组

In a single-process run, and when the matrix is not large, the PETSc matrix can be converted into a numpy array

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 也会被存下来.

The two representations of the same row are put side by side below. getRow reads a row directly from the PETSc matrix and returns (an array of column indices, an array of values); A_numpy[0, :] is the dense row after the conversion. The comparison shows that only the columns belonging to the sparsity pattern of the matrix appear in the sparse representation, and that an entry inside the sparsity pattern is stored even when its value happens to be 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/

在代码中加入如下行

Saving a matrix to a file with MatViewFromOptions

https://petsc.org/main/manualpages/Mat/MatViewFromOptions/

Add the following line to the code

A.petscmat.viewFromOptions('-A_view')

那么在命令行可以通过选项 -A_view binary:A.bin 保存 A 到文件 A.bin.

and then the option -A_view binary:A.bin on the command line saves A to the file A.bin.

5.1.2. 把 PETSc 矩阵转换为 csr_matrix 并计算行列式 Converting a PETSc matrix to a csr_matrix and computing the determinant#

下面是一个完整的脚本: 照常组装并求解 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 的输出详细程度.

The following is a complete script: it assembles and solves the Poisson problem as usual, takes the coefficient matrix back from the KSP of the solver (ksp.getOperators() returns two objects, the operator matrix and the preconditioner matrix), converts it into a scipy sparse matrix and factorizes it with LU, and obtains the absolute value of the determinant from the product of the diagonal entries of L and U. The matrix obtained in this way already has the Dirichlet boundary conditions applied, so it differs from the matrix obtained directly from assemble(a).

There is a pitfall here that is easily overlooked: splu factorizes Pr*A*Pc = L*U, so the product of the diagonal entries of L and U only gives the absolute value of the determinant; the sign further depends on the parities of the two permutations lu.perm_r and lu.perm_c. In this example the two permutations happen to have the same parity, so the signs cancel each other and the plain product happens to look right; for another matrix this need not be the case.

The mpiexec command below demonstrates how to watch the solution process with command-line options: -ksp_view prints the full configuration of the solver, -ksp_monitor prints the residual at every step, -ksp_error_if_not_converged makes a failure to converge raise an error immediately, -mat_mumps_icntl_33 1 makes MUMPS compute the determinant along the way, and -mat_mumps_icntl_4 4 increases the verbosity of the MUMPS output.

# %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 块.

Function and Cofunction store their data in dat, and lend out the underlying PETSc Vec through context managers: dat.vec is read-write, dat.vec_ro is read-only and dat.vec_wo is write-only. On leaving a with block that was entered for writing, Firedrake marks the halo as needing to be exchanged again (it is not synchronized on the spot; the actual exchange is deferred until the halo is used next), so the operations on the Vec have to be written inside the block and it must not be kept for use outside it. In addition, all three properties carry @mpi.collective, so in parallel all processes have to enter the with block together.

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.

The local-to-global map records the correspondence between the local numbering and the global numbering on the current process; during parallel assembly it is what accumulates the local contributions into the right places of the global matrix and vector. One is constructed by hand below: every process owns rank + 3 degrees of freedom, Scan first computes the starting offset of this process within the global numbering, and a contiguous range of global numbers is generated from it. Every line printed by lgmap.view() is “ ”. The notebook is executed on a single process, so there are only 3 degrees of freedom and the global numbering is simply 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:

PETSc.Options() reads the options passed on the command line (or set inside the script), and is very common in Firedrake scripts. Firedrake’s solver_parameters also end up in the same options database, so the solver parameters can be written either in the script or given on the command line, and there is no need to edit the code over and over while tuning. There is a precondition for the command line to take effect, though: by default Firedrake generates an options prefix for every solver automatically (of the form firedrake_1_), and a bare -ksp_type cg on the command line then does not match, and no error is reported; the solver still runs with the default preonly + lu. In order to tune parameters on the command line, the solver has to be created with an explicit options_prefix='' (both scripts later in this chapter do this), or with a prefix of one’s own that is then given on the command line. The example below reads two custom options, -lcs and -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 对象, 下面记录两个常用的做法.

  1. 求解完成需要检查是否收敛

下面的例子故意构造一个无法完成 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 的错误码, 例子中按错误码分别给出提示.

For building a custom KSP to solve a linear system, please refer to the PETSc documentation. The solver that Firedrake configures through solver_parameters is in the end also a PETSc KSP object; two common practices are recorded below.

  1. Checking whether the solve converged

The example below deliberately constructs a 2×2 matrix whose LU factorization cannot be completed, in order to demonstrate how to obtain the reason after a failure: ksp.getConvergedReason() gives the convergence status at the KSP level, and ksp.getPC().getFailedReason() gives the reason for the failure at the preconditioner level. This example prints DIVERGED_PCSETUP_FAILED and FACTOR_NUMERIC_ZEROPIVOT: the preconditioner (an LU factorization here) encountered a zero pivot during the numerical factorization.

The OptionsManager in the code writes the parameter dictionary into the PETSc options database, and with om.inserted_options() makes sure that these options are only in effect inside that block; OptionsManager has been moved to the petsctools package in recent versions, which is why the import here is guarded by a try/except ImportError. Furthermore, PETSc’s own errors are raised as PETSc.Error, whose ierr attribute is the PETSc error code; the example gives a separate message for each error code.

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

  1. 查看特征值和残差变化, 并保存图片Viewing the eigenvalues and the residual history, and saving the pictures

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 状态 Inspecting the ksp status in Firedrake#

在 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 端对它的命名不同, 不是两个不同的失败.

In Firedrake the same KSP object can be obtained from the solver with solver.snes.getKSP(). The conjugate gradient method is used below without a preconditioner on purpose, with the number of iterations limited by ksp_max_it: 4 so that the solve is bound not to converge, and the exception is then caught inside the loop so that the run carries on. The mesh size is determined by the -N read through PETSc.Options(), and defaults to 32*size.

Two of the lines in the output come from different places: Linear solve did not converge due to DIVERGED_ITS iterations 4 is printed by PETSc itself because of the option ksp_converged_reason, whereas DIVERGED_MAX_IT is the name looked up by the converged reason code with KSPReasons. The two correspond to the same reason code (-3); it is only that PETSc’s C side and the petsc4py side name it differently, not two different failures.

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 的例子放在下一节.

DMPlex is the data structure that manages unstructured meshes; internally it represents the topological relations of the mesh with a directed acyclic graph [HKA+21, LMKG16]. The Plex in DMPlex stands for Complex. In a parallel computation the mesh is divided into pieces that are distributed over the processes.

Every ordinary Firedrake mesh has a DMPlex behind it, which is obtained with mesh.topology_dm (VertexOnlyMesh is the exception: it is backed by a DMSwarm). A script that coarsens a mesh directly with DMPlex is given first below; the example that inspects a DMPlex with a Viewer is left to the next section.

5.4.1. 网格粗化 Mesh coarsening#

下面的脚本用 PETSc.DMPlex 建一张立方体网格, 分发到各个进程, 再调用 coarsen() 得到一张更粗的网格. 三处 viewFromOptions 分别在初始、分发后和粗化后输出网格, 输出目标由命令行选项给出 (格式见下一节), 因此不加选项时它们什么都不做. coarsen() 依赖 PETSc 的网格适应模块, 用哪个库由 -dm_adaptor 决定 (末尾注释中的命令用的是 parmmg), 需要 PETSc 构建时带入相应的库. 这段代码以独立脚本的形式给出 (可以用 mpiexec 并行运行), 没有在笔记本中执行.

The script below builds a cube mesh with PETSc.DMPlex, distributes it over the processes, and then calls coarsen() to obtain a coarser mesh. The three viewFromOptions calls output the mesh initially, after the distribution and after the coarsening; the output target is given by command-line options (the format is described in the next section), so without any options they do nothing. coarsen() relies on PETSc’s mesh adaptation module, and which library is used is decided by -dm_adaptor (the command in the trailing comment uses parmmg), which requires PETSc to have been built with the corresponding library. This code is given as a standalone script (it can be run in parallel with mpiexec) and is not executed in the notebook.

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 文件, 甚至画成图.

Viewer is PETSc’s unified interface for outputting objects. Most PETSc objects have two viewing methods: view() prints directly to the screen, while viewFromOptions('-xxx_view') leaves the output target and format to the command-line options. The advantage of the latter is that the same object can be printed to the screen, written as a binary or HDF5 file, or even drawn as a picture, without changing the code.

  1. 选项字符串的语法Syntax of the option string (PetscObjectViewFromOptions): https://petsc.org/main/manualpages/Sys/PetscObjectViewFromOptions/

Viewer 的选项取值:The values the viewer option can take:

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)

  1. 绘图类 viewer 的选项Options of the drawing viewers (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 读取并查看网格文件 Reading and viewing a mesh file with 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 自己维护的拓扑信息.

Below a gmsh mesh is read into a DMPlex directly with petsc4py, without going through Firedrake; dm.view() prints its basic information, and several kinds of files are then written with different Viewers: createHDF5 writes HDF5, which is convenient for reloading later; writing another copy in the HDF5_XDMF format allows PETSc’s own petsc_gen_xdmf.py to generate an .xdmf file for visualization in ParaView; createVTK writes a VTK file directly.

In the output of dm.view(), Number of 0-cells, 1-cells and 2-cells are the numbers of vertices, edges and cells respectively, and the Labels section lists the labels on the mesh, among which depth and celltype are the topological information that DMPlex maintains itself.

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') 都会受它影响.

petsc4py does not expose the PetscDraw object, and its Viewer().createDraw() cannot specify an output format such as tikz or a target file name either, so the only way to use a viewer of draw type is through the options database: insertString first writes -dm_view into PETSc.Options(), and viewFromOptions('-dm_view') then makes the object output according to that option. Here draw:tikz:data/Lshape.tex draws the mesh as TikZ code saved in a .tex file, which can be inserted into a LaTeX document directly. Note that insertString writes into the global options database, so every viewFromOptions('-dm_view') in the same process afterwards is affected by it.

# 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 网格 Viewing a Firedrake mesh with DMPlex#

Firedrake 网格底层的 DMPlex 用 mesh.topology_dm 取到. 网格分区的信息只有在多进程下才看得出来, 因此下面借助 ipyparallel 启动 2 个 MPI 进程, 之后带 %%px 的单元都在这两个进程上执行.

The DMPlex underlying a Firedrake mesh is obtained with mesh.topology_dm. The information about the mesh partitioning is only visible with more than one process, so ipyparallel is used below to start 2 MPI processes; from then on all the cells carrying %%px are executed on these two processes.

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 上的计算可以在通信进行的同时先做.

Below the same 8×8 mesh is built on 2 processes and its DMPlex is printed. The Number of 0-cells per rank line and the ones like it are each followed by two numbers, the numbers of vertices, edges and cells on the local DMPlex of each of the two processes; their sums are larger than the totals of the serial run, because the entities near the partition interface have one copy on each of the neighboring processes. Note that these are the local counts (including the halo), not the counts “owned by this process”: the sums of the counts actually owned should be exactly equal to the totals of the serial run. The three labels pyop2_core, pyop2_owned and pyop2_ghost are added by Firedrake and split the mesh points into three classes: ghost is the layer borrowed from the other processes, owned is the part owned by this process that has to be sent to the others, and core is the part owned by this process that does not have to be sent, so the computation on core can be done while the communication is going on.

%%px --block
from firedrake import *

mesh = RectangleMesh(8, 8, 1, 1)
mesh.topology_dm.view()

[stdout:0] DM Object: DM_0x30024de0_0 2 MPI processes
  type: plex
DM_0x30024de0_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 The star forest PetscSF#

PetscSF 是 PETSc 中用于进程间数据交换的数据结构. 它存储了根节点(当前进程上)和叶子节点(任意进程上)的对应关系, 可以把数据从根节点发送到叶子节点, 也可以把数据从叶子节点收集到根节点 [ZBB+22].

下面的例子在 rank 0 上建一张单位正方形的三角形网格 (9 个顶点, 8 个单元), 在顶点上定义一个每点 1 个自由度的 Section, 再把网格分发到 2 个进程: distribute() 返回的 SF 描述网格点的对应关系, 用它调用 distributeSection 可以把 Section 一并分发, 得到叶子进程上的 Section 与偏移, 最后由 createSectionSF 得到自由度层面的 SF. 函数里反复调用 distributeSection, 是为了验证其中几种调用写法 (是否预先创建 leafSection) 给出的结果一致.

PetscSF is the data structure PETSc uses to exchange data between processes. It stores the correspondence between the root nodes (on the current process) and the leaf nodes (on any process), and can send data from the roots to the leaves as well as gather data from the leaves back to the roots [ZBB+22].

The example below builds a triangular mesh of the unit square on rank 0 (9 vertices, 8 cells), defines on the vertices a Section with 1 degree of freedom per point, and then distributes the mesh over 2 processes: the SF returned by distribute() describes the correspondence of the mesh points, and calling distributeSection with it distributes the Section as well, giving the Section and the offsets on the leaf processes; createSectionSF finally gives the SF at the level of the degrees of freedom. distributeSection is called repeatedly inside the function in order to check that the several calling conventions (whether or not a leafSection is created beforehand) give the same result.

%%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, 多出来的正是分区界面上被两个进程共享的顶点.

The previous cell only defines the function; it is actually called below. The output is the degree-of-freedom SF after the distribution: rank 0 has 9 roots (1 degree of freedom on each of the 9 vertices of the original mesh; before the distribution the whole mesh is on rank 0, so the number of roots on rank 1 is 0), and each of the two processes has 6 leaves; 0 <- (0,3) means that leaf 0 of this process corresponds to root 3 on rank 0. The leaves add up to 12 rather than 9, and the extra ones are exactly the vertices on the partition interface that are shared by the two processes.

%%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 相关源码时才需要了解.

PCPatch is a preconditioner built on patches of cells: the global problem is split into many small local problems that are solved separately and then put back together, and it is commonly used as a multigrid smoother. One can use the patch types provided by PETSc (star, vanka, pardecomp) or construct one’s own.

From the dm and sf information about the solution vector supplied by the user, it constructs the correspondence between the degrees of freedom on each patch and the global degrees of freedom. What is recorded below are the internal variables of PCPatch; they only need to be understood when patches have to be constructed by hand, or when the relevant PETSc source code has to be read.

/* 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 一致)

gtol is the map from the degree-of-freedom numbering of this process to the local degrees of freedom on a single patch. The degrees of freedom of this process are a local numbering relative to the whole solution, but a global numbering relative to a single patch.

dofs is the local degree-of-freedom map on a patch (it can be seen as the cell_node_map of the patch, and agrees with the map that was passed in)

References

[HKA+21] (1,2)

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.

[LMKG16] (1,2)

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.

[ZBB+22] (1,2)

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.