exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Second version of BlockILUCp preconditioner

Browse source
This commit is contained in:
Vuko Vukcevic 2015-07-08 08:38:47 +02:00 committed by Hrvoje Jasak
parent b71eb3340c
commit 114f462896
5 changed files with 198 additions and 31 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

190
src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C
View file

@ -30,7 +30,7 @@ License
template<class Type>
Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
(
const blockLduMatrix<Type>& blockLdum,
const BlockLduMatrix<Type>& blockLdum,
const label extensionLevel,
const polyMesh& polyMesh
)
@ -68,24 +68,74 @@ Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
const unallocLabelList& faceMap = extLduAddr_.faceMap();
// Avoid assuming it's upper if the matrix is symmetric
if (blockLdum.hasUpper())
if (blockLdum.thereIsUpper())
{
// Allocate extended upper only
extendedUpperPtr_ = new TypeCoeffField
(
extLduAddr_.extendedUpperAddr().size(),
pTraits<Type>::zero
extLduAddr_.extendedUpperAddr().size()
);
TypeCoeffField& extUpper = *extendedUpperPtr_;
// Get upper coeffs from underlying lduMatrix
const TypeCoeffField& upper = blockLdum.upper();
if (upper.activeType() == blockCoeffBase::SCALAR)
{
// Helper type definition
typedef typename CoeffField<Type>::scalarTypeField activeType;
// Get references to fields
const activeType& activeUpper = upper.asScalar();
activeType& activeExtUpper = extUpper.asScalar();
// Copy non-zero coeffs from basic lduMatrix into corresponding
// positions
forAll (upper, faceI)
{
extUpper[faceMap[faceI]] = upper[faceI];
activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
}
}
else if (upper.activeType() == blockCoeffBase::LINEAR)
{
// Helper type definition
typedef typename CoeffField<Type>::linearTypeField activeType;
// Get references to fields
const activeType& activeUpper = upper.asLinear();
activeType& activeExtUpper = extUpper.asLinear();
// Copy non-zero coeffs from basic lduMatrix into corresponding
// positions
forAll (upper, faceI)
{
activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
}
}
else if (upper.activeType() == blockCoeffBase::SQUARE)
{
// Helper type definition
typedef typename CoeffField<Type>::squareTypeField activeType;
// Get references to fields
const activeType& activeUpper = upper.asSquare();
activeType& activeExtUpper = extUpper.asSquare();
// Copy non-zero coeffs from basic lduMatrix into corresponding
// positions
forAll (upper, faceI)
{
activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
}
}
else
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
) << "Problem between ordinary block matrix and extended"
<< " block matrix upper coeffs type morphing."
<< abort(FatalError);
}
}
else
@ -93,19 +143,69 @@ Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
// Allocate extended lower only
extendedLowerPtr_ = new TypeCoeffField
(
extLduAddr_.extendedLowerAddr().size(),
pTraits<Type>::zero
extLduAddr_.extendedLowerAddr().size()
);
TypeCoeffField& extLower = *extendedLowerPtr_;
// Get lower coeffs from underlying lduMatrix
const TypeCoeffField& lower = blockLdum.lower();
if (lower.activeType() == blockCoeffBase::SCALAR)
{
// Helper type definition
typedef typename CoeffField<Type>::scalarTypeField activeType;
// Get references to fields
const activeType& activeLower = lower.asScalar();
activeType& activeExtLower = extLower.asScalar();
// Copy non-zero coeffs from basic lduMatrix into corresponding
// positions
forAll (lower, faceI)
{
extLower[faceMap[faceI]] = lower[faceI];
activeExtLower[faceMap[faceI]] = activeLower[faceI];
}
}
else if (lower.activeType() == blockCoeffBase::LINEAR)
{
// Helper type definition
typedef typename CoeffField<Type>::linearTypeField activeType;
// Get references to fields
const activeType& activeLower = lower.asLinear();
activeType& activeExtLower = extLower.asLinear();
// Copy non-zero coeffs from basic lduMatrix into corresponding
// positions
forAll (lower, faceI)
{
activeExtLower[faceMap[faceI]] = activeLower[faceI];
}
}
else if (lower.activeType() == blockCoeffBase::SQUARE)
{
// Helper type definition
typedef typename CoeffField<Type>::squareTypeField activeType;
// Get references to fields
const activeType& activeLower = lower.asSquare();
activeType& activeExtLower = extLower.asSquare();
// Copy non-zero coeffs from basic lduMatrix into corresponding
// positions
forAll (lower, faceI)
{
activeExtLower[faceMap[faceI]] = activeLower[faceI];
}
}
else
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
) << "Problem between ordinary block matrix and extended"
<< " block matrix lower coeffs type morphing."
<< abort(FatalError);
}
}
}
@ -118,22 +218,82 @@ Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
const label nExtFaces = extLduAddr_.extendedUpperAddr().size();
// Allocate extended upper and lower
extendedUpperPtr_ = new TypeCoeffField(nExtFaces, pTraits<Type>::zero);
extendedUpperPtr_ = new TypeCoeffField(nExtFaces);
TypeCoeffField& extUpper = *extendedUpperPtr_;
extendedLowerPtr_ = new TypeCoeffField(nExtFaces, pTraits<Type>::zero);
extendedLowerPtr_ = new TypeCoeffField(nExtFaces);
TypeCoeffField& extLower = *extendedLowerPtr_;
// Get upper and lower coeffs from underlying lduMatrix
const TypeCoeffField& upper = blockLdum.upper();
const TypeCoeffField& lower = blockLdum.lower();
// Assuming lower and upper have the same active type
if (upper.activeType() == blockCoeffBase::SCALAR)
{
// Helper type definition
typedef typename CoeffField<Type>::scalarTypeField activeType;
// Get references to fields
const activeType& activeUpper = upper.asScalar();
activeType& activeExtUpper = extUpper.asScalar();
const activeType& activeLower = lower.asScalar();
activeType& activeExtLower = extLower.asScalar();
// Copy non-zero coeffs from basic lduMatrix into corresponding
// positions
forAll (upper, faceI)
{
extUpper[faceMap[faceI]] = upper[faceI];
extLower[faceMap[faceI]] = lower[faceI];
activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
activeExtLower[faceMap[faceI]] = activeLower[faceI];
}
}
else if (upper.activeType() == blockCoeffBase::LINEAR)
{
// Helper type definition
typedef typename CoeffField<Type>::linearTypeField activeType;
// Get references to fields
const activeType& activeUpper = upper.asLinear();
activeType& activeExtUpper = extUpper.asLinear();
const activeType& activeLower = lower.asLinear();
activeType& activeExtLower = extLower.asLinear();
// Copy non-zero coeffs from basic lduMatrix into corresponding
// positions
forAll (upper, faceI)
{
activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
activeExtLower[faceMap[faceI]] = activeLower[faceI];
}
}
else if (upper.activeType() == blockCoeffBase::SQUARE)
{
// Helper type definition
typedef typename CoeffField<Type>::squareTypeField activeType;
// Get references to fields
const activeType& activeUpper = upper.asSquare();
activeType& activeExtUpper = extUpper.asSquare();
const activeType& activeLower = lower.asSquare();
activeType& activeExtLower = extLower.asSquare();
// Copy non-zero coeffs from basic lduMatrix into corresponding
// positions
forAll (upper, faceI)
{
activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
activeExtLower[faceMap[faceI]] = activeLower[faceI];
}
}
else
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
) << "Problem between ordinary block matrix and extended"
<< " block matrix upper/lower coeffs type morphing."
<< abort(FatalError);
}
}
}
@ -172,8 +332,7 @@ Foam::extendedBlockLduMatrix<Type>::extendedLower()
{
extendedLowerPtr_ = new TypeCoeffField
(
extLduAddr_.extendedLowerAddr().size(),
pTraits<Type>::zero
extLduAddr_.extendedLowerAddr().size()
);
}
}
@ -196,8 +355,7 @@ Foam::extendedBlockLduMatrix<Type>::extendedUpper()
{
extendedUpperPtr_ = new TypeCoeffField
(
extLduAddr_.extendedUpperAddr().size(),
pTraits<Type>::zero
extLduAddr_.extendedUpperAddr().size()
);
}
}

8
src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H
View file

@ -58,7 +58,7 @@ class extendedBlockLduMatrix
{
// Public data type
typedef typename BlockLduMatrix<Type>::TypeCoeffField TypeCoeffField;
typedef CoeffField<Type> TypeCoeffField;
private:
@ -145,6 +145,12 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "extendedBlockLduMatrix.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

5
src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C
View file

@ -54,7 +54,7 @@ void Foam::BlockILUCpPrecon<Type>::calcFactorization
if (!this->matrix_.diagonal())
{
// Create multiplication function object
typename BlockCoeff<Type>::multiply mult;
typename BlockCoeff<LDUType>::multiply mult;
// Get necessary const access to extended ldu addressing
const extendedLduAddressing& addr = extBlockMatrix_.extendedLduAddr();
@ -98,7 +98,6 @@ void Foam::BlockILUCpPrecon<Type>::calcFactorization
// Start and end of k-th row (upper) and k-th column (lower)
fStart = ownStartPtr[rowI];
fEnd = ownStartPtr[rowI + 1];
// Initialize temporary working diagonal
zDiagI = diagPtr[rowI];
@ -172,7 +171,7 @@ void Foam::BlockILUCpPrecon<Type>::calcFactorization
// Update diagonal entry, inverting it for future use
LDUType& diagRowI = diagPtr[rowI];
diagRowI = mult.inv(zDiagI);
diagRowI = mult.inverse(zDiagI);
// Index for updating L and U
register label zwIndex;

4
src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.C
View file

@ -39,9 +39,9 @@ namespace Foam
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// Excluded tensor and vector, need reduced versions. VV, 3/Jul/2015.
// Excluded tensor, need reduced version. VV, 3/Jul/2015.
makeBlockPrecon(blockScalarPrecon, blockILUCpPreconScalar);
//makeBlockPrecon(blockVectorPrecon, blockILUCpPreconVector);
makeBlockPrecon(blockVectorPrecon, blockILUCpPreconVector);
//makeBlockPrecon(blockTensorPrecon, blockILUCpPreconTensor);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

4
src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C
View file

@ -31,6 +31,7 @@ License
#include "blockDiagonalPrecons.H"
#include "blockGaussSeidelPrecons.H"
#include "BlockCholeskyPrecon.H"
#include "BlockILUCpPrecon.H"
#include "blockLduSmoothers.H"
#include "blockGaussSeidelSmoothers.H"
@ -83,6 +84,9 @@ makeBlockPrecon(block##Type##Precon, block##Type##GaussSeidelPrecon); \
typedef BlockCholeskyPrecon<type > block##Type##CholeskyPrecon; \
makeBlockPrecon(block##Type##Precon, block##Type##CholeskyPrecon); \
\
typedef BlockILUCpPrecon<type > block##Type##ILUCpPrecon; \
makeBlockPrecon(block##Type##Precon, block##Type##ILUCpPrecon); \
\
/* Smoothers */ \
typedef BlockLduSmoother<type > block##Type##Smoother; \
defineNamedTemplateTypeNameAndDebug(block##Type##Smoother, 0); \