[转载]OFtutorial05_basicParallelComputing

OpenFOAM中的并行计算基础

原文地址:https://www.cnblogs.com/xubonan/p/15480481.html

代码：

Pout << "Hello from processor " << Pstream::myProcNo() << "! I am working on "
     << mesh.C().size() << " cells" << endl;

 

作用：

如果并行化运行一个算例，整个计算域的网格会按照核数分解。每一个小块在一个单独的核中计算，与串行计算类似，每个核中都会有一个对象的实例，比如：mesh，U，p。这些实例可能有不同的大小，但是也比整个网格小很多。与Info不同，每一个核都能通过’Pout‘输出内容。Info只能在主核内输出（processor0）。

注：

在不同的线程中间交换数据的时候，需要调用OpenMPI的线程

代码：

scalar meshVolume(0.);
forAll(mesh.V(),cellI)
   meshVolume += mesh.V()[cellI];

 

作用：

这会遍历每一个子区域中的cell并累加他们的体积。

代码：

Pout << "Mesh volume on this processor: " << meshVolume << endl;
reduce(meshVolume, sumOp<scalar>());
Info << "Total mesh volume on all processors: " << meshVolume
    // Note how the reudction operation may be done in place without defning
    // a temporary variable, where appropriate.
     << " over " << returnReduce(mesh.C().size(), sumOp<label>()) << " cells" << endl;

 

作用：

把所有核中的值累加在一起。

在reduction阶段，可以有不同的操作。上述代码中用sumOp描述的累加就是其中之一。也有一些其他的有用的操作，比如：minOp, maxOp. 在使用模板的时候必须指定数据的类型，比如这个文件中用了scalar数据类型，在小括号之前要加上<scalar>。

代码：

Pstream::scatter(meshVolume);
Pout << "Mesh volume on this processor is now " << meshVolume << endl;

作用：

通过scatter向每一个核中分发数据

代码：

List<label> nInternalFaces (Pstream::nProcs()), nBoundaries (Pstream::nProcs());
nInternalFaces[Pstream::myProcNo()] = mesh.Cf().size();
nBoundaries[Pstream::myProcNo()] = mesh.boundary().size();

作用：

检测所有的线程里面某些东西的分布。这可以通过一个list实现，其中list中的每一个元素只能被一个核写入。

代码：

Pstream::gatherList(nInternalFaces);
Pstream::gatherList(nBoundaries);
Pstream::scatterList(nInternalFaces);
Pstream::scatterList(nBoundaries);

作用：

• list可以被主线程生成
• 也可以向每个线程分发

代码：

if (Pstream::master())
{
    forAll(nInternalFaces,i)
        Info << "Processor " << i << " has " << nInternalFaces[i]
             << " internal faces and " << nBoundaries[i] << " boundary patches" << endl;
}

作用：

也可以只在主线程内操作（在这个例子中，我们使用Info方法，这个方法只在主线程内进行。）现在生成的list有所有的线程的信息。

代码：

forAll(mesh.boundary(),patchI)
    Pout << "Patch " << patchI << " named " << mesh.boundary()[patchI].name() << endl;

作用：

当网格被分解的时候，线程之间的interface被转化成了patches，这代表每个子区域把线程的边界看成边界条件。

代码：

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;
}

作用：

对于线程的patches，可以检测其类型。这与前文中讲到的如何检测一个patch是不是empty类型相似。

下面是把tutorial 04 调整为并行处理的例子。不同的地方都会有NOTE标识出来。

代码：

#include "createFields.H"

头文件中的代码：

Info << "Reading transportProperties\n" << endl;

IOdictionary transportProperties
(
    IOobject
    (
        "transportProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

dimensionedScalar nu
(
    "nu",
    dimViscosity,
    transportProperties.lookup("nu")
);

Info<< "Reading field p\n" << endl;
volScalarField p
(
    IOobject
    (
        "p",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Reading field U\n" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

作用：

• 在OpenFOAM中，会经常把一大块代码写在一个.H文件中，使得求解器本身更易读。这种做法与标准C++的处理不同，因为头文件一般只涉及声明，不涉及定义。
• createFields 是除了通用头文件setRootCase, createTime, and createMesh之外的一个很常用的头文件，这个头文件一般来说每个solver都不一样。
• 这里我们把所有的处理场数据和输运参数的代码转移到createFields 文件中，并包含在主程序中。

代码：

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>());
scalar f (1.);

作用：

**NOTE: ** 得到全局数值，把有单位的量转化为无单位的量。

代码：

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.);

作用：

向场中赋值。

• sin函数的输入是无单位的数据，所以需要.value()把当前的时间转化为无单位的数值。
• r的单位为长度单位，需要在其中加入一个很小的值以防止除0.
• 最后，计算得到的数值必须符合压力的单位: m^2 * s^-2

**注: ** 此时压力为原来的压力除以密度，不可压流动常用这种处理方式

代码：

p.correctBoundaryConditions();

作用：

**NOTE: **需要更新核边界上的值。如果不更新，求梯度的操作在核边界附近会出问题。

代码：

U = fvc::grad(p)*dimensionedScalar("tmp", dimTime, 1.);
runTime.write();

作用：

通过压力梯度计算速度。

总体代码：

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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"

    // For a case being run in parallel, the domain is decomposed into several
    // processor meshes. Each of them is run in a separate process and holds
    // instances of objects like mesh, U or p just as in a single-threaded (serial)
    // computation. These will have different sizes, of course, as they hold
    // fewer elements than the whole, undecomposed, mesh.
    // Pout is a stream to which each processor can write, unlike Info which only
    // gets used by the head process (processor0)
    Pout << "Hello from processor " << Pstream::myProcNo() << "! I am working on "
         << mesh.C().size() << " cells" << endl;

    // To exchange information between processes, special OpenMPI routines need
    // to be called.

    // This goes over each cell in the subdomain and integrates their volume.
    scalar meshVolume(0.);
    forAll(mesh.V(),cellI)
        meshVolume += mesh.V()[cellI];

    // Add the values from all processes together
    Pout << "Mesh volume on this processor: " << meshVolume << endl;
    reduce(meshVolume, sumOp<scalar>());
    Info << "Total mesh volume on all processors: " << meshVolume
        // Note how the reudction operation may be done in place without defning
        // a temporary variable, where appropriate.
         << " over " << returnReduce(mesh.C().size(), sumOp<label>()) << " cells" << endl;
    // During the reduction stage, different operations may be carried out, summation,
    // described by the sumOp template, being one of them.
    // Other very useful operations are minOp and maxOp.
    // Note how the type
    // of the variable must be added to make an instance of the template, here
    // this is done by adding <scalar> in front of the brackets.
    // Custom reduction operations are easy to implement but need fluency in
    // object-oriented programming in OpenFOAM, so we'll skip this for now.

    // 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;

    // It is often useful to check the distribution of something across all
    // processors. This may be done using a list, with each element of it
    // being written to by only one processor.
    List<label> nInternalFaces (Pstream::nProcs()), nBoundaries (Pstream::nProcs());
    nInternalFaces[Pstream::myProcNo()] = mesh.Cf().size();
    nBoundaries[Pstream::myProcNo()] = mesh.boundary().size();

    // The list may then be gathered on the head node as
    Pstream::gatherList(nInternalFaces);
    Pstream::gatherList(nBoundaries);
    // Scattering a list is also possbile
    Pstream::scatterList(nInternalFaces);
    Pstream::scatterList(nBoundaries);

    // It can also be useful to do things on the head node only
    // (in this case this is meaningless since we are using Info, which already
    // checks this and executes on the head node).
    // Note how the gathered lists hold information for all processors now.
    if (Pstream::master())
    {
        forAll(nInternalFaces,i)
            Info << "Processor " << i << " has " << nInternalFaces[i]
                 << " internal faces and " << nBoundaries[i] << " boundary patches" << endl;
    }

    // As the mesh is decomposed, interfaces between processors are turned
    // into patches, meaning each subdomain sees a processor boundary as a
    // boundary condition.
    forAll(mesh.boundary(),patchI)
        Pout << "Patch " << patchI << " named " << mesh.boundary()[patchI].name() << endl;

    // When looking for processor patches, it is useful to check their type,
    // similarly to how one can check if a patch is of empty type
    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.);

        // NOTE: this is needed to update the values on the processor boundaries.
        // If this is not done, the gradient operator will get confused around the
        // processor patches.
        p.correctBoundaryConditions();

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

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

    return 0;
}

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