请教一个问题，关于Foam::fvMatrix<Type>::A() 函数的

yuan_neu

template<class Type>
Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
{
    tmp<volScalarField> tAphi
    (
        new volScalarField
        (
            IOobject
            (
                "A("+psi_.name()+')',
                psi_.instance(),
                psi_.mesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            psi_.mesh(),
            dimensions_/psi_.dimensions()/dimVol,
            zeroGradientFvPatchScalarField::typeName
        )
    );

    tAphi().internalField() = D()/psi_.mesh().V();
    tAphi().correctBoundaryConditions();

    return tAphi;
}

为什么要除以V()？？？感觉跟李东岳的icoFoam解析对不上了。
tAphi().internalField() = D()/psi_.mesh().V();

yuan_neu

为了跟动网格的处理在代码形式上统一吗？还是为了与边界处理的一致？

youmengtian

试着回答一下。本身A()这个成员函数的目的就是返回代数矩阵的格心系数。而实际当中fvmatrix在构造时都已经包含了网格体积，即完成了体积分或者面积分运算。因此在取回格心系数（central coeffcient）时需要格外剔除网格体积项。否则在内网格部分就和D()这个成员函数重复了。

yuan_neu

@youmengtian :kiss:

yuan_neu

@youmengtian 你还记得在哪些文件或者代码里描述了这种处理的吗？如果可以的话，指点一下。谢谢

youmengtian

simpleFoam

程迪

@yuan_neu

在icoFoam中插入以下代码

fvVectorMatrix UEqn
(
    fvm::ddt(U)
    +fvm::div(phi,U)
    -fvm::laplacian(nu,U)
)
Info<<"fvc::laplacian(nu,U).dimensions="<<fvc::laplacian(nu,U)().dimensions()<<endl;
Info<<"fvm::laplacian(nu,U).dimensions="<<fvm::laplacian(nu,U)().dimensions()<<endl;
Info<<"phi.dimensions="<<phi.dimensions()<<endl;
Info<<"U.dimensions="<<U.dimensions()<<endl;
Info<<"fvc::div(phi,U).dimensions="<<fvc::div(phi,U)().dimensions()<<endl;
Info<<"fvm::div(phi,U).dimensions="<<fvm::div(phi,U)().dimensions()<<endl;
Info << "Dimension of UEqn is "<<UEqn.dimensions()<<endl;
// ...

输出大概是

fvc::laplacian(nu,U):[0 1 -2 0 0 0 0]
fvm::laplacian(nu,U):[0 4 -2 0 0 0 0]
phi:[0 3 -1 0 0 0 0]
U:[0 1 -1 0 0 0 0]
fvc::div(phi,U):[0 1 -2 0 0 0 0]
fvm::div(phi,U):[0 4 -2 0 0 0 0]
UEqn:[0 4 -2 0 0 0 0]
UEqn.A():[0 0 -1 0 0 0 0]
UEqn.H():[0 1 -2 0 0 0 0]

可以看出：

实际的情况是fvm项是乘以了体积的积分方程离散，输出类型是fvMatrix，fvc项是没有乘以体积的微分方程离散，输出类型是vol<Type>Field
从后面的代码UEqn == -fvc::grad(p)来看，二者是可以用operator==连接的。

根据fvMatrix的operator==源代码可以发现，operator==将GeometricField加到fvMatrix中的时候乘以了体积。

/*
/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
*/
//line: 1458
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
(
    const fvMatrix<Type>& A,
    const DimensionedField<Type, volMesh>& su
)
{
    checkMethod(A, su, "==");
    tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
    tC.ref().source() += su.mesh().V()*su.field(); //!!!
    return tC;
}

fvMatrix的部分函数比较诡异，UEqn.A()的量纲和UEqn相差[0 4 -1 0 0 0 0]，而UEqn.A()和H()相差待求变量速度的量纲[0 1 -1 0 0 0 0]，合理的解释是：
- UEqn是积分方程的离散，对于icoFoam，待求变量是速度U，量纲是[0 1 -1 0 0 0 0]；
- UEqn方程系数和phi的量纲一致，是[0 3 -1 0 0 0 0]，而方程右端源项是[0 4 -2 0 0 0 0]
- UEqn.A()是对角线系数除以体积，所以是[0 0 -1 0 0 0 0]
- UEqn.H()是非对角线部分作用于待求变量并被方程右端源项减去，再除以体积得到的，所以是[0 1 -2 0 0 0 0]
如果你详细看源代码，你会发现UEqn.D()是无量纲的，而UEqn.A()是强行把量纲附上去的，我觉得OF这些地方是有偷懒或者不一致的地方的。

yuan_neu

@程迪 :kiss_big: :kiss_big: :kiss_big: 牛

yuan_neu

@程迪如果你详细看源代码，你会发现UEqn.D()是无量纲的，而UEqn.A()是强行把量纲附上去的，我觉得OF这些地方是有偷懒或者不一致的地方的。

这句话我也有所感觉。