Block ILUC0 precon and smoother
This commit is contained in:
Hrvoje Jasak
parent
737dee5705
commit
6c5d254802
13 changed files with 2250 additions and 0 deletions
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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<class Type>
|
||||
template<class LDUType>
|
||||
void Foam::BlockILUC0Precon<Type>::calcActiveTypeFactorization
|
||||
(
|
||||
Field<LDUType>& preconD,
|
||||
Field<LDUType>& preconUpper,
|
||||
Field<LDUType>& preconLower
|
||||
) const
|
||||
{
|
||||
if (!this->matrix_.diagonal())
|
||||
{
|
||||
// Get number of rows
|
||||
const label nRows = preconD.size();
|
||||
|
||||
// Allocate working fields
|
||||
Field<LDUType> z(nRows, pTraits<LDUType>::zero);
|
||||
Field<LDUType> w(nRows, pTraits<LDUType>::zero);
|
||||
LDUType zDiag(pTraits<LDUType>::zero);
|
||||
|
||||
// Create multiplication function object
|
||||
typename BlockCoeff<Type>::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<LDUType>::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<LDUType>::zero;
|
||||
wPtr[uPtr[faceI]] = pTraits<LDUType>::zero;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"template <class Type>\n"
|
||||
"template <class LDUType>\n"
|
||||
"void BlockILUC0Precon<Type>::calcActiveTypeFactorization\n"
|
||||
"(\n"
|
||||
" Field<LDUType>& preconD,\n"
|
||||
" Field<LDUType>& preconUpper,\n"
|
||||
" Field<LDUType>& preconLower,\n"
|
||||
") const"
|
||||
) << "Unnecessary use of BlockILUC0 preconditioner for diagonal "
|
||||
<< "matrix."
|
||||
<< nl
|
||||
<< "Use BlockDiagonal preconditioner instead."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
void Foam::BlockILUC0Precon<Type>::calcFactorization() const
|
||||
{
|
||||
// Get upper and lower matrix factors
|
||||
CoeffField<Type>& Lower = preconLower_;
|
||||
CoeffField<Type>& 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<class Type>
|
||||
template<class LDUType>
|
||||
void Foam::BlockILUC0Precon<Type>::LUSubstitute
|
||||
(
|
||||
Field<Type>& x,
|
||||
const Field<LDUType>& preconD,
|
||||
const Field<LDUType>& upper,
|
||||
const Field<LDUType>& lower,
|
||||
const Field<Type>& b
|
||||
) const
|
||||
{
|
||||
// Create multiplication function object
|
||||
typename BlockCoeff<Type>::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<class Type>
|
||||
template<class LDUType>
|
||||
void Foam::BlockILUC0Precon<Type>::LUSubstituteT
|
||||
(
|
||||
Field<Type>& xT,
|
||||
const Field<LDUType>& preconD,
|
||||
const Field<LDUType>& upper,
|
||||
const Field<LDUType>& lower,
|
||||
const Field<Type>& bT
|
||||
) const
|
||||
{
|
||||
// Create multiplication function object
|
||||
typename BlockCoeff<Type>::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<class Type>
|
||||
Foam::BlockILUC0Precon<Type>::BlockILUC0Precon
|
||||
(
|
||||
const BlockLduMatrix<Type>& matrix,
|
||||
const dictionary& dict
|
||||
)
|
||||
:
|
||||
BlockLduPrecon<Type>(matrix),
|
||||
preconDiag_(matrix.diag()),
|
||||
preconLower_(matrix.lower()),
|
||||
preconUpper_(matrix.upper())
|
||||
{
|
||||
this->calcFactorization();
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
template<class Type>
|
||||
Foam::BlockILUC0Precon<Type>::~BlockILUC0Precon()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
template<class Type>
|
||||
void Foam::BlockILUC0Precon<Type>::precondition
|
||||
(
|
||||
Field<Type>& x,
|
||||
const Field<Type>& b
|
||||
) const
|
||||
{
|
||||
if (!this->matrix_.diagonal())
|
||||
{
|
||||
// Get upper and lower matrix factors
|
||||
const CoeffField<Type>& Lower = preconLower_;
|
||||
const CoeffField<Type>& 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<Type>::precondition\n"
|
||||
"(\n"
|
||||
" Field<Type>& x,\n"
|
||||
" const Field<Type>& T\n"
|
||||
") const"
|
||||
) << "Problem with coefficient type morphing."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"void BlockILUC0Precon<Type>::precondition\n"
|
||||
"(\n"
|
||||
" Field<Type>& x,\n"
|
||||
" const Field<Type>& b\n"
|
||||
") const"
|
||||
) << "Unnecessary use of BlockILUC0 preconditioner for diagonal "
|
||||
<< "matrix. "
|
||||
<< nl
|
||||
<< "Use BlockDiagonal preconditioner instead."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
void Foam::BlockILUC0Precon<Type>::preconditionT
|
||||
(
|
||||
Field<Type>& xT,
|
||||
const Field<Type>& bT
|
||||
) const
|
||||
{
|
||||
if (!this->matrix_.diagonal())
|
||||
{
|
||||
// Get upper and lower matrix factors
|
||||
const CoeffField<Type>& Lower = preconLower_;
|
||||
const CoeffField<Type>& 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<Type>::preconditionT\n"
|
||||
"(\n"
|
||||
" Field<Type>& x,\n"
|
||||
" const Field<Type>& T\n"
|
||||
") const"
|
||||
) << "Problem with coefficient type morphing."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"void BlockILUC0Precon<Type>::preconditionT\n"
|
||||
"(\n"
|
||||
" Field<Type>& x,\n"
|
||||
" const Field<Type>& b\n"
|
||||
") const"
|
||||
) << "Unnecessary use of BlockILUC0 preconditioner for diagonal "
|
||||
<< "matrix."
|
||||
<< nl
|
||||
<< "Use BlockDiagonal preconditioner instead."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
void Foam::BlockILUC0Precon<Type>::initMatrix()
|
||||
{
|
||||
preconDiag_ = this->matrix_.diag();
|
||||
preconLower_ = this->matrix_.lower();
|
||||
preconUpper_ = this->matrix_.upper();
|
||||
|
||||
this->calcFactorization();
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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 Type>
|
||||
class BlockILUC0Precon
|
||||
:
|
||||
public BlockLduPrecon<Type>
|
||||
{
|
||||
// Private Data
|
||||
|
||||
//- Preconditioned diagonal
|
||||
mutable CoeffField<Type> preconDiag_;
|
||||
|
||||
//- Preconditioned lower
|
||||
mutable CoeffField<Type> preconLower_;
|
||||
|
||||
//- Preconditioned upper
|
||||
mutable CoeffField<Type> 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<class LDUType>
|
||||
void calcActiveTypeFactorization
|
||||
(
|
||||
Field<LDUType>& preconD,
|
||||
Field<LDUType>& preconUpper,
|
||||
Field<LDUType>& 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<class LDUType>
|
||||
void LUSubstitute
|
||||
(
|
||||
Field<Type>& x,
|
||||
const Field<LDUType>& preconD,
|
||||
const Field<LDUType>& upper,
|
||||
const Field<LDUType>& lower,
|
||||
const Field<Type>& b
|
||||
) const;
|
||||
|
||||
//- Performs forward/backward LU substitution on transpose system
|
||||
// Note: both lower, diag and upper have to be of same type
|
||||
template<class LDUType>
|
||||
void LUSubstituteT
|
||||
(
|
||||
Field<Type>& xT,
|
||||
const Field<LDUType>& preconD,
|
||||
const Field<LDUType>& upper,
|
||||
const Field<LDUType>& lower,
|
||||
const Field<Type>& bT
|
||||
) const;
|
||||
|
||||
|
||||
// Decoupled operations, used in template specialisation
|
||||
|
||||
//- Execute preconditioning, decoupled version
|
||||
void decoupledPrecondition
|
||||
(
|
||||
Field<Type>& x,
|
||||
const Field<Type>& b
|
||||
) const;
|
||||
|
||||
//- Execute preconditioning with matrix transpose,
|
||||
// decoupled version
|
||||
void decoupledPreconditionT
|
||||
(
|
||||
Field<Type>& xT,
|
||||
const Field<Type>& bT
|
||||
) const;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
//- Runtime type information
|
||||
TypeName("ILUC0");
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from components
|
||||
BlockILUC0Precon
|
||||
(
|
||||
const BlockLduMatrix<Type>& matrix,
|
||||
const dictionary& dict
|
||||
);
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~BlockILUC0Precon();
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Execute preconditioning
|
||||
virtual void precondition
|
||||
(
|
||||
Field<Type>& x,
|
||||
const Field<Type>& b
|
||||
) const;
|
||||
|
||||
//- Execute preconditioning with matrix transpose
|
||||
virtual void preconditionT
|
||||
(
|
||||
Field<Type>& xT,
|
||||
const Field<Type>& 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
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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<class Type>
|
||||
void Foam::BlockILUC0Precon<Type>::decoupledPrecondition
|
||||
(
|
||||
Field<Type>& x,
|
||||
const Field<Type>& b
|
||||
) const
|
||||
{
|
||||
// Helper type definition of DecoupledCoeffField
|
||||
typedef DecoupledCoeffField<Type> 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<Type>::decoupledPrecondition\n"
|
||||
"(\n"
|
||||
" Field<Type>& x,\n"
|
||||
" const Field<Type>& T\n"
|
||||
") const"
|
||||
) << "Problem with coefficient type morphing."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"void BlockILUC0Precon<Type>::decoupledPrecondition\n"
|
||||
"(\n"
|
||||
" Field<Type>& x,\n"
|
||||
" const Field<Type>& b\n"
|
||||
") const"
|
||||
) << "Unnecessary use of BlockILUC0 preconditioner for diagonal "
|
||||
<< "matrix. "
|
||||
<< nl
|
||||
<< "Use BlockDiagonal preconditioner instead."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
void Foam::BlockILUC0Precon<Type>::decoupledPreconditionT
|
||||
(
|
||||
Field<Type>& xT,
|
||||
const Field<Type>& bT
|
||||
) const
|
||||
{
|
||||
// Helper type definition of DecoupledCoeffField
|
||||
typedef DecoupledCoeffField<Type> 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<Type>::decoupledPreconditionT\n"
|
||||
"(\n"
|
||||
" Field<Type>& x,\n"
|
||||
" const Field<Type>& T\n"
|
||||
") const"
|
||||
) << "Problem with coefficient type morphing."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"void BlockILUC0Precon<Type>::decoupledPreconditionT\n"
|
||||
"(\n"
|
||||
" Field<Type>& x,\n"
|
||||
" const Field<Type>& b\n"
|
||||
") const"
|
||||
) << "Unnecessary use of BlockILUC0 preconditioner for diagonal "
|
||||
<< "matrix. "
|
||||
<< nl
|
||||
<< "Use BlockDiagonal preconditioner instead."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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<scalar> blockILUC0PreconScalar;
|
||||
typedef BlockILUC0Precon<vector> blockILUC0PreconVector;
|
||||
typedef BlockILUC0Precon<tensor> blockILUC0PreconTensor;
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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<scalar>::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<scalar>::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<scalar>::calcFactorization() const
|
||||
{
|
||||
calcActiveTypeFactorization
|
||||
(
|
||||
preconDiag_.asScalar(),
|
||||
preconUpper_.asScalar(),
|
||||
preconLower_.asScalar()
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
template<>
|
||||
void BlockILUC0Precon<scalar>::precondition
|
||||
(
|
||||
scalarField& x,
|
||||
const scalarField& b
|
||||
) const
|
||||
{
|
||||
// Decoupled version
|
||||
notImplemented("void Foam::BlockILUC0Precon<scalar>::precondition");
|
||||
}
|
||||
|
||||
|
||||
template<>
|
||||
void BlockILUC0Precon<scalar>::preconditionT
|
||||
(
|
||||
scalarField& xT,
|
||||
const scalarField& bT
|
||||
) const
|
||||
{
|
||||
// Decoupled version
|
||||
notImplemented("void Foam::BlockILUC0Precon<scalar>::preconditionT");
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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<scalar>::calcActiveTypeFactorization
|
||||
(
|
||||
scalarField& preconD,
|
||||
scalarField& preconUpper,
|
||||
scalarField& preconLower
|
||||
) const;
|
||||
|
||||
|
||||
// Calculate factorization (constructor helper)
|
||||
template<>
|
||||
void BlockILUC0Precon<scalar>::calcFactorization() const;
|
||||
|
||||
|
||||
// Precondition
|
||||
template<>
|
||||
void BlockILUC0Precon<scalar>::precondition
|
||||
(
|
||||
scalarField& x,
|
||||
const scalarField& b
|
||||
) const;
|
||||
|
||||
|
||||
// Precondition transpose
|
||||
template<>
|
||||
void BlockILUC0Precon<scalar>::preconditionT
|
||||
(
|
||||
scalarField& xT,
|
||||
const scalarField& bT
|
||||
) const;
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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<tensor>::calcActiveTypeFactorization
|
||||
(
|
||||
tensorField& preconD,
|
||||
tensorField& preconUpper,
|
||||
tensorField& preconLower
|
||||
) const
|
||||
{
|
||||
// Decoupled version
|
||||
notImplemented("void Foam::BlockILUC0Precon<tensor>::calcFactorization");
|
||||
}
|
||||
|
||||
|
||||
template<>
|
||||
void BlockILUC0Precon<tensor>::calcFactorization() const
|
||||
{
|
||||
// Get upper and lower matrix factors
|
||||
CoeffField<tensor>& Lower = preconLower_;
|
||||
CoeffField<tensor>& 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<tensor>::precondition
|
||||
(
|
||||
tensorField& x,
|
||||
const tensorField& b
|
||||
) const
|
||||
{
|
||||
// Decoupled version
|
||||
notImplemented("void Foam::BlockILUC0Precon<tensor>::precondition");
|
||||
}
|
||||
|
||||
|
||||
template<>
|
||||
void BlockILUC0Precon<tensor>::preconditionT
|
||||
(
|
||||
tensorField& xT,
|
||||
const tensorField& bT
|
||||
) const
|
||||
{
|
||||
// Decoupled version
|
||||
notImplemented("void Foam::BlockILUC0Precon<tensor>::preconditionT");
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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<tensor>::calcActiveTypeFactorization
|
||||
(
|
||||
tensorField& preconD,
|
||||
tensorField& preconUpper,
|
||||
tensorField& preconLower
|
||||
) const;
|
||||
|
||||
|
||||
// Calculate factorization (constructor helper)
|
||||
template<>
|
||||
void BlockILUC0Precon<tensor>::calcFactorization() const;
|
||||
|
||||
|
||||
// Precondition
|
||||
template<>
|
||||
void BlockILUC0Precon<tensor>::precondition
|
||||
(
|
||||
tensorField& x,
|
||||
const tensorField& b
|
||||
) const;
|
||||
|
||||
|
||||
// Precondition transpose
|
||||
template<>
|
||||
void BlockILUC0Precon<tensor>::preconditionT
|
||||
(
|
||||
tensorField& xT,
|
||||
const tensorField& bT
|
||||
) const;
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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 Type>
|
||||
class BlockILUC0Smoother
|
||||
:
|
||||
public BlockLduSmoother<Type>
|
||||
{
|
||||
// Private Data
|
||||
|
||||
//- ILUCp preconditioner
|
||||
BlockILUC0Precon<Type> precon_;
|
||||
|
||||
//- Correction array
|
||||
mutable Field<Type> xCorr_;
|
||||
|
||||
//- Residual array
|
||||
mutable Field<Type> 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<Type>& matrix,
|
||||
const dictionary& dict
|
||||
)
|
||||
:
|
||||
BlockLduSmoother<Type>(matrix),
|
||||
precon_(matrix, dict),
|
||||
xCorr_(matrix.lduAddr().size()),
|
||||
residual_(matrix.lduAddr().size())
|
||||
{}
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~BlockILUC0Smoother()
|
||||
{}
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Execute smoothing
|
||||
virtual void smooth
|
||||
(
|
||||
Field<Type>& x,
|
||||
const Field<Type>& 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
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "blockLduMatrices.H"
|
||||
#include "blockLduSmoothers.H"
|
||||
#include "blockILUC0Smoothers.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
makeBlockSmoothers(blockILUC0Smoother);
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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<scalar> blockILUC0SmootherScalar;
|
||||
typedef BlockILUC0Smoother<vector> blockILUC0SmootherVector;
|
||||
typedef BlockILUC0Smoother<tensor> blockILUC0SmootherTensor;
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -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<type > block##Type##CholeskyPrecon; \
|
||||
makeBlockPrecon(block##Type##Precon, block##Type##CholeskyPrecon); \
|
||||
\
|
||||
typedef BlockILUC0Precon<type > block##Type##ILUC0Precon; \
|
||||
makeBlockPrecon(block##Type##Precon, block##Type##ILUC0Precon); \
|
||||
\
|
||||
typedef BlockILUCpPrecon<type > block##Type##ILUCpPrecon; \
|
||||
makeBlockPrecon(block##Type##Precon, block##Type##ILUCpPrecon); \
|
||||
\
|
||||
|
|
@ -101,6 +106,9 @@ makeBlockSmoother(block##Type##Smoother, block##Type##GaussSeidelSmoother); \
|
|||
typedef BlockILUSmoother<type > block##Type##ILUSmoother; \
|
||||
makeBlockSmoother(block##Type##Smoother, block##Type##ILUSmoother); \
|
||||
\
|
||||
typedef BlockILUC0Smoother<type > block##Type##ILUC0Smoother; \
|
||||
makeBlockSmoother(block##Type##Smoother, block##Type##ILUC0Smoother); \
|
||||
\
|
||||
typedef BlockILUCpSmoother<type > block##Type##ILUCpSmoother; \
|
||||
makeBlockSmoother(block##Type##Smoother, block##Type##ILUCpSmoother); \
|
||||
\
|
||||
|
|
|
|||
Reference in a new issue