关于 phi 和 fvVectorMatrix 的两个问题

wwzhao

lduMatrix::diag()和fvMatrix<Type>::D()并不等价。相比lduMatrix::diag()，fvMatrix<Type>::D()考虑了边界对系数矩阵的影响，具体代码对应这句：addCmptAvBoundaryDiag(tdiag.ref());。

CFD中文网

有关diag()和D()严格来说是的。你可以测试一下然后给楼主提供更详细的信息。:sunglasses:

wwzhao

@cfd-china 离散的具体实现可参阅相关源码。

以fvm::ddt(T)为例，调用的是

template<class Type>
tmp<fvMatrix<Type> >
ddt
(
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return fv::ddtScheme<Type>::New
    (
        vf.mesh(),
        vf.mesh().ddtScheme("ddt(" + vf.name() + ')')
    )().fvmDdt(vf);
}

这里调用了fv::ddtScheme<Type>::New，这个函数是一个runtime construction的实现，即在运行时读取字典文件中的离散格式并构造相应的对象：

template<class Type>
tmp<ddtScheme<Type> > ddtScheme<Type>::New
(
    const fvMesh& mesh,
    Istream& schemeData
)
{
    if (fv::debug)
    {
        Info<< "ddtScheme<Type>::New(const fvMesh&, Istream&) : "
               "constructing ddtScheme<Type>"
            << endl;
    }

    if (schemeData.eof())
    {
        FatalIOErrorIn
        (
            "ddtScheme<Type>::New(const fvMesh&, Istream&)",
            schemeData
        )   << "Ddt scheme not specified" << endl << endl
            << "Valid ddt schemes are :" << endl
            << IstreamConstructorTablePtr_->sortedToc()
            << exit(FatalIOError);
    }

    const word schemeName(schemeData);

    typename IstreamConstructorTable::iterator cstrIter =
        IstreamConstructorTablePtr_->find(schemeName);

    if (cstrIter == IstreamConstructorTablePtr_->end())
    {
        FatalIOErrorIn
        (
            "ddtScheme<Type>::New(const fvMesh&, Istream&)",
            schemeData
        )   << "Unknown ddt scheme " << schemeName << nl << nl
            << "Valid ddt schemes are :" << endl
            << IstreamConstructorTablePtr_->sortedToc()
            << exit(FatalIOError);
    }

    return cstrIter()(mesh, schemeData);
}

以常见的Euler格式为例，将会构造EulerDdtScheme类型的对象，并用其fvmDdt()方法返回的fvMatrix，看fvmDdt()的实现：

template<class Type>
tmp<fvMatrix<Type> >
EulerDdtScheme<Type>::fvmDdt
(
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<fvMatrix<Type> > tfvm
    (
        new fvMatrix<Type>
        (
            vf,
            vf.dimensions()*dimVol/dimTime
        )
    );

    fvMatrix<Type>& fvm = tfvm();

    scalar rDeltaT = 1.0/mesh().time().deltaTValue();

    fvm.diag() = rDeltaT*mesh().Vsc();

    if (mesh().moving())
    {
        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc0();
    }
    else
    {
        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc();
    }

    return tfvm;
}

可以看到fvm::ddt(T)产生的fvm.diag()为时间倒数乘以网格单元体积，其产生的源项fvm.source()为时间倒数乘以上一时刻的T乘以网格单元体积。

对于fvm::div及fvm::laplacian可以做类似分析。

插一句，fvm::div离散不会产生source()，fvm::laplacian会产生source()。

mengweilm425

@wwzhao @李东岳 @cfd-china

非常感谢大家的回复！我现在的阶段还看不好代码，但是我又对fvMatrix做了一些测试。

还是如下的fvMatrix

    fvVectorMatrix UEqn_diffusion
    (
       fvm::laplacian(nu, U)
    );

我发现在解这个方程的时候

\begin{equation}
[A][x]=[b]
\end{equation}

矩阵[A] 实际上是 UEqn.D() + UEqn.upper() + UEqn.lower(), 而没有使用UEqn.diag()。UEqn.D()在UEqn.diag()的基础上加上了边界条件对当前网格的影响。上面方程中的[b]，也不是UEqn.source()，而是和边界条件中已知的速度相关的值。我直接求解得到的结果和openfoam算出来的结果是一样的。由于没有真正的去看代码，不知道Openfoam里面有没有相关的function可以直接输出上述$[A][x]=[b]$方程组中的[A]和[b]，还请各位大神指点。

再次感谢大家的帮助！

wwzhao

@mengweilm425 $[A][x]=[b]$ 中的 [A] 可以看成是 internal field 离散后的系数矩阵和 boundary conditions 对系数矩阵的影响两部分之和，其中 internal field 离散之后的系数矩阵表示为 upper()，diag() 和 lower()。boundary condition 对系数矩阵的影响则放在 internalCoeffs_ 里面。[b] 也一样，可以分为 internal field 产生的源项 source()，及 boundary conditions 对源项的影响 boundaryCoeffs_。

OpenFOAM中没有直接输出 [A] 和 [b] 的函数，关于如何将 boundary conditions 的影响加到 fvMatrix 里面，你可以参考 addBoundaryDiag，addCmptAvBoundaryDiag 及 addBoundarySource 这几个函数。

@wwzhao 在关于 phi 和 fvVectorMatrix 的两个问题中说：

@李东岳 fvMatrix的source并不是线性代数方程组Ax=b中的b，而是指加在网格单元中心的力源项。

上面的帖子描述不太准确，确切说 [b] 是 internal field 产生的源项 source() + boundary field 产生的源项 boundaryCoeffs_。

mengweilm425

@wwzhao 谢谢您的回复，这就明朗多了。我再好好看看这部分代码。:happy:

vivian

@mengweilm425 您好看了你的研究过程，觉得很好。我也想这样分步骤运行各个程序，以分清各个量含义。请问您是怎么做到的。我先编译了自己求解器，在文件中加入自己想输出的量，希望在输出中出现想看到的量。结果并没有看到想看的内容。请问您是如何。。

CFD中文网

@vivian
我猜测你的求解器名字没有替换？如果编译后的是icoFoam2，记得输入icoFoam2而不是原来的icoFoam..

vivian

@cfd-china 您好，谢谢您的回复。是替换了的，是我给他速度赋值赋错啦。。

mengweilm425

@vivian @cfd-china 看来我来晚了:happy:

Wayne

厉害了我的哥哥姐姐们

CFD中文网

关于 phi 和 fvVectorMatrix 的两个问题