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

Block AMG matrix selection with 2 norms

Browse source
This commit is contained in:
Hrvoje Jasak 2017-05-10 13:15:45 +01:00
parent f600c863c4
commit f0014817fe
3 changed files with 102 additions and 52 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

102
src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C
View file

@ -25,7 +25,7 @@ Class
BlockMatrixSelection BlockMatrixSelection
Description Description
Selective AMG policy for block matrices Classical AMG coarsening algorithm for block matrices.
Author Author
Tessa Uroic, FMENA Tessa Uroic, FMENA
@ -37,7 +37,6 @@ Author
#include "BlockSAMGInterfaceField.H" #include "BlockSAMGInterfaceField.H"
#include "coarseBlockAMGLevel.H" #include "coarseBlockAMGLevel.H"
#include "PriorityList.H" #include "PriorityList.H"
#include "crMatrix.H"
#include "SortableList.H" #include "SortableList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -47,7 +46,7 @@ template<class Type>
const Foam::debug::tolerancesSwitch const Foam::debug::tolerancesSwitch
Foam::BlockMatrixSelection<Type>::epsilon_ Foam::BlockMatrixSelection<Type>::epsilon_
( (
"samgWeightFactor", "samgStrongConnectionFactor",
0.25 0.25
); );
@ -71,9 +70,13 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
// Note: not taking norm magnitudes. HJ, 28/Feb/2017 // Note: not taking norm magnitudes. HJ, 28/Feb/2017
//------------------------------------------------------------------------------
// CALCULATE COARSENING NORM
//------------------------------------------------------------------------------
// Calculate norm for diagonal coefficients // Calculate norm for diagonal coefficients
scalarField normDiag(matrixSize); scalarField normDiag(matrixSize);
normPtr_->normalize(normDiag, matrix_.diag()); coarseningNormPtr_->normalize(normDiag, matrix_.diag());
// Note: this needs to be untangled for symmetric matrices. // Note: this needs to be untangled for symmetric matrices.
// If the matrix is symmetric and you ask for lower, it will be manufactured // If the matrix is symmetric and you ask for lower, it will be manufactured
@ -84,16 +87,32 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
// Calculate norm for upper triangle coeffs (magUpper) // Calculate norm for upper triangle coeffs (magUpper)
scalarField normUpper(row.size()); scalarField normUpper(row.size());
normPtr_->normalize(normUpper, matrix_.upper()); coarseningNormPtr_->normalize(normUpper, matrix_.upper());
// Calculate norm for lower triangle coeffs (magLower) // Calculate norm for lower triangle coeffs (magLower)
scalarField normLower(column.size()); scalarField normLower(column.size());
normPtr_->normalize(normLower, matrix_.lower()); coarseningNormPtr_->normalize(normLower, matrix_.lower());
// Calculate norm magnitudes // Calculate norm magnitudes
scalarField magNormUpper = mag(normUpper); scalarField magNormUpper = mag(normUpper);
scalarField magNormLower = mag(normLower); scalarField magNormLower = mag(normLower);
//------------------------------------------------------------------------------
// CALCULATE INTERPOLATION NORM
//------------------------------------------------------------------------------
// Calculate norm for diagonal coefficients (store into magDiag)
scalarField interpolationNormDiag(matrixSize);
coarseningNormPtr_->normalize(interpolationNormDiag, matrix_.diag());
// Calculate norm for upper triangle coeffs (magUpper)
scalarField interpolationNormUpper(column.size());
coarseningNormPtr_->normalize(interpolationNormUpper, matrix_.upper());
// Calculate norm for lower triangle coeffs (magLower)
scalarField interpolationNormLower(column.size());
coarseningNormPtr_->normalize(interpolationNormLower, matrix_.lower());
//------------------------------------------------------------------------------ //------------------------------------------------------------------------------
// CRITERIA FOR COARSENING: FIND STRONG CONNECTIONS FOR EACH ROW // CRITERIA FOR COARSENING: FIND STRONG CONNECTIONS FOR EACH ROW
//------------------------------------------------------------------------------ //------------------------------------------------------------------------------
@ -175,8 +194,8 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
// Set addressing for the matrix: // Set addressing for the matrix:
// stongCol and strongElement are arrays needed to create a compressed row // stongCol and strongElement are arrays needed to create a compressed row
// matrix // matrix
const labelList& strongRowStart = strong.crAddr().row(); const labelList& strongRowStart = strong.crAddr().rowStart();
labelList& strongCol = strong.col(); labelList& strongCol = strong.column();
scalarField& strongCoeff = strong.coeffs(); scalarField& strongCoeff = strong.coeffs();
// Counter for counting the strong elements in a row // Counter for counting the strong elements in a row
@ -198,7 +217,9 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
// to count the number of negative strong elements for // to count the number of negative strong elements for
// each row i // each row i
strongCol[strongRowStart[i] + counter[i]] = j; strongCol[strongRowStart[i] + counter[i]] = j;
strongCoeff[strongRowStart[i] + counter[i]] = normUpper[k]; strongCoeff[strongRowStart[i] + counter[i]] =
interpolationNormUpper[k];
counter[i]++; counter[i]++;
influence[j]++; influence[j]++;
} }
@ -207,7 +228,9 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
if (magNormLower[k] > epsLargestNorm[j]) if (magNormLower[k] > epsLargestNorm[j])
{ {
strongCol[strongRowStart[j] + counter[j]] = i; strongCol[strongRowStart[j] + counter[j]] = i;
strongCoeff[strongRowStart[j] + counter[j]] = normLower[k]; strongCoeff[strongRowStart[j] + counter[j]] =
interpolationNormLower[k];
counter[j]++; counter[j]++;
influence[i]++; influence[i]++;
} }
@ -220,8 +243,8 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
// Transpose the compressed row matrix to use for coarsening // Transpose the compressed row matrix to use for coarsening
crAddressing Taddr = strong.crAddr().T(); crAddressing Taddr = strong.crAddr().T();
const labelList& tRow = Taddr.row(); const labelList& tRow = Taddr.rowStart();
const labelList& tCol = Taddr.col(); const labelList& tCol = Taddr.column();
// Label the equations COARSE and FINE based on the number of // Label the equations COARSE and FINE based on the number of
// influences. // influences.
@ -380,20 +403,20 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
// Extract diagonal component to compare the sign with the // Extract diagonal component to compare the sign with the
// off-diagonal coeffs - search for negative and positive // off-diagonal coeffs - search for negative and positive
// connections // connections
scalar diagComponent = normDiag[i]; scalar signDiag = sign(interpolationNormDiag[i]);
// Upper triangle - use rowStart // Upper triangle - use rowStart
for (label k = rowStart[i]; k < rowStart[i + 1]; k++) for (label k = rowStart[i]; k < rowStart[i + 1]; k++)
{ {
// POSITIVE contribution // POSITIVE contribution
if (sign(normUpper[k]) == sign(diagComponent)) if (sign(normUpper[k]) == signDiag)
{ {
positiveElemSum[i] += normUpper[k]; positiveElemSum[i] += interpolationNormUpper[k];
} }
// NEGATIVE contribution // NEGATIVE contribution
else else
{ {
negativeElemSum[i] += normUpper[k]; negativeElemSum[i] += interpolationNormUpper[k];
} }
} }
@ -406,14 +429,14 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
label index = losortAddr[k]; label index = losortAddr[k];
// POSITIVE contribution // POSITIVE contribution
if (sign(normLower[index]) == sign(diagComponent)) if (sign(normLower[index]) == signDiag)
{ {
positiveElemSum[i] += normLower[index]; positiveElemSum[i] += interpolationNormLower[index];
} }
// NEGATIVE contribution // NEGATIVE contribution
else else
{ {
negativeElemSum[i] += normLower[index]; negativeElemSum[i] += interpolationNormLower[index];
} }
} }
} }
@ -426,7 +449,8 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
{ {
if (rowLabel_[i] == FINE) if (rowLabel_[i] == FINE)
{ {
scalar diagSign = sign(normDiag[i]); scalar signDiag = sign(interpolationNormDiag[i]);
for for
( (
label index = strongRowStart[i]; label index = strongRowStart[i];
@ -440,8 +464,8 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
// Sum positive connections // Sum positive connections
if if
( (
sign(strongCoeff[index]) == diagSign && sign(strongCoeff[index]) == signDiag
rowLabel_[j] == COARSE && rowLabel_[j] == COARSE
) )
{ {
posCoarseElemSum[i] += strongCoeff[index]; posCoarseElemSum[i] += strongCoeff[index];
@ -478,7 +502,7 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
if (rowLabel_[i] == FINE) if (rowLabel_[i] == FINE)
{ {
scalar diagonalCoeff = normDiag[i]; scalar diagonalCoeff = interpolationNormDiag[i];
// Check whether this Eqn has COARSE negative strong // Check whether this Eqn has COARSE negative strong
// contributions! // contributions!
@ -509,10 +533,8 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
{ {
prolongationCoeff.append prolongationCoeff.append
( (
(-positiveElemSum[i]) (-positiveElemSum[i])/posCoarseElemSum[i]*
/posCoarseElemSum[i] strongCoeff[k]/(diagonalCoeff)
*strongCoeff[k]
/(diagonalCoeff)
); );
prolongationCol.append(j - fineRowsCount[j]); prolongationCol.append(j - fineRowsCount[j]);
@ -650,7 +672,14 @@ Foam::BlockMatrixSelection<Type>::BlockMatrixSelection
: :
BlockMatrixCoarsening<Type>(matrix, dict, groupSize, minCoarseEqns), BlockMatrixCoarsening<Type>(matrix, dict, groupSize, minCoarseEqns),
matrix_(matrix), matrix_(matrix),
normPtr_(BlockCoeffNorm<Type>::New(dict)), coarseningNormPtr_
(
BlockCoeffNorm<Type>::New(dict.subDict("coarseningNorm"))
),
interpolationNormPtr_
(
BlockCoeffNorm<Type>::New(dict.subDict("interpolationNorm"))
),
nCoarseEqns_(0), nCoarseEqns_(0),
coarsen_(false), coarsen_(false),
Pptr_(NULL), Pptr_(NULL),
@ -723,8 +752,8 @@ Foam::BlockMatrixSelection<Type>::restrictMatrix() const
//------------------------------------------------------------------------------ //------------------------------------------------------------------------------
// Restriction addressing // Restriction addressing
const labelList& rowStartR = crR.row(); const labelList& rowStartR = crR.rowStart();
const labelList& colR = crR.col(); const labelList& colR = crR.column();
// Matrix A addressing // Matrix A addressing
const unallocLabelList& rowStartA = matrix_.lduAddr().ownerStartAddr(); const unallocLabelList& rowStartA = matrix_.lduAddr().ownerStartAddr();
@ -734,8 +763,8 @@ Foam::BlockMatrixSelection<Type>::restrictMatrix() const
const unallocLabelList& losortStart = matrix_.lduAddr().losortStartAddr(); const unallocLabelList& losortStart = matrix_.lduAddr().losortStartAddr();
// Prolongation addressing // Prolongation addressing
const labelList& rowStartP = crP.row(); const labelList& rowStartP = crP.rowStart();
const labelList& colP = crP.col(); const labelList& colP = crP.column();
// coeffLabel is used for checking if a contribution in that address already // coeffLabel is used for checking if a contribution in that address already
// exists // exists
@ -1058,6 +1087,7 @@ Foam::BlockMatrixSelection<Type>::restrictMatrix() const
SAMGInterface::New SAMGInterface::New
( (
coarseAddrPtr(), coarseAddrPtr(),
P,
coarseInterfaces, coarseInterfaces,
fineInterface, fineInterface,
fineInterface.interfaceInternalField(rowLabel_), fineInterface.interfaceInternalField(rowLabel_),
@ -1443,8 +1473,8 @@ void Foam::BlockMatrixSelection<Type>::restrictResidual
// Get restriction addressing // Get restriction addressing
const crAddressing& crR = R.crAddr(); const crAddressing& crR = R.crAddr();
const labelList& rowStartR = crR.row(); const labelList& rowStartR = crR.rowStart();
const labelList& colR = crR.col(); const labelList& colR = crR.column();
// Coefficients of restriction // Coefficients of restriction
const scalarField& coeffR = R.coeffs(); const scalarField& coeffR = R.coeffs();
@ -1483,8 +1513,8 @@ void Foam::BlockMatrixSelection<Type>::prolongateCorrection
// Get prolongation addressing // Get prolongation addressing
const crAddressing& crP = P.crAddr(); const crAddressing& crP = P.crAddr();
const labelList& rowStartP = crP.row(); const labelList& rowStartP = crP.rowStart();
const labelList& colP = crP.col(); const labelList& colP = crP.column();
const label sizeP = crP.nRows(); const label sizeP = crP.nRows();
// Coefficients of prolongation // Coefficients of prolongation

28
src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.H
View file

@ -25,10 +25,24 @@ Class
BlockMatrixSelection BlockMatrixSelection
Description Description
Selective AMG policy for block matrices Classical AMG coarsening algorithm for block matrices. The choice of coarse
equations is based on the largest number of strong connections with other
equations. For optimal performance, strong connections should have sign
opposite to diagonal coefficient. Coarse level matrix is obtained by
multiplying the restriction matrix, fine level matrix and prolongation
matrix. The algorithm closely follows theoretical background from work of
Klaus Stueben.
Author Author
Tessa Uroic, FMENA Tessa Uroic, FMENA, 2017
References
[1] T. Clees:
"AMG strategies for PDE systems with applications in industrial
semiconductor simulation", PhD Thesis, University of Cologne,
Germany, 2004
[2] U. Trottenberg, C. Oosterlee, A. Schueller:
"Multigrid", Elsevier, Academic Press, 2001
SourceFiles SourceFiles
BlockMatrixSelection.C BlockMatrixSelection.C
@ -65,8 +79,14 @@ class BlockMatrixSelection
const BlockLduMatrix<Type>& matrix_; const BlockLduMatrix<Type>& matrix_;
//- Reference to a templated norm calculator //- Reference to a templated norm calculator
//- INFO: Component norm is recommended, based on the pressure equation //- INFO: This norm will be calculated for coarsening process (it can be
autoPtr<BlockCoeffNorm<Type> > normPtr_; //- arbitrary: component norm, max norm, two norm)
autoPtr<BlockCoeffNorm<Type> > coarseningNormPtr_;
//- Reference to a templated norm calculator
//- INFO: This norm will be used for interpolation. Component norm based
//- on a representative equation is recommended.
autoPtr<BlockCoeffNorm<Type> > interpolationNormPtr_;
//- Number of coarse equations //- Number of coarse equations
label nCoarseEqns_; label nCoarseEqns_;

24
src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockSelectiveAMG/BlockSelectiveAMG.C
View file

@ -175,8 +175,8 @@ void Foam::BlockSelectiveAMG<Type>::calcCoarsening()
// Set addressing for the matrix: // Set addressing for the matrix:
// stongCol and strongElement are arrays needed to create a compressed row // stongCol and strongElement are arrays needed to create a compressed row
// matrix // matrix
const labelList& strongRowStart = strong.crAddr().row(); const labelList& strongRowStart = strong.crAddr().rowStart();
labelList& strongCol = strong.col(); labelList& strongCol = strong.column();
scalarField& strongCoeff = strong.coeffs(); scalarField& strongCoeff = strong.coeffs();
// Counter for counting the strong elements in a row // Counter for counting the strong elements in a row
@ -220,8 +220,8 @@ void Foam::BlockSelectiveAMG<Type>::calcCoarsening()
// Transpose the compressed row matrix to use for coarsening // Transpose the compressed row matrix to use for coarsening
crAddressing Taddr = strong.crAddr().T(); crAddressing Taddr = strong.crAddr().T();
const labelList& tRow = Taddr.row(); const labelList& tRow = Taddr.rowStart();
const labelList& tCol = Taddr.col(); const labelList& tCol = Taddr.column();
// Label the equations COARSE and FINE based on the number of // Label the equations COARSE and FINE based on the number of
// influences. // influences.
@ -705,8 +705,8 @@ Foam::BlockSelectiveAMG<Type>::restrictMatrix() const
// GET ADDRESSING // GET ADDRESSING
//------------------------------------------------------------------------------ //------------------------------------------------------------------------------
// Restriction addressing // Restriction addressing
const labelList& rowStartR = crR.row(); const labelList& rowStartR = crR.rowStart();
const labelList& colR = crR.col(); const labelList& colR = crR.column();
// Matrix A addressing // Matrix A addressing
const unallocLabelList& rowStartA = matrix_.lduAddr().ownerStartAddr(); const unallocLabelList& rowStartA = matrix_.lduAddr().ownerStartAddr();
@ -716,8 +716,8 @@ Foam::BlockSelectiveAMG<Type>::restrictMatrix() const
const unallocLabelList& losortStart = matrix_.lduAddr().losortStartAddr(); const unallocLabelList& losortStart = matrix_.lduAddr().losortStartAddr();
// Prolongation addressing // Prolongation addressing
const labelList& rowStartP = crP.row(); const labelList& rowStartP = crP.rowStart();
const labelList& colP = crP.col(); const labelList& colP = crP.column();
//------------------------------------------------------------------------------ //------------------------------------------------------------------------------
// COUNTING COARSE COEFFICIENTS WHICH WILL APPEAR IN TRIPLE PRODUCT // COUNTING COARSE COEFFICIENTS WHICH WILL APPEAR IN TRIPLE PRODUCT
@ -1278,8 +1278,8 @@ void Foam::BlockSelectiveAMG<Type>::restrictResidual
// Get restriction addressing // Get restriction addressing
const crAddressing& crR = R.crAddr(); const crAddressing& crR = R.crAddr();
const labelList& rowStartR = crR.row(); const labelList& rowStartR = crR.rowStart();
const labelList& colR = crR.col(); const labelList& colR = crR.column();
// Coefficients of restriction // Coefficients of restriction
const scalarField& coeffR = R.coeffs(); const scalarField& coeffR = R.coeffs();
@ -1318,8 +1318,8 @@ void Foam::BlockSelectiveAMG<Type>::prolongateCorrection
// Get prolongation addressing // Get prolongation addressing
const crAddressing& crP = P.crAddr(); const crAddressing& crP = P.crAddr();
const labelList& rowStartP = crP.row(); const labelList& rowStartP = crP.rowStart();
const labelList& colP = crP.col(); const labelList& colP = crP.column();
const label sizeP = crP.nRows(); const label sizeP = crP.nRows();
// Coefficients of prolongation // Coefficients of prolongation