OpenFOAM中矩阵系数的管理

本文参考的链接

用有限体积法离散方程最终得到的线性方程，OpenFOAM中使用专门的矩阵求解器(PCG. PBiCG)对其进行求解。

概述

使用有限体积法离散后的方程形式
$$A\varphi =b$$

• A是系数矩阵；
• $\varphi$待求解变量；
• b是离散后的源项；

A是一个$N\times N$方阵，N表示离散后网格单元的个数；OpenFOAM使用一种特殊的方法来求解线性代数方程组，主要使用的类为

• lduMatrix (lower - diagonal - upper - matrix storage class)；
• lduAddressing (是一个索引类)；
• fvMatrix (这个更牛逼，我还没看)

lduMatrix

lduMatrix存储矩阵的系数在下面三个数组(Field就是一种高级点的数组)中

// lduMatrix.H
scalarField *lowerPtr_;
scalarField *diagPtr_;
scalarField *upperPtr_;

这些数组包含矩阵中所有的非零系数，使用以下函数得到具体的数值

// lduMatrix.H
scalarField& upper();
scalarField& lower();
scalarField& diag();

lowerPtr_和upperPtr_(指向数组)的大小等于内部面数；diagPtr_的大小等于网格单元数。

Diagonals

每个单元对应且只对应一个对角系数(都非零)；所以对角线元素为N(N为网格单元数)。

upper and lower triangles

矩阵的上下三角部分总共有$N(N-1)$，但是其中大部分为0(大型稀疏矩阵)，非零系数代表着这两个网格单元之间相互影响，这两个网格单元共享同一个内部面；OpenFOAM的处理方式为：其中一个对应下三角矩阵，另一个对应上三角矩阵。

举一个例子，下面是一个$3\times 3$矩阵

图中，黑色点为网格单元中心，蓝色点为面中心，在OpenFOAM中使用矩阵表示其中的非零系数为

在OpenFOAM中进行某些数学运算时，需要辨别两个系数的符号(比如通量，对某一单元来说为正，但是对其相邻单元来说就应该为负)；为了解决这个问题OpenFOAM关于面的owner和neighour，owenr、neighbour以及对角系数定义如下
$$A_{ij}=\{\begin{array}{rl}& owner,i>j\\ & diagonal,i=j\\ & neighbour,i<j\end{array}$$

owner和neighbour是相对某一个面来说的，内部面相邻的两个单元标号小的称为owner，大的称为neighbour。

下三角矩阵表示owner单元对主单元的影响系数，上三角矩阵表示neighbour单元对主单元的影响系数，那么lowerPtr_，upperPtr_以及diagPtr_的矩阵系数为

值得注意的是lowerPtr_和upperPtr_根据面号存储系数，而数组diagPtr_根据单元编号存储系数。在实际计算通量等物理量时需要知道面法向量，OpenFOAM中正的面法向量方向为由owner单元指向neighbour单元，如下图[^1]

lduAddressing

在OpenFOAM中，lower、diagonal以及upper系数都是有scalarField存储(本质上是list)，但是它们的索引不同。对角系数通过单元进行索引；upper和lower通过面索引，为了使得索引匹配，OpenFOAM为lower和upper提供寻址数组，以便将他们的面索引转换为相应的单元索引。这个索引数组(高级点的数组)称为lduAddressing。

lduAddressing中有两个函数，返回包含每个面的owner(lowerAddr())和neighbour(upperAddr())，例如

lowerAddr() = 1 2 1 2 3 5 4 6 5 4 7 8 
upperAddr() = 2 3 6 5 4 6 5 7 8 9 8 9

根据面编号排序，lowerAddr()为面的owner，upperAddr()指的是面的neighbour。

When do you need it?

在需要对涉及对角和非对角系数的lduMatrix的所有行执行操作时，都需要lduAddressing。

How do you use it?

分为以下几点

• 将涉及上三角矩阵的操作和涉及下三角矩阵的操作分开；
• 遍历网格中的所有面；
• 用面标号访问上下三角矩阵的系数 upper[facei], lower[facei]；
• 使用包装好的面索引访问对角线元素
  • diag[l[facei]] -= lower[facei]，l是下三角矩阵的lduAddressing；
  • diag[u[facei]] -= upper[facei]，u是上三角矩阵的lduAddressing。

Example

举个例子，lduMatrix中有个函数为negSumDiag经常被使用
$$A_p=-\sum _r{A}_r$$
源代码的实现为

// lduMatrix.C
for (label face=0; face<l.size(); face++)
{
   Diag[l[face]] -= Lower[face];
   Diag[u[face]] -= Upper[face];
}

• Diag为对角系数；
• l是下三角矩阵的lduAddressing；
• u是上三角矩阵的lduAddressing；
• Lower是下三角矩阵；
• Upper是上三角矩阵。

Diagonal, symmetric and asymmertric matrices

lduMatrix的矩阵定义中，必须定义对角线(其他的不是必须的)。如果只定义上三角或者下三角部分，那么假定矩阵是对称的(节省存储资源)。如果两部分都被定义，那么它将假定为非对称矩阵(即使这两部分相等)。

[1^]:Moukalled, F., L. Mangani, and M. Darwish. “The finite volume method in computational fluid dynamics.” An Advanced Introduction with OpenFOAM and Matlab (2016)