From 6c5d2548028e4651e22a0ad56501768206a3eb87 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 16 May 2018 11:50:18 +0100 Subject: [PATCH] Block ILUC0 precon and smoother --- .../BlockILUC0Precon/BlockILUC0Precon.C | 842 ++++++++++++++++++ .../BlockILUC0Precon/BlockILUC0Precon.H | 194 ++++ .../BlockILUC0PreconDecoupled.C | 240 +++++ .../BlockILUC0Precon/blockILUC0Precons.C | 48 + .../BlockILUC0Precon/blockILUC0Precons.H | 64 ++ .../BlockILUC0Precon/scalarBlockILUC0Precon.C | 276 ++++++ .../BlockILUC0Precon/scalarBlockILUC0Precon.H | 90 ++ .../BlockILUC0Precon/tensorBlockILUC0Precon.C | 147 +++ .../BlockILUC0Precon/tensorBlockILUC0Precon.H | 90 ++ .../BlockILUC0Smoother/BlockILUC0Smoother.H | 147 +++ .../BlockILUC0Smoother/blockILUC0Smoothers.C | 42 + .../BlockILUC0Smoother/blockILUC0Smoothers.H | 62 ++ .../BlockLduSolvers/blockVectorNSolvers.C | 8 + 13 files changed, 2250 insertions(+) create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0Precon.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0Precon.H create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0PreconDecoupled.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/blockILUC0Precons.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/blockILUC0Precons.H create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/scalarBlockILUC0Precon.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/scalarBlockILUC0Precon.H create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/tensorBlockILUC0Precon.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/tensorBlockILUC0Precon.H create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUC0Smoother/BlockILUC0Smoother.H create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUC0Smoother/blockILUC0Smoothers.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUC0Smoother/blockILUC0Smoothers.H diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0Precon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0Precon.C new file mode 100644 index 000000000..22f631199 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0Precon.C @@ -0,0 +1,842 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +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 . + +Description + Block variant of 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 "BlockILUC0Precon.H" +#include "error.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +template +void Foam::BlockILUC0Precon::calcActiveTypeFactorization +( + Field& preconD, + Field& preconUpper, + Field& preconLower +) 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; + + // Get necessary const access to matrix addressing + const lduAddressing& addr = this->matrix_.lduAddr(); + + // Get upper/lower addressing + const label* const __restrict__ uPtr = addr.upperAddr().begin(); + const label* const __restrict__ lPtr = addr.lowerAddr().begin(); + + // Get owner start addressing + const label* const __restrict__ ownStartPtr = + addr.ownerStartAddr().begin(); + + // Get losort and losort start addressing + const label* const __restrict__ lsrPtr = addr.losortAddr().begin(); + const label* const __restrict__ lsrStartPtr = + addr.losortStartAddr().begin(); + + // Get access to factored matrix entries + LDUType* __restrict__ diagPtr = preconD.begin(); + LDUType* __restrict__ upperPtr = preconUpper.begin(); + LDUType* __restrict__ lowerPtr = preconLower.begin(); + + // Get access to working fields + LDUType* __restrict__ zPtr = z.begin(); + LDUType* __restrict__ wPtr = w.begin(); + + // Define start and end face ("virtual" face when extended addressing + // is used) of this row/column. + 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 and column fields + for (register label faceI = fStart; faceI < fEnd; ++faceI) + { + // Note: z addressed by neighbour of face (column index for + // upper) + 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/upper coeff loop (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 -= mult.activeTypeMultiply + ( + 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]] -= mult.activeTypeMultiply + ( + lowerPtr[losortCoeff], + upperPtr[faceI] + ); + + wPtr[uPtr[faceI]] -= mult.activeTypeMultiply + ( + lowerPtr[faceI], + upperPtr[losortCoeff] + ); + } + } + + // Update diagonal entry, inverting it for future use + LDUType& diagRowI = diagPtr[rowI]; + diagRowI = mult.inverse(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] = mult.activeTypeMultiply + ( + wPtr[zwIndex], + diagRowI + ); + } + + // Reset temporary working fields + zDiag = pTraits::zero; + + // 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]] = pTraits::zero; + wPtr[uPtr[faceI]] = pTraits::zero; + } + } + } + } + else + { + FatalErrorIn + ( + "template \n" + "template \n" + "void BlockILUC0Precon::calcActiveTypeFactorization\n" + "(\n" + " Field& preconD,\n" + " Field& preconUpper,\n" + " Field& preconLower,\n" + ") const" + ) << "Unnecessary use of BlockILUC0 preconditioner for diagonal " + << "matrix." + << nl + << "Use BlockDiagonal preconditioner instead." + << abort(FatalError); + } +} + + +template +void Foam::BlockILUC0Precon::calcFactorization() const +{ + // Get upper and lower matrix factors + CoeffField& Lower = preconLower_; + CoeffField& Upper = preconUpper_; + + // 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::BlockILUC0Precon::LUSubstitute +( + Field& x, + const Field& preconD, + const Field& upper, + const Field& lower, + const Field& b +) const +{ + // Create multiplication function object + typename BlockCoeff::multiply mult; + + // Get matrix addressing + const lduAddressing& addr = this->matrix_.lduAddr(); + const unallocLabelList& upperAddr = addr.upperAddr(); + const unallocLabelList& lowerAddr = addr.lowerAddr(); + const unallocLabelList& losortAddr = addr.losortAddr(); + + // Solve Lz = b with forward substitution in block form. 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]] -= mult + ( + lower[losortCoeffI], + x[lowerAddr[losortCoeffI]] + ); + } + + // Solve Ux = b with back substitution in block form. U is chosen to be + // upper triangular with diagonal entries corresponding to preconD + + // Multiply with inverse diagonal + forAll(x, i) + { + x[i] = mult(preconD[i], x[i]); + } + + // Back substitution loop + forAllReverse (upper, coeffI) + { + // Get row index + rowI = lowerAddr[coeffI]; + + // Subtract already updated upper part from the solution + x[rowI] -= mult + ( + preconD[rowI], + mult + ( + upper[coeffI], + x[upperAddr[coeffI]] + ) + ); + } +} + + +template +template +void Foam::BlockILUC0Precon::LUSubstituteT +( + Field& xT, + const Field& preconD, + const Field& upper, + const Field& lower, + const Field& bT +) const +{ + // Create multiplication function object + typename BlockCoeff::multiply mult; + + // Get matrix addressing + const lduAddressing& addr = this->matrix_.lduAddr(); + const unallocLabelList& upperAddr = addr.upperAddr(); + const unallocLabelList& lowerAddr = addr.lowerAddr(); + const unallocLabelList& losortAddr = addr.losortAddr(); + + // Solve U^T z = b with forward substitution in block form. lower is + // chosen to be unit triangular - U^T (transpose U) "contains" diagonal + // entries. z does not need to be stored + + // Note: transpose should be used for all block coeffs. + + // Initialize x field + forAll(xT, i) + { + xT[i] = mult + ( + mult.transpose(preconD[i]), + bT[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 + xT[rowI] -= mult + ( + mult.transpose(preconD[rowI]), + mult + ( + mult.transpose(upper[losortCoeffI]), + xT[lowerAddr[losortCoeffI]] + ) + ); + } + + // 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 + xT[lowerAddr[coeffI]] -= mult + ( + mult.transpose(lower[coeffI]), + xT[upperAddr[coeffI]] + ); + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::BlockILUC0Precon::BlockILUC0Precon +( + const BlockLduMatrix& matrix, + const dictionary& dict +) +: + BlockLduPrecon(matrix), + preconDiag_(matrix.diag()), + preconLower_(matrix.lower()), + preconUpper_(matrix.upper()) +{ + this->calcFactorization(); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::BlockILUC0Precon::~BlockILUC0Precon() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::BlockILUC0Precon::precondition +( + Field& x, + const Field& b +) const +{ + if (!this->matrix_.diagonal()) + { + // Get upper and lower matrix factors + const CoeffField& Lower = preconLower_; + const CoeffField& Upper = preconUpper_; + + // Execute preconditioning by LU substitution. + // 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) + { + LUSubstitute + ( + x, + preconDiag_.asScalar(), + Upper.asScalar(), + Lower.asScalar(), + b + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + LUSubstitute + ( + x, + preconDiag_.asLinear(), // Promotes to linear + Upper.asLinear(), + Lower.asLinear(), + b + ); + } + else if (Upper.activeType() == blockCoeffBase::SQUARE) + { + LUSubstitute + ( + x, + preconDiag_.asSquare(), // Promotes to square + Upper.asSquare(), + Lower.asSquare(), + b + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) + { + if (Upper.activeType() == blockCoeffBase::SCALAR) + { + LUSubstitute + ( + x, + preconDiag_.asLinear(), + Upper.asLinear(), // Promotes to linear + Lower.asLinear(), // Promotes to linear + b + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + LUSubstitute + ( + x, + preconDiag_.asLinear(), + Upper.asLinear(), + Lower.asLinear(), + b + ); + } + else if (Upper.activeType() == blockCoeffBase::SQUARE) + { + LUSubstitute + ( + x, + preconDiag_.asSquare(), // Promotes to square + Upper.asSquare(), + Lower.asSquare(), + b + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::SQUARE) + { + if (Upper.activeType() == blockCoeffBase::SCALAR) + { + LUSubstitute + ( + x, + preconDiag_.asSquare(), + Upper.asSquare(), // Promotes to square + Lower.asSquare(), // Promotes to square + b + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + LUSubstitute + ( + x, + preconDiag_.asSquare(), + Upper.asSquare(), // Promotes to square + Lower.asSquare(), // Promotes to square + b + ); + } + else if (Upper.activeType() == blockCoeffBase::SQUARE) + { + LUSubstitute + ( + x, + preconDiag_.asSquare(), + Upper.asSquare(), + Lower.asSquare(), + b + ); + } + } + else + { + FatalErrorIn + ( + "void BlockILUC0Precon::precondition\n" + "(\n" + " Field& x,\n" + " const Field& T\n" + ") const" + ) << "Problem with coefficient type morphing." + << abort(FatalError); + } + } + else + { + FatalErrorIn + ( + "void BlockILUC0Precon::precondition\n" + "(\n" + " Field& x,\n" + " const Field& b\n" + ") const" + ) << "Unnecessary use of BlockILUC0 preconditioner for diagonal " + << "matrix. " + << nl + << "Use BlockDiagonal preconditioner instead." + << abort(FatalError); + } +} + + +template +void Foam::BlockILUC0Precon::preconditionT +( + Field& xT, + const Field& bT +) const +{ + if (!this->matrix_.diagonal()) + { + // Get upper and lower matrix factors + const CoeffField& Lower = preconLower_; + const CoeffField& Upper = preconUpper_; + + // Execute transpose preconditioning by transpose LU substitution. + // 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) + { + LUSubstituteT + ( + xT, + preconDiag_.asScalar(), + Upper.asScalar(), + Lower.asScalar(), + bT + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + LUSubstituteT + ( + xT, + preconDiag_.asLinear(), // Promotes to linear + Upper.asLinear(), + Lower.asLinear(), + bT + ); + } + else if (Upper.activeType() == blockCoeffBase::SQUARE) + { + LUSubstituteT + ( + xT, + preconDiag_.asSquare(), // Promotes to square + Upper.asSquare(), + Lower.asSquare(), + bT + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) + { + if (Upper.activeType() == blockCoeffBase::SCALAR) + { + LUSubstituteT + ( + xT, + preconDiag_.asLinear(), + Upper.asLinear(), // Promotes to linear + Lower.asLinear(), // Promotes to linear + bT + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + LUSubstituteT + ( + xT, + preconDiag_.asLinear(), + Upper.asLinear(), + Lower.asLinear(), + bT + ); + } + else if (Upper.activeType() == blockCoeffBase::SQUARE) + { + LUSubstituteT + ( + xT, + preconDiag_.asSquare(), // Promotes to square + Upper.asSquare(), + Lower.asSquare(), + bT + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::SQUARE) + { + if (Upper.activeType() == blockCoeffBase::SCALAR) + { + LUSubstituteT + ( + xT, + preconDiag_.asSquare(), + Upper.asSquare(), // Promotes to square + Lower.asSquare(), // Promotes to square + bT + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + LUSubstituteT + ( + xT, + preconDiag_.asSquare(), + Upper.asSquare(), // Promotes to square + Lower.asSquare(), // Promotes to square + bT + ); + } + else if (Upper.activeType() == blockCoeffBase::SQUARE) + { + LUSubstituteT + ( + xT, + preconDiag_.asSquare(), + Upper.asSquare(), + Lower.asSquare(), + bT + ); + } + } + else + { + FatalErrorIn + ( + "void BlockILUC0Precon::preconditionT\n" + "(\n" + " Field& x,\n" + " const Field& T\n" + ") const" + ) << "Problem with coefficient type morphing." + << abort(FatalError); + } + } + else + { + FatalErrorIn + ( + "void BlockILUC0Precon::preconditionT\n" + "(\n" + " Field& x,\n" + " const Field& b\n" + ") const" + ) << "Unnecessary use of BlockILUC0 preconditioner for diagonal " + << "matrix." + << nl + << "Use BlockDiagonal preconditioner instead." + << abort(FatalError); + } +} + + +template +void Foam::BlockILUC0Precon::initMatrix() +{ + preconDiag_ = this->matrix_.diag(); + preconLower_ = this->matrix_.lower(); + preconUpper_ = this->matrix_.upper(); + + this->calcFactorization(); +} + + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0Precon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0Precon.H new file mode 100644 index 000000000..c79d4d1a0 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0Precon.H @@ -0,0 +1,194 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +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 + BlockILUC0Precon + +Description + Block variant of 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 + BlockILUC0Precon.C + BlockILUC0PreconDecoupled.C + +\*---------------------------------------------------------------------------*/ + +#ifndef BlockILUC0Precon_H +#define BlockILUC0Precon_H + +#include "BlockLduPrecon.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class BlockILUC0Precon Declaration +\*---------------------------------------------------------------------------*/ + +template +class BlockILUC0Precon +: + public BlockLduPrecon +{ + // Private Data + + //- Preconditioned diagonal + mutable CoeffField preconDiag_; + + //- Preconditioned lower + mutable CoeffField preconLower_; + + //- Preconditioned upper + mutable CoeffField preconUpper_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + BlockILUC0Precon(const BlockILUC0Precon&); + + //- Disallow default bitwise assignment + void operator=(const BlockILUC0Precon&); + + //- Calculate active type dependant factorization. + // Note: both lower, diag and upper have to be of same type + template + void calcActiveTypeFactorization + ( + Field& preconD, + Field& preconUpper, + Field& preconLower + ) 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 + template + void LUSubstitute + ( + Field& x, + const Field& preconD, + const Field& upper, + const Field& lower, + const Field& b + ) const; + + //- Performs forward/backward LU substitution on transpose system + // Note: both lower, diag and upper have to be of same type + template + void LUSubstituteT + ( + Field& xT, + const Field& preconD, + const Field& upper, + const Field& lower, + const Field& bT + ) const; + + + // Decoupled operations, used in template specialisation + + //- Execute preconditioning, decoupled version + void decoupledPrecondition + ( + Field& x, + const Field& b + ) const; + + //- Execute preconditioning with matrix transpose, + // decoupled version + void decoupledPreconditionT + ( + Field& xT, + const Field& bT + ) const; + + +public: + + //- Runtime type information + TypeName("ILUC0"); + + + // Constructors + + //- Construct from components + BlockILUC0Precon + ( + const BlockLduMatrix& matrix, + const dictionary& dict + ); + + + //- Destructor + virtual ~BlockILUC0Precon(); + + + // Member Functions + + //- Execute preconditioning + virtual void precondition + ( + Field& x, + const Field& b + ) const; + + //- Execute preconditioning with matrix transpose + virtual void preconditionT + ( + Field& xT, + const Field& bT + ) const; + + //- Re-initialise preconditioner after matrix coefficient update + virtual void initMatrix(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "BlockILUC0Precon.C" +# include "BlockILUC0PreconDecoupled.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0PreconDecoupled.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0PreconDecoupled.C new file mode 100644 index 000000000..4c90c4d77 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/BlockILUC0PreconDecoupled.C @@ -0,0 +1,240 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +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 . + +Description + Block variant of ILU preconditioning with arbitrary level of fill in (p), + based on Crout algorithm. Decoupled version. + + Reference: Saad, Y.: Iterative Methods for Sparse Linear Systems + (2nd Edition), SIAM, 2003. + +Author + Vuko Vukcevic, FMENA Zagreb. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#include "BlockILUC0Precon.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +void Foam::BlockILUC0Precon::decoupledPrecondition +( + Field& x, + const Field& b +) const +{ + // Helper type definition of DecoupledCoeffField + typedef DecoupledCoeffField DecoupledTypeCoeffField; + + if (!this->matrix_.diagonal()) + { + // Get upper and lower matrix factors + const DecoupledTypeCoeffField& Lower = preconLower_; + const DecoupledTypeCoeffField& Upper = preconUpper_; + + // Execute preconditioning by LU substitution. + // 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) + { + LUSubstitute + ( + x, + preconDiag_.asScalar(), + Upper.asScalar(), + Lower.asScalar(), + b + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + LUSubstitute + ( + x, + preconDiag_.asLinear(), // Promotes to linear + Upper.asLinear(), + Lower.asLinear(), + b + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) + { + if (Upper.activeType() == blockCoeffBase::SCALAR) + { + LUSubstitute + ( + x, + preconDiag_.asLinear(), + Upper.asLinear(), // Promotes to linear + Lower.asLinear(), // Promotes to linear + b + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + LUSubstitute + ( + x, + preconDiag_.asLinear(), + Upper.asLinear(), + Lower.asLinear(), + b + ); + } + } + else + { + FatalErrorIn + ( + "void BlockILUC0Precon::decoupledPrecondition\n" + "(\n" + " Field& x,\n" + " const Field& T\n" + ") const" + ) << "Problem with coefficient type morphing." + << abort(FatalError); + } + } + else + { + FatalErrorIn + ( + "void BlockILUC0Precon::decoupledPrecondition\n" + "(\n" + " Field& x,\n" + " const Field& b\n" + ") const" + ) << "Unnecessary use of BlockILUC0 preconditioner for diagonal " + << "matrix. " + << nl + << "Use BlockDiagonal preconditioner instead." + << abort(FatalError); + } +} + + +template +void Foam::BlockILUC0Precon::decoupledPreconditionT +( + Field& xT, + const Field& bT +) const +{ + // Helper type definition of DecoupledCoeffField + typedef DecoupledCoeffField DecoupledTypeCoeffField; + + if (!this->matrix_.diagonal()) + { + // Get upper and lower matrix factors + const DecoupledTypeCoeffField& Lower = preconLower_; + const DecoupledTypeCoeffField& Upper = preconUpper_; + + // Execute preconditioning by LU substitution. + // 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) + { + LUSubstituteT + ( + xT, + preconDiag_.asScalar(), + Upper.asScalar(), + Lower.asScalar(), + bT + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + LUSubstituteT + ( + xT, + preconDiag_.asLinear(), // Promotes to linear + Upper.asLinear(), + Lower.asLinear(), + bT + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) + { + if (Upper.activeType() == blockCoeffBase::SCALAR) + { + LUSubstituteT + ( + xT, + preconDiag_.asLinear(), + Upper.asLinear(), // Promotes to linear + Lower.asLinear(), // Promotes to linear + bT + ); + } + else if (Upper.activeType() == blockCoeffBase::LINEAR) + { + LUSubstituteT + ( + xT, + preconDiag_.asLinear(), + Upper.asLinear(), + Lower.asLinear(), + bT + ); + } + } + else + { + FatalErrorIn + ( + "void BlockILUC0Precon::decoupledPreconditionT\n" + "(\n" + " Field& x,\n" + " const Field& T\n" + ") const" + ) << "Problem with coefficient type morphing." + << abort(FatalError); + } + } + else + { + FatalErrorIn + ( + "void BlockILUC0Precon::decoupledPreconditionT\n" + "(\n" + " Field& x,\n" + " const Field& b\n" + ") const" + ) << "Unnecessary use of BlockILUC0 preconditioner for diagonal " + << "matrix. " + << nl + << "Use BlockDiagonal preconditioner instead." + << abort(FatalError); + } +} + + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/blockILUC0Precons.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/blockILUC0Precons.C new file mode 100644 index 000000000..7cd64f517 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/blockILUC0Precons.C @@ -0,0 +1,48 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +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 . + +Description + ILUCp preconditioners. + +Author + Vuko Vukcevic, FMENA Zagreb. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#include "blockLduMatrices.H" +#include "blockLduPrecons.H" +#include "blockILUC0Precons.H" +#include "addToRunTimeSelectionTable.H" + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +makeBlockPrecons(blockILUC0Precon); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/blockILUC0Precons.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/blockILUC0Precons.H new file mode 100644 index 000000000..e42a1ba41 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/blockILUC0Precons.H @@ -0,0 +1,64 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +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 + BlockILUC0Precon + +Description + Typedefs for ILUCp preconditioning + +Author + Vuko Vukcevic, FMENA Zagreb. All rights reserved. + +SourceFiles + blockILUC0Precons.C + +\*---------------------------------------------------------------------------*/ + +#ifndef blockILUC0Precons_H +#define blockILUC0Precons_H + +// Specialisations not implemented. VV, 3/Jul/2015. +#include "scalarBlockILUC0Precon.H" +#include "tensorBlockILUC0Precon.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +typedef BlockILUC0Precon blockILUC0PreconScalar; +typedef BlockILUC0Precon blockILUC0PreconVector; +typedef BlockILUC0Precon blockILUC0PreconTensor; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/scalarBlockILUC0Precon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/scalarBlockILUC0Precon.C new file mode 100644 index 000000000..6d9de0f6d --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/scalarBlockILUC0Precon.C @@ -0,0 +1,276 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +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 + BlockILUC0Precon + +Description + Template specialisation for scalar block ILUCp preconditioning + +Author + Vuko Vukcevic, FMENA Zagreb. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#ifndef scalarBlockILUC0Precon_H +#define scalarBlockILUC0Precon_H + +#include "BlockILUC0Precon.H" +#include "scalarBlockILUC0Precon.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template<> +template<> +void BlockILUC0Precon::calcActiveTypeFactorization +( + scalarField& preconD, + scalarField& preconUpper, + scalarField& preconLower +) 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 lduAddressing& addr = this->matrix_.lduAddr(); + + // Get upper/lower extended addressing + const label* const __restrict__ uPtr = addr.upperAddr().begin(); + const label* const __restrict__ lPtr = addr.lowerAddr().begin(); + + // Get extended owner start addressing + const label* const __restrict__ ownStartPtr = + addr.ownerStartAddr().begin(); + + // Get extended losort and losort start addressing + const label* const __restrict__ lsrPtr = + addr.losortAddr().begin(); + const label* const __restrict__ lsrStartPtr = + addr.losortStartAddr().begin(); + + // Get access to factored matrix entries + scalar* __restrict__ diagPtr = preconD.begin(); + scalar* __restrict__ upperPtr = preconUpper.begin(); + scalar* __restrict__ lowerPtr = preconLower.begin(); + + // Get access to working fields + scalar* __restrict__ zPtr = z.begin(); + scalar* __restrict__ wPtr = w.begin(); + + // Define start and end face ("virtual" face when extended addressing is + // used) of this row/column. + 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/upper coeff loop (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 + { + FatalErrorIn + ( + "template<>\n" + "template<>\n" + "void BlockILUC0Precon::calcFactorization\n" + "(\n" + " scalarField& preconD,\n" + " scalarField& extUpper,\n" + " scalarField& extLower,\n" + " scalarField& zDiag\n," + " scalarField& z,\n" + " scalarField& w,\n" + ") const" + ) << "Unnecessary use of BlockILUC0 preconditioner for diagonal " + << "matrix." + << nl + << "Use BlockDiagonal preconditioner instead." + << abort(FatalError); + } +} + + +template<> +void BlockILUC0Precon::calcFactorization() const +{ + calcActiveTypeFactorization + ( + preconDiag_.asScalar(), + preconUpper_.asScalar(), + preconLower_.asScalar() + ); +} + + +template<> +void BlockILUC0Precon::precondition +( + scalarField& x, + const scalarField& b +) const +{ + // Decoupled version + notImplemented("void Foam::BlockILUC0Precon::precondition"); +} + + +template<> +void BlockILUC0Precon::preconditionT +( + scalarField& xT, + const scalarField& bT +) const +{ + // Decoupled version + notImplemented("void Foam::BlockILUC0Precon::preconditionT"); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/scalarBlockILUC0Precon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/scalarBlockILUC0Precon.H new file mode 100644 index 000000000..7b0e6e62b --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/scalarBlockILUC0Precon.H @@ -0,0 +1,90 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +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 + BlockILUC0Precon + +Description + Template specialisation for scalar block ILUCp preconditioning + +Author + Vuko Vukcevic, FMENA Zagreb. All rights reserved. + +SourceFiles + scalarBlockILUC0Precon.C + +\*---------------------------------------------------------------------------*/ + +#ifndef scalarBlockILUC0Precon_H +#define scalarBlockILUC0Precon_H + +#include "BlockILUC0Precon.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Calculate active type factorization +template<> +template<> +void BlockILUC0Precon::calcActiveTypeFactorization +( + scalarField& preconD, + scalarField& preconUpper, + scalarField& preconLower +) const; + + +// Calculate factorization (constructor helper) +template<> +void BlockILUC0Precon::calcFactorization() const; + + +// Precondition +template<> +void BlockILUC0Precon::precondition +( + scalarField& x, + const scalarField& b +) const; + + +// Precondition transpose +template<> +void BlockILUC0Precon::preconditionT +( + scalarField& xT, + const scalarField& bT +) const; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/tensorBlockILUC0Precon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/tensorBlockILUC0Precon.C new file mode 100644 index 000000000..b123dc7de --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/tensorBlockILUC0Precon.C @@ -0,0 +1,147 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +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 + BlockILUC0Precon + +Description + Template specialisation for tensor block ILUCp preconditioning + +Author + Vuko Vukcevic, FMENA Zagreb. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#ifndef tensorBlockILUC0Precon_H +#define tensorBlockILUC0Precon_H + +#include "BlockILUC0Precon.H" +#include "tensorBlockILUC0Precon.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template<> +template<> +void BlockILUC0Precon::calcActiveTypeFactorization +( + tensorField& preconD, + tensorField& preconUpper, + tensorField& preconLower +) const +{ + // Decoupled version + notImplemented("void Foam::BlockILUC0Precon::calcFactorization"); +} + + +template<> +void BlockILUC0Precon::calcFactorization() const +{ + // Get upper and lower matrix factors + CoeffField& Lower = preconLower_; + CoeffField& Upper = preconUpper_; + + // 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 BlockILUC0Precon::precondition +( + tensorField& x, + const tensorField& b +) const +{ + // Decoupled version + notImplemented("void Foam::BlockILUC0Precon::precondition"); +} + + +template<> +void BlockILUC0Precon::preconditionT +( + tensorField& xT, + const tensorField& bT +) const +{ + // Decoupled version + notImplemented("void Foam::BlockILUC0Precon::preconditionT"); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/tensorBlockILUC0Precon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/tensorBlockILUC0Precon.H new file mode 100644 index 000000000..a21eb07df --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUC0Precon/tensorBlockILUC0Precon.H @@ -0,0 +1,90 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +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 + BlockILUC0Precon + +Description + Template specialisation for tensor block ILUCp preconditioning + +Author + Vuko Vukcevic, FMENA Zagreb. All rights reserved. + +SourceFiles + tensorBlockILUC0Precon.C + +\*---------------------------------------------------------------------------*/ + +#ifndef tensorBlockILUC0Precon_H +#define tensorBlockILUC0Precon_H + +#include "BlockILUC0Precon.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Calculate active type factorization +template<> +template<> +void BlockILUC0Precon::calcActiveTypeFactorization +( + tensorField& preconD, + tensorField& preconUpper, + tensorField& preconLower +) const; + + +// Calculate factorization (constructor helper) +template<> +void BlockILUC0Precon::calcFactorization() const; + + +// Precondition +template<> +void BlockILUC0Precon::precondition +( + tensorField& x, + const tensorField& b +) const; + + +// Precondition transpose +template<> +void BlockILUC0Precon::preconditionT +( + tensorField& xT, + const tensorField& bT +) const; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUC0Smoother/BlockILUC0Smoother.H b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUC0Smoother/BlockILUC0Smoother.H new file mode 100644 index 000000000..6644cbf07 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUC0Smoother/BlockILUC0Smoother.H @@ -0,0 +1,147 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +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 + BlockILUC0Smoother + +Description + ILU Cp smoother + +SourceFiles + blockILUC0Smoothers.C + +\*---------------------------------------------------------------------------*/ + +#ifndef BlockILUC0Smoother_H +#define BlockILUC0Smoother_H + +#include "BlockLduSmoother.H" +#include "blockILUC0Precons.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class BlockILUC0Smoother Declaration +\*---------------------------------------------------------------------------*/ + +template +class BlockILUC0Smoother +: + public BlockLduSmoother +{ + // Private Data + + //- ILUCp preconditioner + BlockILUC0Precon precon_; + + //- Correction array + mutable Field xCorr_; + + //- Residual array + mutable Field residual_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + BlockILUC0Smoother(const BlockILUC0Smoother&); + + //- Disallow default bitwise assignment + void operator=(const BlockILUC0Smoother&); + + +public: + + //- Runtime type information + TypeName("ILUC0"); + + + // Constructors + + //- Construct from components + BlockILUC0Smoother + ( + const BlockLduMatrix& matrix, + const dictionary& dict + ) + : + BlockLduSmoother(matrix), + precon_(matrix, dict), + xCorr_(matrix.lduAddr().size()), + residual_(matrix.lduAddr().size()) + {} + + + //- Destructor + virtual ~BlockILUC0Smoother() + {} + + + // Member Functions + + //- Execute smoothing + virtual void smooth + ( + Field& x, + const Field& b, + const label nSweeps + ) const + { + for (label sweep = 0; sweep < nSweeps; sweep++) + { + // Calculate residual + this-> matrix_.Amul(residual_, x); + + // residual = b - Ax + forAll (b, i) + { + residual_[i] = b[i] - residual_[i]; + } + + precon_.precondition(xCorr_, residual_); + + // Add correction to x + x += xCorr_; + } + } + + //- Re-initialise preconditioner after matrix coefficient update + virtual void initMatrix() + { + precon_.initMatrix(); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUC0Smoother/blockILUC0Smoothers.C b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUC0Smoother/blockILUC0Smoothers.C new file mode 100644 index 000000000..312b801e6 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUC0Smoother/blockILUC0Smoothers.C @@ -0,0 +1,42 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#include "blockLduMatrices.H" +#include "blockLduSmoothers.H" +#include "blockILUC0Smoothers.H" +#include "addToRunTimeSelectionTable.H" + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +makeBlockSmoothers(blockILUC0Smoother); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUC0Smoother/blockILUC0Smoothers.H b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUC0Smoother/blockILUC0Smoothers.H new file mode 100644 index 000000000..26cd39459 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUC0Smoother/blockILUC0Smoothers.H @@ -0,0 +1,62 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +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 + BlockILUC0Smoother + +Description + Typedefs for Incomplete Lower-Upper (ILU) with Cp fill-in smoother + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved. + +SourceFiles + blockILUC0Smoothers.C + +\*---------------------------------------------------------------------------*/ + +#ifndef blockILUC0Smoothers_H +#define blockILUC0Smoothers_H + +#include "BlockILUC0Smoother.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +typedef BlockILUC0Smoother blockILUC0SmootherScalar; +typedef BlockILUC0Smoother blockILUC0SmootherVector; +typedef BlockILUC0Smoother blockILUC0SmootherTensor; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C index 7d472d67d..3597145b5 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C @@ -31,11 +31,13 @@ License #include "blockDiagonalPrecons.H" #include "blockGaussSeidelPrecons.H" #include "BlockCholeskyPrecon.H" +#include "blockILUC0Precons.H" #include "BlockILUCpPrecon.H" #include "blockLduSmoothers.H" #include "blockGaussSeidelSmoothers.H" #include "BlockILUSmoother.H" +#include "blockILUC0Smoothers.H" #include "BlockILUCpSmoother.H" #include "blockLduSolvers.H" @@ -87,6 +89,9 @@ makeBlockPrecon(block##Type##Precon, block##Type##GaussSeidelPrecon); \ typedef BlockCholeskyPrecon block##Type##CholeskyPrecon; \ makeBlockPrecon(block##Type##Precon, block##Type##CholeskyPrecon); \ \ +typedef BlockILUC0Precon block##Type##ILUC0Precon; \ +makeBlockPrecon(block##Type##Precon, block##Type##ILUC0Precon); \ + \ typedef BlockILUCpPrecon block##Type##ILUCpPrecon; \ makeBlockPrecon(block##Type##Precon, block##Type##ILUCpPrecon); \ \ @@ -101,6 +106,9 @@ makeBlockSmoother(block##Type##Smoother, block##Type##GaussSeidelSmoother); \ typedef BlockILUSmoother block##Type##ILUSmoother; \ makeBlockSmoother(block##Type##Smoother, block##Type##ILUSmoother); \ \ +typedef BlockILUC0Smoother block##Type##ILUC0Smoother; \ +makeBlockSmoother(block##Type##Smoother, block##Type##ILUC0Smoother); \ + \ typedef BlockILUCpSmoother block##Type##ILUCpSmoother; \ makeBlockSmoother(block##Type##Smoother, block##Type##ILUCpSmoother); \ \