本案例演示利用OpenFOAM进行并行计算的基础方法。

CFD计算自然离不开并行计算，在对大规模问题进行计算时，OpenFOAM也可以使用并行模式。然而利用并行模式进行计算时，可能需要修改一些程序代码，尤其是在需要获取全场信息时。本文以一个简单的案例演示并行程序编写时需要注意的一些问题。

1 文件准备

利用下面的命令准备文件：

run
foamNewApp demo5
cd demo5 && touch createFields.H

如下图所示。

2 程序代码

• 头文件createFields.H

// 定义一些物理量
 
Info << "读取transportProperties文件" << endl;
IOdictionary transportProperties
(
    IOobject
    (
        "transportProperties", //文件名
        runTime.constant(),     //文件路径
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);
// 读取物理属性nu
dimensionedScalar nu
(
    "nu",
    dimViscosity,
    transportProperties
);
 
 
Info << "读取速度U" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);
 
Info << "读取压力场p" << endl;
volScalarField p
(
    IOobject
    (
        "p",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);
 

• 源文件demo5.C

#include "fvCFD.H"
// * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
    // 检查案例文件组织结构
    #include "setRootCase.H"
    // 创建runTime对象
    #include "createTime.H"
    // 创建fvMesh对象
    #include "createMesh.H"
 
    /*
    对于并行运行的案例，计算域被分解成若干个processor网格。
    每个子域都在单独的处理器中运行，并保存网格、U或p等对象的实例，就像串行计算一样。
    Pout是每个处理器都可以写入的流，而Info只能由head进程(Processor0)使用。
    在并行计算时，若想要输出计算节点的信息，需要使用Pout流
    */
    Pout << "这是处理器No：" << Pstream::myProcNo() << ",网格数量：" << mesh.C().size() << endl;
 
    // 在处理器之间交换信息，需要调用特殊的OpenMPI进程
    // 计算子域的体积
    scalar meshVolume;
    forAll(mesh.V(), cellI)
    {
        meshVolume += mesh.V()[cellI];
    }
    Pout << "本处理器上的网格体积：" << meshVolume << endl;
 
    // 计算所有的子域的网格体积，使用全局约简函数reduce
    reduce(meshVolume, sumOp());
    // 利用Info在head节点上输出信息
    Info << "所有处理器上的总网格数量：" << returnReduce(mesh.C().size(), sumOp()) << endl
         << "网格总体积：" << meshVolume << endl;
 
    /*
    在并行约简时可以执行各种操作，如前面使用的sumOp模板所描述的求和是其中之一。
    其他非常有用的操作还有minOp和maxOp，分别用于提取最小值与最大值。
    注意必须添加变量类型才能创建模板实例，这里是通过在尖括号中指定变量类型来实现的。
    自定义约简操作很容易实现，但是需要熟练地使用OpenFOAM中的面向对象编程，所以现在先跳过这一步。
    */
    // 在所有处理器发送一个值可以通过scatter函数来实现
    Pstream ::scatter(meshVolume);
    Pout << "现在本处理器上的网格体积为：" << meshVolume << endl;
 
    /*
    检查某些内容在所有处理器上的分布通常很有用。
    可以使用List来完成，列表中的每个元素只由一个处理器写入。
    */
    List nInternalFaces(Pstream::nProcs());
    List nBoundaries(Pstream::nProcs());
    nInternalFaces[Pstream::myProcNo()] = mesh.Cf().size();
    nBoundaries[Pstream::myProcNo()] = mesh.boundary().size();
 
    // 在head节点收集信息
    Pstream::gatherList(nInternalFaces);
    Pstream::gatherList(nBoundaries);
    // 也可以散布数据
    Pstream::scatterList(nInternalFaces);
    Pstream::scatterList(nBoundaries);
 
    /*
    有时候仅在头节点上执行操作也很有用
    */
    if (Pstream::master())
    {
        forAll(nInternalFaces, i)
        {
            // 这里其实可以直接使用Info，不需要进行master判断
            Pout << "处理器[" << i << "]拥有 " << nInternalFaces[i]
                << " 个内部面，以及 " << nBoundaries[i] << " 个边界面" << endl;
        }
    }
 
    /*
    由于区域分解，导致处理器之间的交界面变为了patch，意味着每个子域都将processor边界视为边界
    */
    forAll(mesh.boundary(), patchI)
    {
        Pout << "Patch " << patchI << " 的名称为 " << mesh.boundary()[patchI].name() << endl;
    }
 
    /*
    在搜索processor边界时，可以先判断其类型
    */
    forAll(mesh.boundary(), patchI)
    {
        const polyPatch &pp = mesh.boundaryMesh()[patchI];
        if (isA(pp))
        {
            Pout << "Patch " << patchI << " 的名称为：" << mesh.boundary()[patchI].name()
                << " ，其定义在processor边界上." << endl;
        }
    }
 
    /*
    这是Tutorial 2中的代码的一个示例，该代码已调整为并行运行。
    在OpenFOAM中，通常会将大部分声明代码放到单独的.h文件中，以使求解器本身的代码更具可读性。
    除了setRootCase、createTime和createMesh(它们是通用的)之外，
    一个非常常见的头文件是createFields，其对于每个求解器通常是唯一的。
    这里已经将处理设置物理场和输运常量的所有代码放到了这个头文件中。
    */
    #include "createFields.H"
    const dimensionedVector originVector("x0", dimLength, vector(0.05, 0.05, 0.005));
    volScalarField r(mag(mesh.C() - originVector));
 
    // 计算网格中心到指定点的最大距离，全局参数需要改造成并行代码
    const scalar rFarCell = returnReduce(max(r).value(), maxOp());
 
    scalar f(1.0);
    Info << "n开始时间迭代" << endl;
 
    while (runTime.loop())
    {
        // 计算压力，确保量纲一致
        p = Foam::sin(2. * constant::mathematical::pi * f *
                      runTime.time().value()) / (r / rFarCell + dimensionedScalar("small", dimLength, 1e-12)) *
            dimensionedScalar("tmp", dimensionSet(0, 3, -2, 0, 0), 1.);
 
        // 并行计算必须修正边界
        p.correctBoundaryConditions();
 
        U = fvc::grad(p) * dimensionedScalar("tmp", dimTime, 1.0);
        runTime.write();
    }
    // * * * * * * * * * *  * * * * * * * * * * //
    Foam::Info << nl;
    runTime.printExecutionTime(Info);
    Info << "Endn" << endl;
    return 0;
}

利用命令wmake编译程序，确保没有错误或警告信息，如下图所示。

3 测试程序

这里使用cavity的案例，并行计算需要添加decomposeParDict字典。进入案例根路径，运行下面的命令添加字典文件：

foamGetDict decomposeParDict

利用decomposeParDict字典文件指定并行分区方式。这里使用的字典文件如下：

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      decomposeParDict;
}
// * * * * * * * * * * * * * * * * * //
 
numberOfSubdomains 4;
method hierarchical;
hierarchicalCoeffs
{
    n               (4 1 1);
    delta 0.0001;
    order           xyz;
}

进入案例根路径，运行下面的命令运行案例：

blockMesh
decomposePar
mpirun -np 4 demo5 -parallel

运行结果如下图所示，可以看到并行模式获取了正确的信息。

（完毕）