ENH: improve handling of finiteArea mesh with distributed roots

- in makeFaMesh, the serial fields are now only read on the master
  process and broadcast to the other ranks. The read+distribute is
  almost identical to that used in redistributePar, except that in
  this case entire fields are sent and not a zero-sized subset.

- improved internal faMesh checking for files so that the TryNew
  method works with distributed roots.
This commit is contained in:
Mark Olesen
2022-11-20 15:16:31 +01:00
parent 21e7ce8f42
commit 013f3cccc4
17 changed files with 440 additions and 243 deletions

View File

@ -38,6 +38,7 @@ SourceFiles
#define Foam_fieldsDistributor_H
#include "IOobjectList.H"
#include "bitSet.H"
#include "boolList.H"
#include "PtrList.H"
#include "GeometricField.H"
@ -53,6 +54,22 @@ namespace Foam
class fieldsDistributor
{
// Private Methods
//- Read volume/surface/point/area fields that may or may not exist
//- on all processors
template<class BoolListType, class GeoField, class MeshSubsetter>
static void readFieldsImpl
(
const BoolListType& haveMeshOnProc,
const MeshSubsetter* subsetter,
const typename GeoField::Mesh& mesh,
IOobjectList& allObjects,
PtrList<GeoField>& fields,
const bool deregister
);
public:
// Reading helpers
@ -103,6 +120,32 @@ public:
);
//- Read volume/surface/point/area fields that may or may not exist
//- on all processors
template<class GeoField, class MeshSubsetter>
static void readFields
(
const bitSet& haveMeshOnProc,
const MeshSubsetter* subsetter,
const typename GeoField::Mesh& mesh,
IOobjectList& allObjects,
PtrList<GeoField>& fields,
const bool deregister = false
);
//- Read volume/surface/point/area fields that may or may not exist
//- on all processors
template<class GeoField, class MeshSubsetter>
static void readFields
(
const boolList& haveMeshOnProc,
const MeshSubsetter* subsetter,
const typename GeoField::Mesh& mesh,
IOobjectList& allObjects,
PtrList<GeoField>& fields,
const bool deregister = false
);
//- Read volume/surface/point/area fields that may or may not exist
//- on all processors
template<class GeoField, class MeshSubsetter>
@ -110,7 +153,7 @@ public:
(
const boolList& haveMeshOnProc,
const typename GeoField::Mesh& mesh,
const autoPtr<MeshSubsetter>& subsetterPtr,
const autoPtr<MeshSubsetter>& subsetter,
IOobjectList& allObjects,
PtrList<GeoField>& fields,
const bool deregister = false

View File

@ -101,12 +101,12 @@ void Foam::fieldsDistributor::readFields
}
template<class GeoField, class MeshSubsetter>
void Foam::fieldsDistributor::readFields
template<class BoolListType, class GeoField, class MeshSubsetter>
void Foam::fieldsDistributor::readFieldsImpl
(
const boolList& haveMeshOnProc,
const BoolListType& haveMeshOnProc,
const MeshSubsetter* subsetter,
const typename GeoField::Mesh& mesh,
const autoPtr<MeshSubsetter>& subsetterPtr,
IOobjectList& allObjects,
PtrList<GeoField>& fields,
const bool deregister
@ -122,7 +122,7 @@ void Foam::fieldsDistributor::readFields
wordList masterNames(objectNames);
Pstream::broadcast(masterNames);
if (haveMeshOnProc[Pstream::myProcNo()] && objectNames != masterNames)
if (haveMeshOnProc.test(Pstream::myProcNo()) && objectNames != masterNames)
{
FatalErrorInFunction
<< "Objects not synchronised across processors." << nl
@ -173,9 +173,10 @@ void Foam::fieldsDistributor::readFields
bool decompose = true;
for (const int proci : Pstream::subProcs())
{
if (haveMeshOnProc[proci])
if (haveMeshOnProc.test(proci))
{
decompose = false;
break;
}
}
@ -197,7 +198,7 @@ void Foam::fieldsDistributor::readFields
Pstream::parRun(oldParRun); // Restore any changes
}
else if (haveMeshOnProc[Pstream::myProcNo()])
else if (haveMeshOnProc.test(Pstream::myProcNo()))
{
// Have mesh so just try to load
forAll(masterNames, i)
@ -226,20 +227,18 @@ void Foam::fieldsDistributor::readFields
OPBstream toProcs(UPstream::masterNo()); // worldComm
const label nDicts = (subsetterPtr ? fields.size() : label(0));
const label nDicts = (subsetter ? fields.size() : label(0));
toProcs << nDicts << token::BEGIN_LIST; // Begin list
if (nDicts)
if (nDicts && subsetter)
{
// Disable communication for interpolate() method
const bool oldParRun = Pstream::parRun(false);
const auto& subsetter = subsetterPtr();
forAll(fields, i)
for (const auto& fld : fields)
{
tmp<GeoField> tsubfld = subsetter.interpolate(fields[i]);
tmp<GeoField> tsubfld = subsetter->interpolate(fld);
// Surround each with {} as dictionary entry
toProcs.beginBlock();
@ -258,15 +257,14 @@ void Foam::fieldsDistributor::readFields
IPBstream fromMaster(UPstream::masterNo()); // worldComm
// But only consume where needed...
if (!haveMeshOnProc[Pstream::myProcNo()])
if (!haveMeshOnProc.test(Pstream::myProcNo()))
{
fromMaster >> fieldDicts;
}
}
// Use the received dictionaries to create fields
// (will be empty if we didn't require them)
// Use the received dictionaries (if any) to create missing fields.
// Disable communication when constructing from dictionary
const bool oldParRun = Pstream::parRun(false);
@ -327,4 +325,76 @@ void Foam::fieldsDistributor::readFields
}
template<class GeoField, class MeshSubsetter>
void Foam::fieldsDistributor::readFields
(
const bitSet& haveMeshOnProc,
const MeshSubsetter* subsetter,
const typename GeoField::Mesh& mesh,
IOobjectList& allObjects,
PtrList<GeoField>& fields,
const bool deregister
)
{
readFieldsImpl
(
haveMeshOnProc,
subsetter,
mesh,
allObjects,
fields,
deregister
);
}
template<class GeoField, class MeshSubsetter>
void Foam::fieldsDistributor::readFields
(
const boolList& haveMeshOnProc,
const MeshSubsetter* subsetter,
const typename GeoField::Mesh& mesh,
IOobjectList& allObjects,
PtrList<GeoField>& fields,
const bool deregister
)
{
readFieldsImpl
(
haveMeshOnProc,
subsetter,
mesh,
allObjects,
fields,
deregister
);
}
template<class GeoField, class MeshSubsetter>
void Foam::fieldsDistributor::readFields
(
const boolList& haveMeshOnProc,
const typename GeoField::Mesh& mesh,
const autoPtr<MeshSubsetter>& subsetter,
IOobjectList& allObjects,
PtrList<GeoField>& fields,
const bool deregister
)
{
readFieldsImpl
(
haveMeshOnProc,
subsetter.get(),
mesh,
allObjects,
fields,
deregister
);
}
// ************************************************************************* //

View File

@ -11,6 +11,7 @@ faMesh/faBoundaryMesh/faBoundaryMeshEntries.C
faMesh/faMeshSubset/faMeshSubset.C
faMesh/faMeshTools/faMeshTools.C
faMesh/faMeshTools/faMeshToolsChecks.C
faMesh/faMeshTools/faMeshToolsProcAddr.C
faPatches = faMesh/faPatches

View File

@ -978,7 +978,7 @@ bool Foam::faMesh::movePoints()
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
IOobject::NO_REGISTER
),
S()
);

View File

@ -64,7 +64,7 @@ bool Foam::faMesh::hasSystemFiles(const polyMesh& pMesh)
pMesh,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
IOobject::NO_REGISTER
),
expect // typeName (ununsed?)
)
@ -76,6 +76,7 @@ bool Foam::faMesh::hasSystemFiles(const polyMesh& pMesh)
}
}
// Only needed on master
Pstream::broadcast(looksValid);
return looksValid;
@ -87,8 +88,8 @@ bool Foam::faMesh::hasFiles(const polyMesh& pMesh)
// As well as system/{faSchemes,faSolution}
//
// expect these:
// - timeValue/faMesh/faBoundary
// - timeValue/instance/faMesh/faceLabels
// - instance/faMesh/faceLabels
// - instance/faMesh/faBoundary
bool looksValid = hasSystemFiles(pMesh);
@ -96,15 +97,23 @@ bool Foam::faMesh::hasFiles(const polyMesh& pMesh)
{
const fileOperation& fp = Foam::fileHandler();
fileName subDir(pMesh.dbDir()/faMesh::meshSubDir);
// The geometry instance for faMesh/faceLabels
// Must use READ_IF_PRESENT to avoid aborting if not available
const word instance = pMesh.time().findInstance
(
pMesh.dbDir()/faMesh::meshSubDir,
"faceLabels",
IOobject::READ_IF_PRESENT
);
for
(
const wordPair& expect
: List<wordPair>
({
{"faBoundary", "faBoundaryMesh"},
{"faceLabels", "labelList"}
{"faceLabels", "labelList"},
{"faBoundary", "faBoundaryMesh"}
})
)
{
@ -115,18 +124,18 @@ bool Foam::faMesh::hasFiles(const polyMesh& pMesh)
(
fp.filePath
(
false, // non-global
false, // non-global
IOobject
(
dataFile,
pMesh.time().findInstance(subDir, dataFile),
instance,
faMesh::meshSubDir,
pMesh,
IOobject::MUST_READ,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE,
false
IOobject::NO_REGISTER
),
dataClass // typeName (ununsed?)
dataClass // typeName (ununsed?)
)
);
@ -136,7 +145,8 @@ bool Foam::faMesh::hasFiles(const polyMesh& pMesh)
}
}
Pstream::broadcast(looksValid);
// Everybody needs it, or they all fail
Pstream::reduceAnd(looksValid);
}
return looksValid;

View File

@ -109,7 +109,7 @@ Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::newMesh
io.db(),
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
IOobject::NO_REGISTER
)
);
@ -282,7 +282,7 @@ Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::loadOrCreateMesh
io.db(),
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
IOobject::NO_REGISTER
)
);

View File

@ -31,6 +31,7 @@ Description
SourceFiles
faMeshTools.C
faMeshToolsChecks.C
faMeshToolsProcAddr.C
faMeshToolsTemplates.C
@ -134,6 +135,13 @@ public:
const GeometricField<Type, faePatchField, edgeMesh>& fld,
const bool primitiveOrdering = false
);
//- Report mesh information
static void printMeshChecks
(
const faMesh& mesh,
const int verbose = 1
);
};

View File

@ -0,0 +1,205 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "faMeshTools.H"
#include "areaFields.H"
#include "edgeFields.H"
#include "processorFaPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::faMeshTools::printMeshChecks
(
const faMesh& mesh,
const int verbose
)
{
const faBoundaryMesh& patches = mesh.boundary();
const label nNonProcessor = patches.nNonProcessor();
const label nPatches = patches.size();
label nLocalProcEdges = 0;
if (Pstream::parRun())
{
for (const faPatch& fap : patches)
{
const auto* cpp = isA<processorFaPatch>(fap);
if (cpp)
{
nLocalProcEdges += fap.nEdges();
}
}
}
const labelList nFaces
(
UPstream::listGatherValues<label>(mesh.nFaces())
);
const labelList nPoints
(
UPstream::listGatherValues<label>(mesh.nPoints())
);
const labelList nEdges
(
UPstream::listGatherValues<label>(mesh.nEdges())
);
const labelList nIntEdges
(
UPstream::listGatherValues<label>(mesh.nInternalEdges())
);
// The "real" (non-processor) boundary edges
const labelList nBndEdges
(
UPstream::listGatherValues<label>
(
mesh.nBoundaryEdges() - nLocalProcEdges
)
);
const labelList nProcEdges
(
UPstream::listGatherValues<label>(nLocalProcEdges)
);
// Format output as
// Number of faces: ...
// per-proc: (...)
const auto reporter =
[&,verbose](const char* tag, const labelList& list)
{
Info<< " Number of " << tag << ": " << sum(list) << nl;
if (Pstream::parRun() && verbose)
{
int padding = static_cast<int>
(
// strlen(" Number of ") - strlen("per-proc")
(12 - 8)
+ strlen(tag)
);
do { Info<< ' '; } while (--padding > 0);
Info<< "per-proc: " << flatOutput(list) << nl;
}
};
Info<< "----------------" << nl
<< "Mesh Information" << nl
<< "----------------" << nl
<< " " << "boundingBox: " << boundBox(mesh.points()) << nl;
if (Pstream::master())
{
reporter("faces", nFaces);
reporter("points", nPoints);
reporter("edges", nEdges);
reporter("internal edges", nIntEdges);
reporter("boundary edges", nBndEdges);
if (Pstream::parRun())
{
reporter("processor edges", nProcEdges);
}
}
Info<< "----------------" << nl
<< "Patches" << nl
<< "----------------" << nl;
for (label patchi = 0; patchi < nNonProcessor; ++patchi)
{
const faPatch& p = patches[patchi];
// Report physical size (nEdges) not virtual size
Info<< " " << "patch " << p.index()
<< " (size: " << returnReduce(p.nEdges(), sumOp<label>())
<< ") name: " << p.name()
<< nl;
}
Info<< "----------------" << nl
<< "Used polyPatches: " << flatOutput(mesh.whichPolyPatches()) << nl;
// Geometry information
Info<< nl;
{
scalarMinMax limit(gMinMax(mesh.S().field()));
Info<< "Face area:" << nl
<< " min = " << limit.min() << " max = " << limit.max() << nl;
}
{
scalarMinMax limit(minMax(mesh.magLe().primitiveField()));
// Include processor boundaries into 'internal' edges
if (Pstream::parRun())
{
for (label patchi = nNonProcessor; patchi < nPatches; ++patchi)
{
limit.add(minMax(mesh.magLe().boundaryField()[patchi]));
}
reduce(limit, minMaxOp<scalar>());
}
Info<< "Edge length (internal):" << nl
<< " min = " << limit.min() << " max = " << limit.max() << nl;
// Include (non-processor) boundaries
for (label patchi = 0; patchi < nNonProcessor; ++patchi)
{
limit.add(minMax(mesh.magLe().boundaryField()[patchi]));
}
if (Pstream::parRun())
{
reduce(limit, minMaxOp<scalar>());
}
Info<< "Edge length:" << nl
<< " min = " << limit.min() << " max = " << limit.max() << nl;
}
// Not particularly meaningful
#if 0
{
MinMax<vector> limit(gMinMax(mesh.faceAreaNormals().field()));
Info<< "Face area normals:" << nl
<< " min = " << limit.min() << " max = " << limit.max() << nl;
}
#endif
}
// ************************************************************************* //

View File

@ -238,7 +238,7 @@ Foam::faMeshTools::readProcAddressing
mesh.thisDb(),
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE,
false // no register
IOobject::NO_REGISTER
);
//if (ioAddr.typeHeaderOk<labelIOList>(true))
@ -327,7 +327,7 @@ void Foam::faMeshTools::writeProcAddressing
(procMesh && !decompose ? procMesh->thisDb() : mesh.thisDb()),
IOobject::NO_READ,
IOobject::NO_WRITE,
false // no register
IOobject::NO_REGISTER
);
@ -418,7 +418,7 @@ void Foam::faMeshTools::writeProcAddressing
mesh.thisDb(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false // no register
IOobject::NO_REGISTER
),
map
);

View File

@ -45,6 +45,7 @@ SourceFiles
#define Foam_faFieldDecomposer_H
#include "faMesh.H"
#include "faMeshSubset.H"
#include "faPatchFieldMapper.H"
#include "edgeFields.H"
@ -446,6 +447,26 @@ public:
const IOobjectList& objects
);
//- Read all fields given mesh and objects.
//- Supports reading/sending fields
void readAllFields
(
const bitSet& haveMeshOnProc,
const faMeshSubset* subsetter,
const faMesh& mesh,
IOobjectList& objects
);
//- Read all fields given mesh and objects.
//- Supports reading/sending fields
void readAllFields
(
const boolList& haveMeshOnProc,
const faMeshSubset* subsetter,
const faMesh& mesh,
IOobjectList& objects
);
//- Decompose and write all fields
void decomposeAllFields
(

View File

@ -77,7 +77,11 @@ public:
bool empty() const noexcept { return !size(); }
void readAll(const faMesh& mesh, const IOobjectList& objects)
void readAll
(
const faMesh& mesh,
const IOobjectList& objects
)
{
#undef doLocalCode
#define doLocalCode(Type) \
@ -113,6 +117,45 @@ public:
#undef doLocalCode
}
template<class BoolListType>
void readAll
(
const BoolListType& haveMeshOnProc,
const faMeshSubset*& subsetter,
const faMesh& mesh,
IOobjectList& objects
)
{
#undef doLocalCode
#define doLocalCode(Type) \
{ \
fieldsDistributor::readFields \
( \
haveMeshOnProc, \
subsetter, \
mesh, \
objects, \
Type##AreaFields_ \
); \
fieldsDistributor::readFields \
( \
haveMeshOnProc, \
subsetter, \
mesh, \
objects, \
Type##EdgeFields_ \
); \
}
doLocalCode(scalar);
doLocalCode(vector);
doLocalCode(sphericalTensor);
doLocalCode(symmTensor);
doLocalCode(tensor);
#undef doLocalCode
}
template<class GeoField>
static void decompose
(
@ -205,6 +248,36 @@ void Foam::faFieldDecomposer::fieldsCache::readAllFields
}
void Foam::faFieldDecomposer::fieldsCache::readAllFields
(
const bitSet& haveMeshOnProc,
const faMeshSubset* subsetter,
const faMesh& mesh,
IOobjectList& objects
)
{
if (cache_)
{
cache_->readAll(haveMeshOnProc, subsetter, mesh, objects);
}
}
void Foam::faFieldDecomposer::fieldsCache::readAllFields
(
const boolList& haveMeshOnProc,
const faMeshSubset* subsetter,
const faMesh& mesh,
IOobjectList& objects
)
{
if (cache_)
{
cache_->readAll(haveMeshOnProc, subsetter, mesh, objects);
}
}
void Foam::faFieldDecomposer::fieldsCache::decomposeAllFields
(
const faFieldDecomposer& decomposer,