exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Block AMG solver updates

Browse source
This commit is contained in:
Hrvoje Jasak 2014-05-29 00:18:30 +01:00
parent 75c037bf4a
commit 05dc011f83
8 changed files with 290 additions and 160 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
src/foam/matrices/blockLduMatrix/BlockAmg/BlockAmgCycle.C
View file

@ -133,10 +133,7 @@ void Foam::BlockAmgCycle<Type>::fixedCycle
Field<Type>& bCoarse = coarseLevelPtr_->levelPtr_->b(); Field<Type>& bCoarse = coarseLevelPtr_->levelPtr_->b();
// Zero out coarse x // Zero out coarse x
forAll (xCoarse,i) xCoarse = pTraits<Type>::zero;
{
xCoarse[i] *= 0.0;
}
// Restrict residual: optimisation on number of pre-sweeps // Restrict residual: optimisation on number of pre-sweeps
levelPtr_->restrictResidual levelPtr_->restrictResidual

327
src/foam/matrices/blockLduMatrix/BlockAmg/BlockMatrixCoarsening/BlockMatrixAgglomeration.C
View file

@ -25,7 +25,7 @@ Class
BlockMatrixAgglomeration BlockMatrixAgglomeration
Description Description
Agglomerative block matrix AMG corsening, adjusted for BlockLduMatrix Agglomerative block matrix AMG corsening
Author Author
Klas Jareteg, 2012-12-13 Klas Jareteg, 2012-12-13
@ -33,6 +33,7 @@ Author
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "BlockMatrixAgglomeration.H" #include "BlockMatrixAgglomeration.H"
#include "boolList.H"
#include "coeffFields.H" #include "coeffFields.H"
#include "addToRunTimeSelectionTable.H" #include "addToRunTimeSelectionTable.H"
#include "BlockGAMGInterfaceField.H" #include "BlockGAMGInterfaceField.H"
@ -40,6 +41,12 @@ Author
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
template<class Type>
const Foam::scalar Foam::BlockMatrixAgglomeration<Type>::diagFactor_
(
debug::tolerances("aamgDiagFactor", 1e-8)
);
template<class Type> template<class Type>
const Foam::scalar Foam::BlockMatrixAgglomeration<Type>::weightFactor_ = 0.65; const Foam::scalar Foam::BlockMatrixAgglomeration<Type>::weightFactor_ = 0.65;
@ -47,7 +54,7 @@ const Foam::scalar Foam::BlockMatrixAgglomeration<Type>::weightFactor_ = 0.65;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Type> template<class Type>
void Foam::BlockMatrixAgglomeration<Type>::calcChild() void Foam::BlockMatrixAgglomeration<Type>::calcAgglomeration()
{ {
// Algorithm: // Algorithm:
// 1) Create temporary equation addressing using a double-pass algorithm. // 1) Create temporary equation addressing using a double-pass algorithm.
@ -96,9 +103,9 @@ void Foam::BlockMatrixAgglomeration<Type>::calcChild()
// Create column and coefficient index array // Create column and coefficient index array
{ {
// Use child array to count number of entries per row. // Use agglomIndex to count number of entries per row.
// Reset the list for counting // Reset the list for counting
labelList& nPerRow = child_; labelList& nPerRow = agglomIndex_;
nPerRow = 0; nPerRow = 0;
forAll (upperAddr, coeffI) forAll (upperAddr, coeffI)
@ -123,8 +130,8 @@ void Foam::BlockMatrixAgglomeration<Type>::calcChild()
nPerRow[lowerAddr[coeffI]]++; nPerRow[lowerAddr[coeffI]]++;
} }
// Reset child array // Reset agglomeration index array
child_ = -1; agglomIndex_ = -1;
} }
@ -132,43 +139,122 @@ void Foam::BlockMatrixAgglomeration<Type>::calcChild()
// Get matrix coefficients // Get matrix coefficients
const CoeffField<Type>& diag = matrix_.diag(); const CoeffField<Type>& diag = matrix_.diag();
const CoeffField<Type>& upper = matrix_.upper();
// Coefficient magnitudes are pre-calculated
scalarField magDiag(diag.size());
scalarField magOffDiag(upperAddr.size());
normPtr_->coeffMag(diag, magDiag);
if (matrix_.asymmetric())
{
scalarField magUpper(upperAddr.size());
scalarField magLower(upperAddr.size());
normPtr_->coeffMag(matrix_.upper(), magUpper);
normPtr_->coeffMag(matrix_.lower(), magLower);
magOffDiag = Foam::max(magUpper, magLower);
}
else if (matrix_.symmetric())
{
normPtr_->coeffMag(matrix_.upper(), magOffDiag);
}
else
{
// Diag only matrix. Reset and return
agglomIndex_ = 0;
nCoarseEqns_ = 1;
return;
}
labelList sizeOfGroups(nRows, 0); labelList sizeOfGroups(nRows, 0);
nCoarseEqns_ = 0; nCoarseEqns_ = 0;
label indexU, indexG, colI, curEqn, nextEqn, groupPassI; // Gather disconnected and weakly connected equations into cluster zero
// Weak connection is assumed to be the one where the off-diagonal
// coefficient is smaller than diagFactor_*diag
{
// Algorithm
// Mark all cells to belong to zero cluster
// Go through all upper and lower coefficients and for the ones
// larger than threshold mark the equations out of cluster zero
scalar magRowDiag, magColDiag, weight, weightU, weightG; scalarField magScaledDiag = diagFactor_*magDiag;
// Magnitudes pre-calculated boolList zeroCluster(diag.size(), true);
Field<scalar> magDiag(diag.size());
Field<scalar> magUpper(upper.size());
normPtr_->coeffMag(diag,magDiag); forAll (magOffDiag, coeffI)
normPtr_->coeffMag(upper,magUpper); {
if (magOffDiag[coeffI] > magScaledDiag[upperAddr[coeffI]])
{
zeroCluster[upperAddr[coeffI]] = false;
}
if (magOffDiag[coeffI] > magScaledDiag[lowerAddr[coeffI]])
{
zeroCluster[lowerAddr[coeffI]] = false;
}
}
// Collect solo equations
label nSolo = 0;
forAll (zeroCluster, eqnI)
{
if (zeroCluster[eqnI])
{
// Found solo equation
nSolo++;
agglomIndex_[eqnI] = nCoarseEqns_;
}
}
reduce(nSolo, sumOp<label>());
if (nSolo > 0)
{
// Found solo equations
nCoarseEqns_++;
if (BlockLduMatrix<Type>::debug >= 2)
{
Info<< "Found " << nSolo << " weakly connected equations."
<< endl;
}
}
}
// Go through the off-diagonal and create clusters, marking the child array
label indexUngrouped, indexGrouped;
label colI, curEqn, nextEqn, groupPassI;
scalar magRowDiag, magColDiag;
scalar weight, weightUngrouped, weightGrouped;
for (label eqnI = 0; eqnI < nRows; eqnI++) for (label eqnI = 0; eqnI < nRows; eqnI++)
{ {
if (child_[eqnI] == -1) if (agglomIndex_[eqnI] == -1)
{ {
curEqn = eqnI; curEqn = eqnI;
indexU = -1; indexUngrouped = -1;
indexG = -1; indexGrouped = -1;
child_[curEqn] = nCoarseEqns_; agglomIndex_[curEqn] = nCoarseEqns_;
magRowDiag = magDiag[curEqn]; magRowDiag = magDiag[curEqn];
for (groupPassI = 1; groupPassI < groupSize_; groupPassI++) for (groupPassI = 1; groupPassI < groupSize_; groupPassI++)
{ {
weightU = 0; weightUngrouped = 0;
weightG = 0; weightGrouped = 0;
indexU = -1; indexUngrouped = -1;
indexG = -1; indexGrouped = -1;
for for
( (
@ -181,40 +267,42 @@ void Foam::BlockMatrixAgglomeration<Type>::calcChild()
magColDiag = magDiag[colI]; magColDiag = magDiag[colI];
weight = (magUpper[cIndex[rowCoeffI]] weight = magOffDiag[cIndex[rowCoeffI]]/max(magRowDiag, magColDiag);
/ max(magRowDiag, magColDiag));
if (child_[colI] == -1) if (agglomIndex_[colI] == -1)
{ {
if (indexU == -1 || weight > weightU) if (indexUngrouped == -1 || weight > weightUngrouped)
{ {
indexU = rowCoeffI; indexUngrouped = rowCoeffI;
weightU = weight; weightUngrouped = weight;
} }
} }
else if (child_[curEqn] != child_[colI]) else if (agglomIndex_[curEqn] != agglomIndex_[colI])
{ {
if (indexG == -1 || weight > weightG) if (indexGrouped == -1 || weight > weightGrouped)
{ {
indexG = rowCoeffI; indexGrouped = rowCoeffI;
weightG = weight; weightGrouped = weight;
} }
} }
} }
if if
( (
indexU != -1 indexUngrouped != -1
&& (indexG == -1 || weightU >= weightFactor_*weightG) && (
indexGrouped == -1
|| weightUngrouped >= weightFactor_*weightGrouped
)
) )
{ {
// Found new element of group. Add it and use as // Found new element of group. Add it and use as
// start of next search // start of next search
nextEqn = cols[indexU]; nextEqn = cols[indexUngrouped];
child_[nextEqn] = child_[curEqn]; agglomIndex_[nextEqn] = agglomIndex_[curEqn];
sizeOfGroups[child_[curEqn]]; sizeOfGroups[agglomIndex_[curEqn]];
curEqn = nextEqn; curEqn = nextEqn;
} }
@ -228,17 +316,20 @@ void Foam::BlockMatrixAgglomeration<Type>::calcChild()
if if
( (
groupPassI > 1 groupPassI > 1
|| indexG == -1 || indexGrouped == -1
|| sizeOfGroups[child_[cols[indexG]]] > (groupSize_ + 2) || (
sizeOfGroups[agglomIndex_[cols[indexGrouped]]]
> (groupSize_ + 2)
)
) )
{ {
sizeOfGroups[child_[eqnI]]++; sizeOfGroups[agglomIndex_[eqnI]]++;
nCoarseEqns_++; nCoarseEqns_++;
} }
else else
{ {
child_[eqnI] = child_[cols[indexG]]; agglomIndex_[eqnI] = agglomIndex_[cols[indexGrouped]];
sizeOfGroups[child_[cols[indexG]]]++; sizeOfGroups[agglomIndex_[cols[indexGrouped]]]++;
} }
} }
} }
@ -258,6 +349,19 @@ void Foam::BlockMatrixAgglomeration<Type>::calcChild()
reduce(coarsen_, andOp<bool>()); reduce(coarsen_, andOp<bool>());
if (BlockLduMatrix<Type>::debug >= 2)
{
Pout << "Coarse level size: " << nCoarseEqns_;
if (coarsen_)
{
Pout << ". Accepted" << endl;
}
else
{
Pout << ". Rejected" << endl;
}
}
} }
@ -274,19 +378,13 @@ Foam::BlockMatrixAgglomeration<Type>::BlockMatrixAgglomeration
: :
BlockMatrixCoarsening<Type>(matrix, dict, groupSize, minCoarseEqns), BlockMatrixCoarsening<Type>(matrix, dict, groupSize, minCoarseEqns),
matrix_(matrix), matrix_(matrix),
normPtr_ normPtr_(BlockCoeffNorm<Type>::New(dict)),
( agglomIndex_(matrix_.lduAddr().size()),
BlockCoeffNorm<Type>::New
(
dict
)
),
child_(matrix_.lduAddr().size()),
groupSize_(groupSize), groupSize_(groupSize),
nCoarseEqns_(0), nCoarseEqns_(0),
coarsen_(false) coarsen_(false)
{ {
calcChild(); calcAgglomeration();
} }
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
@ -314,7 +412,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
// Construct the coarse matrix and ldu addressing for the next level // Construct the coarse matrix and ldu addressing for the next level
// Algorithm: // Algorithm:
// 1) Loop through all fine coeffs. If the child labels on two sides are // 1) Loop through all fine coeffs. If the agglom labels on two sides are
// different, this creates a coarse coeff. Define owner and neighbour // different, this creates a coarse coeff. Define owner and neighbour
// for this coeff based on cluster IDs. // for this coeff based on cluster IDs.
// 2) Check if the coeff has been seen before. If yes, add the coefficient // 2) Check if the coeff has been seen before. If yes, add the coefficient
@ -332,14 +430,14 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
const label nFineCoeffs = upperAddr.size(); const label nFineCoeffs = upperAddr.size();
# ifdef FULLDEBUG # ifdef FULLDEBUG
if (child_.size() != matrix_.lduAddr().size()) if (agglomIndex_.size() != matrix_.lduAddr().size())
{ {
FatalErrorIn FatalErrorIn
( (
"autoPtr<BlockLduMatrix<Type> > BlockMatrixAgglomeration<Type>::" "autoPtr<BlockLduMatrix<Type> > BlockMatrixAgglomeration<Type>::"
"restrictMatrix() const" "restrictMatrix() const"
) << "Child array does not correspond to fine level. " << endl ) << "agglomIndex array does not correspond to fine level. " << endl
<< " Child size: " << child_.size() << " Size: " << agglomIndex_.size()
<< " number of equations: " << matrix_.lduAddr().size() << " number of equations: " << matrix_.lduAddr().size()
<< abort(FatalError); << abort(FatalError);
} }
@ -369,8 +467,8 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
// Loop through all fine coeffs // Loop through all fine coeffs
forAll (upperAddr, fineCoeffi) forAll (upperAddr, fineCoeffi)
{ {
label rmUpperAddr = child_[upperAddr[fineCoeffi]]; label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffi]];
label rmLowerAddr = child_[lowerAddr[fineCoeffi]]; label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffi]];
if (rmUpperAddr == rmLowerAddr) if (rmUpperAddr == rmLowerAddr)
{ {
@ -480,6 +578,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
initCoarseNeighb.setSize(0); initCoarseNeighb.setSize(0);
coarseCoeffMap.setSize(0); coarseCoeffMap.setSize(0);
// Create coarse-level coupled interfaces // Create coarse-level coupled interfaces
// Create coarse interfaces, addressing and coefficients // Create coarse interfaces, addressing and coefficients
@ -517,7 +616,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
interfaceFields[intI].interface().initInternalFieldTransfer interfaceFields[intI].interface().initInternalFieldTransfer
( (
Pstream::blocking, Pstream::blocking,
child_ agglomIndex_
); );
} }
} }
@ -541,7 +640,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
fineInterface.internalFieldTransfer fineInterface.internalFieldTransfer
( (
Pstream::blocking, Pstream::blocking,
child_ agglomIndex_
) )
) )
); );
@ -563,7 +662,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
( (
*coarseAddrPtr, *coarseAddrPtr,
fineInterface, fineInterface,
fineInterface.interfaceInternalField(child_), fineInterface.interfaceInternalField(agglomIndex_),
fineInterfaceAddr[intI] fineInterfaceAddr[intI]
).ptr() ).ptr()
); );
@ -664,7 +763,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
squareTypeField& activeCoarseLower = coarseLower.asSquare(); squareTypeField& activeCoarseLower = coarseLower.asSquare();
const squareTypeField& activeFineLower = fineLower.asSquare(); const squareTypeField& activeFineLower = fineLower.asSquare();
restrictResidual(fineDiag, coarseDiag); restrictDiag(fineDiag, coarseDiag);
forAll(coeffRestrictAddr, fineCoeffI) forAll(coeffRestrictAddr, fineCoeffI)
{ {
@ -688,27 +787,41 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
} }
else if (fineUpper.activeType() == blockCoeffBase::LINEAR) else if (fineUpper.activeType() == blockCoeffBase::LINEAR)
{ {
FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const") FatalErrorIn
<< "Matrix diagonal of square type and upper of linear type is not implemented" (
"autoPtr<amgMatrix>"
"BlockMatrixAgglomeration<Type>::restrictMatrix() const"
) << "Matrix diagonal of square type and upper of "
<< "linear type is not implemented"
<< abort(FatalError); << abort(FatalError);
} }
else else
{ {
FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const") FatalErrorIn
<< "Matrix diagonal of square type and upper of scalar type is not implemented" (
"autoPtr<amgMatrix>"
"BlockMatrixAgglomeration<Type>::restrictMatrix() const"
) << "Matrix diagonal of square type and upper of "
<< "scalar type is not implemented"
<< abort(FatalError); << abort(FatalError);
} }
} }
else if (fineDiag.activeType() == blockCoeffBase::LINEAR) else if (fineDiag.activeType() == blockCoeffBase::LINEAR)
{ {
FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const") FatalErrorIn
<< "Matrix diagonal of linear type not implemented" (
"autoPtr<amgMatrix>"
"BlockMatrixAgglomeration<Type>::restrictMatrix() const"
) << "Matrix diagonal of linear type not implemented"
<< abort(FatalError); << abort(FatalError);
} }
else else
{ {
FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const") FatalErrorIn
<< "Matrix diagonal of scalar type not implemented" (
"autoPtr<amgMatrix>"
"BlockMatrixAgglomeration<Type>::restrictMatrix() const"
) << "Matrix diagonal of scalar type not implemented"
<< abort(FatalError); << abort(FatalError);
} }
} }
@ -722,7 +835,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
squareTypeField& activeCoarseDiag = coarseDiag.asSquare(); squareTypeField& activeCoarseDiag = coarseDiag.asSquare();
const squareTypeField& activeFineUpper = fineUpper.asSquare(); const squareTypeField& activeFineUpper = fineUpper.asSquare();
restrictResidual(fineDiag, coarseDiag); restrictDiag(fineDiag, coarseDiag);
forAll(coeffRestrictAddr, fineCoeffI) forAll(coeffRestrictAddr, fineCoeffI)
{ {
@ -730,27 +843,38 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
if (cCoeff >= 0) if (cCoeff >= 0)
{ {
activeCoarseUpper[cCoeff] activeCoarseUpper[cCoeff] +=
+= activeFineUpper[fineCoeffI]; activeFineUpper[fineCoeffI];
} }
else else
{ {
// Add the fine face coefficient into the diagonal. // Add the fine face coefficient into the diagonal
activeCoarseDiag[-1 - cCoeff] // Note: upper and lower coeffs are transpose of
+= 2*activeFineUpper[fineCoeffI]; // each other. HJ, 28/May/2014
activeCoarseDiag[-1 - cCoeff] +=
activeFineUpper[fineCoeffI]
+ activeFineUpper[fineCoeffI].T();
} }
} }
} }
else if (fineUpper.activeType() == blockCoeffBase::LINEAR) else if (fineUpper.activeType() == blockCoeffBase::LINEAR)
{ {
FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const") FatalErrorIn
<< "Matrix diagonal of square type and upper of linear type is not implemented" (
"autoPtr<amgMatrix>"
"BlockMatrixAgglomeration<Type>::restrictMatrix() const"
) << "Matrix diagonal of square type and upper of "
<< "linear type is not implemented"
<< abort(FatalError); << abort(FatalError);
} }
else else
{ {
FatalErrorIn("autoPtr<amgMatrix> BlockMatrixAgglomeration<Type>::restrictMatrix() const") FatalErrorIn
<< "Matrix diagonal of square type and upper of scalar type is not implemented" (
"autoPtr<amgMatrix>"
"BlockMatrixAgglomeration<Type>::restrictMatrix() const"
) << "Matrix diagonal of square type and upper of "
<< "scalar type is not implemented"
<< abort(FatalError); << abort(FatalError);
} }
@ -775,48 +899,47 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
} }
} }
return autoPtr<BlockLduMatrix<Type> > return autoPtr<BlockLduMatrix<Type> >(coarseMatrixPtr);
(
new BlockLduMatrix<Type>
(
coarseMatrix
)
);
} }
template<class Type> template<class Type>
void Foam::BlockMatrixAgglomeration<Type>::restrictResidual void Foam::BlockMatrixAgglomeration<Type>::restrictDiag
( (
const CoeffField<Type>& res, const CoeffField<Type>& Coeff,
CoeffField<Type>& coarseRes CoeffField<Type>& coarseCoeff
) const ) const
{ {
typedef CoeffField<Type> TypeCoeffField; typedef CoeffField<Type> TypeCoeffField;
if (res.activeType() == blockCoeffBase::SQUARE && if
coarseRes.activeType() == blockCoeffBase::SQUARE) (
Coeff.activeType() == blockCoeffBase::SQUARE
&& coarseCoeff.activeType() == blockCoeffBase::SQUARE
)
{ {
typedef typename TypeCoeffField::squareType squareType;
typedef typename TypeCoeffField::squareTypeField squareTypeField; typedef typename TypeCoeffField::squareTypeField squareTypeField;
squareTypeField& activeCoarseRes = coarseRes.asSquare(); squareTypeField& activeCoarseCoeff = coarseCoeff.asSquare();
const squareTypeField& activeRes = res.asSquare(); const squareTypeField& activeCoeff = Coeff.asSquare();
// KRJ: 2013-02-05: To be done by correct pTraits? forAll (coarseCoeff, i)
forAll (coarseRes, i)
{ {
activeCoarseRes[i] *= 0.0; activeCoarseCoeff[i] = pTraits<squareType>::zero;
} }
forAll (res, i) forAll (Coeff, i)
{ {
activeCoarseRes[child_[i]] += activeRes[i]; activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i];
} }
} }
else else
{ {
FatalErrorIn("void BlockMatrixAgglomeration<Type>::restrictResidual() const") FatalErrorIn
<< "Only present for square type coeff type" (
"void BlockMatrixAgglomeration<Type>::restrictDiag() const"
) << "Only present for square type coeff type"
<< abort(FatalError); << abort(FatalError);
} }
} }
@ -833,7 +956,7 @@ void Foam::BlockMatrixAgglomeration<Type>::restrictResidual
forAll (res, i) forAll (res, i)
{ {
coarseRes[child_[i]] += res[i]; coarseRes[agglomIndex_[i]] += res[i];
} }
} }
@ -847,7 +970,7 @@ void Foam::BlockMatrixAgglomeration<Type>::prolongateCorrection
{ {
forAll (x, i) forAll (x, i)
{ {
x[i] += coarseX[child_[i]]; x[i] += coarseX[agglomIndex_[i]];
} }
} }

29
src/foam/matrices/blockLduMatrix/BlockAmg/BlockMatrixCoarsening/BlockMatrixAgglomeration.H
View file

@ -59,9 +59,6 @@ class BlockMatrixAgglomeration
: :
public BlockMatrixCoarsening<Type> public BlockMatrixCoarsening<Type>
{ {
public:
// Private Data // Private Data
//- Reference to matrix //- Reference to matrix
@ -70,8 +67,9 @@ public:
//- Reference to a templated norm calculator //- Reference to a templated norm calculator
autoPtr<BlockCoeffNorm<Type> > normPtr_; autoPtr<BlockCoeffNorm<Type> > normPtr_;
//- Child array //- Child array: for each fine equation give agglomeration cluster
labelField child_; // index
labelField agglomIndex_;
//- Group size //- Group size
label groupSize_; label groupSize_;
@ -82,6 +80,7 @@ public:
//- Can a coarse level be constructed? //- Can a coarse level be constructed?
bool coarsen_; bool coarsen_;
// Private Member Functions // Private Member Functions
//- Disallow default bitwise copy construct //- Disallow default bitwise copy construct
@ -90,12 +89,15 @@ public:
// Disallow default bitwise assignment // Disallow default bitwise assignment
void operator=(const BlockMatrixAgglomeration<Type>&); void operator=(const BlockMatrixAgglomeration<Type>&);
//- Calculate child //- Calculate agglomeration index (child)
void calcChild(); void calcAgglomeration();
// Private Static Data // Private Static Data
//- Diagonal scaling factor
static const scalar diagFactor_;
//- Weighting factor //- Weighting factor
static const scalar weightFactor_; static const scalar weightFactor_;
@ -117,6 +119,7 @@ public:
const label minCoarseEqns const label minCoarseEqns
); );
// Destructor // Destructor
virtual ~BlockMatrixAgglomeration(); virtual ~BlockMatrixAgglomeration();
@ -131,8 +134,7 @@ public:
} }
//- Restrict matrix //- Restrict matrix
virtual autoPtr<BlockLduMatrix<Type> > restrictMatrix virtual autoPtr<BlockLduMatrix<Type> > restrictMatrix() const;
() const;
//- Restrict residual //- Restrict residual
virtual void restrictResidual virtual void restrictResidual
@ -141,11 +143,11 @@ public:
Field<Type>& coarseRes Field<Type>& coarseRes
) const; ) const;
//- Restrict residual using CoeffField //- Restrict CoeffField. Used for diag coefficient
void restrictResidual void restrictDiag
( (
const CoeffField<Type>& res, const CoeffField<Type>& Coeff,
CoeffField<Type>& coarseRes CoeffField<Type>& coarseCoeff
) const; ) const;
//- Prolongate correction //- Prolongate correction
@ -154,7 +156,6 @@ public:
Field<Type>& x, Field<Type>& x,
const Field<Type>& coarseX const Field<Type>& coarseX
) const; ) const;
}; };

29
src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C
View file

@ -198,6 +198,7 @@ void Foam::coarseBlockAmgLevel<Type>::solve
label maxIter = Foam::min(2*coarseningPtr_->minCoarseEqns(), 1000); label maxIter = Foam::min(2*coarseningPtr_->minCoarseEqns(), 1000);
// Create artificial dictionary for top-level solution
dictionary topLevelDict; dictionary topLevelDict;
topLevelDict.add("nDirections", "5"); topLevelDict.add("nDirections", "5");
topLevelDict.add("minIter", 1); topLevelDict.add("minIter", 1);
@ -233,15 +234,10 @@ void Foam::coarseBlockAmgLevel<Type>::solve
matrixPtr_, matrixPtr_,
topLevelDict topLevelDict
).solve(x, b); ).solve(x, b);
if (BlockLduMatrix<Type>::debug >= 2)
{
coarseSolverPerf.print();
}
} }
else else
{ {
topLevelDict.add("preconditioner", "ILU0"); topLevelDict.add("preconditioner", "Cholesky");
coarseSolverPerf = coarseSolverPerf =
//BlockGaussSeidelSolver<Type> //BlockGaussSeidelSolver<Type>
@ -252,6 +248,25 @@ void Foam::coarseBlockAmgLevel<Type>::solve
matrixPtr_, matrixPtr_,
topLevelDict topLevelDict
).solve(x, b); ).solve(x, b);
}
// Escape cases of top-level solver divergence
if
(
coarseSolverPerf.nIterations() == maxIter
&& (
coarseSolverPerf.finalResidual()
>= coarseSolverPerf.initialResidual()
)
)
{
// Top-level solution failed. Attempt rescue
// HJ, 27/Jul/2013
multiply(x, invDiag, b);
// Print top level correction failure as info for user
coarseSolverPerf.print();
}
if (BlockLduMatrix<Type>::debug >= 2) if (BlockLduMatrix<Type>::debug >= 2)
{ {
@ -259,8 +274,6 @@ void Foam::coarseBlockAmgLevel<Type>::solve
} }
} }
}
template<class Type> template<class Type>
void Foam::coarseBlockAmgLevel<Type>::scaleX void Foam::coarseBlockAmgLevel<Type>::scaleX

17
src/foam/matrices/blockLduMatrix/BlockAmg/fineBlockAmgLevel.C
View file

@ -180,17 +180,26 @@ void Foam::fineBlockAmgLevel<Type>::solve
const scalar relTol const scalar relTol
) const ) const
{ {
Info<< "Fine level solver"<<endl; Info<< "Fine level solver"<<endl;
// Create artificial dictionary for finest solution
dictionary finestDict;
// finestDict.add("nDirections", "5");
finestDict.add("minIter", 1);
finestDict.add("maxIter", 1000);
finestDict.add("tolerance", tolerance);
finestDict.add("relTol", relTol);
if (matrix_.symmetric()) if (matrix_.symmetric())
{ {
finestDict.add("preconditioner", "Cholesky");
BlockSolverPerformance<Type> coarseSolverPerf = BlockSolverPerformance<Type> coarseSolverPerf =
BlockCGSolver<Type> BlockCGSolver<Type>
( (
"topLevelCorr", "topLevelCorr",
matrix_, matrix_,
dict_ finestDict
).solve(x, b); ).solve(x, b);
if (BlockLduMatrix<Type>::debug >= 2) if (BlockLduMatrix<Type>::debug >= 2)
@ -200,12 +209,14 @@ void Foam::fineBlockAmgLevel<Type>::solve
} }
else else
{ {
finestDict.add("preconditioner", "Cholesky");
BlockSolverPerformance<Type> coarseSolverPerf = BlockSolverPerformance<Type> coarseSolverPerf =
BlockBiCGStabSolver<Type> BlockBiCGStabSolver<Type>
( (
"topLevelCorr", "topLevelCorr",
matrix_, matrix_,
dict_ finestDict
).solve(x, b); ).solve(x, b);
if (BlockLduMatrix<Type>::debug >= 2) if (BlockLduMatrix<Type>::debug >= 2)

26
src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAmgSolver/BlockAmgSolver.C
View file

@ -31,12 +31,8 @@ Author
#include "BlockAmgSolver.H" #include "BlockAmgSolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
//- Construct from matrix and solver data stream
template<class Type> template<class Type>
Foam::BlockAmgSolver<Type>::BlockAmgSolver Foam::BlockAmgSolver<Type>::BlockAmgSolver
( (
@ -76,34 +72,22 @@ Foam::BlockAmgSolver<Type>::solve
this->fieldName() this->fieldName()
); );
// Create local references to avoid the spread this-> ugliness
const BlockLduMatrix<Type>& matrix = this->matrix_;
scalar norm = this->normFactor(x, b); scalar norm = this->normFactor(x, b);
Field<Type> wA(x.size()); // Calculate initial residual
solverPerf.initialResidual() = gSum(cmptMag(amg_.residual(x, b)))/norm;
// Calculate residual. Note: sign of residual swapped for efficiency
matrix.Amul(wA, x);
wA -= b;
solverPerf.initialResidual() = gSum(cmptMag(wA))/norm;
solverPerf.finalResidual() = solverPerf.initialResidual(); solverPerf.finalResidual() = solverPerf.initialResidual();
if (!solverPerf.checkConvergence(this->tolerance(), this->relTolerance())) if (!this->stop(solverPerf))
{ {
do do
{ {
amg_.cycle(x, b); amg_.cycle(x, b);
// Re-calculate residual. Note: sign of residual swapped solverPerf.finalResidual() =
// for efficiency gSum(cmptMag(amg_.residual(x, b)))/norm;
matrix.Amul(wA, x);
wA -= b;
solverPerf.finalResidual() = gSum(cmptMag(wA))/norm;
solverPerf.nIterations()++; solverPerf.nIterations()++;
} while (!this->stop(solverPerf)); } while (!this->stop(solverPerf));
} }

1
src/lduSolvers/amg/amgPolicy/pamgPolicy.C
View file

@ -29,6 +29,7 @@ Description
Author Author
Hrvoje Jasak, Wikki Ltd. All rights reserved Hrvoje Jasak, Wikki Ltd. All rights reserved
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "pamgPolicy.H" #include "pamgPolicy.H"

2
src/lduSolvers/amg/coarseAmgLevel.C
View file

@ -237,7 +237,7 @@ void Foam::coarseAmgLevel::solve
).solve(x, b, cmpt); ).solve(x, b, cmpt);
} }
// Escape cases of solver divergence // Escape cases of top-level solver divergence
if if
( (
coarseSolverPerf.nIterations() == maxIter coarseSolverPerf.nIterations() == maxIter