From 5dc5fa9d1be4b6a476752b3e31fbe164e977ca4e Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 22 Feb 2013 17:25:33 +0000 Subject: [PATCH] ENH: GAMGSolver: initial version to use any solver on coarsest level --- .../lduMatrix/solvers/GAMG/GAMGSolver.C | 505 +++++++++++++++++- .../lduMatrix/solvers/GAMG/GAMGSolver.H | 51 ++ .../lduMatrix/solvers/GAMG/GAMGSolverSolve.C | 167 +++++- 3 files changed, 721 insertions(+), 2 deletions(-) diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C index f0275284f6..8877ed01cd 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C @@ -24,6 +24,10 @@ License \*---------------------------------------------------------------------------*/ #include "GAMGSolver.H" +#include "processorGAMGInterface.H" +#include "processorLduInterfaceField.H" +#include "GAMGInterfaceField.H" +#include "processorGAMGInterfaceField.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -39,6 +43,145 @@ namespace Foam } +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +void Foam::GAMGSolver::gatherMatrices +( + const labelList& procIDs, + const PtrList& procMeshes, + + const lduMatrix& mat, + const FieldField& interfaceBouCoeffs, + const FieldField& interfaceIntCoeffs, + const lduInterfaceFieldPtrsList& interfaces, + + PtrList& otherMats, + PtrList >& otherBouCoeffs, + PtrList >& otherIntCoeffs, + PtrList& otherInterfaces +) const +{ + const label meshComm = mat.mesh().comm(); + + //lduInterfacePtrsList interfaces(mesh.interfaces()); + + if (Pstream::myProcNo(meshComm) == procIDs[0]) + { + // Master. + otherMats.setSize(procIDs.size()-1); + otherBouCoeffs.setSize(procIDs.size()-1); + otherIntCoeffs.setSize(procIDs.size()-1); + otherInterfaces.setSize(procIDs.size()-1); + + for (label i = 1; i < procIDs.size(); i++) + { + const lduMesh& procMesh = procMeshes[i]; + const lduInterfacePtrsList procInterfaces = procMesh.interfaces(); + + IPstream fromSlave + ( + Pstream::scheduled, + procIDs[i], + 0, // bufSize + Pstream::msgType(), + meshComm + ); + + otherMats.set(i-1, new lduMatrix(procMesh, fromSlave)); + + // Receive number of/valid interfaces + boolList validTransforms(fromSlave); + List validRanks(fromSlave); + + // Size coefficients + otherBouCoeffs.set + ( + i-1, + new FieldField(validTransforms.size()) + ); + otherIntCoeffs.set + ( + i-1, + new FieldField(validTransforms.size()) + ); + otherInterfaces.set + ( + i-1, + new lduInterfaceFieldPtrsList(validTransforms.size()) + ); + + forAll(validTransforms, intI) + { + if (validTransforms[intI]) + { + const processorGAMGInterface& interface = + refCast + ( + procInterfaces[intI] + ); + + + otherBouCoeffs[i-1].set(intI, new scalarField(fromSlave)); + otherIntCoeffs[i-1].set(intI, new scalarField(fromSlave)); + otherInterfaces[i-1].set + ( + intI, + GAMGInterfaceField::New + ( + interface, //procInterfaces[intI], + validTransforms[intI], + validRanks[intI] + ).ptr() + ); + } + } + } + } + else + { + // Send to master + OPstream toMaster + ( + Pstream::scheduled, + procIDs[0], + 0, + Pstream::msgType(), + meshComm + ); + + + // Count valid interfaces + boolList validTransforms(interfaceBouCoeffs.size(), false); + List validRanks(interfaceBouCoeffs.size(), -1); + forAll(interfaces, intI) + { + if (interfaces.set(intI)) + { + const processorLduInterfaceField& interface = + refCast + ( + interfaces[intI] + ); + + validTransforms[intI] = interface.doTransform(); + validRanks[intI] = interface.rank(); + } + } + + toMaster << mat << validTransforms << validRanks; + forAll(validTransforms, intI) + { + if (validTransforms[intI]) + { + toMaster + << interfaceBouCoeffs[intI] + << interfaceIntCoeffs[intI]; + } + } + } +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::GAMGSolver::GAMGSolver @@ -74,6 +217,7 @@ Foam::GAMGSolver::GAMGSolver interpolateCorrection_(false), scaleCorrection_(matrix.symmetric()), directSolveCoarsest_(false), + processorAgglomerate_(false), agglomeration_(GAMGAgglomeration::New(matrix_, controlDict_)), matrixLevels_(agglomeration_.size()), @@ -94,10 +238,15 @@ Foam::GAMGSolver::GAMGSolver if (directSolveCoarsest_) { - label coarseComm = matrixLevels_[coarsestLevel].mesh().comm(); + const lduMesh& coarsestMesh = matrixLevels_[coarsestLevel].mesh(); + + label coarseComm = coarsestMesh.comm(); label oldWarn = UPstream::warnComm; UPstream::warnComm = coarseComm; + Pout<< "Solve direct on coasestmesh (level=" << coarsestLevel + << ") using communicator " << coarseComm << endl; + coarsestLUMatrixPtr_.set ( new LUscalarMatrix @@ -108,6 +257,344 @@ Foam::GAMGSolver::GAMGSolver ) ); + UPstream::warnComm = oldWarn; + } + else if (processorAgglomerate_) + { + const lduMesh& coarsestMesh = matrixLevels_[coarsestLevel].mesh(); + + label coarseComm = coarsestMesh.comm(); + label oldWarn = UPstream::warnComm; + UPstream::warnComm = coarseComm; + + Pout<< "Solve generic on coarsestmesh (level=" << coarsestLevel + << ") using communicator " << coarseComm << endl; + + // All processors to master + const List& procIDs = UPstream::procID(coarseComm); + + // Gather all meshes onto procIDs[0] + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // Note: move into GAMGAgglomeration + + Pout<< "procIDs:" << procIDs << endl; + + PtrList procMeshes; + lduPrimitiveMesh::gather + ( + coarsestMesh, + procIDs, + procMeshes + ); + + forAll(procMeshes, procI) + { + Pout<< "** procMesh " << procI << " " + << procMeshes[procI].info() + << endl; + } + Pout<< endl; + + + // Gather all matrix coefficients onto procIDs[0] + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + PtrList otherMats; + PtrList > otherBouCoeffs; + PtrList > otherIntCoeffs; + PtrList otherInterfaces; + gatherMatrices + ( + procIDs, + procMeshes, + + matrix, + interfaceBouCoeffs, + interfaceIntCoeffs, + interfaces, + + otherMats, + otherBouCoeffs, + otherIntCoeffs, + otherInterfaces + ); + + Pout<< "Own matrix:" << matrix.info() << endl; + + forAll(otherMats, i) + { + Pout<< "** otherMats " << i << " " + << otherMats[i].info() + << endl; + } + Pout<< endl; + + + if (Pstream::myProcNo(coarseComm) == procIDs[0]) + { + // Agglomerate all addressing + // ~~~~~~~~~~~~~~~~~~~~~~~~~~ + + allMeshPtr_.reset + ( + new lduPrimitiveMesh + ( + coarseComm, + procIDs, + procMeshes, + + cellOffsets_, + faceMap_, + boundaryMap_, + boundaryFaceMap_ + ) + ); + const lduMesh& allMesh = allMeshPtr_(); + + + // Agglomerate all matrix + // ~~~~~~~~~~~~~~~~~~~~~~ + + allMatrixPtr_.reset(new lduMatrix(allMesh)); + lduMatrix& allMatrix = allMatrixPtr_(); + + if (matrix.hasDiag()) + { + scalarField& allDiag = allMatrix.diag(); + SubList(allDiag, matrix.diag().size()).assign + ( + matrix.diag() + ); + forAll(otherMats, i) + { + SubList + ( + allDiag, + otherMats[i].diag().size(), + cellOffsets_[i+1] + ).assign + ( + otherMats[i].diag() + ); + } + } + if (matrix.hasLower()) + { + scalarField& allLower = allMatrix.lower(); + UIndirectList(allLower, faceMap_[0]) = + matrix.lower(); + forAll(otherMats, i) + { + UIndirectList(allLower, faceMap_[i+1]) = + otherMats[i].lower(); + } + } + if (matrix.hasUpper()) + { + scalarField& allUpper = allMatrix.upper(); + UIndirectList(allUpper, faceMap_[0]) = + matrix.upper(); + forAll(otherMats, i) + { + UIndirectList(allUpper, 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(boundaryMap_, procI) + { + const FieldField& procBouCoeffs + ( + (procI == 0) + ? interfaceBouCoeffs + : otherBouCoeffs[procI-1] + ); + const FieldField& procIntCoeffs + ( + (procI == 0) + ? interfaceIntCoeffs + : otherIntCoeffs[procI-1] + ); + const lduInterfaceFieldPtrsList procInterfaces + ( + (procI == 0) + ? interfaces + : otherInterfaces[procI-1] + ); + + + const labelList& bMap = boundaryMap_[procI]; + forAll(bMap, procIntI) + { + const scalarField& procBou = procBouCoeffs[procIntI]; + const scalarField& procInt = procIntCoeffs[procIntI]; + + label allIntI = bMap[procIntI]; + + if (allIntI != -1) + { + // So this boundary has been preserved. Copy + // data across. + + if (!allInterfaces_.set(allIntI)) + { + // Construct lduInterfaceField + + const processorGAMGInterfaceField& procInt = + refCast + ( + procInterfaces[procIntI] + ); + + allInterfaces_.set + ( + allIntI, + GAMGInterfaceField::New + ( + refCast + ( + allMeshInterfaces[allIntI] + ), + procInt.doTransform(), + procInt.rank() + ).ptr() + ); + } + + + // Map data from processor to complete mesh + + scalarField& allBou = + allInterfaceBouCoeffs_[allIntI]; + scalarField& allInt = + allInterfaceIntCoeffs_[allIntI]; + + const labelList& map = + boundaryFaceMap_[procI][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 + + // 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& allU = + allMesh.lduAddr().upperAddr(); + const labelList& allL = + allMesh.lduAddr().lowerAddr(); + const label off = cellOffsets_[procI]; + + forAll(map, i) + { + if (map[i] >= 0) + { + FatalErrorIn + ( + "GAMGSolver::GAMGSolver()" + ) << "problem." + << abort(FatalError); + } + + label allFaceI = -(map[i]+1); + 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); + } + + + if (allL[allFaceI] == allCellI) + { + if (matrix.hasUpper()) + { + allMatrix.upper()[allFaceI] = + procBou[i]; + } + if (matrix.hasLower()) + { + allMatrix.lower()[allFaceI] = + procInt[i]; + } + } + else + { + if (matrix.hasUpper()) + { + allMatrix.upper()[allFaceI] = + procInt[i]; + } + if (matrix.hasLower()) + { + allMatrix.lower()[allFaceI] = + procBou[i]; + } + } + } + } + } + } + } + UPstream::warnComm = oldWarn; } } @@ -184,6 +671,22 @@ void Foam::GAMGSolver::readControls() controlDict_.readIfPresent("interpolateCorrection", interpolateCorrection_); controlDict_.readIfPresent("scaleCorrection", scaleCorrection_); controlDict_.readIfPresent("directSolveCoarsest", directSolveCoarsest_); + controlDict_.readIfPresent("processorAgglomerate", processorAgglomerate_); + + Pout<< "GAMGSolver settings :" + << " cacheAgglomeration:" << cacheAgglomeration_ + << " nPreSweeps:" << nPreSweeps_ + << " preSweepsLevelMultiplier:" << preSweepsLevelMultiplier_ + << " maxPreSweeps:" << maxPreSweeps_ + << " nPostSweeps:" << nPostSweeps_ + << " postSweepsLevelMultiplier:" << postSweepsLevelMultiplier_ + << " maxPostSweeps:" << maxPostSweeps_ + << " nFinestSweeps:" << nFinestSweeps_ + << " interpolateCorrection:" << interpolateCorrection_ + << " scaleCorrection:" << scaleCorrection_ + << " directSolveCoarsest:" << directSolveCoarsest_ + << " processorAgglomerate:" << processorAgglomerate_ + << endl; } diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H index f62a9341a2..136e6d7923 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H @@ -109,6 +109,9 @@ class GAMGSolver //- Direct or iteratively solve the coarsest level bool directSolveCoarsest_; + //- Agglomerate processors + bool processorAgglomerate_; + //- The agglomeration const GAMGAgglomeration& agglomeration_; @@ -129,6 +132,28 @@ class GAMGSolver autoPtr coarsestLUMatrixPtr_; + + // Processor addressing agglomeration - to be moved into + // GAMGAgglomeration + + //- Combined coarsest meshes for all processors + autoPtr allMeshPtr_; + + //- Mapping from processor to coarsestAllMesh + labelList cellOffsets_; + labelListList faceMap_; + labelListList boundaryMap_; + labelListListList boundaryFaceMap_; + + // Processor matrix agglomeration + + //- Combined coarsest level matrix for all processors + autoPtr allMatrixPtr_; + FieldField allInterfaceBouCoeffs_; + FieldField allInterfaceIntCoeffs_; + lduInterfaceFieldPtrsList allInterfaces_; + + // Private Member Functions //- Read control parameters from the control dictionary @@ -158,6 +183,32 @@ class GAMGSolver //- Agglomerate coarse matrix void agglomerateMatrix(const label fineLevelIndex); + //- Collect matrices from other processors + void gatherMatrices + ( + const labelList& procIDs, + const PtrList& procMeshes, + + const lduMatrix& mat, + const FieldField& interfaceBouCoeffs, + const FieldField& interfaceIntCoeffs, + const lduInterfaceFieldPtrsList& interfaces, + + PtrList& otherMats, + PtrList >& otherBouCoeffs, + PtrList >& otherIntCoeffs, + PtrList& otherInterfaces + ) const; + + //- Gather field in allMatrix form + void gatherField + ( + const label meshComm, + const labelList& procIDs, + const scalarField& coarsestSource, + scalarField& allSource + ) const; + //- Interpolate the correction after injected prolongation void interpolate ( diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C index edd9fefaad..f9088e0b77 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C @@ -31,6 +31,67 @@ License // * * * * * * * * * * * * * * * 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 + ( + 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 slaveSlots + ( + allSource, + cellOffsets_[i+1]-cellOffsets_[i], + cellOffsets_[i] + ); + + UList& 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, @@ -451,9 +512,113 @@ void Foam::GAMGSolver::solveCoarsestLevel coarsestCorrField = coarsestSource; coarsestLUMatrixPtr_->solve(coarsestCorrField); } + else if (processorAgglomerate_) + { + Pout<< "** Processor agglomeration" << endl; + + const lduMatrix& allMatrix = allMatrixPtr_(); + //const lduMesh& allMesh = allMeshPtr_(); + + const List& procIDs = UPstream::procID(coarseComm); + + scalarField allSource; + gatherField + ( + coarseComm, + procIDs, + coarsestSource, + allSource + ); + + scalarField allCorrField; + solverPerformance coarseSolverPerf; + + if (Pstream::myProcNo(coarseComm) == procIDs[0]) + { + allCorrField.setSize(allSource.size(), 0); + + if (allMatrix.asymmetric()) + { + coarseSolverPerf = BICCG + ( + "coarsestLevelCorr", + allMatrix, + allInterfaceBouCoeffs_, + allInterfaceIntCoeffs_, + allInterfaces_, + tolerance_, + relTol_ + ).solve + ( + allCorrField, + allSource + ); + } + else + { + coarseSolverPerf = ICCG + ( + "coarsestLevelCorr", + allMatrix, + allInterfaceBouCoeffs_, + allInterfaceIntCoeffs_, + allInterfaces_, + tolerance_, + relTol_ + ).solve + ( + allCorrField, + allSource + ); + } + + + // Distribute correction + coarsestCorrField = SubList + ( + allCorrField, + coarsestSource.size(), + 0 + ); + for (label i=1; i < procIDs.size(); i++) + { + OPstream toSlave + ( + Pstream::scheduled, + procIDs[i], + 0, // bufSize + Pstream::msgType(), + coarseComm + ); + + toSlave << SubList + ( + allCorrField, + cellOffsets_[i+1]-cellOffsets_[i], + cellOffsets_[i] + ); + } + } + else + { + IPstream fromMaster + ( + Pstream::scheduled, + 0, + 0, // bufSize + Pstream::msgType(), + coarseComm + ); + fromMaster >> coarsestCorrField; + } + + if (debug >= 2) + { + coarseSolverPerf.print(Info(coarseComm)); + } + } else { - const label coarsestLevel = matrixLevels_.size() - 1; coarsestCorrField = 0; solverPerformance coarseSolverPerf;