diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C index efca1ceab..7207aba0d 100644 --- a/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C +++ b/src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C @@ -37,7 +37,7 @@ Author #include "BlockSAMGInterfaceField.H" #include "coarseBlockAMGLevel.H" #include "PriorityList.H" -#include "SortableList.H" +#include "labelPair.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -53,13 +53,93 @@ Foam::BlockMatrixSelection::epsilon_ // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +template +Foam::autoPtr +Foam::BlockMatrixSelection::filterProlongation +( + const crMatrix& prolongationMatrix, + const labelList& fineFaceCells +) const +{ + // Get the addressing and coefficients of the prolongation matrix + // on my side + const labelList& pRowStart = prolongationMatrix.crAddr().rowStart(); + const labelList& pColumn = prolongationMatrix.crAddr().column(); + const scalarField& pCoeffs = prolongationMatrix.coeffs(); + const label pNCols = prolongationMatrix.crAddr().nCols(); + + // Count how many weighting factors should be sent across - to avoid + // using dynamic lists + label countWeights = 0; + + forAll (fineFaceCells, i) + { + label end = fineFaceCells[i] + 1; + label start = fineFaceCells[i]; + countWeights += pRowStart[end] - pRowStart[start]; + } + + // Create filtered prolongation matrix arrays + labelList filteredRow(fineFaceCells.size() + 1); + labelList filteredCol(countWeights, 0); + scalarField filteredCoeffs(countWeights, 0); + + // Select the part of the prolongation matrix to send (only for equations + // on processor boundary) + + // Mark the start of the compressed row addressing + // Please note: row start addressing and the owner (faceCells) are linked by + // their address, i.e. the beginning of row faceCells[i] will be at + // rowStart[i], and NOT at rowStart[faceCells[i]]! + filteredRow[0] = 0; + + forAll (fineFaceCells, i) + { + label nCoeffs = filteredRow[i]; + const label start = fineFaceCells[i]; + const label end = fineFaceCells[i] + 1; + + // Copy coefficients from the prolongation matrix on my side to the + // processor boundary (filtered) prolongation matrix + for (label k = pRowStart[start]; k < pRowStart[end]; k++) + { + filteredCoeffs[nCoeffs] = pCoeffs[k]; + filteredCol[nCoeffs] = pColumn[k]; + nCoeffs++; + } + + // Update row start addressing for next row + filteredRow[i + 1] = nCoeffs; + } + + // Create filtered prolongation matrix + autoPtr tFilteredP + ( + new crMatrix + ( + fineFaceCells.size(), + pNCols, + filteredRow, + filteredCol + ) + ); + crMatrix& filteredP = tFilteredP(); + + // Transfer coefficients + filteredP.coeffs().transfer(filteredCoeffs); + + return tFilteredP; +} + + + template void Foam::BlockMatrixSelection::calcCoarsening() { //------------------------------------------------------------------------------ // MATRIX DATA: ADDRESSING, COEFFICIENTS, COEFF NORMS //------------------------------------------------------------------------------ - Info<< "Start equation selection" << endl; + // Get addressing const unallocLabelList& rowStart = matrix_.lduAddr().ownerStartAddr(); const unallocLabelList& losortAddr = matrix_.lduAddr().losortAddr(); @@ -78,13 +158,6 @@ void Foam::BlockMatrixSelection::calcCoarsening() scalarField normDiag(matrixSize); coarseningNormPtr_->normalize(normDiag, matrix_.diag()); - // Note: this needs to be untangled for symmetric matrices. - // If the matrix is symmetric and you ask for lower, it will be manufactured - // for you and will double the memory. Therefore, all loops dealing with - // lower need to be under (if matrix.assymetric()) {...} protection. - // Please get the code to work first and then refactor - // HJ, 18/Feb/2017 - // Calculate norm for upper triangle coeffs (magUpper) scalarField normUpper(row.size()); coarseningNormPtr_->normalize(normUpper, matrix_.upper()); @@ -101,17 +174,17 @@ void Foam::BlockMatrixSelection::calcCoarsening() // CALCULATE INTERPOLATION NORM //------------------------------------------------------------------------------ - // Calculate norm for diagonal coefficients (store into magDiag) - scalarField interpolationNormDiag(matrixSize); - coarseningNormPtr_->normalize(interpolationNormDiag, matrix_.diag()); + // Calculate norm for diagonal coefficients (store into magDiag) + scalarField interpolationNormDiag(matrixSize); + interpolationNormPtr_->normalize(interpolationNormDiag, matrix_.diag()); - // Calculate norm for upper triangle coeffs (magUpper) - scalarField interpolationNormUpper(column.size()); - coarseningNormPtr_->normalize(interpolationNormUpper, matrix_.upper()); + // Calculate norm for upper triangle coeffs (magUpper) + scalarField interpolationNormUpper(column.size()); + interpolationNormPtr_->normalize(interpolationNormUpper, matrix_.upper()); - // Calculate norm for lower triangle coeffs (magLower) - scalarField interpolationNormLower(column.size()); - coarseningNormPtr_->normalize(interpolationNormLower, matrix_.lower()); + // Calculate norm for lower triangle coeffs (magLower) + scalarField interpolationNormLower(column.size()); + interpolationNormPtr_->normalize(interpolationNormLower, matrix_.lower()); //------------------------------------------------------------------------------ // CRITERIA FOR COARSENING: FIND STRONG CONNECTIONS FOR EACH ROW @@ -128,58 +201,68 @@ void Foam::BlockMatrixSelection::calcCoarsening() // Create largest norm. It will be multiplied by epsilon later scalarField epsLargestNorm(matrixSize, 0); + register label j = -1; + + for (register label i = 0; i < matrixSize; i++) + { + scalar signDiag = sign(normDiag[i]); + + for (register label k = rowStart[i]; k < rowStart[i + 1]; k++) + { + // Do row coefficient + scalar magAij = mag(min(signDiag*normUpper[k], 0)); + + // Upper triangle + if (magAij > epsLargestNorm[i]) + { + epsLargestNorm[i] = magAij; + } + + // Do col coefficient + scalar magAji = mag(min(signDiag*normLower[k], 0)); + + // Lower triangle + j = column[k]; + + if (magAji > epsLargestNorm[j]) + { + epsLargestNorm[j] = magAji; + } + } + } + + // Multiply largest norm by epsilon. This is now used below + epsLargestNorm *= epsilon_(); + // Count strong elements in each row - for rowStart addressing // Note: checking magnitude of coeff, which accounts for both strong // positive and negative coefficients. HJ, 28/Feb/2017 labelList strongCoeffCounter(matrixSize, 0); - register label j; - for (register label i = 0; i < matrixSize; i++) { - for (register label k = rowStart[i]; k < rowStart[i+1]; k++) + scalar signDiag = sign(normDiag[i]); + + for (register label k = rowStart[i]; k < rowStart [i + 1]; k++) { + // Do row coefficient + scalar magAij = mag(min(signDiag*normUpper[k], 0)); + // Upper triangle - if (magNormUpper[k] > epsLargestNorm[i]) + if (magAij > epsLargestNorm[i]) { - epsLargestNorm[i] = magNormUpper[k]; + strongCoeffCounter[i]++; } - // Lower triangle - j = column[k]; + // Do col coefficient + scalar magAji = mag(min(signDiag*normLower[k], 0)); - if (magNormLower[k] > epsLargestNorm[j]) - { - epsLargestNorm[j] = magNormLower[k]; - } - } - } - - // Multiply largest norm by epsilon. This is now it is used below - epsLargestNorm *= epsilon_(); - - // Count strong elements in the matrix to create addressing - label count = 0; - - for (register label i = 0; i < matrixSize; i++) - { - for (register label k = rowStart[i]; k < rowStart [i+1]; k++) - { // Column of coefficient k in row i j = column[k]; - // Do the upper triangle (rowwise) - if (magNormUpper[k] > epsLargestNorm[i]) - { - strongCoeffCounter[i]++; - count++; - } - - // Do the lower triangle (columnwise) - if (magNormLower[k] > epsLargestNorm[j]) + if (magAji > epsLargestNorm[j]) { strongCoeffCounter[j]++; - count++; } } } @@ -201,17 +284,16 @@ void Foam::BlockMatrixSelection::calcCoarsening() // Counter for counting the strong elements in a row labelList counter(matrixSize, 0); - // Counter for strong elements in a column - labelList influence(matrixSize, 0); - for (register label i = 0; i < matrixSize; i++) { + scalar signDiag = sign(normDiag[i]); + for (register label k = rowStart[i]; k < rowStart [i+1]; k++) { - j = column[k]; - // Check elements in upper triangle - if (magNormUpper[k] > epsLargestNorm[i]) + scalar magAij = mag(min(signDiag*normUpper[k], 0)); + + if (magAij > epsLargestNorm[i]) { // Store the strong elements into crMatrix, use counter // to count the number of negative strong elements for @@ -221,18 +303,21 @@ void Foam::BlockMatrixSelection::calcCoarsening() interpolationNormUpper[k]; counter[i]++; - influence[j]++; } // Check elements in lower triangle - if (magNormLower[k] > epsLargestNorm[j]) + scalar magAji = mag(min(signDiag*normLower[k], 0)); + + // Column of coefficient k in row i + j = column[k]; + + if (magAji > epsLargestNorm[j]) { strongCol[strongRowStart[j] + counter[j]] = i; strongCoeff[strongRowStart[j] + counter[j]] = interpolationNormLower[k]; counter[j]++; - influence[i]++; } } } @@ -260,15 +345,18 @@ void Foam::BlockMatrixSelection::calcCoarsening() { // Set weights for each equation (weight == number of strong // connections in column!) - equationWeight.set(i, influence[i]); + equationWeight.set(i, strongRowStart[i + 1] - strongRowStart[i]); // Label rows without connections as FINE - if (counter[i] == 0) + if (strongRowStart[i + 1] == strongRowStart[i]) { rowLabel_[i] = FINE; } } + // Start counting coarse equations + nCoarseEqns_ = 0; + while (!equationWeight.empty()) { // removeHead = return index with largest weight and remove @@ -276,7 +364,7 @@ void Foam::BlockMatrixSelection::calcCoarsening() if (rowLabel_[topElement] == UNDECIDED) { - rowLabel_[topElement] = COARSE; + rowLabel_[topElement] = nCoarseEqns_; nCoarseEqns_++; // UPPER TRIANGLE - row-wise @@ -342,41 +430,6 @@ void Foam::BlockMatrixSelection::calcCoarsening() } } -//------------------------------------------------------------------------------ -// CALCULATE COARSE MATRIX ADDRESSING -//------------------------------------------------------------------------------ - - // Calculate the number of fine rows with index smaller than i and - // save it - it will be used for prolongation matrix addresing - labelField fineRowsCount(matrixSize, 0); - - for (label i = 0; i < matrixSize; i++) - { - if (rowLabel_[i] == FINE) - { - if (i == 0) - { - fineRowsCount[i] = 1; - } - else - { - fineRowsCount[i] = fineRowsCount[i - 1] + 1; - } - } - else - { - if (i == 0) - { - continue; - } - else - { - fineRowsCount[i] = fineRowsCount[i - 1]; - } - } - - } - //------------------------------------------------------------------------------ // CALCULATING CONTRIBUTIONS TO THE SCALING FACTOR //------------------------------------------------------------------------------ @@ -393,6 +446,49 @@ void Foam::BlockMatrixSelection::calcCoarsening() // Sum of positive COARSE elements in row scalarField posCoarseElemSum(matrixSize, 0); + // From processor boundary + const typename BlockLduInterfaceFieldPtrsList::Type& interfaceFields = + matrix_.interfaces(); + + forAll (interfaceFields, intI) + { + if (interfaceFields.set(intI)) + { + // Get norm of boundary coefficients + scalarField normBou(matrix_.coupleUpper()[intI].size()); + interpolationNormPtr_->normalize + ( + normBou, + matrix_.coupleUpper()[intI] + ); + + // Get addressing + const labelList& owner = + interfaceFields[intI].coupledInterface().faceCells(); + + forAll (normBou, coeffI) + { + // Get row from faceCells + const label i = owner[coeffI]; + + // Get sign of diagonal coeff + scalar signDiag = sign(interpolationNormDiag[i]); + + // Adjust sign of off-diag coeff + // Note: additional minus because the sign of interface + // coeffs is opposite from the normal matrix off-diagonal + + scalar coeff = -signDiag*normBou[coeffI]; + + // Add negative/positive contribution into corresponding field + // which contributes to prolongation weight factor + positiveElemSum += Foam::max(coeff, 0); + negativeElemSum += Foam::min(coeff,0); + } + } + } + + // Adding positive and negative coeffs in each row (numerator of // scaling factor) for (label i = 0; i < matrixSize; i++) @@ -403,21 +499,15 @@ void Foam::BlockMatrixSelection::calcCoarsening() // Extract diagonal component to compare the sign with the // off-diagonal coeffs - search for negative and positive // connections - scalar signDiag = sign(interpolationNormDiag[i]); + scalar signDiag = sign(interpolationNormDiag[i]); // Upper triangle - use rowStart for (label k = rowStart[i]; k < rowStart[i + 1]; k++) { - // POSITIVE contribution - if (sign(normUpper[k]) == signDiag) - { - positiveElemSum[i] += interpolationNormUpper[k]; - } - // NEGATIVE contribution - else - { - negativeElemSum[i] += interpolationNormUpper[k]; - } + scalar coeff = signDiag*interpolationNormUpper[k]; + + positiveElemSum[i] += Foam::max(coeff, 0); + negativeElemSum[i] += Foam::min(coeff, 0); } // Lower triangle - use losort to go row-wise (Why? This @@ -428,16 +518,10 @@ void Foam::BlockMatrixSelection::calcCoarsening() { label index = losortAddr[k]; - // POSITIVE contribution - if (sign(normLower[index]) == signDiag) - { - positiveElemSum[i] += interpolationNormLower[index]; - } - // NEGATIVE contribution - else - { - negativeElemSum[i] += interpolationNormLower[index]; - } + scalar coeff = signDiag*interpolationNormLower[index]; + + positiveElemSum[i] += Foam::max(coeff, 0); + negativeElemSum[i] += Foam::min(coeff, 0); } } } @@ -460,22 +544,12 @@ void Foam::BlockMatrixSelection::calcCoarsening() { // Get column of strong coeff label j = strongCol[index]; + scalar coeff = signDiag*strongCoeff[index]; - // Sum positive connections - if - ( - sign(strongCoeff[index]) == signDiag - && rowLabel_[j] == COARSE - ) + if (rowLabel_[j] != FINE) { - posCoarseElemSum[i] += strongCoeff[index]; - } - else if - ( - rowLabel_[j] == COARSE - ) - { - negCoarseElemSum[i] += strongCoeff[index]; + posCoarseElemSum[i] += Foam::max(coeff, 0); + negCoarseElemSum[i] += Foam::min(coeff, 0); } } } @@ -509,14 +583,14 @@ void Foam::BlockMatrixSelection::calcCoarsening() if (negCoarseElemSum[i] == 0) { -# ifdef FULLDEBUG + //# ifdef FULLDEBUG // Positive connections are used for interpolation! Diagonal // dominance is not conserved! Info << "Interpolation from positive neighbours!" << nl << "Equation " << i << endl; -# endif + //# endif // Row has no negative COARSE contributions. // Interpolating from positive contributions! @@ -529,7 +603,7 @@ void Foam::BlockMatrixSelection::calcCoarsening() { label j = strongCol[k]; - if (rowLabel_[j] == COARSE) + if (rowLabel_[j] != FINE) { prolongationCoeff.append ( @@ -537,7 +611,7 @@ void Foam::BlockMatrixSelection::calcCoarsening() strongCoeff[k]/(diagonalCoeff) ); - prolongationCol.append(j - fineRowsCount[j]); + prolongationCol.append(rowLabel_[j]); rowCount++; } @@ -558,7 +632,7 @@ void Foam::BlockMatrixSelection::calcCoarsening() if ( - rowLabel_[j] == COARSE && + rowLabel_[j] != FINE && sign(strongCoeff[k]) != sign(diagonalCoeff) ) { @@ -570,7 +644,7 @@ void Foam::BlockMatrixSelection::calcCoarsening() /(diagonalCoeff + positiveElemSum[i]) ); - prolongationCol.append(j - fineRowsCount[j]); + prolongationCol.append(rowLabel_[j]); rowCount++; } } @@ -580,7 +654,7 @@ void Foam::BlockMatrixSelection::calcCoarsening() else // The row is COARSE - injection { prolongationCoeff.append(1); - prolongationCol.append(i - fineRowsCount[i]); + prolongationCol.append(rowLabel_[i]); rowCount++; } @@ -604,31 +678,6 @@ void Foam::BlockMatrixSelection::calcCoarsening() prolongation.coeffs().transfer(prolongationCoeff); - // Ranking for smoothing sweeps - first all coarse equations, - // then all fine - labelField eqnRank_(matrixSize, 0); - label rank = 0; - - // COARSE - for (label i = 0; i < matrixSize; i++) - { - if (rowLabel_[i] == COARSE) - { - eqnRank_[rank] = i; - rank++; - } - } - - // FINE - for (label i = 0; i < matrixSize; i++) - { - if (rowLabel_[i] == FINE) - { - eqnRank_[rank] = i; - rank++; - } - } - // 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 @@ -642,6 +691,20 @@ void Foam::BlockMatrixSelection::calcCoarsening() reduce(coarsen_, andOp()); + if (blockLduMatrix::debug >= 2) + { + Pout<< "Coarse level size: " << nCoarseEqns_; + + if (coarsen_) + { + Pout<< ". Accepted" << endl; + } + else + { + Pout<< ". Rejected" << endl; + } + } + if (coarsen_) { // Coarsening OK, make restriction matrix @@ -654,8 +717,6 @@ void Foam::BlockMatrixSelection::calcCoarsening() // Coarsening did not succeed. Delete Pptr deleteDemandDrivenData(Pptr_); } - - Info<< "End equation selection" << endl; } @@ -706,7 +767,6 @@ template Foam::autoPtr > Foam::BlockMatrixSelection::restrictMatrix() const { - Info<< "Start matrix restriction" << endl; if (!coarsen_) { FatalErrorIn("autoPtr samgPolicy::restrictMatrix() const") @@ -771,11 +831,9 @@ Foam::BlockMatrixSelection::restrictMatrix() const labelList coeffLabel(nCoarseEqns_, -1); label nCoarseCoeffs = 0; - // Array for creating rowStart addressing for upper triangle - labelList nUpperCoeffs(nCoarseEqns_, 0); - // Create coarse addressing - owner and neighbour pairs // HJ: give a better guess of size of coarse off-diagonal + HashTable > coarseOwnNbr; DynamicList