首发于有所思
【翻译】OpenFOAM中的空间离散化与矩阵系数

【翻译】OpenFOAM中的空间离散化与矩阵系数

翻译者的话:我在研究icoFoam的过程中,遇到了一些问题,主要是对fvMatrix的理解不够深入。之后在网络上找到了一份很好的文献(OpenFOAM -空間の離散化と係数行列の取り扱い),作者是Fumiya Nozaki,他画了很多漂亮的示意图,讲述了一维导热问题的fvMatrix例子,现在整理一下重点,分享给大家,里面掺杂了一些我自己用代码还原矩阵的过程,帮助大家理解OpenFOAM里面矩阵运算的过程。

OpenFOAM中的数学模型,最终是要解线性方程组Ax=b

矩阵A中非0元素所代表的意义,是相邻两个cell之间的关系。

矩阵A的数据以LDU的方式分三块保存在iduMatrix类里。

但是LDU只保存非0元素,因为矩阵大部分都是0,把那么多0保存下载没意义。

解决的方法是LDU里只保存非0元素,然后用一个IduAddressing的类来保存矩阵的哪些位置有值这样的信息,如此就可以把LDU里保存的非0元素一个个填回矩阵A里,还原这个矩阵的本来面目。

举个例子。5号cell和2号cell之间有接触,所以A[5][2]就不为0,同理A[2][5]也不为0。

fvMatrix里面不止有IduMatrix和IduAddressing,还有其他一些变量。我最关心的b就是source,其他边界条件对矩阵A和b的影响主要是通过internalCoeffs和boundaryCoeffs实现的,二者的区别是internalCoeffs修改A,boundaryCoeffs修改b。

昏昏欲睡的同学们,重头戏到了:

设定一个一维导热问题,分析一下OpenFOAM是怎么运算的。

这个问题是有解析解的。

修改系统自带的laplacianFoam,改成你自己定义的新Foam,把其中的公式改成这个问题里的方程:

fvScalarMatrix TEqn(fvm::laplacian(k, T));
TEqn.solve();

跑一跑,结果与解析解一致。

现在修改代码,令其输出矩阵A和b:

// Assigning diagonal coefficients
for(label i=0; i<Nc; i++)
{
    A[i][i] = TEqn.diag()[i];
}
// Assigning off-diagonal coefficients
for(label faceI=0; faceI<TEqn.lduAddr().lowerAddr().size(); faceI++)
{
    label l = TEqn.lduAddr().lowerAddr()[faceI];
    label u = TEqn.lduAddr().upperAddr()[faceI];
    A[l][u] = TEqn.upper()[faceI];
    A[u][l] = TEqn.upper()[faceI];
}
// Assigning contribution from BC
forAll(T.boundaryField(), patchI)
{
    const fvPatch &pp = T.boundaryField()[patchI].patch();
    forAll(pp, faceI)
    {
        label cellI = pp.faceCells()[faceI];
        A[cellI][cellI] += TEqn.internalCoeffs()[patchI][faceI];
        A.source()[cellI] += TEqn.boundaryCoeffs()[patchI][faceI];
    }
}

可以看出主要给A和b赋值的就是TEqn这个类,这个类下面又有diag(),lduAddr(),upper(),internalCoeffs()和boundaryCoeffs()在起作用。

看起来lower()和source()项都没用上?

想了一下,可能是因为这个case里upper和lower是对称的,所以可以省略lower。真实计算时应该是:

A[l][u] = TEqn.upper()[faceI];
    A[u][l] = TEqn.lower()[faceI];

这样也比较符合上面那张矩阵图里的定义:

source项没有用上应该是因为这个公式里没有source,所以可以省略,真实计算时应该有:

A.source()=TEqn.source()

现在我们推导一下,如何从TEqn类下面的diag(),lduAddr(),lower(), upper(),source(), internalCoeffs()和boundaryCoeffs()算出矩阵A和向量b。

令solver在解方程前后各输出一次变量:

fvScalarMatrix TEqn(fvm::laplacian(k, T));
Info<< "TEqn = " << TEqn << endl;
Info<< "TEqn.lduAddr().lowerAddr() = " << TEqn.lduAddr().lowerAddr() << endl;
Info<< "TEqn.lduAddr().upperAddr() = " << TEqn.lduAddr().upperAddr() << endl;
TEqn.solve();
Info<< "T = " << T << endl;

结果:

//输出lower,diag和upper
TEqn = false true true 5(-100 -200 -200 -200 -100)4(100 100 100 100)
//输出dimension,不重要
[1 2 -3 0 0 0 0]
//输出source,也就是b
5{0}
//输出internalCoeffs
4
(
1(-200)
1(-200)
0()
0()
)

//输出boundaryCoeffs
4
(
1(-20000)
1(-100000)
0()
0()
)

//输出lduAddr的位置信息
TEqn.lduAddr().lowerAddr() = 4(0 1 2 3)
TEqn.lduAddr().upperAddr() = 4(1 2 3 4)
DICPCG:  Solving for T, Initial residual = 1, Final residual = 1.21266e-16, No Iterations 1
//输出更新后的T
T = dimensions      [0 0 0 1 0 0 0];

internalField   nonuniform List<scalar> 5(140 220 300 380 460);

boundaryField
{
    left
    {
        type            fixedValue;
        value           uniform 100;
    }
    right
    {
        type            fixedValue;
        value           uniform 500;
    }
    topAndBottom
    {
        type            empty;
    }
    frontAndBack
    {
        type            empty;
    }
}

我们操作一下:

Fig.1 矩阵操作流程

最后得到一个矩阵A和向量b,带入Matlab里算一算结果:

与代码中的计算结果对比一下:

¥n==Coefficient Matrix==
-300 100 0 0 0 -20000
100 -200 100 0 0 0
0 100 -200 100 0 0
0 0 100 -200 100 0
0 0 0 100 -300 -100000
Solution: 5(140 220 300 380 460)

一模一样,说明OpenFOAM里也是这样操作的。


这样你就满足了吗?不,你还不够满足。因为这里的输出矩阵函数是我自己写的,不是OpenFOAM自带的函数,OpenFOAM自己没有一个函数给你输出这个矩阵的完整样子。所以我写的代码存在着与OpenFOAM不一致的可能。你需要一种方法,确认OpenFOAM中最终就是得到了这样的矩阵A和b。


这个方法就是gdbOF

简单说,gdb是一个linux下调代码用的通用工具,gdbOF是OpenFOAM社区开发的专门针对OpenFOAM的简单版gdb,他主要有两个功能:

  • 显示离散化矩阵的元素值;
  • 显示cell和boundary上的值。

我们现在要利用的就是第一个功能,其实也就第一个功能有用。cell和boundary上的值用别的方法也可以看到。

下载地址在这儿:

Contrib gdbOF - OpenFOAMWikiopenfoamwiki.net

按照网页里附赠的manual安装,这里有个细节,把openfoam调成debug模式后要重新编译一次,manual里忘记说了。这也意味着你必须使用从源码编译的OpenFOAM。使用傻瓜安装的朋友们哭了,但是这也没办法,这就是人生,这就是阶级,无产阶级哭,资产阶级笑。试想,假如你出生自带2亿家产,你还会玩OpenFOAM这种穷人游戏吗?这样想想你就心理平衡了,人呐,要认命。你已经够幸运了,还能看到这篇攻略,写攻略的人浪费的时间比你还多。比起我,你节约了一天的时间,比起这篇文章的作者Fumiya Nozaki,你节约了一个月的时间。

cd $WM_PROJECT_DIR
wclean all
./Allwmake

装好gdbOF之后,记得把刚刚编译好的新solver再编译一次,因为地址改变了,gdb找不到新的solver可执行文件的位置:

cd mylaplacianFoam
wclean && wmake

装好之后:

这一步之后假如出现“No symbol table is loaded. Use the "file" command.”的提示:

exec-file mylaplacianFoam

按照上述操作,可以得到矩阵。假如你不想输出那么多数,只想输出特定位置的矩阵元素,就按照作业6中的命令输:

pfvmatrixfull this matrix.txt 1 1 3 3

最后把等式右边的b也输出了,对照一下,确实A和b都跟我们算出来的一样,这下就放心了,OpenFOAM里确实就是这么运算的。

这样你就满足了吗?不,你还不够满足。因为即使你看到了Matrix,你还是不能100%确定他就是按照你脑海中的equation和boundary condition和scheme算出来的。你需要分析Ueqn中的每一项,看看他们是否是按照你脑海中的设想给Matrix赋值的。


这个case的公式非常简单,就是:

放在OpenFOAM里就是一个项:

TEqn(fvm::laplacian(k,T));

所以我们只需要计算一个项,就可以得到方程组了。这应该是OpenFOAM里最简单的case。我们现在的目的就是用gdbOF反推,看看每个点上是怎么算的。

关于这个fvm::laplacian(k,T)的执行情况,需要持续追踪,具体调用的函数过程如下:

陈与论:用GdbOF追踪laplacian函数在OpenFOAM中的定义过程zhuanlan.zhihu.com图标

如同文中所说,主要干活的只有两个文件:

  • gaussLaplacianSchemes.C (算b)
  • gaussLaplacianScheme.C (算A)

gaussLaplacianSchemes.C里把干活的函数定义成macro(宏)了,这非常讨厌,因为写成宏的function就不能用gdb一步一步测试了。

我对此也没办法,只能绕开宏的部分,把不是宏的部分检测一下。可以看出gaussLaplacianSchemes.C是先调用了gaussLaplacianScheme.C里的fvmLaplacianUncorrected函数来形成一个fvMatrix,然后再给这个fvMatrix的source项做些修改。那么我们直接定位fvmLaplacianUncorrected函数:

gdb mylaplacianFoam
b gaussLaplacianScheme.C:52
r

依次输出这三行命令,就会进入第一次调用fvmLaplacianUncorrected函数的场景:

函数中赋值的语句就是如下两句:

经过测试,这两句之后的矩阵fvm就是“Fig.1 矩阵操作流程”图中的第4个矩阵。

那么这两个语句是如何运作的呢?

第一句:

deltaCoeffs在宏里面,我看不到。我估计就是1/d,d是两个cell中心的间距。

https://www.cfd-online.com/Forums/openfoam/109640-some-questions-about-laplacian-nonorthogonal-correction.htmlwww.cfd-online.com

上网求证了一下,看到Fumiya Nozaki(也就是这个ppt的作者)的回答:

Therefore, the coefficients calculated in the following line are gamma_f * |Delta|/|d|, where Delta is the vector defined in the Dr. Jasak's thesis(Eq. (3.32)).
fvm.upper() = deltaCoeffs.internalField()*gammaMagSf.internalField();

看来我们的意见是一致的:

upper=k*\frac{|\Delta|}{d}

其中gamma就是热导率k, \Delta 是带方向的面积向量,d是两个cell中心的间距。

这里k=1000,delta=0.01,d=0.1,所以upper的值都是100.

第二句:

这里的意思也很明确,就是每个cell都减去跟upper一样的值。

diag=-k*\frac{|\Delta|}{d}

拿笔算一下,没毛病。唯一比较疑惑的是为什么lower也有值。明明前面只对upper赋值了。

找了一下,OpenFOAM guide/Matrices in OpenFOAM

fvm.lower is never mentioned - this is because it is equal to fvm.upper and a symmetric matrix is being created.

查了一下代码,确实是这样,这个对称的操作是在呼叫lower()或upper()时完成的,upper和lower都是:

  • 假如两个都没被赋值,被呼叫时就会赋值0;
  • 假如其中一个被呼叫时自己没被赋值而另一个有值,则直接把另一个link过来。
  • 假如其中一个被呼叫时自己有值,则不做动作。
Foam::scalarField& Foam::lduMatrix::lower()
{
    if (!lowerPtr_)
    {
        if (upperPtr_)
        {
            lowerPtr_ = new scalarField(*upperPtr_);
        }
        else
        {
            lowerPtr_ = new scalarField(lduAddr().lowerAddr().size(), 0.0);
        }
    }

    return *lowerPtr_;
}

感觉微妙,似乎设定死了这个scheme就是要对称的,不对称不行。

接着是设定边界cell的系数。也就是internalfield和boundaryfield的赋值。

InternalCoeffs=-k*\frac{|\Delta|}{d/2}

boundaryCoeffs=-k*\frac{|\Delta|}{d/2}*T_{right}

没毛病,这里的d是d/2。是因为最靠近边界的cell离墙壁只有1/2个cell。

接下来就是宏了。宏里面定义了source,这个明眼人都看得出来。问题是internalCoeffs和boundaryCoeffs没看到执行。这两个参数应该被用来修正矩阵A。我猜测这又是OpenFOAM的“打包”设计,把这个执行函数写进某个动作里去了。大致搜了一下github,感觉像是TEqn.solve()干的。

用GDB查了一下,果然如此:

基本路线是:

solve()->fvMatrix.solve()->solveSegregated()->addBoundaryDiag->addToInternalField

也就是说,在fvMatrix中储存的矩阵A和b,都是没有考虑边界影响的,直到solve解方程的时候再把边界效应考虑进去,这样做有什么好处,我暂时还不清楚。我猜测,可能是为了方便矩阵相加。加好了之后算的时候一起考虑边界层,这样能方便wall function之类的东西(都是我猜的,不一定对)。

关于用gdb调试宏,网上也没什么好方法,老司机告诉我,最好的办法就是对着源代码,一点一点的s和fin过去,不要n。听起来有点蠢,但确实是目前能采用的唯一方法了。


///////////////////////////////////////////////////////////////////////////


言归正传,我们得到了矩阵的两个公式:

upper=k*\frac{|\Delta|}{d}

diag=-k*\frac{|\Delta|}{d}

这两个公式看起来不太像laplacian项,这中间需要一些数学推导:

现在再看这两个公式,是不是顺眼多了?

upper=k*\frac{|\Delta|}{d}

diag=-k*\frac{|\Delta|}{d}

(顺便一提,linkin上显示,做这个ppt的Fumiya老哥从2016年开始已经2年处于“无职”状态了,不知道他现在过得怎么样。经常在twitter上发一些讲座信息,可能在专门搞培训吧。)

还有朋友会问:那你在laplacianSchemes里定义的那个corrected是干什么的呢?

实际计算中,可能会遇到复杂的mesh,比如这样的,当Sf面积法向量跟cell中心连线不平行的时候,就不能那么简单的离散了。

解决方法是将面积向量Sf做分解,一部分跟cell中心连线平行的照常带入矩阵A计算。另一部分跟接触面平行的带入右边的b计算。

Jasak的论文里介绍了3种非正交修正的方法,OpenFOAM里用的是第三种。区别就是delta的定义不同。

-----------------------------------------------------------------

(Jasak,OpenFOAM的主要作者之一(35%),克罗地亚人,1993年在伦敦帝国学院(也叫做帝国理工)读PhD期间开始参与写FOAM。后来回克罗地亚当了教授,又由于理念不合跟另一个主要作者Henry(60%)闹翻了,自己搞了一个Wikki公司。他的博士论文是OpenFOAM的入门经典教程之一。)

http://powerlab.fsb.hr/ped/kturbo/OpenFOAM/docs/HrvojeJasakPhD.pdfpowerlab.fsb.hr

以下链接是一些有趣的创业内幕,我知道你们都喜欢这种东西:

https://www.cfd-online.com/Forums/openfoam/152605-original-openfoam-author-s.htmlwww.cfd-online.com

-----------------------------------------------------------------

现在回答之前的问题:

laplacianSchemes里定义的那个corrected就是用来控制算不算非正交修正的。corrected就表示要算,uncorrected就表示不要算。

还有一种中庸的,limited。

算出来的结果当然是FVM,因为高斯定义已经对体积积分了,配合非常默契。可惜这种优点只能在laplacian运算中体现出来。其他运算一般还是直接乘以体积。

对于每个cell来说,laplacian实际上在分别计算相邻cell对他的影响,然后加到一起,就形成了现在的矩阵。将上面的矩阵拆开,就是5个方程。最上方和最下方的方程实际上还包涵了边界条件,用来表示边界条件对这个cell的影响,非常有趣。

a=k*\frac{|\Delta|}{d}


分析到此为止,提前祝大家春节快乐!

编辑于 2018-01-23

文章被以下专栏收录