From 549d26f3077da311d9dffeafdc6b4412bdd99d81 Mon Sep 17 00:00:00 2001 From: mattijs <mattijs> Date: Fri, 1 Mar 2013 06:03:10 +0000 Subject: [PATCH] ENH: GAMGAgglomeration: moved lduMesh agglomeration into GAMGAgglomeration --- .../GAMGAgglomeration/GAMGAgglomeration.C | 97 ++++++ .../GAMGAgglomeration/GAMGAgglomeration.H | 56 ++++ .../lduMatrix/solvers/GAMG/GAMGSolver.C | 297 ++++++++--------- .../lduMatrix/solvers/GAMG/GAMGSolver.H | 37 +-- .../lduMatrix/solvers/GAMG/GAMGSolverSolve.C | 124 ++----- src/OpenFOAM/meshes/lduMesh/lduMesh.C | 30 ++ .../meshes/lduMesh/lduPrimitiveMesh.C | 306 +++++++----------- .../meshes/lduMesh/lduPrimitiveMesh.H | 35 +- 8 files changed, 488 insertions(+), 494 deletions(-) diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C index a724e513f7a..49c7aad9331 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C @@ -255,4 +255,101 @@ const Foam::lduInterfacePtrsList& Foam::GAMGAgglomeration::interfaceLevel } +void Foam::GAMGAgglomeration::gatherMeshes +( + const labelList& procIDs, + const label allMeshComm +) const +{ + const lduMesh& coarsestMesh = meshLevels_.last(); + + label coarseComm = coarsestMesh.comm(); + + label oldWarn = UPstream::warnComm; + UPstream::warnComm = coarseComm; + + Pout<< "GAMGAgglomeration : gathering coarsestmesh (level=" + << meshLevels_.size()-1 + << ") using communicator " << coarseComm << endl; + + lduPrimitiveMesh::gather(coarsestMesh, procIDs, otherMeshes_); + + Pout<< "** Own Mesh " << coarsestMesh.info() << endl; + forAll(otherMeshes_, i) + { + Pout<< "** otherMesh " << i << " " + << otherMeshes_[i].info() + << endl; + } + Pout<< endl; + + if (Pstream::myProcNo(coarseComm) == procIDs[0]) + { + // Agglomerate all addressing + // ~~~~~~~~~~~~~~~~~~~~~~~~~~ + + allMeshPtr_.reset + ( + new lduPrimitiveMesh + ( + allMeshComm, + procIDs, + coarsestMesh, + otherMeshes_, + + cellOffsets_, + faceMap_, + boundaryMap_, + boundaryFaceMap_ + ) + ); + } + + UPstream::warnComm = oldWarn; +} + + +const Foam::lduPrimitiveMesh& Foam::GAMGAgglomeration::allMesh() const +{ + if (!allMeshPtr_.valid()) + { + FatalErrorIn("GAMGAgglomeration::allMesh() const") + << "Processor-agglomerated mesh not constructed. Did you call" + << " gatherMeshes?" + << exit(FatalError); + } + return allMeshPtr_(); +} + + +const Foam::labelList& Foam::GAMGAgglomeration::cellOffsets() const +{ + return cellOffsets_; +} + + +const Foam::labelListList& Foam::GAMGAgglomeration::faceMap() const +{ + return faceMap_; +} + + +const Foam::labelListList& Foam::GAMGAgglomeration::boundaryMap() const +{ + return boundaryMap_; +} + + +const Foam::labelListListList& Foam::GAMGAgglomeration::boundaryFaceMap() const +{ + return boundaryFaceMap_; +} + + +const Foam::PtrList<Foam::lduMesh>& Foam::GAMGAgglomeration::otherMeshes() const +{ + return otherMeshes_; +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H index 78e9aa536b4..7a89de8736b 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H @@ -99,6 +99,30 @@ protected: // Warning: Needs to be deleted explicitly. PtrList<lduInterfacePtrsList> interfaceLevels_; + + //- Processor agglomeration + + //- Combined coarsest meshes for all processors + mutable autoPtr<lduPrimitiveMesh> allMeshPtr_; + + //- Per remote processor the local equivalent. TBD: hack. Can be + // deleted once lduMatrices get received and merged into single + // one without constructing local ones + mutable PtrList<lduMesh> otherMeshes_; + + //- Mapping from processor to coarsestAllMesh cells + mutable labelList cellOffsets_; + + //- Mapping from processor to coarsestAllMesh face + mutable labelListList faceMap_; + + //- Mapping from processor to coarsestAllMesh boundary + mutable labelListList boundaryMap_; + + //- Mapping from processor to coarsestAllMesh boundary face + mutable labelListListList boundaryFaceMap_; + + //- Assemble coarse mesh addressing void agglomerateLduAddressing(const label fineLevelIndex); @@ -250,6 +274,38 @@ public: const Field<Type>& cf, const label coarseLevelIndex ) const; + + + // Procesor agglomeration + + //- Collect and combine processor meshes into allMesh. Pass in + // communicator for combined mesh. + void gatherMeshes + ( + const labelList& procIDs, + const label allMeshComm + ) const; + + //- Combined mesh + const lduPrimitiveMesh& allMesh() const; + + //- Mapping from processor to coarsestAllMesh cells + const labelList& cellOffsets() const; + + //- Mapping from processor to coarsestAllMesh face + const labelListList& faceMap() const; + + //- Mapping from processor to coarsestAllMesh boundary + const labelListList& boundaryMap() const; + + //- Mapping from processor to coarsestAllMesh boundary face + const labelListListList& boundaryFaceMap() const; + + //- Per remote processor the local equivalent. TBD: hack. Can be + // deleted once lduMatrices get received and merged into single + // one without constructing local ones + const PtrList<lduMesh>& otherMeshes() const; + }; diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C index b068d8e639c..38f51c8a011 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C @@ -45,10 +45,12 @@ namespace Foam // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // +// Gather matrices. +// Note: matrices get constructed with dummy mesh void Foam::GAMGSolver::gatherMatrices ( const labelList& procIDs, - const PtrList<lduMesh>& procMeshes, + const lduMesh& dummyMesh, const label meshComm, const lduMatrix& mat, @@ -59,16 +61,14 @@ void Foam::GAMGSolver::gatherMatrices PtrList<lduMatrix>& otherMats, PtrList<FieldField<Field, scalar> >& otherBouCoeffs, PtrList<FieldField<Field, scalar> >& otherIntCoeffs, - PtrList<lduInterfaceFieldPtrsList>& otherInterfaces + List<boolList>& otherTransforms, + List<List<int> >& otherRanks ) const { Pout<< "GAMGSolver::gatherMatrices :" << " collecting matrices on procs:" << procIDs << " using comm:" << meshComm << endl; - - //lduInterfacePtrsList interfaces(mesh.interfaces()); - if (Pstream::myProcNo(meshComm) == procIDs[0]) { // Master. @@ -76,85 +76,59 @@ void Foam::GAMGSolver::gatherMatrices Pout<< "GAMGSolver::gatherMatrices :" << " master:" << Pstream::myProcNo(meshComm) << endl; - for (label i=0; i < procMeshes.size(); i++) - { - if (meshComm != procMeshes[i].comm()) - { - FatalErrorIn("GAMGSolver::gatherMatrices()") - << "Inconsistent communicator :" - << "master processor:" << procIDs[0] - << " comm:" << meshComm - << " slave processor:" << procIDs[i] - << " comm:" << procMeshes[i].comm() - << abort(FatalError); - } - } - otherMats.setSize(procIDs.size()-1); otherBouCoeffs.setSize(procIDs.size()-1); otherIntCoeffs.setSize(procIDs.size()-1); - otherInterfaces.setSize(procIDs.size()-1); + otherTransforms.setSize(procIDs.size()-1); + otherRanks.setSize(procIDs.size()-1); - for (label i = 1; i < procIDs.size(); i++) + for (label procI = 1; procI < procIDs.size(); procI++) { - const lduMesh& procMesh = procMeshes[i]; - const lduInterfacePtrsList procInterfaces = procMesh.interfaces(); + label otherI = procI-1; IPstream fromSlave ( Pstream::scheduled, - procIDs[i], + procIDs[procI], 0, // bufSize Pstream::msgType(), meshComm ); - otherMats.set(i-1, new lduMatrix(procMesh, fromSlave)); + otherMats.set(otherI, new lduMatrix(dummyMesh, fromSlave)); // Receive number of/valid interfaces - boolList validTransforms(fromSlave); - List<int> validRanks(fromSlave); + boolList& procTransforms = otherTransforms[otherI]; + List<int>& procRanks = otherRanks[otherI]; + + fromSlave >> procTransforms; + fromSlave >> procRanks; // Size coefficients otherBouCoeffs.set ( - i-1, - new FieldField<Field, scalar>(validTransforms.size()) + otherI, + new FieldField<Field, scalar>(procRanks.size()) ); otherIntCoeffs.set ( - i-1, - new FieldField<Field, scalar>(validTransforms.size()) + otherI, + new FieldField<Field, scalar>(procRanks.size()) ); - otherInterfaces.set - ( - i-1, - new lduInterfaceFieldPtrsList(validTransforms.size()) - ); - - forAll(validTransforms, intI) + forAll(procRanks, intI) { - if (validTransforms[intI]) + if (procRanks[intI] != -1) { - const processorGAMGInterface& interface = - refCast<const processorGAMGInterface> - ( - procInterfaces[intI] - ); - - - otherBouCoeffs[i-1].set(intI, new scalarField(fromSlave)); - otherIntCoeffs[i-1].set(intI, new scalarField(fromSlave)); - otherInterfaces[i-1].set + otherBouCoeffs[otherI].set ( intI, - GAMGInterfaceField::New - ( - interface, //procInterfaces[intI], - validTransforms[intI], - validRanks[intI] - ).ptr() + new scalarField(fromSlave) + ); + otherIntCoeffs[otherI].set + ( + intI, + new scalarField(fromSlave) ); } } @@ -169,8 +143,8 @@ void Foam::GAMGSolver::gatherMatrices // Send to master // Count valid interfaces - boolList validTransforms(interfaceBouCoeffs.size(), false); - List<int> validRanks(interfaceBouCoeffs.size(), -1); + boolList procTransforms(interfaceBouCoeffs.size(), false); + List<int> procRanks(interfaceBouCoeffs.size(), -1); forAll(interfaces, intI) { if (interfaces.set(intI)) @@ -181,8 +155,8 @@ void Foam::GAMGSolver::gatherMatrices interfaces[intI] ); - validTransforms[intI] = interface.doTransform(); - validRanks[intI] = interface.rank(); + procTransforms[intI] = interface.doTransform(); + procRanks[intI] = interface.rank(); } } @@ -198,10 +172,10 @@ void Foam::GAMGSolver::gatherMatrices meshComm ); - toMaster << mat << validTransforms << validRanks; - forAll(validTransforms, intI) + toMaster << mat << procTransforms << procRanks; + forAll(procRanks, intI) { - if (validTransforms[intI]) + if (procRanks[intI] != -1) { toMaster << interfaceBouCoeffs[intI] @@ -310,23 +284,22 @@ Foam::GAMGSolver::GAMGSolver // All processors to master const List<int>& procIDs = UPstream::procID(coarseComm); - - // Gather all meshes onto procIDs[0] - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - // Note: move into GAMGAgglomeration - Pout<< "procIDs:" << procIDs << endl; - PtrList<lduMesh> procMeshes; - lduPrimitiveMesh::gather(coarsestMesh, procIDs, procMeshes); + // Allocate a communicator for single processor + label masterComm = UPstream::allocateCommunicator + ( + coarseComm, + labelList(1, 0) + ); + Pout<< "** Allocated communicator " << masterComm + << " for processor " << procIDs[0] + << " out of " << procIDs << endl; + - forAll(procMeshes, procI) - { - Pout<< "** procMesh " << procI << " " - << procMeshes[procI].info() - << endl; - } - Pout<< endl; + // Gather all meshes onto procIDs[0] + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + agglomeration_.gatherMeshes(procIDs, masterComm); // Gather all matrix coefficients onto procIDs[0] @@ -335,11 +308,12 @@ Foam::GAMGSolver::GAMGSolver PtrList<lduMatrix> otherMats; PtrList<FieldField<Field, scalar> > otherBouCoeffs; PtrList<FieldField<Field, scalar> > otherIntCoeffs; - PtrList<lduInterfaceFieldPtrsList> otherInterfaces; + List<boolList> otherTransforms; + List<List<int> > otherRanks; gatherMatrices ( procIDs, - procMeshes, + coarsestMesh, coarseComm, coarsestMatrix, @@ -350,55 +324,28 @@ Foam::GAMGSolver::GAMGSolver otherMats, otherBouCoeffs, otherIntCoeffs, - otherInterfaces - ); - - Pout<< "Own matrix:" << coarsestMatrix.info() << endl; - - forAll(otherMats, i) - { - Pout<< "** otherMats " << i << " " - << otherMats[i].info() - << endl; - } - Pout<< endl; - - - // Allocate a communicator for single processor - label masterComm = UPstream::allocateCommunicator - ( - coarseComm, - labelList(1, 0) + otherTransforms, + otherRanks ); - Pout<< "** Allocated communicator " << masterComm - << " for processor " << procIDs[0] - << " out of " << procIDs << endl; if (Pstream::myProcNo(coarseComm) == procIDs[0]) { - // Agglomerate all addressing - // ~~~~~~~~~~~~~~~~~~~~~~~~~~ + // Agglomerate all matrix + // ~~~~~~~~~~~~~~~~~~~~~~ - allMeshPtr_.reset - ( - new lduPrimitiveMesh - ( - masterComm, - procIDs, - procMeshes, + Pout<< "Own matrix:" << coarsestMatrix.info() << endl; - cellOffsets_, - faceMap_, - boundaryMap_, - boundaryFaceMap_ - ) - ); - const lduMesh& allMesh = allMeshPtr_(); + forAll(otherMats, i) + { + Pout<< "** otherMats " << i << " " + << otherMats[i].info() + << endl; + } + Pout<< endl; - // Agglomerate all matrix - // ~~~~~~~~~~~~~~~~~~~~~~ + const lduMesh& allMesh = agglomeration_.allMesh(); allMatrixPtr_.reset(new lduMatrix(allMesh)); lduMatrix& allMatrix = allMatrixPtr_(); @@ -420,7 +367,7 @@ Foam::GAMGSolver::GAMGSolver ( allDiag, otherMats[i].diag().size(), - cellOffsets_[i+1] + agglomeration_.cellOffsets()[i+1] ).assign ( otherMats[i].diag() @@ -430,23 +377,35 @@ Foam::GAMGSolver::GAMGSolver if (coarsestMatrix.hasLower()) { scalarField& allLower = allMatrix.lower(); - UIndirectList<scalar>(allLower, faceMap_[0]) = - coarsestMatrix.lower(); + UIndirectList<scalar> + ( + allLower, + agglomeration_.faceMap()[0] + ) = coarsestMatrix.lower(); forAll(otherMats, i) { - UIndirectList<scalar>(allLower, faceMap_[i+1]) = - otherMats[i].lower(); + UIndirectList<scalar> + ( + allLower, + agglomeration_.faceMap()[i+1] + ) = otherMats[i].lower(); } } if (coarsestMatrix.hasUpper()) { scalarField& allUpper = allMatrix.upper(); - UIndirectList<scalar>(allUpper, faceMap_[0]) = - coarsestMatrix.upper(); + UIndirectList<scalar> + ( + allUpper, + agglomeration_.faceMap()[0] + ) = coarsestMatrix.upper(); forAll(otherMats, i) { - UIndirectList<scalar>(allUpper, faceMap_[i+1]) = - otherMats[i].upper(); + UIndirectList<scalar> + ( + allUpper, + agglomeration_.faceMap()[i+1] + ) = otherMats[i].upper(); } } @@ -470,7 +429,7 @@ Foam::GAMGSolver::GAMGSolver } labelList nBounFaces(allMeshInterfaces.size()); - forAll(boundaryMap_, procI) + forAll(agglomeration_.boundaryMap(), procI) { const FieldField<Field, scalar>& procBouCoeffs ( @@ -484,15 +443,8 @@ Foam::GAMGSolver::GAMGSolver ? coarsestIntCoeffs : otherIntCoeffs[procI-1] ); - const lduInterfaceFieldPtrsList procInterfaces - ( - (procI == 0) - ? coarsestInterfaces - : otherInterfaces[procI-1] - ); - - const labelList& bMap = boundaryMap_[procI]; + const labelList& bMap = agglomeration_.boundaryMap()[procI]; forAll(bMap, procIntI) { label allIntI = bMap[procIntI]; @@ -506,11 +458,27 @@ Foam::GAMGSolver::GAMGSolver { // Construct lduInterfaceField - const processorGAMGInterfaceField& procInt = - refCast<const processorGAMGInterfaceField> - ( - procInterfaces[procIntI] - ); + 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 ( @@ -521,8 +489,8 @@ Foam::GAMGSolver::GAMGSolver ( allMeshInterfaces[allIntI] ), - procInt.doTransform(), - procInt.rank() + doTransform, + rank ).ptr() ); } @@ -535,8 +503,8 @@ Foam::GAMGSolver::GAMGSolver scalarField& allInt = allInterfaceIntCoeffs_[allIntI]; - const labelList& map = - boundaryFaceMap_[procI][procIntI]; + const labelList& map = agglomeration_. + boundaryFaceMap()[procI][procIntI]; const scalarField& procBou = procBouCoeffs[procIntI]; const scalarField& procInt = @@ -558,27 +526,15 @@ Foam::GAMGSolver::GAMGSolver { // Boundary has become internal face - // Problem: we don't know whether face is flipped - // or not. Either we interrogate the original - // mesh (requires procMeshes which we don't want) - // or we need to store this information somehow - // in the mapping (boundaryMap_, boundaryFaceMap_) - - - //Hack. See above. - //const labelList& fCells = - // procInterfaces[procI].faceCells(); - const labelList& fCells = - procMeshes[procI].lduAddr().patchAddr(procIntI); - - const labelList& map = - boundaryFaceMap_[procI][procIntI]; + const labelList& map = agglomeration_. + boundaryFaceMap()[procI][procIntI]; const labelList& allU = allMesh.lduAddr().upperAddr(); const labelList& allL = allMesh.lduAddr().lowerAddr(); - const label off = cellOffsets_[procI]; + const label off = agglomeration_. + cellOffsets()[procI]; const scalarField& procBou = procBouCoeffs[procIntI]; @@ -594,12 +550,12 @@ Foam::GAMGSolver::GAMGSolver if (coarsestMatrix.hasUpper()) { allMatrix.upper()[allFaceI] = - procBou[i]; + -procBou[i]; } if (coarsestMatrix.hasLower()) { allMatrix.lower()[allFaceI] = - procInt[i]; + -procInt[i]; } } else @@ -609,12 +565,12 @@ Foam::GAMGSolver::GAMGSolver if (coarsestMatrix.hasUpper()) { allMatrix.upper()[allFaceI] = - procInt[i]; + -procInt[i]; } if (coarsestMatrix.hasLower()) { allMatrix.lower()[allFaceI] = - procBou[i]; + -procBou[i]; } } @@ -626,6 +582,15 @@ Foam::GAMGSolver::GAMGSolver ? map[i] : -map[i]-1 ); + + const labelList& fCells = + lduPrimitiveMesh::mesh + ( + coarsestMesh, + agglomeration_.otherMeshes(), + procI + ).lduAddr().patchAddr(procIntI); + label allCellI = off + fCells[i]; if @@ -649,6 +614,8 @@ Foam::GAMGSolver::GAMGSolver } } } + + Pout<< "** Assembled allMatrix:" << allMatrix.info() << endl; } UPstream::warnComm = oldWarn; diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H index e88ac38ffa8..1736bf474ac 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H @@ -131,20 +131,6 @@ class GAMGSolver //- LU decompsed coarsest matrix autoPtr<LUscalarMatrix> coarsestLUMatrixPtr_; - - - // Processor addressing agglomeration - to be moved into - // GAMGAgglomeration - - //- Combined coarsest meshes for all processors - autoPtr<lduPrimitiveMesh> allMeshPtr_; - - //- Mapping from processor to coarsestAllMesh - labelList cellOffsets_; - labelListList faceMap_; - labelListList boundaryMap_; - labelListListList boundaryFaceMap_; - // Processor matrix agglomeration //- Combined coarsest level matrix for all processors @@ -187,7 +173,8 @@ class GAMGSolver void gatherMatrices ( const labelList& procIDs, - const PtrList<lduMesh>& procMeshes, + const lduMesh& myMesh, + const PtrList<lduMesh>& otherMeshes, const label meshComm, const lduMatrix& mat, @@ -200,14 +187,22 @@ class GAMGSolver PtrList<FieldField<Field, scalar> >& otherIntCoeffs, PtrList<lduInterfaceFieldPtrsList>& otherInterfaces ) const; - - //- Gather field in allMatrix form - void gatherField + void gatherMatrices ( - const label meshComm, const labelList& procIDs, - const scalarField& coarsestSource, - scalarField& allSource + 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; //- Interpolate the correction after injected prolongation diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C index f0ef2813778..735fe01dce0 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C @@ -27,71 +27,10 @@ License #include "ICCG.H" #include "BICCG.H" #include "SubField.H" - +#include "globalIndex.H" // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::GAMGSolver::gatherField -( - const label meshComm, - const labelList& procIDs, - const scalarField& coarsestSource, - scalarField& allSource -) const -{ - if (Pstream::myProcNo(meshComm) == procIDs[0]) - { - // Master. - const lduMesh& allMesh = allMeshPtr_(); - - allSource.setSize(allMesh.lduAddr().size()); - - SubList<scalar> - ( - allSource, - coarsestSource.size(), - 0 - ).assign(coarsestSource); - - for (label i = 1; i < procIDs.size(); i++) - { - IPstream fromSlave - ( - Pstream::scheduled, - procIDs[i], - 0, // bufSize - Pstream::msgType(), - meshComm - ); - - SubList<scalar> slaveSlots - ( - allSource, - cellOffsets_[i+1]-cellOffsets_[i], - cellOffsets_[i] - ); - - UList<scalar>& l = slaveSlots; - - fromSlave >> l; - } - } - else - { - // Send to master - OPstream toMaster - ( - Pstream::scheduled, - procIDs[0], - 0, - Pstream::msgType(), - meshComm - ); - toMaster << coarsestSource; - } -} - - Foam::solverPerformance Foam::GAMGSolver::solve ( scalarField& psi, @@ -519,7 +458,9 @@ void Foam::GAMGSolver::solveCoarsestLevel const List<int>& procIDs = UPstream::procID(coarseComm); scalarField allSource; - gatherField + + const globalIndex cellOffsets(agglomeration_.cellOffsets()); + cellOffsets.gather ( coarseComm, procIDs, @@ -578,55 +519,25 @@ void Foam::GAMGSolver::solveCoarsestLevel } UPstream::warnComm = oldWarn; } - - - // Distribute correction - coarsestCorrField = SubList<scalar> - ( - allCorrField, - coarsestSource.size(), - 0 - ); - for (label i=1; i < procIDs.size(); i++) - { - Pout<< "** master sending to " << procIDs[i] << endl; - - OPstream toSlave - ( - Pstream::scheduled, - procIDs[i], - 0, // bufSize - Pstream::msgType(), - coarseComm - ); - - toSlave << SubList<scalar> - ( - allCorrField, - cellOffsets_[i+1]-cellOffsets_[i], - cellOffsets_[i] - ); - } - } - else - { - Pout<< "** slave receiving from " << procIDs[0] << endl; - IPstream fromMaster - ( - Pstream::scheduled, - procIDs[0], - 0, // bufSize - Pstream::msgType(), - coarseComm - ); - fromMaster >> coarsestCorrField; } + // Scatter to all processors + coarsestCorrField.setSize(coarsestSource.size()); + cellOffsets.scatter + ( + coarseComm, + procIDs, + allCorrField, + coarsestCorrField + ); if (debug >= 2) { coarseSolverPerf.print(Info(coarseComm)); } + + Pout<< "procAgglom: coarsestSource :" << coarsestSource << endl; + Pout<< "procAgglom: coarsestCorrField:" << coarsestCorrField << endl; } else { @@ -672,6 +583,9 @@ void Foam::GAMGSolver::solveCoarsestLevel { coarseSolverPerf.print(Info(coarseComm)); } + + Pout<< "GC: coarsestSource :" << coarsestSource << endl; + Pout<< "GC: coarsestCorrField:" << coarsestCorrField << endl; } UPstream::warnComm = oldWarn; diff --git a/src/OpenFOAM/meshes/lduMesh/lduMesh.C b/src/OpenFOAM/meshes/lduMesh/lduMesh.C index 803cd607d60..4435d810744 100644 --- a/src/OpenFOAM/meshes/lduMesh/lduMesh.C +++ b/src/OpenFOAM/meshes/lduMesh/lduMesh.C @@ -91,6 +91,36 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const InfoProxy<lduMesh>& ip) } } + + // Print actual contents + { + const labelList& l = addr.lowerAddr(); + const labelList& u = addr.upperAddr(); + forAll(l, faceI) + { + Pout<< " face:" << faceI << " l:" << l[faceI] + << " u:" << u[faceI] << endl; + } + forAll(interfaces, i) + { + if (interfaces.set(i)) + { + const labelUList& faceCells = addr.patchAddr(i); + if (faceCells.size()) + { + Pout<< " patch:" << i + << " type:" << interfaces[i].type() << endl; + + forAll(faceCells, i) + { + Pout<< " " << i << " own:" << faceCells[i] + << endl; + } + } + } + } + } + os.check("Ostream& operator<<(Ostream&, const lduMesh&"); return os; diff --git a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C index 43a17578322..8eb47828a69 100644 --- a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C +++ b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C @@ -217,22 +217,12 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh } -//Foam::lduPrimitiveMesh::lduPrimitiveMesh(Istream& is) -//: -// lduAddressing(readLabel(is)), -// lowerAddr_(is), -// upperAddr_(is), -// patchAddr_(is), -// interfaces_(is), -// patchSchedule_(is), -// comm_(readLabel(is)) -//{} - Foam::lduPrimitiveMesh::lduPrimitiveMesh ( const label comm, const labelList& procIDs, - const PtrList<lduMesh>& procMeshes, + const lduMesh& myMesh, + const PtrList<lduMesh>& otherMeshes, labelList& cellOffsets, labelListList& faceMap, @@ -240,7 +230,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh labelListListList& boundaryFaceMap ) : - lduAddressing(size(procMeshes)), + lduAddressing(myMesh.lduAddr().size() + size(otherMeshes)), lowerAddr_(0), upperAddr_(0), patchAddr_(0), @@ -261,40 +251,49 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh << procIDs[i-1] << exit(FatalError); } - if (procMeshes[i].comm() != procMeshes[i-1].comm()) + } + + const label currentComm = myMesh.comm(); + + forAll(otherMeshes, i) + { + if (otherMeshes[i].comm() != currentComm) { FatalErrorIn ( "lduPrimitiveMesh::lduPrimitiveMesh(..)" - ) << "Communicator " << procMeshes[i].comm() << " at index " << i + ) << "Communicator " << otherMeshes[i].comm() + << " at index " << i << " differs from that of predecessor " - << procMeshes[i-1].comm() + << currentComm << exit(FatalError); } } - const label currentComm = procMeshes[0].comm(); + const label nMeshes = otherMeshes.size()+1; // Cells get added in order. - cellOffsets.setSize(procMeshes.size()+1); + cellOffsets.setSize(nMeshes+1); cellOffsets[0] = 0; - forAll(procMeshes, procMeshI) + for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++) { - const lduMesh& mesh = procMeshes[procMeshI]; + const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI); + cellOffsets[procMeshI+1] = cellOffsets[procMeshI] - + mesh.lduAddr().size(); + + procMesh.lduAddr().size(); } // Faces initially get added in order, sorted later - labelList internalFaceOffsets(procMeshes.size()+1); + labelList internalFaceOffsets(nMeshes+1); internalFaceOffsets[0] = 0; - forAll(procMeshes, procMeshI) + for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++) { - const lduMesh& mesh = procMeshes[procMeshI]; + const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI); + internalFaceOffsets[procMeshI+1] = internalFaceOffsets[procMeshI] - + mesh.lduAddr().lowerAddr().size(); + + procMesh.lduAddr().lowerAddr().size(); } // Count how faces get added. Interfaces inbetween get merged. @@ -304,7 +303,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh // - interface in procMesh EdgeMap<labelPairList> mergedMap ( - 2*procMeshes[0].interfaces().size() + 2*myMesh.interfaces().size() ); // Unmerged interfaces: map from two processors back to @@ -312,18 +311,18 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh // - interface in procMesh EdgeMap<labelPairList> unmergedMap(mergedMap.size()); - boundaryMap.setSize(procMeshes.size()); - boundaryFaceMap.setSize(procMeshes.size()); + boundaryMap.setSize(nMeshes); + boundaryFaceMap.setSize(nMeshes); label nBoundaryFaces = 0; label nInterfaces = 0; - labelList nCoupledFaces(procMeshes.size(), 0); + labelList nCoupledFaces(nMeshes, 0); - forAll(procMeshes, procMeshI) + for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++) { const lduInterfacePtrsList interfaces = - procMeshes[procMeshI].interfaces(); + mesh(myMesh, otherMeshes, procMeshI).interfaces(); // Inialise all boundaries as merged boundaryMap[procMeshI].setSize(interfaces.size(), -1); @@ -449,7 +448,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh label procMeshI = elems[i][0]; label interfaceI = elems[i][1]; const lduInterfacePtrsList interfaces = - procMeshes[procMeshI].interfaces(); + mesh(myMesh, otherMeshes, procMeshI).interfaces(); const processorLduInterface& pldui = refCast<const processorLduInterface> ( @@ -476,7 +475,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh label procMeshI = elems[i][0]; label interfaceI = elems[i][1]; const lduInterfacePtrsList interfaces = - procMeshes[procMeshI].interfaces(); + mesh(myMesh, otherMeshes, procMeshI).interfaces(); const processorLduInterface& pldui = refCast<const processorLduInterface> ( @@ -495,12 +494,13 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh // Adapt faceOffsets for internal interfaces - labelList faceOffsets(procMeshes.size()+1); + labelList faceOffsets(nMeshes+1); faceOffsets[0] = 0; - faceMap.setSize(procMeshes.size()); - forAll(procMeshes, procMeshI) + faceMap.setSize(nMeshes); + for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++) { - label nInternal = procMeshes[procMeshI].lduAddr().lowerAddr().size(); + const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI); + label nInternal = procMesh.lduAddr().lowerAddr().size(); faceOffsets[procMeshI+1] = faceOffsets[procMeshI] @@ -524,10 +524,12 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh // Old internal faces and resolved coupled interfaces // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - forAll(procMeshes, procMeshI) + for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++) { - const labelUList& l = procMeshes[procMeshI].lduAddr().lowerAddr(); - const labelUList& u = procMeshes[procMeshI].lduAddr().upperAddr(); + const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI); + + const labelUList& l = procMesh.lduAddr().lowerAddr(); + const labelUList& u = procMesh.lduAddr().upperAddr(); // Add internal faces label allFaceI = faceOffsets[procMeshI]; @@ -540,8 +542,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh } // Add merged interfaces - const lduInterfacePtrsList interfaces = - procMeshes[procMeshI].interfaces(); + const lduInterfacePtrsList interfaces = procMesh.interfaces(); forAll(interfaces, intI) { @@ -596,7 +597,13 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh const lduInterfacePtrsList nbrInterfaces = - procMeshes[nbrProcMeshI].interfaces(); + mesh + ( + myMesh, + otherMeshes, + nbrProcMeshI + ).interfaces(); + const labelUList& faceCells = interfaces[intI].faceCells(); @@ -648,7 +655,13 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh label procMeshI = elems[i][0]; label interfaceI = elems[i][1]; const lduInterfacePtrsList interfaces = - procMeshes[procMeshI].interfaces(); + mesh + ( + myMesh, + otherMeshes, + procMeshI + ).interfaces(); + n += interfaces[interfaceI].faceCells().size(); } @@ -661,7 +674,12 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh label procMeshI = elems[i][0]; label interfaceI = elems[i][1]; const lduInterfacePtrsList interfaces = - procMeshes[procMeshI].interfaces(); + mesh + ( + myMesh, + otherMeshes, + procMeshI + ).interfaces(); boundaryMap[procMeshI][interfaceI] = allInterfaceI; labelList& bfMap = boundaryFaceMap[procMeshI][interfaceI]; @@ -755,87 +773,39 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +const Foam::lduMesh& Foam::lduPrimitiveMesh::mesh +( + const lduMesh& myMesh, + const PtrList<lduMesh>& otherMeshes, + const label meshI +) +{ + return (meshI == 0 ? myMesh : otherMeshes[meshI-1]); +} + + void Foam::lduPrimitiveMesh::gather ( const lduMesh& mesh, const labelList& procIDs, - PtrList<lduMesh>& procMeshes + PtrList<lduMesh>& otherMeshes ) { - const label meshComm = mesh.comm(); + // Force calculation of schedule (since does parallel comms) + (void)mesh.lduAddr().patchSchedule(); - lduInterfacePtrsList interfaces(mesh.interfaces()); - const lduAddressing& addressing = mesh.lduAddr(); + const label meshComm = mesh.comm(); if (Pstream::myProcNo(meshComm) == procIDs[0]) { - procMeshes.setSize(procIDs.size()); - - // Master mesh (copy for now. TBD) - { - // Convert to processorGAMGInterfaces in new communicator - lduInterfacePtrsList newInterfaces(interfaces.size()); - labelListList newPatchAddr(interfaces.size()); - - forAll(interfaces, intI) - { - if (interfaces.set(intI)) - { - const processorLduInterface& pldui = - refCast<const processorLduInterface>(interfaces[intI]); - - newPatchAddr[intI] = interfaces[intI].faceCells(); - labelList faceRestrictAddressing - ( - identity(interfaces[intI].faceCells().size()) - ); - - newInterfaces.set - ( - intI, - new processorGAMGInterface - ( - intI, - newInterfaces, - newPatchAddr[intI], - faceRestrictAddressing, - pldui.comm(), - pldui.myProcNo(), - pldui.neighbProcNo(), - pldui.forwardT(), - pldui.tag() - ) - ); - } - } - - - procMeshes.set - ( - 0, - new lduPrimitiveMesh - ( - addressing.size(), - addressing.lowerAddr(), - addressing.upperAddr(), - newPatchAddr, - newInterfaces, - nonBlockingSchedule<processorGAMGInterface> - ( - newInterfaces - ), - meshComm - ) - ); - } + otherMeshes.setSize(procIDs.size()-1); // Slave meshes for (label i = 1; i < procIDs.size(); i++) { //Pout<< "on master :" - // << " receiving from slave " << procIDs[i] - // << endl; + // << " receiving from slave " << procIDs[i] << endl; IPstream fromSlave ( @@ -849,46 +819,35 @@ void Foam::lduPrimitiveMesh::gather label nCells = readLabel(fromSlave); labelList lowerAddr(fromSlave); labelList upperAddr(fromSlave); - labelListList patchAddr(fromSlave); - labelList comm(fromSlave); - labelList myProcNo(fromSlave); - labelList neighbProcNo(fromSlave); - List<tensorField> forwardT(fromSlave); - labelList tag(fromSlave); - - // Convert to processorGAMGInterfaces - lduInterfacePtrsList newInterfaces(patchAddr.size()); - forAll(patchAddr, intI) + boolList validInterface(fromSlave); + + // Construct GAMGInterfaces + lduInterfacePtrsList newInterfaces(validInterface.size()); + labelListList patchAddr(validInterface.size()); + forAll(validInterface, intI) { - if (tag[intI] != -1) + if (validInterface[intI]) { - labelList faceRestrictAddressing - ( - identity(patchAddr[intI].size()) - ); + word coupleType(fromSlave); newInterfaces.set ( intI, - new processorGAMGInterface + GAMGInterface::New ( + coupleType, intI, newInterfaces, - patchAddr[intI], - faceRestrictAddressing, - comm[intI], - myProcNo[intI], - neighbProcNo[intI], - forwardT[intI], - tag[intI] - ) + fromSlave + ).ptr() ); + patchAddr[intI] = newInterfaces[intI].faceCells(); } } - procMeshes.set + otherMeshes.set ( - i, + i-1, new lduPrimitiveMesh ( nCells, @@ -900,7 +859,8 @@ void Foam::lduPrimitiveMesh::gather ( newInterfaces ), - meshComm + meshComm, + true ) ); } @@ -909,28 +869,12 @@ void Foam::lduPrimitiveMesh::gather { // Send to master - //- Extract info from processorGAMGInterface - labelListList patchAddr(interfaces.size()); - labelList comm(interfaces.size(), -1); - labelList myProcNo(interfaces.size(), -1); - labelList neighbProcNo(interfaces.size(), -1); - List<tensorField> forwardT(interfaces.size()); - labelList tag(interfaces.size(), -1); - + const lduAddressing& addressing = mesh.lduAddr(); + lduInterfacePtrsList interfaces(mesh.interfaces()); + boolList validInterface(interfaces.size()); forAll(interfaces, intI) { - if (interfaces.set(intI)) - { - const processorLduInterface& pldui = - refCast<const processorLduInterface>(interfaces[intI]); - - patchAddr[intI] = interfaces[intI].faceCells(); - comm[intI] = pldui.comm(); - myProcNo[intI] = pldui.myProcNo(); - neighbProcNo[intI] = pldui.neighbProcNo(); - forwardT[intI] = pldui.forwardT(); - tag[intI] = pldui.tag(); - } + validInterface[intI] = interfaces.set(intI); } OPstream toMaster @@ -941,38 +885,34 @@ void Foam::lduPrimitiveMesh::gather Pstream::msgType(), meshComm ); + + //Pout<< "sent nCells:" << addressing.size() + // << "sent lowerAddr:" << addressing.lowerAddr().size() + // << "sent upperAddr:" << addressing.upperAddr() + // << "sent validInterface:" << validInterface + // << endl; + toMaster << addressing.size() << addressing.lowerAddr() << addressing.upperAddr() - << patchAddr - << comm - << myProcNo - << neighbProcNo - << forwardT - << tag; - } -} + << validInterface; + forAll(interfaces, intI) + { + if (interfaces.set(intI)) + { + const GAMGInterface& interface = refCast<const GAMGInterface> + ( + interfaces[intI] + ); -// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // - -//Foam::Ostream& Foam::operator<< -//( -// Ostream& os, -// const Foam::lduPrimitiveMesh& mesh -//) -//{ -// os << mesh.size() -// << mesh.lowerAddr_ -// << mesh.upperAddr_ -// << mesh.patchAddr_ -// << mesh.interfaces_ -// << mesh.patchSchedule_ -// << mesh.comm_ -// -// return os; -//} + toMaster << interface.type(); + interface.write(toMaster); + } + } + } +} // ************************************************************************* // diff --git a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H index 01bbd0fe926..9826490b248 100644 --- a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H +++ b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H @@ -43,12 +43,6 @@ SourceFiles namespace Foam { -//// Forward declaration of friend functions and operators -// -//class lduPrimitiveMesh; -// -//Ostream& operator<<(Ostream&, const lduPrimitiveMesh&); - /*---------------------------------------------------------------------------*\ Class lduPrimitiveMesh Declaration \*---------------------------------------------------------------------------*/ @@ -132,9 +126,6 @@ public: bool reUse ); - ////- Construct from Istream - //lduPrimitiveMesh(Istream& is); - //- 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. @@ -149,13 +140,14 @@ public: // 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. + //// TBD: hardcoded for processorGAMGInterface. Should use some + //// 'clone' or New functionality instead. lduPrimitiveMesh ( const label comm, const labelList& procIDs, - const PtrList<lduMesh>& procMeshes, + const lduMesh& myMesh, + const PtrList<lduMesh>& otherMeshes, labelList& cellOffsets, labelListList& faceMap, @@ -219,19 +211,22 @@ public: // Helper - //- Gather meshes from processors. Hardcoded for - // processorGAMGInterface + //- Select either mesh0 (meshI is 0) or otherMeshes[meshI-1] + static const lduMesh& mesh + ( + const lduMesh& mesh0, + const PtrList<lduMesh>& otherMeshes, + const label meshI + ); + + //- Gather meshes from other processors onto procIDs[0]. + // Hardcoded for GAMGInterface. static void gather ( const lduMesh& mesh, const labelList& procIDs, - PtrList<lduMesh>& procMeshes + PtrList<lduMesh>& otherMeshes ); - -// -// // Ostream operator -// -// friend Ostream& operator<<(Ostream&, const procLduMatrix&); }; -- GitLab