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
+
+// ************************************************************************* //