本案例演示在OpenFOAM中编写工具程序的基本过程。案例实现一个小功能，用于输出最大的网格变化率。

1 程序框架

OpenFOAM提供了foamNewApp帮助我们快速地创建程序框架。

执行下面的命令：

cd $FOAM_RUN/Demo
foamNewApp volumeRatioCheck
cd volumeRatioCheck

程序结构如下图所示。

可以看到程序自动创建了一个名为volumeRatioCheck.C的源文件，以及Make文件夹及其内部的files与options文件。

默认情况files文件中的内容如下所示，这样编译后的可执行文件就被放到了$FOAM_USER_APPBIN中。

这里修改一下files文件，以使编译后的可执行文件放在Demo目录下，将其修改为：

volumeRatioCheck.C
EXE = volumeRatioCheck

查看options文件，其内容如下图所示，此文件无需修改。

可以查看源文件volumeRatioCheck.C的内容，如下所示。

foamNewApp已经帮助我们将程序的大致框架搭建起来了。现在只需要修改volumeRatioCheck.C文件即可。

2 原始程序

本程序的目的是找出计算区域内最大的网格体积比率，这可以通过对所有网格遍历来实现，其中需要能够获取网格的数据。编写程序代码如下：

#include "fvCFD.H"
// * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
    // 初始化案例，会检查程序运行当前路径是否为OpenFOAM案例的根路径
    #include "setRootCase.H"
    // 创建一个Time对象runTime
    #include "createTime.H"
    // 创建一个fvMesh对象mesh，后面所有的网格信息获取全部通过mesh对象来完成
    #include "createMesh.H"
 
    // mesh.C()返回所有网格体心坐标列表；size()函数返回网格数量
    const label c = mesh.C().size();
    Info << "拥有网格数量为:" << c << endl;
 
    // 利用函数cellCells()返回网格的neighbour列表
    const labelListList &neighbour = mesh.cellCells();
 
    // 定义两个临时变量，用于存储结果
    // ratios为比率列表，将所有的体积比都放到这个列表中
    // volumeRatio用于临时存储体积比
    List ratios(0);
    scalar volumeRatio = 0;
 
    // 遍历所有的neighbour，注意neightbour的类型
    forAll(neighbour, cellI)
    {
        // 定义一个整数列表n，存储所有的网格的neighbour编号列表
        List n = neighbour[cellI];
        // 利用mesh.V()获取网格的体积
        const scalar cellVolume = mesh.V()[cellI];
 
        // 在n中进行遍历
        forAll(n, i)
        {
            label neighbourIndex = n[i];
            // 得到neighbour网格的体积
            scalar neighbourVolume = mesh.V()[neighbourIndex];
            if (neighbourVolume >= cellVolume)
            {
                // 得到体积比，并将其放入ratios中
                volumeRatio = neighbourVolume / cellVolume;
                ratios.append(volumeRatio);
            }
        }
    }
    // 利用max函数得到列表中的最大值
    Info << "最大体积比为：" << max(ratios) << endl;
 
    // 下面的程序是foamNewApp自动生成的，输出一些时间信息，可以保留
    Info << nl;
    runTime.printExecutionTime(Info);
    Info << "Endn"
         << endl;
 
    return 0;
}

程序代码编辑完毕后，可以使用wmake进行编译，如下图所示。

不过本程序的执行需要基于OpenFOAM案例文件结构以及网格，因此这里从tutorials中拷贝一个算例到当前路径下，并生成网格。

利用下面的命令测试运行：

cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily .
cd pitzDaily
blockMesh
../volumeRatioCheck

运行结果如下图所示。

可以看到程序统计得到了计算域内的最大网格体积比。

3 程序改进

3.1 改进版本1

前面的程序虽然能正常工作，但是在网格数量非常多的时候，程序效率很低。其主要原因在于代码：

ratios.append(volumeRatio);

此行代码将计算得到的volumeRatio的值插入到列表ratios中。然而列表元素的插入操作非常消耗时间，因为其需要先复制列表，当原始列表中元素较多时，复制工作性能非常低下。

这里可以换一种思路，比如一开始就将列表内存开辟完毕，在有数据需要插入时，利用赋值操作代替插入操作，可以节省大量的时间。

将程序代码改造为：

#include "fvCFD.H"
int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
 
    const label c = mesh.C().size();
    Info << "拥有网格数量为:" << c << endl;
    const labelListList &neighbour = mesh.cellCells();
 
    // 定义一个数量明确的List，其数量为计算域内的所有内部面的数量
    // mesh.Cf().size()得到内部面数量
    label len  = mesh.Cf().size();
    scalar initial = 0;
    List ratios(len,initial);
    label counter = 0;
    scalar volumeRatio = 0;
 
    forAll(neighbour, cellI)
    {
        List n = neighbour[cellI];
        const scalar cellVolume = mesh.V()[cellI];
 
        forAll(n, i)
        {
            label neighbourIndex = n[i];
            scalar neighbourVolume = mesh.V()[neighbourIndex];
            if (neighbourVolume >= cellVolume)
            {
                volumeRatio = neighbourVolume / cellVolume;
                // 采用直接复制的方式可以避免列表复制带来的额外开销
                ratios[counter] = volumeRatio;
                counter++;
            }
        }
    }
    Info << "最大体积比为：" << max(ratios) << endl;
 
    Info << nl;
    runTime.printExecutionTime(Info);
 
    Info << "Endn"
         << endl;
    return 0;
}
 

改造后的代码执行结果如下图所示。程序执行时间从0.17s减少为0.08s，运行时间为之前的一半。

改造后的代码的确可以节省大量的时间，尤其对大量网格进行操作时更加明显。

3.2 改进版本2

然而这个代码依然可以继续改进。在本代码中我们定义了一个用于存放体积率的列表ratios，此列表存储的元素数量与网格数量一致，比如有1000万网格，那么此列表元素个数也有1000万个，显然很消耗内存。

我们其实并不需要将所有的volumeRatio保存在列表中，而只需要定义一个临时变量，在计算体积率之后与此变量进行比较，若比此变量值大，则更新临时变量的值为新的体积率。这样的话可以不需要定义一个庞大的列表，从而节省内存。

改造后代码如下所示：

#include "fvCFD.H"
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
 
    const label c = mesh.C().size();
    Info << "拥有网格数量为:" << c << endl;
    const labelListList &neighbour = mesh.cellCells();
 
    scalar volumeRatio = 0;
    // 定义临时变量tmpRatio
    scalar tmpRatio = 0;
 
    forAll(neighbour, cellI)
    {
        List n = neighbour[cellI];
        const scalar cellVolume = mesh.V()[cellI];
 
        forAll(n, i)
        {
            label neighbourIndex = n[i];
            scalar neighbourVolume = mesh.V()[neighbourIndex];
            if (neighbourVolume >= cellVolume)
            {
                volumeRatio = neighbourVolume / cellVolume;
                // 如果volumeRatio的值比tmpRatio的值大，则更新tmpRatio的值
                if (volumeRatio > tmpRatio)
                {
                    tmpRatio = volumeRatio;
                }
            }
        }
    }
    // 直接输出最大体积率
    Info << "最大体积比为：" << tmpRatio << endl;
 
    Info << nl;
    runTime.printExecutionTime(Info);
 
    Info << "Endn"
         << endl;
 
    return 0;
}
 

程序执行结果如下图所示。这里的改进并没有提高计算效率，但是对于减少内存消耗是很有帮助的，尤其是在网格数量很多时表现更为明显。

在编写新代码时，重要的是要考虑要计算哪些数据以及如何存储这些数据。计算流体力学和相关模拟通常非常耗费内存和处理器，对程序代码的微小改进有时可以显著地提高计算效率和节省内存。

4 读取文件数据

现在要对程序功能进行扩展。

功能需求：在案例的constant路径下有一个名为volumeRatioDict的字典文件，其中定义了一个名为maxRatio的字典，需要通过程序判断体积比率超过maxRatio的网格数量。

这里涉及到如何从字典文件中获取数据的问题。

4.1 准备字典文件

OpenFOAM的字典文件是一个文本文件，可以复制constant中的其他字典文件并修改其内容。

文件volumeRatioDict的文件内容为：

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      volumeRatioDict;
}
maxRatio              1.2;

4.2 修改程序代码

修改程序代码。

#include "fvCFD.H"
 
// * * * * * * * * * * * * * * * //
 
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
 
    // 从字典文件volumeRatioDict中读取数据
    IOdictionary volumeRatioDict(
        IOobject(
            "volumeRatioDict",      //文件名
            mesh.time().constant(), //文件位于constant路径下
            mesh,                   //注册类对象
            IOobject::MUST_READ,    //必须读取
            IOobject::NO_WRITE      // 文件只读
            ));
    scalar maxRatio(readScalar(volumeRatioDict.lookup("maxRatio")));
    // 定义变量nFail用于存储不满足要求的网格数量
    label nFail = 0;
 
    const label c = mesh.C().size();
    Info << "拥有网格数量为:" << c << endl;
    const labelListList &neighbour = mesh.cellCells();
 
    scalar volumeRatio = 0;
    // 定义临时变量tmpRatio
    scalar tmpRatio = 0;
 
    forAll(neighbour, cellI)
    {
        List n = neighbour[cellI];
        const scalar cellVolume = mesh.V()[cellI];
 
        forAll(n, i)
        {
            label neighbourIndex = n[i];
            scalar neighbourVolume = mesh.V()[neighbourIndex];
            if (neighbourVolume >= cellVolume)
            {
                volumeRatio = neighbourVolume / cellVolume;
                // 如果volumeRatio的值比tmpRatio的值大，则更新tmpRatio的值
                if (volumeRatio > tmpRatio)
                {
                    tmpRatio = volumeRatio;
                }
 
                if (volumeRatio > maxRatio)
                {
                    nFail++;
                }
            }
        }
    }
    // 直接输出最大体积比率
    Info << "最大体积比率为：" << tmpRatio << endl;
    Info << "超过标准的网格数量：" << nFail << endl;
 
    Info << nl;
    runTime.printExecutionTime(Info);
 
    Info << "Endn"
         << endl;
 
    return 0;
}

程序运行结果如下图所示。

程序顺利地实现了目标。

需要注意的是，上面创建的程序只能用在串行模式下使用，若要在并行模式下使用，则需要进行并行代码改造，我们在后面的教程中再来聊并行的问题。

（本文结束）