From 52d2c0c0996d8f5a7a2678006f24864eef95c65d Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Sat, 27 Jun 2015 19:27:34 +0200 Subject: [PATCH] ILUCp preconditioner - ILU with run-time selectable level of fill in --- .../extendedLduAddressing.C | 2 +- .../extendedLduMatrix/extendedLduMatrix.H | 4 +- src/lduSolvers/Make/files | 1 + src/lduSolvers/lduPrecon/ILUC0/ILUC0.C | 42 +- src/lduSolvers/lduPrecon/ILUCp/ILUCp.C | 419 ++++++++++++++++++ src/lduSolvers/lduPrecon/ILUCp/ILUCp.H | 146 ++++++ 6 files changed, 596 insertions(+), 18 deletions(-) create mode 100644 src/lduSolvers/lduPrecon/ILUCp/ILUCp.C create mode 100644 src/lduSolvers/lduPrecon/ILUCp/ILUCp.H diff --git a/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.C b/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.C index fbf984c66..70f179a44 100644 --- a/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.C +++ b/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.C @@ -80,7 +80,7 @@ Foam::extendedLduAddressing::extendedLduAddressing "extendedLduAddressing::extendedLduAddressing" ) << "Extension level 0 not allowed as it is the same as ordinary " - << " lduAddressing." + << "lduAddressing." << abort(FatalError); } } diff --git a/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.H b/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.H index 4ee138282..7357f603e 100644 --- a/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.H +++ b/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.H @@ -101,13 +101,13 @@ public: // Access //- Const access to underlying basic lduMatrix - const lduMatrix& basicLduMatrix() + const lduMatrix& basicLduMatrix() const { return basicLduMatrix_; }; //- Const access to extendedLduAddressing - const extendedLduAddressing& extendedLduAddr() + const extendedLduAddressing& extendedLduAddr() const { return extLduAddr_; }; diff --git a/src/lduSolvers/Make/files b/src/lduSolvers/Make/files index ecc549233..b79c0d51e 100644 --- a/src/lduSolvers/Make/files +++ b/src/lduSolvers/Make/files @@ -5,6 +5,7 @@ lduPrecon = lduPrecon $(lduPrecon)/CholeskyPrecon/CholeskyPrecon.C $(lduPrecon)/ILU0/ILU0.C $(lduPrecon)/ILUC0/ILUC0.C +$(lduPrecon)/ILUCp/ILUCp.C $(lduPrecon)/symGaussSeidelPrecon/symGaussSeidelPrecon.C $(lduPrecon)/amgPrecon/amgPrecon.C diff --git a/src/lduSolvers/lduPrecon/ILUC0/ILUC0.C b/src/lduSolvers/lduPrecon/ILUC0/ILUC0.C index f98c67305..784476083 100644 --- a/src/lduSolvers/lduPrecon/ILUC0/ILUC0.C +++ b/src/lduSolvers/lduPrecon/ILUC0/ILUC0.C @@ -291,16 +291,14 @@ void Foam::ILUC0::precondition const unallocLabelList& lowerAddr = addr.lowerAddr(); const unallocLabelList& losortAddr = addr.losortAddr(); - // Get upper matrix coefficients - const scalarField& upper = matrix_.upper(); - // Solve Lz = b with forward substitution. preconLower_ is chosen to // be unit triangular. z does not need to be stored // Initialize x field x = b; - label losortCoeffI; + register label losortCoeffI; + register label rowI; // Forward substitution loop forAll (preconLower_, coeffI) @@ -319,20 +317,28 @@ void Foam::ILUC0::precondition // Multiply with inverse diagonal x *= preconDiag_; - label rowI; - // Back substitution loop - forAllReverse (upper, coeffI) + forAllReverse (preconUpper_, coeffI) { // Get row index rowI = lowerAddr[coeffI]; // Subtract already updated upper part from the solution - x[rowI] -= upper[coeffI]*x[upperAddr[coeffI]]*preconDiag_[rowI]; + x[rowI] -= + preconUpper_[coeffI]*x[upperAddr[coeffI]]*preconDiag_[rowI]; } } else { + WarningIn + ( + "void ILUC0::precondition" + "(scalarField& x, const scalarField& b, const direction cmpt)" + ) << "Unnecessary use of ILUC0 preconditioner for diagonal matrix. " + << nl + << "Use diagonal preconditioner instead." + << endl; + // Diagonal preconditioning forAll(x, i) { @@ -357,9 +363,6 @@ void Foam::ILUC0::preconditionT const unallocLabelList& lowerAddr = addr.lowerAddr(); const unallocLabelList& losortAddr = addr.losortAddr(); - // Get upper matrix coefficients - const scalarField& upper = matrix_.upper(); - // Solve U^T z = b with forward substitution. preconLower_ is chosen to // be unit triangular - U^T (transpose U) "contains" diagonal entries. z // does not need to be stored. @@ -370,11 +373,11 @@ void Foam::ILUC0::preconditionT x[i] = b[i]*preconDiag_[i]; } - label losortCoeffI; - label rowI; + register label losortCoeffI; + register label rowI; // Forward substitution loop - forAll (upper, coeffI) + forAll (preconUpper_, coeffI) { // Get current losortCoeff to ensure row by row access losortCoeffI = losortAddr[coeffI]; @@ -384,7 +387,7 @@ void Foam::ILUC0::preconditionT // Subtract already updated lower (upper transpose) part from the // solution - x[rowI] -= upper[losortCoeffI]*x[lowerAddr[losortCoeffI]]* + x[rowI] -= preconUpper_[losortCoeffI]*x[lowerAddr[losortCoeffI]]* preconDiag_[rowI]; } @@ -399,6 +402,15 @@ void Foam::ILUC0::preconditionT } else { + WarningIn + ( + "void ILUC0::preconditionT" + "(scalarField& x, const scalarField& b, const direction cmpt)" + ) << "Unnecessary use of ILUC0 preconditioner for diagonal matrix. " + << nl + << "Use diagonal preconditioner instead." + << endl; + // Diagonal preconditioning forAll(x, i) { diff --git a/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C b/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C new file mode 100644 index 000000000..c80ea1d1e --- /dev/null +++ b/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C @@ -0,0 +1,419 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + ILUCp + +Description + ILU preconditioning with arbitrary level of fill in (p), based on Crout + algorithm. + + Reference: Saad, Y.: Iterative Methods for Sparse Linear Systems (2nd + Edition), SIAM, 2003. + +Author + Vuko Vukcevic, FMENA Zagreb. All rights reserved + +\*---------------------------------------------------------------------------*/ + +#include "ILUCp.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(ILUCp, 0); + + lduPreconditioner:: + addasymMatrixConstructorToTable + addILUCpPreconditionerAsymMatrixConstructorToTable_; + + // Add to symmetric constructor table as well until we implement Choleskyp + // preconditioner. VV, 27/Jun/2015. + lduPreconditioner:: + addsymMatrixConstructorToTable + addILUCpPreconditionerSymMatrixConstructorToTable_; +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + + +void Foam::ILUCp::calcFactorization() +{ + if (!matrix_.diagonal()) + { + // Get necessary const access to extended ldu addressing + const extendedLduAddressing& addr = extMatrix_.extendedLduAddr(); + + // Get upper/lower extended addressing + const label* const __restrict__ uPtr = addr.extendedUpperAddr().begin(); + const label* const __restrict__ lPtr = addr.extendedLowerAddr().begin(); + + // Get extended owner start addressing + const label* const __restrict__ ownStartPtr = + addr.extendedOwnerStartAddr().begin(); + + // Get extended losort and losort start addressing + const label* const __restrict__ lsrPtr = + addr.extendedLosortAddr().begin(); + const label* const __restrict__ lsrStartPtr = + addr.extendedLosortStartAddr().begin(); + + // Get access to factored matrix entries + scalar* __restrict__ diagPtr = preconDiag_.begin(); + scalar* __restrict__ upperPtr = extMatrix_.extendedUpper().begin(); + scalar* __restrict__ lowerPtr = extMatrix_.extendedLower().begin(); + + // Get access to working fields + scalar* __restrict__ zPtr = z_.begin(); + scalar* __restrict__ wPtr = w_.begin(); + + // Get number of rows + const label nRows = preconDiag_.size(); + + // Define start and end face ("virtual" face when extended addressing is + // used) of this row/column, and number of non zero off diagonal entries + register label fStart, fEnd, fLsrStart, fLsrEnd; + + // Crout LU factorization + + // Row by row loop (k - loop). + for (register label rowI = 0; rowI < nRows; ++rowI) + { + // Start and end of k-th row (upper) and k-th column (lower) + fStart = ownStartPtr[rowI]; + fEnd = ownStartPtr[rowI + 1]; + + // Initialize temporary working diagonal + zDiag_ = diagPtr[rowI]; + + // Initialize temporary working row field + for (register label faceI = fStart; faceI < fEnd; ++faceI) + { + // Note: z addressed by neighbour of face (column index for + // upper), w addressed by neighbour of face (row index for + // lower) + zPtr[uPtr[faceI]] = upperPtr[faceI]; + wPtr[uPtr[faceI]] = lowerPtr[faceI]; + } + + // Start and end of k-th row (lower) and k-th column (upper) + fLsrStart = lsrStartPtr[rowI]; + fLsrEnd = lsrStartPtr[rowI + 1]; + + // Lower coeff loop (first i - loop) + for + ( + register label faceLsrI = fLsrStart; + faceLsrI < fLsrEnd; + ++faceLsrI + ) + { + // Get losort coefficient for this face + const register label losortCoeff = lsrPtr[faceLsrI]; + + // Get corresponding row index for upper (i label) + const label i = lPtr[losortCoeff]; + + // Update diagonal + zDiag_ -= lowerPtr[losortCoeff]*upperPtr[losortCoeff]; + + // Get end of row for cell i + const register label fEndRowi = ownStartPtr[i + 1]; + + // Upper coeff loop (additional loop to avoid checking the + // existence of certain upper coeffs) + for + ( + // Diagonal is already updated (losortCoeff + 1 = start) + register label faceI = losortCoeff + 1; + faceI < fEndRowi; + ++faceI + ) + { + zPtr[uPtr[faceI]] -= lowerPtr[losortCoeff]*upperPtr[faceI]; + wPtr[uPtr[faceI]] -= upperPtr[losortCoeff]*lowerPtr[faceI]; + } + } + + // Update diagonal entry, inverting it for future use + scalar& diagRowI = diagPtr[rowI]; + diagRowI = 1.0/zDiag_; + + // Index for updating L and U + register label zwIndex; + + // Update upper and lower coeffs + for (register label faceI = fStart; faceI < fEnd; ++faceI) + { + // Get index for current face + zwIndex = uPtr[faceI]; + + // Update L and U decomposition for this row (column) + upperPtr[faceI] = zPtr[zwIndex]; + lowerPtr[faceI] = wPtr[zwIndex]*diagRowI; + } + + // Reset temporary working fields + zDiag_ = 0; + + // Only reset parts of the working fields that have been updated in + // this step (for this row and column) + for + ( + register label faceLsrI = fLsrStart; + faceLsrI < fLsrEnd; + ++faceLsrI + ) + { + // Get losort coefficient for this face + const register label losortCoeff = lsrPtr[faceLsrI]; + + // Get corresponding row index for upper (i label) + const label i = lPtr[losortCoeff]; + + // Get end of row for cell i + const register label fEndRowi = ownStartPtr[i + 1]; + + for + ( + register label faceI = losortCoeff + 1; + faceI < fEndRowi; + ++faceI + ) + { + zPtr[uPtr[faceI]] = 0.0; + wPtr[uPtr[faceI]] = 0.0; + } + } + } + } + else + { + forAll (preconDiag_, i) + { + preconDiag_[i] = 1.0/preconDiag_[i]; + } + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::ILUCp::ILUCp +( + const lduMatrix& matrix, + const FieldField& coupleBouCoeffs, + const FieldField& coupleIntCoeffs, + const lduInterfaceFieldPtrsList& interfaces, + const dictionary& dict +) +: + lduPreconditioner + ( + matrix, + coupleBouCoeffs, + coupleIntCoeffs, + interfaces + ), + preconDiag_(matrix_.diag()), + p_(readLabel(dict.lookup("fillInLevel"))), + extMatrix_ + ( + matrix, + p_, + matrix.mesh().thisDb().parent().lookupObject + < + polyMesh + >(polyMesh::defaultRegion) + ), + zDiag_(0), + z_(preconDiag_.size(), 0), + w_(preconDiag_.size(), 0) +{ + calcFactorization(); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::ILUCp::~ILUCp() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::ILUCp::precondition +( + scalarField& x, + const scalarField& b, + const direction +) const +{ + if (!matrix_.diagonal()) + { + // Get matrix addressing + const extendedLduAddressing& addr = extMatrix_.extendedLduAddr(); + const unallocLabelList& upperAddr = addr.extendedUpperAddr(); + const unallocLabelList& lowerAddr = addr.extendedLowerAddr(); + const unallocLabelList& losortAddr = addr.extendedLosortAddr(); + + // Get upper and lower matrix factors + const scalarField& lower = extMatrix_.extendedLower(); + const scalarField& upper = extMatrix_.extendedUpper(); + + // Solve Lz = b with forward substitution. lower is chosen to be unit + // triangular. z does not need to be stored + + // Initialize x field + x = b; + + register label losortCoeffI; + register label rowI; + + // Forward substitution loop + forAll (lower, coeffI) + { + // Get current losortCoeff to ensure row by row access + losortCoeffI = losortAddr[coeffI]; + + // Subtract already updated lower part from the solution + x[upperAddr[losortCoeffI]] -= + lower[losortCoeffI]*x[lowerAddr[losortCoeffI]]; + } + + // Solve Ux = b with back substitution. U is chosen to be upper + // triangular with diagonal entries corresponding to preconDiag_ + + // Multiply with inverse diagonal + x *= preconDiag_; + + // Back substitution loop + forAllReverse (upper, coeffI) + { + // Get row index + rowI = lowerAddr[coeffI]; + + // Subtract already updated upper part from the solution + x[rowI] -= upper[coeffI]*x[upperAddr[coeffI]]*preconDiag_[rowI]; + } + } + else + { + WarningIn + ( + "void ILUCp::precondition" + "(scalarField& x, const scalarField& b, const direction cmpt)" + ) << "Unnecessary use of ILUCp preconditioner for diagonal matrix. " + << nl + << "Use diagonal preconditioner instead." + << endl; + + // Diagonal preconditioning + forAll(x, i) + { + x[i] = b[i]*preconDiag_[i]; + } + } +} + + +void Foam::ILUCp::preconditionT +( + scalarField& x, + const scalarField& b, + const direction cmpt +) const +{ + if (!matrix_.diagonal()) + { + // Get matrix addressing + const extendedLduAddressing& addr = extMatrix_.extendedLduAddr(); + const unallocLabelList& upperAddr = addr.extendedUpperAddr(); + const unallocLabelList& lowerAddr = addr.extendedLowerAddr(); + const unallocLabelList& losortAddr = addr.extendedLosortAddr(); + + // Get upper and lower matrix factors + const scalarField& lower = extMatrix_.extendedLower(); + const scalarField& upper = extMatrix_.extendedUpper(); + + // Solve U^T z = b with forward substitution. lower is chosen to + // be unit triangular - U^T (transpose U) "contains" diagonal entries. z + // does not need to be stored. + + // Initialize x field + forAll(x, i) + { + x[i] = b[i]*preconDiag_[i]; + } + + register label losortCoeffI; + register label rowI; + + // Forward substitution loop + forAll (upper, coeffI) + { + // Get current losortCoeff to ensure row by row access + losortCoeffI = losortAddr[coeffI]; + + // Get row index + rowI = upperAddr[losortCoeffI]; + + // Subtract already updated lower (upper transpose) part from the + // solution + x[rowI] -= upper[losortCoeffI]*x[lowerAddr[losortCoeffI]]* + preconDiag_[rowI]; + } + + // Solve L^T x = z with back substitution. L^T is unit upper triangular + + // Back substitution loop + forAllReverse (lower, coeffI) + { + // Subtract already updated upper part from the solution + x[lowerAddr[coeffI]] -= lower[coeffI]*x[upperAddr[coeffI]]; + } + } + else + { + WarningIn + ( + "void ILUCp::preconditionT" + "(scalarField& x, const scalarField& b, const direction cmpt)" + ) << "Unnecessary use of ILUCp preconditioner for diagonal matrix. " + << nl + << "Use diagonal preconditioner instead." + << endl; + + // Diagonal preconditioning + forAll(x, i) + { + x[i] = b[i]*preconDiag_[i]; + } + } +} + + +// ************************************************************************* // diff --git a/src/lduSolvers/lduPrecon/ILUCp/ILUCp.H b/src/lduSolvers/lduPrecon/ILUCp/ILUCp.H new file mode 100644 index 000000000..274118b3c --- /dev/null +++ b/src/lduSolvers/lduPrecon/ILUCp/ILUCp.H @@ -0,0 +1,146 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + ILUCp + +Description + ILU preconditioning with arbitrary level of fill in (p), based on Crout + algorithm. + + Reference: Saad, Y.: Iterative Methods for Sparse Linear Systems (2nd + Edition), SIAM, 2003. + +Author + Vuko Vukcevic, FMENA Zagreb. All rights reserved + +SourceFiles + ILUCp.C + +\*---------------------------------------------------------------------------*/ + +#ifndef ILUCp_H +#define ILUCp_H + +#include "lduMatrix.H" +#include "extendedLduMatrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class ILUCp Declaration +\*---------------------------------------------------------------------------*/ + +class ILUCp +: + public lduPreconditioner +{ + // Private Data + + //- Preconditioned diagonal + scalarField preconDiag_; + + //- Fill in level + const label p_; + + //- Extended lduMatrix + extendedLduMatrix extMatrix_; + + //- Temporary working diagonal + scalar zDiag_; + + //- Temporary working row field + scalarField z_; + + //- Temporary Working column field + scalarField w_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + ILUCp(const ILUCp&); + + //- Disallow default bitwise assignment + void operator=(const ILUCp&); + + //- Calculate LU factorization + void calcFactorization(); + + +public: + + //- Runtime type information + TypeName("ILUCp"); + + + // Constructors + + //- Construct from matrix and dictionary + ILUCp + ( + const lduMatrix& matrix, + const FieldField& coupleBouCoeffs, + const FieldField& coupleIntCoeffs, + const lduInterfaceFieldPtrsList& interfaces, + const dictionary& dict + ); + + + // Destructor + + virtual ~ILUCp(); + + + // Member Functions + + //- Execute preconditioning + virtual void precondition + ( + scalarField& x, + const scalarField& b, + const direction cmpt = 0 + ) const; + + //- Execute preconditioning with matrix transpose + virtual void preconditionT + ( + scalarField& x, + const scalarField& b, + const direction cmpt = 0 + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //