zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

BUG: lduPrimitiveMesh: multi-agglomeration failing to create interfaces; memory leaks

Browse Source
This commit is contained in:
mattijs
2013-03-04 18:32:59 +00:00
parent 549d26f307
commit 62930bb4bf
13 changed files with 1061 additions and 746 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
View File

@ -257,10 +257,18 @@ const Foam::lduInterfacePtrsList& Foam::GAMGAgglomeration::interfaceLevel
void Foam::GAMGAgglomeration::gatherMeshes
(
const labelList& procAgglomMap,
const labelList& procIDs,
const label allMeshComm
) const
{
if (allMeshPtr_.valid())
{
FatalErrorIn("GAMGAgglomeration::gatherMeshes(..)")
<< "Processor-agglomerated mesh already constructed"
<< exit(FatalError);
}
const lduMesh& coarsestMesh = meshLevels_.last();
label coarseComm = coarsestMesh.comm();
@ -293,6 +301,8 @@ void Foam::GAMGAgglomeration::gatherMeshes
new lduPrimitiveMesh
(
allMeshComm,
procAgglomMap,
procIDs,
coarsestMesh,
otherMeshes_,
@ -303,6 +313,8 @@ void Foam::GAMGAgglomeration::gatherMeshes
boundaryFaceMap_
)
);
Pout<< "** Agglomerated Mesh " << allMeshPtr_().info() << endl;
}
UPstream::warnComm = oldWarn;

9
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H
View File

@ -278,10 +278,15 @@ public:
// Procesor agglomeration
//- Collect and combine processor meshes into allMesh. Pass in
// communicator for combined mesh.
//- Collect and combine processor meshes into allMesh.
// - allMeshComm : communicator for combined mesh.
// - procAgglomMap : per processor the new agglomerated processor
// (rank in allMeshComm!). Global information.
// - procIDs : local information: same for all in
// agglomerated processor.
void gatherMeshes
(
const labelList& procAgglomMap,
const labelList& procIDs,
const label allMeshComm
) const;

536
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C
View File

@ -24,10 +24,6 @@ License
\*---------------------------------------------------------------------------*/
#include "GAMGSolver.H"
#include "processorGAMGInterface.H"
#include "processorLduInterfaceField.H"
#include "GAMGInterfaceField.H"
#include "processorGAMGInterfaceField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -43,149 +39,6 @@ namespace Foam
}
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
// Gather matrices.
// Note: matrices get constructed with dummy mesh
void Foam::GAMGSolver::gatherMatrices
(
const labelList& procIDs,
const lduMesh& dummyMesh,
const label meshComm,
const lduMatrix& mat,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const FieldField<Field, scalar>& interfaceIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
PtrList<lduMatrix>& otherMats,
PtrList<FieldField<Field, scalar> >& otherBouCoeffs,
PtrList<FieldField<Field, scalar> >& otherIntCoeffs,
List<boolList>& otherTransforms,
List<List<int> >& otherRanks
) const
{
Pout<< "GAMGSolver::gatherMatrices :"
<< " collecting matrices on procs:" << procIDs
<< " using comm:" << meshComm << endl;
if (Pstream::myProcNo(meshComm) == procIDs[0])
{
// Master.
Pout<< "GAMGSolver::gatherMatrices :"
<< " master:" << Pstream::myProcNo(meshComm) << endl;
otherMats.setSize(procIDs.size()-1);
otherBouCoeffs.setSize(procIDs.size()-1);
otherIntCoeffs.setSize(procIDs.size()-1);
otherTransforms.setSize(procIDs.size()-1);
otherRanks.setSize(procIDs.size()-1);
for (label procI = 1; procI < procIDs.size(); procI++)
{
label otherI = procI-1;
IPstream fromSlave
(
Pstream::scheduled,
procIDs[procI],
0, // bufSize
Pstream::msgType(),
meshComm
);
otherMats.set(otherI, new lduMatrix(dummyMesh, fromSlave));
// Receive number of/valid interfaces
boolList& procTransforms = otherTransforms[otherI];
List<int>& procRanks = otherRanks[otherI];
fromSlave >> procTransforms;
fromSlave >> procRanks;
// Size coefficients
otherBouCoeffs.set
(
otherI,
new FieldField<Field, scalar>(procRanks.size())
);
otherIntCoeffs.set
(
otherI,
new FieldField<Field, scalar>(procRanks.size())
);
forAll(procRanks, intI)
{
if (procRanks[intI] != -1)
{
otherBouCoeffs[otherI].set
(
intI,
new scalarField(fromSlave)
);
otherIntCoeffs[otherI].set
(
intI,
new scalarField(fromSlave)
);
}
}
}
}
else
{
Pout<< "GAMGSolver::gatherMatrices :"
<< " sending from:" << Pstream::myProcNo(meshComm)
<< " to master:" << procIDs[0] << endl;
// Send to master
// Count valid interfaces
boolList procTransforms(interfaceBouCoeffs.size(), false);
List<int> procRanks(interfaceBouCoeffs.size(), -1);
forAll(interfaces, intI)
{
if (interfaces.set(intI))
{
const processorLduInterfaceField& interface =
refCast<const processorLduInterfaceField>
(
interfaces[intI]
);
procTransforms[intI] = interface.doTransform();
procRanks[intI] = interface.rank();
}
}
Pout<< "GAMGSolver::gatherMatrices :"
<< " sending matrix:" << mat.info() << endl;
OPstream toMaster
(
Pstream::scheduled,
procIDs[0],
0,
Pstream::msgType(),
meshComm
);
toMaster << mat << procTransforms << procRanks;
forAll(procRanks, intI)
{
if (procRanks[intI] != -1)
{
toMaster
<< interfaceBouCoeffs[intI]
<< interfaceIntCoeffs[intI];
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::GAMGSolver::GAMGSolver
@ -282,342 +135,99 @@ Foam::GAMGSolver::GAMGSolver
Pout<< "Solve generic on coarsestmesh (level=" << coarsestLevel
<< ") using communicator " << coarseComm << endl;
// All processors to master
const List<int>& procIDs = UPstream::procID(coarseComm);
Pout<< "procIDs:" << procIDs << endl;
//const List<int>& procIDs = UPstream::procID(coarseComm);
// Processor restriction map: per processor the coarse processor
labelList procAgglomMap(UPstream::nProcs(coarseComm));
procAgglomMap[0] = 0;
procAgglomMap[1] = 0;
procAgglomMap[2] = 1;
procAgglomMap[3] = 1;
// Determine the master processors
Map<label> agglomToMaster(procAgglomMap.size());
forAll(procAgglomMap, procI)
{
label coarseI = procAgglomMap[procI];
Map<label>::iterator fnd = agglomToMaster.find(coarseI);
if (fnd == agglomToMaster.end())
{
agglomToMaster.insert(coarseI, procI);
}
else
{
fnd() = max(fnd(), procI);
}
}
labelList masterProcs(agglomToMaster.size());
forAllConstIter(Map<label>, agglomToMaster, iter)
{
masterProcs[iter.key()] = iter();
}
// Allocate a communicator for the linear solver
// Allocate a communicator for single processor
label masterComm = UPstream::allocateCommunicator
(
coarseComm,
labelList(1, 0)
masterProcs
);
Pout<< "** Allocated communicator " << masterComm
<< " for processor " << procIDs[0]
<< " out of " << procIDs << endl;
<< " for indices " << masterProcs
<< " in processor list " << UPstream::procID(coarseComm)
<< endl;
// Gather all meshes onto procIDs[0]
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
agglomeration_.gatherMeshes(procIDs, masterComm);
List<int> agglomProcIDs;
// Collect all the processors in my agglomeration
{
label myProcID = Pstream::myProcNo(coarseComm);
label myAgglom = procAgglomMap[myProcID];
// Get all processors agglomerating to the same coarse processor
agglomProcIDs = findIndices(procAgglomMap, myAgglom);
// Make sure the master is the first element.
label index = findIndex
(
agglomProcIDs,
agglomToMaster[myAgglom]
);
Swap(agglomProcIDs[0], agglomProcIDs[index]);
}
Pout<< "procAgglomMap:" << procAgglomMap << endl;
Pout<< "agglomProcIDs:" << agglomProcIDs << endl;
// Gather all matrix coefficients onto procIDs[0]
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Gather matrix and mesh onto agglomProcIDs[0]
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
PtrList<lduMatrix> otherMats;
PtrList<FieldField<Field, scalar> > otherBouCoeffs;
PtrList<FieldField<Field, scalar> > otherIntCoeffs;
List<boolList> otherTransforms;
List<List<int> > otherRanks;
gatherMatrices
procAgglomerateMatrix
(
procIDs,
// Current mesh and matrix
coarsestMesh,
coarseComm,
coarsestMatrix,
coarsestInterfaces,
coarsestBouCoeffs,
coarsestIntCoeffs,
coarsestInterfaces,
otherMats,
otherBouCoeffs,
otherIntCoeffs,
otherTransforms,
otherRanks
// Agglomeration information
procAgglomMap,
agglomProcIDs,
masterComm,
// Resulting matrix
allMatrixPtr_,
allInterfaceBouCoeffs_,
allInterfaceIntCoeffs_,
allPrimitiveInterfaces_,
allInterfaces_
);
if (Pstream::myProcNo(coarseComm) == procIDs[0])
{
// Agglomerate all matrix
// ~~~~~~~~~~~~~~~~~~~~~~
Pout<< "Own matrix:" << coarsestMatrix.info() << endl;
forAll(otherMats, i)
{
Pout<< "** otherMats " << i << " "
<< otherMats[i].info()
<< endl;
}
Pout<< endl;
const lduMesh& allMesh = agglomeration_.allMesh();
allMatrixPtr_.reset(new lduMatrix(allMesh));
lduMatrix& allMatrix = allMatrixPtr_();
if (coarsestMatrix.hasDiag())
{
scalarField& allDiag = allMatrix.diag();
SubList<scalar>
(
allDiag,
coarsestMatrix.diag().size()
).assign
(
coarsestMatrix.diag()
);
forAll(otherMats, i)
{
SubList<scalar>
(
allDiag,
otherMats[i].diag().size(),
agglomeration_.cellOffsets()[i+1]
).assign
(
otherMats[i].diag()
);
}
}
if (coarsestMatrix.hasLower())
{
scalarField& allLower = allMatrix.lower();
UIndirectList<scalar>
(
allLower,
agglomeration_.faceMap()[0]
) = coarsestMatrix.lower();
forAll(otherMats, i)
{
UIndirectList<scalar>
(
allLower,
agglomeration_.faceMap()[i+1]
) = otherMats[i].lower();
}
}
if (coarsestMatrix.hasUpper())
{
scalarField& allUpper = allMatrix.upper();
UIndirectList<scalar>
(
allUpper,
agglomeration_.faceMap()[0]
) = coarsestMatrix.upper();
forAll(otherMats, i)
{
UIndirectList<scalar>
(
allUpper,
agglomeration_.faceMap()[i+1]
) = otherMats[i].upper();
}
}
// Agglomerate interface fields and coefficients
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
lduInterfacePtrsList allMeshInterfaces = allMesh.interfaces();
allInterfaceBouCoeffs_.setSize(allMeshInterfaces.size());
allInterfaceIntCoeffs_.setSize(allMeshInterfaces.size());
allInterfaces_.setSize(allMeshInterfaces.size());
forAll(allMeshInterfaces, intI)
{
const lduInterface& patch = allMeshInterfaces[intI];
label size = patch.faceCells().size();
allInterfaceBouCoeffs_.set(intI, new scalarField(size));
allInterfaceIntCoeffs_.set(intI, new scalarField(size));
}
labelList nBounFaces(allMeshInterfaces.size());
forAll(agglomeration_.boundaryMap(), procI)
{
const FieldField<Field, scalar>& procBouCoeffs
(
(procI == 0)
? coarsestBouCoeffs
: otherBouCoeffs[procI-1]
);
const FieldField<Field, scalar>& procIntCoeffs
(
(procI == 0)
? coarsestIntCoeffs
: otherIntCoeffs[procI-1]
);
const labelList& bMap = agglomeration_.boundaryMap()[procI];
forAll(bMap, procIntI)
{
label allIntI = bMap[procIntI];
if (allIntI != -1)
{
// So this boundary has been preserved. Copy
// data across.
if (!allInterfaces_.set(allIntI))
{
// Construct lduInterfaceField
bool doTransform = false;
int rank = -1;
if (procI == 0)
{
const processorGAMGInterfaceField& procInt =
refCast
<
const processorGAMGInterfaceField
>
(
coarsestInterfaces[procIntI]
);
doTransform = procInt.doTransform();
rank = procInt.rank();
}
else
{
doTransform =
otherTransforms[procI-1][procIntI];
rank = otherRanks[procI-1][procIntI];
}
allInterfaces_.set
(
allIntI,
GAMGInterfaceField::New
(
refCast<const GAMGInterface>
(
allMeshInterfaces[allIntI]
),
doTransform,
rank
).ptr()
);
}
// Map data from processor to complete mesh
scalarField& allBou =
allInterfaceBouCoeffs_[allIntI];
scalarField& allInt =
allInterfaceIntCoeffs_[allIntI];
const labelList& map = agglomeration_.
boundaryFaceMap()[procI][procIntI];
const scalarField& procBou =
procBouCoeffs[procIntI];
const scalarField& procInt =
procIntCoeffs[procIntI];
forAll(map, i)
{
label allFaceI = map[i];
if (allFaceI < 0)
{
FatalErrorIn("GAMGSolver::GAMGSolver()")
<< "problem." << abort(FatalError);
}
allBou[allFaceI] = procBou[i];
allInt[allFaceI] = procInt[i];
}
}
else if (procBouCoeffs.set(procIntI))
{
// Boundary has become internal face
const labelList& map = agglomeration_.
boundaryFaceMap()[procI][procIntI];
const labelList& allU =
allMesh.lduAddr().upperAddr();
const labelList& allL =
allMesh.lduAddr().lowerAddr();
const label off = agglomeration_.
cellOffsets()[procI];
const scalarField& procBou =
procBouCoeffs[procIntI];
const scalarField& procInt =
procIntCoeffs[procIntI];
forAll(map, i)
{
if (map[i] >= 0)
{
label allFaceI = map[i];
if (coarsestMatrix.hasUpper())
{
allMatrix.upper()[allFaceI] =
-procBou[i];
}
if (coarsestMatrix.hasLower())
{
allMatrix.lower()[allFaceI] =
-procInt[i];
}
}
else
{
label allFaceI = -map[i]-1;
if (coarsestMatrix.hasUpper())
{
allMatrix.upper()[allFaceI] =
-procInt[i];
}
if (coarsestMatrix.hasLower())
{
allMatrix.lower()[allFaceI] =
-procBou[i];
}
}
// Simple check
label allFaceI =
(
map[i] >= 0
? map[i]
: -map[i]-1
);
const labelList& fCells =
lduPrimitiveMesh::mesh
(
coarsestMesh,
agglomeration_.otherMeshes(),
procI
).lduAddr().patchAddr(procIntI);
label allCellI = off + fCells[i];
if
(
allCellI != allL[allFaceI]
&& allCellI != allU[allFaceI]
)
{
FatalErrorIn
(
"GAMGSolver::GAMGSolver()"
) << "problem."
<< " allFaceI:" << allFaceI
<< " local cellI:" << fCells[i]
<< " allCellI:" << allCellI
<< " allLower:" << allL[allFaceI]
<< " allUpper:" << allU[allFaceI]
<< abort(FatalError);
}
}
}
}
}
Pout<< "** Assembled allMatrix:" << allMatrix.info() << endl;
}
UPstream::warnComm = oldWarn;
}
}

41
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
View File

@ -137,6 +137,7 @@ class GAMGSolver
autoPtr<lduMatrix> allMatrixPtr_;
FieldField<Field, scalar> allInterfaceBouCoeffs_;
FieldField<Field, scalar> allInterfaceIntCoeffs_;
PtrList<lduInterfaceField> allPrimitiveInterfaces_;
lduInterfaceFieldPtrsList allInterfaces_;
@ -171,23 +172,6 @@ class GAMGSolver
//- Collect matrices from other processors
void gatherMatrices
(
const labelList& procIDs,
const lduMesh& myMesh,
const PtrList<lduMesh>& otherMeshes,
const label meshComm,
const lduMatrix& mat,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const FieldField<Field, scalar>& interfaceIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
PtrList<lduMatrix>& otherMats,
PtrList<FieldField<Field, scalar> >& otherBouCoeffs,
PtrList<FieldField<Field, scalar> >& otherIntCoeffs,
PtrList<lduInterfaceFieldPtrsList>& otherInterfaces
) const;
void gatherMatrices
(
const labelList& procIDs,
const lduMesh& dummyMesh,
@ -205,6 +189,29 @@ class GAMGSolver
List<List<int> >& otherRanks
) const;
//- Agglomerate processor matrices
void procAgglomerateMatrix
(
// Current mesh and matrix
const lduMesh& coarsestMesh,
const lduMatrix& coarsestMatrix,
const lduInterfaceFieldPtrsList& coarsestInterfaces,
const FieldField<Field, scalar>& coarsestBouCoeffs,
const FieldField<Field, scalar>& coarsestIntCoeffs,
// Agglomeration information
const labelList& procAgglomMap,
const List<int>& agglomProcIDs,
label masterComm,
// Resulting matrix
autoPtr<lduMatrix>& allMatrixPtr,
FieldField<Field, scalar>& allInterfaceBouCoeffs,
FieldField<Field, scalar>& allInterfaceIntCoeffs,
PtrList<lduInterfaceField>& allPrimitiveInterfaces,
lduInterfaceFieldPtrsList& allInterfaces
);
//- Interpolate the correction after injected prolongation
void interpolate
(

515
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C
View File

@ -25,6 +25,8 @@ License
#include "GAMGSolver.H"
#include "GAMGInterfaceField.H"
#include "processorLduInterfaceField.H"
#include "processorGAMGInterfaceField.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -207,4 +209,517 @@ void Foam::GAMGSolver::agglomerateMatrix(const label fineLevelIndex)
}
// Gather matrices.
// Note: matrices get constructed with dummy mesh
void Foam::GAMGSolver::gatherMatrices
(
const labelList& procIDs,
const lduMesh& dummyMesh,
const label meshComm,
const lduMatrix& mat,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const FieldField<Field, scalar>& interfaceIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
PtrList<lduMatrix>& otherMats,
PtrList<FieldField<Field, scalar> >& otherBouCoeffs,
PtrList<FieldField<Field, scalar> >& otherIntCoeffs,
List<boolList>& otherTransforms,
List<List<int> >& otherRanks
) const
{
Pout<< "GAMGSolver::gatherMatrices :"
<< " collecting matrices on procs:" << procIDs
<< " using comm:" << meshComm << endl;
if (Pstream::myProcNo(meshComm) == procIDs[0])
{
// Master.
Pout<< "GAMGSolver::gatherMatrices :"
<< " master:" << Pstream::myProcNo(meshComm) << endl;
otherMats.setSize(procIDs.size()-1);
otherBouCoeffs.setSize(procIDs.size()-1);
otherIntCoeffs.setSize(procIDs.size()-1);
otherTransforms.setSize(procIDs.size()-1);
otherRanks.setSize(procIDs.size()-1);
for (label procI = 1; procI < procIDs.size(); procI++)
{
label otherI = procI-1;
IPstream fromSlave
(
Pstream::scheduled,
procIDs[procI],
0, // bufSize
Pstream::msgType(),
meshComm
);
otherMats.set(otherI, new lduMatrix(dummyMesh, fromSlave));
// Receive number of/valid interfaces
boolList& procTransforms = otherTransforms[otherI];
List<int>& procRanks = otherRanks[otherI];
fromSlave >> procTransforms;
fromSlave >> procRanks;
// Size coefficients
otherBouCoeffs.set
(
otherI,
new FieldField<Field, scalar>(procRanks.size())
);
otherIntCoeffs.set
(
otherI,
new FieldField<Field, scalar>(procRanks.size())
);
forAll(procRanks, intI)
{
if (procRanks[intI] != -1)
{
otherBouCoeffs[otherI].set
(
intI,
new scalarField(fromSlave)
);
otherIntCoeffs[otherI].set
(
intI,
new scalarField(fromSlave)
);
}
}
}
}
else
{
Pout<< "GAMGSolver::gatherMatrices :"
<< " sending from:" << Pstream::myProcNo(meshComm)
<< " to master:" << procIDs[0] << endl;
// Send to master
// Count valid interfaces
boolList procTransforms(interfaceBouCoeffs.size(), false);
List<int> procRanks(interfaceBouCoeffs.size(), -1);
forAll(interfaces, intI)
{
if (interfaces.set(intI))
{
const processorLduInterfaceField& interface =
refCast<const processorLduInterfaceField>
(
interfaces[intI]
);
procTransforms[intI] = interface.doTransform();
procRanks[intI] = interface.rank();
}
}
//Pout<< "GAMGSolver::gatherMatrices :"
// << " sending matrix:" << mat.info() << endl;
OPstream toMaster
(
Pstream::scheduled,
procIDs[0],
0,
Pstream::msgType(),
meshComm
);
toMaster << mat << procTransforms << procRanks;
forAll(procRanks, intI)
{
if (procRanks[intI] != -1)
{
toMaster
<< interfaceBouCoeffs[intI]
<< interfaceIntCoeffs[intI];
}
}
}
}
void Foam::GAMGSolver::procAgglomerateMatrix
(
// Current mesh and matrix
const lduMesh& coarsestMesh,
const lduMatrix& coarsestMatrix,
const lduInterfaceFieldPtrsList& coarsestInterfaces,
const FieldField<Field, scalar>& coarsestBouCoeffs,
const FieldField<Field, scalar>& coarsestIntCoeffs,
// Agglomeration information
const labelList& procAgglomMap,
const List<int>& agglomProcIDs,
label masterComm,
// Resulting matrix
autoPtr<lduMatrix>& allMatrixPtr,
FieldField<Field, scalar>& allInterfaceBouCoeffs,
FieldField<Field, scalar>& allInterfaceIntCoeffs,
PtrList<lduInterfaceField>& allPrimitiveInterfaces,
lduInterfaceFieldPtrsList& allInterfaces
)
{
label coarseComm = coarsestMesh.comm();
label oldWarn = UPstream::warnComm;
UPstream::warnComm = coarseComm;
// Construct (on the agglomeration) a complete mesh with mapping
agglomeration_.gatherMeshes(procAgglomMap, agglomProcIDs, masterComm);
// Gather all matrix coefficients onto agglomProcIDs[0]
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
PtrList<lduMatrix> otherMats;
PtrList<FieldField<Field, scalar> > otherBouCoeffs;
PtrList<FieldField<Field, scalar> > otherIntCoeffs;
List<boolList> otherTransforms;
List<List<int> > otherRanks;
gatherMatrices
(
agglomProcIDs,
coarsestMesh,
coarseComm,
coarsestMatrix,
coarsestBouCoeffs,
coarsestIntCoeffs,
coarsestInterfaces,
otherMats,
otherBouCoeffs,
otherIntCoeffs,
otherTransforms,
otherRanks
);
if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
{
// Agglomerate all matrix
// ~~~~~~~~~~~~~~~~~~~~~~
Pout<< "Own matrix:" << coarsestMatrix.info() << endl;
forAll(otherMats, i)
{
Pout<< "** otherMats " << i << " "
<< otherMats[i].info()
<< endl;
}
Pout<< endl;
const lduMesh& allMesh = agglomeration_.allMesh();
allMatrixPtr.reset(new lduMatrix(allMesh));
lduMatrix& allMatrix = allMatrixPtr();
if (coarsestMatrix.hasDiag())
{
scalarField& allDiag = allMatrix.diag();
SubList<scalar>
(
allDiag,
coarsestMatrix.diag().size()
).assign
(
coarsestMatrix.diag()
);
forAll(otherMats, i)
{
SubList<scalar>
(
allDiag,
otherMats[i].diag().size(),
agglomeration_.cellOffsets()[i+1]
).assign
(
otherMats[i].diag()
);
}
}
if (coarsestMatrix.hasLower())
{
scalarField& allLower = allMatrix.lower();
UIndirectList<scalar>
(
allLower,
agglomeration_.faceMap()[0]
) = coarsestMatrix.lower();
forAll(otherMats, i)
{
UIndirectList<scalar>
(
allLower,
agglomeration_.faceMap()[i+1]
) = otherMats[i].lower();
}
}
if (coarsestMatrix.hasUpper())
{
scalarField& allUpper = allMatrix.upper();
UIndirectList<scalar>
(
allUpper,
agglomeration_.faceMap()[0]
) = coarsestMatrix.upper();
forAll(otherMats, i)
{
UIndirectList<scalar>
(
allUpper,
agglomeration_.faceMap()[i+1]
) = otherMats[i].upper();
}
}
// Agglomerate interface fields and coefficients
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
lduInterfacePtrsList allMeshInterfaces = allMesh.interfaces();
allInterfaceBouCoeffs.setSize(allMeshInterfaces.size());
allInterfaceIntCoeffs.setSize(allMeshInterfaces.size());
allPrimitiveInterfaces.setSize(allMeshInterfaces.size());
allInterfaces.setSize(allMeshInterfaces.size());
forAll(allMeshInterfaces, intI)
{
const lduInterface& patch = allMeshInterfaces[intI];
label size = patch.faceCells().size();
allInterfaceBouCoeffs.set(intI, new scalarField(size));
allInterfaceIntCoeffs.set(intI, new scalarField(size));
}
labelList nBounFaces(allMeshInterfaces.size());
forAll(agglomeration_.boundaryMap(), procI)
{
const FieldField<Field, scalar>& procBouCoeffs
(
(procI == 0)
? coarsestBouCoeffs
: otherBouCoeffs[procI-1]
);
const FieldField<Field, scalar>& procIntCoeffs
(
(procI == 0)
? coarsestIntCoeffs
: otherIntCoeffs[procI-1]
);
const labelList& bMap = agglomeration_.boundaryMap()[procI];
forAll(bMap, procIntI)
{
label allIntI = bMap[procIntI];
if (allIntI != -1)
{
// So this boundary has been preserved. Copy
// data across.
if (!allInterfaces.set(allIntI))
{
// Construct lduInterfaceField
bool doTransform = false;
int rank = -1;
if (procI == 0)
{
const processorGAMGInterfaceField& procInt =
refCast
<
const processorGAMGInterfaceField
>
(
coarsestInterfaces[procIntI]
);
doTransform = procInt.doTransform();
rank = procInt.rank();
}
else
{
doTransform =
otherTransforms[procI-1][procIntI];
rank = otherRanks[procI-1][procIntI];
}
allPrimitiveInterfaces.set
(
allIntI,
GAMGInterfaceField::New
(
refCast<const GAMGInterface>
(
allMeshInterfaces[allIntI]
),
doTransform,
rank
).ptr()
);
allInterfaces.set
(
allIntI,
&allPrimitiveInterfaces[allIntI]
);
}
// Map data from processor to complete mesh
scalarField& allBou = allInterfaceBouCoeffs[allIntI];
scalarField& allInt = allInterfaceIntCoeffs[allIntI];
const labelList& map =
agglomeration_.boundaryFaceMap()[procI][procIntI];
Pout<< " from proc:" << procI << " interface:" << procIntI
<< " mapped to faces:" << map << endl;
const scalarField& procBou = procBouCoeffs[procIntI];
const scalarField& procInt = procIntCoeffs[procIntI];
forAll(map, i)
{
label allFaceI = map[i];
if (allFaceI < 0)
{
FatalErrorIn("GAMGSolver::GAMGSolver()")
<< "problem." << abort(FatalError);
}
allBou[allFaceI] = procBou[i];
allInt[allFaceI] = procInt[i];
}
}
else if (procBouCoeffs.set(procIntI))
{
// Boundary has become internal face
const labelList& map =
agglomeration_.boundaryFaceMap()[procI][procIntI];
const labelList& allU = allMesh.lduAddr().upperAddr();
const labelList& allL = allMesh.lduAddr().lowerAddr();
const label off = agglomeration_.cellOffsets()[procI];
const scalarField& procBou = procBouCoeffs[procIntI];
const scalarField& procInt = procIntCoeffs[procIntI];
forAll(map, i)
{
if (map[i] >= 0)
{
label allFaceI = map[i];
if (coarsestMatrix.hasUpper())
{
allMatrix.upper()[allFaceI] = -procBou[i];
}
if (coarsestMatrix.hasLower())
{
allMatrix.lower()[allFaceI] = -procInt[i];
}
}
else
{
label allFaceI = -map[i]-1;
if (coarsestMatrix.hasUpper())
{
allMatrix.upper()[allFaceI] = -procInt[i];
}
if (coarsestMatrix.hasLower())
{
allMatrix.lower()[allFaceI] = -procBou[i];
}
}
// Simple check
label allFaceI =
(
map[i] >= 0
? map[i]
: -map[i]-1
);
const labelList& fCells =
lduPrimitiveMesh::mesh
(
coarsestMesh,
agglomeration_.otherMeshes(),
procI
).lduAddr().patchAddr(procIntI);
label allCellI = off + fCells[i];
if
(
allCellI != allL[allFaceI]
&& allCellI != allU[allFaceI]
)
{
FatalErrorIn
(
"GAMGSolver::GAMGSolver()"
) << "problem."
<< " allFaceI:" << allFaceI
<< " local cellI:" << fCells[i]
<< " allCellI:" << allCellI
<< " allLower:" << allL[allFaceI]
<< " allUpper:" << allU[allFaceI]
<< abort(FatalError);
}
}
}
}
}
Pout<< "** Assembled allMatrix:" << allMatrix.info() << endl;
forAll(allInterfaces, intI)
{
if (allInterfaces.set(intI))
{
Pout<< " patch:" << intI
<< " type:" << allInterfaces[intI].type()
<< " size:"
<< allInterfaces[intI].interface().
faceCells().size()
<< endl;
const scalarField& bouCoeffs = allInterfaceBouCoeffs[intI];
const scalarField& intCoeffs = allInterfaceIntCoeffs[intI];
forAll(bouCoeffs, faceI)
{
Pout<< " " << faceI
<< "\tbou:" << bouCoeffs[faceI]
<< "\tint:" << intCoeffs[faceI]
<< endl;
}
}
}
}
UPstream::warnComm = oldWarn;
}
//XXXX
// ************************************************************************* //

55
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
View File

@ -453,9 +453,49 @@ void Foam::GAMGSolver::solveCoarsestLevel
}
else if (processorAgglomerate_)
{
Pout<< "** Processor agglomeration" << endl;
const List<int>& procIDs = UPstream::procID(coarseComm);
//XXXX
labelList procAgglomMap(UPstream::nProcs(coarseComm));
procAgglomMap[0] = 0;
procAgglomMap[1] = 0;
procAgglomMap[2] = 1;
procAgglomMap[3] = 1;
// Determine the master processors
Map<label> agglomToMaster(procAgglomMap.size());
forAll(procAgglomMap, procI)
{
label coarseI = procAgglomMap[procI];
Map<label>::iterator fnd = agglomToMaster.find(coarseI);
if (fnd == agglomToMaster.end())
{
agglomToMaster.insert(coarseI, procI);
}
else
{
fnd() = max(fnd(), procI);
}
}
List<int> agglomProcIDs;
// Collect all the processors in my agglomeration
{
label myProcID = Pstream::myProcNo(coarseComm);
label myAgglom = procAgglomMap[myProcID];
// Get all processors agglomerating to the same coarse processor
agglomProcIDs = findIndices(procAgglomMap, myAgglom);
// Make sure the master is the first element.
label index = findIndex
(
agglomProcIDs,
agglomToMaster[myAgglom]
);
Swap(agglomProcIDs[0], agglomProcIDs[index]);
}
//XXXX
scalarField allSource;
@ -463,7 +503,7 @@ void Foam::GAMGSolver::solveCoarsestLevel
cellOffsets.gather
(
coarseComm,
procIDs,
agglomProcIDs,
coarsestSource,
allSource
);
@ -471,7 +511,7 @@ void Foam::GAMGSolver::solveCoarsestLevel
scalarField allCorrField;
solverPerformance coarseSolverPerf;
if (Pstream::myProcNo(coarseComm) == procIDs[0])
if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
{
const lduMatrix& allMatrix = allMatrixPtr_();
@ -481,7 +521,8 @@ void Foam::GAMGSolver::solveCoarsestLevel
label solveComm = allMatrix.mesh().comm();
label oldWarn = UPstream::warnComm;
UPstream::warnComm = solveComm;
Pout<< "** Master:Solving on comm:" << solveComm << endl;
Pout<< "** Master:Solving on comm:" << solveComm
<< " with procs:" << UPstream::procID(solveComm) << endl;
if (allMatrix.asymmetric())
{
@ -521,12 +562,14 @@ void Foam::GAMGSolver::solveCoarsestLevel
}
}
Pout<< "done master solve." << endl;
// Scatter to all processors
coarsestCorrField.setSize(coarsestSource.size());
cellOffsets.scatter
(
coarseComm,
procIDs,
agglomProcIDs,
allCorrField,
coarsestCorrField
);

13
src/OpenFOAM/meshes/lduMesh/lduMesh.C
View File

@ -111,6 +111,19 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const InfoProxy<lduMesh>& ip)
Pout<< " patch:" << i
<< " type:" << interfaces[i].type() << endl;
if (isA<processorLduInterface>(interfaces[i]))
{
const processorLduInterface& pi = refCast
<
const processorLduInterface
>(interfaces[i]);
Pout<< " myProcNo:" << pi.myProcNo()
<< " neighbProcNo:" << pi.neighbProcNo()
<< " comm:" << pi.comm()
<< endl;
}
forAll(faceCells, i)
{
Pout<< " " << i << " own:" << faceCells[i]

508
src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C
View File

@ -29,6 +29,31 @@ License
#include "labelPair.H"
#include "processorGAMGInterface.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
//- Less operator for pairs of <processor><index>
class procLess
{
const labelPairList& lst_;
public:
procLess(const labelPairList& lst)
:
lst_(lst)
{}
bool operator()(const label a, const label b)
{
return lst_[a].first() < lst_[b].first();
}
};
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::lduPrimitiveMesh::checkUpperTriangular
@ -135,34 +160,34 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
patchSchedule_(ps),
comm_(comm)
{
Pout<< "lduPrimitiveMesh :"
<< " nCells:" << nCells
<< " l:" << lowerAddr_.size()
<< " u:" << upperAddr_.size()
<< " pa:" << patchAddr_.size()
<< " interfaces:" << interfaces_.size()
<< " comm:" << comm_
<< endl;
forAll(interfaces_, i)
{
if (interfaces_.set(i))
{
if (isA<processorLduInterface>(interfaces_[i]))
{
const processorLduInterface& pi = refCast
<
const processorLduInterface
>(interfaces_[i]);
Pout<< " patch:" << i
<< " size:" << patchAddr_[i].size()
<< " myProcNo:" << pi.myProcNo()
<< " neighbProcNo:" << pi.neighbProcNo()
<< " comm:" << pi.comm()
<< endl;
}
}
}
//Pout<< "lduPrimitiveMesh :"
// << " nCells:" << nCells
// << " l:" << lowerAddr_.size()
// << " u:" << upperAddr_.size()
// << " pa:" << patchAddr_.size()
// << " interfaces:" << interfaces_.size()
// << " comm:" << comm_
// << endl;
//forAll(interfaces_, i)
//{
// if (interfaces_.set(i))
// {
// if (isA<processorLduInterface>(interfaces_[i]))
// {
// const processorLduInterface& pi = refCast
// <
// const processorLduInterface
// >(interfaces_[i]);
//
// Pout<< " patch:" << i
// << " size:" << patchAddr_[i].size()
// << " myProcNo:" << pi.myProcNo()
// << " neighbProcNo:" << pi.neighbProcNo()
// << " comm:" << pi.comm()
// << endl;
// }
// }
//}
}
@ -186,40 +211,42 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
patchSchedule_(ps),
comm_(comm)
{
Pout<< "lduPrimitiveMesh :"
<< " nCells:" << nCells
<< " l:" << lowerAddr_.size()
<< " u:" << upperAddr_.size()
<< " pa:" << patchAddr_.size()
<< " interfaces:" << interfaces_.size()
<< " comm:" << comm_
<< endl;
forAll(interfaces_, i)
{
if (interfaces_.set(i))
{
if (isA<processorLduInterface>(interfaces_[i]))
{
const processorLduInterface& pi = refCast
<
const processorLduInterface
>(interfaces_[i]);
Pout<< " patch:" << i
<< " size:" << patchAddr_[i].size()
<< " myProcNo:" << pi.myProcNo()
<< " neighbProcNo:" << pi.neighbProcNo()
<< " comm:" << pi.comm()
<< endl;
}
}
}
//Pout<< "lduPrimitiveMesh :"
// << " nCells:" << nCells
// << " l:" << lowerAddr_.size()
// << " u:" << upperAddr_.size()
// << " pa:" << patchAddr_.size()
// << " interfaces:" << interfaces_.size()
// << " comm:" << comm_
// << endl;
//forAll(interfaces_, i)
//{
// if (interfaces_.set(i))
// {
// if (isA<processorLduInterface>(interfaces_[i]))
// {
// const processorLduInterface& pi = refCast
// <
// const processorLduInterface
// >(interfaces_[i]);
//
// Pout<< " patch:" << i
// << " size:" << patchAddr_[i].size()
// << " myProcNo:" << pi.myProcNo()
// << " neighbProcNo:" << pi.neighbProcNo()
// << " comm:" << pi.comm()
// << endl;
// }
// }
//}
}
Foam::lduPrimitiveMesh::lduPrimitiveMesh
(
const label comm,
const labelList& procAgglomMap,
const labelList& procIDs,
const lduMesh& myMesh,
const PtrList<lduMesh>& otherMeshes,
@ -238,21 +265,6 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
patchSchedule_(0),
comm_(comm)
{
// Sanity check.
for (label i = 1; i < procIDs.size(); i++)
{
if (procIDs[i] <= procIDs[i-1])
{
FatalErrorIn
(
"lduPrimitiveMesh::lduPrimitiveMesh(..)"
) << "Processor " << procIDs[i] << " at index " << i
<< " should be higher numbered than its predecessor "
<< procIDs[i-1]
<< exit(FatalError);
}
}
const label currentComm = myMesh.comm();
forAll(otherMeshes, i)
@ -272,6 +284,30 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
const label nMeshes = otherMeshes.size()+1;
const label myAgglom = procAgglomMap[UPstream::myProcNo(currentComm)];
Pout<< "I am " << UPstream::myProcNo(currentComm)
<< " agglomerating into " << myAgglom
<< " as are " << findIndices(procAgglomMap, myAgglom)
<< endl;
forAll(procIDs, i)
{
if (procAgglomMap[procIDs[i]] != procAgglomMap[procIDs[0]])
{
FatalErrorIn
(
"lduPrimitiveMesh::lduPrimitiveMesh(..)"
) << "Processor " << procIDs[i]
<< " agglomerates to " << procAgglomMap[procIDs[i]]
<< " whereas other processors " << procIDs
<< " agglomerate to "
<< UIndirectList<label>(procAgglomMap, procIDs)
<< exit(FatalError);
}
}
// Cells get added in order.
cellOffsets.setSize(nMeshes+1);
cellOffsets[0] = 0;
@ -298,15 +334,13 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
// Count how faces get added. Interfaces inbetween get merged.
// Merged interfaces: map from two processors back to
// Merged interfaces: map from two coarse processors back to
// - procMeshes
// - interface in procMesh
EdgeMap<labelPairList> mergedMap
(
2*myMesh.interfaces().size()
);
// (estimate size from number of patches of mesh0)
EdgeMap<labelPairList> mergedMap(2*myMesh.interfaces().size());
// Unmerged interfaces: map from two processors back to
// Unmerged interfaces: map from two coarse processors back to
// - procMeshes
// - interface in procMesh
EdgeMap<labelPairList> unmergedMap(mergedMap.size());
@ -315,8 +349,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
boundaryFaceMap.setSize(nMeshes);
label nBoundaryFaces = 0;
label nInterfaces = 0;
label nOtherInterfaces = 0;
labelList nCoupledFaces(nMeshes, 0);
for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
@ -333,73 +366,46 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
{
if (interfaces.set(intI))
{
if (isA<processorLduInterface>(interfaces[intI]))
const lduInterface& ldui = interfaces[intI];
if (isA<processorLduInterface>(ldui))
{
const processorLduInterface& pldui =
refCast<const processorLduInterface>
(
interfaces[intI]
);
if (pldui.myProcNo() != procIDs[procMeshI])
refCast<const processorLduInterface>(ldui);
label agglom0 = procAgglomMap[pldui.myProcNo()];
label agglom1 = procAgglomMap[pldui.neighbProcNo()];
const edge procEdge(agglom0, agglom1);
if (agglom0 != myAgglom && agglom1 != myAgglom)
{
FatalErrorIn
(
"lduPrimitiveMesh::lduPrimitiveMesh(..)"
) << "proc:" << procIDs[procMeshI]
<< " myProcNo:" << pldui.myProcNo()
<< abort(FatalError);
FatalErrorIn("lduPrimitiveMesh::lduPrimitiveMesh(..)")
<< "At mesh from processor " << procIDs[procMeshI]
<< " have interface " << intI
<< " with myProcNo:" << pldui.myProcNo()
<< " with neighbProcNo:" << pldui.neighbProcNo()
<< exit(FatalError);
}
const edge procEdge
(
min(pldui.myProcNo(), pldui.neighbProcNo()),
max(pldui.myProcNo(), pldui.neighbProcNo())
);
label index = findIndex(procIDs, pldui.neighbProcNo());
if (index == -1)
{
// Still external interface
Pout<< "external interface: myProcNo:"
<< pldui.myProcNo()
<< " nbr:" << pldui.neighbProcNo()
<< " size:" << interfaces[intI].faceCells().size()
<< endl;
nBoundaryFaces += interfaces[intI].faceCells().size();
EdgeMap<labelPairList>::iterator iter =
unmergedMap.find(procEdge);
if (iter != unmergedMap.end())
{
iter().append(labelPair(procMeshI, intI));
}
else
{
unmergedMap.insert
(
procEdge,
labelPairList(1, labelPair(procMeshI, intI))
);
}
nInterfaces++;
}
else
else if (agglom0 == myAgglom && agglom1 == myAgglom)
{
// Merged interface
Pout<< "merged interface: myProcNo:" << pldui.myProcNo()
Pout<< "merged interface: myProcNo:"
<< pldui.myProcNo()
<< " nbr:" << pldui.neighbProcNo()
<< " size:" << interfaces[intI].faceCells().size()
<< " size:" << ldui.faceCells().size()
<< endl;
if (pldui.myProcNo() < pldui.neighbProcNo())
label nbrProcMeshI = findIndex
(
procIDs,
pldui.neighbProcNo()
);
if (procMeshI < nbrProcMeshI)
{
nCoupledFaces[procMeshI] +=
interfaces[intI].faceCells().size();
// I am 'master' since get lowest numbered cells
nCoupledFaces[procMeshI] += ldui.faceCells().size();
}
EdgeMap<labelPairList>::iterator iter =
@ -418,6 +424,30 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
);
}
}
else
{
Pout<< "external interface: myProcNo:"
<< pldui.myProcNo()
<< " nbr:" << pldui.neighbProcNo()
<< " size:" << ldui.faceCells().size()
<< endl;
EdgeMap<labelPairList>::iterator iter =
unmergedMap.find(procEdge);
if (iter != unmergedMap.end())
{
iter().append(labelPair(procMeshI, intI));
}
else
{
unmergedMap.insert
(
procEdge,
labelPairList(1, labelPair(procMeshI, intI))
);
}
}
}
else
{
@ -428,8 +458,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
<< " of unhandled type " << interfaces[intI].type()
<< exit(FatalError);
nBoundaryFaces += interfaces[intI].faceCells().size();
nInterfaces++;
nOtherInterfaces++;
}
}
}
@ -437,11 +466,12 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
//if (debug)
{
Pout<< "Remaining interfaces:" << endl;
forAllConstIter(EdgeMap<labelPairList>, unmergedMap, iter)
{
Pout<< " procEdge:" << iter.key() << endl;
Pout<< " agglom procEdge:" << iter.key() << endl;
const labelPairList& elems = iter();
forAll(elems, i)
{
@ -449,6 +479,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
label interfaceI = elems[i][1];
const lduInterfacePtrsList interfaces =
mesh(myMesh, otherMeshes, procMeshI).interfaces();
const processorLduInterface& pldui =
refCast<const processorLduInterface>
(
@ -464,12 +495,14 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
Pout<< endl;
}
}
//if (debug)
{
Pout<< "Merged interfaces:" << endl;
forAllConstIter(EdgeMap<labelPairList>, mergedMap, iter)
{
Pout<< " procEdge:" << iter.key() << endl;
Pout<< " agglom procEdge:" << iter.key() << endl;
const labelPairList& elems = iter();
forAll(elems, i)
{
label procMeshI = elems[i][0];
@ -548,46 +581,46 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
{
if (interfaces.set(intI))
{
if (isA<processorLduInterface>(interfaces[intI]))
const lduInterface& ldui = interfaces[intI];
if (isA<processorLduInterface>(ldui))
{
const processorLduInterface& pldui =
refCast<const processorLduInterface>
(
interfaces[intI]
);
refCast<const processorLduInterface>(ldui);
// Look up corresponding interfaces
label myP = pldui.myProcNo();
label nbrP = pldui.neighbProcNo();
label nbrProcMeshI = findIndex(procIDs, nbrP);
if (myP < nbrP)
if (procMeshI < nbrProcMeshI)
{
// I am 'master' since my cell numbers will be lower
// since cells get added in procMeshI order.
label agglom0 = procAgglomMap[myP];
label agglom1 = procAgglomMap[nbrP];
EdgeMap<labelPairList>::const_iterator fnd =
mergedMap.find(edge(myP, nbrP));
mergedMap.find(edge(agglom0, agglom1));
if (fnd != mergedMap.end())
{
const labelPairList& elems = fnd();
// Find nbrP in elems
label nbrProcMeshI = -1;
label nbrIntI = -1;
if (procIDs[elems[0][0]] == nbrP)
forAll(elems, i)
{
nbrProcMeshI = elems[0][0];
nbrIntI = elems[0][1];
}
else
{
nbrProcMeshI = elems[1][0];
nbrIntI = elems[1][1];
if (elems[i][0] == nbrProcMeshI)
{
nbrIntI = elems[i][1];
break;
}
}
if
(
elems.size() != 2
|| procIDs[nbrProcMeshI] != nbrP
)
if (nbrIntI == -1)
{
FatalErrorIn
(
@ -596,17 +629,16 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
}
const lduInterfacePtrsList nbrInterfaces =
mesh
(
myMesh,
otherMeshes,
nbrProcMeshI
).interfaces();
const lduInterfacePtrsList nbrInterfaces = mesh
(
myMesh,
otherMeshes,
nbrProcMeshI
).interfaces();
const labelUList& faceCells =
interfaces[intI].faceCells();
ldui.faceCells();
const labelUList& nbrFaceCells =
nbrInterfaces[nbrIntI].faceCells();
@ -615,8 +647,8 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
labelList& nbrBfMap =
boundaryFaceMap[nbrProcMeshI][nbrIntI];
bfMap.setSize(faceCells.size(), -1);
nbrBfMap.setSize(faceCells.size(), -1);
bfMap.setSize(faceCells.size());
nbrBfMap.setSize(faceCells.size());
forAll(faceCells, pfI)
{
@ -639,57 +671,71 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
// Kept interfaces
// ~~~~~~~~~~~~~~~
interfaces_.setSize(nInterfaces);
interfaces_.setSize(unmergedMap.size() + nOtherInterfaces);
primitiveInterfaces_.setSize(interfaces_.size());
label allInterfaceI = 0;
forAllConstIter(EdgeMap<labelPairList>, unmergedMap, iter)
{
Pout<< "procEdge:" << iter.key() << endl;
Pout<< "agglom procEdge:" << iter.key() << endl;
const labelPairList& elems = iter();
// Sort processors in increasing order
labelList order(identity(elems.size()));
Foam::sort(order, procLess(elems));
// Count
label n = 0;
forAll(elems, i)
forAll(order, i)
{
label procMeshI = elems[i][0];
label interfaceI = elems[i][1];
const lduInterfacePtrsList interfaces =
mesh
(
myMesh,
otherMeshes,
procMeshI
).interfaces();
const labelPair& elem = elems[order[i]];
label procMeshI = elem[0];
label interfaceI = elem[1];
const lduInterfacePtrsList interfaces = mesh
(
myMesh,
otherMeshes,
procMeshI
).interfaces();
Pout<< " adding interface:" << " procMeshI:" << procMeshI
<< " interface:" << interfaceI
<< " size:" << interfaces[interfaceI].faceCells().size()
<< endl;
n += interfaces[interfaceI].faceCells().size();
}
// Size
labelField allFaceCells(n);
labelField allFaceRestrictAddressing(n);
n = 0;
forAll(elems, i)
// Fill
forAll(order, i)
{
label procMeshI = elems[i][0];
label interfaceI = elems[i][1];
const lduInterfacePtrsList interfaces =
mesh
(
myMesh,
otherMeshes,
procMeshI
).interfaces();
const labelPair& elem = elems[order[i]];
label procMeshI = elem[0];
label interfaceI = elem[1];
const lduInterfacePtrsList interfaces = mesh
(
myMesh,
otherMeshes,
procMeshI
).interfaces();
boundaryMap[procMeshI][interfaceI] = allInterfaceI;
labelList& bfMap = boundaryFaceMap[procMeshI][interfaceI];
const labelUList& l = interfaces[interfaceI].faceCells();
bfMap.setSize(l.size());
forAll(l, faceI)
{
allFaceCells[n] = cellOffsets[procMeshI]+l[faceI];
allFaceRestrictAddressing[n] = n;
bfMap[faceI] = faceI;
bfMap[faceI] = n;
n++;
}
}
@ -697,42 +743,38 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
// Find out local and remote processor in new communicator
label myProcNo = -1;
label neighbProcNo = -1;
if (findIndex(procIDs, iter.key()[0]) != -1)
// See what the two processors map onto
if (iter.key()[0] == myAgglom)
{
myProcNo = UPstream::procNo
(
comm_,
currentComm,
iter.key()[0]
);
neighbProcNo = UPstream::procNo
(
comm_,
currentComm,
iter.key()[1]
);
if (iter.key()[1] == myAgglom)
{
FatalErrorIn
(
"lduPrimitiveMesh::lduPrimitiveMesh(..)"
) << "problem procEdge:" << iter.key()
<< exit(FatalError);
}
neighbProcNo = iter.key()[1];
}
else
{
myProcNo = UPstream::procNo
(
comm_,
currentComm,
iter.key()[1]
);
neighbProcNo = UPstream::procNo
(
comm_,
currentComm,
iter.key()[0]
);
if (iter.key()[1] != myAgglom)
{
FatalErrorIn
(
"lduPrimitiveMesh::lduPrimitiveMesh(..)"
) << "problem procEdge:" << iter.key()
<< exit(FatalError);
}
neighbProcNo = iter.key()[0];
}
interfaces_.set
primitiveInterfaces_.set
(
allInterfaceI,
new processorGAMGInterface
@ -742,31 +784,39 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
allFaceCells,
allFaceRestrictAddressing,
comm_,
myProcNo,
myAgglom,
neighbProcNo,
tensorField(), // forwardT
Pstream::msgType() // tag
)
);
interfaces_.set(allInterfaceI, &primitiveInterfaces_[allInterfaceI]);
Pout<< "Created " << interfaces_[allInterfaceI].type()
<< " interface at " << allInterfaceI
<< " comm:" << comm_
<< " myProcNo:" << myAgglom
<< " neighbProcNo:" << neighbProcNo
<< " nFaces:" << allFaceCells.size()
<< endl;
allInterfaceI++;
}
// Extract faceCells from interfaces_
labelListList patchAddr_(interfaces_.size());
patchAddr_.setSize(interfaces_.size());
forAll(interfaces_, coarseIntI)
{
if (interfaces_.set(coarseIntI))
{
patchAddr_[coarseIntI] =
interfaces_[coarseIntI].faceCells();
patchAddr_[coarseIntI] = interfaces_[coarseIntI].faceCells();
}
}
patchSchedule_ = nonBlockingSchedule<processorGAMGInterface>(interfaces_);
Pout<< "lowerAddr_:" << lowerAddr_ << endl;
Pout<< "upperAddr_:" << upperAddr_ << endl;
checkUpperTriangular(cellOffsets.last(), lowerAddr_, upperAddr_);
}

31
src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H
View File

@ -67,6 +67,9 @@ class lduPrimitiveMesh
// with only those pointing to interfaces being set
lduInterfacePtrsList interfaces_;
//- Concrete interfaces if above interfaces were created locally
PtrList<lduInterface> primitiveInterfaces_;
//- Patch field evaluation schedule
lduSchedule patchSchedule_;
@ -127,25 +130,21 @@ public:
);
//- Construct by combining multiple meshes. The meshes come from
// processors procIDs but these are only used to detect the interfaces
// that are to be merged.
//
// Sets mapping information:
// cellOffsets : per mesh the starting cell
// faceMap : per mesh, per internal face, the combined face
// boundaryMap : per mesh, per boundary, the combined boundary
// or -1 for merged boundaries
// boundaryFaceMap : per mesh, per boundary face:
// kept boundaries: the new boundary face
// merged boundaries: the new internal face (on the owner) or
// -faceI-1 (on the neighbour)
//
//// TBD: hardcoded for processorGAMGInterface. Should use some
//// 'clone' or New functionality instead.
// processors procIDs:
// procIDs[0] : local processor (myMesh)
// procIDs[i] : processor where otherMeshes[i-1] comes from
// procAgglomMap : for every processor which processor it agglomerates
// onto. The new processor numbers are in compact
// numbering (so ranks in communicator comm), i.e.
// similar to cell-restrict-addressing.
// We need this information to be able to map
// inter-processor interfaces
lduPrimitiveMesh
(
const label comm,
const labelList& procIDs,
const labelList& procAgglomMap,
const labelList& procIDs, // procIDs[0] = myMesh
const lduMesh& myMesh,
const PtrList<lduMesh>& otherMeshes,

50
src/Pstream/mpi/UPstream.C
View File

@ -359,10 +359,15 @@ void Foam::UPstream::allocatePstreamCommunicator
}
else
{
//std::cout
// << "MPI : Allocating new communicator at index " << index
// << " from parent " << parentIndex
// << std::endl;
std::cout
<< "MPI : Allocating new communicator at index " << index
<< " from parent " << parentIndex
<< std::endl;
for (label i=0; i < procIDs_[index].size(); i++)
{
std::cout<< " " << i << " rank:" << procIDs_[index][i]
<< std::endl;
}
// Create new group
MPI_Group_incl
@ -395,19 +400,38 @@ void Foam::UPstream::allocatePstreamCommunicator
}
else
{
//std::cout
// << "MPI : New comm "
// << long(PstreamGlobals::MPICommunicators_[index])
// << std::endl;
MPI_Comm_rank
std::cout
<< "MPI : New comm "
<< long(PstreamGlobals::MPICommunicators_[index])
<< std::endl;
if
(
PstreamGlobals::MPICommunicators_[index],
&myProcNo_[index]
);
MPI_Comm_rank
(
PstreamGlobals::MPICommunicators_[index],
&myProcNo_[index]
)
)
{
FatalErrorIn
(
"UPstream::allocatePstreamCommunicator\n"
"(\n"
" const label,\n"
" const labelList&\n"
")\n"
) << "Problem :"
<< " when allocating communicator at " << index
<< " from ranks " << procIDs_[index]
<< " of parent " << parentIndex
<< " cannot find my own rank"
<< Foam::exit(FatalError);
}
}
}
//std::cout<< "MPI : I am rank " << myProcNo_[index] << std::endl;
std::cout<< "MPI : in communicator " << index
<< " I am rank " << myProcNo_[index] << std::endl;
}

11
src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.H
View File

@ -146,6 +146,17 @@ public:
{
return fineCyclicAMIInterface_.reverseT();
}
// I/O
//- Write to stream
virtual void write(Ostream&) const
{
//TBD. How to serialise the AMI such that we can stream
// cyclicAMIGAMGInterface.
notImplemented("cyclicAMIGAMGInterface::write(Ostream&) const");
}
};

12
src/meshTools/regionCoupled/GAMG/interfaces/regionCoupledGAMGInterface/regionCoupledGAMGInterface.H
View File

@ -85,6 +85,18 @@ public:
virtual ~regionCoupledGAMGInterface();
// I/O
//- Write to stream
virtual void write(Ostream&) const
{
//TBD. How to serialise the AMI such that we can stream
// regionCoupledGAMGInterface.
notImplemented
(
"regionCoupledGAMGInterface::write(Ostream&) const"
);
}
};

14
src/meshTools/regionCoupled/GAMG/interfaces/regionCoupledGAMGInterface/regionCoupledWallGAMGInterface.H
View File

@ -86,6 +86,20 @@ public:
//- Destructor
virtual ~regionCoupledWallGAMGInterface();
// I/O
//- Write to stream
virtual void write(Ostream&) const
{
//TBD. How to serialise the AMI such that we can stream
// regionCoupledWallGAMGInterface.
notImplemented
(
"regionCoupledWallGAMGInterface::write(Ostream&) const"
);
}
};