SemiImplicitSource in fvOptions
本文最后更新于：2021年1月7日 下午
本篇主要讲述具体的源项类： 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);
}
}
}上述设置，相当于求解如下方程
∂T∂t=Su+Sp⋅T∂T∂t=Su+Sp⋅T
其中 Su=0.05,Sp=0Su=0.05,Sp=0 。
反观上面对代码的分析，可知对于当前的设置，SuSu 的量纲为 TEqn 的量纲除以体积量纲，SpSp 的量为 SuSu 量纲除以 T 的量纲，这与上面给出的方程是一致的。
但是，要注意， SemiImplicitSource 有两个模式：absolute 和 specific，区别在于代码里的 VDash 的取值不一样。对于 absolute 模式，VDash = V，即所选的 cellZone 的体积；对于 specific 模式，VDash = 1.0。而代码里的 Su 和 Sp 的值都是用在fvOptions 词典文件设置的值除以 VDash 。
所以，确切地说，求解的应该是如下积分方程：
∫V∂T∂tdV−∫V∇⋅(DT∇T)dV=∫V(SuVDash+SpVDashT)dV∫V∂T∂tdV−∫V∇⋅(DT∇T)dV=∫V(SuVDash+SpVDashT)dV
其中 VDashVDash 是选定的 cellZone 的体积。
为了验证以上的推演，作了如下两个测试：
源项参数设置为
T (0.01 0)，源项作用的区域为一个体积为VDash=0.001VDash=0.001 的 cellZone，为了消除热量传递，设置 DT=0DT=0，初始整个区域的 T 均为0， 模拟时间为1 s。
按照上述的推演，如果是absolute模式，最终 1 s 时选定的 cellZone 的温度将是 T=t⋅SuVDash=1s⋅0.01k/s0.001=10kT=t⋅SuVDash=1s⋅0.01k/s0.001=10k；如果是specific模式，那么最终1 s 时选定的 cellZone 的温度将是 T=1s⋅0.01k/s1.0=0.01kT=1s⋅0.01k/s1.0=0.01k。以上结果在测试算例中得到了证实。源项参数设置为
T (1.0 2.0)，同样设置
DT=0DT=0
。这种情况下，可以先求微分方程
∂T∂t=x+yT∂T∂t=x+yT
的解，经简单计算得到
T=eyt⋅eyc−xyT=eyt⋅eyc−xy
cc
为任意常数。
若
T|t=0=0kT|t=0=0k
，则可以得到定解为
T=eyt⋅x−xyT=eyt⋅x−xy
若使用
specific模式，则根据当前的设置得到1 s时的解为
T=e2⋅1⋅1−12=3.194528kT=e2⋅1⋅1−12=3.194528k
算例测试结果为：
时间离散格式:Euler，ΔT=0.001sΔT=0.001s，T=3.20193kT=3.20193k；
时间离散格式：Euler，ΔT=0.0001sΔT=0.0001s，T=3.19527kT=3.19527k，可见减小时间步后结果与解析解吻合度提高了很多。
时间离散格式：CrankNicolson，ΔT=0.001sΔT=0.001s，T=3.19454kT=3.19454k。可见用高阶的时间离散格式，在同样时间步下能得到误差更小的结果。
同样，如果用absolute 模式，且源项设置为 T (0.001 0.002) ，应该能得到一样的结果，实际上算例测试正是如此。由此可以认为推演得到了证实。
至此，SemiImplicitSource 这个源项的核心部分就算是明了了。可是，现在测试的是非常简单的情形，如果在求解多个方程，且有多个方程里都加入了源项的情况下，该怎么给不同的方程设置不一样的源项呢？要解决这个问题，需要先理解清楚 fvOptions 的调用过程。
fvOptions 源项的调用过程
下面至下而上地来看看 fvOptions 的调用过程。
TEqn里，有一个调用 fvOptions 的语句： fvOptions(T)，上一篇讲过， fvOptions 的定义为
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 函数。为了说明这一点，先从构造函数看起
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 源项的调用途径就打通了，接下来可以继续具体分析特定的源项了。
本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议 ，转载请注明出处！