DMPlex 和 Mesh DMPlex and Mesh

6. DMPlex 和 Mesh DMPlex and Mesh#

DMPlex 的基本用法参见 DMPlex.

Overview

How Firedrake meshes are represented by PETSc DMPlex: creating a DMPlex, PETSc.Section, the relation between DMPlex and Firedrake Mesh/FunctionSpace, and locating geometric entities and their DoF numbering.

本章的示例大多会跑两遍: 一遍在 notebook 自己的串行内核里, 一遍通过 ipyparallel 在 2 个 MPI 进程上跑 (后者的代码单元以 %%px 开头). 这样安排是因为 DMPlex 的实体编号、core/owned/ghost 分类和 node_set 的大小只有在网格被分区之后才显出差别, 串行的结果只是并行结果的退化情形.

先启动一个 2 进程的 MPI 集群.

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'>

%%px 是 ipyparallel 提供的 cell magic, 表示这个单元的代码在集群的每个引擎 (即每个 MPI 进程) 上各执行一遍, --block 让 notebook 等所有引擎跑完再继续. 引擎有自己独立的命名空间, 所以要单独 import 一次.

%%px --block
from firedrake import *
from firedrake.petsc import PETSc
from mpi4py import MPI
import numpy as np

在串行环境中导入必要的包

from firedrake import *
from firedrake.petsc import PETSc
from mpi4py import MPI
import numpy as np

6.1. 创建 DMPlex#

DMPlex 把网格里所有维度的几何实体 (顶点、边、面、单元) 统一编号成”点” (point), 再用 cone (某个点的直接边界) 和 support (直接包含某个点的实体) 记录它们之间的覆盖关系. 整体是一张有向无环图, 也就是下面第二张图里的 Hasse 图.

下面两张 PETSc 官方插图画的是 doublet 网格: 两个共享一条边的三角形, 共 2 个单元、 5 条边、4 个顶点, 一共 11 个点. 图中同一维度的实体占据连续的一段编号: 0-1 是单元, 2-5 是顶点, 6-10 是边.

A 2D doublet mesh, two triangles sharing an edge

The Hasse diagram for our 2D doublet mesh, expressed as a DAG

6.1.1. 底层创建方法#

手工构造一个 DMPlex 的步骤是: 先用 setChart 声明一共有多少个点, 再用 setConeSize 给每个点声明它的 cone 有多大 (三角形有 3 条边, 边有 2 个顶点, 顶点没有 cone, 保持默认的 0 即可), 最后 setUp 分配存储. 下面按上图的编号构造出那张 doublet 网格.

plex = PETSc.DMPlex().create()
plex.setDimension(2)
plex.setChart(0, 11)
# plex.setConeSize(point, number of points that cover the point)
plex.setConeSize(0, 3)
plex.setConeSize(1, 3)
plex.setConeSize(6, 2)
plex.setConeSize(7, 2)
plex.setConeSize(8, 2)
plex.setConeSize(9, 2)
plex.setConeSize(10, 2)
plex = plex.setUp()  # plex.setUp() return self

接着用 setCone 逐点填入覆盖关系: 每个三角形给出它的 3 条边, 每条边给出它的 2 个顶点. symmetrize() 由 cone 反算出 support, stratify() 计算每个点的 depth (顶点为 0、边为 1、单元为 2). view() 的输出里应看到 4 个 0-cells、5 个 1-cells、 2 个 2-cells, 与上图一致.

# plex.setCone(point, [points that cover the point])
plex.setCone(0, [6, 7, 8])
plex.setCone(1, [7, 9, 10])
plex.setCone(6, [2, 3])
plex.setCone(7, [3, 4])
plex.setCone(8, [4, 2])
plex.setCone(9, [4, 5])
plex.setCone(10, [5, 3])

plex.symmetrize()
plex.stratify()

plex.view()

DM Object: 1 MPI process
  type: plex
DM_0x2dc0dcd0_0 in 2 dimensions:
  Number of 0-cells per rank: 4
  Number of 1-cells per rank: 5
  Number of 2-cells per rank: 2
Labels:
  depth: 3 strata with value/size (0 (4), 1 (5), 2 (2))
  celltype: 3 strata with value/size (0 (4), 1 (5), 3 (2))

6.1.2. 使用高级接口#

下面用 createFromCellList 重建同一张 doublet 网格: 只需给出 2 个三角形各自的顶点编号和 4 个顶点的坐标, interpolate=True 让 PETSc 自动补出中间的边. view() 报出的点数和各维实体个数与上面手工构造的完全一致 (两段输出在 DM Object: 的名字和 Labels: 的排列顺序上有出入, 不必在意).

cells = [
    [0, 1, 2],
    [1, 3, 2]
]
coords = [
    [-1, 0],
    [0, -1],
    [0, 1],
    [1, 0]
]

plex = PETSc.DMPlex().createFromCellList(
    dim=2, cells=cells, coords=coords, interpolate=True, comm=None)

plex.view()

DM Object: DM_0x2dc0dcd0_1 1 MPI process
  type: plex
DM_0x2dc0dcd0_1 in 2 dimensions:
  Number of 0-cells per rank: 4
  Number of 1-cells per rank: 5
  Number of 2-cells per rank: 2
Labels:
  celltype: 3 strata with value/size (0 (4), 1 (5), 3 (2))
  depth: 3 strata with value/size (0 (4), 1 (5), 2 (2))

6.2. `PETSc.Section`#

PETSc.Section 描述的是”每个点上放几个自由度、放在数据数组的什么位置”, 是从拓扑 (DMPlex) 到数据布局之间的桥梁. 下面建两个 Section 做对照: section1 用默认顺序, section2 额外设一个排列 (permutation), 看同样的自由度在两种顺序下 offset 有什么不同.

section1 = PETSc.Section().create()
section1.setChart(*plex.getChart())
section2 = PETSc.Section().create()
section2.setChart(*plex.getChart())

这里用 RCM (反向 Cuthill-McKee) 算一个重排. 注意方向: getOrdering 返回的是 old_to_new, 而 setPermutation 要的是反过来的 reordering[new]->[old], 因此先用一次索引赋值把它求逆.

old_to_new = plex.getOrdering(PETSc.Mat.OrderingType.RCM).indices
reordering = np.empty_like(old_to_new)  # reordering[new] -> old
reordering[old_to_new] = np.arange(old_to_new.size, dtype=old_to_new.dtype)

perm = PETSc.IS().createGeneral(reordering)
section2.setPermutation(perm)

只给顶点 (depth 为 0 的点) 各放 1 个自由度, 两个 Section 放得完全一样; setUp() 之后 offset 才被算出来.

ps, pe = plex.getDepthStratum(0)
for p in range(ps, pe):
    section1.addDof(p, 1)
    section2.addDof(p, 1)

section1.setUp()
section2.setUp()

对比两个 Section 的输出: dof 都是点 2-5 上各 1 个, 其余点为 0; 但 section1 的 offset 依次是 0、1、2、3, section2 因为设了排列变成 3、0、2、1. 可见 permutation 改变的是数据在数组里的存放次序, 不改变每个点上有多少自由度.

section1.view()
section2.view()

PetscSection Object: 1 MPI process
  type not yet set
Process 0:
  (   0) dof  0 offset   0
  (   1) dof  0 offset   0
  (   2) dof  1 offset   0
  (   3) dof  1 offset   1
  (   4) dof  1 offset   2
  (   5) dof  1 offset   3
  (   6) dof  0 offset   4
  (   7) dof  0 offset   4
  (   8) dof  0 offset   4
  (   9) dof  0 offset   4
  (  10) dof  0 offset   4
PetscSection Object: 1 MPI process
  type not yet set
Process 0:
  (   0) dof  0 offset   0
  (   1) dof  0 offset   0
  (   2) dof  1 offset   3
  (   3) dof  1 offset   0
  (   4) dof  1 offset   2
  (   5) dof  1 offset   1
  (   6) dof  0 offset   4
  (   7) dof  0 offset   4
  (   8) dof  0 offset   4
  (   9) dof  0 offset   4
  (  10) dof  0 offset   4

6.3. Firedrake 中的 Mesh#

Firedrake 的 Mesh 内部持有一个 DMPlex, 用 mesh.topology_dm 就能拿到. 但 Firedrake 并不直接沿用 DMPlex 的点编号, 而是重新编一套: 一方面按 RCM 重排以改善数据局部性, 另一方面要把本进程可见的点按 PyOP2 的三个类别分成连续的三段.

三个类别的判据是 (以 dmcommon.mark_entity_classes 的实现为准):

ghost: DMPlex point SF 的 leaf, 即本进程的影子点, 实际归别的进程所有;
owned: 本身不是 ghost, 但它所在的某个本地单元的闭包里含有 ghost 点;
core: 其余的点, 即本进程拥有、且周围一圈单元的闭包里都没有影子点.

这样分是为了让计算和通信重叠: core 上的计算不碰任何影子数据, 可以在 halo 交换还没完成时就先算; owned 上的计算必须等交换结束. core + owned 就是”本进程持有”的部分, ghost 不算在内.

顺带一提, mark_entity_classes 的 docstring 把 core/owned 描述成”是否在 send halo 里”, 和上面的实现并不一致: 实测存在被标为 core、却同时落在别的进程影子区里的点. 读源码时以实现为准.

这一节就是看这套编号是如何得来的.

mesh = RectangleMesh(4, 4, 1, 1)
plex = mesh.topology_dm
rank, size = mesh.comm.rank, mesh.comm.size

同一张网格在 2 个进程上再建一次. 引擎的命名空间和 notebook 内核是分开的, 因此 import 和建网格都要重做一遍; 这一次网格会被分成两块, 每个进程只拿到自己那块外加一圈影子区.

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

mesh = RectangleMesh(4, 4, 1, 1)
plex = mesh.topology_dm
rank, size = mesh.comm.rank, mesh.comm.size

mesh._dm_renumbering 就是上面说的那套编号, 类型是 PETSc.IS, 含义为 indices[新编号] = 原 DMPlex 点号. 每个进程只列出自己看得见的点 (含影子区), 因此两个进程各自的长度都小于串行时的总点数; 数组从前到后依次是 core、owned、ghost 三段.

%%px --block
mesh._dm_renumbering.indices

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

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

和重新编号相关的函数 (位于 firedrake/mesh.py 和 firedrake/cython/dmcommon.pyx)

MeshTopology._default_reordering: RCM 重排, _default_reordering[new]->[old]; 由 topology_dm.getOrdering(RCM) 求逆得到, 只管数据局部性, 不管并行分类.
MeshTopology._dm_renumbering: 最终使用的实体编号, _dm_renumbering[new]->[old]; 在 RCM 的基础上再按 core/owned/ghost 分块.
MeshTopology._renumber_entities: 上面两步的入口; reorder=False 时跳过 RCM, 直接按 DMPlex 原始顺序遍历.
dmcommon.mark_entity_classes: 检查 DMPlex 的 point SF, 给每个点打上 pyop2_core、pyop2_owned、pyop2_ghost 之一的 label.
dmcommon.get_entity_classes: 把这三个 label 统计成每个维度的累积分界数组, 即 mesh._entity_classes.
dmcommon.plex_renumbering: 真正生成排列的地方, 返回的 PETSc.IS 满足 indices[PyOP2 实体号] = DMPlex 点号, 与第 2 条同向. 它只按 RCM 顺序把单元遍历一遍, 对每个单元闭包里尚未编号的点查它属于哪一类, 再用 core/owned/ghost 三个各自独立的写游标写进对应的一段 —— 三段因此各自连续, 并不是把网格遍历三遍.

下面把 mesh._default_reordering 和手工算出来的 RCM 排列并排打印, 两者应当完全一致, 说明 _default_reordering 就是对本进程的 DMPlex 调用 getOrdering(RCM) 再求逆. 它和 _dm_renumbering 不同 —— 后者在此基础上还按 core/owned/ghost 分了块.

%%px --block
# https://petsc.org/release/manualpages/DMPlex/DMPlexGetOrdering
old_to_new = plex.getOrdering(PETSc.Mat.OrderingType.RCM).indices
reordering = np.empty_like(old_to_new)
reordering[old_to_new] = np.arange(old_to_new.size, dtype=old_to_new.dtype)
mesh._default_reordering, reordering

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

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

这个单元先 view() 出 DMPlex. 并行时它只打印一份: Number of ...-cells per rank: 那几行按进程逐个列出实体个数, 而紧随其后 Labels: 里的 pyop2_core: 1 strata with value/size (1 (...)) 等只是 0 号进程的数字, 不是全局总量, 也不是别的进程的数字. 这三个 pyop2_* label 就是 mark_entity_classes 打上去的.

接着打印 mesh._entity_classes, 它是逐进程的. 每行对应一个维度的实体 (第 0 行是顶点, 第 1 行是边, 第 2 行是单元), 三个数是累积分界, 即 [core 数, core+owned 数, 三类总数]. 串行时没有影子区, 三个数相等.

%%px --block
mesh.topology_dm.view()
entity_classes = mesh._entity_classes
ec_msg = ', '.join([f'{i}-cells: {ec}' for i, ec in enumerate(entity_classes)])
PETSc.Sys.syncPrint(f'[{rank}/{size}] mesh._entity_classes (core, owned, ghost): {ec_msg}')
PETSc.Sys.syncFlush()

[stdout:0] DM Object: DM_0x1ab381b0_0 2 MPI processes
  type: plex
DM_0x1ab381b0_0 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))
[0/2] mesh._entity_classes (core, owned, ghost): 0-cells: [ 5 10 19], 1-cells: [13 26 38], 2-cells: [ 8 16 20]
[1/2] mesh._entity_classes (core, owned, ghost): 0-cells: [10 15 19], 1-cells: [26 30 38], 2-cells: [16 16 20]

计算 node_set 相关的函数 (位于 firedrake/functionspacedata.py、firedrake/mesh.py、 firedrake/cython/dmcommon.pyx 和 firedrake/halo.py)

functionspacedata.get_node_set: 汇合下面几项, 组装出 PyOP2 的节点集合 node_set.
functionspacedata.get_global_numbering: 取得 (必要时新建) 描述节点布局的 PETSc.Section, 也就是后面用到的 V.global_numbering. 注意它返回的是 (Section, 受约束节点数) 二元组, V.global_numbering 是其中的第一项.
AbstractMeshTopology.create_section: 网格一侧的入口, 转手调用下一项, 同样返回上述二元组 (它自己的 docstring 只写了 Section, 不准).
dmcommon.create_section: 按 nodes_per_entity 给每个 DMPlex 点分配自由度个数和 offset, 生成上述 Section.
firedrake.Halo: 用 DMPlex 的 point SF 和这个 Section 构造影子区的通信结构.
make_dofs_per_plex_entity: MeshTopology 的方法, 把有限元的 entity_dofs 压成 “每个维度的实体各带几个节点”, 即 nodes_per_entity.

计算 node_set 相关的概念

entity_dofs (FIAT FIAT/finite_element.py): 参考单元上”哪个实体带哪几个局部自由度”的映射, 通过 V.finat_element.entity_dofs() 取到.
nodes_per_entity (firedrake/functionspacedata.py): 把 entity_dofs 压成按维度计的节点数元组, 是 create_section 的直接输入.

把这两个概念打印出来对照. CG1 只在顶点上有自由度, 所以 nodes_per_entity 是 (1, 0, 0); CG2 在顶点和边上各有一个, 是 (1, 1, 0). 元组的三个位置依次对应顶点、边、单元.

V1 = FunctionSpace(mesh, 'CG', 1)
V2 = FunctionSpace(mesh, 'CG', 2)
PETSc.Sys.Print(f'V1 entity_dofs: {V1.finat_element.entity_dofs()}')
PETSc.Sys.Print(f'V2 entity_dofs: {V2.finat_element.entity_dofs()}')
nodes_per_entity_V1 = tuple(mesh.make_dofs_per_plex_entity(V1.finat_element.entity_dofs()))
nodes_per_entity_V2 = tuple(mesh.make_dofs_per_plex_entity(V2.finat_element.entity_dofs()))
PETSc.Sys.Print(f'V1 nodes_per_entity: {nodes_per_entity_V1}')
PETSc.Sys.Print(f'V2 nodes_per_entity: {nodes_per_entity_V2}')

V1 entity_dofs: {0: {0: [0], 1: [1], 2: [2]}, 1: {2: [], 1: [], 0: []}, 2: {0: []}}
V2 entity_dofs: {0: {0: [0], 1: [3], 2: [5]}, 1: {2: [1], 1: [2], 0: [4]}, 2: {0: []}}
V1 nodes_per_entity: (1, 0, 0)
V2 nodes_per_entity: (1, 1, 0)

同一段代码在 2 个进程上跑一遍. entity_dofs 和 nodes_per_entity 只取决于有限元本身, 与网格如何分区无关, 所以并行的输出和串行完全一样; 这里用的 PETSc.Sys.Print 只在 0 号进程打印, 因此每行仍只出现一次.

%%px --block
V1 = FunctionSpace(mesh, 'CG', 1)
V2 = FunctionSpace(mesh, 'CG', 2)
PETSc.Sys.Print(f'V1 entity_dofs: {V1.finat_element.entity_dofs()}')
PETSc.Sys.Print(f'V2 entity_dofs: {V2.finat_element.entity_dofs()}')
nodes_per_entity_V1 = tuple(mesh.make_dofs_per_plex_entity(V1.finat_element.entity_dofs()))
nodes_per_entity_V2 = tuple(mesh.make_dofs_per_plex_entity(V2.finat_element.entity_dofs()))
PETSc.Sys.Print(f'V1 nodes_per_entity: {nodes_per_entity_V1}')
PETSc.Sys.Print(f'V2 nodes_per_entity: {nodes_per_entity_V2}')

[stdout:0] V1 entity_dofs: {0: {0: [0], 1: [1], 2: [2]}, 1: {2: [], 1: [], 0: []}, 2: {0: []}}
V2 entity_dofs: {0: {0: [0], 1: [3], 2: [5]}, 1: {2: [1], 1: [2], 0: [4]}, 2: {0: []}}
V1 nodes_per_entity: (1, 0, 0)
V2 nodes_per_entity: (1, 1, 0)

章首说过 node_set 的大小要分区之后才有意义, 这里补上. node_set 是 PyOP2 用来描述 “本进程有哪些节点”的集合, 它的分类正是从 mesh._entity_classes 来的. 下面打印两个容易混淆的尺寸: node_set.size 是本进程持有的节点数, 也就是 core + owned; 而 V.node_count 还把影子区算在内, 两者不要混用. 两个进程的 size 加起来等于全局节点数; 对 CG1 而言, size 恰好等于上面 _entity_classes 第 0 行的第 2 个数.

%%px --block
PETSc.Sys.syncPrint(f'[{rank}/{size}] V1 node_set.size (owned): {V1.node_set.size}, '
                    f'node_count (with halo): {V1.node_count}')
PETSc.Sys.syncFlush()

[stdout:0] [0/2] V1 node_set.size (owned): 10, node_count (with halo): 19
[1/2] V1 node_set.size (owned): 15, node_count (with halo): 19

6.4. `FunctionSpace`#

把前面几节串起来: 一个 FunctionSpace 的自由度布局由 FunctionSpaceData (即 V._shared_data) 缓存, 其中 global_numbering 就是上一节那种 PETSc.Section, node_set 是 PyOP2 的节点集合; dof_dset 在其上再包一层 DMShell, 供 PETSc 一侧使用. 这几个对象的引用关系如下.

下图使用 Sphinx / Jupyter Book 图形插件 plantuml 制作.

skinparam monochrome true
skinparam defaultFontSize 14
skinparam defaultFontName Aapex

class "FunctionSpace" as fs {
*dm: DMShell, cached_property
*_shared_data: FunctionSpaceData
*node_set (_shared_data.node_set): Set
*dof_dset: DataSet
*cell_node_list
*finat_element
--
*_dm(): { dm = dof_dset.dm; attach_hooks(dm, level, sf, section); }
}

class "FunctionSpaceData" as fsd {
*entity_node_lists: dict
*node_set: Set
*global_numbering: PETSc.Section
}

class "DataSet" as dset {
*dm: DMShell, cached_property
*layout_vec: Vec, cached_property
--
*dm()
}
note left of dset::dm()
dm = PETSc.DMShell().create(comm=self.comm)
dm.setGlobalVector(self.layout_vec)
end note

class "DMShell" as dm {
* PetscSF
* PetscSection
--
}

fs o-- dset
fs o-- fsd
dset *-- dm — Fig. 6.1 class FunctionSpace#

6.5. 示例: 寻找特定的几何实体以及节点编号#

前面几节讲的是机制, 这一节解决一个具体问题: 已知 gmsh 打的物理标签, 如何在 DMPlex 里定位到对应的边或面, 再把 DMPlex 的点号翻译成 Firedrake 的节点编号, 从而能直接去 Function 的数据数组里读取这些点上的值. 翻译这一步靠的就是 V.global_numbering 这个 PETSc.Section: getOffset(点号) 给出该点上第一个节点在本进程数据数组中的下标.

并行时还多一件事: 同一个实体可能同时出现在几个进程的影子区里, 因此要用 pyop2_core 和 pyop2_owned 两个 label 把它筛给唯一的持有者, 否则会被重复处理.

6.5.1. 二维示例#

使用 global_numbering 寻找某条边界上端点, 以及在这条边界上与它相邻的点 (2D mesh)

网格 gmsh/rectangle.msh 的四条边界在 gmsh/rectangle.geo 里各打了一个物理标签: 1 是下边界, 2 是上边界, 3 是左边界, 4 是右边界. 载入后这些标签落在 DMPlex 的 Face Sets label 上 —— 在下面 topology_dm.view() 的输出里可以看到 Face Sets: 4 strata.

from firedrake import *
from py.intro_utils import triplot
import matplotlib.pyplot as plt
import numpy as np

rectangle = Mesh("gmsh/rectangle.msh")
fig, axes = plt.subplots(figsize=[4, 3])
triplot(rectangle, axes=axes)
axes.set_aspect('equal')
rectangle.topology_dm.view()

DM Object: firedrake_default_topology 1 MPI process
  type: plex
firedrake_default_topology in 2 dimensions:
  Number of 0-cells per rank: 98
  Number of 1-cells per rank: 259
  Number of 2-cells per rank: 162
Labels:
  celltype: 3 strata with value/size (0 (98), 1 (259), 3 (162))
  depth: 3 strata with value/size (0 (98), 1 (259), 2 (162))
  Cell Sets: 1 strata with value/size (1 (162))
  Face Sets: 4 strata with value/size (1 (17), 2 (17), 3 (17), 4 (17))
  exterior_facets: 1 strata with value/size (1 (64))
  interior_facets: 1 strata with value/size (1 (325))
  pyop2_core: 1 strata with value/size (1 (519))
  pyop2_owned: 0 strata with value/size ()
  pyop2_ghost: 0 strata with value/size ()

../_images/dba6bf00836f8f1b866b72ee2b19e41255851493af2761cb4719e4e9da163bfd.png

并行版本. 网格被分成两块, 每个进程只画自己分到的那部分; 这里额外固定了坐标轴范围, 方便和串行的图对照.

%%px --block
from firedrake import *
from py.intro_utils import triplot
import matplotlib.pyplot as plt
import numpy as np

rectangle = Mesh("gmsh/rectangle.msh")
fig, axes = plt.subplots(figsize=[4, 3])
triplot(rectangle, axes=axes)
axes.set_aspect('equal')
axes.set_xlim([-0.1, 1.1])
axes.set_ylim([-0.1, 1.1])
rectangle.topology_dm.view()

[stdout:0] DM Object: firedrake_default_topology 2 MPI processes
  type: plex
firedrake_default_topology in 2 dimensions:
  Number of 0-cells per rank: 61 60
  Number of 1-cells per rank: 150 149
  Number of 2-cells per rank: 90 90
Labels:
  depth: 3 strata with value/size (0 (61), 1 (150), 2 (90))
  celltype: 3 strata with value/size (0 (61), 1 (150), 3 (90))
  Cell Sets: 1 strata with value/size (1 (90))
  Face Sets: 3 strata with value/size (2 (17), 3 (11), 4 (9))
  exterior_facets: 1 strata with value/size (1 (35))
  interior_facets: 1 strata with value/size (1 (194))
  pyop2_core: 1 strata with value/size (1 (190))
  pyop2_owned: 1 strata with value/size (1 (60))
  pyop2_ghost: 1 strata with value/size (1 (51))

[output:1]

../_images/ff2a1ed9a8a2a8d1c424c5b931f2ccbe1a707674e6a8e2b3f596d8548c55d86a.png

[output:0]

../_images/e494f0115004b1b96a3135aa21a977e27a5d8a45f3f5fe9ed2ac6b2e3d32f67a.png

下面这个函数做三件事. 第一, 用 Face Sets label 分别取出两条边界上的点 —— Firedrake 从 gmsh 载入时会把物理曲线上的边连同它的端点一起打上标签, 因此两个集合的交集 (np.intersect1d) 恰好是这两条边界的公共角点. 第二, 在这个角点的 support (所有以它为端点的边) 里挑一条既属于目标边界、又是本进程 core 或 owned 的边, 这条边的 cone 给出另一个端点. 第三, 用 V.global_numbering.getOffset 把两个 DMPlex 点号换成 Firedrake 的节点号, 拿去 coordinates.dat 里取坐标.

串行这一次取的是 tag 1 (下边界) 和 tag 3 (左边界), 得到的是角点 (0, 0) 以及它在下边界上的邻点. 注意取坐标的 rectangle.coordinates.dat.data_ro_with_halos 写在 if 外面: 它是集合操作, 只让一部分进程调用会把程序挂住.

还有一个容易被忽略的前提: getOffset 给出的是该函数空间自己的节点下标, 拿它去索引坐标数组, 只在 V 是 CG1 时才成立 —— CG1 的节点与网格顶点一一对应, 编号也一致. 换成 CG2 就不再成立: 实测同一个顶点在 CG2 下的 offset 是 134, 而这张网格的坐标数组只有 98 个元素, 直接索引会 IndexError. 要按几何位置取值, 应当用坐标场自己的函数空间的编号.

def get_interface_element_with_contact_point(mesh, V, interface_tag, adj_line_tag):
    dm = mesh.topology_dm
    edge_label = dm.getLabel("Face Sets")  # 2D Face is Edge
    core_label = dm.getLabel("pyop2_core")
    owned_label = dm.getLabel("pyop2_owned")

    edge_label_values = edge_label.getValueIS().indices
    interface_indices = []
    if interface_tag in edge_label_values:
        interface_indices = edge_label.getStratumIS(interface_tag).indices

    adj_line = []
    if adj_line_tag in edge_label_values:
        adj_line = edge_label.getStratumIS(adj_line_tag).indices

    points = np.intersect1d(interface_indices, adj_line)

    plex_element = []
    if len(points) > 0:
        point = points[0]
        support = dm.getSupport(point)
        for edge in support:
            if edge_label.getValue(edge) == interface_tag and \
                (core_label.getValue(edge) == 1 or owned_label.getValue(edge) == 1):
                cone = dm.getCone(edge)
                adj_point = cone[1] if cone[0] == point else cone[0]
                plex_element = [point, adj_point]
                break
            
    local_section = V.global_numbering  # global_numbering is a local section
    element = [local_section.getOffset(_) for _ in plex_element]

    return element

# Face Sets tags of gmsh/rectangle.msh: 1 lower, 2 upper, 3 left, 4 right
V = FunctionSpace(rectangle, 'CG', 1)
element = get_interface_element_with_contact_point(rectangle, V, interface_tag=1, adj_line_tag=3)
coords_data = rectangle.coordinates.dat.data_ro_with_halos  # This must be outside the if condition (mpi collective)

rank, size = rectangle.comm.rank, rectangle.comm.size
if len(element) > 0:
    coords = [coords_data[_] for _ in element]
    PETSc.Sys.syncPrint(f"[{rank}/{size}] node {element[0]}: {coords[0]}), node {element[1]}: {coords[1]}")
PETSc.Sys.syncFlush()

[0/1] node 41: [0. 0.]), node 47: [0.125 0.   ]

并行版本换成 tag 2 (上边界) 和 tag 4 (右边界), 找的是角点 (1, 1). 2 进程下只有一个进程打印出一行, 另一个进程的 element 是空的.

原因不在 core/owned 判断, 而在分区: 这一次的分区把上边界整块分给了其中一个进程, 另一个进程的 Face Sets 里压根没有 tag 2 这个值, interface_indices 是空的, np.intersect1d 直接得到空集, 于是 if len(points) > 0 之下的代码一次都没执行. 分区是会变的 —— 改进程数, 甚至只是改变本次会话中建网格的先后顺序, 有输出的就可能换成另一个进程. 因此不要把”哪个进程有输出”当成固定结论.

那么 core_label / owned_label 那个判断起什么作用? 它解决的是另一个问题: 同一条边可能同时出现在几个进程的影子区里, 有了这个条件才能保证它只被唯一的持有者处理, 结果不会随进程数翻倍. 本例进程数少, 角点周围那几条边恰好全是 core, 因此这个条件一条也没筛掉 —— 在这个单元里真正起筛选作用的是同一行的 edge_label.getValue(edge) == interface_tag, 它把角点 support 里的内部边和另一条边界上的边排除掉了.

%%px --block
def get_interface_element_with_contact_point(mesh, V, interface_tag, adj_line_tag):
    dm = mesh.topology_dm
    edge_label = dm.getLabel("Face Sets")  # 2D Face is Edge
    core_label = dm.getLabel("pyop2_core")
    owned_label = dm.getLabel("pyop2_owned")

    edge_label_values = edge_label.getValueIS().indices
    interface_indices = []
    if interface_tag in edge_label_values:
        interface_indices = edge_label.getStratumIS(interface_tag).indices

    adj_line = []
    if adj_line_tag in edge_label_values:
        adj_line = edge_label.getStratumIS(adj_line_tag).indices

    points = np.intersect1d(interface_indices, adj_line)

    plex_element = []
    if len(points) > 0:
        point = points[0]
        support = dm.getSupport(point)
        for edge in support:
            if edge_label.getValue(edge) == interface_tag and \
                (core_label.getValue(edge) == 1 or owned_label.getValue(edge) == 1):
                cone = dm.getCone(edge)
                adj_point = cone[1] if cone[0] == point else cone[0]
                plex_element = [point, adj_point]
                break
            
    local_section = V.global_numbering  # global_numbering is a local section
    element = [local_section.getOffset(_) for _ in plex_element]

    return element

# Face Sets tags of gmsh/rectangle.msh: 1 lower, 2 upper, 3 left, 4 right

V = FunctionSpace(rectangle, 'CG', 1)
element = get_interface_element_with_contact_point(rectangle, V, interface_tag=2, adj_line_tag=4)
coords_data = rectangle.coordinates.dat.data_ro_with_halos  # This must be outside the if condition (mpi collective)

rank, size = rectangle.comm.rank, rectangle.comm.size
if len(element) > 0:
    coords = [coords_data[_] for _ in element]
    PETSc.Sys.syncPrint(f"[{rank}/{size}] node {element[0]}: {coords[0]}), node {element[1]}: {coords[1]}")
PETSc.Sys.syncFlush()

[stdout:0] [0/2] node 32: [1. 1.]), node 29: [0.875 1.   ]

6.5.2. 三维示例#

使用 global_numbering 寻找界面上与接触线相邻的三角形 (3D mesh)

本示例网格文件 cylinder.msh 由几何文件 cylinder.geo 生成, 如图:

这是上下两块拼起来的圆柱, 交界面在 z = 0.25 处. 交界面的物理标签是 3 (interface), 它的圆周边界线标签是 1 (contact_line). 载入后, 一维的物理曲线落在 DMPlex 的 Edge Sets label 上, 二维的物理曲面落在 Face Sets 上.

import matplotlib.pyplot as plt
import numpy as np
from py.intro_utils import triplot

cylinder = Mesh("gmsh/cylinder.msh")
fig, axes = plt.subplots(figsize=[4, 3], subplot_kw={'projection': '3d'})
triplot(cylinder, axes=axes)
axes.set_aspect('equal')
cylinder.topology_dm.view()

DM Object: firedrake_default_topology 1 MPI process
  type: plex
firedrake_default_topology in 3 dimensions:
  Number of 0-cells per rank: 157
  Number of 1-cells per rank: 774
  Number of 2-cells per rank: 1106
  Number of 3-cells per rank: 488
Labels:
  celltype: 4 strata with value/size (0 (157), 1 (774), 3 (1106), 6 (488))
  depth: 4 strata with value/size (0 (157), 1 (774), 2 (1106), 3 (488))
  Cell Sets: 2 strata with value/size (1 (244), 2 (244))
  Face Sets: 5 strata with value/size (3 (209), 4 (209), 5 (209), 6 (230), 7 (230))
  Edge Sets: 1 strata with value/size (1 (16))
  exterior_facets: 1 strata with value/size (1 (782))
  interior_facets: 1 strata with value/size (1 (1745))
  pyop2_core: 1 strata with value/size (1 (2525))
  pyop2_owned: 0 strata with value/size ()
  pyop2_ghost: 0 strata with value/size ()

../_images/fe1f76ab0650b8507347256b6331cb1a3c4ccb4ef2c3997ac2d8c52241436b1d.png

并行版本. 每个进程只持有网格的一部分, view() 输出的 Number of 0-cells per rank: ... 会按进程逐个列出本进程可见的实体个数.

%%px --block
import matplotlib.pyplot as plt
import numpy as np
from py.intro_utils import triplot

cylinder = Mesh("gmsh/cylinder.msh")
fig, axes = plt.subplots(figsize=[4, 3], subplot_kw={'projection': '3d'})
triplot(cylinder, axes=axes)
axes.set_aspect('equal')
cylinder.topology_dm.view()

[stdout:0] DM Object: firedrake_default_topology 2 MPI processes
  type: plex
firedrake_default_topology in 3 dimensions:
  Number of 0-cells per rank: 104 109
  Number of 1-cells per rank: 468 476
  Number of 2-cells per rank: 637 640
  Number of 3-cells per rank: 272 272
Labels:
  depth: 4 strata with value/size (0 (104), 1 (468), 2 (637), 3 (272))
  celltype: 4 strata with value/size (0 (104), 1 (468), 3 (637), 6 (272))
  Cell Sets: 2 strata with value/size (1 (136), 2 (136))
  Face Sets: 5 strata with value/size (3 (125), 4 (129), 5 (129), 6 (129), 7 (129))
  Edge Sets: 1 strata with value/size (1 (9))
  exterior_facets: 1 strata with value/size (1 (459))
  interior_facets: 1 strata with value/size (1 (1047))
  pyop2_core: 1 strata with value/size (1 (742))
  pyop2_owned: 1 strata with value/size (1 (470))
  pyop2_ghost: 1 strata with value/size (1 (269))

[output:1]

../_images/28dde51f7697a01a9e27dc10775fcd805d4a584f0bbcbc23255238f2da1fcf88.png

[output:0]

../_images/23aa64d6bd0cebb974c4ff03c9f856458ca636aeba7c8fef67d84d52237a0852.png

三维的做法和二维同构, 只是各维度升一级. 对接触线上属于本进程的每一段 (Edge Sets 里 tag 为 1 的那些边), 在它的 support (所有含这条边的面) 里找那个既属于界面 (Face Sets 的 tag 3)、又是本进程 core 或 owned 的三角形; 这个三角形的三个顶点就是该段的两个端点, 加上用 np.setdiff1d 从面的闭包里挑出来的第三个顶点. 最后同样用 global_numbering 把点号换成节点号, 得到一张形状为 (三角形数, 3) 的 cell_node_map. 中间那句 assert ... getDof(_) > 0 是在确认这些顶点上确实有 CG1 的自由度.

并行时这段代码有个隐患值得指出: 本进程的 contact_line 里可能混进影子段, 它们找不到 core 或 owned 的界面三角形, 于是 faces 会比 contact_line 短, 而下一步的 zip(contact_line, faces) 会静默截断, 不报任何错. 实测 2 进程下就有一个进程是 9 段对 8 个面, 目前没出乱子只是因为那个没配上的影子段恰好排在最后一位. 稳妥的写法是在同一个循环里把 (段, 面) 成对收集, 而不是先分别建两个列表再 zip.

def get_interface_element_include_contact_line(mesh, V, interface_tag, contact_line_tag):
    dm = mesh.topology_dm

    edge_label = dm.getLabel("Edge Sets")
    face_label = dm.getLabel("Face Sets")
    core_label = dm.getLabel("pyop2_core")
    owned_label = dm.getLabel("pyop2_owned")

    edge_label_values = edge_label.getValueIS().indices
    contact_line = []
    if contact_line_tag in edge_label_values:
        contact_line = edge_label.getStratumIS(contact_line_tag).indices

    faces = []
    for seg in contact_line:
        for face in dm.getSupport(seg):
            if face_label.getValue(face) == interface_tag and \
                (core_label.getValue(face) == 1 or owned_label.getValue(face) == 1):
                faces.append(int(face))
                break

    plex_cell_node_map = np.zeros((len(faces), 3), dtype=np.int32)
    for i, (seg, face) in enumerate(zip(contact_line, faces)):
        seg_nodes = dm.getCone(seg)
        plex_cell_node_map[i, :2] = seg_nodes 
        plex_cell_node_map[i, 2:] = np.setdiff1d(
            np.unique(np.array([dm.getCone(_) for _ in dm.getCone(face)]).flatten()),
            seg_nodes)

    local_section = V.global_numbering  # global_numbering is a local section
    cell_node_map = np.zeros_like(plex_cell_node_map)
    for i, cell in enumerate(plex_cell_node_map):
        assert np.all(np.array([local_section.getDof(_) for _ in cell]) > 0)
        cell_node_map[i, :] = [local_section.getOffset(_) for _ in cell]

    return cell_node_map

CONTACT_LINE = 1
INTERFACE = 3
V = FunctionSpace(cylinder, 'CG', 1)
cell_node_map = get_interface_element_include_contact_line(cylinder, V, INTERFACE, CONTACT_LINE)

并行版本. 这张网格的接触线一共 16 段, 这一次的分区下由两个进程分别持有其中一部分; 注意每个进程的 contact_line 里还可能带上邻居的影子段, 因此两边的段数之和会比 16 大 (实测是 8 和 9). 两个进程的 cell_node_map 合起来才是完整的一圈.

%%px --block
def get_interface_element_include_contact_line(mesh, V, interface_tag, contact_line_tag):
    dm = mesh.topology_dm

    edge_label = dm.getLabel("Edge Sets")
    face_label = dm.getLabel("Face Sets")
    core_label = dm.getLabel("pyop2_core")
    owned_label = dm.getLabel("pyop2_owned")

    edge_label_values = edge_label.getValueIS().indices
    contact_line = []
    if contact_line_tag in edge_label_values:
        contact_line = edge_label.getStratumIS(contact_line_tag).indices

    faces = []
    for seg in contact_line:
        for face in dm.getSupport(seg):
            if face_label.getValue(face) == interface_tag and \
                (core_label.getValue(face) == 1 or owned_label.getValue(face) == 1):
                faces.append(int(face))
                break

    plex_cell_node_map = np.zeros((len(faces), 3), dtype=np.int32)
    for i, (seg, face) in enumerate(zip(contact_line, faces)):
        seg_nodes = dm.getCone(seg)
        plex_cell_node_map[i, :2] = seg_nodes  # set the seg nodes first
        plex_cell_node_map[i, 2:] = np.setdiff1d(
            np.unique(np.array([dm.getCone(_) for _ in dm.getCone(face)]).flatten()),
            seg_nodes)

    local_section = V.global_numbering  # global_numbering is a local section
    cell_node_map = np.zeros_like(plex_cell_node_map)
    for i, cell in enumerate(plex_cell_node_map):
        assert np.all(np.array([local_section.getDof(_) for _ in cell]) > 0)
        cell_node_map[i, :] = [local_section.getOffset(_) for _ in cell]

    return cell_node_map

CONTACT_LINE = 1
INTERFACE = 3
V = FunctionSpace(cylinder, 'CG', 1)
cell_node_map = get_interface_element_include_contact_line(cylinder, V, INTERFACE, CONTACT_LINE)

把找到的三角形画出来验证: 先断言它们的顶点 z 坐标都等于 0.25, 也就是确实落在界面上; 再在 xy 平面里画出这些三角形, 用黑色虚线标出其中落在接触线上的那条边. 看到的应是沿着圆周排开的一圈三角形.

# plot the triangle to check if they are on the interface
coords = cylinder.coordinates.dat.data_ro_with_halos
if len(cell_node_map) > 0:
    assert np.allclose(coords[:, 2][cell_node_map], 0.25)
    fig, axes = plt.subplots(figsize=[4, 3])
    c = axes.triplot(coords[:, 0], coords[:, 1], triangles=cell_node_map)
    lines = [[(coords[_, 0], coords[_, 1]) for _ in __[:2] ] for __ in cell_node_map]
    from matplotlib.collections import LineCollection
    line_collection = LineCollection(lines, colors='k', linestyles=':')
    axes.add_collection(line_collection)
    axes.set_xlim([-0.52, 0.52])
    axes.set_ylim([-0.52, 0.52])
    axes.set_aspect("equal")
    axes.grid("on")

../_images/fd3dc116f96f482a48ed531e10a2a0e52362efc1e9b60a20f4cbcfabc0a2a287.png

并行版本, 每个进程画自己持有的那一部分.

%%px --block
# plot the triangle to check if they are on the interface
coords = cylinder.coordinates.dat.data_ro_with_halos
if len(cell_node_map) > 0:
    assert np.allclose(coords[:, 2][cell_node_map], 0.25)
    fig, axes = plt.subplots(figsize=[4, 3])
    c = axes.triplot(coords[:, 0], coords[:, 1], triangles=cell_node_map)
    lines = [[(coords[_, 0], coords[_, 1]) for _ in __[:2] ] for __ in cell_node_map]
    from matplotlib.collections import LineCollection
    line_collection = LineCollection(lines, colors='k', linestyles=':')
    axes.add_collection(line_collection)
    axes.set_xlim([-0.52, 0.52])
    axes.set_ylim([-0.52, 0.52])
    axes.set_aspect("equal")
    axes.grid("on")

[output:1]

../_images/0ab4aa73b6aede69218340768bce83ccb6b90125fe10497855aad3c9afb67139.png

[output:0]

../_images/0898430e1bca0dadb319e0553c068c204b75f2bb827bc884f388c75efcfb7e71.png