2016-03-20

fvOptions 之 semiImplicitSource

上篇浅析了 fvOptions 框架的结构，这篇来看一个具体的源项类： semiImplicitSource 。

先来看看这个源项代码中的关键部分。

SemiImplicitSource.C

#include "SemiImplicitSource.H"
#include "fvMesh.H"
#include "fvMatrices.H"
#include "DimensionedField.H"
#include "fvmSup.H"

template<class Type>
const Foam::wordList Foam::fv::SemiImplicitSource<Type>::volumeModeTypeNames_
(
    IStringStream("(absolute specific)")()  // 初始化 volumeModeTypeNames_，这里有两种模式， `absolute` 和 `specific` ，具体含义下面会解释。 
);

template<class Type>
typename Foam::fv::SemiImplicitSource<Type>::volumeModeType
Foam::fv::SemiImplicitSource<Type>::wordToVolumeModeType // 将字符串转换成 volumeModeType
(
    const word& vmtName
) const
{
    forAll(volumeModeTypeNames_, i)
    {
        if (vmtName == volumeModeTypeNames_[i])
        {
            return volumeModeType(i);
        }
    }

    FatalErrorIn
    (
        "SemiImplicitSource<Type>::volumeModeType"
        "SemiImplicitSource<Type>::wordToVolumeModeType(const word&)"
    )   << "Unknown volumeMode type " << vmtName
        << ". Valid volumeMode types are:" << nl << volumeModeTypeNames_
        << exit(FatalError);

    return volumeModeType(0);
}

template<class Type>
Foam::word Foam::fv::SemiImplicitSource<Type>::volumeModeTypeToWord
(
    const volumeModeType& vmtType
) const
{
    if (vmtType > volumeModeTypeNames_.size())
    {
        return "UNKNOWN";
    }
    else
    {
        return volumeModeTypeNames_[vmtType];
    }
}

template<class Type>
void Foam::fv::SemiImplicitSource<Type>::setFieldData(const dictionary& dict)
{
    fieldNames_.setSize(dict.toc().size());
    injectionRate_.setSize(fieldNames_.size());

    applied_.setSize(fieldNames_.size(), false);

    label i = 0;
    forAllConstIter(dictionary, dict, iter)
    {
        fieldNames_[i] = iter().keyword();
        dict.lookup(iter().keyword()) >> injectionRate_[i];
        i++;
    }

    // Set volume normalisation
    if (volumeMode_ == vmAbsolute)
    {
        VDash_ = V_;
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class Type>
Foam::fv::SemiImplicitSource<Type>::SemiImplicitSource
(
    const word& name,
    const word& modelType,
    const dictionary& dict,
    const fvMesh& mesh
)
:
    option(name, modelType, dict, mesh),
    volumeMode_(vmAbsolute), volumeMode 初始值为 vmAbsolute，也就是字典里的 absolute
    VDash_(1.0), // VDash 初始值为 1.0
    injectionRate_()
{
    read(dict);
}


template<class Type>
void Foam::fv::SemiImplicitSource<Type>::addSup // 这个是最关键的函数，求解器里的 fvOptions(T)，最终就是转换为调用这个函数。
(
    fvMatrix<Type>& eqn,
    const label fieldI
)
{
    if (debug)
    {
        Info<< "SemiImplicitSource<" << pTraits<Type>::typeName
            << ">::addSup for source " << name_ << endl;
    }

    // psi 表示方程中的未知量，比如，eqn(fvm::ddt(T))，则psi其实就相当于T，其量纲也与T的量纲一致。
    const GeometricField<Type, fvPatchField, volMesh>& psi = eqn.psi();

// 创建一个场 Su，其量纲为方程eqn的量纲除以体积的量纲。注意，经测试，假设eqn为（fvm::ddt(T)），则eqn的量纲为k.m^3/s。
    DimensionedField<Type, volMesh> Su
    (
        IOobject
        (
            name_ + fieldNames_[fieldI] + "Su",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensioned<Type>
        (
            "zero",
            eqn.dimensions()/dimVolume,
            pTraits<Type>::zero
        ),
        false
    );

    // 这一句的意思是，将属于cells_这个集合的网格的Su赋值为fvoptions里所设置的第一个参数的值除以体积VDash。这里的VDash,如果模式为absolute，则值为cells_这个集合的网格体积之和，如果模式为specific，则其值为1. 
    // UUIndirectList<Type>(Su, cells_)这一句是利用Su和cells为参数，构建一个UUIndirectList类的临时对象，并调用这个类的重载的“=”操作符对Su进行重新赋值。
    
    UIndirectList<Type>(Su, cells_) = injectionRate_[fieldI].first()/VDash_;

    DimensionedField<scalar, volMesh> Sp
    (
        IOobject
        (
            name_ + fieldNames_[fieldI] + "Sp",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensioned<scalar>
        (
            "zero",
            Su.dimensions()/psi.dimensions(),
            0.0
        ),
        false
    );
    
    UIndirectList<scalar>(Sp, cells_) = injectionRate_[fieldI].second()/VDash_;

    // fvMatrix<Type> 类中对“+=”操作符进行了重载，所以，eqn与Su的相加，相当于eqn+Su*mesh.V()，要不然eqn与Su的量纲不一致。
    eqn += Su + fvm::SuSp(Sp, psi);
}

下面用一个例子来说明 semiImplicitSource 的作用。前提到 scalarTransportFoam 是使用 fvOptions 的一个最简单的求解器，这里对该求解器进一步简化，只保留瞬变项，对流和扩散项都删去，来验证 semiImplicitSource 的作用。
修改之后的 TEqn 为

tmp TEqn
(
    fvm::ddt(T)  == fvOptions(T)
);

然后，fvOptions 词典文件的设置如下：

firstHeatSource
{
    type scalarSemiImplicitSource;
    active          true;
    selectionMode   cellZone;
    cellZone        boxSourceZone;
    scalarSemiImplicitSourceCoeffs 
    {
        volumeMode absolute;
        injectionRateSuSp 
        {
            T (0.05 0);
        }
    }
}

上述设置，相当于求解如下方程
$$
\frac{\partial T}{\partial t}=S_u + S_p\cdot T
$$
其中 $S_u=0.05, S_p = 0$ 。
反观上面对代码的分析，可知对于当前的设置，$S_u$ 的量纲为 TEqn 的量纲除以体积量纲，$S_p$ 的量为 $S_u$ 量纲除以 T 的量纲，这与上面给出的方程是一致的。
但是，要注意， SemiImplicitSource 有两个模式：absolute 和 specific，区别在于代码里的 VDash 的取值不一样。对于 absolute 模式，VDash = V，即所选的 cellZone 的体积；对于 specific 模式，VDash = 1.0。而代码里的 Su 和 Sp 的值都是用在fvOptions 词典文件设置的值除以 VDash 。
所以，确切地说，求解的应该是如下积分方程：
$$
\int_V \frac{\partial T}{\partial t} dV - \int_V \nabla \cdot (D_T \nabla T) dV = \int_V (\frac{S_u}{V_{Dash}}+\frac{S_p}{V_{Dash}}T )dV
$$
其中 $V_{Dash}$ 是选定的 cellZone 的体积。
为了验证以上的推演，作了如下两个测试：

源项参数设置为T (0.01 0)，源项作用的区域为一个体积为$V_{Dash}=0.001$ 的 cellZone，为了消除热量传递，设置 $D_T = 0$，初始整个区域的 T 均为0，模拟时间为1 s。
按照上述的推演，如果是 absolute 模式，最终 1 s 时选定的 cellZone 的温度将是 $T=t\cdot \frac{S_u}{V_{Dash}}=1s\cdot\frac{0.01 k/s}{0.001}=10k$；如果是 specific 模式，那么最终1 s 时选定的 cellZone 的温度将是 $T=1s\cdot\frac{0.01 k/s}{1.0}=0.01k$。以上结果在测试算例中得到了证实。
源项参数设置为T (1.0 2.0)，同样设置 $D_T = 0$。这种情况下，可以先求微分方程
$$\frac{\partial T}{\partial t}=x+yT$$
的解，经简单计算得到
$$T=\frac{e^{yt}\cdot e^{yc}-x}{y}$$
$c$ 为任意常数。
若 $T|\,_{t=0}=0k$，则可以得到定解为
$$T=\frac{e^{yt}\cdot x-x}{y}$$
若使用specific 模式，则根据当前的设置得到1 s时的解为
$$T=\frac{e^{2\cdot 1}\cdot 1 -1}{2}=3.194528 k$$
算例测试结果为：
时间离散格式:Euler，$\Delta T=0.001s$，$T=3.20193k$；
时间离散格式：Euler，$\Delta T=0.0001s$，$T=3.19527k$，可见减小时间步后结果与解析解吻合度提高了很多。
时间离散格式：CrankNicolson，$\Delta T=0.001s$，$T=3.19454k$。可见用高阶的时间离散格式，在同样时间步下能得到误差更小的结果。

同样，如果用absolute 模式，且源项设置为 T (0.001 0.002) ，应该能得到一样的结果，实际上算例测试正是如此。由此可以认为推演得到了证实。

至此，SemiImplicitSource 这个源项的核心部分就算是明了了。可是，现在测试的是非常简单的情形，如果在求解多个方程，且有多个方程里都加入了源项的情况下，该怎么给不同的方程设置不一样的源项呢？要解决这个问题，需要先理解清楚 fvOptions 的调用过程。

fvOptions 源项的调用过程

下面至下而上地来看看 fvOptions 的调用过程。

TEqn里，有一个调用 fvOptions 的语句： fvOptions(T)，上一篇讲过， fvOptions 的定义为

1	fv::IOoptionList fvOptions(mesh);

这就很明显了：建立一个 IOoptionList 类的对象 fvOptions。由此可知，求解器里的 fvOptions 是一个对象的名字，因此 fvOptions(T) 这种用法也只可能是对象调用类中重载的 () 运算符了。从前面的分析，可知 IOoptionList 本身很简单，仅是作为一个接口来用的，所以 () 与算符的重载要去其父类中去找。 IOoptionList 的作用是，从 “constant”（优先）或者 “system” 目录读取 fvOptions 文件，并作为 IOobject 类的对象传递给父类 optionList （从构造函数的成员初始化列表 optionList(mesh, *this) ）

接下来，就该进入 optionList 类了。这个类里，所有可能出现在求解器代码里函数都有了，包括 correct ， constrain ， makeRelative ， makeAbsolute ， relative 以及 () 运算符的重载。但是，注意这里并没有具体的代码实现，而是通过类似

forAll(*this, i)
{
   this->operator[](i).makeAbsolute(phi);
}

调用其他地方的函数。 optionList 的一个重要使命是，统计 fvOptions 文件里定义了多少个源项，并将每一个源项都作为一个储存起来，然后再根据词典的内容创建特定的源项。核心在于 reset 函数。为了说明这一点，先从构造函数看起

1	reset(optionsDict(dict));

可见构造函数里调用了 reset 函数，并且用 optionsDict 函数的返回值作为 reset 函数的参数。前文讲到， IOoptionList 类将 fvOptions 文件的内容以 IOobject 的形式传递给父类 optionList，所以，这里的参数 dict 可以理解为就是fvOptions 文件的内容。optionsDict 函数的代码如下：

const Foam::dictionary& Foam::fv::optionList::optionsDict
(
    const dictionary& dict
) const
{
    if (dict.found("options"))
    {
        return dict.subDict("options");
    }
    else
    {
        return dict;
    }
}

可见，这个函数去 fvOptions 文件里查找关键字 options，如果找到，就将 options 关键字对应的 subDict 内容返回，否则直接返回 fvOptions 文件的内容。举例说，形如

firstHeatSource
{
    type scalarSemiImplicitSource;
    active          true;
    selectionMode   cellZone;
    cellZone        boxSourceZone;
    scalarSemiImplicitSourceCoeffs 
    {
        volumeMode absolute;
        injectionRateSuSp 
        {
            T (0.001 0);
        }
    }
}

的，直接返回，因为这就构成了一个 dictionary；而形如

options
{
    massSource1
    {
        type            scalarSemiImplicitSource;
        $injector1;

        scalarSemiImplicitSourceCoeffs
        {
            volumeMode      absolute;
            injectionRateSuSp
            {
                thermo:rho.air     (1e-3 0); // kg/s
            }
        }
    }

    momentumSource1
    {
        type            vectorSemiImplicitSource;
        $injector1;

        vectorSemiImplicitSourceCoeffs
        {
            volumeMode      absolute;
            injectionRateSuSp
            {
                U.air           ((0 -1e-2 0) 0); // kg*m/s^2
            }
        }
    }

    energySource1
    {
        type            scalarSemiImplicitSource;
        $injector1;

        scalarSemiImplicitSourceCoeffs
        {
            volumeMode      absolute;
            injectionRateSuSp
            {
                e.air      (500 0); // kg*m^2/s^3
            }
        }
    }
}

的，则将options下的每一个 subDict 作为 dictionary 返回。注意，这里的 dictionary 类可以理解为一个容器，每一个

xxxxx
{
    ......
}

都可以作为容器里的一个成员，容器的容量（size）等于总的成员数，每一个成员，其实就对应一个源项。
于是，我们知道 optionsDict 返回了一个有一定数目成员的容器。再来看 reset 函数

void Foam::fv::optionList::reset(const dictionary& dict)
{
    // Count number of active fvOptions
    label count = 0;
    
    // 遍历 dict 容器，确定其成员的数目，即确定定义了几个源项。
    forAllConstIter(dictionary, dict, iter) 
    {
        if (iter().isDict())
        {
            count++;
        }
    }
    this->setSize(count); // setSize 是 PtrList 类的成员，顾名思义，PtrList 是一个 List。PtrList<option> 类的成员是 options 类的对象。
    
    label i = 0;
    // 遍历 dict 容器，根据每一个 dict 容器的成员来建立对应的 option 类的对象，这通过调用 option 类的 New 函数来实现。这是使用 RuntimeSelection 机制的类的很常规的做法。
    forAllConstIter(dictionary, dict, iter)
    {
        if (iter().isDict())
        {
            const word& name = iter().keyword(); // keyword 返回的是类似 energySource1 这样的，是这个源项的一个名字
            const dictionary& sourceDict = iter().dict();

            this->set // set 函数，肯定毫无疑问也是从 PtrList 类中继承而来的
            (
                i++,
                option::New(name, sourceDict, mesh_) // 调用 option 类的 New 函数
            );
        }
    }
}

注意，最重要的是， reset 里实现了对父类 PtrList<option> 的初始化。
理解了这些，再来看 correct 以及 () 操作符重载中的代码，就好理解了：

// 本来 *this 应该是 optionList 类的指针，这里先作个隐式转换，转成 PtrList<option> 类的指针。前面reset函数已经对 PtrList<option> 进行了初始化，使其读入了每一个源项。所以，这里的循环就是对每一个源项进行循环，然后调用对应的函数。operator[]肯定也是在PtrList类中定义的，i 指的是 PtrList 类的第 i 个成员，这里 PtrList 的每一个成员都是一个 option 类的对象，所以，makeAbsolute 函数是定义在 option 类中的函数。
forAll(*this, i)
{
   this->operator[](i).makeAbsolute(phi);
}

根据上面的理解，一个很自然的推论是，定义在 fvOptions 文件中的源项，其作用是叠加的。也就是说，上述的 fvOptions(T)，对定义在 fvOptions 文件中的每一个源项，都会调用一次。经测试，

 firstHeatSource
{
    type scalarSemiImplicitSource;
    active          true;
    selectionMode   cellZone;
    cellZone        boxSourceZone;
    scalarSemiImplicitSourceCoeffs 
    {
        volumeMode absolute;
        injectionRateSuSp 
        {
            T (0.001 0);
        }
    }
}
        
secondHeatSource
{
    type scalarSemiImplicitSource;
    active          true;
    selectionMode   cellZone;
    cellZone        boxSourceZone;
    scalarSemiImplicitSourceCoeffs 
    {
        volumeMode absolute;
        injectionRateSuSp 
        {
            T (0.0 0.002);
        }
    }
}

与

secondHeatSource
{
    type scalarSemiImplicitSource;
    active          true;
    selectionMode   cellZone;
    cellZone        boxSourceZone;
    scalarSemiImplicitSourceCoeffs 
    {
        volumeMode absolute;
        injectionRateSuSp 
        {
            T (0.0001 0.002);
        }
    }
}

的作用是一样的，这证实了源项的作用确实是叠加的。
此外，还要注意一点，那就是 optionList 类中重载的 () 运算符，其实是在调用 option 类中的 addSup 函数，并且其返回值是 fvMatrix 类的对象。

template<class Type>
Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
(
    GeometricField<Type, fvPatchField, volMesh>& fld,
    const word& fieldName
)
{
    checkApplied();

    const dimensionSet ds = fld.dimensions()/dimTime*dimVolume;

    tmp<fvMatrix<Type> > tmtx(new fvMatrix<Type>(fld, ds));
    fvMatrix<Type>& mtx = tmtx();

    forAll(*this, i)
    {
        option& source = this->operator[](i);

        label fieldI = source.applyToField(fieldName);

        if (fieldI != -1)
        {
            source.setApplied(fieldI);

            if (source.isActive())
            {
                if (debug)
                {
                    Info<< "Applying source " << source.name() << " to field "
                        << fieldName << endl;
                }

                source.addSup(mtx, fieldI);
            }
        }
    }
    return tmtx;
}

理解了以上这些，就可以进入 option 类了。
option 类是所有具体的源项类的基类，这个类里处理了所有源项都需要处理的部分，比如确定起始时间，选择源项作用的区域，这些都是在构造函数里完成的，主要是通过调用 setSelection 和 setCellSet 两个函数。这里需要注意的是数据成员 cells_， V_ 以及 fieldNames_，分别定义源项作用区域的网格id （这里源项的作用区域是以 cell 为基础来指定的，即便是 points 模式，实际上源项作用的区域仍然是 points 所在的 cell。），选定区域的体积以及需要开启源项作用的场名（有时候，求解器的多个 Eqn 里有fvOptions，但是实际算例中只想针对特定的场开启源项，这可以通过指定 fieldNames_ 来实现）。

注意 fieldNames 在 option 类中并没有初始化，需要在具体的源项类中指定。此外，makeRelative 等在求解器中实际调用的函数，在 option 类中也并没有进行具体的实现（但不是声明为纯虚函数，仅仅是函数体为空的而已，这里的结果不适合用纯虚函数）

以 SemiImplicitSource 为例，主要去看 addSup 函数，这个函数，关键的一个参数是fieldI，这个参数，指的是 fieldName_ 这个List的 applyToField 函数的返回值，用来判断一个 field 是否要启用源项。fieldName_ 这个List，是在 SemiImplicitSource 类的 setFieldData 函数中赋值的，通过读取 SemiImplicitSourceCoeffs 中的参数来决定。但是这个 fieldName_ 的确定方法根据不同的类有不同的做法，要具体分析。
到此，fvOpptions 源项的调用途径就打通了，接下来可以继续具体分析特定的源项了。

Giskard's CFD Learning Tricks

CFD and Scientific Computing

fvOptions 之 semiImplicitSource

fvOptions 源项的调用过程