一个关于函数fvMatrix<Type>::H()函数的疑惑
-
代码如下，主要的疑惑点是这个for循环中的几行代码，我做了一些注释；首先调用函数
addBoundaryDiag(boundaryDiagCmpt, cmpt);将边界条件对矩阵的影响系数取相反数并赋值给boundaryDiagCmpt场，然后调用函数addCmptAvBoundaryDiag(boundaryDiagCmpt);，不是很理解这几行代码做了一件什么事？希望各位大佬解惑！template<class Type> Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>> Foam::fvMatrix<Type>::H() const { tmp<GeometricField<Type, fvPatchField, volMesh>> tHphi ( new GeometricField<Type, fvPatchField, volMesh> ( IOobject ( "H("+psi_.name()+')', psi_.instance(), psi_.mesh(), IOobject::NO_READ, IOobject::NO_WRITE ), psi_.mesh(), dimensions_/dimVol, extrapolatedCalculatedFvPatchScalarField::typeName ) ); GeometricField<Type, fvPatchField, volMesh>& Hphi = tHphi.ref(); // Loop over field components for (direction cmpt=0; cmpt<Type::nComponents; cmpt++) { scalarField psiCmpt(psi_.primitiveField().component(cmpt)); scalarField boundaryDiagCmpt(psi_.size(), Zero); // 定义一个全 0 的scalarField addBoundaryDiag(boundaryDiagCmpt, cmpt); // 边界条件对矩阵系数的影响（系数此时并未加到矩阵上） boundaryDiagCmpt.negate(); // 对boundaryDiagCmpt场赋值为其相反数 addCmptAvBoundaryDiag(boundaryDiagCmpt); Hphi.primitiveFieldRef().replace(cmpt, boundaryDiagCmpt*psiCmpt); } ...... return tHphi; }对应的
addCmptAvBoundaryDiag(boundaryDiagCmpt);函数和其调用的cmptAv(internalCoeffs_[patchi])函数的代码如下template<class Type> void Foam::fvMatrix<Type>::addCmptAvBoundaryDiag(scalarField& diag) const { for (label fieldi = 0; fieldi < nMatrices(); fieldi++) { const auto& bpsi = this->psi(fieldi).boundaryField(); forAll(bpsi, ptfi) { const label patchi = globalPatchID(fieldi, ptfi); if (patchi != -1) { addToInternalField ( lduAddr().patchAddr(patchi), cmptAv(internalCoeffs_[patchi]), diag ); } } } } template<class Form, class Cmpt, direction Ncmpts> inline Cmpt cmptAv ( const VectorSpace<Form, Cmpt, Ncmpts>& vs ) { return cmptSum(vs)/Ncmpts; // 好像是 (vs[0] + vs[1] + vs[2])/3 } -
考虑对流项$\nabla\cdot(\bfU\bfU)$，离散后为$\sum_f (F_f\bfU_f)$，如果是内部面，这个$f$包含了所有的面。如果是边界面，对流项$\nabla\cdot(\bfU\bfU)$离散后为$\sum_{f_{内部}} (F_{f_{内部}}\bfU_f)+F_{边界}\bfU_{边界}$。
首先对于$\sum_{f_{内部}} (F_{f_{内部}}\bfU_f)$，不管对于哪个速度分量，$F_{f_{内部}}$是固定的，所以针对每个分量离散出来的是这样：
$$
\sum_{f_{内部}} (F_{f_{内部}}\bfU_f)=
\begin{pmatrix}
\sum_{f_{内部}} (F_{f_{内部}} u_f)\\
\sum_{f_{内部}} (F_{f_{内部}} v_f)
\end{pmatrix}
$$主要是其矩阵系数是一样的。这部分系数进入的是OpenFOAM中的lower以及upper元素。
对于$F_{边界}\bfU_{边界}$。这部分系数进入的是internalCoeff以及boundaryCoeff，有多少个patch，internalCoeff里面就有多少个场数据。比如网格有3个边界，那么internalCoeff可能就是这样：
3 ( (第1组数据 对应第1个边界的数据) (第2组数据 对应第2个边界的数据) (第3组数据 对应第3个边界的数据) )boundaryCoeff同样如此。
继续考虑这个$F_{边界}\bfU_{边界}$。一般在给求解器的时候，是类似这种：
fvm::div(phi, U)，很明显phi已经是显性的值。因此对于$F_{边界}\bfU_{边界}$：-
如果给定的是固定值速度边界，$F_{边界}$已知，同时$\bfU_{边界}$已知，因此$F_{边界}\bfU_{边界}$的影响全部进入到boundaryCoeff。internalCoeff没有变化。
-
如果给的是零法向梯度，$F_{边界}$是从U的网格点赋值相等过来，也是已知。$\bfU_{边界}$与内部网格点的值相同。具体多少取决于内部网格点。因此可以将其写为$\bfU_{边界}=\bfU_{内网格点}$，这样离散后，就存在内网格点的$\bfU$，因此会改变矩阵的对角线系数internalCoeff。但对边界boundaryCoeff没有影响。
因为固定值边界对internalCoeff不影响，现在考虑零法向梯度$F_{边界}\bfU_{边界}$对矩阵的影响。假定对于某个边界面的内网格速度$u=1,v=0$。在这种情况下，$F_{边界}\bfU_{边界}$将使得u方程的矩阵以及v方程的矩阵发生区别。如果u变量与v变量的矩阵不一样，如果有代码
fvVectorMatrix UEqn，那么理论上应该区分为u矩阵以及v矩阵。但是OpenFOAM里面只有一个矩阵对角线。首先，对于一套网格系统的内部面，$\sum_f (F_f\bfU_f)$是固定的，因此u与v的矩阵系数一样。因为内部面是大多数。因此可以认为u矩阵以及v矩阵大部分元素相同。也就是这个矩阵的lower以及upper是相同的。
但是边界的印象会导致不同。那么是否需要组建3个矩阵？目前来看OpenFOAM没有这么做。为了将u分量v分量的矩阵相同化，他把边界导致的系数的影响做了个平均。加到了矩阵的对角线元素上。也就是
D()函数。深入看一下
D()函数，代码这样：void fvMatrix<Type>::addCmptAvBoundaryDiag(scalarField& diag) const { forAll(internalCoeffs_, patchI) //对每一个边界，都处理这个边界的internalCoeffs_ { addToInternalField//假如针对某个零法向梯度速度边界 ( lduAddr().patchAddr(patchI), // 这个边界的网格单元的面对应的内部网格单元编号 cmptAv(internalCoeffs_[patchI]), // 对这个patch上的internalCoeffs_，进行三个方向的平均，形成一个。 //比如原来internalCoeffs_可能是(1,2,0), (1,2,0), (1,2,0)，做完了之后成为了(1.5), (1.5), (1.5) diag//原本的对角线元素的贡献。加上考虑内部网格单元编号，与cmptAv的操作，形成一个新的internalCoeffs_的贡献 //也就是原本对角线元素 + 平均后internalCoeffs_的贡献 ); } }所以主要是这个
D()函数，里面有一个平均后的操作，这样不同分量的D()都是一套。进一步的，A()函数就是简单的D()/V，因此对于UEqn.A()，其实只有一套矩阵系数，而不是3套矩阵系数。同理，针对
H()函数，首先要进行类似D()函数的操作，比如下面的代码，就是简单的把平均后internalCoeffs_的贡献，乘以psi的值，就表示边界internalCoeffs_对H的贡献。在这个阶段，不同分量的H都是一样的。// Loop over field components for (direction cmpt=0; cmpt<Type::nComponents; cmpt++) { scalarField psiCmpt = psi_.internalField().component(cmpt); scalarField boundaryDiagCmpt(psi_.size(), 0.0); addBoundaryDiag(boundaryDiagCmpt, cmpt); boundaryDiagCmpt.negate(); addCmptAvBoundaryDiag(boundaryDiagCmpt); tHphi().internalField().replace(cmpt, boundaryDiagCmpt*psiCmpt); }下面这个代码，就是好理解的内网格+源项，就是 http://dyfluid.com/piso.html 里面的$\mathbf{HbyA}$还没除A的时候。
tHphi().internalField() += lduMatrix::H(psi_.internalField()) + source_;另外，边界不是也会影响boundaryCoeffs_么（只不过咱们这个讨论没变化）。但对于其他的情况，也要把boundaryCoeffs_对
H()的影响考虑进来。就有下面的代码：addBoundarySource(tHphi().internalField());好。那这样写，是不是就不准了？也不是。你看看这个东西影响的是$\mathbf{HbyA}$。作为一个迭代算法，中间变量多多少少问题都不大。重要的是最后的phi是不是通量守恒了。以及最终的速度分量，是不是满足动量方程。如果phi守恒，并且我们的动量方程也根本没有调用$\mathbf{HbyA}$（UEqn是怎么写的？）。所以没问题。
不知不觉竟然写了这么多。
-
-
谢谢李老师，打了这么多字！我好好消化下
-
@李东岳 老师，看了你的回复，我重新理了一下，我利用一个小栗子并结合代码来看哈
template<class Type> Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>> Foam::fvMatrix<Type>::H() const { tmp<GeometricField<Type, fvPatchField, volMesh>> tHphi ( ... ); GeometricField<Type, fvPatchField, volMesh>& Hphi = tHphi.ref(); // Loop over field components for (direction cmpt=0; cmpt<Type::nComponents; cmpt++) { scalarField psiCmpt(psi_.primitiveField().component(cmpt)); scalarField boundaryDiagCmpt(psi_.size(), Zero); // 定义一个全 0 的scalarField addBoundaryDiag(boundaryDiagCmpt, cmpt); // 边界条件对矩阵系数的影响（系数此时并未加到矩阵上） boundaryDiagCmpt.negate(); // 对boundaryDiagCmpt场赋值为其相反数 addCmptAvBoundaryDiag(boundaryDiagCmpt); Hphi.primitiveFieldRef().replace(cmpt, boundaryDiagCmpt*psiCmpt); } ...... return tHphi; }假设Patch1的internalCoeffs_是(1 2 3)；那么通过for循环中的这两行代码计算出来
addBoundaryDiag(boundaryDiagCmpt, cmpt); // 边界条件对矩阵系数的影响（系数此时并未加到矩阵上） boundaryDiagCmpt.negate(); // 对boundaryDiagCmpt场赋值为其相反数计算出来的boundaryDiagCmpt（对应for循环中cmpt分别为0,1,2）中Patch1索引对应的值为-1，-2，-3.
再通过代码addCmptAvBoundaryDiag(boundaryDiagCmpt);此时boundaryDiagCmpt（对应for循环中cmpt分别为0,1,2）中Patch1索引对应的值为1，0，-1
然后再执行代码Hphi.primitiveFieldRef().replace(cmpt, boundaryDiagCmpt*psiCmpt);是这样对吗？
-
Patch1的internalCoeffs_是(1 2 3)；
针对某个具体的patch，比如patch1，如果是速度的话，internalCoeffs_应该是
( (1 2 3) (2 3 4) (3 4 5) )每个都有3个分量，这才能跟速度3个分量对上。
我仔细处理这个问题，就是因为如果对速度离散乘fvVectorMatrix的话，应该是3套系数（边界的影响）。但是OpenFOAM明显用的是一套矩阵系数，因为UEqn.A()、UEqn.D()都是一个。原因就在你主贴的疑问里面。之前这方面我还没细看。后来发现就是internalCoeffs_给平均了。你主贴里面问了个平均的操作
cmptAv。比如上面这个internalCoeffs_明显3个分量不一样，这就需要三套系数。但是OpenFOAM给平均了：( (1+2+3)/3 (2+3+4)/3 (3+4+5)/3 )就变成了一套系数。
目前这个平均操作，我并不觉得有理论依据。就是简单平均话了。我认为的原因在上个帖子里面也讨论了。我认为就是这样（不影响收敛）。
其他你写的我也没太看懂..
-
@李东岳 老师，您说的这些我大概都清楚了（关于
cmptAv函数部分），其实我的目前更想知道的是为什么在addCmptAvBoundaryDiag(boundaryDiagCmpt);这行代码之前需要先执行addBoundaryDiag(boundaryDiagCmpt, cmpt); boundaryDiagCmpt.negate();这两行代码呢？
-
addBoundaryDiag(boundaryDiagCmpt, cmpt);
请参考：
void fvMatrix<Type>::addCmptAvBoundaryDiag(scalarField& diag) const { forAll(internalCoeffs_, patchI) //对每一个边界，都处理这个边界的internalCoeffs_ { addToInternalField//假如针对某个零法向梯度速度边界 ( lduAddr().patchAddr(patchI), // 这个边界的网格单元的面对应的内部网格单元编号 cmptAv(internalCoeffs_[patchI]), // 对这个patch上的internalCoeffs_，进行三个方向的平均，形成一个。 //比如原来internalCoeffs_可能是(1,2,0), (1,2,0), (1,2,0)，做完了之后成为了(1.5), (1.5), (1.5) diag//原本的对角线元素的贡献。加上考虑内部网格单元编号，与cmptAv的操作，形成一个新的internalCoeffs_的贡献 //也就是原本对角线元素 + 平均后internalCoeffs_的贡献 ); } }boundaryDiagCmpt.negate();
H操作符是把矩阵系数移到右边。下面的上下对角线移到右边都取负。边界导致的对角线移到右边，也要取负。
for (register label face=0; face<l.size(); face++) { Hphi[u[face]] -= Lower[face]*sf[l[face]];//取负 Hphi[l[face]] -= Upper[face]*sf[u[face]];//取负 }
