From 9c0d7df0bc6619f7db35e8d4d59aa2c4372225a5 Mon Sep 17 00:00:00 2001 From: luofq Date: Thu, 11 Feb 2021 12:29:41 +0800 Subject: [PATCH] =?UTF-8?q?add=2057.Under=5Fthe=5Fhood=5Fof=5FOpenFOAM.md.?= =?UTF-8?q?=20=E6=8F=90=E4=BA=A457=E7=AB=A0=E7=9A=84=E5=88=9D=E7=A8=BF?= =?UTF-8?q?=EF=BC=8Cluofq5?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- 57.Under_the_hood_of_OpenFOAM.md | 1625 ++++++++++++++++++++++++++++++ 1 file changed, 1625 insertions(+) create mode 100644 57.Under_the_hood_of_OpenFOAM.md diff --git a/57.Under_the_hood_of_OpenFOAM.md b/57.Under_the_hood_of_OpenFOAM.md new file mode 100644 index 0000000..7f73dcb --- /dev/null +++ b/57.Under_the_hood_of_OpenFOAM.md @@ -0,0 +1,1625 @@ +### 57 Under the hood of OpenFOAM +本章提供了一些简短的代码示例，用以解释OpenFOAM在特定场景下的行为方式。这些例子均来源于本手册的其他部分，其中某些程序例子的源代码已经在其他地方解释过了。 + +#### 57.1 求解算法 +请参照第Ⅴ部分的第41章、42章及44章。 + +#### 57.2 命名空间 +##### 57.2.1 常量 +物理场中包含了很多常量。因此，在一个特定区域对物理或数学的常量进行定义是一个不错的选择。OpenFOAM在命名空间Foam：：constant中提供了常量。这些预先定义的常量按类别进行分组，例如： +关于电磁场的 + -mu0表示真空中的电磁渗透率 + -epsilon0表示真空中的电子渗透率 +关于物理化学的 + -R表示广义气体常态 +数学的 + -pi表示π + -e表示欧拉数 + +在列表388中展示了如何从源代码中获取常量pi。列表389展示了OF-2.2.x中定义的所有常量。从提高计算性能的角度出发，预先定义像π这样经常使用的常量具有重要的意义。需要注意的是，除以2.0通常采用乘以0.5来替代。从数学上来说，这些操作是相等的，但从计算成本的角度上出发，相比浮点数的除法操作，浮点数乘法操作具有更快的速度。 +此外，需要注意的是OpenFOAM没有额外自定义e和π，它采用的是系统库提供的常量。GNU C库提供的数学变量可参见： +http://www.gnu.org/software/libc/manual/html_node/Mathematical-Constants.html。因此，e和π是通过访问M_E和M_PI进行定义的。 +需要进一步注意的是常量在定义时会带有const描述，这也是C和C++语言中定义场量的唯一稳健方式。 + +``` + scalar foo = constant :: mathematical ::pi; +``` +列表388：一个无用的代码示例，在OpenFOAM源代码中访问π + +``` +const scalar e(M_E); +const scalar pi(M_PI); +const scalar twoPi (2*pi); +const scalar piByTwo (0.5*pi); +``` +列表389：mathematicalConstants.H提供的数学常量 + +在OpenFOAM-extend中，数学常量的访问方式也相似。唯一的不同是相应的命名空间不是 constant::mathematical，而是mathematicalConstants。这是因为OpenFOAM-extend版本大部分是基于OpenFOAM-1.6版本发展起来的。 + +#### 57.3 从字典中查找关键字 +字典中一般存在两种关键字，分别是强制关键字和可选关键字。 +##### 57.3.1 强制关键字 +当强制关键字在字典中未被搜索到时，OpenFOAM会报出一段错误信息并中止运行。 +列表390展示了读取三个强制关键字的操作。其中的lookup（）函数将在列表391中解释。 + +``` +#include "readTimeControls.H" +int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr"))); int nAlphaSubCycles(readInt(pimple.dict().lookup("nAlphaSubCycles"))); +Switch correctAlpha(pimple.dict().lookup("correctAlpha")); +``` +列表390：readTwoPhaseEulerFoamControls.H的内容 + +###### 代码 +列表391的第32行显示，lookup（）函数直接调用了lookupEntry（）函数的返回值。该函数同样调用了另一个函数（lookupEntryPtr（）），并执行错误处理。错误处理机制清楚地显示了OpenFOAM将在找不到关键字时中止运行（第19行）。 + +``` +const Foam:: entry& Foam:: dictionary :: lookupEntry +( +const word& keyword , +bool recursive , +bool patternMatch +) const +{ +const entry* entryPtr = lookupEntryPtr(keyword , recursive , patternMatch); +if (entryPtr == NULL) +{ +FatalIOErrorIn ( +"dictionary :: lookupEntry(const word&, bool , bool) const", +*this +) +<< "keyword " << keyword << " is undefined in dictionary " +<< name() +<< exit(FatalIOError); +} +return *entryPtr; +} +Foam:: ITstream& Foam:: dictionary :: lookup +( +const word& keyword , +bool recursive , +bool patternMatch +) const +{ +return lookupEntry(keyword , recursive , patternMatch).stream (); +} +``` +列表391：dictionary.C的部分内容 + +##### 57.3.2 可选关键字 +从字典中读取可选关键字的函数通常会指定一个默认值，当这个关键字在字典中不存在时则该默认值生效。 +列表392展示了读取三个可选关键字的操作。该读取函数被调用时输入了两个参数，第一个是关键字，第二个是相应的默认值。如果函数lookupOrDefault（）没有找到接口，默认值将被返回。 +``` +const bool adjustTimeStep = +runTime.controlDict ().lookupOrDefault("adjustTimeStep", false); +scalar maxCo = +runTime.controlDict ().lookupOrDefault ("maxCo", 1.0); +scalar maxDeltaT = +runTime.controlDict ().lookupOrDefault ("maxDeltaT", GREAT); +``` +列表392：readTimeControls.H的内容。 + +###### 代码： +列表393展示了函数lookupOrDefault（）的定义。这个函数同样调用了另一个函数来查找关键字，实际上它查找的是字典上分配到相应关键字的值，并且进入了条件判断分支。如果该关键字被找到，对应的值将被返回；如果该关键字没有找到，默认值将被返回（第18行）。 +在列表393中，该函数定义了四个输入参数。然后，在列表392中该函数只被两个输入参数调用。 +该矛盾的解决方案在dictionary.H中的函数定义中可以看到。函数在列表394中被声明，在第6行和第7行中，两个输入变量的默认值被具体定义了。因此，该函数也可以只被另外两个没有相应默认值的输入参数进行调用。如果该函数被它的所有输入参数调用，传值进来的参数将默认值覆盖掉。 +当声明一个使用默认值参数的函数，那些不使用默认值的参数必须排在使用默认值的参数之前。否则，将产生歧义。 + +``` +template +T Foam:: dictionary :: lookupOrDefault +( +const word& keyword , +const T& deflt , +bool recursive , +bool patternMatch +) const +{ +const entry* entryPtr = lookupEntryPtr(keyword , recursive , patternMatch); +if (entryPtr) +{ return pTraits (entryPtr ->stream ()); } +else +{ +return deflt; +} +} +``` +列表393：dictionaryTemplates.C的部分内容 + +``` +template T lookupOrDefault +( +const word&, +const T&, +bool recursive=false , +bool patternMatch=true +) +const; +``` +列表394：dictionary.H的部分内容 + +##### 57.3.3 取代强制关键字 +在一些例子中，一个基类可能要求关键字的存在，它的一部分子类中使用了该关键字，另一些子类则不使用。因此产生这样一种情况，每一组的关键字-默认值都同时被描述，但它是否被使用将取决于使用者的决定。 +然而，不管这些关键字在子类中的使用性质，他们在基类中都是必须定义的，所以将他们设置为可选参数。一个继承而来的子类无法改变相应基类的性质。在某些极端例子中，即使一组关键字-值完全不起作用，但还是要求提供定义。 +如何在基类及其子类中对必要数据进行合理分配，不总是那么简单。除了以上两种明显局限性的案例之外（一是数据只在某个特定的派生类中需要，二是数据在所有派生类中都被需要），其他情况中的数据如何分布将取决于软件设计人员。 +如果有足够数量的相似派生类，那么可以建立一种中间基类。有了中间基类之后，由其衍生出来的派生类所共有的数据可以移交给中间基类。然而，特定模型的类继承特性也一定总能准确反应该类模型在数据使用上的内在逻辑。 +在图130可以看到这样一个中间基类例子。类图中的数据成员分别被各自的类视为强制入口。基类InjectionModel读取了一个标量massTotal，该标量决定了有多少个粒子被注入。然而，在一些派生的注入模型中允许直接对粒子数量进行描述，因此该强制值massTotal将不被使用。在这些例子中，injection模型会报错，提示用户massTotal不会有效果。 + +图130：一些拉格朗日粒子注入模型的类图。采用中间基类来减少代码的重复性，即紧密相关的注入模型代码。 + + +#### 57.4 OpenFOAM 中的特定数据类型 +##### 57.4.1 Switch数据类型 +在字典中拥有大量的开关设置，用以激活或取消激活某一个特征，列表395展示了定义所有有效值的部分源代码。该代码中每一个Switch的值只允许是true或false，即Switch类是一种用作布尔值判断的数据类型。然而，在字典中的Switch可以拥有更多的值，假设他们共享同一种决策。人类语言通常拥有多种方式来回答一个yes-no问题，这可能是允许通过一系列范围值来定义Switch的原因。 + +``` +// NB: values chosen such that bitwise ’&’ 0x1 yields the bool value +// INVALID is also evaluates to false , but don’t rely on that +const char* Foam:: Switch ::names[Foam:: Switch :: INVALID +1] = { +"false", "true", +"off", "on", +"no", "yes", +"n","y", +"f","t", +"none","true", // is there a reasonable counterpart to "none"? +"invalid" +}; +``` +列表395：Switch.C的部分内容 + +列表396展示了一个关于Switch数据类型如何在代码中被使用的例子，该代码读取transportProperties字典。如果没有有效的testSwitch数据入口被接入，那么该switch的值将被设置为false。注意该函数lookupOrDefault()的第二个参数，它读取Switch(false)。它的含义是，一个类型是Switch的对象被创建，并且是通过将布尔值false传给类Switch的构造函数来创建的。之后，该Switch类型的新对象在必要时被用作名为testSwitch的switch默认值。 + +``` +Switch testSwitch(transportProperties.lookupOrDefault ("testSwitch", witch(false))); +``` +列表396：Switch数据类型的使用例子 + +##### 57.4.2 label数据类型 +当查阅类似42.2章节的求解算法时，可以看到计数器的存在。OpenFOAM使用了一种名叫label的数据类型来定义计数器，请看列表266. + +###### 以下内容应用在OpenFOAM-3.0之前的版本 +计数器中最重要的数据类型是整型数据int，列表397包含了文件label.H的部分代码，这正是int数据类型的定义。取决于系统或编译参数，label可以是int、long、long long类型的某一种。 +列表397展示了label的定义，此时int被用作是隐含的数据类型。 + +``` +namespace Foam +{ +typedef int label; +static const label labelMin = INT_MIN; +static const label labelMax = INT_MAX; +inline label readLabel(Istream& is) +{ +return readInt(is); +} +} // End namespace Foam +``` +列表397 label.H 的部分内容 + +###### 以下内容应用在OpenFOAM-3.0之后的版本 +OpenFOAM在数据写入过程中提供了两种数据长度的label类型：32位和64位。从底层的代码来说，这归结为使用int32_还是int64_t。由于int，long和long long数据类型仅保证一个最小的数据长度，因此采用固定宽度整数int32_t或int64_t用作label数据类型的基础。 +label数据类型的长度可以在编译之前通过编译选项WM_LABEL_SIZE进行指定，即可以指定为32或64。借助这里的编译器选项和之前的一些预处理器宏，可以实现整数类型的构建。 + +``` +#define INT_ADD_SIZE(x,s,y) x ## s ## y +#define INT_ADD_DEF_SIZE(x,s,y) INT_ADD_SIZE(x,s,y) +#define INT_SIZE(x,y) INT_ADD_DEF_SIZE(x,WM_LABEL_SIZE ,y) +#if WM_LABEL_SIZE != 32 && WM_LABEL_SIZE != 64 +#error "label.H: WM_LABEL_SIZE must be set to either 32 or 64" +#endif +namespace Foam +{ +typedef INT_SIZE(int , _t) label; +/* code removed for brevity */ +} +``` +列表398：label.H的部分内容 + +在列表398中，可以看到一些帮助宏以及健全性检查。如果label的数据长度不是32或64，则健全性检查会报错。列表398的最后一行是label数据类型的实际类型定义。宏INT_SIZE通过多个帮助程序宏和编译器选项WM_LABEL_SIZE，以实现固定宽度的整数数据类型int32_t或int64_t的构建。INT_ADD_SIZE宏中的##运算符用于连接宏的参数。 + +##### 57.4.3 scalar数据类型 +与整数数据类型相似，OpenFOAM针对浮点数也有相应的特定数据类型，即scalar。Scalar通过编译设置变成float或double类型。在2018年，一种拓展精度的double类型被增加为另一种选项。 +其中基本的浮点数数据类型float是32位宽度，即它的二进制形式会占用32bits的内存。而Double类型会使用两倍的bits，因此也叫做双精度。计算中使用float或double经常被称为单精度或双精度。 +OpenFOAM在为相应编译器选项WM_PRECISION_OPTIOON命名时遵循下列的准则。对单精度和双精度分别采用SP和DP取值。对于最近版本的long double类型，可以通过值LP进行选择。 +采用整数或浮点数等固定宽度的数据类型，可以使得应用更加便捷，即某个计算机的二进制结果文件可以被另一台计算机打开和处理，只要求两台计算机在安装OpenFOAM时采用相同的数据宽度，如32位整数或双精度浮点数等。 +列表399展示了floatScalar类型的定义，其中提及了标准数据类型float。类似地，同样有一个叫做doubleFloat的数据类型，它与标准数据类型double相关。通过引入类型scalar， +OpenFOAM的源代码与实际使用的浮点类型无关。 而floatScalar和doubleScalar类型是中间媒介，用来帮助用户在编译之前进行选择。 +floatScalar和doubleScalar类型主要由负责输入和输出的低层次代码调用。没有应用程序（求解器和工具）或湍流模型等高级库会直接使用这些数据类型。所有高级代码均从实际的浮点数据类型抽象而来。 + +``` +namespace Foam +{ +typedef float floatScalar; +/* code removed for brevity */ +``` +列表399：floatScalar.H的部分内容 + +列表400展示了scalar.H的相关代码行，其中根据WPM_SP或WM_DP的赋值来定义类型scalar。取决于定义为WM_SP或WM_DP，类型scalar将引用floatScalar或doubleScalar。 + +``` +#if defined(WM_SP) +// Define scalar as a float +namespace Foam +{ +typedef floatScalar scalar; +/* code removed for brevity */ +} +#elif defined(WM_DP) +// Define scalar as a double +namespace Foam +{ +typedef doubleScalar scalar; +/* code removed for brevity */ +} +#endif +``` +列表400：scalar.H的部分内容 + +##### 57.4.4 tmp<>数据类型 +有一种针对所有临时数据的特殊类。由于在C++中没有内存管理，程序员必须删除未使用的变量。作者认为所有存储临时数据的tmp类都是旨在区分临时变量与其他变量。 +tmp类使用一了种称为泛型编程的技术 + +##### 57.4.5 IOobject数据类型 +IOobject是处理各种数据结构行为的类。尽管没有IOobject类型的变量，但是理解此类的某些内容将有助于理解OpenFOAM某些特定的内容。 +列表401和402展示了来自求解器twoPhaseEulerFoam源代码的一些例子。如IOobject类如何用于创建场以及字典对象。 +在列表401中，两个volScalarField类型的变量被创建。 volScalarField类的构造函数接收了两个参数。在这两种情况下，接收的第一个参数都是IOobject。 +继续查看IOobject构造函数调用的参数。第一个参数是IOobject的名称，最后两个参数是读和写标志。 +对于场变量alpha1和alpha2，读和写标志是不同的。alpha1是在应用程序开始时进行读取的场，并且由于写标志的存在，无论是否进入数据写入阶段场alpha1都会被自动写入磁盘。相反，场alpha2不会被写入磁盘，并且应用程序也不会尝试读取它。 +IOobject的名称也是应用操作时的文件名称。因此场alpha1将以名为alpha1的文件形式写入磁盘。同样，当应用程序尝试读取场alpha1时，它也会尝试从文件alpha1中进行读取。 + +``` +volScalarField alpha1 ( +IOobject +( +"alpha1", +runTime.timeName (), +mesh , +IOobject ::MUST_READ , +IOobject :: AUTO_WRITE +), +mesh +); +volScalarField alpha2 +( +IOobject +( +"alpha2", +runTime.timeName (), +mesh , +IOobject ::NO_READ , +IOobject :: NO_WRITE +), +scalar (1) - alpha1 +); +``` +列表401：createFields.H中关于体积分数场的定义 + +列表402显示了IOdictionary是如何被定义的。Iodictionary类的构造函数同样接收一个Ioobject类作为输入参数。并且，Ioobject的名字同样作为应用程序进行读入操作时相应文件的名字。注意到其中的读入标志。该标志使得应用程序可以检查文件在运行时间内是否发生修改，如果发生则会再次读入。 + +``` +IOdictionary ppProperties +( +IOobject ( +"ppProperties", +runTime.constant (), +mesh , +IOobject :: MUST_READ_IF_MODIFIED , +IOobject :: NO_WRITE +) +); +``` +列表402：readPPProperties.H中关于字典的定义 + +###### 读&写标志 +在构造函数中，所谓的读写标志是作为输入参数被提供，可参见列表402的第8和9行。 +列表403展示了可用的读/写标志。OpenFOAM-2.0.0引入了标记MUST_READ_IF_MODIFIED。 +目前可用的读取标志有相当大的灵活性。 + +``` +//- Enumeration defining the valid states of an IOobject +enum objectState +{ +GOOD , +BAD +}; +//- Enumeration defining the read options +enum readOption +{ +MUST_READ , +MUST_READ_IF_MODIFIED , +READ_IF_PRESENT , +NO_READ +}; +//- Enumeration defining the write options +enum writeOption +{ +AUTO_WRITE = 0, +NO_WRITE = 1 +}; +``` +列表403：IOobject.H中关于IOobject的对象状态和读/写标志的定义 + +###### Pitfall: Solving for a NO_READ field +作者在修改求解器时无意中曾发现到一个有意思的错误。该错误可以归结为复制和粘贴的错误。但是，作者还是希望分享这一经验。 +如果想将标量输运方程式植入到现有求解器，例如在想要求解一个volScalarField场，则必须先创建该场。可采用复制相关代码的方式进行创建，也有大量的代码文件可供选择，列表404展示了一个示例。该场的名称和写入标志都已相应更改。由于我们要创可视化的图像，因此需要将写标志设置为AUTO_WRITE。但是，没有注意到读入标志的情况。 + +``` +volScalarField T +( +IOobject +( +"T", +runTime.timeName (), +mesh , +IOobject ::NO_READ , +IOobject :: AUTO_WRITE +), +mesh , +dimensionedScalar("zero", dimensionSet (0, 0, 0, 0, 0), 0.0) +); +``` +列表404：通过IOobject::NO_READ读入标志来创建一个场 + +当创建自己的场T并组成该场的输运方程（TEqn）后，想求解这个输运方程式。但是，调用Teqn.solve（）会导致一些意外结果。列表405展示了OpenFOAM发出的错误消息。 + +``` +--> FOAM FATAL ERROR: +valueInternalCoeffs cannot be called for a calculatedFvPatchField +on patch inlet of field T in file "/home/user/OpenFOAM/user -2.3.x/run/foo/case /0/T" +You are probably trying to solve for a field with a default boundary condition. +From function calculatedFvPatchField :: valueInternalCoeffs(const tmp&) +const +in file fields/fvPatchFields/basic/calculated/calculatedFvPatchField.C at line 154. +FOAM exiting +``` +列表405：当尝试求解一个非读入场时OpenFOAM提示的错误信息 + +咋一看，该消息似乎不可思议，因为我们反复检查了文件T中的边界条件。同样，改变边界条件也不会出现不一样的结果。 +该错误消息提示，想要求解的是一个具有默认边界条件的场。这毫无疑问是正确的，但我们需要找出原因。最终发现是由于使用NO_READ标志创建该场，相当于没有边界条件被提供。因此，OpenFOAM会指定默认边界条件。如果从磁盘文件中读取的boundaryField字典未作处理，得出结果也不变。 + +###### 进一步的问题 +仅改变列表404中的读取标志并不能完全解决问题。将读取标志从NO_READ更改为MUST_READ会产生与列表405相同的错误消息。 +其原因是列表404中构造函数所调用的参数。如果要从磁盘中读取场，就不能传递场值（列表404中的第12行）。 +为了使修改后的求解器正常工作，需要删除列表404第12行中传递的参数。OpenFOAM的开发人员已经预见到了这种情况，因此当对带有MUST_READ或MUST_READ_IF_MODIFIED读取标志的构造函数进行传值时，OpenFOAM会发出警告消息，请参见列表406。 + +``` +--> FOAM Warning : From function GeometricField :: readIfPresent () +in file /home/user/OpenFOAM/OpenFOAM-2.3.x/src/OpenFOAM/lnInclude/GeometricField.C at line 108 +read option IOobject :: MUST_READ or MUST_READ_IF_MODIFIED suggests that +a read constructor for field T would be more appropriate. +``` +列表406：由于构造函数参数中的读取标志或初始值不当引起的OpenFOAM警告消息 + +##### 57.4.6 关于随机的东西 +OpenFOAM具有随机数生成器（RNG）。所生成的数字在序列内几乎是随机的，但也取决于算法的质量。由于计算机上的随机数生成器通常是确定性的，因此也称为伪随机数生成器。 +否则，没人能够为真正的随机数生成器编写程序。 +随机数生成器的随机性以其初始状态作为展现形式，这就是众所周知的种子。选择非常数的种子值是获得良好随机数的关键，而使用常定的种子值（每次运行应用程序时使用相同的值）会导致不断出现的随机数是一种循环序列，即对于相同的初始条件，RNG会生成相同的数字序列。 + +###### 好的，差的和糟糕的-以下为相反的顺序 +最糟糕的事情是使用一个常定值作为RNG的种子。在列表407中，我们使用零点种子值， +每次我们运行应用程序时，该值都一样。因此，当我们发现打印到终端的随机数始终相同时，就不需要感到奇怪。即在每次运行包含列表407中代码的应用程序时，会在20至100之间打印相同的包含20个数字的序列。 + +``` +// random stuff +#include "Random.H" +Random ranGen (0); +for (int j = 0; j < 20; j++) +{ +Info << ranGen.integer(1, 100) << endl; +} +``` +列表407：一个简单的随机数测试，糟糕的情况 + +为了获得不同的序列，需要选用一个更好的种子值。实际上，需要选择一个每次运行程序时提供不同值的种子。对于这样的种子值，时间将是一个绝佳的选择。然而，在学习新东西的过程往往会犯错误。在源代码中看到osRandomInteger（）函数，这是一个使用随机数作为种子的随机数生成器，听起来不错。但仔细一想，这更像是个鸡和鸡蛋的问题，继续前进。 +因此我们植入了列表408的代码，这仅是一个拥有不同值的种子。然而，当运行该代码时，发现所得到的数字序列还是一样的循环序列，与前面的种子一样。 +深入研究该代码，发现其中的osRandomInteger()使用了POSIX提供的随机数生成器，但似乎没有为POSIX随机数生成器提供合适的种子。 + +``` +// random stuff +#include "Random.H" +Random ranGen(osRandomInteger ()); +for (int j = 0; j < 20; j++) +{ +Info << ranGen.integer(1, 100) << endl; +} +``` +列表408：一个简单的随机数测试，差的情况 + +按照之前的提法，时间将是绝佳的种子值选择。由于要求得到一种好的情况，需要一些新的东西而不是时间。在列表409中，使用了PID作为RNG的种子生成程序。即使重复运行多次，PID也不太可能产生相同的值。实质上，OS的内核在一定的整数范围内依次分配PIDs，例如在作者的Linux机器上某个处理器的PID值范围是1到32768。如果到达数值范围的末尾，内核就会重新气启动，并且自动跳过仍在使用的数字。此外，当程序是并行运行时，PID依然可以保证不同，即所有子进程都具有唯一的PID值。 + +``` +// random stuff +#include "Random.H" +Random ranGen(pid()); +for (int j = 0; j < 20; j++) +{ +Info << ranGen.integer(1, 100) << endl; +} +``` +列表409：一个简单的随机数测试，好的情况 + +###### 更好的情况 +如前所述，借助时间可以在每次运行程序都提供不同的种子值。函数getTime()会返回自1970年1月1日以来经过的秒数。每次运行应用程序时，列表410的代码都会产生不同的数字序列。并且，PID的重复使用也不再是问题，因为每次重新使用PID时，时间秒数肯定是不同的。当使用时间来播种RNG时，2038年的问题对我们来说不再是问题，因为只需对时间的唯一值感兴趣，而不是时间的正确表示值。 + +``` +// random stuff +#include "Random.H" +#include "clock.H" +Random ranGen(clock:: getTime ()); +for (int j = 0; j < 20; j++) +{ +Info << ranGen.integer(1, 100) << endl; +} +``` +列表410：一个简单的随机数测试，更好的情况 + +###### 绝佳的情况 +以上的方案已经接近完美，剩下唯一的问题是只针对并行运行情况。当为只需串行运行的程序实现随机数生成时，这似乎不会有难度。确实，技巧也很简单。 +可以使用当前时间作为种子值并添加PID。这样可以确保当并行运行、多个进程同时启动时，由于PID的贡献，每个进程都有其独特的种子值。 + +``` +// random stuff +#include "Random.H" +#include "clock.H" +Random ranGen(clock:: getTime ()+pid()); +for (int j = 0; j < 20; j++) +{ +Info << ranGen.integer(1, 100) << endl; +} +``` +列表411：一个简单的随机数测试，绝佳的情况 + + +#### 57.5 OpenFOAM便捷编程使用的特定宏 +##### 57.5.1 For循环 +###### …for lists +UList类定义了一些用于遍历列表内容的宏。我们经常遇到使用forAll，forAllIter或forAllConstIter等语句的循环体。其实这些语句不是C或C ++编程语言的一部分，而是Ulist类提供的宏。 +在列表412中可以看到forAll宏的定义。 + +``` +#define forAll(list , i) \ +for (Foam::label i=0; i<(list).size(); i++) +``` +列表412：UList.H中关于forAll宏的定义 + +当forAll宏被使用时，一个for循环会出现，如列表413所示。 + +``` +forAll(someList , i) +{ +// list body +} +``` +列表413：使用forAll宏来遍历列表 + +列表413中的代码在编译时会被展开成列表414所示的代码。 + +``` +for (Foam::label i=0; i<( someList).size(); i++) +{ +// list body +} +``` +列表414：将forAll 展开，这是C++预处理程序对列表413所做的处理 + +###### ... for containers +对于具有相关迭代器数据类型的容器数据类型，相应有forAllIter和forAllConstIter这两个宏。容器是一种包含顺序容器（列表）和关联容器（哈希表）的抽象概念。迭代器是用于将容器元素访问与容器内部组织进行分离的概念。 +以下，在列表415中可以看到forAllIter宏的定义，该宏使用了迭代器遍历容器，请注意函数begin()和end()的用法，以及访问单个元素的方式。列表416中的示例展示该方式。 + +``` +#define forAllIter(Container ,container ,iter) +for +( +Container :: iterator iter = (container).begin(); +iter != (container).end(); +++iter +) +``` +列表415：UList.H文件中关于AllIter宏的定义 + +在如下的列表416中，该实例遍历了拉格朗日粒子云，并且每个粒子均可以通过迭代器访问。这段代码，或者说是这段代码的开发人员，完全忽略了Lagrangian云类所使用的基础数据结构。 +借助迭代器遍历云并访问其元素，OpenFOAM开发人员可能会改变云的内部数据存储方式，但不会破坏高级代码。 + +``` +forAllIter(Cloud , c, iter) +{ +solidParticle& p = iter(); +// ... +} +``` +列表416：使用AllIter宏获取拉格朗日粒子云的所有粒子 + +###### ... for everything else +在撰写本文时，当前版本是OpenFOAM-5，此类for循环宏仅为类UList定义。因此，该宏只能用于从UList类派生的数据类型。 +如果想对不是列表或UList派生的容器数据类型等进行遍历操作，则需要以传统方式编写for循环体语句。 + +#### 57.6 时间管理 +##### 57.6.1 时间步进 +瞬态求解器需要在每个时间步至少求解一次控制方程。根据求解算法的不同，每一次外部迭代中所进行的内部迭代次数（时间步内的迭代）也不同。 + +###### pimpleFoam +列表417展示了pimpleFoam主循环的开始部分。在三个包含include指示器代码之后，一个runTime对象被创建。这表示当前时间增加，转到下一个时间步。 + +``` +/* code removed for the sake of brevity */ +Info << "\nStarting time loop\n" << endl; +while (runTime.run()) +{ +#include "readTimeControls.H" +#include "CourantNo.H" +#include "setDeltaT.H" +runTime ++; +Info << "Time = " << runTime.timeName () << nl << endl; +/* code continues */ +``` +列表417：pimpleFoam.C中关于pimpleFoam主循环的开始部分 + +###### pisoFoam +列表418展示了pisoFoam主循环的初始部分。 + +``` +/* code removed for the sake of brevity */ +Info << "\nStarting time loop\n" << endl; +while (runTime.loop()) +{ +Info << "Time = " << runTime.timeName () << nl << endl; +#include "readPISOControls.H" +#include "CourantNo.H" +// Pressure -velocity PISO corrector +{ +/* code continues */ +``` +列表418：pisoFoam.C中关于pisoFoam主循环的开始部分 + +可以看到没有任何runTime对象的递增。对此的解释是，它隐含在while语句的条件中。在pisoFoam中，while语句是通过调用函数runTime.loop()的返回值来控制的。而在pimpleFoam中，while语句由函数runTime.run()的返回值进行控制。 +进一步查看runTime.loop()。列表419显示：函数loop（）调用了函数run（），然后通过调用operator++()来递增runTime对象。 + +###### Time类的++操作符 +列表420展示了时间类中++操作符定义代码的第一行。列表420的最后一个指示器将时间值设置为当前时间加上时间步长。 + +``` +bool Foam::Time::loop() +{ +bool running = run(); +if (running) +{ +operator ++(); +} +return running; +``` +列表419：Time.C中关于函数loop()的定义 + +``` +Foam::Time& Foam::Time:: operator ++() +{ +deltaT0_ = deltaTSave_; +deltaTSave_ = deltaT_; +// Save old time name +const word oldTimeName = dimensionedScalar ::name(); +setTime(value () + deltaT_ , timeIndex_ + 1); +/* code removed for the sake of brevity */ +``` +列表420：Time.C中关于函数操作符++的定义 + +##### 57.6.2 设置新的时间步 +在瞬态模拟中，可以通过固定时间步或可变时间步来运行。在固定时间步的模拟中，时间步长是恒定的，并且必须在开始运行之前设置时间步长的值。时间步长会影响模拟的准确性和稳定性，时间步长的值还决定了计算中所能解析的时间尺度大小。时间步长还与时间积分方法的稳定性相关，这是通过Courant-Friedrichs-Lewy（CFL）准则进行关联的。 +OpenFOAM中的大多数瞬态求解器都提供了通过可变时间步长进行瞬态模拟的选项。用户提供一些关于时间步长范围限制的参数，最常用的限制是最大时间步长maxDeltaT， +这是每个新求解时间步长值的上限。同时也是用户控制所求解时间尺度大小的参数。 +决定时间步长的第二个限制参数是最大库郎数，此参数是用于保证数值解的稳定性。 +列表421展示了读取时间控制的相关代码。第一条指令读取controlDict中的内容，指定是否采用可变时间步长，这条代码的含义是很清晰的。如果controlDict中没有相应内容，则默认使用固定时间步长。其他两个语句读取最大Courant数和最大时间步长的值。最大Courant数的默认值为1.0，这是出于显式欧拉时间积分方法的限制考虑。 + +``` +const bool adjustTimeStep = +runTime.controlDict ().lookupOrDefault("adjustTimeStep", false); +scalar maxCo = +runTime.controlDict ().lookupOrDefault ("maxCo", 1.0); +scalar maxDeltaT = +runTime.controlDict ().lookupOrDefault ("maxDeltaT", GREAT); +``` +列表421：readTimeControls.H文件中的内容 + +###### 确定新的时间步长 +新的时间步长必须同时满足上述限制，即最大时间步长和最大库郎数。为了防止数值振荡，需要缓慢增加时间步长。列表422显示了如何计算每个时间步的步长值。 + +``` +if (adjustTimeStep) +{ +scalar maxDeltaTFact = maxCo/( CoNum + SMALL); +scalar deltaTFact = min(min(maxDeltaTFact , 1.0 + 0.1* maxDeltaTFact), 1.2); +runTime.setDeltaT ( +min +( +deltaTFact*runTime.deltaTValue (), +maxDeltaT +) +); +Info << "deltaT = " << runTime.deltaTValue () << endl; +} +``` +列表422：setDeltaT.H文件中的内容 + +进一步查看该代码的实际作用。 +标量maxDeltaTFact（列表422中的第3行，等式（163））表示的是最大库郎数与当前库郎数之间的关系（有关如何确定Courant数的建议，请查看第57.6.4节）。常量SMALL的作用是防止maxCo被零除，这将导致求解器崩溃。 +标量deltaTFact由maxDeltaTFact计算得出。该代码（第4行，公式（164））实现了数值缓和，即时间步的增加速率受到限制。两个min（）函数的嵌套使用可以得出三个值中的最小值。这三个数值中最明显的是最后一个参数。如果此值最小，则下一个时间步长相比上一个时间步的增长是20％。 +等式（164）以数学方式显示了前两个参数中的最小值。图131显示了等式（164）的三个参数。在标量maxDeltaTFact中使用了符号x。在图131中，x的值大于1。等式（166）阐述了相应的原因。x是最大库郎数Comax与当前库郎数Co的比值。由于当前库郎数始终小于最大库郎数，因此我们用f*Comax来替换Co，其中f<1。消去*Comax后，f以倒数的形式存在。因此，得到x始终大于1。 + +函数setDeltaT（）的参数遵循第一个限制条件，即最大时间步长。新计算得到的最小值和最大时间步长将继续传递. + +##### 57.6.3 关于时间传递的注解 +在本节中，我们将仔细研究Time类的实现。 + +###### 类的设计 +快速浏览文件Time.H，发现一些有关时间特性的有趣信息，或更准确地说，是有关Time类的特性。列表423展示了Time类是如何从五个基类中继承而来的。 + +``` +class Time +: +public clock , +public cpuTime , +public TimePaths , +public objectRegistry , +public TimeState +{ +/* class definition */ +} +``` +列表423：Time类的继承信息，Time.H.文件的部分部分内容。 + +###### TimeState类 +从列表423可以看到，由于继承的特性，Time类也是一个TimeState。从列表424可以看到有关TimeState类的继承信息，实质上，TimeState是一个dimensionedScalar。 + +``` +class TimeState +: +public dimensionedScalar +{ +/* class definition */ +} +``` +列表424：TimeState类的继承信息，TimeState.H文件的部分部分内容。 + +###### 区分时间步长 +Time是一种TimeState，而TImeState又是一种dimensionedScalar，了解这一事实有助于理解列表425中的第6、13和21行。其中调用了dimensionedScalar命名空间中的name()函数。 + +``` +Foam::Time& Foam::Time:: operator ++() +{ +// some code removed for brevity +// Save old time name +const word oldTimeName = dimensionedScalar ::name(); +setTime(value () + deltaT_ , timeIndex_ + 1); +// some code removed for brevity +// Check that new time representation differs from old one +if (dimensionedScalar ::name() == oldTimeName) +{ +int oldPrecision = precision_; +do +{ +precision_ ++; +setTime(value (), timeIndex ()); +} +while (precision_ < 100 && dimensionedScalar ::name() == oldTimeName); WarningIn("Time:: operator ++()") +<< "Increased the timePrecision from " << oldPrecision +<< " to " << precision_ +<< " to distinguish between timeNames at time " << value() +<< endl; +if (precision_ == 100 && precision_ != oldPrecision) +{ +// Reached limit. +WarningIn("Time:: operator ++()") +<< "Current time name " << dimensionedScalar ::name() +<< " is the old as the previous one " << oldTimeName +<< endl +<< " This might result in overwriting old results." +<< endl; } +} +// some code removed for brevity +} +``` +列表425：Time类中的递增运算符（++），Time.C文件的部分部分内容。 + +从列表425的第23至27行，可以看到之前章节11.3.2的列表41出现的警告消息的相应产生代码。 +在第21和29行中，可以找到关于时间精度的硬编码限制。如果时间精度值达到100，它将不再增加。 + +###### 精确命名时间 +在11.3.2节中，可以看到timePrecision值可能是误差的来源。现在，将详细说明此误差的实际原因。 +列表426展示了函数timeName（const scalar）是如何定义的。该函数用于从指定标量格式中创建相应正确格式的时间名称。在这个函数中，时间精度以数据成员precision_的形式发挥作用，该成员是Time类的静态数据场量，具有受保护的可见性。在这个函数中，高精度的时间值将被转换为有限精度的字符串形式（时间名称）。 + +``` +//- Return time name of given scalar time +Foam::word Foam::Time:: timeName(const scalar t) +{ +std:: ostringstream buf; +buf.setf(ios_base :: fmtflags(format_), ios_base :: floatfield); +buf.precision(precision_); +buf << t; +return buf.str(); +} +``` +列表426：Time类中的函数timeName（const scalar）；来自Time.C的部分内容。（请注意，描述性注释来自头文件Time.H。） + +当时间前进时，例如使用Time类的增量运算符，函数setTime（）会被调用。列表427显示了此函数的定义。新的时间值将传递给此函数。在第二条指令中，可以看到时间名称是如何被更新为新值的。 + +``` +void Foam::Time:: setTime(const scalar newTime , const label newIndex) +{ +value() = newTime; +dimensionedScalar ::name() = timeName(timeToUserTime(newTime)); +timeIndex_ = newIndex; +} +``` +列表427：Time类中的函数setTime( )，来自Time.C文件的部分内容。 + +setTime（）函数被类似Time类的运算符*调用，参见列表425的第8行。时间索引会增加1。从TimeState类的头文件可以看到，时间索引是属于数据类型label，它实际上是整数数据类型。因此，时间索引是系列连续的数字，可以用来计算时间步长。 + +##### 57.6.4 库郎数 +库仑数Co是时间步长∆t与特征对流时间尺度u / ∆x的比值。等式（168）显示了库郎数是如何定义的。但在实际CFD代码中，库郎数的计算方式会稍有不同。等式（169）说明等式（168）可以通过A / A变形展开得到一个新的等式，该等式的特征参数是通量和控制体的体积而不是速度和离散长度。等式（170）展示了等式（169）的一维有限体积扩展形式。通过E和W面的平均通量值定义对流时间尺度。在一维情况下，此定义显而易见以某种方式体现。但对于二维或三维情况，如何定义特征通量并非是能直截得出的。 + +（168） +（169） +（170） + +###### OpenFOAM中的库郎数 +在OpenFOAM中，所有单元格的库郎数都被计算。但实际上，OpenFOAM最后得出的是最大库朗数（即所有单元格中最大的库郎数），以及平均库郎数（即所有单元格中的平均库郎数）。 +列表428展示了负责计算库郎数的代码。列表428的第8行对应的是等式（171）。sumPhi是一个标量场，包含了每个网格单元的所有面通量大小之和，即将每个网格单元的面通量大小相加。等式（171）对每个单元格都成立。 + +等式（172）是代码中第11行的数学表示。主要是确定sumPhi值与单元体积比值的最大值。其中变量sumPhi和mesh.V（）均包含所有网格单元的值。因此，gMax（）函数将返回最大值。 +等式（173）表示代码第14行。 + +``` +scalar CoNum = 0.0; +scalar meanCoNum = 0.0; +if (mesh.nInternalFaces ()) +{ +scalarField sumPhi +( +fvc:: surfaceSum(mag(phi))().internalField () +); +CoNum = 0.5* gMax(sumPhi/mesh.V().field ())*runTime.deltaTValue (); +meanCoNum = +0.5*( gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue (); +} +Info << "Courant Number mean: " << meanCoNum +<< " max: " << CoNum << endl; +``` +列表428：文件CourantNo.H的内容 + +等式（171） +等式（172） +等式（173） + +###### 讨论 +如上所述，在三维情况下计算Courant数的方法并不简单。本部分只是说明作者的理解， +因此不能保证有效性。比值1/2和φfi的求和将解释如下。 +以二维控制体积为基础，等式（175）展示了求和的完全形式。然后重新排列该方程式得到等式（176）。在等式（176）中，求和被减少为两项，分别是在N-S和W-E主方向上的面通量算术平均值。之后将该总和作为主要方向上平均面部通量的L1范数。 +选择L1范数的原因很简单，相比Euklidian或L2范数，L1范数在计算上更容易。并且使用L1范数看起来是合理的，因为它可以测量运动所覆盖的距离，具体请参见：http://en.wikipedia.org/wiki/Taxicab_geometry。 + +等式（174） +等式（175） +等式（176） +等式（177） +等式（178） +等式（179） +等式（178） +继续引入以下符号 + +等式（179） +等式（180） + +粗略一看，平均库郎数的计算方法似乎不正确，但事实并非如此。 + +等式（180） + +其中x的平均值定义为： + +等式（181） + +接下来，写出库郎数的平均值。其中未标记的求和是针对所有单元进行的。 + +等式（182） +等式（183） +等式（184） + +现在等式（184）类似于（173），两者之间唯一的差别是术语X，我们将进一步讨论它。 + +等式（185） +等式（186） +等式（187） +等式（188） + +假设Vcell/Vcell≈1 +等式（189） +等式（190） + +因此，平均库郎数meanCoNum的计算方法实际上是对库郎数求均值Co。但是，以上是基于一些假设进行的测试和验证。 +首先，根据作者的解释，面通量之和的意义取决于六面体单元，但该论点似乎不适用于四面体单元。其次，Vcell≈1的假设只对同类网格有效。对于均一网格，此假设可以很好地得到满足。但如果最大和最小单元的体积相差很大，则此假设不成立。 + +###### 关于计算成本的一些思考 +为什么将平均库郎数的计算公式从等式（191）转换成等式（192），作者其实也不确定。 +作者认为这样做是出于计算成本的原因。对所有场的值进行求和以及一个相除的计算，相比对两个场量进行按元素相除以及后续所有的相应场元素进行求和，在数学上要节省一半的时间。 +如果除法操作比加法操作花费更多的时间，那么将能实现节省时间的操作，并且该情况是很有可能的。根据系统的不同，浮点除法运算所花的时间可能比浮点乘法要长几倍。 +在第一种情况下，需要进行n次除法和一次加法，其中场量的数量为n。在第二种情况下，相加2n次，然后只进行一次除法。 + +等式（193） + +引入因子，这是Td和Ts的比值。 + +等式（194） +等式（195） +等式（196） +接下来假设 n是足够大的 + +等式（197） +因此，第一个公式执行1 +δ次运算，而第二个公式执行约2n运算。如果δ大于1，则第二个公式将花费较少的时间进行计算。 δ小于1的可能性很小，甚至不可能，因为加法是非常简单的操作。请记住，δ是除法时间与加法时间之间的比率。实际比率的值取决于系统架构、编译器和编程实现形式而有所不同，例如[1]报告单精度和双精度浮点除法的系数为5到6。该参数不考虑所涉及操作的内存使用情况，仅关注浮点操作的数量。 +由于Courant数是在每个时间步迭代之后进行计算的，所以计算Courant数所需的时间会影响数值模拟所需的时间。 + +##### 57.6.5 两相库郎数 +在两相流模拟中，有几种可选方式来计算Courant数。总共有4个速度场（U1，U2，U和Ur）。这些是相1和相2的速度、混合物的速度以及相对速度。求解器twoPhaseEulerFoam将计算混合物和相对速度的库郎数。 +列表429显示了CourantNos.H的内容，该文件是此求解器源代码的一部分。第一行通过预文件CourantNo.H计算混合物的库郎数。这是第57.6.4节中已经进行描述过的文件。当此代码在场phi上运行时，它恰好是混合物的通量，因此将计算得到混合物的库仑数。 +下一行根据相对相通量计算库郎数。在第11行中，这两个库郎数中的最大值将被确定并存储到变量CoNum中。 +CoNum是时间步进机制使用的库郎数。因此，TwoPhaseEulerFoam求解器的可变时间步长是基于混合物和相对速度的这两者中的最大库郎数进行确定的。 + +``` +#include "CourantNo.H" +{ +scalar UrCoNum = 0.5* gMax +( +fvc:: surfaceSum(mag(phi1 - phi2))().internalField ()/mesh.V().field() )*runTime.deltaTValue (); Info << "Max Ur Courant Number = " << UrCoNum << endl; +CoNum = max(CoNum , UrCoNum); } +``` +列表429：文件CourantNos.H的内容 + +#### 57.7 注册机制 +通过浏览用户文档或互联网学习OpenFOAM源代码的某个时刻，，我们可能会遇到诸如注册对象的名词或类似的表达式。本节试图对此主题进行一些阐述，或者至少介绍作者的想法和发现。本节与57.8节密切相关。 + +##### 57.7.1 涉及的类 +这是在各个类的头文件中能找到的描述语句。 +###### IOobject +IOobject 定义了支持objectRegistry隐式管理机制的对象属性，并提供用于执行I/O流的基本结构。 +###### regIOobject +regIOobject是从IOobject类继承过来的抽象类，它通过regIOobject处理自动对象注册。 +###### objectRegistry +objectRegistry为regIOobject的注册表 + +图132展示了围绕类regIOobject的类继承层次结构详细信息。 + +图132 关于regIOobject类层次结构的局部展示； +注意，该图仅对IOobject和regIOobject类是完整的——意味着IOobject没有从任何其他类继承，而regIOobject仅是从IOobject继承过来的；但其他类拥有比上图更多的基类 + +###### IOobject +该类提供了I/O的一些基本特性。在章节57.4.5中展示了该类的实际用法或典型用法。 +###### regIOobject +该类是一个抽象类，正如头文件中所描述的。在图132中，纯虚函数的名字以斜体字形式显示，正是纯虚函数的存在使得该类变成抽象类。这意味着所有从regIOobject继承而来的类都必须重载该纯虚函数。同时，这也表示无法直接在regIOobject类中创建对象。因此，在所有OpenFOAM的源代码中，只能从regIOobject派生类的初始化列表中找到调用regIOobject类的构造函数。 + +###### objectRegistry +objectRegistry的命名与本节相同。实际上，OpenFOAM中没有特定的注册表，而是几种不同的。其中，Time类、cloud类和polyMesh类都是从objectRegistry派生的。图133显示了从objectRegistry继承的类。 + +图133：类objectRegistry的基类；该类从regIOobject和HashTable派生而来； +注意，HashTable的模板参数是regIOobject的指针。因此，objectRegistry是一个regIOobject，也是一个由regIOobject类型指针组成的的HashTable——这正是C ++的模板和继承机制的疯狂和奇异之处。 + +###### objectRegistries：空间和时间 +在OpenFOAM中，通常只有一个Time类型的对象，通常称为runTime。在作者的知识范畴内，没有求解器会同时使用多个Time类的实例。由于大多数求解器也仅具有一套网格，因此作为注册表的Time和polyMesh之间的分隔似乎太大了。 +但是，有些求解器具有多个网格，例如共轭传热求解器。其中最简单的配置是，求解域的固体区域有一套网格，流体区域有另一套网格。但是，求解器chtMultiRegionFoam支持任意数量的流体域和固体域。在这种情况下，固体区域i的温度场T需要选择合适的注册机制，即与固体区域i的网格进行结合。由于流体区域中的流体也具有温度场，因此也需要将这个温度场注册到相应的流体网格中。在有两个温度场T对象的情况下，注册表（即网格和时间）的分隔使我们避免了潜在的名称冲突。因此，网格需要在runTime对象注册表中进行注册，而场是在相应的网格中进行注册。 +OpenFOAM求解器有大量正在使用的对象注册表，最主要的是runTime和mesh对象。对于场量来说，重要的是要知道它们属于哪套网格，这是因为field场量仅仅是一系列数值的列表，只有与网格相联接才赋予field实际的含义，即列表中位置i处的场量是单元格i的中心值。此外，该场量还需要与模拟的实际时间状态关联，否则将缺乏有实际含义的方式来定义或计算时间导数。 + +##### 57.7.2 使用注册表 +一个对象注册表能有什么用？好吧，查看代码。 +在章节66.3中展示了一种搜索特定类型文件的方式。现在，在所有文件中通过文件拓展名.C来搜索lookupObject类型的文件，并计算匹配的总数。列表430显示了可以使用的命令。首先，使用find查找具有指定文件名的所有文件，然后将结果传递到grep命令，后者会在上述文件中继续搜索指定形式的文件。最后，将grep的结果传输到wc命令，再计算相应的行数、字数和字节数。因此，此命令序列返回的第一个数字对应的是匹配的次数。实际匹配数大约是显示数字的一半，这是因为从源代码构建OpenFOAM执行文件的过程中，在lnInclude文件夹中创建了符号链接。 + +``` +find $FOAM_SRC -name ’*.C’ | xargs grep ’lookupObject ’ | wc +``` +列表430：查找和搜索扩展名为.C、并且类型为lookupObject的文件，并计算匹配的次数 + +在作者安装的OpenFOAM-2.3.x中显示了列表430代码运行的结果，总计1068次匹配数，其中有537个文件来自lnInclude目录的符号链接。这也意味着lookupObject（）被大量使用，那么lookupObject（）有什么用？ + +###### 需要知道vs.想知道 +封装或信息隐藏的思想是面向对象编程的基本原理184。其总体思路是将某些东西的实际实现形式隐藏在可公开访问的界面后面。因此，可以对类内部的工作原理进行改变，而不会影响其实际使用。迭代器的概念就是一个关于信息隐藏优越性的示例。典型的容器类实现了被称为迭代器的功能，该功能用于遍历迭代容器的所有元素。通过使用迭代器的公共接口，实际使用的容器可以是任何类型的数据结构（比如链接列表，向量，哈希表等）。 +除了能提供和使用类数据的访问接口之外，这也是一种常见且较好的数据范围限制方法，例如临时数据对于实际使用的类或函数来说是局部的。在植入类的实现设计中，可以将包含在类中和/或传递给类的数据限制为必要的最小值，比如在求解器中使用的粘度定律不需要知道使用了哪个求解来求解离散化方程系统。但是，它们可能需要访问数据，而某些类的原始设计者并没有想到这一点。 + +###### 命名空间 & 作用域 +另一方面是源代码中的名称空间和变量作用域。变量在命名空间和它声明的范围内是可见的。如果查看求解器的底层代码，例如twoPhaseEulerFoam，可以看到很多#include语句和主函数main（）。尽管没有直接找到涉及命名空间的声明语句，但在文件fvCFD.H中有一条声明被隐藏，使得编译器可以使用名称空间Foam。这也是为什么可以在之后例如createFields.H的文件中使用在名称空间Foam中所定义类型名的原因，例如volScalarField。否则，也需要明确指定名称空间，例如Foam :: volScalarField。因此，所有由求解器创建的对象（例如网格，运行时等）在命名空间Foam内都是可见的。 +然而模型也有自己的名称空间。列表431显示了这样具有自己命名空间的一个模型。在命名空间Foam中，一个新的命名空间DiameterModels将被创建。在此命名空间中，定义了一个isothermal类。因此直径模型类的实现不会污染Foam名称空间。 +尽管diameterModels命名空间是Foam名称空间的子集，并且在Foam中声明的所有内容均在Foam :: diameterModels中可见，但是该直径模型也必须与其他模型一起编译到共享库中。因此，在编译这些文件时，编译器完全不知道由例如createFields.H文件所创建的属于Foam命名空间的对象。 + +``` +namespace Foam +{ namespace diameterModels +{ +class isothermal +: +public diameterModel +{ +// code removed +} +``` +列表431： isothermal类的类定义文件，继承自isothermalDiameter.H中的类diameterModel + +###### 查找内容 +列表432显示了isothermal.类的函数d（）是如何定义的。出于上述isothermal.C和reateFields.H位于不同编译单元中的原因，即使压力场p是命名空间Foam的一部分，也无法直接从方法内部内访问压力场p。然而，在其他直径模型不需要访问压力场，例如实现恒定直径的constant. + +``` +Foam::tmp Foam:: diameterModels :: isothermal ::d() const +{ +const volScalarField& p = phase_.U().db().lookupObject +( +"p" +); +return d0_*pow(p0_/p, 1.0/3.0); +} +列表432：isothermalDiameter.C中diameterModel类的函数d（）的定义 +``` +上面的示例显示了lookup机制的价值。由于某些子模型在某些场上运行，很容易从该场中获取对网格的引用，正如在phase_.U（）.db（）是如何进行的。Phase_是直径模型所对应基类的成员。调用phase_.U（）将返回一个对相关相位速度场的引用。所幸速度场是在网格中注册的，否则我们将不知道哪个速度值属于之前应用的某个特定单元格，并且该引用是通过调用db（）来获得对实现的，这是IOobject类的函数方法。这种方便的机制避免对相关子模型的污染，由于这些子模型引用了网格、时间以及某些时候可能需要的场量，或者在某些特殊情况下可能需要的派生类。 +因此，lookupObject（）方法提供了一个工具，以获取在编译时可能未被声明场的引用，并因此该场变为可用。请记住，压力场是在求解器的createFields.H中进行声明的，该文件与直径模型所需编译的库处于不同的编译单元中。如果直径模型和求解器的代码位于同一编译单元（求解器的可执行文件）中，则不需要查找机制。但是，由于OpenFOAM的开发人员追求模块化，因此将所有内容放到一个编译单元中违反了模块化和可重用性的设计原则。 +lookupObject（）方法是模板化的，由于可以在网格物体上注册任何东西，实际上是任何从regIOobject派生出来的东西，并且objectRegistry是regIOobject指针类型的HashTable。因此，在编译时，函数和编译器并不确切知道它将要处理的数据类型。这是模板起作用的地方，模板化方法只被template参数执行一次，当使用该方法时，只需将模板参数替换为实际类型，如lookupObject （“ p”）一样。之后编译器完成其余工作并生成合适的代码。也可以通过使用函数重载来解决此问题，这样无需使用模板，但使用函数重载会以大量代码重复和可维护性差为代价。 + +##### 57.7.3打印注册表 +如果您感到好奇，可以将以下代码行添加到您自己的实用工具进行测试，主要检查在网格和RunTime对象注册表中具体注册了什么。注意，必须从放置代码的位置访问mesh和runTime。此外，在某些情况下，对象的名称可能会有所不同。 +``` +Info << "mesh.names() " << mesh.names() << nl << endl; +Info << "runTime.names() " << runTime.names() << endl; +``` +列表433：将mesh和runTime注册对象的内容打印到终端 + +#### 57.8 I/O-输入&输出 +关于I/O的一些内容已经在57.4.5和57.3章节中介绍过。但是，这方面的材料在组织设计上比较零散或者没有完全设计，这里用一种更普遍的方式进行介绍。 + +##### 57.8.1 输出到终端——OpenFOAM特有的pritf() +在编程中通常需要将一些材料打印到终端，也就是pritf() 调试。在C++中，一般的I/O是基于I/O流实现的，C++的I/O流为内置和用户自定义数据类型都提供了一种类型安全且统一的方式来实现。有关C语言pritf()函数和C++流使用方面的内容，参见列表434和435。 + +``` +#include +int main(int argc , char** argv) +{ +printf("Hello , World!\n"); +return 0; +} +``` +列表434：C的Hello World案例 + +``` +#include +int main() { +std::cout << "Hello World!" << std::endl; +return 0; +} +``` +列表435：C++的Hello World案例 + +OpenFOAM 植入了自己的流库，其通用流库是基于Iostream类，该类的头文件描述中说明了这样做的原因： + IOstream是所有输入/输出系统的抽象基类，包括流、文件以及令牌列表等。 +其中的基本操作包括构造，关闭，读取标记，读取基元和读取二进制块。此外，在附加版本中合并了控制和行计数。通常，人们会使用读取初始成员函数，但是如果正在读取的是未知数据序列上的流，则可以逐个读取，然后进行分析。 + +OpenFOAM可以处理各种流方面的通信操作，其中包括：与用户交互终端I/O，文件I/O和用于并行处理的进程通信。列表436显示了OpenFOAM范畴的Hello World示例，看起来与C ++的示例非常相似。 + +``` +#include "Istream.H" +using namespace Foam; +int main(int argc , char *argv []) +{ +Info << "Hello OpenFOAM!" << endl; +return 0; +} +``` +列表436：OpenFOAM范畴的Hello World案例 + +###### 条件（调试）输出 +printf（）调试是一种非常方便的底层技术，可以用于对代码段进行故障排除。在实际调试时，一旦成功debug将会删除所有打印到终端的代码行。然而，人们也许想要开发软件，该软件可以是健谈的也可以是静默的187。在这种情况下，我们需要条件Info语句。 +列表437显示了Hello World！的带条件输出示例。该列表很长，这是因为决定不采用简单的布尔值来控制条件输出。取而代之的是一种实际情况，其中的详细程度由命令行选项控制。但是需要更多的代码行来处理命令行参数。 + +``` +#include "argList.H" +bool verbose(false); +using namespace Foam; +int main(int argc , char *argv []) +{ +argList :: addNote +( +"This is a \"Hello World !\" program for the OpenFOAM world." +); +argList :: noBanner (); +argList :: noParallel (); +argList :: removeOption("noFunctionObjects"); +argList :: removeOption("case"); + +argList :: addBoolOption ( +"verbose", +"control the chatty -ness of me" +); +Foam:: argList args(argc , argv); +if (args.optionFound("verbose")) +{ +verbose = true; +} +Info << "Hello OpenFOAM!" << endl; +if (verbose) Info << "... and hello to all other non -OpenFOAM worlds!" << endl; +return 0; +} +``` +列表437：用OpenFOAM编写带条件输出的Hello World！示例 + +除了布尔命令行开关之外，还添加了一条注释来告知用户有关可执行文件的信息。当调用可执行文件的命令行选项-help来显示使用情况消息时，将显示此注释。 OpenFOAM默认情况下会添加许多命令行参数，这里删除了其中的一些参数（对于Hello World！程序没有意义的参数，例如parallel选项）。 +代码的倒数第二行是实际控制条件输出的代码，采用了古老但好用的if语句。 +在OpenFOAM-2.3.x的函数对象源代码中，可以观察到了另一种定义有条件输出的方式。在那里，可以将一个参数传递给Info。当使用OpenFOAM-2.4.x和更高版本时，该文件不再编译。 + +``` +// OpenFOAM -2.3.x +Info(log_)<< "Including porosity effects" << endl; + +// OpenFOAM -2.4.x and higher +if (log_) Info << "Including porosity effects" << endl; +``` +列表438：在不同版本OpenFOAM实现由Switch log_控制的条件输出。 +本示例取自力函数对象，参见文件force.C。 + +2016年5月，一种使用Info（log_）条件语句的变体被恢复，其定义如下： + +``` +//- Report write to Foam::Info if the local log switch is true +#define Log +if (log) Info +``` +列表439：Log宏的定义，最初用于功能对象的条件输出。参见文件messageStream.H + +##### 57.8.2注册机制、I / O或runTime.write（）背后的真相 +当您想将当前的模拟状态写入磁盘时，通过runTime对象注册机制进行场注册会更轻松。在很多求解器甚至可以说在所有求解器中，都在main函数的主循环内找到了类似runTime.write（）的指令。该指令通过调用函数write（）实现将场写入磁盘。当每个求解器都向磁盘写入不同的系列场数据时，可能需要问自己，当调用runTime对象的write（）函数时，求解器或OpenFOAM如何知道要写入的时哪些场？这里，其实是Time类的注册机制特性开始发挥作用。由于最终需要读取或写入的所有场都是被注册到runTime对象，runTime会拥有一系列对象（实际上是regIOobjects）需要（或可能要）被写入。实际上，由于objectRegistry是从HashTable类型派生而来的，对象注册表其实是一系列最终要写入（或可能写入）的对象。这里对Time类中write函数的调用会导致Time自身执行迭代（runTime其实是继承自 regIOobject的列表），并且调用列表中每个单独项目的write（）函数。而函数write（）是在regIOobject类中完成定义的。 + +如果进一步考虑C ++的某些规则，对源代码的研究就更深入。列表440展示了当runTime时调用write（）是如何实现的数，请记住Time是由regIOobject类的二代衍生类，其中的中间派生类是objectRegistry。该列表显示了对函数writeObject（）的调用。 + +``` +bool Foam:: regIOobject ::write() const +{ +return writeObject +( +time().writeFormat (), +IOstream :: currentVersion , +time().writeCompression () +); +} +``` +列表440：regIOobjectWrite.C中类regIOobject的函数write（） + +``` +bool Foam::Time:: writeObject +( +IOstream :: streamFormat fmt , +IOstream :: versionNumber ver , +IOstream :: compressionType cmp +) const +{ +if (outputTime ()) +{ +// some code removed +timeDict.regIOobject :: writeObject(fmt , ver , cmp); +bool writeOK = objectRegistry :: writeObject(fmt , ver , cmp); +// further code removed +``` +列表441：TimeIO.C中Time类的函数writeObject（）的部分内容 + +如果搜索Time类及其所有基类的源代码，则会发现Time，regIOobject和objectRegistry +都定义了一个称为writeObject（）的函数。其实这三种函数都共享了相同的声明，即它们接收的是相同的函数变量。由于对writeObject（）的调用没有进一步指定具体的名称空间，实质上它是Time类的writeObject（）函数，并且runTime的类型为Time，因此当调用runTime.write（）时其实调用的是该函数。 + +在列表441中，可以看到Time类的函数writeObject（）定义的部分内容。在那里可以显式调用regIOobject和objectRegistry类的函数writeObject（）。 + +因此，调当用runTime.write（）时，所有三个类（Time，regIOobject和objectRegistry）的writeObject（）函数都将被调用。值得注意的是对regIOobject :: writeObject（）的调用其实是在timeDict对象上实施的。并且，此对象的定义是在调用之前的部分移除代码。查看源代码可以发现，timeDict是一个IOdictionary，它也是一个从regIOobject派生的类，具体参见图132。对timeDict.writeObject（）进行调用的具体执行代码是在时间目录中创建统一文件夹。 + +类objectRegistry的writeObject（）函数实际上通过注册表对其中的所有元素进行迭代。列表442展示了regIOobject指针对应的哈希表是如何进行迭代的。只要写入标记未被设置为NO_WRITE，都会为每个元素调用writeObject（）。实际上调用的是类regIOobject的函数writeObject（），这是因为迭代是通过regIOobject指针进行的。在列表442第16行上的调用将导致已注册的场被写入磁盘。 + +``` +bool Foam:: objectRegistry :: writeObject +( +IOstream :: streamFormat fmt , +IOstream :: versionNumber ver , +IOstream :: compressionType cmp +) const +{ +bool ok = true; +forAllConstIter(HashTable , *this , iter) +{ +// code removed handling debug output +if (iter()->writeOpt () != NO_WRITE) +{ +ok = iter()->writeObject(fmt , ver , cmp) && ok; +} +} +return ok; +} +``` +列表442：objectRegistry.C中类objectRegistry的函数writeObject（）的部分代码 + +总而言之，通过对OpenFOAM源代码的挖掘学习到调用runTime.write（）对应的内部工作原理。首先，Time类将当前状态写入到磁盘的统一文件夹，然后runTime对象的objectRegistry部分会被写入到所有已注册场上。在第57.6节中已经提到，时间类具有倍数除法的特性，并且其中一些甚至具有继承特性。这就强调了必须对C ++有一定的了解，以便能够从OpenFOAM的源代码中推断出具体发生了什么，由于OpenFOAM使用了大量的C ++语言功能，例如多重继承，多态性和模板。在涉及编程案例的上下文中，OpenFOAM使用了（以下并不限于此）的功能：面向对象和泛型编程。 + +#### 57.9形成论点–传递论点 +在表格437和章节57.8中，通过命令行参数来控制应用程序的行为。在章节11.1中讨论了在应用程序之外实施控制的不同方法。需要指出的是，命令行参数是对应用程序的最低控制级别，在每次运行程序时都需要指定命令行参数，并且其中仅有可选参数的默认值是永久性的。本节讨论有关命令行参数的一些要点。 + +##### 57.9.1 事物的顺序 +如果对使用自定义命令行参数的应用程序进行研究，可以发现事物是有一定顺序的。 +列表443显示了如何定义命令行参数的一个简单示例。 + +首先，对静态函数addBoolOption的调用会被添加到自定义的命令行参数。然后，setRootCase.H文件会被包含到头文件。列表444显示了该文件的内容。具体可以看到一个Foam :: argList类型的变量被创建，并将所有命令行参数传递给Foam :: argList的构造函数，并且该构造函数仅对已定义选项的有效性开展检查。实际上，在对Foam :: argList的构造函数进行调用之后，无法继续添加其他选项。一旦构造函数被调用，一个名为args的对象会出现，该对象可用于查找选项并提取信息。 + +``` +argList :: addBoolOption +( +"verbose", +"be more talkative" +); +#include "setRootCase.H" +const bool verbose = args.optionFound("verbose"); +``` +列表443：定义命令行参数相应源代码的执行顺序 + +``` +Foam:: argList args(argc , argv); +if (!args.checkRootCase ()) +{ +Foam:: FatalError.exit(); +} +``` +列表444：文件setRootCase.H的内容 + +##### 57.9.2 处理SPAAAACE! +在参数名称中使用空格并不是一个好主意，因为终端通常会把空格理解为为输入参数的结尾。在常见的UNIX / Linux工具中，多字符的参数通常用连字符届分隔，例如—auto-compress。在OpenFOAM范畴内，可以看到camel案例的使用，即借助-noFunctionObjects处理多字符参数名称。 +但是，如果确实希望在参数名称中包含空格，OpenFOAM也预留了相应的方案。列表445演示了如何定义名为search point的参数。 + +``` +argList :: addOption ( +"search point", +"vector", +"find the cell containing the specified coords at - eg, ’(1 0 0)’" +) +/* no comment */ + +vector v; + if (args.optionReadIfPresent("search point", v)) +{ +Info << "Searching cell at point: " << v << endl; +} +``` +列表445：从命令行参数读取一个点的坐标 + +使用带有空格的参数需要特别注意，如列表446所示。通过引号防止终端将参数名称中的空格理解为参数名称的结尾。 +列表447显示了使用空格定义参数的危害。在这种情况下，引号的使用不再像常规的情况。 + +``` +findCellByPointCoords -’search point ’ ’( -0.993389 -1.90411 12.4942) ’ +``` +列表446：传递一个点的坐标到应用程序，注意其中引用的参数名 + +``` +user@host:∼$ findCellByPointCoords -search point ’( -0.993389 -1.90411 12.4942) ’ +... +--> FOAM FATAL ERROR: +Wrong number of arguments , expected 0 found 2 +Invalid option: -search +FOAM exiting +``` +列表447：传递一个点的坐标到应用程序，省略了参数名的引用将导致报错 + +#### 57.10 湍流模型 +在章节32.3中，已经强调了用户可以在以下选项进行选择； +1.层流模型 +2.RANS湍流模型 +３.LES湍流模型 + +该语句反映在OpenFOAM中实现湍流模型相应的类关系。面向对象编程允许程序员对从人类语言到源代码的关系进行直接转换。关于湍流模型，有两种说法： +１.所有的RAS湍流模型都是湍流模型，但不是所有的湍流模型都是RANS模型； +２.RAS模型与LES模型不同，但他们都是湍流模型。 + +两种声明都可以使用湍流模型类图来反映。最上面的是一个抽象类turbulenceModel。该抽象类提供了所有延伸类的框架。同时，所有湍流模型类的可能共有泛函特征都可以在该类进行定义。之后所有延伸类都继承自该泛函。也就是，每个湍流模型都从该抽象基类继承而来，并且每个湍流模型都会进一步植入特有的函数特征。 + +图片134：湍流模型类的继承类图 + +##### 57.10.1 turbulenceModel抽象基类 +基类turbulenceModel是抽象类，它包含几个纯虚函数。为了能够调用这种函数，这些函数必须在基类派生的子类中进行重载。不能直接调用纯虚拟类。列表448显示了纯虚函数或抽象函数的声明。 其中= 0表示函数是抽象的。 + +``` +//- Return the turbulence viscosity +virtual tmp nut() const = 0; +//- Return the effective viscosity +virtual tmp nuEff() const = 0; +``` +列表448：turbulenceModel.H中虚函数的定义 + +基类不仅包含虚函数，它还包含了所有派生函数所共有的函数类。因此，此函数是在基类声明的。列表449显示了函数nu（）的实现，该函数能用于访问层流或分子粘度。层流粘度是流体本身的特性，与是否湍流无关。因此，湍流模型也需要获取层流粘度的权限。 + +``` +//- Return the laminar viscosity +inline tmp nu() const +{ +return transportModel_.nu(); +} +``` +列表449：turbulenceModel.H中nu()的实现 + +每个从抽象基类派生的子类都至少对抽象函数进行一次重载。而基类中的非抽象函数——比如列表449中的nu( )，可以在派生类中直接使用。无论采用RAS还是LES湍流模型，层流粘度永远保持不变。 + +##### 57.10.2 RASModel类 +RASModel类由抽象类turbulenceModel派生，类RASModel本身是所有RAS湍流模型的基类。他同样是一个抽象类，这是因为他没有对turbulenceModel继承过来的虚函数进行重载。 +然而，RASModel实现了所有RAS湍流模型所共有的所有函数。列表450显示了RASModel类中nuEff（）函数的具体实现。 + +``` +//- Return the effective viscosity +virtual tmp nuEff() const +{ +return tmp +( +new volScalarField("nuEff", nut() + nu()) +); +} +``` +列表450：RASModel.H中nuEff()的实现 + +有效粘度nuEff由流体特性参数中的层流粘度和湍流粘度计算而来。其中，湍流粘度是湍流模型的属性。列表450中的函数nu（）是在turbulenceModel类中实现的，参见列表449。而函数nut（）不是由RASModel类实现的，因此此方法必须在RASModel派生的类中实现。 + +##### 57.10.3 RAS湍流模型 +所有RAS湍流模型都是从RASModel类派生而来的。每个派生类都必须实现剩下的所有抽象函数。图135显示了简化的类图– OpenFOAM中提供了许多的RAS湍流模型。 + +图135：RAS湍流模型的继承类图 + +##### 57.10.4 kEpsilon类 +kEpsilon类是从RASModel继承过来的。 + +``` +class kEpsilon : +public RASModel +{ +/* class definition */ +} +``` +列表451 kEpsilon.H中kEpsilon类的定义 + +函数nut（）必须在kEpsilon中具体实现。列表452显示了如何实现函数nut（）。该函数仅返回类成员变量nut_。 + +``` +//- Return the turbulence viscosity +virtual tmp nut() const +{ +return nut_; +} +``` +列表452 kEpsilon.H中nut()的实现 + +在所有RAS湍流模型之间，nut_的计算方式也有所不同。具体请参阅第61.2.2节中的列表486。 + +#### 57.11 调试机制 +OpenFOAM提供了便捷的调试机制。创建附加模型库时可以调用此机制。 OpenFOAM Wiki的一节介绍了内置调试机制。 +全局调试标志（控制整个系统的调试行为）在\ $ FOAM_SRC /../ etc / con中被指定。从OpenFOAM-2.2.0开始，可以在案例的controlDict中选定调试标志，并把全局调试标志进行覆盖。 +由于这种调试机制依赖于内部变量，在使用这种类型的调试机制时不会有再次编译发生。这种调试有时也被称为printf调试。 +默认情况下，所有调试开关均被初始化为零值，因此特定类的调试功能已被禁用。但当求解器设置案例时，对全局和本地条目的检查操作被执行。列表453显示了controlDict中用于覆盖调试开关的条目。列表454显示了求解器的输出信息，说明了controlDict中的本地设置。 + +``` +DebugSwitches { +DefaultStability 0； +YoonLuttrellAttachment 1； +} +``` +列表453：在案例的controlDict中指定调试开关 + +``` +Overriding DebugSwitches according to controlDict +DefaultStability 0; +YoonLuttrellAttachment 1; +``` +列表454：在案例的controlDict中指定调试开关时的求解输出信息 + +##### 57.11.1 使用调试机制 +如果为c某个类启用调试机制，列表455显示了实际将会如何调用它。该代码是非常简单的，在该场景下一个名为debug的变量被提供。我们仅在if语句中使用此变量。 + +``` +// print debug information +if (debug) +{ +// debug action +} +``` +列表455：在类中使用调试机制 + +##### 57.11.2 使用案例：写入中间场 +列表456显示了名为Ea的函数是如何定义的。出于调试目的，希望将中间场写入磁盘。在列表456的第7行中，雷诺数被计算并将存储到ReB中，ReB被用于生成函数的返回值。在正常操作中，只有返回值是被关注店。但在调试时，中间结果可能也令人感兴趣。默认情况下，场量ReB不会被写入到磁盘，并且在超出范围时该函数就不再存在，即当函数到达其范围的末尾时，变量ReB会被自动删除。 +请注意在第7行中所传递的参数。第一个是场量的名称。我们可以忽略这个论点，但当变量ReB被写入磁盘时，文件名将被第一个参数确定。如果省略此参数，则会使用基于场量生成方式的自动命名名称。在这种情况下，写入的文件将被命名为max（（（mag（（U1-U2））* d）| nu），0.001）。第7行的公式可以很清楚地识别。但包含特殊字符（非字母数字字符）的文件名通常是不可见的。 +在第14行中，write（）函数被手动调用。该函数适用于所有已注册的输入/输出对象204。当从I/O注册的对象Ur构造局部变量ReB时，可以轻松假设ReB也会是这种类型。 + +``` +Foam::tmp Foam:: YoonLuttrellAttachment ::Ea +( +const volScalarField& Ur, const dimensionedScalar& dP +) const +{ +// do stuff +volScalarField ReB("ReB", max( Ur*dB/phase2_.nu(), scalar (1.0e-3) )); +// debug instructions +if (debug) +{ +if (Ur.time().outputTime ()) +{ +ReB.write(); +} +} +// do more stuff } +``` +列表456：手动写入调试相应的中间场 + +#### 57.12 浏览运行时选择和调试魔力 +OpenFOAM提供了一些惊人的功能，例如对流体求解器进行编译时，没有人知道与求解器一起使用的会是哪种湍流模型。实际上，它可能根本没有湍流模型，也可能是任何一种可用的模型。拖曳力模型和两相欧拉求解器也是如此，只是不能使用无拖曳力定律。 +但是，运行时选择机制的背后机理会比本节中介绍的更为复杂。这里重点介绍可以在SchillerNaumann拖曳模型类的源文件中找到的相应宏。我们知道，此拖曳模型是从基类dragModel派生的。为了使运行时选择机制起作用，基类还需要做一些准备。请参阅对运行时选择机制的讨论信息。 +http://openfoamwiki.net/index。 php / OpenFOAM_guide/runTimeSelection_mechanism。希望本节对理解运行时选择机制的一些内部原理有所启发。 +现在借助SchillerNaumann阻力模型作为案例，对OpenFOAM的神奇力量进行分析。列表457和458（第10行和第3行）显示了两行看上去很普通的代码行，正是这些代码将所有的魔力工作都实现了。 + +``` +namespace Foam +{ +class SchillerNaumann +: +public dragModel +{ +public: +//- Runtime type information +TypeName("SchillerNaumann"); +} +} +``` +列表457：SchillerNaumann.H中的相关行 + +``` +namespace Foam +{ +defineTypeNameAndDebug(SchillerNaumann , 0); +} +``` +列表458：SchillerNaumann.C中的相关行 + +##### 57.12.1 part1-TypeName +首先我们检查列表457的第10行 +看起来该调用函数的实际上是一个带有参数的预处理器宏。 宏TypeName其实是在文件typeInfo.H中进行定义。列表459显示了它的定义。 +\ #define宏至少由两部分组成。首先是标识符，然后直到第207行的结尾都是括号中的可选参数列表，或者至少是替换令牌列表。当宏由预处理器进行扩展时，标识符（在这种情况下为TypeName）将被替换为替换标记（参数列表之后的所有指令）。宏不能覆盖多于一行，但借助使用反斜杠（\\）可以将当前行继续延伸到下一行。 + +``` +//- Declare a ClassName () with extra virtual type info +#define TypeName(TypeNameString) +ClassName(TypeNameString); +virtual const word& type() const { return typeName; } +``` +列表459 typeInfo.H中定义的宏 + +因此，该行TypeName（“ SchillerNaumann”）会扩展到—— +``` +ClassName("SchillerNaumann"); +virtual const word& type() const { return typeName; } +``` + +第二行是函数定义。由于此函数是由宏定义的，该函数其实是为了所有在类定义中声明了TypeName宏的类。这正说明了使用预处理器宏的主要原因之一-只需为重复的代码执行一次写入。 +上面列表的第一行本身就是一个宏。列表460显示了扩展ClassName宏所必需的宏定义。 + +``` +//- Add typeName information from argument \a TypeNameString to a class. +// Also declares debug information. +#define ClassName(TypeNameString) +ClassNameNoDebug(TypeNameString); +static int debug +//- Add typeName information from argument \a TypeNameString to a class. +// Without debug information +#define ClassNameNoDebug(TypeNameString) +static const char* typeName_ () { return TypeNameString; } +static const ::Foam::word typeName +``` +列表460：className.H中的两个宏定义 + +然后，进一步扩展宏TypeName("SchillerNaumann") 。 +``` +static const char* typeName_ () { return "SchillerNaumann"; } +static const ::Foam::word typeName +static int debug +virtual const word& type() const { return typeName; } +``` + +由于将TypeName（“ SchillerNaumann”）宏放置到动词+ SchillerNaumann +类的定义中。宏添加了两个函数定义（第一行和最后一行），其中一个是静态函数，以及两个静态变量（两个中心线）。 +类的静态元素（变量或函数）其实是对类的所有实例都只存储在一次的元素。在两阶段欧拉求解器中，可能存在SchillerNaumann类的两个实例-如果为两相都指定了此模型。无论该类的两个实例中是哪一个调用了函数typeName（），它调用的始终同一个函数。在这种情况下（返回类的名称），使用静态函数是很有意义的，并且是实现任务的唯一明智方法。 +TypeName（“ SchillerNaumann”）宏用于创建一个用于返回类名称的函数，以及一个返回类型名称的函数。显然，在设计OpenFOAM时，类名和类型名并不等效。由TypeName（“ SchillerNaumann”）宏创建的变量是包含类型名称的静态变量和名为debug的静态变量。该调试变量控制了第57.11节中介绍的调试机制。 + +##### 57.12.2 Part 2 – defineTypeNameAndDebug +现在检查列表458的第三行，重复如下。 +``` +defineTypeNameAndDebug(SchillerNaumann , 0); +``` + +宏defineTypeNameAndDebug在文件className.H中定义。 + +``` +//- Define the typeName and debug information +#define defineTypeNameAndDebug(Type , DebugSwitch) +defineTypeName(Type); +defineDebugSwitch(Type , DebugSwitch) +``` +列表461：className.H中宏的定义 + +因此该宏被扩展为两个宏 +``` +defineTypeName(SchillerNaumann); +defineDebugSwitch(SchillerNaumann , 0); +``` +列表462展示了扩展为两个宏的相应的必须宏是如何定义的。 + +``` +//- Define the typeName , with alternative lookup as \a Name +#define defineTypeNameWithName(Type , Name) +const ::Foam::word Type:: typeName(Name) + +//- Define the typeName +#define defineTypeName(Type) +defineTypeNameWithName(Type , Type:: typeName_ ()) +//- Define the debug information , lookup as \a Name +#define defineDebugSwitchWithName(Type , Name , DebugSwitch) +int Type::debug (:: Foam::debug:: debugSwitch(Name , DebugSwitch)) + +//- Define the debug information +#define defineDebugSwitch(Type , DebugSwitch) +defineDebugSwitchWithName(Type , Type:: typeName_ (), DebugSwitch); +registerDebugSwitchWithName(Type , Type , Type:: typeName_ ()) +``` +列表462：debugName.H中定义的四个宏 + +因此，该宏扩展为 + +``` +const ::Foam::word SchillerNaumann :: typeName(SchillerNaumann :: typeName_ ()); +int SchillerNaumann :: debug (:: Foam::debug :: debugSwitch(SchillerNaumann , 0)); registerDebugSwitchWithName(SchillerNaumann , SchillerNaumann , SchillerNaumann :: typeName_ ()); +``` + +defineTypeNameAndDebug（SchillerNaumann，0）相应扩展宏的第一行将函数typeName_（）的返回值分配静态变量typeName。这起到了使类名和类型名具有相等值的效果。但此框架的设置方式允许使用不同的名称。 +第二行将调用函数:: Foam :: debug :: debugSwitch（SchillerNaumann，0）调用的返回值分配给静态变量SchillerNaumann :: debug。其中不能将该值直接用于分配到静态变量的原因是，被调用的函数将调试开关添加到字典上了，请参见列表463。 +扩展宏的最后一行调用了另一个宏。列表464显示了registerDebugSwitchWithName的宏定义。 + +``` +int Foam::debug :: debugSwitch(const char* name , const int defaultValue) +{ +return debugSwitches ().lookupOrAddDefault +( +name , defaultValue , false , false +); +} +``` +列表463：将调试开关增加到debug.C中的字典 + +``` +//- Define the debug information , lookup as \a Name +#define registerDebugSwitchWithName(Type ,Tag ,Name) +class add##Tag## ToDebug +: +public ::Foam:: simpleRegIOobject +{ +public: +add##Tag## ToDebug(const char* name) +: +::Foam:: simpleRegIOobject(Foam::debug :: addDebugObject , name) +{ } +virtual ∼add##Tag## ToDebug () +{ } +virtual void readData(Foam:: Istream& is) +{ +Type::debug = readLabel(is); +} +virtual void writeData(Foam:: Ostream& os) const +{ +os << Type::debug; +} +}; +add##Tag## ToDebug add##Tag## ToDebug_(Name) +``` +列表464：debugName.H中对于宏registerDebugSwitchWithName的定义 + +##### 57.12.3 在公园散步：展示其中的一些魔力 +在以上各节中，我们了解了两个非常强大的预处理器宏。那么，它们的存在是为什么呢？ +对于运行时选择机制的实用性，湍流模型是非常突出的示例。在编译时（用户或OpenFOAM开发人员编译求解器时），没人知道最终要在模拟中使用哪种具体的湍流模型。因此，需要在运行时（求解器读取所有案例信息时）决定使用哪种湍流模型。为了避免为每种湍流模型编写求解器，选择以通用方式编写求解器。即在编译求解器时，没有人在乎实际使用的湍流模型，甚至说求解器也不在乎。基类turbulenceModel会告诉编译器和求解器湍流模型是如何工作的，这是用户在编译时需要知道的。 +但是，在运行时需要决定使用哪种湍流模型。幸运的是，OpenFOAM很仔细地设计了相关机制，用户无需担心。但是，在某些情况下，用户想知道当前使用的是哪种湍流模型，可以通过读取案例数据或使用运行时机制的魔术来实现。 +列表465显示了三行代码。该行的目的是打印typeName_（）和type（）函数的返回值。这两种函数分别由第57.12.1节和57.12.2节中介绍的两个宏提供。 + +``` +Info << "Happy printf () debugging:" << endl; +Info << turbulence ->typeName_ () << endl; +Info << turbulence ->type() << endl; +``` +列表465：一些奇妙应用，相应的源代码 + +列表466显示了列表465的三行代码相应的执行结果。列表465中的代码假定湍流建模以其通用形式使用，例如在pimpleFoam的案例中。在此示例中变量turbulence的类型为autoPtr turbulence。 +从输出信息可以看到，turbulence变量确实属于turbulenceModel类型。然而，由于turbulenceModel类是一个抽象基类，没有求解器会实际使用turbulenceModel本身。在实际情况下，求解器使用的是kOmega湍流模型。因此，type（）函数返回的是实际湍流模型的名称。正如前文某些段落所述，这里同样可以看到类名和类型名之间的隐含区别含义。对于具体类的案例来说，他们都是一样的。但对于基类来说，其中的区别有重要的意义。 + +``` +Happy printf () debugging: +turbulenceModel +kOmega +``` +列表466：一些奇妙的应用，相应的输出 + + +#### 57.13；一些关于OpenFOAM并行运行的笔记 +OpenFOAM完全有能力在多个处理器上运行，然而，即使OpenFOAM通过其非常精巧的设计隐藏了程序员在并行运行所需解决的部分或大部分技术问题，仍然有一些问题或可能偶然遇到的陷阱进行考虑。 +首先，花点时间感谢OpenFOAM在并行运行背景下花费大量代价设计的通用性。如果您像任何明智的人一样将自己限制在OpenFOAM的高级数据结构上，那么您基本上可以很轻松地实现并行运行。您可能不需要深入了解拉格朗日粒子是如何遍历子域，或者如何计算子域最边界处场的梯度。 +本节专门讨论一些主题，这些主题为我们这些普通人提供了一次不会搞砸OpenFOAM带来的魔力的机会。出于好奇，还介绍了其他主题，例如打印到Terminal，以及如何实现？ + +##### 57.13.1 Printing to Terminal输出到终端 +当并行运行一个应用程序时，我们可能希望每个并行进程都能够被打印到终端，或者我们可能希望只有主进程会被打印到终端。 OpenFOAM允许您同时执行。 + +###### Hush! The master is printing嘘！主进程正在打印 +如果只希望将主进程打印到终端，则一切正常。只需继续使用我们都熟悉的Info语句即可 + +``` +// this should be printed only once +Info << "Hello OpenFOAM World!" << endl; +``` +列表467：使用Info语句打印到终端 + +列表469显示了列表467的输出，它是并行运行4个进程的结果。 + +###### Parallel print, the choir of processes并行打印，处理器的大合唱 + +如果我们希望所有并行进程都将其输出打印到Terminal，则可以使用Pout语句。该语句允许每个正在运行的进程打印到终端。 + +``` +// this should be printed only once +Info << "Hello OpenFOAM World!" << endl; +// this should be printed once for each running process +Pout << " ... this is process " << pid() << " speaking!" << endl; +``` +列表468：使用Pout打印到终端 + +Pout语句在输出消息时会附上处理器编号，编号从0开始计数，并且0是主进程。注意在列表469中，各个进程打印到终端的顺序不是固定的。 + +``` +Hello OpenFOAM World! +[1] ... this is process 2315 speaking! +[2] ... this is process 2316 speaking! [0] ... this is process 2314 speaking! +[3] ... this is process 2317 speaking! +Finalising parallel run +``` +列表469：列表468的输出结果 + +##### 57.13.2 Condensing values (sums, extrema) +如果想对一个场量的所有值求和，或者想要确定极值（最小值和最大值），我们很可能希望在全局范围内执行这些操作，即在所有子域中并行执行。 +其中需要做的就是执行具体操作，例如在子域中进行求和、通信（例如，所有处理器均向主处理器进行报告），并继续操作（主处理器确定最终结果）。 +由于OpenFOAM的开发人员设计了对用户十分友好的代码库，包括名为gSum，gMin，gMax的语句都可以做到这一点。 + +``` +// correct in serial , however , not so in parallel +Info << "Sum of field U : " << sum(mag(X)*mesh.V()) << endl; +// correct summing -up +Info << "Sum of field U : " << gSum(mag(X)*mesh.V()) << endl; +``` +列表470：计算并打印场的总数到终端 + +如果要对mesh.V（）求和，它是由单元格体积组成的标量场量，当使用gSum（mesh.V（））时，则得出总体积。如果我们使用sum（mesh.V（）），那么得到的将是大致的总体积除以进程数。结果大约是域体积的N分之一，这是因为区域划分并不总是产生大小恰好是全局域N分之一的子域，比如当网格单元总数为偶数、并且划分区域为三部分时。 + +#### 57.14 Math-like syntax in OpenFOAM +OpenFOAM以某种方式在其底层代码设计了集成数学符号的方法，或引用原始论文[66]： + +目的是开发一个C ++类库，能够将复杂的数学和物理模型植入成为高级数学表达式。通过使高级代码尽可能接近标准矢量和张量表示法，可以简化此操作。 + +##### 57.14.1 Pitfalls +###### Operator precedence: use caution with vector products +在数学中，某些运算的优先级高于其他运算216。大多数编程语言也是如此。 C ++中运算符优先级的一般规则严格遵循数学上已有的规则。 +但是，随着OpenFOAM将一些诸如向量、张量的概念扩展到标准C ++的基本数学功能中，运算符优先级的问题变得很有趣。由于C ++中没有向量乘积，并且无法定义全新的运算符，因此必须重载一些现有的运算符以实现向量积。此外，不能更改现有运算符的运算符优先级。 +事实证明，通过按位AND（＆）和按位XOR（^）已实现对内部和外部向量积的重载。使用这些运算符的副作用是，内部和外部向量积的优先级低于基本运算符（例如加法，减法和乘法），这是因为按位AND和按位XOR运算符的优先级很低。 + +对于毫无戒心的用户，当以下公式导致编译错误时，可能会感到惊讶： + +如上所述，意外的运算符优先级导致 + +尝试从上面编译公式，将其转换为OpenFOAM源代码： + +``` +const vector a(1.0, 2.0, 3.0); +const vector b(0.0, 42.0, 666.0); +const scalar c(0.815); +vector d( a & b - c); +``` +列表471：这个简单的向量操作会导致编译错误 + +在潜在的大量编译器警告和错误消息行中，如下的一行是最相关的。它告知我们当前的目标是从向量中减去标量，但相应的代码未被定义。 + +``` + error: no match for ’operator -’ (operand types are ’const vector {aka const Foam::Vector < double >}’ and ’const scalar {aka const double}’) +``` +列表472：编译错误信息 + +因此，建议用户使用大量的括号来避免此类问题。 Hrv Jasak在论坛post218中很好地解释了运算符优先级的问题： + 这是因为operator＆和operator ^实际上是C ++中的二进制C ++运算符，并且该语言不允许您定义自己的运算符或更改运算符优先级。 + 我们选择运算符重载来获得漂亮的语法-使代码看起来更漂亮。因此，我只推荐矢量/张量积中的括号。 +以上情况同样适用于双内积，它采用两个至少二阶张量。通过重载逻辑AND（&&）运算符可实现此运算，但该运算符的优先级甚至比按位二进制运算符还要低。 -- Gitee