From 63506b5c81dbf9fa9f8fdb4253f21361603dcef3 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Tue, 14 Jul 2015 10:03:19 +0200 Subject: [PATCH] Third version of BlockILUCp preconditioner - still not working properly --- src/foam/Make/files | 4 + .../extendedBlockLduMatrices.H | 4 +- .../extendedBlockLduMatrix.C | 64 ++-- .../extendedBlockLduMatrix.H | 4 + .../sphericalTensorExtendedBlockLduMatrix.C | 249 +++++++++++++++ .../sphericalTensorExtendedBlockLduMatrix.H | 63 ++++ .../symmTensorExtendedBlockLduMatrix.C | 247 +++++++++++++++ .../symmTensorExtendedBlockLduMatrix.H | 63 ++++ .../tensorExtendedBlockLduMatrix.C | 243 +++++++++++++++ .../tensorExtendedBlockLduMatrix.H | 63 ++++ .../BlockILUCpPrecon/BlockILUCpPrecon.C | 291 +++++++++--------- .../BlockILUCpPrecon/BlockILUCpPrecon.H | 27 +- .../BlockILUCpPrecon/blockILUCpPrecons.C | 5 +- .../BlockILUCpPrecon/blockILUCpPrecons.H | 2 - .../BlockILUCpPrecon/scalarBlockILUCpPrecon.C | 41 ++- .../BlockILUCpPrecon/scalarBlockILUCpPrecon.H | 17 +- .../BlockILUCpPrecon/tensorBlockILUCpPrecon.C | 65 +++- .../BlockILUCpPrecon/tensorBlockILUCpPrecon.H | 18 +- src/foam/primitives/BlockCoeff/BlockCoeff.H | 31 ++ .../BlockCoeff/DecoupledBlockCoeff.H | 22 ++ 20 files changed, 1285 insertions(+), 238 deletions(-) create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.H create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.H create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.H diff --git a/src/foam/Make/files b/src/foam/Make/files index fa3f384e8..caaa327fe 100644 --- a/src/foam/Make/files +++ b/src/foam/Make/files @@ -663,7 +663,11 @@ matrices/blockLduMatrix/BlockLduMatrix/symmTensorBlockLduMatrix.C matrices/blockLduMatrix/BlockLduMatrix/tensorBlockLduMatrix.C matrices/blockLduMatrix/BlockLduMatrix/blockVectorNMatrices.C matrices/blockLduMatrix/BlockLduMatrix/BlockConstraint/scalarBlockConstraint.C + matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrices.C +matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.C +matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.C +matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.C matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockVectorNMatrices.C BlockLduInterfaceField = matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrices.H b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrices.H index 0a64ae2fc..b647b618b 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrices.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrices.H @@ -35,7 +35,9 @@ Author #ifndef extendedBlockLduMatrices_H #define extendedBlockLduMatrices_H -#include "extendedBlockLduMatrix.H" +#include "sphericalTensorExtendedBlockLduMatrix.H" +#include "symmTensorExtendedBlockLduMatrix.H" +#include "tensorExtendedBlockLduMatrix.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C index 02893eaf4..da4a82a11 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C @@ -25,36 +25,14 @@ License #include "extendedBlockLduMatrix.H" -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // template -Foam::extendedBlockLduMatrix::extendedBlockLduMatrix +void Foam::extendedBlockLduMatrix::mapOffDiagCoeffs ( - const BlockLduMatrix& blockLdum, - const label extensionLevel, - const polyMesh& polyMesh + const BlockLduMatrix& blockLdum ) -: - basicBlockLduMatrix_(blockLdum), - extLduAddr_ - ( - extendedLduAddressing::New - ( - polyMesh, - blockLdum.lduAddr(), - extensionLevel - ) - ), - extendedLowerPtr_(NULL), - extendedUpperPtr_(NULL) { - if (debug) - { - Info<< "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&) :" - "Constructing extendedBlockLduMatrix." - << endl; - } - if (blockLdum.diagonal()) { WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)") @@ -299,6 +277,42 @@ Foam::extendedBlockLduMatrix::extendedBlockLduMatrix } +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::extendedBlockLduMatrix::extendedBlockLduMatrix +( + const BlockLduMatrix& blockLdum, + const label extensionLevel, + const polyMesh& polyMesh +) +: + basicBlockLduMatrix_(blockLdum), + extLduAddr_ + ( + extendedLduAddressing::New + ( + polyMesh, + blockLdum.lduAddr(), + extensionLevel + ) + ), + extendedLowerPtr_(NULL), + extendedUpperPtr_(NULL) +{ + if (debug) + { + Info<< "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&) :" + "Constructing extendedBlockLduMatrix." + << endl; + } + + // Map off diag coeffs from original block matrix to this extended block + // matrix + mapOffDiagCoeffs(blockLdum); +} + + // * * * * * * * * * * * * * * * * Destructor* * * * * * * * * * * * * * * * // template diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H index 54804f7da..ded28183d 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H @@ -86,6 +86,10 @@ private: //- Disallow default bitwise assignement void operator=(const extendedBlockLduMatrix&); + //- Map upper and lower coeffs from ordinary block matrix to extended + // block matrix + void mapOffDiagCoeffs(const BlockLduMatrix&); + public: diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.C new file mode 100644 index 000000000..720f8f9ec --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.C @@ -0,0 +1,249 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#ifndef sphericalTensorExtendedBlockLduMatrix_H +#define sphericalTensorExtendedBlockLduMatrix_H + +#include "extendedBlockLduMatrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template<> +void Foam::extendedBlockLduMatrix::mapOffDiagCoeffs +( + const BlockLduMatrix& blockLdum +) +{ + if (blockLdum.diagonal()) + { + WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)") + << "Attempted to create extended lower/upper coeffs for block " + << "matrix that is diagonal." + << nl << endl; + } + else if (blockLdum.symmetric()) + { + // 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()) + { + // Allocate 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) + { + // 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 + { + FatalErrorIn + ( + "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" + ) << "Problem between ordinary block matrix and extended" + << " block matrix upper coeffs type morphing." + << abort(FatalError); + } + } + else + { + // Allocate extended lower only + extendedLowerPtr_ = new TypeCoeffField + ( + 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 + { + FatalErrorIn + ( + "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" + ) << "Problem between ordinary block matrix and extended" + << " block matrix lower coeffs type morphing." + << abort(FatalError); + } + } + } + else + { + // Get reference to faceMap in extended addressing + const unallocLabelList& faceMap = extLduAddr_.faceMap(); + + // Get number of extended faces + const label nExtFaces = extLduAddr_.extendedUpperAddr().size(); + + // Allocate extended upper and lower + extendedUpperPtr_ = new TypeCoeffField(nExtFaces); + TypeCoeffField& extUpper = *extendedUpperPtr_; + + 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(); + + // Assuming lower and upper have the same active type + if (upper.activeType() == blockCoeffBase::SCALAR) + { + // 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 + { + FatalErrorIn + ( + "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" + ) << "Problem between ordinary block matrix and extended" + << " block matrix upper/lower coeffs type morphing." + << abort(FatalError); + } + } +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.H b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.H new file mode 100644 index 000000000..18e9f55ab --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.H @@ -0,0 +1,63 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +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 + extendedBlockLduMatrix + +Description + Template specialisation for sphericalTensor extended block matrix + +Author + Vuko Vukcevic. All rights reserved. + +SourceFiles + sphericalTensorExtendedBlockLduMatrix.C + +\*---------------------------------------------------------------------------*/ + +#ifndef sphericalTensorExtendedBlockLduMatrix_H +#define sphericalTensorExtendedBlockLduMatrix_H + +#include "extendedBlockLduMatrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template<> +void extendedBlockLduMatrix::mapOffDiagCoeffs +( + const BlockLduMatrix& blockLdum +); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.C new file mode 100644 index 000000000..63715b2d3 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.C @@ -0,0 +1,247 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#ifndef symmTensorExtendedBlockLduMatrix_H +#define symmTensorExtendedBlockLduMatrix_H + +#include "extendedBlockLduMatrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template<> +void Foam::extendedBlockLduMatrix::mapOffDiagCoeffs +( + const BlockLduMatrix& blockLdum +) +{ + if (blockLdum.diagonal()) + { + WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)") + << "Attempted to create extended lower/upper coeffs for block " + << "matrix that is diagonal." + << nl << endl; + } + else if (blockLdum.symmetric()) + { + // 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()) + { + // Allocate 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) + { + // 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 + { + FatalErrorIn + ( + "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" + ) << "Problem between ordinary block matrix and extended" + << " block matrix upper coeffs type morphing." + << abort(FatalError); + } + } + else + { + // Allocate extended lower only + extendedLowerPtr_ = new TypeCoeffField + ( + 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 + { + FatalErrorIn + ( + "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" + ) << "Problem between ordinary block matrix and extended" + << " block matrix lower coeffs type morphing." + << abort(FatalError); + } + } + } + else + { + // Get reference to faceMap in extended addressing + const unallocLabelList& faceMap = extLduAddr_.faceMap(); + + // Get number of extended faces + const label nExtFaces = extLduAddr_.extendedUpperAddr().size(); + + // Allocate extended upper and lower + extendedUpperPtr_ = new TypeCoeffField(nExtFaces); + TypeCoeffField& extUpper = *extendedUpperPtr_; + + 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(); + + // Assuming lower and upper have the same active type + if (upper.activeType() == blockCoeffBase::SCALAR) + { + // 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 + { + FatalErrorIn + ( + "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" + ) << "Problem between ordinary block matrix and extended" + << " block matrix upper/lower coeffs type morphing." + << abort(FatalError); + } + } +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.H b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.H new file mode 100644 index 000000000..617fcb091 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.H @@ -0,0 +1,63 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +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 + extendedBlockLduMatrix + +Description + Template specialisation for symmTensor extended block matrix + +Author + Vuko Vukcevic. All rights reserved. + +SourceFiles + symmTensorExtendedBlockLduMatrix.C + +\*---------------------------------------------------------------------------*/ + +#ifndef symmTensorExtendedBlockLduMatrix_H +#define symmTensorExtendedBlockLduMatrix_H + +#include "extendedBlockLduMatrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template<> +void extendedBlockLduMatrix::mapOffDiagCoeffs +( + const BlockLduMatrix& blockLdum +); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.C new file mode 100644 index 000000000..25efd56ea --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.C @@ -0,0 +1,243 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#ifndef tensorExtendedBlockLduMatrix_H +#define tensorExtendedBlockLduMatrix_H + +#include "extendedBlockLduMatrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template<> +void Foam::extendedBlockLduMatrix::mapOffDiagCoeffs +( + const BlockLduMatrix& blockLdum +) +{ + if (blockLdum.diagonal()) + { + WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)") + << "Attempted to create extended lower/upper coeffs for block " + << "matrix that is diagonal." + << nl << endl; + } + else if (blockLdum.symmetric()) + { + // 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()) + { + // Allocate 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) + { + // 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 + { + FatalErrorIn + ( + "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" + ) << "Problem between ordinary block matrix and extended" + << " block matrix upper coeffs type morphing." + << abort(FatalError); + } + } + else + { + // Allocate extended lower only + extendedLowerPtr_ = new TypeCoeffField + ( + 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 + { + FatalErrorIn + ( + "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" + ) << "Problem between ordinary block matrix and extended" + << " block matrix lower coeffs type morphing." + << abort(FatalError); + } + } + } + else + { + // Get reference to faceMap in extended addressing + const unallocLabelList& faceMap = extLduAddr_.faceMap(); + + // Get number of extended faces + const label nExtFaces = extLduAddr_.extendedUpperAddr().size(); + + // Allocate extended upper and lower + extendedUpperPtr_ = new TypeCoeffField(nExtFaces); + TypeCoeffField& extUpper = *extendedUpperPtr_; + + 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(); + + // Assuming lower and upper have the same active type + if (upper.activeType() == blockCoeffBase::SCALAR) + { + // 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 + { + FatalErrorIn + ( + "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)" + ) << "Problem between ordinary block matrix and extended" + << " block matrix upper/lower coeffs type morphing." + << abort(FatalError); + } + } +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.H b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.H new file mode 100644 index 000000000..3393c7a36 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.H @@ -0,0 +1,63 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +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 + extendedBlockLduMatrix + +Description + Template specialisation for tensor extended block matrix + +Author + Vuko Vukcevic. All rights reserved. + +SourceFiles + tensorExtendedBlockLduMatrix.C + +\*---------------------------------------------------------------------------*/ + +#ifndef tensorExtendedBlockLduMatrix_H +#define tensorExtendedBlockLduMatrix_H + +#include "extendedBlockLduMatrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template<> +void extendedBlockLduMatrix::mapOffDiagCoeffs +( + const BlockLduMatrix& blockLdum +); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C index d7dd651fd..732833923 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C @@ -41,20 +41,25 @@ Author template template -void Foam::BlockILUCpPrecon::calcFactorization +void Foam::BlockILUCpPrecon::calcActiveTypeFactorization ( Field& preconD, Field& extUpper, - Field& extLower, - Field& zDiag, - Field& z, - Field& w -) + Field& extLower +) const { if (!this->matrix_.diagonal()) { + // Get number of rows + const label nRows = preconD.size(); + + // Allocate working fields + Field z(nRows, pTraits::zero); + Field w(nRows, pTraits::zero); + LDUType zDiag(pTraits::zero); + // 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(); @@ -81,10 +86,6 @@ void Foam::BlockILUCpPrecon::calcFactorization // Get access to working fields LDUType* __restrict__ zPtr = z.begin(); LDUType* __restrict__ wPtr = w.begin(); - LDUType& zDiagI = zDiag[0]; - - // Get number of rows - const label nRows = preconD.size(); // Define start and end face ("virtual" face when extended addressing is // used) of this row/column. @@ -98,8 +99,9 @@ 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]; + zDiag = diagPtr[rowI]; // Initialize temporary working row and column fields for (register label faceI = fStart; faceI < fEnd; ++faceI) @@ -132,7 +134,7 @@ void Foam::BlockILUCpPrecon::calcFactorization // Update diagonal // WARNING: Not sure about order of multiplication. // Check it. VV, 3/Jul/2015 - zDiagI -= mult + zDiag -= mult.activeTypeMultiply ( lowerPtr[losortCoeff], upperPtr[losortCoeff] @@ -153,7 +155,7 @@ void Foam::BlockILUCpPrecon::calcFactorization { // WARNING: Not sure about order of multiplication. // Check it. VV, 3/Jul/2015 - zPtr[uPtr[faceI]] -= mult + zPtr[uPtr[faceI]] -= mult.activeTypeMultiply ( lowerPtr[losortCoeff], upperPtr[faceI] @@ -161,17 +163,22 @@ void Foam::BlockILUCpPrecon::calcFactorization // WARNING: Not sure about order of multiplication. // Check it. VV, 3/Jul/2015 - wPtr[uPtr[faceI]] -= mult + wPtr[uPtr[faceI]] -= mult.activeTypeMultiply ( - upperPtr[losortCoeff], - lowerPtr[faceI] + lowerPtr[faceI], + mult.activeTypeMultiply + ( +// diagPtr[lPtr[losortCoeff]], // Note: diag already inverted + mult.inverse(zDiag), // Note: diag already inverted + upperPtr[losortCoeff] + ) ); } } // Update diagonal entry, inverting it for future use LDUType& diagRowI = diagPtr[rowI]; - diagRowI = mult.inverse(zDiagI); + diagRowI = mult.inverse(zDiag); // Index for updating L and U register label zwIndex; @@ -187,15 +194,16 @@ void Foam::BlockILUCpPrecon::calcFactorization // WARNING: Not sure about order of multiplication. // Check it. VV, 3/Jul/2015 - lowerPtr[faceI] = mult - ( - diagRowI, - wPtr[zwIndex] - ); +// lowerPtr[faceI] = mult.activeTypeMultiply +// ( +// diagRowI, +// wPtr[zwIndex] +// ); + lowerPtr[faceI] = wPtr[zwIndex]; } // Reset temporary working fields - zDiagI = pTraits::zero; + zDiag = pTraits::zero; // Only reset parts of the working fields that have been updated in // this step (for this row and column) @@ -234,7 +242,7 @@ void Foam::BlockILUCpPrecon::calcFactorization ( "template \n" "template \n" - "void BlockILUCpPrecon::calcFactorization\n" + "void BlockILUCpPrecon::calcActiveTypeFactorization\n" "(\n" " Field& preconD,\n" " Field& extUpper,\n" @@ -252,6 +260,109 @@ void Foam::BlockILUCpPrecon::calcFactorization } +template +void Foam::BlockILUCpPrecon::calcFactorization() const +{ + // Get upper and lower matrix factors + CoeffField& Lower = extBlockMatrix_.extendedLower(); + CoeffField& Upper = extBlockMatrix_.extendedUpper(); + + // Calculate factorization + // Note: lower, diag and upper must have same type as required by the + // algorithm. This is handled by lowest possible promotion + if (preconDiag_.activeType() == blockCoeffBase::SCALAR) + { + if (Upper.activeType() == blockCoeffBase::SCALAR) + { + calcActiveTypeFactorization + ( + preconDiag_.asScalar(), + Upper.asScalar(), + Lower.asScalar() + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + calcActiveTypeFactorization + ( + preconDiag_.asLinear(), // Promotes to linear + Upper.asLinear(), + Lower.asLinear() + ); + } + else if (Upper.activeType() == blockCoeffBase::SQUARE) + { + calcActiveTypeFactorization + ( + preconDiag_.asSquare(), // Promotes to square + Upper.asSquare(), + Lower.asSquare() + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) + { + if (Upper.activeType() == blockCoeffBase::SCALAR) + { + calcActiveTypeFactorization + ( + preconDiag_.asLinear(), + Upper.asLinear(), // Promotes to linear + Lower.asLinear() // Promotes to linear + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + calcActiveTypeFactorization + ( + preconDiag_.asLinear(), + Upper.asLinear(), + Lower.asLinear() + ); + } + else if (Upper.activeType() == blockCoeffBase::SQUARE) + { + calcActiveTypeFactorization + ( + preconDiag_.asSquare(), // Promotes to square + Upper.asSquare(), + Lower.asSquare() + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::SQUARE) + { + if (Upper.activeType() == blockCoeffBase::SCALAR) + { + calcActiveTypeFactorization + ( + preconDiag_.asSquare(), + Upper.asSquare(), // Promotes to square + Lower.asSquare() // Promotes to square + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + calcActiveTypeFactorization + ( + preconDiag_.asSquare(), + Upper.asSquare(), // Promotes to square + Lower.asSquare() // Promotes to square + ); + } + else if (Upper.activeType() == blockCoeffBase::SQUARE) + { + calcActiveTypeFactorization + ( + preconDiag_.asSquare(), + Upper.asSquare(), + Lower.asSquare() + ); + } + } +} + + template template void Foam::BlockILUCpPrecon::LUSubstitute @@ -420,135 +531,9 @@ Foam::BlockILUCpPrecon::BlockILUCpPrecon ( polyMesh::defaultRegion ) - ), - zDiag_(1), - z_(preconDiag_.size()), - w_(preconDiag_.size()) + ) { - // Get upper and lower matrix factors - CoeffField& Lower = extBlockMatrix_.extendedLower(); - CoeffField& Upper = extBlockMatrix_.extendedUpper(); - - // Calculate factorization - // Note: lower, diag and upper must have same type as required by the - // algorithm. This is handled by lowest possible promotion - if (preconDiag_.activeType() == blockCoeffBase::SCALAR) - { - if (Upper.activeType() == blockCoeffBase::SCALAR) - { - calcFactorization - ( - preconDiag_.asScalar(), - Upper.asScalar(), - Lower.asScalar(), - zDiag_.asScalar(), - z_.asScalar(), - w_.asScalar() - ); - } - else if (Upper.activeType() == blockCoeffBase::LINEAR) - { - calcFactorization - ( - preconDiag_.asLinear(), // Promotes to linear - Upper.asLinear(), - Lower.asLinear(), - zDiag_.asLinear(), - z_.asLinear(), - w_.asLinear() - ); - } - else if (Upper.activeType() == blockCoeffBase::SQUARE) - { - calcFactorization - ( - preconDiag_.asSquare(), // Promotes to square - Upper.asSquare(), - Lower.asSquare(), - zDiag_.asSquare(), - z_.asSquare(), - w_.asSquare() - ); - } - } - else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) - { - if (Upper.activeType() == blockCoeffBase::SCALAR) - { - calcFactorization - ( - preconDiag_.asLinear(), - Upper.asLinear(), // Promotes to linear - Lower.asLinear(), // Promotes to linear - zDiag_.asLinear(), // Promotes to linear - z_.asLinear(), // Promotes to linear - w_.asLinear() // Promotes to linear - ); - } - else if (Upper.activeType() == blockCoeffBase::LINEAR) - { - calcFactorization - ( - preconDiag_.asLinear(), - Upper.asLinear(), - Lower.asLinear(), - zDiag_.asLinear(), - z_.asLinear(), - w_.asLinear() - ); - } - else if (Upper.activeType() == blockCoeffBase::SQUARE) - { - calcFactorization - ( - preconDiag_.asSquare(), // Promotes to square - Upper.asSquare(), - Lower.asSquare(), - zDiag_.asSquare(), - z_.asSquare(), - w_.asSquare() - ); - } - } - else if (preconDiag_.activeType() == blockCoeffBase::SQUARE) - { - if (Upper.activeType() == blockCoeffBase::SCALAR) - { - calcFactorization - ( - preconDiag_.asSquare(), - Upper.asSquare(), // Promotes to square - Lower.asSquare(), // Promotes to square - zDiag_.asSquare(), // Promotes to square - z_.asSquare(), // Promotes to square - w_.asSquare() // Promotes to square - ); - } - else if (Upper.activeType() == blockCoeffBase::LINEAR) - { - calcFactorization - ( - preconDiag_.asSquare(), - Upper.asSquare(), // Promotes to square - Lower.asSquare(), // Promotes to square - zDiag_.asSquare(), // Promotes to square - z_.asSquare(), // Promotes to square - w_.asSquare() // Promotes to square - ); - } - else if (Upper.activeType() == blockCoeffBase::SQUARE) - { - calcFactorization - ( - preconDiag_.asSquare(), - Upper.asSquare(), - Lower.asSquare(), - zDiag_.asSquare(), - z_.asSquare(), - w_.asSquare() - ); - } - } + calcFactorization(); } diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H index c59dee611..f56a371f9 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H @@ -63,22 +63,13 @@ class BlockILUCpPrecon // Private Data //- Preconditioned diagonal - CoeffField preconDiag_; + mutable CoeffField preconDiag_; //- Fill in level const label p_; //- Extended lduMarix - extendedBlockLduMatrix extBlockMatrix_; - - //- Temporary working diagonal - CoeffField zDiag_; - - //- Temporary working row field - CoeffField z_; - - //- Temporary working column field - CoeffField w_; + mutable extendedBlockLduMatrix extBlockMatrix_; // Private Member Functions @@ -89,18 +80,18 @@ class BlockILUCpPrecon //- Disallow default bitwise assignment void operator=(const BlockILUCpPrecon&); - //- Calculate LU factorization. + //- Calculate active type dependant factorization. // Note: both lower, diag and upper have to be of same type template - void calcFactorization + void calcActiveTypeFactorization ( Field& preconD, Field& extUpper, - Field& extLower, - Field& zDiag, - Field& z, - Field& w - ); + Field& extLower + ) const; + + //- Calculate LU factorization (constructor helper) + void calcFactorization() const; //- Performs forward/backward LU substitution // Note: both lower, diag and upper have to be of same type diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.C index 4afe80483..83e302df3 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.C @@ -39,10 +39,7 @@ namespace Foam // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -// Excluded tensor, need reduced version. VV, 3/Jul/2015. -makeBlockPrecon(blockScalarPrecon, blockILUCpPreconScalar); -makeBlockPrecon(blockVectorPrecon, blockILUCpPreconVector); -//makeBlockPrecon(blockTensorPrecon, blockILUCpPreconTensor); +makeBlockPrecons(blockILUCpPrecon); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.H index 6e0716277..c942c0118 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.H @@ -42,8 +42,6 @@ SourceFiles #include "scalarBlockILUCpPrecon.H" #include "tensorBlockILUCpPrecon.H" -#include "BlockILUCpPrecon.H" - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.C index ad1a9d923..4a90575ec 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.C @@ -45,18 +45,23 @@ namespace Foam template<> template<> -void BlockILUCpPrecon::calcFactorization +void BlockILUCpPrecon::calcActiveTypeFactorization ( scalarField& preconD, scalarField& extUpper, - scalarField& extLower, - scalarField& zDiag, - scalarField& z, - scalarField& w -) + scalarField& extLower +) const { if (!matrix_.diagonal()) { + // Get number of rows + const label nRows = preconD.size(); + + // Allocate working fields + scalarField z(nRows, 0.0); + scalarField w(nRows, 0.0); + scalar zDiag = 0.0; + // Get necessary const access to extended ldu addressing const extendedLduAddressing& addr = extBlockMatrix_.extendedLduAddr(); @@ -82,10 +87,6 @@ void BlockILUCpPrecon::calcFactorization // Get access to working fields scalar* __restrict__ zPtr = z.begin(); scalar* __restrict__ wPtr = w.begin(); - scalar& zDiagI = zDiag[0]; - - // Get number of rows - const label nRows = preconD.size(); // Define start and end face ("virtual" face when extended addressing is // used) of this row/column. @@ -101,7 +102,7 @@ void BlockILUCpPrecon::calcFactorization fEnd = ownStartPtr[rowI + 1]; // Initialize temporary working diagonal - zDiagI = diagPtr[rowI]; + zDiag = diagPtr[rowI]; // Initialize temporary working row field for (register label faceI = fStart; faceI < fEnd; ++faceI) @@ -132,7 +133,7 @@ void BlockILUCpPrecon::calcFactorization const label i = lPtr[losortCoeff]; // Update diagonal - zDiagI -= lowerPtr[losortCoeff]*upperPtr[losortCoeff]; + zDiag -= lowerPtr[losortCoeff]*upperPtr[losortCoeff]; // Get end of row for cell i const register label fEndRowi = ownStartPtr[i + 1]; @@ -154,7 +155,7 @@ void BlockILUCpPrecon::calcFactorization // Update diagonal entry, inverting it for future use scalar& diagRowI = diagPtr[rowI]; - diagRowI = 1.0/zDiagI; + diagRowI = 1.0/zDiag; // Index for updating L and U register label zwIndex; @@ -171,7 +172,7 @@ void BlockILUCpPrecon::calcFactorization } // Reset temporary working fields - zDiagI = 0; + zDiag = 0; // Only reset parts of the working fields that have been updated in // this step (for this row and column) @@ -228,6 +229,18 @@ void BlockILUCpPrecon::calcFactorization } +template<> +void BlockILUCpPrecon::calcFactorization() const +{ + calcActiveTypeFactorization + ( + preconDiag_.asScalar(), + extBlockMatrix_.extendedUpper().asScalar(), + extBlockMatrix_.extendedLower().asScalar() + ); +} + + template<> void BlockILUCpPrecon::precondition ( diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.H index 47d93814a..3469d9cb0 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.H @@ -45,18 +45,21 @@ SourceFiles namespace Foam { -// Calculate factorization +// Calculate active type factorization template<> template<> -void BlockILUCpPrecon::calcFactorization +void BlockILUCpPrecon::calcActiveTypeFactorization ( scalarField& preconD, scalarField& extUpper, - scalarField& extLower, - scalarField& zDiag, - scalarField& z, - scalarField& w -); + scalarField& extLower +) const; + + +// Calculate factorization (constructor helper) +template<> +void BlockILUCpPrecon::calcFactorization() const; + // Precondition template<> diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/tensorBlockILUCpPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/tensorBlockILUCpPrecon.C index ba88b4d80..9f62308f2 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/tensorBlockILUCpPrecon.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/tensorBlockILUCpPrecon.C @@ -43,24 +43,75 @@ Author namespace Foam { -// Calculate factorization template<> template<> -void BlockILUCpPrecon::calcFactorization +void BlockILUCpPrecon::calcActiveTypeFactorization ( tensorField& preconD, tensorField& extUpper, - tensorField& extLower, - tensorField& zDiag, - tensorField& z, - tensorField& w -) + tensorField& extLower +) const { // Decoupled version notImplemented("void Foam::BlockILUCpPrecon::calcFactorization"); } +template<> +void BlockILUCpPrecon::calcFactorization() const +{ + // Get upper and lower matrix factors + CoeffField& Lower = extBlockMatrix_.extendedLower(); + CoeffField& Upper = extBlockMatrix_.extendedUpper(); + + // Calculate factorization + // Note: lower, diag and upper must have same type as required by the + // algorithm. This is handled by lowest possible promotion + if (preconDiag_.activeType() == blockCoeffBase::SCALAR) + { + if (Upper.activeType() == blockCoeffBase::SCALAR) + { + calcActiveTypeFactorization + ( + preconDiag_.asScalar(), + Upper.asScalar(), + Lower.asScalar() + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + calcActiveTypeFactorization + ( + preconDiag_.asLinear(), // Promotes to linear + Upper.asLinear(), + Lower.asLinear() + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) + { + if (Upper.activeType() == blockCoeffBase::SCALAR) + { + calcActiveTypeFactorization + ( + preconDiag_.asLinear(), + Upper.asLinear(), // Promotes to linear + Lower.asLinear() // Promotes to linear + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + calcActiveTypeFactorization + ( + preconDiag_.asLinear(), + Upper.asLinear(), + Lower.asLinear() + ); + } + } +} + + template<> void BlockILUCpPrecon::precondition ( diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/tensorBlockILUCpPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/tensorBlockILUCpPrecon.H index 60ac4e0fc..7cd0e9cc1 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/tensorBlockILUCpPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/tensorBlockILUCpPrecon.H @@ -45,20 +45,23 @@ SourceFiles namespace Foam { -// Calculate factorization +// Calculate active type factorization template<> template<> -void BlockILUCpPrecon::calcFactorization +void BlockILUCpPrecon::calcActiveTypeFactorization ( tensorField& preconD, tensorField& extUpper, - tensorField& extLower, - tensorField& zDiag, - tensorField& z, - tensorField& w -); + tensorField& extLower +) const; +// Calculate factorization (constructor helper) +template<> +void BlockILUCpPrecon::calcFactorization() const; + + +// Precondition template<> void BlockILUCpPrecon::precondition ( @@ -67,6 +70,7 @@ void BlockILUCpPrecon::precondition ) const; +// Precondition transpose template<> void BlockILUCpPrecon::preconditionT ( diff --git a/src/foam/primitives/BlockCoeff/BlockCoeff.H b/src/foam/primitives/BlockCoeff/BlockCoeff.H index d2aa1df59..c0207adaf 100644 --- a/src/foam/primitives/BlockCoeff/BlockCoeff.H +++ b/src/foam/primitives/BlockCoeff/BlockCoeff.H @@ -126,6 +126,37 @@ public: } + // Coefficient times coefficient multiplication. Needed for BlockILUCp + // preconditioner. VV, 12/Jul/2015. + + scalarType activeTypeMultiply + ( + const scalarType& a, + const scalarType& b + ) const + { + return a*b; + } + + linearType activeTypeMultiply + ( + const linearType& a, + const linearType& b + ) const + { + return cmptMultiply(a, b); + } + + squareType activeTypeMultiply + ( + const squareType& a, + const squareType& b + ) const + { + return (a & b); + } + + // Transpose functions scalarType transpose(const scalarType& c) const diff --git a/src/foam/primitives/BlockCoeff/DecoupledBlockCoeff.H b/src/foam/primitives/BlockCoeff/DecoupledBlockCoeff.H index ea9614175..46e4c7f09 100644 --- a/src/foam/primitives/BlockCoeff/DecoupledBlockCoeff.H +++ b/src/foam/primitives/BlockCoeff/DecoupledBlockCoeff.H @@ -113,6 +113,28 @@ public: } + // Coefficient times coefficient multiplication. Needed for BlockILUCp + // preconditioner. VV, 12/Jul/2015. + + scalarType activeTypeMultiply + ( + const scalarType& a, + const scalarType& b + ) const + { + return a*b; + } + + linearType activeTypeMultiply + ( + const linearType& a, + const linearType& b + ) const + { + return cmptMultiply(a, b); + } + + // Inverse functions scalarType inverse(const scalarType& c) const