关于源项处理方法fvm::su和fvm::sp的疑问

wh3296

一般形式的对流扩散方程为

将对流项和扩散项的体积分转为面积分：

源项可以写作

在《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

OpenFOAM中的方程一般写成这种形式

   fvm::ddt(U)
 + fvm::div(phi, U)
... 
==
   fvc::Su(...)
   fvc::Sp(...)
   fvc::SuSp(...)

fvMatrix 的 == 操作符执行优先级最低。先将 == 左右的各项相加，最后执行 ==。

而 == 操作符被重载了，实际上是 - ，所以当对源项进行操作的时候是需要减去 SuVp，同理，对系数矩阵的修改需要用加号。

wh3296

@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();
问题比较基础，还需多多请教。谢谢！

史浩

我是这样理解的，在fvMatrix中，对角线元素本省就是带了一个负号。在组装矩阵的时候，先确定好非对角线元素的值，然后对角线元素的值为该行非对角线元素和的相反数，也就是对角线元素一开始就带了一个负号。关于符号问题，对于每个项的符号也略有差别，所以会出现一些困惑。
这一部分你可以参考fvm::laplacian中关于拉普拉斯项的隐式离散和fvm::ddt的离散来了解一下矩阵组装中的符号问题
加油~

wwzhao

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 神奇的方程式代码与原始控制方程的一致性才这么设计其内部实现的。

李东岳

多谢各位高手分享讨论！

wh3296

@wwzhao
举个例子的话，我的理解是：
对于方程 fvm::ddt(U)+fvm::div(phi, U)== fvc::Su(...) ，fvc::Su(...)在等号的源项一侧，fvc::Su(...)项的存在使矩阵方程 [A][x]=[b] 的 [b] 项增加；
而 fvm.source() -= mesh.V()*su.field();这一行的意思是在源项中减去SuVp；
这个地方还是理解不了...请指教。

wwzhao

@wh3296 == 相当于减号，负负得正。。。

浪迹天大

@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();
}

wh3296

@史浩我好好看一下，谢谢！

wwzhao

@浪迹天大

因为源项通常写在控制方程等号的右边，为了保证源代码与控制方程形式统一， fvm::Su fvm::Sp 以及 fvm::SuSp 必须写在 == 右边。

你可以看一下 OpenFOAM 中湍流模型中的用法。

wwzhao

@史浩在关于源项处理方法fvm::su和fvm::sp的疑问中说：

我是这样理解的，在fvMatrix中，对角线元素本省就是带了一个负号。在组装矩阵的时候，先确定好非对角线元素的值，然后对角线元素的值为该行非对角线元素和的相反数，也就是对角线元素一开始就带了一个负号。关于符号问题，对于每个项的符号也略有差别，所以会出现一些困惑。
这一部分你可以参考fvm::laplacian中关于拉普拉斯项的隐式离散和fvm::ddt的离散来了解一下矩阵组装中的符号问题
加油~

对角元素是正数，而不是负数。非对角元素可正可负。

对角元素占优才能迭代求解，否则就成了病态矩阵。

这点可以通过输出 fvm.diag() 确认，参考这篇博文。

浪迹天大

@wwzhao
谢谢您更新了我的认知！
现在我消化一下：
假如我们要求解方程 $\frac{\partial \rho k}{\partia 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 请问我这里的理解正确吗？

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 的方法进行处理。

以上方法的最终目的都是增加矩阵对角占优，使线性方程组的迭代更稳定，更容易收敛。

wh3296

@wwzhao 明白了，谢谢老师！

浪迹天大

@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 有人提出过，但是我到现在为止仍看不明白。

wwzhao

@浪迹天大

你这个问题很有意思，从你的推导分析来看，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 是个深坑，用起来要非常小心才行。