关于negSumDiag()的一个问题
-
在看关于negSumDiag()的代码的时候遇到一个问题。diag上面的项是由该“行”上面的的off-diag的项相加后取负值得到。为什么下面的代码感觉像是把该“列”的数相加了呢？
template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::negSumDiag()
{
const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
Field<DType>& Diag = diag();const unallocLabelList& l = lduAddr().lowerAddr(); const unallocLabelList& u = lduAddr().upperAddr(); for (label face=0; face<l.size(); face++) { Diag[l[face]] -= Lower[face]; Diag[u[face]] -= Upper[face]; }}
以这个网格为例
lowerAdd 为 (0 1 2 1 2 3 4 3 )
upperAdd 为 (2 2 3 6 4 5 6 5 )矩阵的形式为
那么Diag[0] 应该是由 Upper[1]得到。这个和上面的代码不符。如果矩阵是对称的，这个计算没有问题。而如果是 非对称的，是不是就有问题了？
谢谢大家！
-
为了进一步测试，建立一个3X3的网格，并创建一个fvMatrix。
fvVectorMatrix UEqn_lap ( fvm::laplacian(nu, U) );网格编号如下：
|6|7|8|
|3|4|5|
|0|1|2|输出 upper 和lower 的address，没有问题
Info<< "UEqn_lap upper add before sum " << UEqn_lap.lduAddr().upperAddr() << nl << endl; Info<< "UEqn_lap lower add before sum " << UEqn_lap.lduAddr().lowerAddr() << nl << endl;UEqn_lap upper add before sum
12
(
1
3
2
4
5
4
6
5
7
8
7
8
)UEqn_lap lower add before sum
12
(
0
0
1
1
2
3
3
4
4
5
6
7
)下面把upper lower 和diag重新赋值 并输出
UEqn_lap.upper() = 100; UEqn_lap.diag() = 0; UEqn_lap.lower() = 1; Info<< "UEqn_lap diag after change " << UEqn_lap.diag() << nl << endl; Info<< "UEqn_lap upper after change " << UEqn_lap.upper() << nl << endl; Info<< "UEqn_lap lower after change " << UEqn_lap.lower() << nl << endl;输出为：
UEqn_lap diag after change 9{0}
UEqn_lap upper after change 12{100}
UEqn_lap lower after change 12{1}
调用
UEqn_lap.negSumDiag(); Info<< "UEqn_lap diag after negSumDiag " << UEqn_lap.diag() << nl << endl;输出结果为：
UEqn_lap diag after negSumDiag 9(-2 -102 -101 -102 -202 -201 -101 -201 -200)
我特意把upper里面的数值设的很大，从输出的结果看，diag 是把lower里面的数加和了，也就是把列中的数加和了。
这个结果让我有点儿弄不明白。
-
@mengweilm425 这个问题cfd-online上也有人讨论过，不过没有结果。我个人觉得这确实是OpenFOAM的bug，不过由于fvm离散后的系数矩阵的特殊性质，列元素相加和行元素相加得到的结果是一样的，因此这个bug不影响计算。有机会细说。
-