关于源项处理方法fvm::su和fvm::sp的疑问
-
一般形式的对流扩散方程为
将对流项和扩散项的体积分转为面积分：
源项可以写作
在《OpenFOAM中的有限体积离散》（原文链接：http://marinecfd.xyz/post/openfoam-finite-volume-discretization/）中，有如下表述：
fvm::su 将 su 的值乘以网格体积后施加到 source() 上：// OpenFOAM/src/finiteVolume/finiteVolume/fvm/fvmSup.C#L32 template<class Type> Foam::tmp<Foam::fvMatrix<Type>> Foam::fvm::Su ( const DimensionedField<Type, volMesh>& su, const GeometricField<Type, fvPatchField, volMesh>& vf ) { const fvMesh& mesh = vf.mesh(); tmp<fvMatrix<Type>> tfvm ( new fvMatrix<Type> ( vf, dimVol*su.dimensions() ) ); fvMatrix<Type>& fvm = tfvm.ref(); fvm.source() -= mesh.V()*su.field(); return tfvm; }fvm::sp 将 sp 的值乘以网格体积后施加到对角元素 diag() 上：
// OpenFOAM/src/finiteVolume/finiteVolume/fvm/fvmSup.C#L98 template<class Type>template<class Type> Foam::tmp<Foam::fvMatrix<Type>> Foam::fvm::Sp ( const volScalarField::Internal& sp, const GeometricField<Type, fvPatchField, volMesh>& vf ) { const fvMesh& mesh = vf.mesh(); tmp<fvMatrix<Type>> tfvm ( new fvMatrix<Type> ( vf, dimVol*sp.dimensions()*vf.dimensions() ) ); fvMatrix<Type>& fvm = tfvm.ref(); fvm.diag() += mesh.V()*sp.field(); return tfvm; }请问，在fvmSup.C文件中，
fvm.source() -= mesh.V()*su.field(); fvm.diag() += mesh.V()*sp.field();SuVp项不是本来就在等号右边吗？为什么代码中是减号？
SpVpϕp要移项至等号左边合并成系数矩阵，为什么代码中是加号？
谢谢！ -
@wwzhao 谢谢您的解答！
不过还有个小问题，根据您之前在其他地方的回答，在矩阵方程 [A][x]=[b] 中，[b] 可以分为 internal field 产生的源项 source()，及 boundary conditions 对源项的影响 boundaryCoeffs_。
即：在[A][x]=[b]中，fvm.source() 是等号右边的项。
在方程fvm::ddt(U) +fvm::div(phi, U) == fvc::Su(...) fvc::Sp(...) fvc::SuSp(...) （==是-）中，（1）fvc::Su也是和源项在一侧的，还是不能理解 fvm.source() -= mesh.V()*su.field();
（2）同时，既然Sp项前面有减号，为什么在fvm.diag()上可以直接加Sp项呢？ fvm.diag() += mesh.V()*sp.field();
问题比较基础，还需多多请教。谢谢！ -
fvm.source() 并不是[A][x]=[b]等号右边的项，而只是[b]中的一部分，另一部分是boundaryCoeffs_，在求解线性方程时才会组装成完整的[b]。
关于SuSp的问题，你看的应该是
fvm::Su和fvm::Sp的代码，而不是fvc::。假设方程右端的源项需要隐式离散，进入[A]中，那么移项到左边后将使
fvm.diag()减小，前面的==实现了负号的功能，所以fvm::Su需要增加fvm.diag()。同理，若方程右端的源项需要显示离散，进入[b]中，那么将使
fvm.source()增大，前面的==实现了负号的功能，所以fvm::Sp需要减小fvm::source()。总而言之，
fvm::Su、fvm::Sp和==一起，是为了实现OpenFOAM 神奇的方程式代码与原始控制方程的一致性才这么设计其内部实现的。 -
@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 是个深坑，用起来要非常小心才行。
