From 22a2e59c3b460c2038ded20d299c9889defe0666 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Thu, 10 Sep 2015 09:33:02 +0200 Subject: [PATCH] Minor extendedLdu matrices consistency rewrite compared to ldu matrices --- .../extendedBlockLduMatrix.C | 177 +++++------------- .../extendedLduMatrix/extendedLduMatrix.C | 55 ++---- src/lduSolvers/lduPrecon/ILUCp/ILUCp.C | 5 +- 3 files changed, 71 insertions(+), 166 deletions(-) diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C index da4a82a11..42981b44b 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C @@ -45,146 +45,73 @@ void Foam::extendedBlockLduMatrix::mapOffDiagCoeffs // Get reference to faceMap in extended addressing const unallocLabelList& faceMap = extLduAddr_.faceMap(); - // Avoid assuming it's upper if the matrix is symmetric - if (blockLdum.thereIsUpper()) + // Matrix is considered symmetric if the upper is allocated and lower + // is not allocated. Allocating extended upper only. + extendedUpperPtr_ = new TypeCoeffField + ( + extLduAddr_.extendedUpperAddr().size() + ); + TypeCoeffField& extUpper = *extendedUpperPtr_; + + // Get upper coeffs from underlying lduMatrix + const TypeCoeffField& upper = blockLdum.upper(); + + if (upper.activeType() == blockCoeffBase::SCALAR) { - // Allocate extended upper only - extendedUpperPtr_ = new TypeCoeffField - ( - extLduAddr_.extendedUpperAddr().size() - ); - TypeCoeffField& extUpper = *extendedUpperPtr_; + // Helper type definition + typedef typename CoeffField::scalarTypeField activeType; - // Get upper coeffs from underlying lduMatrix - const TypeCoeffField& upper = blockLdum.upper(); + // Get references to fields + const activeType& activeUpper = upper.asScalar(); + activeType& activeExtUpper = extUpper.asScalar(); - if (upper.activeType() == blockCoeffBase::SCALAR) + // Copy non-zero coeffs from basic lduMatrix into corresponding + // positions + forAll (upper, faceI) { - // Helper type definition - typedef typename CoeffField::scalarTypeField activeType; - - // Get references to fields - const activeType& activeUpper = upper.asScalar(); - activeType& activeExtUpper = extUpper.asScalar(); - - // Copy non-zero coeffs from basic lduMatrix into corresponding - // positions - forAll (upper, faceI) - { - activeExtUpper[faceMap[faceI]] = activeUpper[faceI]; - } + activeExtUpper[faceMap[faceI]] = activeUpper[faceI]; } - else if (upper.activeType() == blockCoeffBase::LINEAR) + } + else if (upper.activeType() == blockCoeffBase::LINEAR) + { + // Helper type definition + typedef typename CoeffField::linearTypeField activeType; + + // Get references to fields + const activeType& activeUpper = upper.asLinear(); + activeType& activeExtUpper = extUpper.asLinear(); + + // Copy non-zero coeffs from basic lduMatrix into corresponding + // positions + forAll (upper, faceI) { - // Helper type definition - typedef typename CoeffField::linearTypeField activeType; - - // Get references to fields - const activeType& activeUpper = upper.asLinear(); - activeType& activeExtUpper = extUpper.asLinear(); - - // Copy non-zero coeffs from basic lduMatrix into corresponding - // positions - forAll (upper, faceI) - { - activeExtUpper[faceMap[faceI]] = activeUpper[faceI]; - } + activeExtUpper[faceMap[faceI]] = activeUpper[faceI]; } - else if (upper.activeType() == blockCoeffBase::SQUARE) - { - // Helper type definition - typedef typename CoeffField::squareTypeField activeType; + } + else if (upper.activeType() == blockCoeffBase::SQUARE) + { + // Helper type definition + typedef typename CoeffField::squareTypeField activeType; - // Get references to fields - const activeType& activeUpper = upper.asSquare(); - activeType& activeExtUpper = extUpper.asSquare(); + // Get references to fields + const activeType& activeUpper = upper.asSquare(); + activeType& activeExtUpper = extUpper.asSquare(); - // Copy non-zero coeffs from basic lduMatrix into corresponding - // positions - forAll (upper, faceI) - { - activeExtUpper[faceMap[faceI]] = activeUpper[faceI]; - } - } - else + // Copy non-zero coeffs from basic lduMatrix into corresponding + // positions + forAll (upper, faceI) { - FatalErrorIn - ( - "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" - ) << "Problem between ordinary block matrix and extended" - << " block matrix upper coeffs type morphing." - << abort(FatalError); + activeExtUpper[faceMap[faceI]] = activeUpper[faceI]; } } else { - // Allocate extended lower only - extendedLowerPtr_ = new TypeCoeffField + FatalErrorIn ( - extLduAddr_.extendedLowerAddr().size() - ); - TypeCoeffField& extLower = *extendedLowerPtr_; - - // Get lower coeffs from underlying lduMatrix - const TypeCoeffField& lower = blockLdum.lower(); - - if (lower.activeType() == blockCoeffBase::SCALAR) - { - // Helper type definition - typedef typename CoeffField::scalarTypeField activeType; - - // Get references to fields - const activeType& activeLower = lower.asScalar(); - activeType& activeExtLower = extLower.asScalar(); - - // Copy non-zero coeffs from basic lduMatrix into corresponding - // positions - forAll (lower, faceI) - { - activeExtLower[faceMap[faceI]] = activeLower[faceI]; - } - } - else if (lower.activeType() == blockCoeffBase::LINEAR) - { - // Helper type definition - typedef typename CoeffField::linearTypeField activeType; - - // Get references to fields - const activeType& activeLower = lower.asLinear(); - activeType& activeExtLower = extLower.asLinear(); - - // Copy non-zero coeffs from basic lduMatrix into corresponding - // positions - forAll (lower, faceI) - { - activeExtLower[faceMap[faceI]] = activeLower[faceI]; - } - } - else if (lower.activeType() == blockCoeffBase::SQUARE) - { - // Helper type definition - typedef typename CoeffField::squareTypeField activeType; - - // Get references to fields - const activeType& activeLower = lower.asSquare(); - activeType& activeExtLower = extLower.asSquare(); - - // Copy non-zero coeffs from basic lduMatrix into corresponding - // positions - forAll (lower, faceI) - { - activeExtLower[faceMap[faceI]] = activeLower[faceI]; - } - } - else - { - FatalErrorIn - ( - "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" - ) << "Problem between ordinary block matrix and extended" - << " block matrix lower coeffs type morphing." - << abort(FatalError); - } + "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" + ) << "Problem between ordinary block matrix and extended" + << " block matrix upper coeffs type morphing." + << abort(FatalError); } } else diff --git a/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.C b/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.C index bd2f2c085..c2f76e2ab 100644 --- a/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.C +++ b/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.C @@ -74,46 +74,23 @@ Foam::extendedLduMatrix::extendedLduMatrix // Get reference to faceMap in extended addressing const unallocLabelList& faceMap = extLduAddr_.faceMap(); - // Avoid assuming it's upper if the matrix is symmetric - if (ldum.hasUpper()) + // Matrix is considered symmetric if the upper is allocated and lower + // is not allocated. Allocating extended upper only. + extendedUpperPtr_ = new scalarField + ( + extLduAddr_.extendedUpperAddr().size(), + 0.0 + ); + scalarField& extUpper = *extendedUpperPtr_; + + // Get upper coeffs from underlying lduMatrix + const scalarField& upper = ldum.upper(); + + // Copy non-zero coeffs from basic lduMatrix into corresponding + // positions + forAll (upper, faceI) { - // Allocate extended upper only - extendedUpperPtr_ = new scalarField - ( - extLduAddr_.extendedUpperAddr().size(), - 0.0 - ); - scalarField& extUpper = *extendedUpperPtr_; - - // Get upper coeffs from underlying lduMatrix - const scalarField& upper = ldum.upper(); - - // Copy non-zero coeffs from basic lduMatrix into corresponding - // positions - forAll (upper, faceI) - { - extUpper[faceMap[faceI]] = upper[faceI]; - } - } - else - { - // Allocate extended lower only - extendedLowerPtr_ = new scalarField - ( - extLduAddr_.extendedLowerAddr().size(), - 0.0 - ); - scalarField& extLower = *extendedLowerPtr_; - - // Get lower coeffs from underlying lduMatrix - const scalarField& lower = ldum.lower(); - - // Copy non-zero coeffs from basic lduMatrix into corresponding - // positions - forAll (lower, faceI) - { - extLower[faceMap[faceI]] = lower[faceI]; - } + extUpper[faceMap[faceI]] = upper[faceI]; } } else diff --git a/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C b/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C index 1704f1ce0..969cc0896 100644 --- a/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C +++ b/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C @@ -49,8 +49,9 @@ namespace Foam addasymMatrixConstructorToTable addILUCpPreconditionerAsymMatrixConstructorToTable_; - // Add to symmetric constructor table as well until we implement Choleskyp - // preconditioner. VV, 27/Jun/2015. + // Add to symmetric constructor table as well. Cholesky with fill in would + // yield the same sparseness pattern as the original matrix, hence it is not + // implemented. VV, 10/Sep/2015. lduPreconditioner:: addsymMatrixConstructorToTable addILUCpPreconditionerSymMatrixConstructorToTable_;