From ca64bab21d3021ffa6760004a8432e62053f59a4 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Wed, 8 Jul 2015 08:38:47 +0200 Subject: [PATCH] Second version of BlockILUCp preconditioner --- .../extendedBlockLduMatrix.C | 208 +++++++++++++++--- .../extendedBlockLduMatrix.H | 8 +- .../BlockILUCpPrecon/BlockILUCpPrecon.C | 5 +- .../BlockILUCpPrecon/blockILUCpPrecons.C | 4 +- .../BlockLduSolvers/blockVectorNSolvers.C | 4 + 5 files changed, 198 insertions(+), 31 deletions(-) diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C index 2d5a284a2..02893eaf4 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C @@ -30,7 +30,7 @@ License template Foam::extendedBlockLduMatrix::extendedBlockLduMatrix ( - const blockLduMatrix& blockLdum, + const BlockLduMatrix& blockLdum, const label extensionLevel, const polyMesh& polyMesh ) @@ -68,24 +68,74 @@ Foam::extendedBlockLduMatrix::extendedBlockLduMatrix const unallocLabelList& faceMap = extLduAddr_.faceMap(); // Avoid assuming it's upper if the matrix is symmetric - if (blockLdum.hasUpper()) + if (blockLdum.thereIsUpper()) { // Allocate extended upper only extendedUpperPtr_ = new TypeCoeffField ( - extLduAddr_.extendedUpperAddr().size(), - pTraits::zero + extLduAddr_.extendedUpperAddr().size() ); TypeCoeffField& extUpper = *extendedUpperPtr_; // Get upper coeffs from underlying lduMatrix const TypeCoeffField& upper = blockLdum.upper(); - // Copy non-zero coeffs from basic lduMatrix into corresponding - // positions - forAll (upper, faceI) + if (upper.activeType() == blockCoeffBase::SCALAR) { - extUpper[faceMap[faceI]] = 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]; + } + } + 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) + { + activeExtUpper[faceMap[faceI]] = activeUpper[faceI]; + } + } + 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(); + + // Copy non-zero coeffs from basic lduMatrix into corresponding + // positions + forAll (upper, faceI) + { + activeExtUpper[faceMap[faceI]] = activeUpper[faceI]; + } + } + else + { + FatalErrorIn + ( + "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" + ) << "Problem between ordinary block matrix and extended" + << " block matrix upper coeffs type morphing." + << abort(FatalError); } } else @@ -93,19 +143,69 @@ Foam::extendedBlockLduMatrix::extendedBlockLduMatrix // Allocate extended lower only extendedLowerPtr_ = new TypeCoeffField ( - extLduAddr_.extendedLowerAddr().size(), - pTraits::zero + extLduAddr_.extendedLowerAddr().size() ); TypeCoeffField& extLower = *extendedLowerPtr_; // Get lower coeffs from underlying lduMatrix const TypeCoeffField& lower = blockLdum.lower(); - // Copy non-zero coeffs from basic lduMatrix into corresponding - // positions - forAll (lower, faceI) + if (lower.activeType() == blockCoeffBase::SCALAR) { - extLower[faceMap[faceI]] = lower[faceI]; + // 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); } } } @@ -118,22 +218,82 @@ Foam::extendedBlockLduMatrix::extendedBlockLduMatrix const label nExtFaces = extLduAddr_.extendedUpperAddr().size(); // Allocate extended upper and lower - extendedUpperPtr_ = new TypeCoeffField(nExtFaces, pTraits::zero); + extendedUpperPtr_ = new TypeCoeffField(nExtFaces); TypeCoeffField& extUpper = *extendedUpperPtr_; - extendedLowerPtr_ = new TypeCoeffField(nExtFaces, pTraits::zero); + extendedLowerPtr_ = new TypeCoeffField(nExtFaces); TypeCoeffField& extLower = *extendedLowerPtr_; // Get upper and lower coeffs from underlying lduMatrix const TypeCoeffField& upper = blockLdum.upper(); const TypeCoeffField& lower = blockLdum.lower(); - // Copy non-zero coeffs from basic lduMatrix into corresponding - // positions - forAll (upper, faceI) + // Assuming lower and upper have the same active type + if (upper.activeType() == blockCoeffBase::SCALAR) { - extUpper[faceMap[faceI]] = upper[faceI]; - extLower[faceMap[faceI]] = lower[faceI]; + // Helper type definition + typedef typename CoeffField::scalarTypeField activeType; + + // Get references to fields + const activeType& activeUpper = upper.asScalar(); + activeType& activeExtUpper = extUpper.asScalar(); + const activeType& activeLower = lower.asScalar(); + activeType& activeExtLower = extLower.asScalar(); + + // Copy non-zero coeffs from basic lduMatrix into corresponding + // positions + forAll (upper, faceI) + { + activeExtUpper[faceMap[faceI]] = activeUpper[faceI]; + activeExtLower[faceMap[faceI]] = activeLower[faceI]; + } + } + 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(); + const activeType& activeLower = lower.asLinear(); + activeType& activeExtLower = extLower.asLinear(); + + // Copy non-zero coeffs from basic lduMatrix into corresponding + // positions + forAll (upper, faceI) + { + activeExtUpper[faceMap[faceI]] = activeUpper[faceI]; + activeExtLower[faceMap[faceI]] = activeLower[faceI]; + } + } + 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(); + const activeType& activeLower = lower.asSquare(); + activeType& activeExtLower = extLower.asSquare(); + + // Copy non-zero coeffs from basic lduMatrix into corresponding + // positions + forAll (upper, faceI) + { + activeExtUpper[faceMap[faceI]] = activeUpper[faceI]; + activeExtLower[faceMap[faceI]] = activeLower[faceI]; + } + } + else + { + FatalErrorIn + ( + "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" + ) << "Problem between ordinary block matrix and extended" + << " block matrix upper/lower coeffs type morphing." + << abort(FatalError); } } } @@ -172,8 +332,7 @@ Foam::extendedBlockLduMatrix::extendedLower() { extendedLowerPtr_ = new TypeCoeffField ( - extLduAddr_.extendedLowerAddr().size(), - pTraits::zero + extLduAddr_.extendedLowerAddr().size() ); } } @@ -196,8 +355,7 @@ Foam::extendedBlockLduMatrix::extendedUpper() { extendedUpperPtr_ = new TypeCoeffField ( - extLduAddr_.extendedUpperAddr().size(), - pTraits::zero + extLduAddr_.extendedUpperAddr().size() ); } } diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H index 3147d6209..54804f7da 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H @@ -58,7 +58,7 @@ class extendedBlockLduMatrix { // Public data type - typedef typename BlockLduMatrix::TypeCoeffField TypeCoeffField; + typedef CoeffField TypeCoeffField; private: @@ -145,6 +145,12 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#ifdef NoRepository +# include "extendedBlockLduMatrix.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + #endif // ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C index a16f5bd27..d7dd651fd 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C @@ -54,7 +54,7 @@ void Foam::BlockILUCpPrecon::calcFactorization if (!this->matrix_.diagonal()) { // Create multiplication function object - typename BlockCoeff::multiply mult; + typename BlockCoeff::multiply mult; // Get necessary const access to extended ldu addressing const extendedLduAddressing& addr = extBlockMatrix_.extendedLduAddr(); @@ -98,7 +98,6 @@ void Foam::BlockILUCpPrecon::calcFactorization // Start and end of k-th row (upper) and k-th column (lower) fStart = ownStartPtr[rowI]; fEnd = ownStartPtr[rowI + 1]; - // Initialize temporary working diagonal zDiagI = diagPtr[rowI]; @@ -172,7 +171,7 @@ void Foam::BlockILUCpPrecon::calcFactorization // Update diagonal entry, inverting it for future use LDUType& diagRowI = diagPtr[rowI]; - diagRowI = mult.inv(zDiagI); + diagRowI = mult.inverse(zDiagI); // Index for updating L and U register label zwIndex; diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.C index f4dace338..4afe80483 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.C @@ -39,9 +39,9 @@ namespace Foam // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -// Excluded tensor and vector, need reduced versions. VV, 3/Jul/2015. +// Excluded tensor, need reduced version. VV, 3/Jul/2015. makeBlockPrecon(blockScalarPrecon, blockILUCpPreconScalar); -//makeBlockPrecon(blockVectorPrecon, blockILUCpPreconVector); +makeBlockPrecon(blockVectorPrecon, blockILUCpPreconVector); //makeBlockPrecon(blockTensorPrecon, blockILUCpPreconTensor); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C index 38cd3d559..d962948ec 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C @@ -31,6 +31,7 @@ License #include "blockDiagonalPrecons.H" #include "blockGaussSeidelPrecons.H" #include "BlockCholeskyPrecon.H" +#include "BlockILUCpPrecon.H" #include "blockLduSmoothers.H" #include "blockGaussSeidelSmoothers.H" @@ -83,6 +84,9 @@ makeBlockPrecon(block##Type##Precon, block##Type##GaussSeidelPrecon); \ typedef BlockCholeskyPrecon block##Type##CholeskyPrecon; \ makeBlockPrecon(block##Type##Precon, block##Type##CholeskyPrecon); \ \ +typedef BlockILUCpPrecon block##Type##ILUCpPrecon; \ +makeBlockPrecon(block##Type##Precon, block##Type##ILUCpPrecon); \ + \ /* Smoothers */ \ typedef BlockLduSmoother block##Type##Smoother; \ defineNamedTemplateTypeNameAndDebug(block##Type##Smoother, 0); \