关于源项处理方法fvm::su和fvm::sp的疑问
-
@wwzhao 我也有同样的疑惑，不知道为什么
Su这个函数里边是负号。
如果我把这一项放在方程左边呢？它的作用是使整个fvMatrix的source减去一个数啊。fvm::ddt(U) + fvm::div(phi, U) + fvm::Su(sourceField, U)我之前是这样的理解的：
[A]x = [b]对应是[A]x - [b]。
OF 中的fvMatrix的diag对应A，source对应-b。
这样就能解释为什么对一个矩阵+::Su(sourceField, x)，后果却是这个矩阵的source被减去了一个数mesh.V()*sourceField.field()。
因为source减去一个数，相当于[A]x - [b]里边的 b 增加了这个数：[A]x - [b + mesh.V()*sourceField.field()]。
我这个理解很大概率是错的。所以求大神指正！Foam::tmp<Foam::fvMatrix<Type>> Foam::fvm::Su ( const DimensionedField<Type, volMesh>& su, const GeometricField<Type, fvPatchField, volMesh>& vf ) { fvm.source() -= mesh.V()*su.field(); } -
@史浩 在 关于源项处理方法fvm::su和fvm::sp的疑问 中说：
我是这样理解的，在fvMatrix中，对角线元素本省就是带了一个负号。在组装矩阵的时候，先确定好非对角线元素的值，然后对角线元素的值为该行非对角线元素和的相反数，也就是对角线元素一开始就带了一个负号。关于符号问题，对于每个项的符号也略有差别，所以会出现一些困惑。
这一部分你可以参考fvm::laplacian中关于拉普拉斯项的隐式离散和fvm::ddt的离散来了解一下矩阵组装中的符号问题
加油~对角元素是正数，而不是负数。非对角元素可正可负。
对角元素占优才能迭代求解，否则就成了病态矩阵。
这点可以通过输出 fvm.diag() 确认，参考这篇博文。
-
@wwzhao
谢谢您更新了我的认知！
现在我消化一下：
假如我们要求解方程 $\frac{\ \rho k}{\p t} = -30 + 50$。
左边这项当然变成了：fvm::ddt(rho, k)
右边第一项可以变成：- fvm::Sp(30/k, k)，移动这一项到左边，即把负号去掉了，然后把Sp的定义代入，即变成了fvm.diag += V*30/k。这项加到对角元素上去了，和我们要求解的方程是一致的，并且使对角更占优了。
右边第二项可以变成：fvm::Su(50, k)，移动这一项到左边，并且把Su的定义代入，即变成了fvm.source += V*50，与我们要求解的方程是一致的。如果右边第一项写成这样：
- fvm::Su(30, k)，那么它移动到左边后，将会使fvm.source -= V*30，与我们要求解的方程是一致的。
如果右边第二项写成这样：fvm::Sp(50/k, k)，那么它移到左边后，将会使fvm.diag -= V*50/k。这不利于对角占优，因此不推荐这样做。
结论：对于方程右边的项：
负数，写成Sp——推荐这样做，因为可以使对角更占优。
负数，写成Su——可以这样做。
正数，写成Sp——不推荐这样做，因为会损害对角占优。
正数，写成Su——应该这样做。
还有一种没讨论的情况：正负未知就写成SuSp。
@wwzhao 请问我这里的理解正确吗？ -
总结的不错。
Su Sp 是灵活处理源项的方法，正确使用方法如下：
1、若源项为负，如你举的例子中的
-30，一般写作 $\frac{-30}{k^{old}}k$，$k^{old}$ 代表上一时刻的值，采用fvm::Sp离散后，源项就进入 [A] 中对角元素，使对角元素更占优。2、若源项为正，如你举的例子中的
50，则直接用fvm::Su离散，源项进入右端项 [b] 中。3、若源项是一个有正有负的 volScalarField，则利用
fvm::SuSp离散，它会判断源项每个网格单元的正负号，然后选择用 Su 还是 Sp 的方法进行处理。以上方法的最终目的都是增加矩阵对角占优，使线性方程组的迭代更稳定，更容易收敛。
-
@wwzhao 感谢您的答疑解惑！现在已经清楚
Su和Sp的作用了。但是对于SuSp仍有疑问：
它的代码简化如下：Foam::fvm::SuSp ( const DimensionedField<scalar, volMesh>& susp, const GeometricField<Type, fvPatchField, volMesh>& vf ) { fvm.diag() += mesh.V()*max(susp.field(), scalar(0));//susp为正，这句话起作用 fvm.source() -= mesh.V()*min(susp.field(), scalar(0))*vf.internalField();//susp为负，这句话起作用 }仍然考虑方程右边有两个源项，假设是 $\frac{\partial \rho k}{\partial t} = A+B$。
先假设 A B 正负未知，代码可以写成：== SuSp(A/k, k) + SuSp(B/k, k)。
实际计算时，假设 A=-30，B=50。
那么根据SuSp的定义，代码变成了：fvm.source() -= V*30 fvm.diag() += V*50/k移到方程左边，则变成了：
fvm.source() += V*30 fvm.diag() -= V*50/k这里却削弱了对角占优啊。
不知道哪里出现了问题。
该问题在 CFDonline 有人提出过，但是我到现在为止仍看不明白。 -
你这个问题很有意思，从你的推导分析来看，SuSp 确实是削弱了对角占优。
但其中有一处错误：V30 应该为 V-30。不过这改变的是源项，而与对角元素无关。我们一起分析下是怎么回事。
如果将方程按照这种形式写：$\frac{\partial k}{\partial t}=-(-A)-(-B)$
那么代码可以写成：
== -fvm::SuSp(-A/k, k) - fvm::SuSp(-B/k, k)。实际计算时，假设 A=-30，B=50。
那么 $-A/k>0$，$-B/k<0$，fvm前面的负号与==抵消，代码实际执行的是fvm.diag() += V*30 fvm.source() -= V*(-50)*k这又的确增加了对角占优，显然
fvm::SuSp不能随心所欲地用！-fvm::SuSp(-A/k, k)和fvm::SuSp(A/k, k)的结果不一样，而前一种写法才是我们想要的结果。我翻了 OpenFOAM 的代码后发现，凡是在
==右边使用的 Sp 和 SuSp 前面必有负号，凡是在==左边使用的 Sp 前面必为正号 (代码中基本没有出现fvm::Su)。看来 SuSp 是个深坑，用起来要非常小心才行。