OpenFOAM中的并行计算

OpenFOAM并行计算的特点

(1) 单指令多数据(同一套代码，不同网格数据)

(2) 相对独立地运行n个进程，进程间的通信由OpenMPI实现。如果不调用MPI函数
（reduce，returnReduce，scatter，gatherList，scatterList）那么这些进程始终是独立的。

(3) 网格区块之间使用特殊的边界条件

代码示例

可参见OpenFOAM教程https://gitee.com/fr13ndsdp/OpenFOAM_TUTORIAL中的OFtutorial05_basicParallelComputing

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"

    //对于OpenFOAM并行程序，求解区域被分割为若干各子区域。每个进程会分配独立地一个区域，拥有自己
    //的变量。不同进程的网格是不同的，所以U或p场也是不同的。

    // Pout 是一个输出流，每一个进程都可以输出内容
    // Info 也是一个输出流，只由头进程输出(processor0)
    Pout << "Hello from processor " << Pstream::myProcNo() << "! I am working on "
         << mesh.C().size() << " cells" << endl;

    // 需要调用openMPI实现进程间的通信

    // 在当前进程中对单元体积求和，存入meshVolume
    scalar meshVolume(0.);
    forAll(mesh.V(),cellI)
        meshVolume += mesh.V()[cellI];

    // 使用reduce函数将各进程的meshVolume变量汇总求和
    //reduce函数中的运算还可有minOp和maxOp等
    Pout << "Mesh volume on this processor: " << meshVolume << endl;
    reduce(meshVolume, sumOp<scalar>());
    Info << "Total mesh volume on all processors: " << meshVolume
        //可以使用returnReduce代替reduce函数，直接返回结果而不改变原始变量
         << " over " << returnReduce(mesh.C().size(), sumOp<label>()) << " cells" << endl;


    // Spreading a value across all processors is done using a scatter operation.
    Pstream::scatter(meshVolume);
    Pout << "Mesh volume on this processor is now " << meshVolume << endl;

    // 通常可以定义一个列表（数组），长度为Pstream::nProcs()，每个元素储存该进程对应的数据
    // Pstream::myProcNo()给出当前进程的编号
    List<label> nInternalFaces (Pstream::nProcs()), nBoundaries (Pstream::nProcs());
    nInternalFaces[Pstream::myProcNo()] = mesh.Cf().size();
    nBoundaries[Pstream::myProcNo()] = mesh.boundary().size();
    // 列表可以汇总到头进程
    Pstream::gatherList(nInternalFaces);
    Pstream::gatherList(nBoundaries);
    // 列表也可以分散到所有进程
    Pstream::scatterList(nInternalFaces);
    Pstream::scatterList(nBoundaries);

    // 使用Pstream::master()仅在头进程中执行代码
    if (Pstream::master())
    {
        forAll(nInternalFaces,i)
            Pout << "Processor " << i << " has " << nInternalFaces[i]
                 << " internal faces and " << nBoundaries[i] << " boundary patches" << endl;
    }

    // 网格被划分为若干块之后，块与块之间则会形成新的边界
    forAll(mesh.boundary(),patchI)
        Pout << "Patch " << patchI << " named " << mesh.boundary()[patchI].name() << endl;

    // 检查边界是否是processor boundary
    forAll(mesh.boundary(),patchI)
    {
        const polyPatch& pp = mesh.boundaryMesh()[patchI];
        if (isA<processorPolyPatch>(pp))
            Pout << "Patch " << patchI << " named " << mesh.boundary()[patchI].name()
                 << " is definitely a processor boundary!" << endl;
    }

    // ---
    // this is an example implementation of the code from tutoral 2 which
    // has been adjusted to run in parallel. Each difference is highlighted
    // as a NOTE.

    // It is conventional in OpenFOAM to move large parts of code to separate
    // .H files to make the code of the solver itself more readable. This is not
    // a standard C++ practice, as header files are normally associated with
    // declarations rather than definitions.
    // A very common include, apart from the setRootCase, createTime, and createMesh,
    // which are generic, is createFields, which is often unique for each solver.
    // Here we've moved all of the parts of the code dealing with setting up the fields
    // and transport constants into this include file.
    #include "createFields.H"

    // pre-calculate geometric information using field expressions rather than
    // cell-by-cell assignment.
	const dimensionedVector originVector("x0", dimLength, vector(0.05,0.05,0.005));
    volScalarField r (mag(mesh.C()-originVector));
    // NOTE: we need to get a global value; convert from dimensionedScalar to scalar
	const scalar rFarCell = returnReduce(max(r).value(), maxOp<scalar>());
    scalar f (1.);

	Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        // assign values to the field;
        // sin function expects a dimensionless argument, hence need to convert
        // current time using .value().
        // r has dimensions of length, hence the small value being added to it
        // needs to match that.
        // Finally, the result has to match dimensions of pressure, which are
        // m^2 / s^-2/
		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.);

        // 注意：correctBoundaryConditions函数可完成网格块之间的数据交换
        // 如果没有这一行代码，界面处就会出错
        p.correctBoundaryConditions();

        // calculate velocity from gradient of pressure
		U = fvc::grad(p)*dimensionedScalar("tmp", dimTime, 1.);
		runTime.write();
	}

    Info<< "End\n" << endl;

    return 0;
}

// ************************************************************************* //