ENH: replace listCombineReduce with equivalent

- a few places still used listCombineReduce instead of
  newer constructs (eg allGatherList) with fewer MPI calls.

- align triangulation handling of turbulentDFSEMInlet and
  patchInjectionBase with meshTools/triangulatedPatch
  (will ease future code refactoring)
This commit is contained in:
Mark Olesen 2024-02-05 09:55:41 +01:00
parent d89ecc74be
commit f6825c7952
7 changed files with 188 additions and 204 deletions

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2016-2022 OpenCFD Ltd.
Copyright (C) 2016-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -132,58 +132,76 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::initialisePatch()
patchNormal_ /= mag(patchNormal_) + ROOTVSMALL;
// Decompose the patch faces into triangles to enable point search
const polyPatch& patch = this->patch().patch();
const pointField& points = patch.points();
// Triangulate the patch faces and create addressing
DynamicList<label> triToFace(2*patch.size());
DynamicList<scalar> triMagSf(2*patch.size());
DynamicList<face> triFace(2*patch.size());
DynamicList<face> tris(5);
// Set zero value at the start of the tri area list
triMagSf.append(0.0);
forAll(patch, faceI)
{
const face& f = patch[faceI];
tris.clear();
f.triangles(points, tris);
forAll(tris, i)
label nTris = 0;
for (const face& f : patch)
{
triToFace.append(faceI);
triFace.append(tris[i]);
triMagSf.append(tris[i].mag(points));
nTris += f.nTriangles();
}
DynamicList<labelledTri> dynTriFace(nTris);
DynamicList<face> tris(8); // work array
forAll(patch, facei)
{
const face& f = patch[facei];
tris.clear();
f.triangles(points, tris);
for (const auto& t : tris)
{
dynTriFace.emplace_back(t[0], t[1], t[2], facei);
}
}
// Transfer to persistent storage
triFace_.transfer(dynTriFace);
}
const label myProci = UPstream::myProcNo();
const label numProc = UPstream::nProcs();
// Calculate the cumulative triangle weights
{
triCumulativeMagSf_.resize_nocopy(triFace_.size()+1);
auto iter = triCumulativeMagSf_.begin();
// Set zero value at the start of the tri area/weight list
scalar patchArea = 0;
*iter++ = patchArea;
// Calculate cumulative and total area
for (const auto& t : triFace_)
{
patchArea += t.mag(points);
*iter++ = patchArea;
}
sumTriMagSf_.resize_nocopy(numProc+1);
sumTriMagSf_[0] = 0;
{
scalarList::subList slice(sumTriMagSf_, numProc, 1);
slice[myProci] = patchArea;
Pstream::allGatherList(slice);
}
// Convert to cumulative sum of areas per proc
for (label i = 1; i < sumTriMagSf_.size(); ++i)
{
sumTriMagSf_[i] += sumTriMagSf_[i-1];
}
}
sumTriMagSf_ = Zero;
sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
Pstream::listCombineReduce(sumTriMagSf_, maxEqOp<scalar>());
for (label i = 1; i < triMagSf.size(); ++i)
{
triMagSf[i] += triMagSf[i-1];
}
// Transfer to persistent storage
triFace_.transfer(triFace);
triToFace_.transfer(triToFace);
triCumulativeMagSf_.transfer(triMagSf);
// Convert sumTriMagSf_ into cumulative sum of areas per proc
for (label i = 1; i < sumTriMagSf_.size(); ++i)
{
sumTriMagSf_[i] += sumTriMagSf_[i-1];
}
// Global patch area (over all processors)
patchArea_ = sumTriMagSf_.last();
patchArea_ = sumTriMagSf_.back();
// Local patch bounds (this processor)
patchBounds_ = boundBox(patch.localPoints(), false);
@ -239,33 +257,32 @@ Foam::pointIndexHit Foam::turbulentDFSEMInletFvPatchVectorField::setNewPosition
const bool global
)
{
// Initialise to miss
pointIndexHit pos(false, vector::max, -1);
const polyPatch& patch = this->patch().patch();
const pointField& points = patch.points();
label triI = -1;
if (global)
{
const scalar areaFraction =
rndGen_.globalPosition<scalar>(0, patchArea_);
// Determine which processor to use
label procI = 0;
label proci = 0;
forAllReverse(sumTriMagSf_, i)
{
if (areaFraction >= sumTriMagSf_[i])
{
procI = i;
proci = i;
break;
}
}
if (Pstream::myProcNo() == procI)
if (UPstream::myProcNo() == proci)
{
// Find corresponding decomposed face triangle
label triI = 0;
const scalar offset = sumTriMagSf_[procI];
triI = 0;
const scalar offset = sumTriMagSf_[proci];
forAllReverse(triCumulativeMagSf_, i)
{
if (areaFraction > triCumulativeMagSf_[i] + offset)
@ -274,20 +291,13 @@ Foam::pointIndexHit Foam::turbulentDFSEMInletFvPatchVectorField::setNewPosition
break;
}
}
// Find random point in triangle
const face& tf = triFace_[triI];
const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
pos.hitPoint(tri.randomPoint(rndGen_));
pos.setIndex(triToFace_[triI]);
}
}
else
{
// Find corresponding decomposed face triangle on local processor
label triI = 0;
const scalar maxAreaLimit = triCumulativeMagSf_.last();
triI = 0;
const scalar maxAreaLimit = triCumulativeMagSf_.back();
const scalar areaFraction = rndGen_.position<scalar>(0, maxAreaLimit);
forAllReverse(triCumulativeMagSf_, i)
@ -298,16 +308,22 @@ Foam::pointIndexHit Foam::turbulentDFSEMInletFvPatchVectorField::setNewPosition
break;
}
}
// Find random point in triangle
const face& tf = triFace_[triI];
const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
pos.hitPoint(tri.randomPoint(rndGen_));
pos.setIndex(triToFace_[triI]);
}
return pos;
if (triI >= 0)
{
return pointIndexHit
(
true,
// Find random point in triangle
triFace_[triI].tri(points).randomPoint(rndGen_),
triFace_[triI].index()
);
}
// No hit
return pointIndexHit(false, vector::max, -1);
}
@ -573,7 +589,6 @@ turbulentDFSEMInletFvPatchVectorField
patchArea_(-1),
triFace_(),
triToFace_(),
triCumulativeMagSf_(),
sumTriMagSf_(Pstream::nProcs() + 1, Zero),
patchNormal_(Zero),
@ -615,7 +630,6 @@ turbulentDFSEMInletFvPatchVectorField
patchArea_(ptf.patchArea_),
triFace_(ptf.triFace_),
triToFace_(ptf.triToFace_),
triCumulativeMagSf_(ptf.triCumulativeMagSf_),
sumTriMagSf_(ptf.sumTriMagSf_),
patchNormal_(ptf.patchNormal_),
@ -656,7 +670,6 @@ turbulentDFSEMInletFvPatchVectorField
patchArea_(-1),
triFace_(),
triToFace_(),
triCumulativeMagSf_(),
sumTriMagSf_(Pstream::nProcs() + 1, Zero),
patchNormal_(Zero),
@ -684,10 +697,11 @@ turbulentDFSEMInletFvPatchVectorField
Foam::turbulentDFSEMInletFvPatchVectorField::
turbulentDFSEMInletFvPatchVectorField
(
const turbulentDFSEMInletFvPatchVectorField& ptf
const turbulentDFSEMInletFvPatchVectorField& ptf,
const DimensionedField<vector, volMesh>& iF
)
:
fixedValueFvPatchField<vector>(ptf),
fixedValueFvPatchField<vector>(ptf, iF),
U_(ptf.U_.clone(patch().patch())),
R_(ptf.R_.clone(patch().patch())),
L_(ptf.L_.clone(patch().patch())),
@ -702,7 +716,6 @@ turbulentDFSEMInletFvPatchVectorField
patchArea_(ptf.patchArea_),
triFace_(ptf.triFace_),
triToFace_(ptf.triToFace_),
triCumulativeMagSf_(ptf.triCumulativeMagSf_),
sumTriMagSf_(ptf.sumTriMagSf_),
patchNormal_(ptf.patchNormal_),
@ -723,40 +736,10 @@ turbulentDFSEMInletFvPatchVectorField
Foam::turbulentDFSEMInletFvPatchVectorField::
turbulentDFSEMInletFvPatchVectorField
(
const turbulentDFSEMInletFvPatchVectorField& ptf,
const DimensionedField<vector, volMesh>& iF
const turbulentDFSEMInletFvPatchVectorField& ptf
)
:
fixedValueFvPatchField<vector>(ptf, iF),
U_(ptf.U_.clone(patch().patch())),
R_(ptf.R_.clone(patch().patch())),
L_(ptf.L_.clone(patch().patch())),
delta_(ptf.delta_),
d_(ptf.d_),
kappa_(ptf.kappa_),
Uref_(ptf.Uref_),
Lref_(ptf.Lref_),
scale_(ptf.scale_),
m_(ptf.m_),
nCellPerEddy_(ptf.nCellPerEddy_),
patchArea_(ptf.patchArea_),
triFace_(ptf.triFace_),
triToFace_(ptf.triToFace_),
triCumulativeMagSf_(ptf.triCumulativeMagSf_),
sumTriMagSf_(ptf.sumTriMagSf_),
patchNormal_(ptf.patchNormal_),
patchBounds_(ptf.patchBounds_),
eddies_(ptf.eddies_),
v0_(ptf.v0_),
rndGen_(ptf.rndGen_),
sigmax_(ptf.sigmax_),
maxSigmaX_(ptf.maxSigmaX_),
nEddy_(ptf.nEddy_),
curTimeIndex_(-1),
singleProc_(ptf.singleProc_),
writeEddies_(ptf.writeEddies_)
turbulentDFSEMInletFvPatchVectorField(ptf, ptf.internalField())
{}

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2016-2021 OpenCFD Ltd.
Copyright (C) 2016-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -146,6 +146,7 @@ SourceFiles
#include "eddy.H"
#include "pointIndexHit.H"
#include "PatchFunction1.H"
#include "labelledTri.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -204,13 +205,11 @@ class turbulentDFSEMInletFvPatchVectorField
//- Patch area - total across all processors
scalar patchArea_;
//- Decomposed patch faces as a list of triangles
faceList triFace_;
//- The polyPatch faces as triangles, the index of each corresponds
//- to the undecomposed patch face index.
List<labelledTri> triFace_;
//- Addressing from per triangle to patch face
labelList triToFace_;
//- Cumulative triangle area per triangle face
//- Cumulative triangle area per triangle face (processor-local)
scalarList triCumulativeMagSf_;
//- Cumulative area fractions per processor

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017, 2020 OpenFOAM Foundation
Copyright (C) 2017-2023 OpenCFD Ltd.
Copyright (C) 2017-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -101,11 +101,9 @@ void Foam::Cloud<ParticleType>::writeCloudUniformProperties() const
)
);
labelList np(Pstream::nProcs(), Zero);
np[Pstream::myProcNo()] = ParticleType::particleCount_;
// FIXME: replace with Pstream::allGatherList(np);
Pstream::listCombineReduce(np, maxEqOp<label>());
labelList np(UPstream::nProcs(), Foam::zero{});
np[UPstream::myProcNo()] = ParticleType::particleCount_;
Pstream::allGatherList(np);
uniformPropsDict.add
(

View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2013-2017 OpenFOAM Foundation
Copyright (C) 2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -43,13 +44,12 @@ Foam::patchInjectionBase::patchInjectionBase
:
patchName_(patchName),
patchId_(mesh.boundaryMesh().findPatchID(patchName_)),
patchArea_(0.0),
patchArea_(0),
patchNormal_(),
cellOwners_(),
triFace_(),
triToFace_(),
triCumulativeMagSf_(),
sumTriMagSf_(Pstream::nProcs() + 1, Zero)
sumTriMagSf_()
{
if (patchId_ < 0)
{
@ -71,7 +71,6 @@ Foam::patchInjectionBase::patchInjectionBase(const patchInjectionBase& pib)
patchNormal_(pib.patchNormal_),
cellOwners_(pib.cellOwners_),
triFace_(pib.triFace_),
triToFace_(pib.triToFace_),
triCumulativeMagSf_(pib.triCumulativeMagSf_),
sumTriMagSf_(pib.sumTriMagSf_)
{}
@ -88,48 +87,68 @@ void Foam::patchInjectionBase::updateMesh(const polyMesh& mesh)
cellOwners_ = patch.faceCells();
// Triangulate the patch faces and create addressing
DynamicList<label> triToFace(2*patch.size());
DynamicList<scalar> triMagSf(2*patch.size());
DynamicList<face> triFace(2*patch.size());
DynamicList<face> tris(5);
// Set zero value at the start of the tri area list
triMagSf.append(0.0);
forAll(patch, facei)
{
const face& f = patch[facei];
tris.clear();
f.triangles(points, tris);
forAll(tris, i)
label nTris = 0;
for (const face& f : patch)
{
triToFace.append(facei);
triFace.append(tris[i]);
triMagSf.append(tris[i].mag(points));
nTris += f.nTriangles();
}
DynamicList<labelledTri> dynTriFace(nTris);
DynamicList<face> tris(8); // work array
forAll(patch, facei)
{
const face& f = patch[facei];
tris.clear();
f.triangles(points, tris);
for (const auto& t : tris)
{
dynTriFace.emplace_back(t[0], t[1], t[2], facei);
}
}
// Transfer to persistent storage
triFace_.transfer(dynTriFace);
}
sumTriMagSf_ = Zero;
sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
Pstream::listCombineReduce(sumTriMagSf_, maxEqOp<scalar>());
const label myProci = UPstream::myProcNo();
const label numProc = UPstream::nProcs();
for (label i = 1; i < triMagSf.size(); i++)
// Calculate the cumulative triangle weights
{
triMagSf[i] += triMagSf[i-1];
}
triCumulativeMagSf_.resize_nocopy(triFace_.size()+1);
// Transfer to persistent storage
triFace_.transfer(triFace);
triToFace_.transfer(triToFace);
triCumulativeMagSf_.transfer(triMagSf);
auto iter = triCumulativeMagSf_.begin();
// Convert sumTriMagSf_ into cumulative sum of areas per proc
for (label i = 1; i < sumTriMagSf_.size(); i++)
{
sumTriMagSf_[i] += sumTriMagSf_[i-1];
// Set zero value at the start of the tri area/weight list
scalar patchArea = 0;
*iter++ = patchArea;
// Calculate cumulative and total area
for (const auto& t : triFace_)
{
patchArea += t.mag(points);
*iter++ = patchArea;
}
sumTriMagSf_.resize_nocopy(numProc+1);
sumTriMagSf_[0] = 0;
{
scalarList::subList slice(sumTriMagSf_, numProc, 1);
slice[myProci] = patchArea;
Pstream::allGatherList(slice);
}
// Convert to cumulative
for (label i = 1; i < sumTriMagSf_.size(); ++i)
{
sumTriMagSf_[i] += sumTriMagSf_[i-1];
}
}
const scalarField magSf(mag(patch.faceAreas()));
@ -152,12 +171,12 @@ Foam::label Foam::patchInjectionBase::setPositionAndCell
{
label facei = -1;
if (cellOwners_.size() > 0)
if (!cellOwners_.empty())
{
// Determine which processor to inject from
const label proci = whichProc(fraction01);
if (Pstream::myProcNo() == proci)
if (UPstream::myProcNo() == proci)
{
const scalar areaFraction = fraction01*patchArea_;
@ -174,15 +193,13 @@ Foam::label Foam::patchInjectionBase::setPositionAndCell
}
// Set cellOwner
facei = triToFace_[trii];
facei = triFace_[trii].index();
cellOwner = cellOwners_[facei];
// Find random point in triangle
const polyPatch& patch = mesh.boundaryMesh()[patchId_];
const pointField& points = patch.points();
const face& tf = triFace_[trii];
const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
const point pf(tri.randomPoint(rnd));
const point pf = triFace_[trii].tri(points).randomPoint(rnd);
// Position perturbed away from face (into domain)
const scalar a = rnd.position(scalar(0.1), scalar(0.5));
@ -239,17 +256,9 @@ Foam::label Foam::patchInjectionBase::setPositionAndCell
tetPti = cellTetIs[teti].tetPt();
}
}
else
{
cellOwner = -1;
tetFacei = -1;
tetPti = -1;
// Dummy position
position = pTraits<vector>::max;
}
}
else
if (facei == -1)
{
cellOwner = -1;
tetFacei = -1;

View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2013-2016 OpenFOAM Foundation
Copyright (C) 2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -39,21 +40,22 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef patchInjectionBase_H
#define patchInjectionBase_H
#ifndef Foam_patchInjectionBase_H
#define Foam_patchInjectionBase_H
#include "word.H"
#include "labelList.H"
#include "scalarList.H"
#include "vectorList.H"
#include "faceList.H"
#include "labelledTri.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declarations
// Forward Declarations
class polyMesh;
class fvMesh;
class Random;
@ -66,7 +68,7 @@ class patchInjectionBase
{
protected:
// Protected data
// Protected Data
//- Patch name
const word patchName_;
@ -83,13 +85,11 @@ protected:
//- List of cell labels corresponding to injector positions
labelList cellOwners_;
//- Decomposed patch faces as a list of triangles
faceList triFace_;
//- The polyPatch faces as triangles, the index of each corresponds
//- to the undecomposed patch face index.
List<labelledTri> triFace_;
//- Addressing from per triangle to patch face
labelList triToFace_;
//- Cumulative triangle area per triangle face
//- Cumulative triangle area per triangle face (processor-local)
scalarList triCumulativeMagSf_;
//- Cumulative area fractions per processor

View File

@ -171,31 +171,32 @@ void Foam::triangulatedPatch::update()
// Set zero value at the start of the tri area/weight list
scalar patchArea = 0;
*iter = patchArea;
++iter;
*iter++ = patchArea;
// Calculate cumulative and total area (processor-local at this point)
for (const auto& t : triFace_)
{
patchArea += t.mag(points);
*iter = patchArea;
++iter;
*iter++ = patchArea;
}
// FIXME: use allGatherList of subslice
scalarList procSumWght(numProc+1, Foam::zero{});
procSumWght[myProci+1] = patchArea;
Pstream::listCombineReduce(procSumWght, maxEqOp<scalar>());
scalarList procSumArea(numProc+1);
procSumArea[0] = 0;
{
scalarList::subList slice(procSumArea, numProc, 1);
slice[myProci] = patchArea;
Pstream::allGatherList(slice);
}
// Convert to cumulative
for (label i = 1; i < procSumWght.size(); ++i)
for (label i = 1; i < procSumArea.size(); ++i)
{
procSumWght[i] += procSumWght[i-1];
procSumArea[i] += procSumArea[i-1];
}
const scalar offset = procSumWght[myProci];
const scalar totalArea = procSumWght.back();
const scalar offset = procSumArea[myProci];
const scalar totalArea = procSumArea.back();
// Apply processor offset and normalise - for a global 0-1 interval
for (scalar& w : triWght_)

View File

@ -669,19 +669,13 @@ void Foam::topOVariablesBase::writeFluidSolidInterface
// invertOneToMany(surfFaces_.size(), changedFaceToCutFace);
// Transform origin cut faces to a global numbering
labelList cuttingFacesPerProc(Pstream::nProcs(), Zero);
cuttingFacesPerProc[Pstream::myProcNo()] = changedFaces.size();
Pstream::listCombineReduce(cuttingFacesPerProc, plusEqOp<label>());
labelList passedFaces(Pstream::nProcs(), Zero);
for (label i = 1; i < Pstream::nProcs(); ++i)
{
passedFaces[i] = passedFaces[i - 1] + cuttingFacesPerProc[i - 1];
}
const label changedFacesOffset =
globalIndex::calcOffset(changedFaces.size());
forAll(changedFacesInfo, facei)
for (auto& info : changedFacesInfo)
{
changedFacesInfo[facei].data() += passedFaces[Pstream::myProcNo()];
info.data() += changedFacesOffset;
}
}