Code of the Week #20: Sp(A,B)

OpenFOAM
1

#201907

今天，我们来聊一聊源项那些事儿。

源项在方程中是很难避免的。比如，在动量方程中，

\frac{\partial U}{\partial t} + \nabla\cdot(UU) = -\nabla p + \nabla\cdot\tau + S

S 就是源项，源项代表那些不能归类于时间项、对流项、扩散项的那些，如果源项可以线性化为

S=A\times U

则在写方程时，把源项加入

fvVectorMatrix UEqn
(
    fvm::ddt(U)
  + fvm::div(phi, U)
  - fvm::laplacian(nu,U)
  - fvm::Sp(A,U)
);

使用 fvm:: 代表隐式源项， U 是待求项，系数A 将会对待求矩阵的系数的主元产生影响，如果使用 fvc::Sp(A,U) 就会使用已知的U，然后这一项就会放在形成的矩阵方程等式的右边。

A 的类型应是场或一个标量，而 U 应是场。

fvc::Sp() 位于 src/finiteVolume/finiteVolume/fvc/fvcSup.C，

template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> >
Sp
(
    const volScalarField& sp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return sp*vf;  // 直接返回两者相乘结果
}

template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> >
Sp
(
    const tmp<volScalarField>& tsp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return tsp*vf;
}


template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> >
Sp
(
    const dimensionedScalar& sp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return sp*vf;
}

隐式离散的 fvm::Sp() 位于src/finiteVolume/finiteVolume/fvm/fvmSup.C

template<class Type>
Foam::tmp<Foam::fvMatrix<Type> >
Foam::fvm::Sp
(
    const DimensionedField<scalar, volMesh>& sp,
    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();

    fvm.diag() += mesh.V()*sp.field(); // 返回系数矩阵对角线

    return tfvm;
}

template<class Type>
Foam::tmp<Foam::fvMatrix<Type> >
Foam::fvm::Sp
(
    const tmp<DimensionedField<scalar, volMesh> >& tsp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<fvMatrix<Type> > tfvm = fvm::Sp(tsp(), vf);
    tsp.clear();
    return tfvm;
}

template<class Type>
Foam::tmp<Foam::fvMatrix<Type> >
Foam::fvm::Sp
(
    const tmp<volScalarField>& tsp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<fvMatrix<Type> > tfvm = fvm::Sp(tsp(), vf);
    tsp.clear();
    return tfvm;
}


template<class Type>
Foam::tmp<Foam::fvMatrix<Type> >
Foam::fvm::Sp
(
    const dimensionedScalar& sp,
    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();

    fvm.diag() += mesh.V()*sp.value(); // 返回系数矩阵对角线

    return tfvm;
}

template<class Type>
Foam::zeroField
Foam::fvm::Sp
(
    const zeroField&,
    GeometricField<Type, fvPatchField, volMesh>&
)
{
    return zeroField();
}

由此可以看到，fvm::Sp() 返回的类型是矩阵类型，而 fvc::Sp() 返回的是场。

大家也可能在源代码中看到 SuSp()，其中 fvc::SuSp() 的定义位于 src/finiteVolume/finiteVolume/fvc/fvcSup.C

template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> >
SuSp
(
    const volScalarField& sp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return sp*vf;
}

template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> >
SuSp
(
    const tmp<volScalarField>& tsp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return tsp*vf;
}

:question:1. 可以尝试比较一下 fvc::Sp() 与 fvc::SuSp() 有哪些不同？
:question:2. 下面一段程序是 fvm::SuSP()，你能说一说它是什么意思吗？

template<class Type>
Foam::tmp<Foam::fvMatrix<Type> >
Foam::fvm::SuSp
(
    const DimensionedField<scalar, volMesh>& susp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    tmp<fvMatrix<Type> > tfvm
    (
        new fvMatrix<Type>
        (
            vf,
            dimVol*susp.dimensions()*vf.dimensions()
        )
    );
    fvMatrix<Type>& fvm = tfvm();

    fvm.diag() += mesh.V()*max(susp.field(), scalar(0));

    fvm.source() -= mesh.V()*min(susp.field(), scalar(0))
        *vf.internalField();

    return tfvm;
}
1 个赞
4

我觉得，fvm::SuSp(A,B) 是一个智能源项。
如果 A>0 则，将系数加到矩阵的对角线上；
如果 A<0 则，将 A*B 放到矩阵方程的等号右边。