diff --git a/src/foam/Make/files b/src/foam/Make/files index ef723aefb..70a89b2f0 100644 --- a/src/foam/Make/files +++ b/src/foam/Make/files @@ -718,9 +718,7 @@ $(BlockSAMGInterfaceFields)/GGIBlockSAMGInterfaceField/GGIBlockSAMGInterfaceFiel BlockMatrixCoarsening = $(BlockAMG)/BlockMatrixCoarsening $(BlockMatrixCoarsening)/BlockMatrixCoarsening/blockMatrixCoarsenings.C -$(BlockMatrixCoarsening)/BlockMatrixAgglomeration/blockMatrixAgglomerations.C $(BlockMatrixCoarsening)/BlockMatrixClustering/blockMatrixClusterings.C -$(BlockMatrixCoarsening)/BlockMatrixReorderedClustering/blockMatrixReorderedClusterings.C $(BlockMatrixCoarsening)/BlockMatrixSelection/blockMatrixSelections.C BlockLduPrecons = matrices/blockLduMatrix/BlockLduPrecons diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGCycle.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGCycle.C index 101d1b9e6..8a8145ac3 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGCycle.C +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGCycle.C @@ -107,6 +107,22 @@ void Foam::BlockAMGCycle::makeCoarseLevels(const label nMaxLevels) } +template +Foam::BlockLduMatrix& Foam::BlockAMGCycle::coarseMatrix() +{ + if (!coarseLevelPtr_) + { + FatalErrorIn + ( + "BlockLduMatrix& BlockAMGCycle::coarseMatrix()" + ) << "Coarse level not available" + << abort(FatalError); + } + + return coarseLevelPtr_->matrix(); +} + + template void Foam::BlockAMGCycle::fixedCycle ( @@ -193,7 +209,21 @@ void Foam::BlockAMGCycle::fixedCycle else { // Call direct solver - levelPtr_->solve(x, b, 1e-9, 0); + levelPtr_->solve(x, b, 1e-7, 0); + } +} + + +template +void Foam::BlockAMGCycle::initMatrix() +{ + // If present, update coarse levels recursively + if (coarseLevelPtr_) + { + // Update current level + levelPtr_->initLevel(coarseLevelPtr_->levelPtr_); + + coarseLevelPtr_->initMatrix(); } } diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGCycle.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGCycle.H index 56da581f5..517d7f90c 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGCycle.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGCycle.H @@ -143,6 +143,9 @@ public: return nLevels_; } + //- Return coarse matrix + BlockLduMatrix& coarseMatrix(); + //- Calculate residual void residual ( @@ -166,6 +169,9 @@ public: const label nPostSweeps, const bool scale ) const; + + //- Re-initialise preconditioner after matrix coefficient update + void initMatrix(); }; diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGLevel.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGLevel.H index 9db62dc0e..b1a119fcd 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGLevel.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGLevel.H @@ -86,6 +86,9 @@ public: // Member Functions + //- Return reference to matrix + virtual BlockLduMatrix& matrix() = 0; + //- Return reference to x virtual Field& x() = 0; @@ -144,6 +147,12 @@ public: //- Create next level from current level virtual autoPtr makeNextLevel() const = 0; + + //- Re-initialise AMG level after matrix coefficient update + virtual void initLevel + ( + autoPtr >& coarseLevelPtr + ) = 0; }; diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.C deleted file mode 100644 index 13c0d280e..000000000 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.C +++ /dev/null @@ -1,1312 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Class - BlockMatrixAgglomeration - -Description - Agglomerative block matrix AMG corsening - -Author - Klas Jareteg, 2012-12-13 - -\*---------------------------------------------------------------------------*/ - -#include "BlockMatrixAgglomeration.H" -#include "boolList.H" -#include "tolerancesSwitch.H" -#include "coeffFields.H" -#include "BlockAMGInterfaceField.H" -#include "coarseBlockAMGLevel.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -template -const Foam::debug::tolerancesSwitch -Foam::BlockMatrixAgglomeration::weightFactor_ -( - "aamgWeightFactor", - 0.65 -); - - -template -const Foam::debug::tolerancesSwitch -Foam::BlockMatrixAgglomeration::diagFactor_ -( - "aamgDiagFactor", - 1e-8 -); - - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -template -void Foam::BlockMatrixAgglomeration::calcAgglomeration() -{ - // Algorithm: - // 1) Create temporary equation addressing using a double-pass algorithm. - // to create the offset table. - // 2) Loop through all equations and for each equation find the best fit - // neighbour. If all neighbours are grouped, add equation to best group - - // Create row-based addressing - - const label nRows = matrix_.lduAddr().size(); - - const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); - const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); - - // For each equation get number of coefficients in a row - labelList cols(upperAddr.size() + lowerAddr.size()); - labelList cIndex(upperAddr.size() + lowerAddr.size()); - labelList rowOffset(nRows + 1, 0); - - // Count the number of coefficients - forAll (upperAddr, coeffI) - { - rowOffset[upperAddr[coeffI]]++; - } - - forAll (lowerAddr, coeffI) - { - rowOffset[lowerAddr[coeffI]]++; - } - - label nCoeffs = 0; - - forAll (rowOffset, eqnI) - { - nCoeffs += rowOffset[eqnI]; - } - - rowOffset[nRows] = nCoeffs; - - for (label eqnI = nRows - 1; eqnI >= 0; --eqnI) - { - rowOffset[eqnI] = rowOffset[eqnI + 1] - rowOffset[eqnI]; - } - - rowOffset[0] = 0; - - // Create column and coefficient index array - { - // Use agglomIndex to count number of entries per row. - // Reset the list for counting - labelList& nPerRow = agglomIndex_; - nPerRow = 0; - - forAll (upperAddr, coeffI) - { - cols[rowOffset[upperAddr[coeffI]] + nPerRow[upperAddr[coeffI]]] = - lowerAddr[coeffI]; - - cIndex[rowOffset[upperAddr[coeffI]] + nPerRow[upperAddr[coeffI]]] = - coeffI; - - nPerRow[upperAddr[coeffI]]++; - } - - forAll (lowerAddr, coeffI) - { - cols[rowOffset[lowerAddr[coeffI]] + nPerRow[lowerAddr[coeffI]]] = - upperAddr[coeffI]; - - cIndex[rowOffset[lowerAddr[coeffI]] + nPerRow[lowerAddr[coeffI]]] = - coeffI; - - nPerRow[lowerAddr[coeffI]]++; - } - - // Reset agglomeration index array - agglomIndex_ = -1; - } - - - // Calculate agglomeration - - // Get matrix coefficients - const CoeffField& diag = matrix_.diag(); - - // Coefficient magnitudes are pre-calculated - scalarField magDiag(diag.size()); - scalarField magOffDiag(upperAddr.size()); - - // Calculate norm of block coefficient. Note: the norm - // may be signed, ie. it will take the sign of the coefficient - normPtr_->normalize(magDiag, diag); - - // Take magnitude of the norm - mag(magDiag, magDiag); - - if (matrix_.asymmetric()) - { - scalarField magUpper(upperAddr.size()); - scalarField magLower(upperAddr.size()); - - // Take max of norm magnitude in upper and lower triangle - normPtr_->normalize(magUpper, matrix_.upper()); - mag(magUpper, magUpper); - - normPtr_->normalize(magLower, matrix_.lower()); - mag(magLower, magLower); - - Foam::max(magOffDiag, magUpper, magLower); - } - else if (matrix_.symmetric()) - { - // Take max of norm magnitude - normPtr_->normalize(magOffDiag, matrix_.upper()); - mag(magOffDiag, magOffDiag); - } - else - { - // Diag only matrix. Reset and return - agglomIndex_ = 0; - nCoarseEqns_ = 1; - - return; - } - - labelList sizeOfGroups(nRows, 0); - - nCoarseEqns_ = 0; - - // 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 - - scalarField magScaledDiag = diagFactor_()*magDiag; - - boolList zeroCluster(diag.size(), true); - - forAll (magOffDiag, coeffI) - { - if (magOffDiag[coeffI] > magScaledDiag[upperAddr[coeffI]]) - { - zeroCluster[upperAddr[coeffI]] = false; - } - - if (magOffDiag[coeffI] > magScaledDiag[lowerAddr[coeffI]]) - { - zeroCluster[lowerAddr[coeffI]] = false; - } - } - - // Collect solo equations - forAll (zeroCluster, eqnI) - { - if (zeroCluster[eqnI]) - { - // Found solo equation - nSolo_++; - - agglomIndex_[eqnI] = nCoarseEqns_; - } - } - - if (nSolo_ > 0) - { - // Found solo equations - sizeOfGroups[nCoarseEqns_] = nSolo_; - - nCoarseEqns_++; - } - } - - // 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++) - { - if (agglomIndex_[eqnI] == -1) - { - curEqn = eqnI; - - indexUngrouped = -1; - indexGrouped = -1; - - agglomIndex_[curEqn] = nCoarseEqns_; - - magRowDiag = magDiag[curEqn]; - - for (groupPassI = 1; groupPassI < groupSize_; groupPassI++) - { - weightUngrouped = 0; - weightGrouped = 0; - - indexUngrouped = -1; - indexGrouped = -1; - - for - ( - label rowCoeffI = rowOffset[curEqn]; - rowCoeffI < rowOffset[curEqn + 1]; - rowCoeffI++ - ) - { - colI = cols[rowCoeffI]; - - magColDiag = magDiag[colI]; - - weight = magOffDiag[cIndex[rowCoeffI]]/ - max(magRowDiag, magColDiag); - - if (agglomIndex_[colI] == -1) - { - if (indexUngrouped == -1 || weight > weightUngrouped) - { - indexUngrouped = rowCoeffI; - weightUngrouped = weight; - } - } - else if (agglomIndex_[curEqn] != agglomIndex_[colI]) - { - if (indexGrouped == -1 || weight > weightGrouped) - { - indexGrouped = rowCoeffI; - weightGrouped = weight; - } - } - } - - if - ( - indexUngrouped != -1 - && ( - indexGrouped == -1 - || weightUngrouped >= weightFactor_()*weightGrouped - ) - ) - { - // Found new element of group. Add it and use as - // start of next search - - nextEqn = cols[indexUngrouped]; - - agglomIndex_[nextEqn] = agglomIndex_[curEqn]; - sizeOfGroups[agglomIndex_[curEqn]]++; - - curEqn = nextEqn; - } - else - { - // Group full or cannot be extended - break; - } - } - - if - ( - groupPassI > 1 - || indexGrouped == -1 - || ( - sizeOfGroups[agglomIndex_[cols[indexGrouped]]] - > (groupSize_ + 2) - ) - ) - { - sizeOfGroups[agglomIndex_[eqnI]]++; - nCoarseEqns_++; - } - else - { - agglomIndex_[eqnI] = agglomIndex_[cols[indexGrouped]]; - sizeOfGroups[agglomIndex_[cols[indexGrouped]]]++; - } - } - } - - // Resize group size count - sizeOfGroups.setSize(nCoarseEqns_); - - // The decision on parallel agglomeration needs to be made for the - // whole gang of processes; otherwise I may end up with a different - // number of agglomeration levels on different processors. - - // If the number of coarse equations is less than minimum and - // if the matrix has reduced in size by at least 1/3, coarsen - if - ( - nCoarseEqns_ > BlockMatrixCoarsening::minCoarseEqns() - && 3*nCoarseEqns_ <= 2*nRows - ) - { - coarsen_ = true; - } - - reduce(coarsen_, andOp()); - - if (blockLduMatrix::debug >= 3) - { - // Count singleton clusters - label nSingleClusters = 0; - - // Adjust start based on solo cells cluster status - label start = 0; - - if (nSolo_ > 0) - { - start = 1; - } - - for (label i = start; i < sizeOfGroups.size(); i++) - { - if (sizeOfGroups[i] == 1) - { - nSingleClusters++; - } - } - - Pout<< "Coarse level size: " << nCoarseEqns_ - << " nSolo = " << nSolo_ - << " cluster (" << min(sizeOfGroups) << " " - << max(sizeOfGroups) - << "). N singleton clusters = " << nSingleClusters; - - if (coarsen_) - { - Pout << ". Accepted" << endl; - } - else - { - Pout << ". Rejected" << endl; - } - } -} - - -template -void Foam::BlockMatrixAgglomeration::restrictDiag -( - const CoeffField& Coeff, - CoeffField& coarseCoeff -) const -{ - typedef CoeffField TypeCoeffField; - - if - ( - Coeff.activeType() == blockCoeffBase::SQUARE - && coarseCoeff.activeType() == blockCoeffBase::SQUARE - ) - { - typedef typename TypeCoeffField::squareType squareType; - typedef typename TypeCoeffField::squareTypeField squareTypeField; - - squareTypeField& activeCoarseCoeff = coarseCoeff.asSquare(); - const squareTypeField& activeCoeff = Coeff.asSquare(); - - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } - - forAll (Coeff, i) - { - activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; - } - } - else if - ( - Coeff.activeType() == blockCoeffBase::LINEAR - && coarseCoeff.activeType() == blockCoeffBase::LINEAR - ) - { - typedef typename TypeCoeffField::linearType linearType; - typedef typename TypeCoeffField::linearTypeField linearTypeField; - - linearTypeField& activeCoarseCoeff = coarseCoeff.asLinear(); - const linearTypeField& activeCoeff = Coeff.asLinear(); - - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } - - forAll (Coeff, i) - { - activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; - } - } - else if - ( - Coeff.activeType() == blockCoeffBase::SCALAR - && coarseCoeff.activeType() == blockCoeffBase::SCALAR - ) - { - typedef typename TypeCoeffField::scalarType scalarType; - typedef typename TypeCoeffField::scalarTypeField scalarTypeField; - - scalarTypeField& activeCoarseCoeff = coarseCoeff.asScalar(); - const scalarTypeField& activeCoeff = Coeff.asScalar(); - - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } - - forAll (Coeff, i) - { - activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; - } - } - else - { - FatalErrorIn - ( - "void BlockMatrixAgglomeration::restrictDiag() const" - ) << "Problem in coeff type morphing" - << abort(FatalError); - } -} - - -template -template -void Foam::BlockMatrixAgglomeration::agglomerateCoeffs -( - const labelList& coeffRestrictAddr, - Field& activeCoarseDiag, - Field& activeCoarseUpper, - const Field& activeFineUpper, - const Field& activeFineUpperTranspose -) const -{ - // Does the matrix have solo equations - bool soloEqns = nSolo_ > 0; - - // Get addressing - const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); - const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); - - forAll(coeffRestrictAddr, fineCoeffI) - { - label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffI]]; - label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffI]]; - - // If the coefficient touches block zero and - // solo equations are present, skip it - if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) - { - continue; - } - - label cCoeff = coeffRestrictAddr[fineCoeffI]; - - if (cCoeff >= 0) - { - activeCoarseUpper[cCoeff] += activeFineUpper[fineCoeffI]; - } - else - { - // Add the fine face coefficient into the diagonal - // Note: upper and lower coeffs are transpose of - // each other. HJ, 28/May/2014 - activeCoarseDiag[-1 - cCoeff] += - activeFineUpper[fineCoeffI] - + activeFineUpperTranspose[fineCoeffI]; - } - } -} - - -template -template -void Foam::BlockMatrixAgglomeration::agglomerateCoeffs -( - const labelList& coeffRestrictAddr, - Field& activeCoarseDiag, - Field& activeCoarseUpper, - const Field& activeFineUpper, - Field& activeCoarseLower, - const Field& activeFineLower -) const -{ - // Does the matrix have solo equations - bool soloEqns = nSolo_ > 0; - - // Get addressing - const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); - const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); - - forAll(coeffRestrictAddr, fineCoeffI) - { - label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffI]]; - label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffI]]; - - // If the coefficient touches block zero and - // solo equations are present, skip it - if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) - { - continue; - } - - label cCoeff = coeffRestrictAddr[fineCoeffI]; - - if (cCoeff >= 0) - { - activeCoarseUpper[cCoeff] += activeFineUpper[fineCoeffI]; - activeCoarseLower[cCoeff] += activeFineLower[fineCoeffI]; - } - else - { - // Add the fine face coefficients into the diagonal. - activeCoarseDiag[-1 - cCoeff] += - activeFineUpper[fineCoeffI] - + activeFineLower[fineCoeffI]; - } - } -} - - -template -void Foam::BlockMatrixAgglomeration::restrictDiagDecoupled -( - const CoeffField& Coeff, - CoeffField& coarseCoeff -) const -{ - typedef CoeffField TypeCoeffField; - - if - ( - Coeff.activeType() == blockCoeffBase::LINEAR - && coarseCoeff.activeType() == blockCoeffBase::LINEAR - ) - { - typedef typename TypeCoeffField::linearType linearType; - typedef typename TypeCoeffField::linearTypeField linearTypeField; - - linearTypeField& activeCoarseCoeff = coarseCoeff.asLinear(); - const linearTypeField& activeCoeff = Coeff.asLinear(); - - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } - - forAll (Coeff, i) - { - activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; - } - } - else if - ( - Coeff.activeType() == blockCoeffBase::SCALAR - && coarseCoeff.activeType() == blockCoeffBase::SCALAR - ) - { - typedef typename TypeCoeffField::scalarType scalarType; - typedef typename TypeCoeffField::scalarTypeField scalarTypeField; - - scalarTypeField& activeCoarseCoeff = coarseCoeff.asScalar(); - const scalarTypeField& activeCoeff = Coeff.asScalar(); - - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } - - forAll (Coeff, i) - { - activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; - } - } - else - { - FatalErrorIn - ( - "void BlockMatrixAgglomeration::restrictDiagDecoupled()" - " const" - ) << "Problem in coeff type morphing" - << abort(FatalError); - } -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -template -Foam::BlockMatrixAgglomeration::BlockMatrixAgglomeration -( - const BlockLduMatrix& matrix, - const dictionary& dict, - const label groupSize, - const label minCoarseEqns -) -: - BlockMatrixCoarsening(matrix, dict, groupSize, minCoarseEqns), - matrix_(matrix), - normPtr_(BlockCoeffNorm::New(dict)), - agglomIndex_(matrix_.lduAddr().size()), - groupSize_(groupSize), - nSolo_(0), - nCoarseEqns_(0), - coarsen_(false), - lTime_() -{ - calcAgglomeration(); -} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -template -Foam::BlockMatrixAgglomeration::~BlockMatrixAgglomeration() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -template -Foam::autoPtr > -Foam::BlockMatrixAgglomeration::restrictMatrix() const -{ - if (!coarsen_) - { - FatalErrorIn - ( - "autoPtr > " - "BlockMatrixAgglomeration::restrictMatrix() const" - ) << "Requesting coarse matrix when it cannot be created" - << abort(FatalError); - } - - // Construct the coarse matrix and ldu addressing for the next level - // Algorithm: - // 1) Loop through all fine coeffs. If the agglom labels on two sides are - // different, this creates a coarse coeff. Define owner and neighbour - // for this coeff based on cluster IDs. - // 2) Check if the coeff has been seen before. If yes, add the coefficient - // to the appropriate field (stored with the equation). If no, create - // a new coeff with neighbour ID and add the coefficient - // 3) Once all the coeffs have been created, loop through all clusters and - // insert the coeffs in the upper order. At the same time, collect the - // owner and neighbour addressing. - // 4) Agglomerate the diagonal by summing up the fine diagonal - - // Get addressing - const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); - const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); - - const label nFineCoeffs = upperAddr.size(); - -# ifdef FULLDEBUG - if (agglomIndex_.size() != matrix_.lduAddr().size()) - { - FatalErrorIn - ( - "autoPtr >" - "BlockMatrixAgglomeration::restrictMatrix() const" - ) << "agglomIndex array does not correspond to fine level. " << endl - << " Size: " << agglomIndex_.size() - << " number of equations: " << matrix_.lduAddr().size() - << abort(FatalError); - } -# endif - - // Does the matrix have solo equations - bool soloEqns = nSolo_ > 0; - - // Storage for block neighbours and coefficients - - // Guess initial maximum number of neighbours in block - label maxNnbrs = 10; - - // Number of neighbours per block - labelList blockNnbrs(nCoarseEqns_, 0); - - // Setup initial packed storage for neighbours and coefficients - labelList blockNbrsData(maxNnbrs*nCoarseEqns_); - - // Create face-restriction addressing - labelList coeffRestrictAddr(nFineCoeffs); - - // Initial neighbour array (not in upper-triangle order) - labelList initCoarseNeighb(nFineCoeffs); - - // Counter for coarse coeffs - label nCoarseCoeffs = 0; - - // Loop through all fine coeffs - forAll (upperAddr, fineCoeffi) - { - label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffi]]; - label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffi]]; - - // If the coefficient touches block zero and solo equations are - // present, skip it - if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) - { - continue; - } - - if (rmUpperAddr == rmLowerAddr) - { - // For each fine coeff inside of a coarse cluster keep the address - // of the cluster corresponding to the coeff in the - // coeffRestrictAddr as a negative index - coeffRestrictAddr[fineCoeffi] = -(rmUpperAddr + 1); - } - else - { - // This coeff is a part of a coarse coeff - - label cOwn = rmUpperAddr; - label cNei = rmLowerAddr; - - // Get coarse owner and neighbour - if (rmUpperAddr > rmLowerAddr) - { - cOwn = rmLowerAddr; - cNei = rmUpperAddr; - } - - // Check the neighbour to see if this coeff has already been found - bool nbrFound = false; - label& ccnCoeffs = blockNnbrs[cOwn]; - - for (int i = 0; i < ccnCoeffs; i++) - { - if (initCoarseNeighb[blockNbrsData[maxNnbrs*cOwn + i]] == cNei) - { - nbrFound = true; - coeffRestrictAddr[fineCoeffi] = - blockNbrsData[maxNnbrs*cOwn + i]; - break; - } - } - - if (!nbrFound) - { - if (ccnCoeffs >= maxNnbrs) - { - // Double the size of list and copy data - label oldMaxNnbrs = maxNnbrs; - maxNnbrs *= 2; - - // Resize and copy list - const labelList oldBlockNbrsData = blockNbrsData; - blockNbrsData.setSize(maxNnbrs*nCoarseEqns_); - - forAll (blockNnbrs, i) - { - for (int j = 0; j < blockNnbrs[i]; j++) - { - blockNbrsData[maxNnbrs*i + j] = - oldBlockNbrsData[oldMaxNnbrs*i + j]; - } - } - } - - blockNbrsData[maxNnbrs*cOwn + ccnCoeffs] = nCoarseCoeffs; - initCoarseNeighb[nCoarseCoeffs] = cNei; - coeffRestrictAddr[fineCoeffi] = nCoarseCoeffs; - ccnCoeffs++; - - // New coarse coeff created - nCoarseCoeffs++; - } - } - } // End for all fine coeffs - - - // Renumber into upper-triangular order - - // All coarse owner-neighbour storage - labelList coarseOwner(nCoarseCoeffs); - labelList coarseNeighbour(nCoarseCoeffs); - labelList coarseCoeffMap(nCoarseCoeffs); - - label coarseCoeffi = 0; - - forAll (blockNnbrs, cci) - { - label* cCoeffs = &blockNbrsData[maxNnbrs*cci]; - label ccnCoeffs = blockNnbrs[cci]; - - for (int i = 0; i < ccnCoeffs; i++) - { - coarseOwner[coarseCoeffi] = cci; - coarseNeighbour[coarseCoeffi] = initCoarseNeighb[cCoeffs[i]]; - coarseCoeffMap[cCoeffs[i]] = coarseCoeffi; - coarseCoeffi++; - } - } - - forAll(coeffRestrictAddr, fineCoeffi) - { - label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffi]]; - label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffi]]; - - // If the coefficient touches block zero and solo equations are - // present, skip it - if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) - { - continue; - } - - if (coeffRestrictAddr[fineCoeffi] >= 0) - { - coeffRestrictAddr[fineCoeffi] = - coarseCoeffMap[coeffRestrictAddr[fineCoeffi]]; - } - } - - // Clear the temporary storage for the coarse matrix data - blockNnbrs.setSize(0); - blockNbrsData.setSize(0); - initCoarseNeighb.setSize(0); - coarseCoeffMap.setSize(0); - - - // Create coarse-level coupled interfaces - - // Create coarse interfaces, addressing and coefficients - const typename BlockLduInterfaceFieldPtrsList::Type& - interfaceFields = - const_cast&>(matrix_).interfaces(); - - // Set the coarse interfaces and coefficients - lduInterfacePtrsList coarseInterfaces(interfaceFields.size()); - - labelListList coarseInterfaceAddr(interfaceFields.size()); - - // Add the coarse level - - // Set the coarse ldu addressing onto the list - autoPtr coarseAddrPtr - ( - new lduPrimitiveMesh - ( - nCoarseEqns_, - coarseOwner, - coarseNeighbour, - Pstream::worldComm, //HJ, AMG Comm fineMesh.comm(), - true - ) - ); - - // Initialise transfer of restrict addressing on the interface - // HJ, must use blocking comms. HJ, 9/Jun/2016 - forAll (interfaceFields, intI) - { - if (interfaceFields.set(intI)) - { - interfaceFields[intI].coupledInterface().initInternalFieldTransfer - ( - Pstream::blocking, - agglomIndex_ - ); - } - } - - // Store coefficients to avoid tangled communications - // HJ, 1/Apr/2009 - FieldField fineInterfaceAddr(interfaceFields.size()); - - forAll (interfaceFields, intI) - { - if (interfaceFields.set(intI)) - { - const lduInterface& fineInterface = - interfaceFields[intI].coupledInterface(); - - fineInterfaceAddr.set - ( - intI, - new labelField - ( - fineInterface.internalFieldTransfer - ( - Pstream::blocking, - agglomIndex_ - ) - ) - ); - } - } - - // Create AMG interfaces - forAll (interfaceFields, intI) - { - if (interfaceFields.set(intI)) - { - const lduInterface& fineInterface = - interfaceFields[intI].coupledInterface(); - - coarseInterfaces.set - ( - intI, - AMGInterface::New - ( - coarseAddrPtr(), - coarseInterfaces, - fineInterface, - fineInterface.interfaceInternalField(agglomIndex_), - fineInterfaceAddr[intI] - ).ptr() - ); - } - } - - forAll (interfaceFields, intI) - { - if (interfaceFields.set(intI)) - { - const AMGInterface& coarseInterface = - refCast(coarseInterfaces[intI]); - - coarseInterfaceAddr[intI] = coarseInterface.faceCells(); - } - } - - // Add interfaces - coarseAddrPtr->addInterfaces - ( - coarseInterfaces, - coarseInterfaceAddr, - matrix_.patchSchedule() - ); - - // Set the coarse level matrix - autoPtr > coarseMatrixPtr - ( - new BlockLduMatrix(coarseAddrPtr()) - ); - BlockLduMatrix& coarseMatrix = coarseMatrixPtr(); - - // Get interfaces from coarse matrix - typename BlockLduInterfaceFieldPtrsList::Type& - coarseInterfaceFieldsTransfer = coarseMatrix.interfaces(); - - // Aggolmerate the upper and lower coupled coefficients - forAll (interfaceFields, intI) - { - if (interfaceFields.set(intI)) - { - const AMGInterface& coarseInterface = - refCast(coarseInterfaces[intI]); - - coarseInterfaceFieldsTransfer.set - ( - intI, - BlockAMGInterfaceField::New - ( - coarseInterface, - interfaceFields[intI] - ).ptr() - ); - - // Since the type of agglomeration is now templated, agglomeration - // of block coefficients must be done by a FIELD (not interface) - // via a new set of virtual functions - // HJ, 16/Mar/2016 - - // Note: in the scalar AMG, agglomeration is done by the interface - // (always scalar) but in the block matrix it is done by a - // templated block interface field - // HJ, 16/Mar/2016 - - // Cast the interface into AMG type - const BlockAMGInterfaceField& coarseField = - refCast > - ( - coarseInterfaceFieldsTransfer[intI] - ); - - coarseMatrix.coupleUpper().set - ( - intI, - coarseField.agglomerateBlockCoeffs - ( - matrix_.coupleUpper()[intI] - ) - ); - - coarseMatrix.coupleLower().set - ( - intI, - coarseField.agglomerateBlockCoeffs - ( - matrix_.coupleLower()[intI] - ) - ); - } - } - - // Matrix restriction done! - - typedef CoeffField TypeCoeffField; - - typedef typename TypeCoeffField::squareTypeField squareTypeField; - typedef typename TypeCoeffField::linearTypeField linearTypeField; - typedef typename TypeCoeffField::scalarTypeField scalarTypeField; - - TypeCoeffField& coarseUpper = coarseMatrix.upper(); - TypeCoeffField& coarseDiag = coarseMatrix.diag(); - const TypeCoeffField& fineUpper = matrix_.upper(); - const TypeCoeffField& fineDiag = matrix_.diag(); - - // KRJ: 2013-01-31: Many cases needed as there are different combinations - - // Note: - // In coarsening, the off-diagonal coefficient type should be preserved - // and the case of taking out of off-diag from diag may need to be handled - // separately (expand the off-diag coeff to diag type before removing it - // from the diag coefficient). Since this has not been encountered yet - // only matching diag/off-diag types are handled. - // HJ, 15/Feb/2016 - if (matrix_.symmetric()) - { - if - ( - fineDiag.activeType() == blockCoeffBase::SQUARE - || fineUpper.activeType() == blockCoeffBase::SQUARE - ) - { - squareTypeField& activeCoarseDiag = coarseDiag.asSquare(); - - squareTypeField& activeCoarseUpper = coarseUpper.asSquare(); - const squareTypeField& activeFineUpper = fineUpper.asSquare(); - - // Use lower as transpose of upper - squareTypeField activeFineUpperTranspose = - activeFineUpper.T(); - - restrictDiag(fineDiag, coarseDiag); - - agglomerateCoeffs - ( - coeffRestrictAddr, - activeCoarseDiag, - activeCoarseUpper, - activeFineUpper, - activeFineUpperTranspose - ); - } - else if - ( - fineDiag.activeType() == blockCoeffBase::LINEAR - || fineUpper.activeType() == blockCoeffBase::LINEAR - ) - { - linearTypeField& activeCoarseDiag = coarseDiag.asLinear(); - - linearTypeField& activeCoarseUpper = coarseUpper.asLinear(); - const linearTypeField& activeFineUpper = fineUpper.asLinear(); - - restrictDiag(fineDiag, coarseDiag); - - agglomerateCoeffs - ( - coeffRestrictAddr, - activeCoarseDiag, - activeCoarseUpper, - activeFineUpper, - activeFineUpper - ); - } - else if - ( - fineDiag.activeType() == blockCoeffBase::SCALAR - || fineUpper.activeType() == blockCoeffBase::SCALAR - ) - { - scalarTypeField& activeCoarseDiag = coarseDiag.asScalar(); - - scalarTypeField& activeCoarseUpper = coarseUpper.asScalar(); - const scalarTypeField& activeFineUpper = fineUpper.asScalar(); - - restrictDiag(fineDiag, coarseDiag); - - agglomerateCoeffs - ( - coeffRestrictAddr, - activeCoarseDiag, - activeCoarseUpper, - activeFineUpper, - activeFineUpper - ); - } - else - { - FatalErrorIn - ( - "autoPtr >" - "BlockMatrixAgglomeration::restrictMatrix() const" - ) << "Matrix coeff type morphing error, symmetric matrix" - << abort(FatalError); - } - } - else // asymmetric matrix - { - TypeCoeffField& coarseLower = coarseMatrix.lower(); - const TypeCoeffField& fineLower = matrix_.lower(); - - if - ( - fineDiag.activeType() == blockCoeffBase::SQUARE - || fineUpper.activeType() == blockCoeffBase::SQUARE - ) - { - squareTypeField& activeCoarseDiag = coarseDiag.asSquare(); - - squareTypeField& activeCoarseUpper = coarseUpper.asSquare(); - const squareTypeField& activeFineUpper = fineUpper.asSquare(); - - squareTypeField& activeCoarseLower = coarseLower.asSquare(); - const squareTypeField& activeFineLower = fineLower.asSquare(); - - restrictDiag(fineDiag, coarseDiag); - - agglomerateCoeffs - ( - coeffRestrictAddr, - activeCoarseDiag, - activeCoarseUpper, - activeFineUpper, - activeCoarseLower, - activeFineLower - ); - } - else if - ( - fineDiag.activeType() == blockCoeffBase::LINEAR - || fineUpper.activeType() == blockCoeffBase::LINEAR - ) - { - linearTypeField& activeCoarseDiag = coarseDiag.asLinear(); - - linearTypeField& activeCoarseUpper = coarseUpper.asLinear(); - const linearTypeField& activeFineUpper = fineUpper.asLinear(); - - linearTypeField& activeCoarseLower = coarseLower.asLinear(); - const linearTypeField& activeFineLower = fineLower.asLinear(); - - restrictDiag(fineDiag, coarseDiag); - - agglomerateCoeffs - ( - coeffRestrictAddr, - activeCoarseDiag, - activeCoarseUpper, - activeFineUpper, - activeCoarseLower, - activeFineLower - ); - } - else if - ( - fineDiag.activeType() == blockCoeffBase::SCALAR - || fineUpper.activeType() == blockCoeffBase::SCALAR - ) - { - scalarTypeField& activeCoarseDiag = coarseDiag.asScalar(); - - scalarTypeField& activeCoarseUpper = coarseUpper.asScalar(); - const scalarTypeField& activeFineUpper = fineUpper.asScalar(); - - scalarTypeField& activeCoarseLower = coarseLower.asScalar(); - const scalarTypeField& activeFineLower = fineLower.asScalar(); - - restrictDiag(fineDiag, coarseDiag); - - agglomerateCoeffs - ( - coeffRestrictAddr, - activeCoarseDiag, - activeCoarseUpper, - activeFineUpper, - activeCoarseLower, - activeFineLower - ); - } - else - { - FatalErrorIn - ( - "autoPtr >" - "BlockMatrixAgglomeration::restrictMatrix() const" - ) << "Matrix coeff type morphing error, asymmetric matrix" - << abort(FatalError); - } - } - - // Create and return BlockAMGLevel - return autoPtr > - ( - new coarseBlockAMGLevel - ( - coarseAddrPtr, - coarseMatrixPtr, - this->dict(), - this->type(), - this->groupSize(), - this->minCoarseEqns() - ) - ); -} - - -template -void Foam::BlockMatrixAgglomeration::restrictResidual -( - const Field& res, - Field& coarseRes -) const -{ - coarseRes = pTraits::zero; - - forAll (res, i) - { - coarseRes[agglomIndex_[i]] += res[i]; - } -} - - -template -void Foam::BlockMatrixAgglomeration::prolongateCorrection -( - Field& x, - const Field& coarseX -) const -{ - forAll (x, i) - { - x[i] += coarseX[agglomIndex_[i]]; - } -} - - -// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.H deleted file mode 100644 index e75d523b3..000000000 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.H +++ /dev/null @@ -1,212 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Class - BlockMatrixAgglomeration - -Description - Agglomerative block matrix AMG coarsening - -Author - Klas Jareteg, 2013-04-15 - -SourceFiles - BlockMatrixAgglomeration.C - -\*---------------------------------------------------------------------------*/ - -#ifndef BlockMatrixAgglomeration_H -#define BlockMatrixAgglomeration_H - -#include "BlockMatrixCoarsening.H" -#include "BlockLduMatrix.H" -#include "BlockCoeffNorm.H" -#include "BlockCoeff.H" -#include "tolerancesSwitch.H" -#include "cpuTime.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -/*---------------------------------------------------------------------------*\ - Class BlockMatrixAgglomeration Declaration -\*---------------------------------------------------------------------------*/ - -template -class BlockMatrixAgglomeration -: - public BlockMatrixCoarsening -{ - // Private Data - - //- Reference to matrix - const BlockLduMatrix& matrix_; - - //- Norm calculator - autoPtr > normPtr_; - - //- Child array: for each fine equation give agglomeration cluster - // index - labelField agglomIndex_; - - //- Group size - label groupSize_; - - //- Number of solo cells - label nSolo_; - - //- Number of coarse equations - label nCoarseEqns_; - - //- Can a coarse level be constructed? - bool coarsen_; - - cpuTime lTime_; - - - // Private Member Functions - - //- Disallow default bitwise copy construct - BlockMatrixAgglomeration(const BlockMatrixAgglomeration&); - - // Disallow default bitwise assignment - void operator=(const BlockMatrixAgglomeration&); - - - //- Calculate agglomeration index (child) - void calcAgglomeration(); - - //- Restrict CoeffField. Used for diag coefficient - void restrictDiag - ( - const CoeffField& Coeff, - CoeffField& coarseCoeff - ) const; - - //- Agglomerate coeffs, symmetric matrix - template - void agglomerateCoeffs - ( - const labelList& coeffRestrictAddr, - Field& activeCoarseDiag, - Field& activeCoarseUpper, - const Field& activeFineUpper, - const Field& activeFineUpperTranspose - ) const; - - //- Agglomerate coeffs, assymmetric matrix - template - void agglomerateCoeffs - ( - const labelList& coeffRestrictAddr, - Field& activeCoarseDiag, - Field& activeCoarseUpper, - const Field& activeFineUpper, - Field& activeCoarseLower, - const Field& activeFineLower - ) const; - - //- Restrict CoeffField, decoupled version. Used for diag coefficient - void restrictDiagDecoupled - ( - const CoeffField& Coeff, - CoeffField& coarseCoeff - ) const; - - // Private Static Data - - //- Weighting factor - static const debug::tolerancesSwitch weightFactor_; - - //- Diagonal scaling factor - static const debug::tolerancesSwitch diagFactor_; - - -public: - - //- Runtime type information - TypeName("AAMG"); - - - // Constructors - - //- Construct from matrix and group size - BlockMatrixAgglomeration - ( - const BlockLduMatrix& matrix, - const dictionary& dict, - const label groupSize, - const label minCoarseEqns - ); - - - //- Destructor - virtual ~BlockMatrixAgglomeration(); - - - // Member Functions - - //- Can a coarse level be constructed? - virtual bool coarsen() const - { - return coarsen_; - } - - //- Restrict matrix - virtual autoPtr > restrictMatrix() const; - - //- Restrict residual - virtual void restrictResidual - ( - const Field& res, - Field& coarseRes - ) const; - - //- Prolongate correction - virtual void prolongateCorrection - ( - Field& x, - const Field& coarseX - ) const; -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#ifdef NoRepository -# include "BlockMatrixAgglomeration.C" -#endif - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/blockMatrixAgglomerations.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/blockMatrixAgglomerations.C deleted file mode 100644 index d90e8d9db..000000000 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/blockMatrixAgglomerations.C +++ /dev/null @@ -1,43 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -\*---------------------------------------------------------------------------*/ - -#include "blockMatrixAgglomerations.H" -#include "blockMatrixCoarsenings.H" -#include "coarseBlockAMGLevel.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -makeBlockMatrixCoarsenings(blockMatrixAgglomeration); - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/blockMatrixAgglomerations.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/blockMatrixAgglomerations.H deleted file mode 100644 index b0d6c08fb..000000000 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/blockMatrixAgglomerations.H +++ /dev/null @@ -1,65 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Class - BlockMatrixAgglomeration - -Description - Typedefs for block matrix AMG coarsening. - -Author - Klas Jareteg, 2012-12-15 - -SourceFiles - blockMatrixAgglomeration.C - -\*---------------------------------------------------------------------------*/ - -#ifndef blockMatrixAgglomerations_H -#define blockMatrixAgglomerations_H - -#include "scalarBlockMatrixAgglomeration.H" -#include "tensorBlockMatrixAgglomeration.H" -#include "BlockMatrixAgglomeration.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -typedef BlockMatrixAgglomeration blockMatrixAgglomerationScalar; -typedef BlockMatrixAgglomeration blockMatrixAgglomerationVector; -typedef BlockMatrixAgglomeration blockMatrixAgglomerationTensor; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/scalarBlockMatrixAgglomeration.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/scalarBlockMatrixAgglomeration.H deleted file mode 100644 index 05477172e..000000000 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/scalarBlockMatrixAgglomeration.H +++ /dev/null @@ -1,86 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Class - BlockMatrixAgglomeration - -Description - Specialisation of the BlockMatrixAgglomeration for scalars. - -Author - Klas Jareteg, 2013-01-31 - -\*---------------------------------------------------------------------------*/ - -#ifndef scalarBlockMatrixAgglomeration_H -#define scalarBlockMatrixAgglomeration_H - -#include "blockMatrixCoarsenings.H" -#include "blockMatrixAgglomerations.H" -#include "BlockMatrixAgglomeration.H" -#include "BlockMatrixCoarsening.H" -#include "runTimeSelectionTables.H" -#include "scalarBlockLduMatrix.H" -#include "scalarBlockConstraint.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * Forward declaration of template friend fuctions * * * * * * * // - - -/*---------------------------------------------------------------------------*\ - Class BlockMatrixAgglomeration Declaration -\*---------------------------------------------------------------------------*/ - -// Disable restrict matrix: no square type -template<> -inline autoPtr > -BlockMatrixAgglomeration::restrictMatrix() const -{ - FatalErrorIn - ( - "autoPtr > " - "BlockMatrixAgglomeration::restrictMatrix() const" - ) << "Function not implemented for Type=scalar. " << endl - << abort(FatalError); - - // Dummy return to keep compiler happy - return autoPtr >(NULL); -} - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/tensorBlockMatrixAgglomeration.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/tensorBlockMatrixAgglomeration.H deleted file mode 100644 index 6a19fdeae..000000000 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/tensorBlockMatrixAgglomeration.H +++ /dev/null @@ -1,83 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Class - BlockMatrixAgglomeration - -Description - Specialisation of the BlockMatrixAgglomeration for tensors. - -Author - Klas Jareteg, 2013-01-31 - -\*---------------------------------------------------------------------------*/ - -#ifndef tensorBlockMatrixAgglomeration_H -#define tensorBlockMatrixAgglomeration_H - -#include "blockMatrixCoarsenings.H" -#include "blockMatrixAgglomerations.H" -#include "BlockMatrixAgglomeration.H" -#include "runTimeSelectionTables.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -/*---------------------------------------------------------------------------*\ - Class BlockMatrixAgglomeration Declaration -\*---------------------------------------------------------------------------*/ - -// Disable restrict matrix: no square type -template<> -inline autoPtr > -BlockMatrixAgglomeration::restrictMatrix -() const -{ - FatalErrorIn - ( - "autoPtr > " - "BlockMatrixAgglomeration::" - "restrictMatrix() const" - ) << "Function not implemented for Type=tensor. " << endl - << abort(FatalError); - - // Dummy return to keep compiler happy - return autoPtr >(NULL); -} - - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C index c5752bba9..980d46cb9 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C @@ -555,10 +555,8 @@ void Foam::BlockMatrixClustering::restrictDiag squareTypeField& activeCoarseCoeff = coarseCoeff.asSquare(); const squareTypeField& activeCoeff = Coeff.asSquare(); - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } + // Reset coefficients to zero + activeCoarseCoeff = pTraits::zero; forAll (Coeff, i) { @@ -577,10 +575,8 @@ void Foam::BlockMatrixClustering::restrictDiag linearTypeField& activeCoarseCoeff = coarseCoeff.asLinear(); const linearTypeField& activeCoeff = Coeff.asLinear(); - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } + // Reset coefficients to zero + activeCoarseCoeff = pTraits::zero; forAll (Coeff, i) { @@ -599,10 +595,8 @@ void Foam::BlockMatrixClustering::restrictDiag scalarTypeField& activeCoarseCoeff = coarseCoeff.asScalar(); const scalarTypeField& activeCoeff = Coeff.asScalar(); - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } + // Reset coefficients to zero + activeCoarseCoeff = pTraits::zero; forAll (Coeff, i) { @@ -624,7 +618,6 @@ template template void Foam::BlockMatrixClustering::agglomerateCoeffs ( - const labelList& coeffRestrictAddr, Field& activeCoarseDiag, Field& activeCoarseUpper, const Field& activeFineUpper, @@ -638,7 +631,10 @@ void Foam::BlockMatrixClustering::agglomerateCoeffs const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); - forAll(coeffRestrictAddr, fineCoeffI) + // Reset coefficients to zero. Cannot touch the diagonal + activeCoarseUpper = pTraits::zero; + + forAll(coeffRestrictAddr_, fineCoeffI) { label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffI]]; label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffI]]; @@ -650,7 +646,7 @@ void Foam::BlockMatrixClustering::agglomerateCoeffs continue; } - label cCoeff = coeffRestrictAddr[fineCoeffI]; + label cCoeff = coeffRestrictAddr_[fineCoeffI]; if (cCoeff >= 0) { @@ -673,7 +669,6 @@ template template void Foam::BlockMatrixClustering::agglomerateCoeffs ( - const labelList& coeffRestrictAddr, Field& activeCoarseDiag, Field& activeCoarseUpper, const Field& activeFineUpper, @@ -688,7 +683,11 @@ void Foam::BlockMatrixClustering::agglomerateCoeffs const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); - forAll(coeffRestrictAddr, fineCoeffI) + // Reset coefficients to zero. Cannot touch the diagonal + activeCoarseUpper = pTraits::zero; + activeCoarseLower = pTraits::zero; + + forAll(coeffRestrictAddr_, fineCoeffI) { label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffI]]; label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffI]]; @@ -700,7 +699,7 @@ void Foam::BlockMatrixClustering::agglomerateCoeffs continue; } - label cCoeff = coeffRestrictAddr[fineCoeffI]; + label cCoeff = coeffRestrictAddr_[fineCoeffI]; if (cCoeff >= 0) { @@ -739,10 +738,8 @@ void Foam::BlockMatrixClustering::restrictDiagDecoupled linearTypeField& activeCoarseCoeff = coarseCoeff.asLinear(); const linearTypeField& activeCoeff = Coeff.asLinear(); - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } + // Reset coefficients to zero + activeCoarseCoeff = pTraits::zero; forAll (Coeff, i) { @@ -761,10 +758,8 @@ void Foam::BlockMatrixClustering::restrictDiagDecoupled scalarTypeField& activeCoarseCoeff = coarseCoeff.asScalar(); const scalarTypeField& activeCoeff = Coeff.asScalar(); - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } + // Reset coefficients to zero + activeCoarseCoeff = pTraits::zero; forAll (Coeff, i) { @@ -800,6 +795,7 @@ Foam::BlockMatrixClustering::BlockMatrixClustering maxGroupSize_(readLabel(dict.lookup("maxGroupSize"))), normPtr_(BlockCoeffNorm::New(dict)), agglomIndex_(matrix_.lduAddr().size()), + coeffRestrictAddr_(), groupSize_(groupSize), nSolo_(0), nCoarseEqns_(0), @@ -850,8 +846,6 @@ Foam::BlockMatrixClustering::restrictMatrix() const const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); - const label nFineCoeffs = upperAddr.size(); - # ifdef FULLDEBUG if (agglomIndex_.size() != matrix_.lduAddr().size()) { @@ -866,7 +860,8 @@ Foam::BlockMatrixClustering::restrictMatrix() const } # endif - // Does the matrix have solo equations + // If the matrix will be coarsened, create off-diagonal agglomeration + // Does the matrix have solo equations bool soloEqns = nSolo_ > 0; // Storage for block neighbours and coefficients @@ -880,8 +875,10 @@ Foam::BlockMatrixClustering::restrictMatrix() const // Setup initial packed storage for neighbours and coefficients labelList blockNbrsData(maxNnbrs*nCoarseEqns_); + const label nFineCoeffs = upperAddr.size(); + // Create face-restriction addressing - labelList coeffRestrictAddr(nFineCoeffs); + coeffRestrictAddr_.setSize(nFineCoeffs); // Initial neighbour array (not in upper-triangle order) labelList initCoarseNeighb(nFineCoeffs); @@ -906,8 +903,8 @@ Foam::BlockMatrixClustering::restrictMatrix() const { // For each fine coeff inside of a coarse cluster keep the address // of the cluster corresponding to the coeff in the - // coeffRestrictAddr as a negative index - coeffRestrictAddr[fineCoeffi] = -(rmUpperAddr + 1); + // coeffRestrictAddr_ as a negative index + coeffRestrictAddr_[fineCoeffi] = -(rmUpperAddr + 1); } else { @@ -932,7 +929,7 @@ Foam::BlockMatrixClustering::restrictMatrix() const if (initCoarseNeighb[blockNbrsData[maxNnbrs*cOwn + i]] == cNei) { nbrFound = true; - coeffRestrictAddr[fineCoeffi] = + coeffRestrictAddr_[fineCoeffi] = blockNbrsData[maxNnbrs*cOwn + i]; break; } @@ -962,7 +959,7 @@ Foam::BlockMatrixClustering::restrictMatrix() const blockNbrsData[maxNnbrs*cOwn + ccnCoeffs] = nCoarseCoeffs; initCoarseNeighb[nCoarseCoeffs] = cNei; - coeffRestrictAddr[fineCoeffi] = nCoarseCoeffs; + coeffRestrictAddr_[fineCoeffi] = nCoarseCoeffs; ccnCoeffs++; // New coarse coeff created @@ -995,7 +992,7 @@ Foam::BlockMatrixClustering::restrictMatrix() const } } - forAll(coeffRestrictAddr, fineCoeffi) + forAll(coeffRestrictAddr_, fineCoeffi) { label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffi]]; label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffi]]; @@ -1007,10 +1004,10 @@ Foam::BlockMatrixClustering::restrictMatrix() const continue; } - if (coeffRestrictAddr[fineCoeffi] >= 0) + if (coeffRestrictAddr_[fineCoeffi] >= 0) { - coeffRestrictAddr[fineCoeffi] = - coarseCoeffMap[coeffRestrictAddr[fineCoeffi]]; + coeffRestrictAddr_[fineCoeffi] = + coarseCoeffMap[coeffRestrictAddr_[fineCoeffi]]; } } @@ -1019,8 +1016,6 @@ Foam::BlockMatrixClustering::restrictMatrix() const blockNbrsData.setSize(0); initCoarseNeighb.setSize(0); coarseCoeffMap.setSize(0); - - // Create coarse-level coupled interfaces // Create coarse interfaces, addressing and coefficients @@ -1028,8 +1023,7 @@ Foam::BlockMatrixClustering::restrictMatrix() const const_cast& >(matrix_).interfaces().size(); const typename BlockLduInterfaceFieldPtrsList::Type& - interfaceFields = - const_cast&>(matrix_).interfaces(); + interfaceFields = matrix_.interfaces(); // Set the coarse interfaces and coefficients lduInterfacePtrsList coarseInterfaces(interfaceSize); @@ -1046,7 +1040,7 @@ Foam::BlockMatrixClustering::restrictMatrix() const nCoarseEqns_, coarseOwner, coarseNeighbour, - Pstream::worldComm, //HJ, AMG Comm fineMesh.comm(), + Pstream::worldComm, true ) ); @@ -1140,7 +1134,39 @@ Foam::BlockMatrixClustering::restrictMatrix() const ); BlockLduMatrix& coarseMatrix = coarseMatrixPtr(); + // Build matrix coefficients + this->updateMatrix(coarseMatrix); + + // Create and return BlockAMGLevel + return autoPtr > + ( + new coarseBlockAMGLevel + ( + coarseAddrPtr, + coarseMatrixPtr, + this->dict(), + this->type(), + this->groupSize(), + this->minCoarseEqns() + ) + ); +} + + +template +void Foam::BlockMatrixClustering::updateMatrix +( + BlockLduMatrix& coarseMatrix +) const +{ + // Get interfaces from fine matrix + const typename BlockLduInterfaceFieldPtrsList::Type& + interfaceFields = matrix_.interfaces(); + // Get interfaces from coarse matrix + lduInterfacePtrsList coarseInterfaces = coarseMatrix.mesh().interfaces(); + + // Get interfaces fields from coarse matrix typename BlockLduInterfaceFieldPtrsList::Type& coarseInterfaceFieldsTransfer = coarseMatrix.interfaces(); @@ -1242,7 +1268,6 @@ Foam::BlockMatrixClustering::restrictMatrix() const agglomerateCoeffs ( - coeffRestrictAddr, activeCoarseDiag, activeCoarseUpper, activeFineUpper, @@ -1264,7 +1289,6 @@ Foam::BlockMatrixClustering::restrictMatrix() const agglomerateCoeffs ( - coeffRestrictAddr, activeCoarseDiag, activeCoarseUpper, activeFineUpper, @@ -1286,7 +1310,6 @@ Foam::BlockMatrixClustering::restrictMatrix() const agglomerateCoeffs ( - coeffRestrictAddr, activeCoarseDiag, activeCoarseUpper, activeFineUpper, @@ -1326,7 +1349,6 @@ Foam::BlockMatrixClustering::restrictMatrix() const agglomerateCoeffs ( - coeffRestrictAddr, activeCoarseDiag, activeCoarseUpper, activeFineUpper, @@ -1352,7 +1374,6 @@ Foam::BlockMatrixClustering::restrictMatrix() const agglomerateCoeffs ( - coeffRestrictAddr, activeCoarseDiag, activeCoarseUpper, activeFineUpper, @@ -1378,7 +1399,6 @@ Foam::BlockMatrixClustering::restrictMatrix() const agglomerateCoeffs ( - coeffRestrictAddr, activeCoarseDiag, activeCoarseUpper, activeFineUpper, @@ -1396,20 +1416,6 @@ Foam::BlockMatrixClustering::restrictMatrix() const << abort(FatalError); } } - - // Create and return BlockAMGLevel - return autoPtr > - ( - new coarseBlockAMGLevel - ( - coarseAddrPtr, - coarseMatrixPtr, - this->dict(), - this->type(), - this->groupSize(), - this->minCoarseEqns() - ) - ); } diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.H index 804d323c2..915e5a651 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.H @@ -76,7 +76,10 @@ class BlockMatrixClustering autoPtr > normPtr_; //- Child array: for each fine equation give a clustering index - labelField agglomIndex_; + labelList agglomIndex_; + + //- Face-restriction addressing + mutable labelList coeffRestrictAddr_; //- Group size label groupSize_; @@ -116,7 +119,6 @@ class BlockMatrixClustering template void agglomerateCoeffs ( - const labelList& coeffRestrictAddr, Field& activeCoarseDiag, Field& activeCoarseUpper, const Field& activeFineUpper, @@ -127,7 +129,6 @@ class BlockMatrixClustering template void agglomerateCoeffs ( - const labelList& coeffRestrictAddr, Field& activeCoarseDiag, Field& activeCoarseUpper, const Field& activeFineUpper, @@ -154,7 +155,7 @@ class BlockMatrixClustering public: //- Runtime type information - TypeName("cluster"); + TypeName("AAMG"); // Constructors @@ -197,6 +198,9 @@ public: Field& x, const Field& coarseX ) const; + + //- Update coarse matrix using same coefficients + virtual void updateMatrix(BlockLduMatrix& coarseMatrix) const; }; diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/scalarBlockMatrixClustering.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/scalarBlockMatrixClustering.H index a274599fe..37841e3a4 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/scalarBlockMatrixClustering.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/scalarBlockMatrixClustering.H @@ -57,18 +57,18 @@ namespace Foam // Disable restrict matrix: no square type template<> -inline autoPtr > -BlockMatrixClustering::restrictMatrix() const +inline void +BlockMatrixClustering::updateMatrix +( + BlockLduMatrix& coarseMatrix +) const { FatalErrorIn ( "autoPtr > " - "BlockMatrixClustering::restrictMatrix() const" + "BlockMatrixClustering::updateMatrix(...) const" ) << "Function not implemented for Type=scalar. " << endl << abort(FatalError); - - // Dummy return to keep compiler happy - return autoPtr >(NULL); } diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/tensorBlockMatrixClustering.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/tensorBlockMatrixClustering.H index 5d3f31946..37ba4cf2e 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/tensorBlockMatrixClustering.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/tensorBlockMatrixClustering.H @@ -49,26 +49,23 @@ namespace Foam Class BlockMatrixClustering Declaration \*---------------------------------------------------------------------------*/ -// Disable restrict matrix: no square type +// Disable updateMatrix matrix: no square type template<> -inline autoPtr > -BlockMatrixClustering::restrictMatrix -() const +inline void +BlockMatrixClustering::updateMatrix +( + BlockLduMatrix& coarseMatrix +) const { FatalErrorIn ( - "autoPtr > " - "BlockMatrixClustering::" - "restrictMatrix() const" + "void BlockMatrixClustering::" + "updateMatrix(...) const" ) << "Function not implemented for Type=tensor. " << endl << abort(FatalError); - - // Dummy return to keep compiler happy - return autoPtr >(NULL); } - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixCoarsening/BlockMatrixCoarsening.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixCoarsening/BlockMatrixCoarsening.H index 74b11a148..e957bec1b 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixCoarsening/BlockMatrixCoarsening.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixCoarsening/BlockMatrixCoarsening.H @@ -172,8 +172,7 @@ public: virtual bool coarsen() const = 0; //- Restrict matrix - virtual autoPtr > restrictMatrix - () const = 0; + virtual autoPtr > restrictMatrix() const = 0; //- Restrict residual virtual void restrictResidual @@ -189,6 +188,8 @@ public: const Field& coarseX ) const = 0; + //- Update coarse matrix using same coefficients + virtual void updateMatrix(BlockLduMatrix& coarseMatrix) const = 0; }; @@ -203,15 +204,15 @@ public: #endif #define makeBlockMatrixCoarsening(BlockMatrixCoarseningType, typeBlockMatrixCoarseningType) \ - \ -defineNamedTemplateTypeNameAndDebug(typeBlockMatrixCoarseningType, 0); \ - \ + \ +defineNamedTemplateTypeNameAndDebug(typeBlockMatrixCoarseningType, 0); \ + \ addToRunTimeSelectionTable(BlockMatrixCoarseningType, typeBlockMatrixCoarseningType, matrix); -#define makeBlockMatrixCoarsenings(blockAmgCoarseningType) \ - \ -makeBlockMatrixCoarsening(blockScalarMatrixCoarsening, blockAmgCoarseningType##Scalar); \ -makeBlockMatrixCoarsening(blockVectorMatrixCoarsening, blockAmgCoarseningType##Vector); \ +#define makeBlockMatrixCoarsenings(blockAmgCoarseningType) \ + \ +makeBlockMatrixCoarsening(blockScalarMatrixCoarsening, blockAmgCoarseningType##Scalar); \ +makeBlockMatrixCoarsening(blockVectorMatrixCoarsening, blockAmgCoarseningType##Vector); \ makeBlockMatrixCoarsening(blockTensorMatrixCoarsening, blockAmgCoarseningType##Tensor); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/BlockMatrixReorderedClustering.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/BlockMatrixReorderedClustering.C deleted file mode 100644 index 5126e480a..000000000 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/BlockMatrixReorderedClustering.C +++ /dev/null @@ -1,1505 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Class - BlockMatrixReorderedClustering - -Description - Block matrix AMG coarsening by Jasak clustering algorithm - -Author - Hrvoje Jasak, Wikki Ltd. All rights reserved - -\*---------------------------------------------------------------------------*/ - -#include "BlockMatrixReorderedClustering.H" -#include "boolList.H" -#include "tolerancesSwitch.H" -#include "coeffFields.H" -#include "BlockAMGInterfaceField.H" -#include "coarseBlockAMGLevel.H" -#include "PriorityList.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -template -const Foam::debug::tolerancesSwitch -Foam::BlockMatrixReorderedClustering::weightFactor_ -( - "aamgWeightFactor", - 0.65 -); - - -template -const Foam::debug::tolerancesSwitch -Foam::BlockMatrixReorderedClustering::diagFactor_ -( - "aamgDiagFactor", - 1e-8 -); - - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -template -void Foam::BlockMatrixReorderedClustering::calcClustering() -{ - if (matrix_.diagonal()) - { - // Diag only matrix. Reset and return - agglomIndex_ = 0; - nCoarseEqns_ = 1; - - return; - } - - // Algorithm: - // 1) Calculate appropriate connection strength for symmetric and asymmetric - // matrix - // 2) Collect solo (disconnected) cells and remove them from agglomeration - // by placing them into cluster zero (solo group) - // 3) Loop through all equations and for each equation find the best fit - // neighbour. If all neighbours are grouped, add equation - // to best group - - // Initialise child array - agglomIndex_ = -1; - - const label nRows = matrix_.lduAddr().size(); - - // Get matrix addressing - const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); - const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); - const unallocLabelList& losortAddr = matrix_.lduAddr().losortAddr(); - - const unallocLabelList& ownerStartAddr = - matrix_.lduAddr().ownerStartAddr(); - const unallocLabelList& losortStartAddr = - matrix_.lduAddr().losortStartAddr(); - - - // Calculate clustering - - // Get matrix coefficients and norms Note: the norm - // may be signed, ie. it will take the sign of the coefficient - - const CoeffField& diag = matrix_.diag(); - scalarField normDiag(diag.size()); - normPtr_->normalize(normDiag, diag); - - scalarField normUpper(upperAddr.size(), 0); - scalarField normLower(upperAddr.size(), 0); - - // Note: - // Matrix properties are no longer assumed, eg. the diag and off-diag - // sign is checked - // HJ, 29/Mar/2017 - - // Note: negative connections are eliminated in max(...) below - // HJ, 30/Mar/2017 - - if (matrix_.thereIsUpper()) - { - normPtr_->normalize(normUpper, matrix_.upper()); - - // Owner: upper triangle - forAll (lowerAddr, coeffI) - { - // Sign of strong positive upper is opposite of the sign of - // its diagonal coefficient - normUpper[coeffI] = Foam::max - ( - -1*sign(normDiag[lowerAddr[coeffI]])*normUpper[coeffI], - 0 - ); - } - } - - if (matrix_.thereIsLower()) - { - normPtr_->normalize(normLower, matrix_.lower()); - - // Neighbour: lower triangle - forAll (lowerAddr, coeffI) - { - // Sign of strong positive upper is opposite of the sign of - // its diagonal coefficient - normLower[coeffI] = Foam::max - ( - -1*sign(normDiag[upperAddr[coeffI]])*normLower[coeffI], - 0 - ); - } - } - else - { - normLower = normUpper; - } - - // Take magnitude of the diagonal norm after the normalisation - // of upper and lower is complete - // Note: upper and lower are already normalised - normDiag = mag(normDiag); - - labelList sizeOfGroups(nRows, 0); - - nCoarseEqns_ = 0; - - // 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 - - scalarField scaledNormDiag = diagFactor_()*normDiag; - - boolList zeroCluster(normDiag.size(), true); - - if (matrix_.symmetric()) - { - // Owner: upper triangle - forAll (lowerAddr, coeffI) - { - if (normUpper[coeffI] > scaledNormDiag[lowerAddr[coeffI]]) - { - zeroCluster[lowerAddr[coeffI]] = false; - } - } - - // Neighbour: lower triangle with symm coefficients - forAll (upperAddr, coeffI) - { - if (normUpper[coeffI] > scaledNormDiag[upperAddr[coeffI]]) - { - zeroCluster[upperAddr[coeffI]] = false; - } - } - } - else if (matrix_.asymmetric()) - { - // Owner: upper triangle - forAll (lowerAddr, coeffI) - { - if (normUpper[coeffI] > scaledNormDiag[lowerAddr[coeffI]]) - { - zeroCluster[lowerAddr[coeffI]] = false; - } - } - - // Neighbour: lower triangle with lower coeffs - forAll (upperAddr, coeffI) - { - if (normLower[coeffI] > scaledNormDiag[upperAddr[coeffI]]) - { - zeroCluster[upperAddr[coeffI]] = false; - } - } - } - - // Collect solo equations - forAll (zeroCluster, eqnI) - { - if (zeroCluster[eqnI]) - { - // Found solo equation - nSolo_++; - - agglomIndex_[eqnI] = nCoarseEqns_; - } - } - - if (nSolo_ > 0) - { - // Found solo equations - sizeOfGroups[nCoarseEqns_] = nSolo_; - - nCoarseEqns_++; - } - } - - - // Loop through cells - // - if the cell is not already grouped, open a new group (seed) - // - find stroungest grouped and ungrouped connection: - // - grouped connection is towards an already grouped cell - // - ungrouped connection is towards an ungrouped cell - // - decide if a grouped or ungrouped connection is better - - label curEqn, indexUngrouped, indexGrouped, colI, groupPassI, - nextUngrouped, nextGrouped; - - scalar curWeight, weightUngrouped, weightGrouped; - - for (label rowI = 0; rowI < nRows; rowI++) - { - if (agglomIndex_[rowI] == -1) - { - // Found new ungrouped equation - curEqn = rowI; - - // Reset grouped and upgrouped index - indexUngrouped = -1; - indexGrouped = -1; - - // Make next group (coarse equation) and search for neighbours - agglomIndex_[curEqn] = nCoarseEqns_; - - // Work on the group until the min group size is satisfied - // As each new element of the group is found, group search starts - // from this element until group cluster size is satisfied - for - ( - groupPassI = 1; - groupPassI < minGroupSize_; - groupPassI++ - ) - { - weightUngrouped = 0; - weightGrouped = 0; - - indexUngrouped = -1; - indexGrouped = -1; - - nextUngrouped = -1; - nextGrouped = -1; - - // Visit all neighbour equations - - // Visit upper triangle coefficients from the equation - for - ( - label rowCoeffI = ownerStartAddr[curEqn]; - rowCoeffI < ownerStartAddr[curEqn + 1]; - rowCoeffI++ - ) - { - // Get column index, upper triangle - colI = upperAddr[rowCoeffI]; - - // curWeight = normUpper[rowCoeffI]/ - // // normDiag[curEqn]; - // max(normDiag[curEqn], normDiag[colI]); - - curWeight = Foam::min - ( - normUpper[rowCoeffI]/normDiag[curEqn], - normLower[rowCoeffI]/normDiag[colI] - ); - - if (agglomIndex_[colI] == -1) - { - // Found ungrouped neighbour - if - ( - indexUngrouped == -1 - || curWeight > weightUngrouped - ) - { - // Found or updated ungrouped neighbour - indexUngrouped = rowCoeffI; - weightUngrouped = curWeight; - - // Record possible next equation for search - nextUngrouped = colI; - } - } - else if (agglomIndex_[curEqn] != agglomIndex_[colI]) - { - // Check for neighbour in solo group - if (nSolo_ == 0 || agglomIndex_[colI] != 0) - { - // Found neighbour belonging to other group - if (indexGrouped == -1 || curWeight > weightGrouped) - { - // Found or updated grouped neighbour - indexGrouped = rowCoeffI; - weightGrouped = curWeight; - - // Record possible next equation for search - nextGrouped = colI; - } - } - } - } - - // Visit lower triangle coefficients from the equation - for - ( - label rowCoeffI = losortStartAddr[curEqn]; - rowCoeffI < losortStartAddr[curEqn + 1]; - rowCoeffI++ - ) - { - // Get column index, lower triangle - colI = lowerAddr[losortAddr[rowCoeffI]]; - - // curWeight = normLower[losortAddr[rowCoeffI]]/ - // // normDiag[curEqn]; - // max(normDiag[curEqn], normDiag[colI]); - - curWeight = Foam::min - ( - normLower[losortAddr[rowCoeffI]]/normDiag[curEqn], - normUpper[losortAddr[rowCoeffI]]/normDiag[colI] - ); - - if (agglomIndex_[colI] == -1) - { - // Found first or better ungrouped neighbour - if - ( - indexUngrouped == -1 - || curWeight > weightUngrouped - ) - { - // Found or updated ungrouped neighbour - indexUngrouped = rowCoeffI; - weightUngrouped = curWeight; - - // Record possible next equation for search - nextUngrouped = colI; - } - } - else if (agglomIndex_[curEqn] != agglomIndex_[colI]) - { - // Check for neighbour in solo group - if (nSolo_ == 0 || agglomIndex_[colI] != 0) - { - // Found first of better neighbour belonging to - // other group - if - ( - indexGrouped == -1 - || curWeight > weightGrouped - ) - { - // Found or updated grouped neighbour - indexGrouped = rowCoeffI; - weightGrouped = curWeight; - - // Record possible next equation for search - nextGrouped = colI; - } - } - } - } - - // Decide what to do with current equation - - // If ungrouped candidate exists and it is stronger - // than best weighted grouped candidate, use it - if - ( - indexUngrouped != -1 - && ( - indexGrouped == -1 - || weightUngrouped >= weightFactor_()*weightGrouped - ) - ) - { - // Found new element of group. Add it and use as - // start of next search - - agglomIndex_[nextUngrouped] = agglomIndex_[curEqn]; - sizeOfGroups[agglomIndex_[curEqn]]++; - - // Search from nextEqn - curEqn = nextUngrouped; - } - else - { - // Group full or cannot be extended with new - // ungrouped candidates - break; - } - } - - - // Finished group passes - if - ( - groupPassI > 1 - || indexGrouped == -1 - || sizeOfGroups[agglomIndex_[nextGrouped]] >= maxGroupSize_ - ) - { - // If this is a solo cell and a group is available - // force it into the group irrespective of size - if - ( - sizeOfGroups[agglomIndex_[rowI]] == 0 - && indexGrouped != -1 - - ) - { - // Group exists, but it's too big. Add it anyway - agglomIndex_[rowI] = agglomIndex_[nextGrouped]; - sizeOfGroups[agglomIndex_[nextGrouped]]++; - } - else - { - // No group and no solo group. Make its own group - sizeOfGroups[agglomIndex_[rowI]]++; - nCoarseEqns_++; - } - } - else - { - // Dump current cell into the best group available - agglomIndex_[rowI] = agglomIndex_[nextGrouped]; - sizeOfGroups[agglomIndex_[nextGrouped]]++; - } - } - } - - // Resize group size count - sizeOfGroups.setSize(nCoarseEqns_); - - // Go through faces belonging to a single cluster and sum up their - // off-Diagonals - this will go into Priority list for reordering - - // Weight array - scalarField clusterWeight(nCoarseEqns_, 0); - - for (label eqn = 0; eqn < nRows; eqn++) - { - // Find neighbours of cells which are shared with a cell in another - // cluster - this is a new face on the coarse level - for (label i = ownerStartAddr[eqn]; i < ownerStartAddr[eqn + 1]; i++) - { - // column of upper coefficient (i.e. neighbour) - label j = upperAddr[i]; - - // Check whether the neighbour is in my cluster - if (agglomIndex_[eqn] != agglomIndex_[j]) - { - // Neighbouring cluster found - new face on coarse level - - // Upper triangle - give weight to me - clusterWeight[agglomIndex_[eqn]] += normUpper[i]; - - // Lower triangle - give weight to neighbour - clusterWeight[agglomIndex_[j]] += normLower[i]; - } - } - } - - // Calculate the sum of off-diags for coarse eqns and put them into priority - // array - // Be careful to match the index of the equation to the coarse level] - - // Create the PriorityList, where the weights are sumMagOffDiagCoeffs - PriorityList cardinality(nCoarseEqns_); - - for (label eqnI = 0; eqnI < nCoarseEqns_; ++eqnI) - { - cardinality.set(eqnI, clusterWeight[eqnI]); - } - - // Reorder my coarse equations according to priority array and save into - // renumberedCluster - labelList renumberedCluster(nCoarseEqns_, -1); - - label count = 0; - while (!cardinality.empty()) - { - renumberedCluster[count] = cardinality.removeHead(); - count++; - } - - // Write the new ordering of clusters into agglomIndex_ array - forAll (agglomIndex_, i) - { - agglomIndex_[i] = renumberedCluster[agglomIndex_[i]]; - } - - // The decision on parallel agglomeration needs to be made for the - // whole gang of processes; otherwise I may end up with a different - // number of agglomeration levels on different processors. - - // If the number of coarse equations is less than minimum and - // if the matrix has reduced in size by at least 1/3, coarsen - if - ( - nCoarseEqns_ > BlockMatrixCoarsening::minCoarseEqns() - && 3*nCoarseEqns_ <= 2*nRows - ) - { - coarsen_ = true; - } - - reduce(coarsen_, andOp()); - - if (blockLduMatrix::debug >= 3) - { - // Count singleton clusters - label nSingleClusters = 0; - - // Adjust start based on solo cells cluster status - label start = 0; - - if (nSolo_ > 0) - { - start = 1; - } - - for (label i = start; i < sizeOfGroups.size(); i++) - { - if (sizeOfGroups[i] == 1) - { - nSingleClusters++; - } - } - - Pout<< "Coarse level size: " << nCoarseEqns_ - << " nSolo = " << nSolo_ - << " cluster (" << min(sizeOfGroups) << " " - << max(sizeOfGroups) - << "). N singleton clusters = " << nSingleClusters; - - if (coarsen_) - { - Pout << ". Accepted" << endl; - } - else - { - Pout << ". Rejected" << endl; - } - } -} - - -template -void Foam::BlockMatrixReorderedClustering::restrictDiag -( - const CoeffField& Coeff, - CoeffField& coarseCoeff -) const -{ - typedef CoeffField TypeCoeffField; - - if - ( - Coeff.activeType() == blockCoeffBase::SQUARE - && coarseCoeff.activeType() == blockCoeffBase::SQUARE - ) - { - typedef typename TypeCoeffField::squareType squareType; - typedef typename TypeCoeffField::squareTypeField squareTypeField; - - squareTypeField& activeCoarseCoeff = coarseCoeff.asSquare(); - const squareTypeField& activeCoeff = Coeff.asSquare(); - - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } - - forAll (Coeff, i) - { - activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; - } - } - else if - ( - Coeff.activeType() == blockCoeffBase::LINEAR - && coarseCoeff.activeType() == blockCoeffBase::LINEAR - ) - { - typedef typename TypeCoeffField::linearType linearType; - typedef typename TypeCoeffField::linearTypeField linearTypeField; - - linearTypeField& activeCoarseCoeff = coarseCoeff.asLinear(); - const linearTypeField& activeCoeff = Coeff.asLinear(); - - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } - - forAll (Coeff, i) - { - activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; - } - } - else if - ( - Coeff.activeType() == blockCoeffBase::SCALAR - && coarseCoeff.activeType() == blockCoeffBase::SCALAR - ) - { - typedef typename TypeCoeffField::scalarType scalarType; - typedef typename TypeCoeffField::scalarTypeField scalarTypeField; - - scalarTypeField& activeCoarseCoeff = coarseCoeff.asScalar(); - const scalarTypeField& activeCoeff = Coeff.asScalar(); - - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } - - forAll (Coeff, i) - { - activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; - } - } - else - { - FatalErrorIn - ( - "void BlockMatrixReorderedClustering::restrictDiag() const" - ) << "Problem in coeff type morphing" - << abort(FatalError); - } -} - - -template -template -void Foam::BlockMatrixReorderedClustering::agglomerateCoeffs -( - const labelList& coeffRestrictAddr, - Field& activeCoarseDiag, - Field& activeCoarseUpper, - const Field& activeFineUpper, - const Field& activeFineUpperTranspose -) const -{ - // Does the matrix have solo equations - bool soloEqns = nSolo_ > 0; - - // Get addressing - const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); - const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); - - forAll(coeffRestrictAddr, fineCoeffI) - { - label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffI]]; - label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffI]]; - - // If the coefficient touches block zero and - // solo equations are present, skip it - if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) - { - continue; - } - - label cCoeff = coeffRestrictAddr[fineCoeffI]; - - if (cCoeff >= 0) - { - activeCoarseUpper[cCoeff] += activeFineUpper[fineCoeffI]; - } - else - { - // Add the fine face coefficient into the diagonal - // Note: upper and lower coeffs are transpose of - // each other. HJ, 28/May/2014 - activeCoarseDiag[-1 - cCoeff] += - activeFineUpper[fineCoeffI] - + activeFineUpperTranspose[fineCoeffI]; - } - } -} - - -template -template -void Foam::BlockMatrixReorderedClustering::agglomerateCoeffs -( - const labelList& coeffRestrictAddr, - Field& activeCoarseDiag, - Field& activeCoarseUpper, - const Field& activeFineUpper, - Field& activeCoarseLower, - const Field& activeFineLower -) const -{ - // Does the matrix have solo equations - bool soloEqns = nSolo_ > 0; - - // Get addressing - const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); - const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); - - forAll(coeffRestrictAddr, fineCoeffI) - { - label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffI]]; - label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffI]]; - - // If the coefficient touches block zero and - // solo equations are present, skip it - if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) - { - continue; - } - - label cCoeff = coeffRestrictAddr[fineCoeffI]; - - if (cCoeff >= 0) - { - activeCoarseUpper[cCoeff] += activeFineUpper[fineCoeffI]; - activeCoarseLower[cCoeff] += activeFineLower[fineCoeffI]; - } - else - { - // Add the fine face coefficients into the diagonal. - activeCoarseDiag[-1 - cCoeff] += - activeFineUpper[fineCoeffI] - + activeFineLower[fineCoeffI]; - } - } -} - - -template -void Foam::BlockMatrixReorderedClustering::restrictDiagDecoupled -( - const CoeffField& Coeff, - CoeffField& coarseCoeff -) const -{ - typedef CoeffField TypeCoeffField; - - if - ( - Coeff.activeType() == blockCoeffBase::LINEAR - && coarseCoeff.activeType() == blockCoeffBase::LINEAR - ) - { - typedef typename TypeCoeffField::linearType linearType; - typedef typename TypeCoeffField::linearTypeField linearTypeField; - - linearTypeField& activeCoarseCoeff = coarseCoeff.asLinear(); - const linearTypeField& activeCoeff = Coeff.asLinear(); - - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } - - forAll (Coeff, i) - { - activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; - } - } - else if - ( - Coeff.activeType() == blockCoeffBase::SCALAR - && coarseCoeff.activeType() == blockCoeffBase::SCALAR - ) - { - typedef typename TypeCoeffField::scalarType scalarType; - typedef typename TypeCoeffField::scalarTypeField scalarTypeField; - - scalarTypeField& activeCoarseCoeff = coarseCoeff.asScalar(); - const scalarTypeField& activeCoeff = Coeff.asScalar(); - - forAll (coarseCoeff, i) - { - activeCoarseCoeff[i] = pTraits::zero; - } - - forAll (Coeff, i) - { - activeCoarseCoeff[agglomIndex_[i]] += activeCoeff[i]; - } - } - else - { - FatalErrorIn - ( - "void BlockMatrixReorderedClustering::restrictDiagDecoupled()" - " const" - ) << "Problem in coeff type morphing" - << abort(FatalError); - } -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -template -Foam::BlockMatrixReorderedClustering::BlockMatrixReorderedClustering -( - const BlockLduMatrix& matrix, - const dictionary& dict, - const label groupSize, - const label minCoarseEqns -) -: - BlockMatrixCoarsening(matrix, dict, groupSize, minCoarseEqns), - matrix_(matrix), - minGroupSize_(readLabel(dict.lookup("minGroupSize"))), - maxGroupSize_(readLabel(dict.lookup("maxGroupSize"))), - normPtr_(BlockCoeffNorm::New(dict)), - agglomIndex_(matrix_.lduAddr().size()), - groupSize_(groupSize), - nSolo_(0), - nCoarseEqns_(0), - coarsen_(false), - lTime_() -{ - calcClustering(); -} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -template -Foam::BlockMatrixReorderedClustering::~BlockMatrixReorderedClustering() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -template -Foam::autoPtr > -Foam::BlockMatrixReorderedClustering::restrictMatrix() const -{ - if (!coarsen_) - { - FatalErrorIn - ( - "autoPtr > " - "BlockMatrixReorderedClustering::restrictMatrix() const" - ) << "Requesting coarse matrix when it cannot be created" - << abort(FatalError); - } - - // Construct the coarse matrix and ldu addressing for the next level - // Algorithm: - // 1) Loop through all fine coeffs. If the agglom labels on two sides are - // different, this creates a coarse coeff. Define owner and neighbour - // for this coeff based on cluster IDs. - // 2) Check if the coeff has been seen before. If yes, add the coefficient - // to the appropriate field (stored with the equation). If no, create - // a new coeff with neighbour ID and add the coefficient - // 3) Once all the coeffs have been created, loop through all clusters and - // insert the coeffs in the upper order. At the same time, collect the - // owner and neighbour addressing. - // 4) Agglomerate the diagonal by summing up the fine diagonal - - // Get addressing - const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); - const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); - - const label nFineCoeffs = upperAddr.size(); - -# ifdef FULLDEBUG - if (agglomIndex_.size() != matrix_.lduAddr().size()) - { - FatalErrorIn - ( - "autoPtr >" - "BlockMatrixReorderedClustering::restrictMatrix() const" - ) << "agglomIndex array does not correspond to fine level. " << endl - << " Size: " << agglomIndex_.size() - << " number of equations: " << matrix_.lduAddr().size() - << abort(FatalError); - } -# endif - - // Does the matrix have solo equations - bool soloEqns = nSolo_ > 0; - - // Storage for block neighbours and coefficients - - // Guess initial maximum number of neighbours in block - label maxNnbrs = 10; - - // Number of neighbours per block - labelList blockNnbrs(nCoarseEqns_, 0); - - // Setup initial packed storage for neighbours and coefficients - labelList blockNbrsData(maxNnbrs*nCoarseEqns_); - - // Create face-restriction addressing - labelList coeffRestrictAddr(nFineCoeffs); - - // Initial neighbour array (not in upper-triangle order) - labelList initCoarseNeighb(nFineCoeffs); - - // Counter for coarse coeffs - label nCoarseCoeffs = 0; - - // Loop through all fine coeffs - forAll (upperAddr, fineCoeffi) - { - label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffi]]; - label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffi]]; - - // If the coefficient touches block zero and solo equations are - // present, skip it - if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) - { - continue; - } - - if (rmUpperAddr == rmLowerAddr) - { - // For each fine coeff inside of a coarse cluster keep the address - // of the cluster corresponding to the coeff in the - // coeffRestrictAddr as a negative index - coeffRestrictAddr[fineCoeffi] = -(rmUpperAddr + 1); - } - else - { - // This coeff is a part of a coarse coeff - - label cOwn = rmUpperAddr; - label cNei = rmLowerAddr; - - // Get coarse owner and neighbour - if (rmUpperAddr > rmLowerAddr) - { - cOwn = rmLowerAddr; - cNei = rmUpperAddr; - } - - // Check the neighbour to see if this coeff has already been found - bool nbrFound = false; - label& ccnCoeffs = blockNnbrs[cOwn]; - - for (int i = 0; i < ccnCoeffs; i++) - { - if (initCoarseNeighb[blockNbrsData[maxNnbrs*cOwn + i]] == cNei) - { - nbrFound = true; - coeffRestrictAddr[fineCoeffi] = - blockNbrsData[maxNnbrs*cOwn + i]; - break; - } - } - - if (!nbrFound) - { - if (ccnCoeffs >= maxNnbrs) - { - // Double the size of list and copy data - label oldMaxNnbrs = maxNnbrs; - maxNnbrs *= 2; - - // Resize and copy list - const labelList oldBlockNbrsData = blockNbrsData; - blockNbrsData.setSize(maxNnbrs*nCoarseEqns_); - - forAll (blockNnbrs, i) - { - for (int j = 0; j < blockNnbrs[i]; j++) - { - blockNbrsData[maxNnbrs*i + j] = - oldBlockNbrsData[oldMaxNnbrs*i + j]; - } - } - } - - blockNbrsData[maxNnbrs*cOwn + ccnCoeffs] = nCoarseCoeffs; - initCoarseNeighb[nCoarseCoeffs] = cNei; - coeffRestrictAddr[fineCoeffi] = nCoarseCoeffs; - ccnCoeffs++; - - // New coarse coeff created - nCoarseCoeffs++; - } - } - } // End for all fine coeffs - - - // Renumber into upper-triangular order - - // All coarse owner-neighbour storage - labelList coarseOwner(nCoarseCoeffs); - labelList coarseNeighbour(nCoarseCoeffs); - labelList coarseCoeffMap(nCoarseCoeffs); - - label coarseCoeffi = 0; - - forAll (blockNnbrs, cci) - { - label* cCoeffs = &blockNbrsData[maxNnbrs*cci]; - label ccnCoeffs = blockNnbrs[cci]; - - for (int i = 0; i < ccnCoeffs; i++) - { - coarseOwner[coarseCoeffi] = cci; - coarseNeighbour[coarseCoeffi] = initCoarseNeighb[cCoeffs[i]]; - coarseCoeffMap[cCoeffs[i]] = coarseCoeffi; - coarseCoeffi++; - } - } - - forAll(coeffRestrictAddr, fineCoeffi) - { - label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffi]]; - label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffi]]; - - // If the coefficient touches block zero and solo equations are - // present, skip it - if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0)) - { - continue; - } - - if (coeffRestrictAddr[fineCoeffi] >= 0) - { - coeffRestrictAddr[fineCoeffi] = - coarseCoeffMap[coeffRestrictAddr[fineCoeffi]]; - } - } - - // Clear the temporary storage for the coarse matrix data - blockNnbrs.setSize(0); - blockNbrsData.setSize(0); - initCoarseNeighb.setSize(0); - coarseCoeffMap.setSize(0); - - - // Create coarse-level coupled interfaces - - // Create coarse interfaces, addressing and coefficients - const label interfaceSize = - const_cast& >(matrix_).interfaces().size(); - - const typename BlockLduInterfaceFieldPtrsList::Type& - interfaceFields = - const_cast&>(matrix_).interfaces(); - - // Set the coarse interfaces and coefficients - lduInterfacePtrsList coarseInterfaces(interfaceSize); - - labelListList coarseInterfaceAddr(interfaceSize); - - // Add the coarse level - - // Set the coarse ldu addressing onto the list - autoPtr coarseAddrPtr - ( - new lduPrimitiveMesh - ( - nCoarseEqns_, - coarseOwner, - coarseNeighbour, - Pstream::worldComm, //HJ, AMG Comm fineMesh.comm(), - true - ) - ); - - // Initialise transfer of restrict addressing on the interface - // HJ, must use blocking comms. HJ, 9/Jun/2016 - forAll (interfaceFields, intI) - { - if (interfaceFields.set(intI)) - { - interfaceFields[intI].coupledInterface().initInternalFieldTransfer - ( - Pstream::blocking, - agglomIndex_ - ); - } - } - - // Store coefficients to avoid tangled communications - // HJ, 1/Apr/2009 - FieldField fineInterfaceAddr(interfaceFields.size()); - - forAll (interfaceFields, intI) - { - if (interfaceFields.set(intI)) - { - const lduInterface& fineInterface = - interfaceFields[intI].coupledInterface(); - - fineInterfaceAddr.set - ( - intI, - new labelField - ( - fineInterface.internalFieldTransfer - ( - Pstream::blocking, - agglomIndex_ - ) - ) - ); - } - } - - // Create AMG interfaces - forAll (interfaceFields, intI) - { - if (interfaceFields.set(intI)) - { - const lduInterface& fineInterface = - interfaceFields[intI].coupledInterface(); - - coarseInterfaces.set - ( - intI, - AMGInterface::New - ( - coarseAddrPtr(), - coarseInterfaces, - fineInterface, - fineInterface.interfaceInternalField(agglomIndex_), - fineInterfaceAddr[intI] - ).ptr() - ); - } - } - - forAll (interfaceFields, intI) - { - if (interfaceFields.set(intI)) - { - const AMGInterface& coarseInterface = - refCast(coarseInterfaces[intI]); - - coarseInterfaceAddr[intI] = coarseInterface.faceCells(); - } - } - - // Add interfaces - coarseAddrPtr->addInterfaces - ( - coarseInterfaces, - coarseInterfaceAddr, - matrix_.patchSchedule() - ); - - // Set the coarse level matrix - autoPtr > coarseMatrixPtr - ( - new BlockLduMatrix(coarseAddrPtr()) - ); - BlockLduMatrix& coarseMatrix = coarseMatrixPtr(); - - // Get interfaces from coarse matrix - typename BlockLduInterfaceFieldPtrsList::Type& - coarseInterfaceFieldsTransfer = coarseMatrix.interfaces(); - - // Aggolmerate the upper and lower coupled coefficients - forAll (interfaceFields, intI) - { - if (interfaceFields.set(intI)) - { - const AMGInterface& coarseInterface = - refCast(coarseInterfaces[intI]); - - coarseInterfaceFieldsTransfer.set - ( - intI, - BlockAMGInterfaceField::New - ( - coarseInterface, - interfaceFields[intI] - ).ptr() - ); - - // Since the type of clustering is now templated, clustering - // of block coefficients must be done by a FIELD (not interface) - // via a new set of virtual functions - // HJ, 16/Mar/2016 - - // Note: in the scalar AMG, clustering is done by the interface - // (always scalar) but in the block matrix it is done by a - // templated block interface field - // HJ, 16/Mar/2016 - - // Cast the interface into AMG type - const BlockAMGInterfaceField& coarseField = - refCast > - ( - coarseInterfaceFieldsTransfer[intI] - ); - - coarseMatrix.coupleUpper().set - ( - intI, - coarseField.agglomerateBlockCoeffs - ( - matrix_.coupleUpper()[intI] - ) - ); - - coarseMatrix.coupleLower().set - ( - intI, - coarseField.agglomerateBlockCoeffs - ( - matrix_.coupleLower()[intI] - ) - ); - } - } - - // Matrix restriction done! - - typedef CoeffField TypeCoeffField; - - typedef typename TypeCoeffField::squareTypeField squareTypeField; - typedef typename TypeCoeffField::linearTypeField linearTypeField; - typedef typename TypeCoeffField::scalarTypeField scalarTypeField; - - TypeCoeffField& coarseUpper = coarseMatrix.upper(); - TypeCoeffField& coarseDiag = coarseMatrix.diag(); - const TypeCoeffField& fineUpper = matrix_.upper(); - const TypeCoeffField& fineDiag = matrix_.diag(); - - // KRJ: 2013-01-31: Many cases needed as there are different combinations - - // Note: - // In coarsening, the off-diagonal coefficient type should be preserved - // and the case of taking out of off-diag from diag may need to be handled - // separately (expand the off-diag coeff to diag type before removing it - // from the diag coefficient). Since this has not been encountered yet - // only matching diag/off-diag types are handled. - // HJ, 15/Feb/2016 - if (matrix_.symmetric()) - { - if - ( - fineDiag.activeType() == blockCoeffBase::SQUARE - || fineUpper.activeType() == blockCoeffBase::SQUARE - ) - { - squareTypeField& activeCoarseDiag = coarseDiag.asSquare(); - - squareTypeField& activeCoarseUpper = coarseUpper.asSquare(); - const squareTypeField& activeFineUpper = fineUpper.asSquare(); - - // Use lower as transpose of upper - squareTypeField activeFineUpperTranspose = - activeFineUpper.T(); - - restrictDiag(fineDiag, coarseDiag); - - agglomerateCoeffs - ( - coeffRestrictAddr, - activeCoarseDiag, - activeCoarseUpper, - activeFineUpper, - activeFineUpperTranspose - ); - } - else if - ( - fineDiag.activeType() == blockCoeffBase::LINEAR - || fineUpper.activeType() == blockCoeffBase::LINEAR - ) - { - linearTypeField& activeCoarseDiag = coarseDiag.asLinear(); - - linearTypeField& activeCoarseUpper = coarseUpper.asLinear(); - const linearTypeField& activeFineUpper = fineUpper.asLinear(); - - restrictDiag(fineDiag, coarseDiag); - - agglomerateCoeffs - ( - coeffRestrictAddr, - activeCoarseDiag, - activeCoarseUpper, - activeFineUpper, - activeFineUpper - ); - } - else if - ( - fineDiag.activeType() == blockCoeffBase::SCALAR - || fineUpper.activeType() == blockCoeffBase::SCALAR - ) - { - scalarTypeField& activeCoarseDiag = coarseDiag.asScalar(); - - scalarTypeField& activeCoarseUpper = coarseUpper.asScalar(); - const scalarTypeField& activeFineUpper = fineUpper.asScalar(); - - restrictDiag(fineDiag, coarseDiag); - - agglomerateCoeffs - ( - coeffRestrictAddr, - activeCoarseDiag, - activeCoarseUpper, - activeFineUpper, - activeFineUpper - ); - } - else - { - FatalErrorIn - ( - "autoPtr >" - "BlockMatrixReorderedClustering::restrictMatrix() const" - ) << "Matrix coeff type morphing error, symmetric matrix" - << abort(FatalError); - } - } - else // asymmetric matrix - { - TypeCoeffField& coarseLower = coarseMatrix.lower(); - const TypeCoeffField& fineLower = matrix_.lower(); - - if - ( - fineDiag.activeType() == blockCoeffBase::SQUARE - || fineUpper.activeType() == blockCoeffBase::SQUARE - ) - { - squareTypeField& activeCoarseDiag = coarseDiag.asSquare(); - - squareTypeField& activeCoarseUpper = coarseUpper.asSquare(); - const squareTypeField& activeFineUpper = fineUpper.asSquare(); - - squareTypeField& activeCoarseLower = coarseLower.asSquare(); - const squareTypeField& activeFineLower = fineLower.asSquare(); - - restrictDiag(fineDiag, coarseDiag); - - agglomerateCoeffs - ( - coeffRestrictAddr, - activeCoarseDiag, - activeCoarseUpper, - activeFineUpper, - activeCoarseLower, - activeFineLower - ); - } - else if - ( - fineDiag.activeType() == blockCoeffBase::LINEAR - || fineUpper.activeType() == blockCoeffBase::LINEAR - ) - { - linearTypeField& activeCoarseDiag = coarseDiag.asLinear(); - - linearTypeField& activeCoarseUpper = coarseUpper.asLinear(); - const linearTypeField& activeFineUpper = fineUpper.asLinear(); - - linearTypeField& activeCoarseLower = coarseLower.asLinear(); - const linearTypeField& activeFineLower = fineLower.asLinear(); - - restrictDiag(fineDiag, coarseDiag); - - agglomerateCoeffs - ( - coeffRestrictAddr, - activeCoarseDiag, - activeCoarseUpper, - activeFineUpper, - activeCoarseLower, - activeFineLower - ); - } - else if - ( - fineDiag.activeType() == blockCoeffBase::SCALAR - || fineUpper.activeType() == blockCoeffBase::SCALAR - ) - { - scalarTypeField& activeCoarseDiag = coarseDiag.asScalar(); - - scalarTypeField& activeCoarseUpper = coarseUpper.asScalar(); - const scalarTypeField& activeFineUpper = fineUpper.asScalar(); - - scalarTypeField& activeCoarseLower = coarseLower.asScalar(); - const scalarTypeField& activeFineLower = fineLower.asScalar(); - - restrictDiag(fineDiag, coarseDiag); - - agglomerateCoeffs - ( - coeffRestrictAddr, - activeCoarseDiag, - activeCoarseUpper, - activeFineUpper, - activeCoarseLower, - activeFineLower - ); - } - else - { - FatalErrorIn - ( - "autoPtr >" - "BlockMatrixReorderedClustering::restrictMatrix() const" - ) << "Matrix coeff type morphing error, asymmetric matrix" - << abort(FatalError); - } - } - - // Create and return BlockAMGLevel - return autoPtr > - ( - new coarseBlockAMGLevel - ( - coarseAddrPtr, - coarseMatrixPtr, - this->dict(), - this->type(), - this->groupSize(), - this->minCoarseEqns() - ) - ); -} - - -template -void Foam::BlockMatrixReorderedClustering::restrictResidual -( - const Field& res, - Field& coarseRes -) const -{ - coarseRes = pTraits::zero; - - forAll (res, i) - { - coarseRes[agglomIndex_[i]] += res[i]; - } -} - - -template -void Foam::BlockMatrixReorderedClustering::prolongateCorrection -( - Field& x, - const Field& coarseX -) const -{ - forAll (x, i) - { - x[i] += coarseX[agglomIndex_[i]]; - } -} - - -// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/BlockMatrixReorderedClustering.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/BlockMatrixReorderedClustering.H deleted file mode 100644 index 4b3058a39..000000000 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/BlockMatrixReorderedClustering.H +++ /dev/null @@ -1,217 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Class - BlockMatrixReorderedClustering - -Description - Block matrix AMG coarsening by Jasak clustering algorithm - -Author - Hrvoje Jasak, Wikki Ltd. All rights reserved. - -SourceFiles - BlockMatrixReorderedClustering.C - -\*---------------------------------------------------------------------------*/ - -#ifndef BlockMatrixReorderedClustering_H -#define BlockMatrixReorderedClustering_H - -#include "BlockMatrixCoarsening.H" -#include "BlockLduMatrix.H" -#include "BlockCoeffNorm.H" -#include "BlockCoeff.H" -#include "tolerancesSwitch.H" -#include "cpuTime.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -/*---------------------------------------------------------------------------*\ - Class BlockMatrixReorderedClustering Declaration -\*---------------------------------------------------------------------------*/ - -template -class BlockMatrixReorderedClustering -: - public BlockMatrixCoarsening -{ - // Private Data - - //- Reference to matrix - const BlockLduMatrix& matrix_; - - //- Min group size - const label minGroupSize_; - - //- Max group size - const label maxGroupSize_; - - //- Reference to a templated norm calculator - autoPtr > normPtr_; - - //- Child array: for each fine equation give a clustering index - labelField agglomIndex_; - - //- Group size - label groupSize_; - - //- Number of solo cells - label nSolo_; - - //- Number of coarse equations - label nCoarseEqns_; - - //- Can a coarse level be constructed? - bool coarsen_; - - cpuTime lTime_; - - - // Private Member Functions - - //- Disallow default bitwise copy construct - BlockMatrixReorderedClustering(const BlockMatrixReorderedClustering&); - - // Disallow default bitwise assignment - void operator=(const BlockMatrixReorderedClustering&); - - - //- Calculate clustering index (child) - void calcClustering(); - - //- Restrict CoeffField. Used for diag coefficient - void restrictDiag - ( - const CoeffField& Coeff, - CoeffField& coarseCoeff - ) const; - - //- Agglomerate coeffs, symmetric matrix - template - void agglomerateCoeffs - ( - const labelList& coeffRestrictAddr, - Field& activeCoarseDiag, - Field& activeCoarseUpper, - const Field& activeFineUpper, - const Field& activeFineUpperTranspose - ) const; - - //- Agglomerate coeffs, assymmetric matrix - template - void agglomerateCoeffs - ( - const labelList& coeffRestrictAddr, - Field& activeCoarseDiag, - Field& activeCoarseUpper, - const Field& activeFineUpper, - Field& activeCoarseLower, - const Field& activeFineLower - ) const; - - //- Restrict CoeffField, decoupled version. Used for diag coefficient - void restrictDiagDecoupled - ( - const CoeffField& Coeff, - CoeffField& coarseCoeff - ) const; - - // Private Static Data - - //- Weighting factor - static const debug::tolerancesSwitch weightFactor_; - - //- Diagonal scaling factor - static const debug::tolerancesSwitch diagFactor_; - - -public: - - //- Runtime type information - TypeName("reorderedCluster"); - - - // Constructors - - //- Construct from matrix and group size - BlockMatrixReorderedClustering - ( - const BlockLduMatrix& matrix, - const dictionary& dict, - const label groupSize, - const label minCoarseEqns - ); - - - //- Destructor - virtual ~BlockMatrixReorderedClustering(); - - - // Member Functions - - //- Can a coarse level be constructed? - virtual bool coarsen() const - { - return coarsen_; - } - - //- Restrict matrix - virtual autoPtr > restrictMatrix() const; - - //- Restrict residual - virtual void restrictResidual - ( - const Field& res, - Field& coarseRes - ) const; - - //- Prolongate correction - virtual void prolongateCorrection - ( - Field& x, - const Field& coarseX - ) const; -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#ifdef NoRepository -# include "BlockMatrixReorderedClustering.C" -#endif - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/blockMatrixReorderedClusterings.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/blockMatrixReorderedClusterings.C deleted file mode 100644 index 55f7ab0f6..000000000 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/blockMatrixReorderedClusterings.C +++ /dev/null @@ -1,43 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -\*---------------------------------------------------------------------------*/ - -#include "blockMatrixReorderedClusterings.H" -#include "blockMatrixCoarsenings.H" -#include "coarseBlockAMGLevel.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -makeBlockMatrixCoarsenings(blockMatrixReorderedClustering); - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/blockMatrixReorderedClusterings.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/blockMatrixReorderedClusterings.H deleted file mode 100644 index 3bdf7054a..000000000 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/blockMatrixReorderedClusterings.H +++ /dev/null @@ -1,65 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Class - BlockMatrixReorderedClustering - -Description - Typedefs for block matrix AMG coarsening by Jasak clustering algorithm. - -Author - Hrvoje Jasak, Wikki Ltd. All rights reserved. - -SourceFiles - blockMatrixReorderedClustering.C - -\*---------------------------------------------------------------------------*/ - -#ifndef blockMatrixReorderedClusterings_H -#define blockMatrixReorderedClusterings_H - -#include "scalarBlockMatrixReorderedClustering.H" -#include "tensorBlockMatrixReorderedClustering.H" -#include "BlockMatrixReorderedClustering.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -typedef BlockMatrixReorderedClustering blockMatrixReorderedClusteringScalar; -typedef BlockMatrixReorderedClustering blockMatrixReorderedClusteringVector; -typedef BlockMatrixReorderedClustering blockMatrixReorderedClusteringTensor; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/scalarBlockMatrixReorderedClustering.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/scalarBlockMatrixReorderedClustering.H deleted file mode 100644 index 84b313a5f..000000000 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/scalarBlockMatrixReorderedClustering.H +++ /dev/null @@ -1,86 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Class - BlockMatrixReorderedClustering - -Description - Specialisation of the BlockMatrixReorderedClustering for scalars. - -Author - Klas Jareteg, 2013-01-31 - -\*---------------------------------------------------------------------------*/ - -#ifndef scalarBlockMatrixReorderedClustering_H -#define scalarBlockMatrixReorderedClustering_H - -#include "blockMatrixCoarsenings.H" -#include "blockMatrixReorderedClusterings.H" -#include "BlockMatrixReorderedClustering.H" -#include "BlockMatrixCoarsening.H" -#include "runTimeSelectionTables.H" -#include "scalarBlockLduMatrix.H" -#include "scalarBlockConstraint.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * Forward declaration of template friend fuctions * * * * * * * // - - -/*---------------------------------------------------------------------------*\ - Class BlockMatrixReorderedClustering Declaration -\*---------------------------------------------------------------------------*/ - -// Disable restrict matrix: no square type -template<> -inline autoPtr > -BlockMatrixReorderedClustering::restrictMatrix() const -{ - FatalErrorIn - ( - "autoPtr > " - "BlockMatrixReorderedClustering::restrictMatrix() const" - ) << "Function not implemented for Type=scalar. " << endl - << abort(FatalError); - - // Dummy return to keep compiler happy - return autoPtr >(NULL); -} - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/tensorBlockMatrixReorderedClustering.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/tensorBlockMatrixReorderedClustering.H deleted file mode 100644 index 53894cb22..000000000 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixReorderedClustering/tensorBlockMatrixReorderedClustering.H +++ /dev/null @@ -1,83 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -Class - BlockMatrixReorderedClustering - -Description - Specialisation of the BlockMatrixReorderedClustering for tensors. - -Author - Klas Jareteg, 2013-01-31 - -\*---------------------------------------------------------------------------*/ - -#ifndef tensorBlockMatrixReorderedClustering_H -#define tensorBlockMatrixReorderedClustering_H - -#include "blockMatrixCoarsenings.H" -#include "blockMatrixReorderedClusterings.H" -#include "BlockMatrixReorderedClustering.H" -#include "runTimeSelectionTables.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -/*---------------------------------------------------------------------------*\ - Class BlockMatrixReorderedClustering Declaration -\*---------------------------------------------------------------------------*/ - -// Disable restrict matrix: no square type -template<> -inline autoPtr > -BlockMatrixReorderedClustering::restrictMatrix -() const -{ - FatalErrorIn - ( - "autoPtr > " - "BlockMatrixReorderedClustering::" - "restrictMatrix() const" - ) << "Function not implemented for Type=tensor. " << endl - << abort(FatalError); - - // Dummy return to keep compiler happy - return autoPtr >(NULL); -} - - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C index e151d4e75..ce3ce69f0 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C @@ -727,7 +727,6 @@ Foam::BlockMatrixSelection::restrictMatrix() const } // Get references to restriction and prolongation matrix - const crMatrix& R = *Rptr_; const crMatrix& P = *Pptr_; @@ -989,8 +988,7 @@ Foam::BlockMatrixSelection::restrictMatrix() const nCoarseEqns_, coarseOwner, coarseNeighbour, - // true - false // DO NOT REUSE STORAGE! HJ, HERE + false // DO NOT REUSE STORAGE! ) ); @@ -1110,12 +1108,77 @@ Foam::BlockMatrixSelection::restrictMatrix() const ); BlockLduMatrix& coarseMatrix = coarseMatrixPtr(); + // Build matrix coefficients + this->updateMatrix(coarseMatrix); + + // Create and return BlockAMGLevel + return autoPtr > + ( + new coarseBlockAMGLevel + ( + coarseAddrPtr, + coarseMatrixPtr, + this->dict(), + this->type(), + 0, // Group size is not used in SAMG + this->minCoarseEqns() + ) + ); +} + + +template +void Foam::BlockMatrixSelection::updateMatrix +( + BlockLduMatrix& coarseMatrix +) const +{ + // Get references to restriction and prolongation matrix + const crMatrix& R = *Rptr_; + const crMatrix& P = *Pptr_; + + // Get connectivity + const crAddressing& crR = R.crAddr(); + const crAddressing& crP = P.crAddr(); + + // Restriction addressing + const labelList& rowR = crR.rowStart(); + const labelList& colR = crR.column(); + + // Matrix A addressing + const unallocLabelList& rowA = matrix_.lduAddr().ownerStartAddr(); + const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); + + // Get interfaces from fine matrix + const typename BlockLduInterfaceFieldPtrsList::Type& + interfaceFields = matrix_.interfaces(); + + // Addressing for lower triangle loop + const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); + const unallocLabelList& losortAddr = matrix_.lduAddr().losortAddr(); + const unallocLabelList& losortStart = matrix_.lduAddr().losortStartAddr(); + + // Prolongation addressing + const labelList& rowP = crP.rowStart(); + const labelList& colP = crP.column(); + + // In order to avoid searching for the off-diagonal coefficient, + // a mark array is used for each row's assembly. + // Mark records the index of the off-diagonal coefficient in the row + // for each neighbour entry. It is reset after completing each row of the + // coarse matrix. + // HJ, 28/Apr/2017 + labelList coeffLabel(nCoarseEqns_, -1); + typedef CoeffField TypeCoeffField; TypeCoeffField& coarseUpper = coarseMatrix.upper(); TypeCoeffField& coarseDiag = coarseMatrix.diag(); TypeCoeffField& coarseLower = coarseMatrix.lower(); + // Get the coarse interfaces and coefficients + lduInterfacePtrsList coarseInterfaces = coarseMatrix.mesh().interfaces(); + //------------------------------------------------------------------------------ // GET COEFFICIENTS //------------------------------------------------------------------------------ @@ -1413,20 +1476,6 @@ Foam::BlockMatrixSelection::restrictMatrix() const ) << "Matrix diagonal of scalar or linear type not implemented" << abort(FatalError); } - - // Create and return BlockAMGLevel - return autoPtr > - ( - new coarseBlockAMGLevel - ( - coarseAddrPtr, - coarseMatrixPtr, - this->dict(), - this->type(), - 0, // Group size is not used in SAMG - this->minCoarseEqns() - ) - ); } diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.H index d0881e308..96b23c2be 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.H @@ -185,6 +185,9 @@ public: Field& x, const Field& coarseX ) const; + + //- Update coarse matrix using same coefficients + virtual void updateMatrix(BlockLduMatrix& coarseMatrix) const; }; diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/scalarBlockMatrixSelection.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/scalarBlockMatrixSelection.H index 368b8bffc..165ea319b 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/scalarBlockMatrixSelection.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/scalarBlockMatrixSelection.H @@ -49,18 +49,17 @@ namespace Foam // Disable restrict matrix // HJ: remove inline in refactoring template<> -inline autoPtr > -BlockMatrixSelection::restrictMatrix() const +inline void BlockMatrixSelection::updateMatrix +( + BlockLduMatrix& coarseMatrix +) const { FatalErrorIn ( "autoPtr > " - "BlockMatrixSelection::restrictMatrix() const" + "BlockMatrixSelection::updateMatrix(...) const" ) << "Function not implemented for Type=scalar. " << endl << abort(FatalError); - - // Dummy return to keep compiler happy - return autoPtr >(NULL); } diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/tensorBlockMatrixSelection.H b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/tensorBlockMatrixSelection.H index 270845f5a..23fc53760 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/tensorBlockMatrixSelection.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/tensorBlockMatrixSelection.H @@ -49,19 +49,18 @@ namespace Foam // Disable restrict matrix: no square type // HJ: remove inline in refactoring template<> -inline autoPtr > -BlockMatrixSelection::restrictMatrix() const +inline void BlockMatrixSelection::updateMatrix +( + BlockLduMatrix& coarseMatrix +) const { FatalErrorIn ( "autoPtr > " "BlockMatrixSelection::" - "restrictMatrix() const" + "updateMatrix(...) const" ) << "Function not implemented for Type=tensor. " << endl << abort(FatalError); - - // Dummy return to keep compiler happy - return autoPtr >(NULL); } diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.C b/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.C index bafb06dc4..694f2a152 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.C +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.C @@ -104,6 +104,13 @@ Foam::coarseBlockAMGLevel::~coarseBlockAMGLevel() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template +Foam::BlockLduMatrix& Foam::coarseBlockAMGLevel::matrix() +{ + return matrixPtr_(); +} + + template Foam::Field& Foam::coarseBlockAMGLevel::x() { @@ -334,4 +341,23 @@ Foam::coarseBlockAMGLevel::makeNextLevel() const } +template +void Foam::coarseBlockAMGLevel::initLevel +( + autoPtr >& coarseLevelPtr +) +{ + // Update matrix by repeating coarsening + + // Update smoother for new matrix + smootherPtr_->initMatrix(); + + // Update coarse matrix if it exists + if (coarseLevelPtr.valid()) + { + coarseningPtr_->updateMatrix(coarseLevelPtr->matrix()); + } +} + + // ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.H b/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.H index fbeb70d75..ade3cc715 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.H @@ -125,6 +125,9 @@ public: return dict_; } + //- Return reference to matrix + virtual BlockLduMatrix& matrix(); + //- Return reference to x virtual Field& x(); @@ -184,6 +187,13 @@ public: //- Create next level from current level virtual autoPtr > makeNextLevel() const; + + //- Re-initialise AMG level after matrix coefficient update + // and update matrix coefficients on coarse level using same coarsening + virtual void initLevel + ( + autoPtr >& coarseLevelPtr + ); }; diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/fineBlockAMGLevel.C b/src/foam/matrices/blockLduMatrix/BlockAMG/fineBlockAMGLevel.C index fc22a49f6..de915e45b 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/fineBlockAMGLevel.C +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/fineBlockAMGLevel.C @@ -81,11 +81,23 @@ Foam::fineBlockAMGLevel::fineBlockAMGLevel // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template +Foam::BlockLduMatrix& Foam::fineBlockAMGLevel::matrix() +{ + FatalErrorIn("Field& Foam::fineBlockAMGLevel::matrix()") + << "Access to matrix is not available on fine level." + << abort(FatalError); + + // Dummy return + return const_cast&>(matrix_); +} + + template Foam::Field& Foam::fineBlockAMGLevel::x() { FatalErrorIn("Field& Foam::fineBlockAMGLevel::x()") - << "x is not available." + << "x is not available on fine level." << abort(FatalError); // Dummy return @@ -97,7 +109,7 @@ template Foam::Field& Foam::fineBlockAMGLevel::b() { FatalErrorIn("Field& Foam::fineBlockAMGLevel::b()") - << "b is not available." + << "b is not available on fine level." << abort(FatalError); // Dummy return @@ -279,4 +291,23 @@ Foam::fineBlockAMGLevel::makeNextLevel() const } +template +void Foam::fineBlockAMGLevel::initLevel +( + autoPtr >& coarseLevelPtr +) +{ + // Fine matrix has been updated externally + + // Update smoother for new matrix + smootherPtr_->initMatrix(); + + // Update coarse matrix if it exists + if (coarseLevelPtr.valid()) + { + coarseningPtr_->updateMatrix(coarseLevelPtr->matrix()); + } +} + + // ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/fineBlockAMGLevel.H b/src/foam/matrices/blockLduMatrix/BlockAMG/fineBlockAMGLevel.H index 2ce5b4638..fa78b44ef 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/fineBlockAMGLevel.H +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/fineBlockAMGLevel.H @@ -125,6 +125,9 @@ public: return dict_; } + //- Return reference to matrix + virtual BlockLduMatrix& matrix(); + //- Return reference to x virtual Field& x(); @@ -183,6 +186,13 @@ public: //- Create next level from current level virtual autoPtr > makeNextLevel() const; + + //- Re-initialise AMG level after matrix coefficient update + // and update matrix coefficients on coarse level using same coarsening + virtual void initLevel + ( + autoPtr >& coarseLevelPtr + ); }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C index d764f1738..9aa243faf 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C @@ -67,7 +67,6 @@ Foam::BlockLduMatrix::BlockLduMatrix(const lduMesh& ldu) coupleUpper_.set(i, new CoeffField(addr.patchAddr(i).size())); coupleLower_.set(i, new CoeffField(addr.patchAddr(i).size())); } - } diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockAMGPrecon/BlockAMGPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockAMGPrecon/BlockAMGPrecon.C index b56392cc1..6457d8a10 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockAMGPrecon/BlockAMGPrecon.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockAMGPrecon/BlockAMGPrecon.C @@ -141,4 +141,11 @@ void Foam::BlockAMGPrecon::precondition } +template +void Foam::BlockAMGPrecon::initMatrix() +{ + amgPtr_->initMatrix(); +} + + // ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockAMGPrecon/BlockAMGPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockAMGPrecon/BlockAMGPrecon.H index 9ab026a3e..f93762632 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockAMGPrecon/BlockAMGPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockAMGPrecon/BlockAMGPrecon.H @@ -25,7 +25,7 @@ Class BlockAMGPrecon Description - AMG preconditioning + Algebraic Multigrid preconditioning for block matrices Author Klas Jareteg, 2012-12-13 @@ -55,7 +55,7 @@ template class BlockAMGCycle; /*---------------------------------------------------------------------------*\ - Class BlockAMGPrecon Declaration + Class BlockAMGPrecon Declaration \*---------------------------------------------------------------------------*/ template @@ -141,6 +141,9 @@ public: Field& x, const Field& b ) const; + + //- Re-initialise preconditioner after matrix coefficient update + virtual void initMatrix(); }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C index ea48e7475..7b6423575 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C @@ -914,7 +914,7 @@ Foam::BlockCholeskyPrecon::BlockCholeskyPrecon BlockLduPrecon(matrix), preconDiag_(matrix.diag()) { - calcPreconDiag(); + this->calcPreconDiag(); } @@ -1301,4 +1301,13 @@ void Foam::BlockCholeskyPrecon::preconditionT } +template +void Foam::BlockCholeskyPrecon::initMatrix() +{ + preconDiag_ = this->matrix_.diag(); + + this->calcPreconDiag(); +} + + // ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.H index 035f5bf75..0e2b91c97 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.H @@ -25,7 +25,7 @@ Class BlockCholeskyPrecon Description - Incomplete Cholesky preconditioning with no fill-in. + Incomplete Cholesky preconditioning with no fill-in for block matrices. Author Hrvoje Jasak, Wikki Ltd. All rights reserved. @@ -215,6 +215,9 @@ public: Field& xT, const Field& bT ) const; + + //- Re-initialise preconditioner after matrix coefficient update + virtual void initMatrix(); }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagonalPrecon/BlockDiagonalPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagonalPrecon/BlockDiagonalPrecon.H index 3071824a0..baf8384de 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagonalPrecon/BlockDiagonalPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagonalPrecon/BlockDiagonalPrecon.H @@ -107,6 +107,10 @@ public: { return precondition(xT, bT); } + + //- Re-initialise preconditioner after matrix coefficient update + virtual void initMatrix() + {} }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C index 3ea12a25a..62c6bbfff 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C @@ -946,4 +946,14 @@ void Foam::BlockGaussSeidelPrecon::preconditionT } +template +void Foam::BlockGaussSeidelPrecon::initMatrix() +{ + invDiag_.clear(); + LUDiag_.clear(); + + this->calcInvDiag(); +} + + // ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H index df85068b3..dd7657998 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H @@ -172,6 +172,9 @@ public: Field& xT, const Field& bT ) const; + + //- Re-initialise preconditioner after matrix coefficient update + virtual void initMatrix(); }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C index bb15a2df4..d5ed14181 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C @@ -840,4 +840,17 @@ void Foam::BlockILUCpPrecon::preconditionT } +template +void Foam::BlockILUCpPrecon::initMatrix() +{ + // Clear extended matrix + extBlockMatrix_.extendedLower().clear(); + extBlockMatrix_.extendedLower().clear(); + + preconDiag_ = this->matrix_.diag(); + + this->calcFactorization(); +} + + // ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H index 1b318a087..508ced666 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H @@ -168,6 +168,9 @@ public: Field& xT, const Field& bT ) const; + + //- Re-initialise preconditioner after matrix coefficient update + virtual void initMatrix(); }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockLduPrecon/BlockLduPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockLduPrecon/BlockLduPrecon.H index 0c42df87b..aca2beb19 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockLduPrecon/BlockLduPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockLduPrecon/BlockLduPrecon.H @@ -140,6 +140,9 @@ public: "(Field& xT, const Field& bT) const" ); } + + //- Re-initialise preconditioner after matrix coefficient update + virtual void initMatrix() = 0; }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockNoPrecon/BlockNoPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockNoPrecon/BlockNoPrecon.H index 0f7646867..df53588aa 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockNoPrecon/BlockNoPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockNoPrecon/BlockNoPrecon.H @@ -109,6 +109,10 @@ public: { precondition(xT, bT); } + + //- Re-initialise preconditioner after matrix coefficient update + virtual void initMatrix() + {} }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockGaussSeidelSmoother/BlockGaussSeidelSmoother.H b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockGaussSeidelSmoother/BlockGaussSeidelSmoother.H index 903c00b1e..ac1f7e22e 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockGaussSeidelSmoother/BlockGaussSeidelSmoother.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockGaussSeidelSmoother/BlockGaussSeidelSmoother.H @@ -111,6 +111,12 @@ public: gs_.precondition(x, b); } } + + //- Re-initialise preconditioner after matrix coefficient update + virtual void initMatrix() + { + gs_.initMatrix(); + } }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUCpSmoother/BlockILUCpSmoother.H b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUCpSmoother/BlockILUCpSmoother.H index 1de496dd1..7721fc5f9 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUCpSmoother/BlockILUCpSmoother.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUCpSmoother/BlockILUCpSmoother.H @@ -127,6 +127,12 @@ public: x += xCorr_; } } + + //- Re-initialise preconditioner after matrix coefficient update + virtual void initMatrix() + { + precon_.initMatrix(); + } }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUSmoother/BlockILUSmoother.H b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUSmoother/BlockILUSmoother.H index 49abb4d4e..9d7cdebd9 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUSmoother/BlockILUSmoother.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUSmoother/BlockILUSmoother.H @@ -133,6 +133,12 @@ public: x += xCorr_; } } + + //- Re-initialise preconditioner after matrix coefficient update + virtual void initMatrix() + { + precon_.initMatrix(); + } }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockLduSmoother/BlockLduSmoother.H b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockLduSmoother/BlockLduSmoother.H index 8fa047bd4..c8c687fe1 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockLduSmoother/BlockLduSmoother.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockLduSmoother/BlockLduSmoother.H @@ -125,6 +125,9 @@ public: const Field& b, const label nSweeps ) const = 0; + + //- Re-initialise preconditioner after matrix coefficient update + virtual void initMatrix() = 0; }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAMGSolver/BlockAMGSolver.H b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAMGSolver/BlockAMGSolver.H index cc20a27c0..a0623705a 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAMGSolver/BlockAMGSolver.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAMGSolver/BlockAMGSolver.H @@ -51,7 +51,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class BlockAMGSolver Declaration + Class BlockAMGSolver Declaration \*---------------------------------------------------------------------------*/ template @@ -93,10 +93,9 @@ public: ); - // Destructor - - virtual ~BlockAMGSolver() - {} + //- Destructor + virtual ~BlockAMGSolver() + {} // Member Functions @@ -107,6 +106,13 @@ public: Field& x, const Field& b ); + + //- Re-initialise solver after matrix coefficient update + virtual void initMatrix() + { + Info<< "REINITIALISE AMG" << endl; + amg_.initMatrix(); + } }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.H b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.H index 1a316b9fc..59cb04519 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.H @@ -102,6 +102,12 @@ public: Field& x, const Field& b ); + + //- Re-initialise solver after matrix coefficient update + virtual void initMatrix() + { + preconPtr_->initMatrix(); + } }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.H b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.H index 3a8ccab84..6f7c73be1 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.H @@ -88,10 +88,9 @@ public: ); - // Destructor - - virtual ~BlockCGSolver() - {} + //- Destructor + virtual ~BlockCGSolver() + {} // Member Functions @@ -102,6 +101,12 @@ public: Field& x, const Field& b ); + + //- Re-initialise solver after matrix coefficient update + virtual void initMatrix() + { + preconPtr_->initMatrix(); + } }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockDiagonal/BlockDiagonalSolver.H b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockDiagonal/BlockDiagonalSolver.H index 69f7f540a..27fc677ca 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockDiagonal/BlockDiagonalSolver.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockDiagonal/BlockDiagonalSolver.H @@ -98,6 +98,10 @@ public: Field& x, const Field& b ); + + //- Re-initialise solver after matrix coefficient update + virtual void initMatrix() + {} }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.H b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.H index f05378da6..158fe9b8f 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.H @@ -115,6 +115,12 @@ public: Field& x, const Field& b ); + + //- Re-initialise solver after matrix coefficient update + virtual void initMatrix() + { + preconPtr_->initMatrix(); + } }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.H b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.H index 02e9f58c2..f99ae7986 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.H @@ -107,6 +107,12 @@ public: Field& x, const Field& b ); + + //- Re-initialise solver after matrix coefficient update + virtual void initMatrix() + { + gs_.initMatrix(); + } }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockILU/BlockILUSolver.H b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockILU/BlockILUSolver.H index 0f5d25b54..22b88b6a8 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockILU/BlockILUSolver.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockILU/BlockILUSolver.H @@ -105,6 +105,13 @@ public: Field& x, const Field& b ); + + + //- Re-initialise solver after matrix coefficient update + virtual void initMatrix() + { + ilu_.initMatrix(); + } }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.H b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.H index b748adb23..5d4cf0490 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.H @@ -208,6 +208,9 @@ public: TypeField& x, const TypeField& b ) = 0; + + //- Re-initialise solver after matrix coefficient update + virtual void initMatrix() = 0; }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/Segregated/SegregatedSolver.H b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/Segregated/SegregatedSolver.H index f7b19298b..59eead92f 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/Segregated/SegregatedSolver.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/Segregated/SegregatedSolver.H @@ -95,10 +95,9 @@ public: ); - // Destructor - - virtual ~SegregatedSolver() - {} + //- Destructor + virtual ~SegregatedSolver() + {} // Member Functions @@ -110,6 +109,9 @@ public: const Field& b ); + //- Re-initialise solver after matrix coefficient update + virtual void initMatrix() + {} }; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C index a58a50e3c..7d472d67d 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C @@ -50,9 +50,7 @@ License #include "blockAMGSolvers.H" #include "blockAMGPrecons.H" #include "blockMatrixCoarsenings.H" -#include "blockMatrixAgglomerations.H" #include "blockMatrixClusterings.H" -#include "blockMatrixReorderedClusterings.H" #include "blockMatrixSelections.H" #include "blockCoeffNorms.H" #include "blockCoeffTwoNorms.H" @@ -153,15 +151,9 @@ typedef BlockMatrixCoarsening block##Type##MatrixCoarsening; \ defineNamedTemplateTypeNameAndDebug(block##Type##MatrixCoarsening, 0); \ defineTemplateRunTimeSelectionTable(block##Type##MatrixCoarsening, matrix); \ \ -typedef BlockMatrixAgglomeration block##Type##MatrixAgglomeration; \ -makeBlockMatrixCoarsening(block##Type##MatrixCoarsening, block##Type##MatrixAgglomeration); \ - \ typedef BlockMatrixClustering block##Type##MatrixClustering; \ makeBlockMatrixCoarsening(block##Type##MatrixCoarsening, block##Type##MatrixClustering); \ \ -typedef BlockMatrixReorderedClustering block##Type##MatrixReorderedClustering; \ -makeBlockMatrixCoarsening(block##Type##MatrixCoarsening, block##Type##MatrixReorderedClustering); \ - \ typedef BlockMatrixSelection block##Type##MatrixSelection; \ makeBlockMatrixCoarsening(block##Type##MatrixCoarsening, block##Type##MatrixSelection); \ \