OpenFOAM中进程之间的信息交流

应用场景：

在OpenFOAM中利用Gussian插值等插值方法将欧拉网格中流场信息插值到某一点过程中需要寻找该点周围网格信息以及流场信息（物理量），如果该点落在进程交接面附近，那么就需要进程进行信息传递，让计算进程可以看到相邻进程部分网格信息以及流场信息。

解决方案：

以如上图所示矩形计算区域为例，计算区域被划分为4个进程，以processor0为例，要想processor0 获得processor1-3的流场信息，主要解决思路为：
由于其他进程（processor1-3）的流场信息以及网格拓扑信息在当前进程(processor0)是看不到的，所以需要一个二维数组来保存其他进程发送的信息。processor0中的二维数组如下表所示：

processor_Id	cells data
0	/
1	…
2	…
3	…

由于当前进程为processor0,所以不会接受到processor0发送的信息（二维数组第一行大小为0），但是其他相邻的进程（processor1-3）会发送数据给当前进程（processor0）,所以数组内会保存数据。OpenFOAM中进程交换过程就是维护这个二维数组的过程。值得注意的是每一个进程都会有一个这样的数组来保存其他进程发送的信息。

1. 获取当前进程有多少个边界网格，并设置选取网格半径大小（procRMax）

// dertermine the pro cells
forAll(fluidCV4DEM, cellI)
{
    if (b_useful[cellI] > SMALL)
    {
        const labelList &faces = meshCells[cellI];
        forAll(faces, faceI)
        {
            label patchID =
                mesh.boundaryMesh().whichPatch(faces[faceI]);

            if (patchID != -1)
            {
                if (
                    isA<processorPolyPatch>(
                        mesh.boundaryMesh()[patchID]))
                {
                    procCellsSet.insert(cellI);
                }
            }
        }
    }
}
labelList procCells = procCellsSet.toc();
sort(procCells);
scalarField procRMax(procCells.size(), 0.0);
forAll(procRMax, i)
{
    scalar volume = mesh.V()[procCells[i]];
    procRMax[i] = Foam::pow(volume, 1.0 / 3.0) + 3.0 * length;
}

得到进程边界处网格，并标记，processor0中边界cells为蓝色标记网格，效果如上图所示。

2. 将网格中心点坐标发送到邻近进程中（processor1-3发送到processor0）

• 首先需要将每个进程边界网格size发送到其他进程（processor1-3发送到processor0），数组sizes记录了相邻进程边界网格个数

// Send and receive num of cells
// next to processor boundaries
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
    if (procI != Pstream::myProcNo())
    {
        // Parallel data exchange
        {
            OPstream toProc(
                Pstream::blocking,
                procI,
                sizeof(label));

            toProc << procCells.size();
        }
    }
}
labelList sizes(Pstream::nProcs(), 0);
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
    if (procI != Pstream::myProcNo())
    {
        // Parallel data exchange
        {
            IPstream fromProc(
                Pstream::blocking,
                procI,
                sizeof(label));

            fromProc >> sizes[procI];
        }
    }
}

• 有了sizes之后将相邻进程网格centers，以及procRMax发送到当前进程(processor1-3发送到processor0)

// Send and receive centres and radii
vectorField centres(procCells.size(), vector::zero);
forAll(centres, cellI)
{
    centres[cellI] = C[procCells[cellI]];
}

for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
    if (procI != Pstream::myProcNo())
    {
        // Parallel data exchange
        {
            OPstream toProc(
                Pstream::blocking,
                procI,
                centres.size() * sizeof(vector) + procRMax.size() * sizeof(scalar));

            toProc << centres << procRMax;
        }
    }
}

FieldField<Field, vector> ctrs(Pstream::nProcs());
FieldField<Field, scalar> rMax(Pstream::nProcs());

for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
    if (procI != Pstream::myProcNo())
    {
        ctrs.set(
            procI,
            new vectorField(sizes[procI], vector::zero));

        rMax.set(
            procI,
            new scalarField(sizes[procI], 0));
        // Parallel data exchange
        {
            IPstream fromProc(
                Pstream::blocking,
                procI,
                sizes[procI] * sizeof(vector) + sizes[procI] * sizeof(scalar));

            fromProc >> ctrs[procI] >> rMax[procI];
        }
    }
    else
    {
        ctrs.set(procI, new vectorField(0));
        rMax.set(procI, new scalarField(0));
    }
}

从这段代码可以看出只有知道了相邻进程（processor1-3）网格大小sizes才可以通过sizes[procI]大小来接受相邻进程（processor1-3）centers，并将其存储到二维数组ctrs中；

3 寻找本进程内（processor0）距离相邻进程（processor1-3）网格中心procRMax距离内的网格

for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
    labelHashSet procCellSet;

    if (procI != Pstream::myProcNo())
    {
        forAll(ctrs[procI], cellI)
        {
            label nearestCellID =
                findNearestCell(mesh, ctrs[procI][cellI], b_useful);

            if (nearestCellID == -1)
            {
                continue;
                FatalErrorIn(
                    "find_cell_cells")
                    << "Can't find nearest cell."
                    << abort(FatalError);
            }

            scalar R = mag(C[nearestCellID] - ctrs[procI][cellI]);

            if (R < rMax[procI][cellI])
            {
                if (!procCellSet.found(nearestCellID) && b_useful[nearestCellID] > SMALL)
                {
                    procCellSet.insert(nearestCellID);
                }

                labelList tmpCellList;
                findCellCells(
                    mesh,
                    ctrs[procI][cellI],
                    nearestCellID,
                    tmpCellList,
                    b_useful);
                forAll(tmpCellList, cI)
                {
                    scalar r =
                        mag(
                            C[tmpCellList[cI]] - ctrs[procI][cellI]);
                    if (r <= rMax[procI][cellI])
                    {
                        if (!procCellSet.found(tmpCellList[cI]))
                        {
                            procCellSet.insert(tmpCellList[cI]);
                        }
                    }
                }
            }
        }
    }

    procCellCells[procI] = procCellSet.toc();
}

如上图所示，processor0中红色网格即是被其他进程processor1-3需要的网格，现进程processor0中的二维数组procCellCells每一行分别包括了距离processor1-3进程边界距离小于procRMax的所有网格ID。

4. 每个进程需要接受各自的procCellCells（数据为ID）对应的网格中心，以processor0为例，需要接收processor1-3中二维数组procCellCells第一行procCellCells[0]对应的网格中心，然后将数据存储在二维数组procCentres中；具体代码如下：

-以processor0为例，获取每个进程的procCellCells[0]数组的大小：

// Send and receive sizes
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
    if (procI != Pstream::myProcNo())
    {
        // Parallel data exchange
        {
            OPstream toProc(
                Pstream::blocking,
                procI,
                sizeof(label));

            toProc << procCellCells[procI].size();
        }
    }
}

labelList procSizes(Pstream::nProcs(), 0);
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
    if (procI != Pstream::myProcNo())
    {
        // Parallel data exchange
        {
            IPstream fromProc(
                Pstream::blocking,
                procI,
                sizeof(label));

            fromProc >> procSizes[procI];
        }
    }
}

• 知道了sizes的大小便可以将网格中心点坐标以同样的方式发送到对应进程：

// Send cell centres
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
    if (procI != Pstream::myProcNo())
    {
        vectorField centres(C, procCellCells[procI]);

        // Parallel data exchange
        {
            OPstream toProc(
                Pstream::blocking,
                procI,
                centres.size() * sizeof(vector));

            toProc << centres;
        }
    }
}

// Receive cell centres
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
    if (procI != Pstream::myProcNo())
    {
        procCentres.set(
            procI,
            new vectorField(procSizes[procI], vector::zero));

        // Parallel data exchange
        {
            IPstream fromProc(
                Pstream::blocking,
                procI,
                procSizes[procI] * sizeof(vector));

            fromProc >> procCentres[procI];
        }
    }
    else
    {
        procCentres.set(procI, new vectorField(0));
    }
}

至此processor1-3中红色网格内所有网格中心对于processor0都是可见的，其数据存储在procCentres二维数组中，再次需要提醒的是每一个进程都包含一个procCentres数组。

如果你只需要网格的几何拓扑信息，做到这一步就可以了，但是如果你还需要物理场的信息，那就需要在每一个时间步内进行信息的发送与接受，下一步介绍这个过程

5.收发的信息

template <class Type>
void sendAndReceive(
    const Field<Type> &psi,
    const labelListList &procCells,
    const FieldField<Field, vector> &procCenters,
    FieldField<Field, Type> &tprocPsi)
{

    forAll(tprocPsi, procI)
    {
        tprocPsi.set(
            procI,
            new Field<Type>(
                procCenters[procI].size(),
                pTraits<Type>::zero));
    }

    if (Pstream::parRun())
    {
        // Send
        for (label procI = 0; procI < Pstream::nProcs(); procI++)
        {
            if (procI != Pstream::myProcNo())
            {
                // Do not send empty lists
                if (!procCells[procI].empty())
                {
                    Field<Type> curPsi(psi, procCells[procI]);

                    // Parallel data exchange
                    OPstream toProc(
                        Pstream::blocking,
                        procI,
                        curPsi.size() * sizeof(Type));
                    toProc << curPsi;
                }
            }
        }
        // Receive
        for (label procI = 0; procI < Pstream::nProcs(); procI++)
        {
            if (procI != Pstream::myProcNo())
            {
                // Do not receive empty lists
                if (!tprocPsi[procI].empty())
                {
                    // Parallel data exchange
                    IPstream fromProc(
                        Pstream::blocking,
                        procI,
                        tprocPsi[procI].size() * sizeof(Type));

                    fromProc >> tprocPsi[procI];
                }
            }
        }
    }
}

该部分代码与第四步收发网格中心点坐标代码类似，由于可能传输标量，矢量或者张量场的信息，所以将此函数写为模板函数的形式。

6. main函数

由于要与课题组基于C语言的程序兼容，所以生成virtual cells的结果是保存在结构中的，main函数首先定义了该结构：

t_virtual_cells virtual_cells;

然后函数：

virtual_cells_init(fluidCV4DEM, &virtual_cells, length);

作用为步骤1-4过程生成procCenters二维数组。
通过调用：

sendAndReceive(labelFluidCell, procCellCells, procCentres, tprocPsi);

将进程边界周围一定范围内的labelFluidCell场的值传送送到当前进程并保存在二维数组tprocPs中。

    forAll(tprocPsi, procI)
    {
        forAll(tprocPsi[procI], cellI)
        {
            Info << tprocPsi[procI][cellI] << endl;
        }
    }

通过上述语句就可输出与processor0相邻的processor1-3边界上的labelFluidCell流场值。
完整的main函数代码如下所示。

#include <math.h>
#include "fvCFD.H"
#include "mathematicalConstants.H"
#include "demandDrivenData.H"
#include "meshSearch.H"
#include "virtual_cells.h"

int main(int argc, char *argv[])
{

#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include <fstream>
#include <iostream>

    t_virtual_cells virtual_cells;

    volScalarField fluidCV4DEM(
        IOobject(
            "fluidCV4DEM",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE),
        mesh,
        dimensionedScalar("zero", dimless, 0));
    volScalarField labelFluidCell(
        IOobject(
            "labelFluidCell",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE),
        mesh,
        dimensionedScalar("0", dimless, 0));

    scalar length = 0.003;
    virtual_cells_init(fluidCV4DEM, &virtual_cells, length);
    FieldField<Field, scalar> *tprocPsi_ = new FieldField<Field, scalar>(Pstream::nProcs());
    FieldField<Field, scalar> &tprocPsi = *tprocPsi_;
    // field communacate
    labelListList &procCellCells = *(virtual_cells.procCellsPtr_);
    FieldField<Field, vector> &procCentres = *(virtual_cells.procCentresPtr_);
    forAll(labelFluidCell, cellI)
    {
        labelFluidCell[cellI] = Pstream::myProcNo();
    }

    sendAndReceive(labelFluidCell, procCellCells, procCentres, tprocPsi);

    forAll(procCentres, procI)
    {
        forAll(procCentres[procI], cellI)
        {
            Info << procCentres[procI][cellI] << endl;
        }
    }
    forAll(tprocPsi, procI)
    {
        forAll(tprocPsi[procI], cellI)
        {
            Info << tprocPsi[procI][cellI] << endl;
        }
    }

    fluidCV4DEM.write();
    labelFluidCell.write();
    Info << "the virtual length is " << length << endl;
    Info << "End!" << endl;
    return 0;
}

进程之间信息交互本身比较抽象，虽然我进行了图文总结，但是不一定一下子就能看明白，如果需要将这段代码使用到自己程序中可能需要看很多遍才能理解，希望我的解读能起到抛砖引玉的效果，完整的代码献上，希望大家能玩出更多花样。