在上一节中,由于讨论了用户在运行模拟时接触到的OpenFOAM的部分内容,所以只简要介绍了一些类。以下各节提供了OpenFOAM中一些选定和常用类的更深入的见解和最佳实践。

5.3.1 Dictionary

dictionary类是OpenFOAM用户最早与之交互的类之一。虽然dictionary类接口不是很复杂,但对于新手用户来说,它的某些方面可能并不明显。

从字典中读取数据

读取字典条目值是使用dictionary类的最基本形式。dictionary类接口提供了多个可用于读取数据的方法。getOrDefault方法就是最常用的例子,因为它不仅提供对指定字典条目的读访问,而且如果没有找到条目,其还能够指定一个默认值。这消除了因缺少非关键输入值而导致的运行时错误。

创建的调用形式如:

const auto & solution(mesh.time().solutionDict());
const auto name(solution.getOrDefault<word> ("parameters1"));
const auto vector1(solution.getOrDefault<vector> ("vector1"));

如上面的代码所示,getOrDefault需要的不仅仅是参数的名称才能读取。它还需要一个带有要在字典中查找的数据类型名称的模板参数。

访问目录

目录包含字典的子字典的名称。在下面的示例中,可以访问由Foam::Time类读取的fvSolution文件的目录:

const dictionary &fvSolutionDict(mesh.time().solutionDict());
Info << fvSolutionDict.toc() << endl;

访问子字典

在开发OpenFOAM时,经常会遇到访问子字典的问题,因为在OpenFOAM中开发的自定义类或求解算法由字典文件配置。假设constant目录中的字典文件A包含以下数据:

axis (1 0 0);
origin (5 10 15);
type "modelA";
modelA
{
name "uniform";
}

为了访问A文件中的子字典modelA,可以使用下面的代码:

const dictionary dict(fileName("constant/A"));
const auto& subDict(dict.subDict("modelA"));

5.3.2 量纲类型

Dimensioned类型将单位附加到标量、矢量和张量上。其扩展了张量运算以包括度量单位:OpenFOAM中的量纲。如速度与动量都是矢量,但它们不能进行代数相加,因为他们的量纲不同。OpenFOAM中的量纲检查是通过类模板dimensioned<Type> 来实现的。

dimensioned<Type> 是一个包装器(适配器)类,其将张量算法的计算委托给经过包装的tensor类,并将量纲检查工作委托给经过包装的dimensionSet对象。通过研究dimensioned<Type> 类模板的算术运算符,可以实现如下所示的+=运算符:

template<class Type> 
void Foam::dimensioned<Type>  :: operator+=
(
    const dimensioned<Type>  &dt
)
{
    dimensions+ += dt.dimensions_;
    value_ += dt.value_;
}

可以看出,有两个算术运算正在执行:一个是对量纲(单位)的运算,另一个是对张量的数值的运算。dimensionSet类的算术运算符负责量纲的检查过程。dimensionSet类的+=算术运算符的源代码如下所示:

bool Foam::dimensionSet :: operator+= (const dimensionSet &ds) const
{
    if(dimensionSet::debug && *this != ds)
    {
        FatalErrorIn("dimensionSet::operator+=(
            const dimensionSet&) const")
            << "Different dimensions for +=" << endl
            << " dimensions : " << *this << " = "
            << ds << endl << abort(FatalError);
    }
}

这提供了足够的信息来准确总结量纲检查是如何执行的:如果量纲检查已打开,则当对不同(!=)量纲的集合执行加法操作时,将生成一个致命错误,该错误将中止程序执行。dimensionSet类将量纲实现为一组物理度量的整数指数,如下所示:

//- Define an enumeration for the names of the dimension exponents
enum dimensionType
{
    MASS, // kilogram kg
    LENGTH, // metre m
    TIME, // second s
    TEMPERATURE, // Kelvin K
    MOLES, // mole mol
    CURRENT, // Ampere A
    LUMINOUS_INTENSITY // Candela Cd
};

如下面的代码所示,量纲检查运算符是根据相等运算符==来执行的。检查过程会迭代操作dimensionSet的量纲。通过循环测试,如果量纲指数之间的差异的大小足够大(> smallExponent),可以认为集合不相等。

bool Foam::dimensionSet::operator==(const dimensionSet& ds) const
{
    for (int Dimension=0; Dimension < nDimensions; ++Dimension)
    {
        if(mag(exponents_[Dimension] - ds.exponents_[Dimension])>  smallExponent)
        {
            return false;
        }
    }
    return true;
}

SmallExponent的值是类静态变量:

static const scalar smallExponent;

其被初始化为SMALL,对于双精度scalar,其值为

通过将$WM_PROJECT_DIR/etc/control Dict中dimsionedSet类的debug标志设置为ON,可以打开量纲检查:

DebugSwitches
{
    Analytical 0;
    APIdiffCoefFunc 0;
    ...
    dictionary 0;
    dimensionSet 1;
    mappedBase 0;
    ...

对于+=运算符,给出的结论与其他量纲张量算术运算完全相同。默认情况下,量纲检查处于激活状态,通常不应在OpenFOAM中将其取消激活。即使自定义应用程序以无量纲的形式实现公式,这些方程式也会使用应该是无量纲的数字进行缩放。量纲检查系统可以检查例如雷诺数:

auto Re = Foam::mag(U) * L/mu;
Info << Re.dimensions() << endl;

雷诺数实际上是无量纲的,在这种情况下,输出的是仅包含零的dimensionSet。然而,如果在计算(例如L)时出错,依赖SI测量单位系统的尺寸检查将发现这一点。因此,即使方程是以无量纲形式编写的,也应该使用量纲检查来检查缩放错误。

OpenFOAM中的许多具有复杂名称的类模板,都提供了typedef。在c中,typedef关键字允许程序员定义较短和更简洁的类型名称。尽管这个typedef可能不会比原来短很多,但它至少保存了部分的输入。在dimensioned的情况下,重点不是名称长度,而是代码样式——dimensionedVector是camel case中的一个名称,用于OpenFOAM应用程序级代码中的类型。

dimensioned<type> 最常用的同义词是dimensionedScalardimensionedVector,在以下示例中使用。此外,请注意,dimensioned类型的构造方式比目前描述的方式稍微复杂一些。标注类型不需要张量值和dimension集合,而是需要一个附加参数:name。例如,如果二维矢量对象按以下方式构造:

dimensionedVector velocity
(
    "velocity",
    dimLength/dimTime,
    vector(1,0,0)
);

dimensionVector momentum
(
    "momentum",
    dimMass * (dimLength / dimTime),
    vector(1,0,0)
);

请注意,有一些预定义的dimensionSet对象用于初始化速度和动量向量对象:dimLength、dimTime和dimMass。它们是常量和全局对象,可以在源文件DimensionSets.C中找到:

const dimensionSet dimless(0,0,0,0,0,0,0);
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0);
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0);
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0);
const dimensionSet dimTemperature(0, 0, 0, 1, 0, 0, 0);
const dimensionSet dimMoles(0, 0, 0, 0, 1, 0, 0);
const dimensionSet dimCurrent(0, 0, 0, 0, 0, 1, 0);
const dimensionSet dimLuminousIntensity(0, 0, 0, 0, 0, 0, 1);

可以使用基本的物理量纲来构造复杂的量纲(例如N=kgm/s2),所以在定义自己的dimensionSet时,使用全局预定义的尺寸集对象来提高代码的可读性。

继续讨论量纲类型的算术,并将速度添加到动量中,如下所示:

momentum += velocity;

将会导致下面的错误提示:

-->  FOAM FATAL ERROR:
Different dimensions for +=
dimensions : [1 1 -1 0 0 0 0] = [0 1 -1 0 0 0 0]
From function dimensionSet::operator+=(const dimensionSet&) const
in file dimensionSet/dimensionSet.C at line 179.
FOAM aborting

OpenFOAM中的量纲检查过程在运行时执行。因此,在量纲操作中包含错误的代码仍然会被编译,但不会运行。OpenFOAM中的量纲检查错误可以通过调试器(例如gdb)轻松调试,前提是OpenFOAM和自定义代码是在调试模式下编译的。

练习:使用前面提到的调试标志关闭量纲检查,并运行动量和速度的算术加法

5.3.3 智能指针

指针是一种特殊的变量,其存储对象的内存地址,可以用来引用该对象。在c中,可以通过值或引用来传递对象,后者对于较大对象来说明显更快。这不仅更快,而且在内存使用方面更加保守,因为没有数据被临时复制。与只使用指针相比,按值返回较大的对象并将其作为函数参数传递在计算上要昂贵得多。

如下面的代码:

result = function(input);

可以改成以非常量引用传递给函数的result参数:

function(result,input);

离散化算法(运算符)首选第一种选择,因为它们通常由算术表达式组成,用于建立数学模型。例如,考虑取自interFoam求解器的动量守恒方程代码:

fvVectorMatrix UEqn
(
    fvm::ddt(rho,U)
    + fvm::div(rhoPhi,U)
    + turbulence-> divDevRhoReff(rho,U)
);

作用于场rho、U及rhoPhi的算子之和将是一个系数矩阵(fvVectorMatrix)。因此,上述代码必须满足以下几点:

  • fvVectorMatrix必须是可拷贝构造的
  • 所有操作符必须返回一个fvVectorMatrix
  • fvVectorMatrix的加法运算符必须返回一个fvVectorMatrix

如果运算符ddt和div被实现为具有可修改参数的函数,则不可能轻松地编写数学模型(通常称为方程模拟)。函数返回的矩阵非常大,因此按值返回它们会带来创建临时对象的代价。借助于OpenFOAM智能指针,OpenFOAM避免以不依赖编译器优化的显式方式创建不必要的复制操作。智能指针在插值函数中初始化,并通过值返回。涉及方程离散化的其他操作,如对流格式,也是如此。例如,高斯对流格式的fvmDiv散度算子初始化了这样一个智能指针(tmp<fvMatrix<Type> > ),如下所示 :

tmp<fvMatrix<Type> >  tfvm
(
    new fvMatrix<Type> 
    (
        vf,
        faceFlux.dimensions()* vf.dimensions()
    )
);

然后它执行计算,负责定义所示系数矩阵中的元素.

fvm.lower() = -weights.internalField()*faceFlux.internalField();
fvm.upper() = fvm.lower() + faceFlux.internalField();
fvm.negSumDiag();
forAll(vf.boundaryField(), patchI)
{
     const fvPatchField<Type> & psf = vf.boundaryField()[patchI];
    const fvsPatchScalarField& patchFlux = faceFlux.boundaryField()[patchI];
    const fvsPatchScalarField& pw = weights.boundaryField()[patchI];
    fvm.internalCoeffs()[patchI] = patchFlux*psf.valueInternalCoeffs(pw);
    fvm.boundaryCoeffs()[patchI] = -patchFlux*psf.valueBoundaryCoeffs(pw);
}

if (tinterpScheme_().corrected())
{
    fvm += fvc::surfaceIntegrate
    (
        faceFlux*tinterpScheme_().correction(vf)
    );
}

然后,在函数开始时初始化的智能指针被按值返回。

return tfvm;

按值返回有限体积矩阵(tfvm)的智能指针的结果是对智能指针进行复制操作。然而,智能指针的拷贝和整个矩阵对象的拷贝之间有一个显著的区别:指针的值只是矩阵的地址,而不是整个矩阵本身。这种方法大大增加了执行的效率。

避免不必要的复制操作有一个常用的简短名称:复制消除。在c编程语言中,可以通过不同的方式来执行拷贝消除:编译器优化(返回值优化,命名为RVO),表达式模板(ET),或者使用c 11语言标准提供的右值引用和移动语义。

通常,当使用指针时,它们指向使用运算符new在堆上创建的对象:

someType *ptr = new someType(argments...);

由于c编程语言不支持自动垃圾收集,因此只能由程序员负责释放资源。

因此,每次调用new操作符之后都需要相应地调用delete操作符。当然,对delete操作符的调用必须放在适当的位置,例如类析构函数。这就可能导致程序员忘记删除指针,从而导致内存泄漏。作为错误的另一个来源,访问已删除指针引用的内存部分将导致未定义的行为。这两个问题都会发生在运行时,通常是难以发现和调试的错误源。为了避免这两个问题,需要避免直接处理原始指针。在c中,原始指针的处理已被称为资源获取初始化(RAII)的习惯用法所取代。

RAII习惯用法指出,原始指针需要封装到一个类中,该类的析构函数负责删除指针并在代码中的适当位置释放资源。在这样的类中调整原始指针,并为调整后的原始指针提供不同的功能,导致了所谓的智能指针的发展。不同的智能指针存在并提供不同的功能。OpenFOAM实现了两个这样的智能指针:autoPtrtmp

只要环境变量是通过源码库的etc/ bashrc配置脚本来设置的,例子的应用程序代码就可以在$PRIMER_EXAMPLES_SRC/applications/test/文件夹下获得。那些有兴趣按照接下来的章节介绍的教程一步步学习的读者,需要创建一个新的可执行的应用程序:testSmartPointers。在OpenFOAM中创建新的应用程序是很简单的,因为有脚本可以为应用程序生成文件框架。要创建一个新的应用程序,需要选择一个目录来放置应用程序的代码,并执行以下命令:

mkdir testSmartPointers
cd testSmartPointers
foamNew source App testSmartPointers
sed -i 's/FOAM/APPBIN/ FOAM_USER_APPBIN/g' Make/options

最后一行将应用程序二进制文件的放置目录从平台目录替换为用户应用程序二进制文件目录。

将你自己的应用程序构建到FOAM_USER_APPBIN中是一个很好的做法,这需要在Make/options构建配置文件中指定。在本节中,使用的是gcc编译器。在阅读文本中关于特定编译器标志的内容时,要注意到这一点

5.3.3.1 使用autoPtr智能指针

为了说明如何使用autoPtr,需要为使用的示例定义一个新类。该类应该在每次调用其构造函数和析构函数时通知用户。在本例中,它继承自基本字段类模板字段,并命名为infoField。可以使用任何其他类,因为编译器对复制操作执行的优化不依赖于对象的大小。

首先,可以如下所示定义infoField类模板:

template<typename Type> 
class infoField : class Field<Type> 
{
    public:
        infoField():Field<type> ()
        {
            Info << "empty constructor" <<endl;
        }

        infoField(const infoField& other):Field<Type> (other)
        {
            Info << "copy constructor" << endl;
        }

        infoField(int size, Type value) : Field<Type> (size,value)
        {
            Info << "size,value constructor" << endl;
        }

        ~infoField()
        {
            Info << "destructor" << endl;
        }

        void operator=(const infoField & other)
        {
            if(this != &other)
            {
                Field<type>  :: operator=(other);
                Info << "assignment operator" << endl;
            }
        }
};

该类模板继承自Field ,并使用以下内容:

  • 空构造函数
  • 拷贝构造函数
  • 析构函数
  • 赋值运算符

每次使用这些函数时,一个Info语句会向标准输出流发出特定调用的信号。这样做的目的只是为了获得哪个函数被执行的信息。

为了继续这个例子,需要定义一个函数模板,其按值返回一个对象。

template<typename Type> 
Type valueReturn(Type const & t)
{
     // One copy construction for the temporary.
    Type temp = t;
    // ... operations (e.g. interpolation) on the temporary variable.
    return temp;
}

为了减少不必要的输入,例子中使用的类型名称被缩短了。

// 简写类型名
typedef infoField<scalar>  infoScalarField;

在主函数中,实现了以下几行:

Info << "Vlaue construction:";
infoScalarField valueConstructed(1e7, 5);
Info << "Empty construction:";
infoScalarField assignedTo;
Info << "Function call" << endl;
assignedTo = valueReturn(valueConstructed);
Info << "Function exit" << endl;

用Debug或Opt选项编译和执行应用程序,将产生完全相同的结果,尽管Debug选项可以打开编译器的优化。正如本节开头所提到的:编译器非常善于识别这样一个事实,即返回的临时对象在赋值后就被丢弃了。不使用几何场的另一个好处是,这个小的示例程序不需要在OpenFOAM仿真案例目录下执行。它可以直接从存储代码的目录中调用。执行该程序将产生以下输出:

?>  testSmartPointers
Value construction : size, value constructor
Empty construction : empty constructor
Function call
copy constructor
assignment operator
destructor
Function exit
destructor
destructor

单独考虑输出的每一行,同时将其与valueReturn函数代码进行比较,有趣的是发现临时返回变量从未被构造。在函数的返回语句中缺少一个拷贝构造,因为该函数是通过值返回的,同时也缺少一个相应的析构器调用。这种行为的原因是编译器自动进行的拷贝消除优化,即使在调试模式下也是如此。为了消除这种优化,可以在Make/options中添加一个额外的编译器标志:

EXE_INC = \
    -I(LIB_SRC)/finiteVolume/lnInclude \
    -fno-elide-constructors
EXE_LIBS = \
    -lfiniteVolume

gcc编译器标志-fno-elide-constructors将阻止编译器进行优化,这对跟踪valueReturn中发生的事情很重要,否则会被编译器优化掉。

注意:记得在执行wmake之前调用wclean。对Make/options文件的修改不会被wmake构建系统识别为需要重新编译的源代码修改。

用上面定义的选项文件编译和运行应用程序,结果是以下输出:

?>  testSmartPointers
Value construction : size, value constructor
Empty construction : empty constructor
Function call
copy constructor # Copy construct tmp
copy constructor # Copy construct the temporary object
destructor # Destruct the tmp - exiting function scope
assignment operator # Assign the temporary object to assignedTo
destructor # Destruct the temporary object
Function exit
destructor
destructor

输出显示了临时对象的不必要的创建和删除。通过在函数返回的地方执行对象的就地构建来避免临时对象的拷贝,已经成为编译器的一个标准选项。经常发生的情况是,即使在Debug模式下用编译器标记-O0 -DFULLDEBUG抑制优化,也不能禁用这个功能。

为了禁用构造函数复制的优化,必须在Make/options中明确传递编译器标志-fno-elide-constructors。

上面的例子有一个主要目的:表明在表达式中没有必要使用autoPtr智能指针,因为在表达式中无论如何都会返回一个命名的临时对象。在现代编译器上,即使在调试模式下编译(对于OpenFOAM来说,这意味着将$WM_COMPILE_OPTION设置为调试),这种优化也是默认开启的。

一些最突出的使用autoPtr的例子是在采用RTS的模型中,例如湍流模型和fvOptions。特定的基类是基于用户在字典中指定的关键字而实例化的。以turbulenceModel类为例,基类被用作autoPtr的模板参数,例如在求解器应用程序中,RTS可以将特定的湍流模型实例化到autoPtr中。RTS允许类用户在运行时将类层次中的特定类的对象进行实例化。以这种方式实例化对象使c能够通过基类指针或引用访问派生类对象。这通常被称为动态多态性。

本书不可能涵盖OpenFOAM使用的c语言标准的所有细节。只要遇到听起来不熟悉的c结构,就应该到其他地方去查找。

湍流模型通常在特定求解器的createFields.H中实例化,它被包含在时间循环开始之前。相关的代码行如下所示:

autoPtr <incompressible::RASModel>  turbulence
(
    incompressible::RASModel::New(U,phi,laminarTransport)
);

很明显,一个autoPtr被用来存储湍流模型。如果使用原始指针来访问RASModel对象,那么在求解器的代码中还需要另一行来删除原始指针,以便在求解器结束时去分配内存。由于autoPtr是一个智能指针,这种删除是在autoPtr析构器中自动执行的。原始指针方法的代码(幸运的是在OpenFOAM中没有任何地方使用)看起来是这样的。

incompressible::RASModel* turbulence = incompressible::RASModel::New(U, phi, laminarTransport)

在求解器代码的末尾,需要有以下几行:

delete turbulence;
turbulence = NULL;

当然,对于一个单一的指针来说,添加行来释放资源可能不是问题。然而,RTS用于:传输模型、边界条件、fvOptions、离散格式、插值格式、梯度格式等等。所有这些对象与场和网格这样的东西相比都很小,所以它们可以通过值来返回吗?从效率的角度来看--是的,特别是考虑到所有现代编译器都实现了返回值优化(RVO)。从灵活性的角度来看--没有机会这样做,因为动态多态性依赖于通过指针或引用的访问。保持运行时的高灵活性,意味着要依赖指针或引用。

当RTS被用来在运行时选择对象时,使用autoPtr或tmp是必要的。

下面的例子研究了 autoPtr 接口和它基于所有权的复制语义。autoPtr拥有它所指向的对象。这是预期的,因为RAII要求智能指针处理资源释放,所以程序员不必这样做。因此,创建autoPtr的副本是很复杂的--复制一个autoPtr会使原来的autoPtr失效,并将对象的所有权转移到副本中。要看这是如何工作的,请考虑下面代码段的主函数。

int main()
{
    // Construct the infoField pointer
    autoPtr<infoScalarField>  ifPtr (new infoScalarField(1e06, 0));
    // Output the pointer data by accessing the reference -
    // using the operator T const & autoPtr<T> ::operator()
    Info << ifPtr() << endl;
    // Create a copy of the ifPtr and transfer the object ownership
    // to ifPtrCopy.
    autoPtr<infoScalarField>  ifPtrCopy (ifPtr);
    Info << ifPtrCopy() << endl;
    // Segmentation fault - accessing a deleted pointer.
    Info << ifPtr() << endl;
    return 0;
}

需要注意的是,一旦对autoPtr对象进行了拷贝,所指向的对象的所有权就会被转移。注释掉导致段故障的那一行,结果应用程序的输出看起来像这样。

?>  testSmartPointers
size, value constructor
1000000{0}
1000000{0}
destructor

很明显,尽管autoPtr对象是通过值传递,也只会引用infoScalarField类的一个构造函数和析构函数。因此,autoPtr可以被用来节省不必要的复制操作。

5.3.3.2 使用tmp智能指针

tmp智能指针通过执行引用计数来防止对象的不必要的复制。引用计数是指同一对象被传递的过程。在这种情况下,引用计数是通过refCount类进行的。这个数据结构是任何应该存储在tmp智能指针中的类的重要基类。它被智能指针包裹着,每次对智能指针进行复制或赋值时,这个对象的引用数量都会增加。

遵照RAII,临时工指针的析构器负责销毁被包装的对象。该析构器检查该对象在当前作用域中的refCount数量。如果这个数字大于0,析构器就简单地将refCount减少1,并允许该对象继续生存。一旦在被包装对象的 refCount 达到 0 的情况下调用析构器,析构器就会删除该对象。

注意:tmp智能指针的定义可以在$FOAM_SRC/OpenFOAM/memory/tmp中找到。

使用tmp智能指针时有一个问题:指针类模板不负责计算引用。引用计数是由封装的类型来实现的。在检查类的析构器时,可以很容易地检查这一点。

template<class T> 
inline Foam::tmp<T> ::~tmp()
{
     if (isTmp_ && ptr_)
    {
         if (ptr_-> okToDelete())
        {
            delete ptr_;
            ptr_ = 0;
        }
        else
        {
             ptr_-> operator--();
        }
    }
}

ptr_属性是一个类型为T的对象的包装好的原始指针,解构器试图访问两个成员函数。

  • T::okToDelete()
  • T::operator--()

因此,被包装的对象必须遵守一个特定的类接口。测试这个问题的另一个方法是尝试用一个琐碎的类来使用tmp,这个类可以定义如下。

class testClass {};

然后试图将这个类包装成一个tmp智能指针。

tmp<testClass>  t1(new testClass());

上述代码导致以下错误。

tmpI.H:108:9: error: ‘class testClass’
has no member named ‘okToDelete’ if (ptr_-> okToDelete())
tmpI.H:115:13: error: ‘class testClass’
has no member named ‘operator--’ ptr_-> operator--();

通过这个错误,编译器抱怨上述成员函数没有被testClass实现的事实。

由于OpenFOAM使对象执行引用计数,所以它被封装在一个此类对象继承的类中,名为refCount。refCount类实现了引用计数器和相关的成员函数。

为了在OpenFOAM中使用tmp ,被包装对象的类型T应该继承于refCount。

如果testClass被修改为继承于refCount。

class testClass : public refCount {};

然后可以用tmp来包装。为了了解引用计数是如何工作的,以及如何避免不必要的对象构造,tmp可以与infoScalarField类一起使用,该类在描述autoScalarField的例子中使用过,然后可以用tmp包装。为了了解引用计数是如何工作的,以及如何避免不必要的对象构造,tmp可以与infoScalarField类一起使用,该类在描述autoPtr的例子中使用。refCount类允许用户通过成员函数refCount::count()获得当前引用计数的信息,这在下面的例子中使用。人工作用域被用来减少临时工对象的寿命,以便它们的析构器被调用。在普通程序代码中,当出现嵌套函数调用或循环时,会出现这种情况。下面是例子的代码。

tmp<infoScalarField>  t1(new infoScalarField(1e06, 0));
Info << "reference count = " << t1-> count() << endl;
{
    tmp<infoScalarField>  t2 (t1);
    Info << "reference count = " << t1-> count() << endl;
    {
        tmp<infoScalarField>  t3(t2);
        Info << "reference count = " << t1-> count() << endl;
    } // t3 destructor called
    Info << "reference count = " << t1-> count() << endl;
} // t2 destructor called
Info << "reference count = " << t1-> count() << endl;

这将导致以下命令行输出。

> ? testSmartPointers
size, value constructor
reference count = 0
reference count = 1
reference count = 2
reference count = 1
reference count = 0
destructor

来自infoField类的单个构造器和相应的析构器输出显示,尽管tmp对象被作为函数参数按值传递,但只进行了一次构造和破坏。

5.3.4 体积场

第10章深入讨论了边界场和边界条件,包括它们背后的理论。

volume field是那些通常用于存储单元中心的场变量值。根据该场变量存储的张量类型,可以使用volScalarField、volVectorField或volTensorField。也有表面场和点场,它们分别在面中心和单元角点存储物理场值。请注意,所有提到的场的构造都是相似的,与它们的类型无关。

在OpenFOAM中,将数值映射到非结构化网格单元的场由通用的GeometricField类模板实现的。具体的场,如体积场、表面场和类似的其他物理场,则通过实例化GeometricField类模板和特定的模板参数,以具体类的形式生成。正如在关于量纲类型的章节中所指出的,OpenFOAM中用于场的类型名称为了方便使用typedef关键字而被缩短。

涉及volScalarField等场变量的编译错误会导致c模板错误,这些错误相当长,读起来不直观。因为这些字段是GeometricField模板的实例,可以调查GeometricField类的模板源代码,以获得更多关于错误的信息。

在能够使用一个物理场的对象之前,它必须被构造,就像任何其他对象一样。类接口有几个重载构造函数,根据特定的情况可能会派上用场。简单的构造函数是复制构造函数。

const volScalarField pOld(p);

顾名思义,它构造了一个原始数据的副本,其类型是相同的。然而,为了使用这个构造函数,首先必须有一个volScalarField。有两种方法可用于此:要么从文件中读取场变量的数据,要么从头开始生成它。基于输入文件的场的初始化是相当直接的,在OpenFOAM的每个求解器或处理程序的createFields.H文件中可以找到这样的例子。为此,IOobject被用来执行访问该文件的文件系统级输入输出,并作为一个封装类来使用。它是物理场所需要的,以便从文件中构建。

volScalarField T
(
    IOobject
    (
        "T",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

上面的代码根据T文件的内容初始化volScalarField,T文件必须存在于时间目录中(或者至少是0-目录)。当然,T文件必须有适当的格式,以便在IOobject的构造函数中使用,然后被用来构造volScalarField。如果T文件不存在,由于IOobject::MUST_READ指令,代码的执行将停止。其他指令也可以选择,这取决于特定的需要。剩下的两个选项是READ_IF_PRESENT和NO_READ。第一个选项只在文件存在于任何一个时间目录中时才读取,否则就通过提供的默认值构建。NO_READ从不从文件中读取任何东西,只是构建物理场。

在某些情况下,需要在不从文件中读取数据的情况下构建物理场。通量场就是一个很好的例子。

surfaceScalarField phi
(
     IOobject
    (
        "phi",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    linearInterpolate(rho*U) & mesh.Sf()
);

如果文件存在,则phi由字段数据本身构建。否则直接从速度场计算通量。

这显示了T和phi的构造的不同;两者都以IOobject作为第一个参数,但T的构造以polyMesh作为第二个参数,而phi则以surfaceScalarField作为第二个参数。如果你仔细看看$FOAM_SRC/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H,你会发现GeometricField类有大量的重载构造函数,它是OpenFOAM中所有场的基类。你可以根据特定的需要使用所有的构造函数。

5.3.4.1 访问网格数据

特定的单元可以通过调用特定场的访问操作符Type&operator[](const label cellI)并将所需的网格标签作为参数来处理。当对任何物理场进行算术运算时,不建议在网格之间进行运算。相反,应该使用物理场算术运算符。

在应用级程序代码中使用场变量的循环会降低代码的可读性,并可能导致计算效率的显著下降。

总之,这种选择特定网格的方法只能在需要使用网格子集进行计算时使用。下面的代码显示了一个如何选择压力场和速度场的网格的例子。

labelList cellIDs(3);
cellIDs[0] = 1;
cellIDs[1] = 42;
cellIDs[3] = 39220;

forAll(cellIDs, cI)
{
    Info << U[cellIDs[cI]] << tab << p[cellIDs[cI]] << endl;
}

5.3.4.2 访问边界数据

正如第1章所述,内部场值与边界场值是分开的。这种网格中心和边界(面中心)值的逻辑分离是由支持FVM的数值插值方式所定义的。

将边界场值与内部场值分开,对OpenFOAM中算法的并行化方式有重要影响。数值操作在OpenFOAM中使用数据并行化,域被分解成子域,数值操作在每个独立的子域上执行。因此,许多并行进程被执行,需要跨越分解(处理器)的边界相互通信。将进程边界建模为边界条件的结果是,在OpenFOAM中基于owner-neighbor寻址的所有数值操作的自动并行化。非结构化的FVM离散化的自动并行化是OpenFOAM的一个值得注意的特点。

下面的例子显示了如何访问边界补丁出口的体积标量场压力的边界场值。关于边界场及其与内部场的区别的更多信息,请参考第10章。出口处的边界场可以根据映射到边界网格patch(边界网格的一个子集)名称的边界ID来找到。

const label outletID(mesh.boundaryMesh().findPatchID("outlet"));

然后使用GeometricField:: boundaryField成员函数访问边界场。

const scalarField& coutletPressure = pressure.boundaryField()[outletID];

体积场pressure有一个成员函数boundaryField,它返回一个指向边界场的指针列表。出口边界字段在该指针列表中的位置由出口ID标签(索引)定义。上面的代码片段将边界物理场设置为常数引用outletPressure,然而对边界物理场的非常数访问是由成员函数提供的。这可以通过观察GeometricField类模板的声明来确认。

GeometricBoundaryField &boundaryFieldRef();

GeometricBoundaryField是一个类模板,它的参数与GeometricField相同,这个类模板的定义被放置在GeometricField类接口的public部分。