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

#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 Like

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