ENH: GAMGSolver: initial version to use any solver on coarsest level

This commit is contained in:
mattijs 2013-02-22 17:25:33 +00:00
parent 2a896ded9e
commit 5dc5fa9d1b
3 changed files with 721 additions and 2 deletions

View File

@ -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<lduMesh>& procMeshes,
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
{
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<int> validRanks(fromSlave);
// Size coefficients
otherBouCoeffs.set
(
i-1,
new FieldField<Field, scalar>(validTransforms.size())
);
otherIntCoeffs.set
(
i-1,
new FieldField<Field, scalar>(validTransforms.size())
);
otherInterfaces.set
(
i-1,
new lduInterfaceFieldPtrsList(validTransforms.size())
);
forAll(validTransforms, intI)
{
if (validTransforms[intI])
{
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
(
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<int> validRanks(interfaceBouCoeffs.size(), -1);
forAll(interfaces, intI)
{
if (interfaces.set(intI))
{
const processorLduInterfaceField& interface =
refCast<const processorLduInterfaceField>
(
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<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
);
forAll(procMeshes, procI)
{
Pout<< "** procMesh " << procI << " "
<< procMeshes[procI].info()
<< endl;
}
Pout<< endl;
// Gather all matrix coefficients onto procIDs[0]
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
PtrList<lduMatrix> otherMats;
PtrList<FieldField<Field, scalar> > otherBouCoeffs;
PtrList<FieldField<Field, scalar> > otherIntCoeffs;
PtrList<lduInterfaceFieldPtrsList> 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<scalar>(allDiag, matrix.diag().size()).assign
(
matrix.diag()
);
forAll(otherMats, i)
{
SubList<scalar>
(
allDiag,
otherMats[i].diag().size(),
cellOffsets_[i+1]
).assign
(
otherMats[i].diag()
);
}
}
if (matrix.hasLower())
{
scalarField& allLower = allMatrix.lower();
UIndirectList<scalar>(allLower, faceMap_[0]) =
matrix.lower();
forAll(otherMats, i)
{
UIndirectList<scalar>(allLower, faceMap_[i+1]) =
otherMats[i].lower();
}
}
if (matrix.hasUpper())
{
scalarField& allUpper = allMatrix.upper();
UIndirectList<scalar>(allUpper, faceMap_[0]) =
matrix.upper();
forAll(otherMats, i)
{
UIndirectList<scalar>(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<Field, scalar>& procBouCoeffs
(
(procI == 0)
? interfaceBouCoeffs
: otherBouCoeffs[procI-1]
);
const FieldField<Field, scalar>& 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<const processorGAMGInterfaceField>
(
procInterfaces[procIntI]
);
allInterfaces_.set
(
allIntI,
GAMGInterfaceField::New
(
refCast<const GAMGInterface>
(
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;
}

View File

@ -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<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
autoPtr<lduMatrix> allMatrixPtr_;
FieldField<Field, scalar> allInterfaceBouCoeffs_;
FieldField<Field, scalar> 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<lduMesh>& procMeshes,
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;
//- 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
(

View File

@ -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<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,
@ -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<int>& 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<scalar>
(
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<scalar>
(
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;