HbyA,phiHbyA,fvc::div(phiHbyA)计算错误问题

yhdthu

最近做pEqn.H中的演算，输出了HbyA，phiHbyA和fvc::div(phiHbyA)，通过手动计算进行验证

首先测试了一个简单的331的立方体模型，由HbyA计算phiHbyA，再进行求和并除以体积求div，计算结果和输出的结果一致，没问题。

换了个模型，吊诡的事情出现了，由HbyA计算phiHbyA没问题，但是由phiHbyA进行求和算div居然是错的？！我很奇怪难道连加法也能算错？

以下是我的一个结果，测试的cell编号是47794：

手动计算的div：-6.3169e-3
输出的div：1.90343e-6

且不说正负号，cell体积的量级就是10^-8，手动求和phiHbyA的量级是10^-10，check了一下代码，也没有其它干扰了，实在想不通内部是怎么算的

李东岳

以下是我的一个结果，测试的cell编号是47794

不太清楚你怎么计算的。你要测试什么？计算$\nabla\cdot\mathbf{U}=0$？

phiHbyA进行求和算div居然是错的？

在OpenFOAM里面，方程$\nabla\cdot\mathbf{U}=0$的植入就是fvc::div(phi)，对应：fvc::surfaceIntegrate(phi)，对应面加和。

我对这个也感兴趣，能否提供更详细信息。

yhdthu

@李东岳前辈好，是这样，我打算测试压力泊松方程源项对计算结果的影响

fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)

以上是压力方程，我想看看**源项fvc::div(phiHbyA)**的计算结果和我自己手动计算的是否一样

所以下一步我需要求phiHbyA，我看了fvc::surfaceIntegrate的源码，fvc::div(phiHbyA)的计算就是将一个cell中不同face的通量phiHbyA相加，之后除以体积mesh.V()

关于phiHbyA的计算，查看代码如下，

surfaceScalarField phiHbyA
(
    "phiHbyA",
    fvc::flux(HbyA)
  + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);

有两个部分组成，一个是通过volVectorField类型的HbyA中心插值到face上再与面矢量做点积，第二个修正项是通过上一时间步对本时间的修正，具体算法代码中有

目前采用的都是2D的算例，一个cell有4个neighbour，具体细节如下：

对于简单算例，几何模型如下，是一个3 * 3 * 1的一个立方体，选择4号单元进行验证，具体数据如图

由HbyA手动计算得得phiHbyA与输出的结果一致，再由phiHbyA计算fvc::div(phiHbyA)也没问题；

对于网格量比较大的复杂问题

由HbyA手动计算得得phiHbyA与输出的结果一致但是，由phiHbyA求和除以mesh.V()所得的结果，与输出的fvc::div(phiHbyA)相差了好几个量级！也就是之前列出来的结果了。

求前辈指教，我现在不懂为啥这个求和还能求错了？或者是求这个fvc::div(phiHbyA)时根本用的就不是phiHbyA求和来的？

李东岳

@yhdthu

#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"

    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

	surfaceScalarField phi
	(
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
       	fvc::interpolate(U) & mesh.Sf()
    );

	Info << "fvc::interpolate(U) " << fvc::interpolate(U) << nl;

	Info << "phi is " << phi << nl;

	Info << "mesh surface vector is " << mesh.Sf() << nl;

	Info << "fvc::div(phi) is " << fvc::div(phi) << nl;

    return 0;
}

请运行上述代码，并把U初始化为：

dimensions      [0 1 -1 0 0 0 0];

internalField   nonuniform List<vector> 
9
(
(0 1 0)
(0 2 0)
(0 3 0)
(0 4 0)
(0 5 0)
(0 6 0)
(0 7 0)
(0 8 0)
(0 9 0)
)
;


//internalField   uniform (0 1 0);

boundaryField
{
    wall
    {
        type            fixedValue;
        value           uniform (0 1 0);
    }
    defaultFaces
    {
        type            empty;
    }
}

我没看出计算错误的问题。fvc::div就是对速度做体积分，通过加和的方式计算。

李东岳

Yes, 我觉得这个算例非常好，因为OKSS2上来就是高斯积分计算，也算是我对你这个感兴趣的动力，我会把这个放在课上分析。

由phiHbyA求和除以mesh.V()所得的结果，与输出的fvc::div(phiHbyA)相差了好几个量级！

用上面那个简单试试，可能你那个太复杂了，小数点很多，数也不工整。客观来讲，不会出现这个现象。

yhdthu

@李东岳前辈您说的这个算例我没异议，这个就和我的那个简单例子一样的，我目前已经掌握HbyA到fvc::div(phiHbyA)的计算方法。

但是，问题来了，我用相同的计算方法去计算我的第二个例子的某个cell，发现求和的时候出现了错误。

即，我手动计算的fvc::div(phiHbyA)为：-6.3169e-3
但是输出的结果却是1.90343e-6

如果只看phiHbyA求和，其结果是10^-10量级，cell的体积是10^-7量级，两个相除就是10^-3次方量级了；

输出的fvc::div(phiHbyA)的结果居然是10^-6量级，这可是除完体积之后的值，也就是说之前的phiHbyA的求和量级是10^-13量级

前辈，这就是我想表达的意思，为什么从phiHbyA就算不到fvc::div(phiHbyA)了？？

李东岳

所以就是说：你实际上已经知道手动计算的方式，但是你算的和代码结果不一样。如果可以的话，请把4号网格4个面的phiHbyA输出，我看看你是怎样 手动计算的fvc::div(phiHbyA) 的，就是说：把你计算-6.3169e-3的过程写出来。

输出的fvc::div(phiHbyA)的结果居然是10^-6量级

这个非常正常，匹配你的压力方程残差。
你可以看出来，这只不过是加减法问题，可能你手动计算存在一个非常小的bug。

yhdthu

@李东岳前辈，我的意思是：

我计算那个方形的4号cell全都没问题，以上第一个方形网格的算例给出了我所查找的所有参数；

有问题的是，用相同的计算方法，去计算第二个算例的47794号cell，从phiHbyA到fvc::div(phiHbyA)的计算出了问题。

我明白量级什么的无所谓，最重要的是数得算的对，但是现在手动算的和程序算的根本就对不上了，而且差了n个量级，这就说不过去了:crying:

yhdthu

@李东岳前辈具体的计算过程如下，这里只说有错误的第二个算例（cell_47794）

需要进行计算的场输出为：
HbyAVol = HbyA;
phiHbyASuf = phiHbyA;
sourcePoissionU = fvc::div(phiHbyA);

根据输出的phiHbyA，4个面的值如下：

E：0.000125002
W：0.000124669
N：-1.94136e-05
S：-1.908e-05
F 和 B 为empty，0

加和：E - W + N - S = 6e-10

mesh.V() = 9.49835e-8

fvc:div(phiHbyA) = (E - W + N - S) / mesh.V() = -0.0063169

但是，输出的fvc::div(phiHbyA) = 1.90343e-6

李东岳

E：0.000125002
W：0.000124669
N：-1.94136e-05
S：-1.908e-05

这个找的正确么

yhdthu

@李东岳前辈，这个百分百没问题，我是用paraview一个一个面点出来看的

李东岳

你是怎么在paraview查看face的值的？

yhdthu

@李东岳通过paraview找到cell编号，查找polymesh里面的owner和neighbor文件定位行号，打开surfacefield对应找就好啦～

李东岳

方法是对的，不过为什么错了？我也很奇怪。

yhdthu

@李东岳更奇怪的是，程序前半段计算phiHbyA都是对的，但是之后求div的时候（加法）却错了，难道加法也能算错:upset:

yhdthu

@李东岳前辈，解决了，原来是输出精度不够导致的，我的锅:lol: