盒子
盒子
文章目录
  1. KSP and PC

PETSc[04] petsc4py中的KSP和PC

2026.01.01 Zhuoyu Chen petsc  热度 ℃

KSP and PC

KSP(Krylov Subspace Methods)负责求解线性系统 $Ax=b$, 它只关心在给定算子$A$和右端项$b$的情况下，如何迭代得到近似解$x \approx A^{-1}b$，并且KSP 不关心矩阵是如何装配，因此KSP是所有PETSc的高层求解器共享的线性核心，其中：

  • SNES 在每一步非线性迭代中反复调用 KSP
  • TS 在每一个时间步中反复调用 SNES / KSP
  • DM 负责数据布局和装配回调，但不参与线性求解本身

在KSP.pyx源文件和petscksp.pxi源文件中，给出了petsc4py中的接口定义和一些封装好的API实现。在KSPType文档中，给出了Krylov求解器的类型。

创建KSP后，需要绑定线性算子$A$，以及一个预处理子$P$，如果没有给出$P$的形式，那么会使用$A$作为预处理子。代码结构如下：

1
2
3
4
5
ksp = PETSc.KSP().create(comm) # create ksp on `comm`
ksp.setOperators(A)   # set linear operator A, and preconditioner A
# or you can use `ksp.setOperators(A, P)` to set a preconditioner
ksp.setFromOptions()  # set ksp from options
ksp.solve(b, x)       # solve Ax=b

如果需要进一步设置ksp的预处理子，可以通过下面代码获取：

1
2
pc = ksp.getPC()
pc.setFromOptions()

更多的内容在PC.pyx源文件和petscpc.pxi源文件中，预处理子的类型可以从PCType文档中获取。通过setPCSide函数，可以设置预处理子的方向（左预处理子或者右预处理子），在PCSide文档里给出了详细内容：

  • PC_SIDE_DEFAULT = -1：由KSP内部决定，或者通过参数-ksp_pc_side传入
  • ``PC_LEFT = 0`：applied after the operator is applied，左预处理子 $M^{-1}Ax=M^{-1}b$
  • PC_RIGHT = 1：applied before the operator is applied，右预处理子 $AM^{-1}y=b, x=M^{-1}y$
  • PC_SYMMETRIC = 2：a portion of the preconditioner is applied before the operator and the transpose of this portion is applied after the operator is applied，对称预处理子 $M=BB^T,B^{-1}AB^{-T}y=B^{-1}b, x=B^{-T}y$

KSP的常见配置方案有：

  • 直接法（LU）： -ksp_type preonly -pc_type lu
  • 迭代法（GMRES + ILU）： -ksp_type gmres -pc_type ilu
  • CG（要求对称正定）： -ksp_type cg -pc_type jacobi

如果需要创建多个KSP，可以使用setOptionsPrefix(prefix)、appendOptionsPrefix(prefix)、getOptionsPrefix()等函数来通过prefix进行区分：

1
2
3
4
5
6
7
8
9
ksp_outer = PETSc.KSP().create(comm)
ksp_outer.setOptionsPrefix("outer_")
ksp_outer.setOperators(A)
ksp_outer.setFromOptions()

ksp_inner = PETSc.KSP().create(comm)
ksp_inner.setOptionsPrefix("inner_")
ksp_inner.setOperators(A2)
ksp_inner.setFromOptions()

之后命令行参数可以通过prefix来区分：

1
2
-outer_ksp_type gmres -outer_pc_type ilu
-inner_ksp_type cg    -inner_pc_type jacobi

为了判断KSP求解器是否成功计算出数值解，以及控制KSP的敛散条件，可以使用setTolerances(rtol, atol, divtol, max_it)来设置，并通过getTolerances来获取这四个设置参数，其中：

  • rtol表示相对误差，如果有预处理子则会返回含预处理子的相对误差，即$|r_k|<rtol \cdot |r_0|$
  • atol表示绝对误差，如果有预处理子则会返回含预处理子的相对误差，即$|r_k|<atol$
  • divtol表示发散条件，即$|r_k|>divtol \cdot |r_0|$时停止迭代
  • max_it表示最大迭代次数

还可以通过下面三个函数，来获取迭代次数、残差范数、停止原因：

  • ksp.getIterationNumber()：返回迭代次数，也常被当作 ksp.its 的来源

  • ksp.getResidualNorm()：返回残差范数，也常被当作 ksp.norm 的来源

  • ksp.getConvergedReason()：返回停止原因，即返回一个 reason 枚举

每次求解完毕后，可以通过这三个函数来输出结束原因：

1
2
3
4
5
ksp.solve(b, x)
reason = ksp.getConvergedReason()
its    = ksp.getIterationNumber()
rnorm  = ksp.getResidualNorm()
PETSc.Sys.Print(f"{reason=}, {its=}, {rnorm=}")

如果需要在求解失败时程序立刻崩溃，可以采用下面两种方案：

1
2
ksp.setErrorIfNotConverged(True)
# or set `-ksp_error_if_not_converged` as option

因为在内层求解（如某些预条件器内部又跑一个 KSP）里，DIVERGED_MAX_IT 不一定被当作错误，因为有些嵌套算法允许内层不完全收敛。还可以通过KSPNormType文档中的值来设置敛散性测试使用的范数类型：

  • KSP_NORM_DEFAULT：使用默认的KSPType
  • KSP_NORM_NONE：不使用范数计算
  • KSP_NORM_PRECONDITIONED：使用预处理残差范数
  • KSP_NORM_UNPRECONDITIONED：使用未预处理的残差范数
  • KSP_NORM_NATURAL：使用自然范数（由线性算子诱导的范数）

在KSP中，其还支持monitor进行求解过程的监控，可以通过-ksp_monitor选项打开。而setPCSide() 和 setNormType()的设置，会影响monitor的输出结果，因此需要设置好这两个类型，才能得知monitor输出的真实含义。

可以通过setMonitor来添加一个monitor，注意KSP的内部会维护一个ksp.__monitor__列表，存储了目前所有注册的monitor信息，使用方法getMonitor会返回这个列表。仿照KSP源文件中monitor(ksp, its, rnorm)的函数签名，可以写出一个类似的实例：

1
2
3
4
5
6
def my_monitor(ksp, its, rnorm):
    if its == 0: PETSc.Sys.Print("start KSP")
    PETSc.Sys.Print(f"[{its}]: rnorm={rnorm:.2e}")
ksp.setMonitor(my_monitor)
ksp.solve(b, x)
ksp.monitorCancel()

monitor中可以传入一些额外参数，比如*args, **kargs，这里可以结合KSP_Monitor和setMonitor中的实现，在ksp.__monitor__中，实际存储的是[(monitor, args, kargs), ...]这样的三元组，而在KSP_Monitor调用的时候，是执行的monitor(Ksp, toInt(its), toReal(rnm), *args, **kargs)语句，因此可以如此添加额外参数（这里的ext1, ext2也可以写成*args）：

1
2
3
4
5
6
7
def my_monitor(ksp, its, rnorm, ext1, ext2, **kargs):
# or: `def my_monitor(ksp, its, rnorm, *args, **kargs):`
    PETSc.Sys.Print(f"{ext1=}, {ext2=}, {kargs=}")
ksp.setMonitor(my_monitor, args=(1, 2), kargs={"x": "3"})
ksp.solve(b, x)
ksp.monitorCancel()
# Output: ext1=1, ext2=2, kargs={'x': '3'}

最后，总结一下比较常用的参数选项：

  • -ksp_view：打印 KSP/PC 的最终配置
  • -ksp_converged_reason：每次 solve 后给出 reason
  • -ksp_monitor：打开monitor控制
  • -ksp_monitor_true_residual：同时打印 true residual 和 preconditioned residual，避免被 “预条件残差下降很快但真实残差不动” 迷惑
  • -ksp_initial_guess_nonzero <true|false>：选择是否使用初始非0解假设