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

Feature ILUCp and BlockILUCp preconditioning, extended addr update and clean-up

Browse source
This commit is contained in:
Hrvoje Jasak 2015-10-28 10:08:07 +00:00
parent b1a012087d
commit 86fc338046
30 changed files with 632 additions and 254 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
src/foam/fields/Fields/VectorNFields/expandContract/ExpandTensorNField.C
View file

@ -77,7 +77,11 @@ UNARY_FUNCTION(sphericalTensorType, cmptType, expandScalar) \
\
UNARY_FUNCTION(tensorType, vectorType, expandLinear) \
UNARY_FUNCTION(diagTensorType, vectorType, expandLinear) \
UNARY_FUNCTION(sphericalTensorType, vectorType, expandLinear)
UNARY_FUNCTION(sphericalTensorType, vectorType, expandLinear) \
\
UNARY_FUNCTION(vectorType, tensorType, sumToDiag) \
UNARY_FUNCTION(vectorType, tensorType, sumMagToDiag)
namespace Foam
{

5
src/foam/fields/Fields/VectorNFields/expandContract/ExpandTensorNField.H
View file

@ -65,7 +65,10 @@ UNARY_FUNCTION(sphericalTensorType, cmptType, expandScalar) \
\
UNARY_FUNCTION(tensorType, vectorType, expandLinear) \
UNARY_FUNCTION(diagTensorType, vectorType, expandLinear) \
UNARY_FUNCTION(sphericalTensorType, vectorType, expandLinear)
UNARY_FUNCTION(sphericalTensorType, vectorType, expandLinear) \
\
UNARY_FUNCTION(vectorType, tensorType, sumToDiag) \
UNARY_FUNCTION(vectorType, tensorType, sumMagToDiag)
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //

19
src/foam/fields/expandContract/expandTensorField.C
View file

@ -85,6 +85,25 @@ void expandLinear(Field<tensor>& res, const UList<vector>& f)
}
}
void sumToDiag(Field<vector>& res, const UList<tensor>& f)
{
forAll (res, i)
{
sumToDiag(res[i], f[i]);
}
}
void sumMagToDiag(Field<vector>& res, const UList<tensor>& f)
{
forAll (res, i)
{
sumMagToDiag(res[i], f[i]);
}
}
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

3
src/foam/fields/expandContract/expandTensorField.H
View file

@ -54,6 +54,9 @@ void expandScalar(Field<vector>& res, const UList<scalar>& f);
void expandScalar(Field<tensor>& res, const UList<scalar>& f);
void expandLinear(Field<tensor>& res, const UList<vector>& f);
void sumToDiag(Field<vector>& res, const UList<tensor>& f);
void sumMagToDiag(Field<vector>& res, const UList<tensor>& f);
} // End namespace Foam
#endif

9
src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C
View file

@ -76,12 +76,15 @@ Foam::coarseBlockAmgLevel<Type>::coarseBlockAmgLevel
BlockLduSmoother<Type>::New
(
matrixPtr_,
dict,
"coarseSmoother"
dict
// ,"coarseSmoother"
)
),
Ax_()
{}
{
Info<< "Coarse AMG level check" << endl;
matrixPtr_->check();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //

9
src/foam/matrices/blockLduMatrix/BlockAmg/fineBlockAmgLevel.C
View file

@ -71,12 +71,15 @@ Foam::fineBlockAmgLevel<Type>::fineBlockAmgLevel
BlockLduSmoother<Type>::New
(
matrix,
dict,
"fineSmoother"
dict
// ,"fineSmoother"
)
),
Ax_()
{}
{
Info<< "Fine AMG level check" << endl;
matrix_.check();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //

5
src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.H
View file

@ -227,9 +227,8 @@ public:
BlockLduMatrix(BlockLduMatrix<Type>&, bool reUse);
// Destructor
virtual ~BlockLduMatrix();
//- Destructor
virtual ~BlockLduMatrix();
// Member functions

207
src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixOperations.C
View file

@ -307,9 +307,14 @@ void Foam::BlockLduMatrix<Type>::check() const
typedef typename TypeCoeffField::linearTypeField linearTypeField;
typedef typename TypeCoeffField::squareTypeField squareTypeField;
// Copy the diagonal
TypeCoeffField DiagCopy(this->diag().size());
// TypeCoeffField DiagCopy = this->diag();
// TODO: Account for coupled boundary patch coeffs
// HJ, 29/Oct/2015
// Get original diag
const TypeCoeffField& d = this->diag();
// Copy the diagonal. It is initialised to zero
TypeCoeffField SumOffDiag(d.size());
const unallocLabelList& l = lduAddr().lowerAddr();
const unallocLabelList& u = lduAddr().upperAddr();
@ -323,67 +328,110 @@ void Foam::BlockLduMatrix<Type>::check() const
if
(
Upper.activeType() == blockCoeffBase::SQUARE
|| DiagCopy.activeType() == blockCoeffBase::SQUARE
|| d.activeType() == blockCoeffBase::SQUARE
)
{
// Note: For a square coefficient, the matrix needs to be analysed
// row-by-row, including the contribution of the off-diagonal
// elements in the diagonal matrix
// Get result as linear
linearTypeField& activeSumOffDiag = SumOffDiag.asLinear();
// Do diagonal elements
// Create the diagonal matrix element without the diagonal
const squareTypeField& activeDiag = d.asSquare();
linearTypeField diagDiag(activeDiag.size());
squareTypeField luDiag(activeDiag.size());
// Take out the diagonal from the diag as a linear type
contractLinear(diagDiag, activeDiag);
// Expand diagonal only to full square type and store into luDiag
expandLinear(luDiag, diagDiag);
// Keep only off-diagonal in luDiag
luDiag = activeDiag - luDiag;
sumMagToDiag(activeSumOffDiag, luDiag);
// Do the off-diagonal elements by collapsing each row
// into the linear form
const squareTypeField& activeUpper = Upper.asSquare();
squareTypeField& activeDiagCopy = DiagCopy.asSquare();
for (register label coeffI = 0; coeffI < l.size(); coeffI++)
{
activeDiagCopy[l[coeffI]] += activeUpper[coeffI].T();
activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
activeSumOffDiag[l[coeffI]] +=
sumMagToDiag(activeUpper[coeffI].T());
activeSumOffDiag[u[coeffI]] +=
sumMagToDiag(activeUpper[coeffI]);
}
Info<< "void BlockLduMatrix<Type>::check() const : "
<< "Symmetric matrix: raw matrix difference: "
<< sum(mag(activeDiagCopy))
<< " scaled: "
<< sum(mag(activeDiagCopy))/sum(mag(this->diag().asSquare()))
// Check diagonal dominance
diagDiag = cmptMag(diagDiag);
// Divide diagonal with sum of off-diagonals
cmptDivide(diagDiag, diagDiag, activeSumOffDiag);
InfoIn("void BlockLduMatrix<Type>::check() const)")
<< "Symmetric matrix: " << activeDiag.size()
<< " diagonal dominance sym square: "
<< Foam::min(diagDiag)
<< endl;
}
else if
(
Upper.activeType() == blockCoeffBase::LINEAR
|| DiagCopy.activeType() == blockCoeffBase::LINEAR
|| d.activeType() == blockCoeffBase::LINEAR
)
{
const linearTypeField& activeDiag = d.asLinear();
const linearTypeField& activeUpper = Upper.asLinear();
linearTypeField& activeDiagCopy = DiagCopy.asLinear();
linearTypeField& activeSumOffDiag = SumOffDiag.asLinear();
for (register label coeffI = 0; coeffI < l.size(); coeffI++)
{
activeDiagCopy[l[coeffI]] += activeUpper[coeffI];
activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
activeSumOffDiag[l[coeffI]] += cmptMag(activeUpper[coeffI]);
activeSumOffDiag[u[coeffI]] += cmptMag(activeUpper[coeffI]);
}
Info<< "void BlockLduMatrix<Type>::check() const : "
<< "Symmetric matrix: raw matrix difference: "
<< sum(mag(activeDiagCopy))
<< " scaled: "
<< sum(mag(activeDiagCopy))/sum(mag(this->diag().asLinear()))
linearTypeField diagDiag = cmptMag(activeDiag);
cmptDivide(diagDiag, diagDiag, activeSumOffDiag);
InfoIn("void BlockLduMatrix<Type>::check() const)")
<< "Symmetric matrix: " << activeDiag.size()
<< " diagonal dominance sym linear: "
<< Foam::min(diagDiag)
<< endl;
}
else if
(
Upper.activeType() == blockCoeffBase::SCALAR
|| DiagCopy.activeType() == blockCoeffBase::SCALAR
|| d.activeType() == blockCoeffBase::SCALAR
)
{
const scalarTypeField& activeDiag = d.asScalar();
const scalarTypeField& activeUpper = Upper.asScalar();
scalarTypeField& activeDiagCopy = DiagCopy.asScalar();
scalarTypeField& activeSumOffDiag = SumOffDiag.asScalar();
for (register label coeffI = 0; coeffI < l.size(); coeffI++)
{
activeDiagCopy[l[coeffI]] += activeUpper[coeffI];
activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
activeSumOffDiag[l[coeffI]] += cmptMag(activeUpper[coeffI]);
activeSumOffDiag[u[coeffI]] += cmptMag(activeUpper[coeffI]);
}
Info<< "void BlockLduMatrix<Type>::check() const : "
<< "Symmetric matrix: raw matrix difference: "
<< sum(mag(activeDiagCopy))
<< " scaled: "
<< sum(mag(activeDiagCopy))/sum(mag(this->diag().asScalar()))
InfoIn("void BlockLduMatrix<Type>::check() const)")
<< "Symmetric matrix: " << activeDiag.size()
<< " diagonal dominance sym scalar: "
<< Foam::min(mag(activeDiag)/activeSumOffDiag)
<< endl;
}
}
@ -398,78 +446,121 @@ void Foam::BlockLduMatrix<Type>::check() const
(
Lower.activeType() == blockCoeffBase::SQUARE
|| Upper.activeType() == blockCoeffBase::SQUARE
|| DiagCopy.activeType() == blockCoeffBase::SQUARE
|| d.activeType() == blockCoeffBase::SQUARE
)
{
// Note: For a square coefficient, the matrix needs to be analysed
// row-by-row, including the contribution of the off-diagonal
// elements in the diagonal matrix
// Get result as linear
linearTypeField& activeSumOffDiag = SumOffDiag.asLinear();
// Do diagonal elements
// Create the diagonal matrix element without the diagonal
const squareTypeField& activeDiag = d.asSquare();
linearTypeField diagDiag(activeDiag.size());
squareTypeField luDiag(activeDiag.size());
// Take out the diagonal from the diag as a linear type
contractLinear(diagDiag, activeDiag);
// Expand diagonal only to full square type and store into luDiag
expandLinear(luDiag, diagDiag);
// Keep only off-diagonal in luDiag
luDiag = activeDiag - luDiag;
sumMagToDiag(activeSumOffDiag, luDiag);
// Do the off-diagonal elements by collapsing each row
// into the linear form
const squareTypeField& activeLower = Lower.asSquare();
const squareTypeField& activeUpper = Upper.asSquare();
squareTypeField& activeDiagCopy = DiagCopy.asSquare();
for (register label coeffI = 0; coeffI < l.size(); coeffI++)
{
activeDiagCopy[l[coeffI]] += activeLower[coeffI];
activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
activeSumOffDiag[l[coeffI]] +=
sumMagToDiag(activeLower[coeffI]);
activeSumOffDiag[u[coeffI]] +=
sumMagToDiag(activeUpper[coeffI]);
}
Info<< "void BlockLduMatrix<Type>::check() const : "
<< "Asymmetric matrix: raw matrix difference: "
<< sum(mag(activeDiagCopy))
<< " scaled: "
<< sum(mag(activeDiagCopy))/sum(mag(this->diag().asSquare()))
// Check diagonal dominance
diagDiag = cmptMag(diagDiag);
// Divide diagonal with sum of off-diagonals
cmptDivide(diagDiag, diagDiag, activeSumOffDiag);
InfoIn("void BlockLduMatrix<Type>::check() const)")
<< "Asymmetric matrix: " << activeDiag.size()
<< " diagonal dominance assym square: "
<< Foam::min(diagDiag)
<< endl;
}
else if
(
Lower.activeType() == blockCoeffBase::LINEAR
|| Upper.activeType() == blockCoeffBase::LINEAR
|| DiagCopy.activeType() == blockCoeffBase::LINEAR
|| d.activeType() == blockCoeffBase::LINEAR
)
{
const linearTypeField& activeDiag = d.asLinear();
const linearTypeField& activeLower = Lower.asLinear();
const linearTypeField& activeUpper = Upper.asLinear();
linearTypeField& activeDiagCopy = DiagCopy.asLinear();
linearTypeField& activeSumOffDiag = SumOffDiag.asLinear();
for (register label coeffI = 0; coeffI < l.size(); coeffI++)
{
activeDiagCopy[l[coeffI]] += activeLower[coeffI];
activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
activeSumOffDiag[l[coeffI]] += cmptMag(activeLower[coeffI]);
activeSumOffDiag[u[coeffI]] += cmptMag(activeUpper[coeffI]);
}
Info<< "void BlockLduMatrix<Type>::check() const : "
<< "Asymmetric matrix: raw matrix difference: "
<< sum(mag(activeDiagCopy))
<< " scaled: "
<< sum(mag(activeDiagCopy))/sum(mag(this->diag().asLinear()))
linearTypeField diagDiag = cmptMag(activeDiag);
cmptDivide(diagDiag, diagDiag, activeSumOffDiag);
InfoIn("void BlockLduMatrix<Type>::check() const)")
<< "Asymmetric matrix: " << activeDiag.size()
<< " diagonal dominance assym linear: "
<< Foam::min(diagDiag)
<< endl;
}
else if
(
Lower.activeType() == blockCoeffBase::SCALAR
|| Upper.activeType() == blockCoeffBase::SCALAR
|| DiagCopy.activeType() == blockCoeffBase::SCALAR
|| d.activeType() == blockCoeffBase::SCALAR
)
{
const scalarTypeField& activeDiag = d.asScalar();
const scalarTypeField& activeLower = Lower.asScalar();
const scalarTypeField& activeUpper = Upper.asScalar();
scalarTypeField& activeDiagCopy = DiagCopy.asScalar();
scalarTypeField& activeSumOffDiag = SumOffDiag.asScalar();
for (register label coeffI = 0; coeffI < l.size(); coeffI++)
{
activeDiagCopy[l[coeffI]] += activeLower[coeffI];
activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
activeSumOffDiag[l[coeffI]] += cmptMag(activeLower[coeffI]);
activeSumOffDiag[u[coeffI]] += cmptMag(activeUpper[coeffI]);
}
Info<< "void BlockLduMatrix<Type>::check() const : "
<< "Asymmetric matrix: raw matrix difference: "
<< sum(mag(activeDiagCopy))
<< " scaled: "
<< sum(mag(activeDiagCopy))/sum(mag(this->diag().asScalar()))
InfoIn("void BlockLduMatrix<Type>::check() const)")
<< "Asymmetric matrix: " << activeDiag.size()
<< " diagonal dominance assym scalar: "
<< Foam::min(mag(activeDiag)/activeSumOffDiag)
<< endl;
}
}
else
{
Info<< "void BlockLduMatrix<Type>::check() const : "
InfoIn("void BlockLduMatrix<Type>::check() const)")
<< "Diagonal matrix" << endl;
}
}

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

@ -35,8 +35,13 @@ void Foam::extendedBlockLduMatrix<Type>::mapOffDiagCoeffs
{
if (blockLdum.diagonal())
{
WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)")
<< "Attempted to create extended lower/upper coeffs for block "
WarningIn
(
"void extendedBlockLduMatrix<Type>::mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<Type>& blockLdum\n"
")"
) << "Attempted to create extended lower/upper coeffs for block "
<< "matrix that is diagonal."
<< nl << endl;
}
@ -108,7 +113,10 @@ void Foam::extendedBlockLduMatrix<Type>::mapOffDiagCoeffs
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
"void extendedBlockLduMatrix<Type>::mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<Type>& blockLdum\n"
")"
) << "Problem between ordinary block matrix and extended"
<< " block matrix upper coeffs type morphing."
<< abort(FatalError);
@ -195,7 +203,10 @@ void Foam::extendedBlockLduMatrix<Type>::mapOffDiagCoeffs
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
"void extendedBlockLduMatrix<Type>::mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<Type>& blockLdum\n"
")"
) << "Problem between ordinary block matrix and extended"
<< " block matrix upper/lower coeffs type morphing."
<< abort(FatalError);
@ -210,27 +221,24 @@ template<class Type>
Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
(
const BlockLduMatrix<Type>& blockLdum,
const label extensionLevel,
const polyMesh& polyMesh
const extendedLduAddressing& extLduAddr
)
:
basicBlockLduMatrix_(blockLdum),
extLduAddr_
(
extendedLduAddressing::New
(
polyMesh,
blockLdum.lduAddr(),
extensionLevel
)
),
extLduAddr_(extLduAddr),
extendedLowerPtr_(NULL),
extendedUpperPtr_(NULL)
{
if (debug)
{
Info<< "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&) :"
"Constructing extendedBlockLduMatrix."
InfoIn
(
"extendedBlockLduMatrix<Type>::extendedBlockLduMatrix\n"
"(\n"
" const BlockLduMatrix<Type>& blockLdum,\n"
" const extendedLduAddressing& extLduAddr\n"
")"
) << "Constructing extendedBlockLduMatrix."
<< endl;
}

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

@ -39,9 +39,8 @@ SourceFiles
#define extendedBlockLduMatrix_H
#include "extendedLduAddressing.H"
#include "polyMesh.H"
#include "className.H"
#include "BlockLduMatrix.H"
#include "className.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -103,8 +102,7 @@ public:
extendedBlockLduMatrix
(
const BlockLduMatrix<Type>&,
const label,
const polyMesh&
const extendedLduAddressing&
);

27
src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.C
View file

@ -41,8 +41,13 @@ void Foam::extendedBlockLduMatrix<Foam::sphericalTensor>::mapOffDiagCoeffs
{
if (blockLdum.diagonal())
{
WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)")
<< "Attempted to create extended lower/upper coeffs for block "
WarningIn
(
"void extendedBlockLduMatrix<sphericalTensor>::mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<sphericalTensor>& blockLdum\n"
")"
) << "Attempted to create extended lower/upper coeffs for block "
<< "matrix that is diagonal."
<< nl << endl;
}
@ -102,7 +107,11 @@ void Foam::extendedBlockLduMatrix<Foam::sphericalTensor>::mapOffDiagCoeffs
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
"void extendedBlockLduMatrix<sphericalTensor>::"
"mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<sphericalTensor>& blockLdum\n"
")"
) << "Problem between ordinary block matrix and extended"
<< " block matrix upper coeffs type morphing."
<< abort(FatalError);
@ -158,7 +167,11 @@ void Foam::extendedBlockLduMatrix<Foam::sphericalTensor>::mapOffDiagCoeffs
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
"void extendedBlockLduMatrix<sphericalTensor>::"
"mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<sphericalTensor>& blockLdum\n"
")"
) << "Problem between ordinary block matrix and extended"
<< " block matrix lower coeffs type morphing."
<< abort(FatalError);
@ -229,7 +242,11 @@ void Foam::extendedBlockLduMatrix<Foam::sphericalTensor>::mapOffDiagCoeffs
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
"void extendedBlockLduMatrix<sphericalTensor>::"
"mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<sphericalTensor>& blockLdum\n"
")"
) << "Problem between ordinary block matrix and extended"
<< " block matrix upper/lower coeffs type morphing."
<< abort(FatalError);

32
src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.C
View file

@ -41,8 +41,13 @@ void Foam::extendedBlockLduMatrix<Foam::symmTensor>::mapOffDiagCoeffs
{
if (blockLdum.diagonal())
{
WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)")
<< "Attempted to create extended lower/upper coeffs for block "
WarningIn
(
"void extendedBlockLduMatrix<symmTensor>::mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<symmTensor>& blockLdum\n"
")"
) << "Attempted to create extended lower/upper coeffs for block "
<< "matrix that is diagonal."
<< nl << endl;
}
@ -102,7 +107,11 @@ void Foam::extendedBlockLduMatrix<Foam::symmTensor>::mapOffDiagCoeffs
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
"void extendedBlockLduMatrix<symmTensor>::"
"mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<symmTensor>& blockLdum\n"
")"
) << "Problem between ordinary block matrix and extended"
<< " block matrix upper coeffs type morphing."
<< abort(FatalError);
@ -158,7 +167,11 @@ void Foam::extendedBlockLduMatrix<Foam::symmTensor>::mapOffDiagCoeffs
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
"void extendedBlockLduMatrix<symmTensor>::"
"mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<symmTensor>& blockLdum\n"
")"
) << "Problem between ordinary block matrix and extended"
<< " block matrix lower coeffs type morphing."
<< abort(FatalError);
@ -188,7 +201,8 @@ void Foam::extendedBlockLduMatrix<Foam::symmTensor>::mapOffDiagCoeffs
if (upper.activeType() == blockCoeffBase::SCALAR)
{
// Helper type definition
typedef typename CoeffField<symmTensor>::scalarTypeField activeType;
typedef typename CoeffField<symmTensor>::scalarTypeField
activeType;
// Get references to fields
const activeType& activeUpper = upper.asScalar();
@ -207,7 +221,8 @@ void Foam::extendedBlockLduMatrix<Foam::symmTensor>::mapOffDiagCoeffs
else if (upper.activeType() == blockCoeffBase::LINEAR)
{
// Helper type definition
typedef typename CoeffField<symmTensor>::linearTypeField activeType;
typedef typename CoeffField<symmTensor>::linearTypeField
activeType;
// Get references to fields
const activeType& activeUpper = upper.asLinear();
@ -227,7 +242,10 @@ void Foam::extendedBlockLduMatrix<Foam::symmTensor>::mapOffDiagCoeffs
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
"void extendedBlockLduMatrix<symmTensor>::mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<symmTensor>& blockLdum\n"
")"
) << "Problem between ordinary block matrix and extended"
<< " block matrix upper/lower coeffs type morphing."
<< abort(FatalError);

36
src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.C
View file

@ -41,8 +41,13 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
{
if (blockLdum.diagonal())
{
WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)")
<< "Attempted to create extended lower/upper coeffs for block "
WarningIn
(
"void extendedBlockLduMatrix<tensor>::mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<tensor>& blockLdum\n"
")"
) << "Attempted to create extended lower/upper coeffs for block "
<< "matrix that is diagonal."
<< nl << endl;
}
@ -67,7 +72,8 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
if (upper.activeType() == blockCoeffBase::SCALAR)
{
// Helper type definition
typedef typename CoeffField<tensor>::scalarTypeField activeType;
typedef typename CoeffField<tensor>::scalarTypeField
activeType;
// Get references to fields
const activeType& activeUpper = upper.asScalar();
@ -83,7 +89,8 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
else if (upper.activeType() == blockCoeffBase::LINEAR)
{
// Helper type definition
typedef typename CoeffField<tensor>::linearTypeField activeType;
typedef typename CoeffField<tensor>::linearTypeField
activeType;
// Get references to fields
const activeType& activeUpper = upper.asLinear();
@ -100,7 +107,10 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
"void extendedBlockLduMatrix<tensor>::mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<tensor>& blockLdum\n"
")"
) << "Problem between ordinary block matrix and extended"
<< " block matrix upper coeffs type morphing."
<< abort(FatalError);
@ -121,7 +131,8 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
if (lower.activeType() == blockCoeffBase::SCALAR)
{
// Helper type definition
typedef typename CoeffField<tensor>::scalarTypeField activeType;
typedef typename CoeffField<tensor>::scalarTypeField
activeType;
// Get references to fields
const activeType& activeLower = lower.asScalar();
@ -137,7 +148,8 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
else if (lower.activeType() == blockCoeffBase::LINEAR)
{
// Helper type definition
typedef typename CoeffField<tensor>::linearTypeField activeType;
typedef typename CoeffField<tensor>::linearTypeField
activeType;
// Get references to fields
const activeType& activeLower = lower.asLinear();
@ -154,7 +166,10 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
"void extendedBlockLduMatrix<tensor>::mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<tensor>& blockLdum\n"
")"
) << "Problem between ordinary block matrix and extended"
<< " block matrix lower coeffs type morphing."
<< abort(FatalError);
@ -223,7 +238,10 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
{
FatalErrorIn
(
"extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
"void extendedBlockLduMatrix<tensor>::mapOffDiagCoeffs\n"
"(\n"
" const BlockLduMatrix<tensor>& blockLdum\n"
")"
) << "Problem between ordinary block matrix and extended"
<< " block matrix upper/lower coeffs type morphing."
<< abort(FatalError);

3
src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C
View file

@ -67,7 +67,6 @@ void Foam::BlockGaussSeidelPrecon<Type>::calcInvDiag()
// For square diagonal invert diagonal only and store the rest
// info LUDiag coefficient. This avoids full inverse of the
// diagonal matrix. HJ, 20/Aug/2015
Info<< "Square diag inverse" << endl;
// Get reference to active diag
const squareTypeField& activeDiag = d.asSquare();
@ -83,7 +82,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::calcInvDiag()
// Expand diagonal only to full square type and store into luDiag
expandLinear(luDiag, lf);
// Keep only off-diagonal in ldDiag.
// Keep only off-diagonal in luDiag.
// Note change of sign to avoid multiplication with -1 when moving
// to the other side. HJ, 20/Aug/2015
luDiag -= activeDiag;

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

@ -65,8 +65,10 @@ void Foam::BlockILUCpPrecon<Type>::calcActiveTypeFactorization
const extendedLduAddressing& addr = extBlockMatrix_.extendedLduAddr();
// Get upper/lower extended addressing
const label* const __restrict__ uPtr = addr.extendedUpperAddr().begin();
const label* const __restrict__ lPtr = addr.extendedLowerAddr().begin();
const label* const __restrict__ uPtr =
addr.extendedUpperAddr().begin();
const label* const __restrict__ lPtr =
addr.extendedLowerAddr().begin();
// Get extended owner start addressing
const label* const __restrict__ ownStartPtr =
@ -87,8 +89,8 @@ void Foam::BlockILUCpPrecon<Type>::calcActiveTypeFactorization
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.
// 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
@ -505,17 +507,15 @@ Foam::BlockILUCpPrecon<Type>::BlockILUCpPrecon
)
:
BlockLduPrecon<Type>(matrix),
preconDiag_(matrix.diag()),
p_(readLabel(dict.lookup("fillInLevel"))),
extBlockMatrix_
(
matrix,
p_,
matrix.mesh().thisDb().time().template lookupObject<polyMesh>
matrix.mesh().lduAddr().extendedAddr
(
polyMesh::defaultRegion
readLabel(dict.lookup("fillInLevel"))
)
)
),
preconDiag_(matrix.diag())
{
calcFactorization();
}

9
src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H
View file

@ -62,15 +62,12 @@ class BlockILUCpPrecon
{
// Private Data
//- Extended ldu matrix
mutable extendedBlockLduMatrix<Type> extBlockMatrix_;
//- Preconditioned diagonal
mutable CoeffField<Type> preconDiag_;
//- Fill in level
const label p_;
//- Extended lduMarix
mutable extendedBlockLduMatrix<Type> extBlockMatrix_;
// Private Member Functions

79
src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.C
View file

@ -25,6 +25,7 @@ License
#include "extendedLduAddressing.H"
#include "mapPolyMesh.H"
#include "demandDrivenData.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -38,21 +39,10 @@ namespace Foam
Foam::extendedLduAddressing::extendedLduAddressing
(
const polyMesh& mesh,
const lduAddressing& lduAddr,
const label extensionLevel
)
:
MeshObject<polyMesh, extendedLduAddressing>(mesh),
lowerAddr_
(
labelList::subList
(
mesh.faceOwner(),
mesh.nInternalFaces()
)
),
upperAddr_(mesh.faceNeighbour()),
lduAddr_(lduAddr),
p_(extensionLevel),
extendedLowerPtr_(NULL),
@ -146,6 +136,44 @@ void Foam::extendedLduAddressing::markNeighbours
}
void Foam::extendedLduAddressing::calcCellCells
(
labelListList& cellCellAddr
) const
{
// Algorithm taken from primiteMeshCellCells.C in calcCellCells() member
// function.
labelList ncc(lduAddr_.size(), 0);
const unallocLabelList& own = lowerAddr();
const unallocLabelList& nei = upperAddr();
forAll (nei, faceI)
{
ncc[own[faceI]]++;
ncc[nei[faceI]]++;
}
// Size the list
cellCellAddr.setSize(ncc.size());
forAll (cellCellAddr, cellI)
{
cellCellAddr[cellI].setSize(ncc[cellI]);
}
ncc = 0;
forAll (nei, faceI)
{
label ownCellI = own[faceI];
label neiCellI = nei[faceI];
cellCellAddr[ownCellI][ncc[ownCellI]++] = neiCellI;
cellCellAddr[neiCellI][ncc[neiCellI]++] = ownCellI;
}
}
void Foam::extendedLduAddressing::calcExtendedLowerUpper() const
{
if (extendedLowerPtr_ && extendedUpperPtr_)
@ -168,19 +196,16 @@ void Foam::extendedLduAddressing::calcExtendedLowerUpper() const
<< abort(FatalError);
}
// Get polyMesh reference
const polyMesh& mesh = this->mesh();
// Allocate dynamic lists for extended owner/neighbour, which are later used
// to define ordinary labelLists (extendedLower, extendedUpper)
// Allocate dynamic lists for extended owner/neighbour, which are later
// used to define ordinary labelLists (extendedLower, extendedUpper)
// Helper type definition
typedef DynamicList<label> DynamicLabelList;
// Get the number of faces
const label nFaces = upperAddr_.size();;
const label nFaces = upperAddr().size();
// Allocate extended owner/neighbour with 4*nFaces*extensionLevel as a guess
// (based on hex mesh assumption) to prevent excessive resizing
// Allocate extended owner/neighbour with 4*nFaces*extensionLevel as a
// guess (based on hex mesh assumption) to prevent excessive resizing
DynamicLabelList dExtOwn(4*nFaces*p_);
DynamicLabelList dExtNei(4*nFaces*p_);
@ -188,10 +213,12 @@ void Foam::extendedLduAddressing::calcExtendedLowerUpper() const
HashSet<label> extCellNbrs(32*p_);
// Get number of cells in the mesh
const label nCells = mesh.nCells();
const label nCells = lduAddr_.size();
// Get a list of cells neighbouring each cell
const labelListList& cellCells = mesh.cellCells();
// Get a list of cells neighbouring each cell.
// Note optimisation to avoid a copy
labelListList cellCells;
calcCellCells(cellCells);
// Loop through cells
for (label cellI = 0; cellI < nCells; ++cellI)
@ -234,8 +261,8 @@ void Foam::extendedLduAddressing::calcFaceMap() const
}
// Get reference to ordinary owner/neighbour addressing
const unallocLabelList& own = lduAddr().lowerAddr();
const unallocLabelList& nbr = lduAddr().upperAddr();
const unallocLabelList& own = lowerAddr();
const unallocLabelList& nbr = upperAddr();
// Allocate memory for faceMap
faceMapPtr_ = new labelList(own.size(), -1);
@ -289,7 +316,7 @@ void Foam::extendedLduAddressing::calcExtendedLosort() const
// Scan the extended neighbour list to find out how many times the cell
// appears as a neighbour of the face. Done this way to avoid guessing
// and resizing list
const label matrixSize = mesh().nCells();
const label matrixSize = lduAddr_.size();
labelList nNbrOfFace(matrixSize, 0);
const unallocLabelList& nbr = extendedUpperAddr();
@ -429,7 +456,7 @@ void Foam::extendedLduAddressing::calcExtendedLosortStart() const
}
// Set up last lookup by hand
lsrtStart[this->mesh().nCells()] = nbr.size();
lsrtStart[lduAddr_.size()] = nbr.size();
}

50
src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.H
View file

@ -25,15 +25,15 @@ Class
Foam::extendedLduAddressing
Description
The class contains data and calculates demand driven extended ldu addressing
of arbitrary level. Level 0 is not allowed as it is the same as ordinary
lduAddressing. Level 1 considers neighbours of neighbours for a given cell,
etc.
The class contains data and calculates demand driven extended ldu
addressing of arbitrary level. Level 0 is not allowed as it is the same
as ordinary lduAddressing. Level 1 considers neighbours of neighbours
for a given cell, etc.
The addressing is completely defined with extended variants of:
owner/neighbour, losort, ownerStart and losortStart: all being demand-driven
data. In order to define extended owner/neighbour lists, additional
"virtual" faces are assumed between cells not sharing a face.
owner/neighbour, losort, ownerStart and losortStart: all being
demand-driven data. In order to define extended owner/neighbour lists,
additional "virtual" faces are assumed between cells not sharing a face.
Author
Vuko Vukcevic, FMENA Zagreb. All rights reserved
@ -46,9 +46,8 @@ SourceFiles
#ifndef extendedLduAddressing_H
#define extendedLduAddressing_H
#include "MeshObject.H"
#include "polyMesh.H"
#include "lduAddressing.H"
#include "HashSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -62,17 +61,9 @@ class mapPolyMesh;
\*---------------------------------------------------------------------------*/
class extendedLduAddressing
:
public MeshObject<polyMesh, extendedLduAddressing>
{
// Private data
//- Lower as a subList of allOwner
labelList::subList lowerAddr_;
//- Upper as a reference to neighbour
const labelList& upperAddr_;
//- Reference to basic lduAddressing
const lduAddressing& lduAddr_;
@ -121,7 +112,7 @@ class extendedLduAddressing
// Parameters:
// masterCellI - searching for extended neighbours of this cell
// nbrCells - neighbouring cells of a cell (based on level)
// cellCells - mesh data to avoid function calls
// cellCells - passing as a parameter for minor optimisation
// extendedNbrs - HashSet containing labels of extended neighbours
// curLevel - current level needed for recursion control
void markNeighbours
@ -136,6 +127,10 @@ class extendedLduAddressing
// Calculation of demand driven data
//- Calculate cells cells using lower/upper addressing
// into provided storage
void calcCellCells(labelListList& cc) const;
//- Calculate extended lower/upper
void calcExtendedLowerUpper() const;
@ -160,10 +155,9 @@ public:
// Constructors
//- Construct given polyMesh
//- Construct given lduAddressing
extendedLduAddressing
(
const polyMesh& mesh,
const lduAddressing& lduAddr,
const label extensionLevel
);
@ -178,22 +172,22 @@ public:
// Access
//- Return basic lduAddressing
const lduAddressing& lduAddr() const
{
return lduAddr_;
}
//- Return lower addressing (i.e. lower label = upper triangle)
const unallocLabelList& lowerAddr() const
{
return lowerAddr_;
return lduAddr_.lowerAddr();
}
//- Return upper addressing (i.e. upper label)
const unallocLabelList& upperAddr() const
{
return upperAddr_;
}
//- Return basic lduAddressing
const lduAddressing& lduAddr() const
{
return lduAddr_;
return lduAddr_.upperAddr();
}
//- Return extended lower addressing

36
src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.C
View file

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "lduAddressing.H"
#include "extendedLduAddressing.H"
#include "demandDrivenData.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -169,6 +170,18 @@ void Foam::lduAddressing::calcLosortStart() const
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::lduAddressing::lduAddressing(const label nEqns)
:
size_(nEqns),
losortPtr_(NULL),
ownerStartPtr_(NULL),
losortStartPtr_(NULL),
extendedAddr_(5)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::lduAddressing::~lduAddressing()
@ -248,4 +261,27 @@ Foam::label Foam::lduAddressing::triIndex(const label a, const label b) const
}
const Foam::extendedLduAddressing&
Foam::lduAddressing::extendedAddr(const label p) const
{
if (p == 0 || p > 4)
{
FatalErrorIn
(
"const Foam::extendedLduAddressing& "
"Foam::lduAddressing::extendedAddr(const label p) const"
) << "Currently supported extended addressing fill-in only "
<< "between order 1 and 4"
<< abort(FatalError);
}
if (!extendedAddr_.set(p))
{
extendedAddr_.set(p, new extendedLduAddressing(*this, p));
}
return extendedAddr_[p];
}
// ************************************************************************* //

24
src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.H
View file

@ -104,8 +104,11 @@ SourceFiles
namespace Foam
{
// Forward declaration of classes
class extendedLduAddressing;
/*---------------------------------------------------------------------------*\
Class lduAddressing Declaration
Class lduAddressing Declaration
\*---------------------------------------------------------------------------*/
class lduAddressing
@ -127,6 +130,9 @@ class lduAddressing
//- Losort start addressing
mutable labelList* losortStartPtr_;
//- Extended addressing with p-order fill-in
mutable PtrList<extendedLduAddressing> extendedAddr_;
// Private Member Functions
@ -149,18 +155,11 @@ class lduAddressing
public:
// Constructor
lduAddressing(const label nEqns)
:
size_(nEqns),
losortPtr_(NULL),
ownerStartPtr_(NULL),
losortStartPtr_(NULL)
{}
lduAddressing(const label nEqns);
// Destructor
virtual ~lduAddressing();
//- Destructor
virtual ~lduAddressing();
// Member Functions
@ -200,6 +199,9 @@ public:
//- Return off-diagonal index given owner and neighbour label
label triIndex(const label a, const label b) const;
//- Return extended addressing given p index
const extendedLduAddressing& extendedAddr(const label p) const;
};

21
src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.C
View file

@ -38,20 +38,11 @@ namespace Foam
Foam::extendedLduMatrix::extendedLduMatrix
(
const lduMatrix& ldum,
const label extensionLevel,
const polyMesh& polyMesh
const extendedLduAddressing& extLduAddr
)
:
basicLduMatrix_(ldum),
extLduAddr_
(
extendedLduAddressing::New
(
polyMesh,
ldum.lduAddr(),
extensionLevel
)
),
extLduAddr_(extLduAddr),
extendedLowerPtr_(NULL),
extendedUpperPtr_(NULL)
{
@ -102,10 +93,10 @@ Foam::extendedLduMatrix::extendedLduMatrix
const label nExtFaces = extLduAddr_.extendedUpperAddr().size();
// Allocate extended upper and lower
extendedUpperPtr_ = new scalarField(nExtFaces, 0.0);
extendedUpperPtr_ = new scalarField(nExtFaces, scalar(0));
scalarField& extUpper = *extendedUpperPtr_;
extendedLowerPtr_ = new scalarField(nExtFaces, 0.0);
extendedLowerPtr_ = new scalarField(nExtFaces, scalar(0));
scalarField& extLower = *extendedLowerPtr_;
// Get upper and lower coeffs from underlying lduMatrix
@ -154,7 +145,7 @@ Foam::scalarField& Foam::extendedLduMatrix::extendedLower()
extendedLowerPtr_ = new scalarField
(
extLduAddr_.extendedLowerAddr().size(),
0.0
scalar(0)
);
}
}
@ -176,7 +167,7 @@ Foam::scalarField& Foam::extendedLduMatrix::extendedUpper()
extendedUpperPtr_ = new scalarField
(
extLduAddr_.extendedUpperAddr().size(),
0.0
scalar(0)
);
}
}

14
src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.H
View file

@ -39,9 +39,8 @@ SourceFiles
#define extendedLduMatrix_H
#include "extendedLduAddressing.H"
#include "polyMesh.H"
#include "className.H"
#include "lduMatrix.H"
#include "className.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -88,12 +87,15 @@ public:
// Constructors
//- Construct given lduMatrix, extension level and polyMesh
extendedLduMatrix(const lduMatrix&, const label, const polyMesh&);
extendedLduMatrix
(
const lduMatrix&,
const extendedLduAddressing&
);
// Destructor
~extendedLduMatrix();
//- Destructor
~extendedLduMatrix();
// Member functions

42
src/foam/primitives/VectorN/TensorNI.H
View file

@ -104,13 +104,14 @@ inline TensorN<Cmpt, length> TensorN<Cmpt, length>::T() const
return transpose;
}
//- Return tensor diagonal
template <class Cmpt, int length>
inline DiagTensorN<Cmpt, length> TensorN<Cmpt, length>::diag() const
{
DiagTensorN<Cmpt, length> dt;
label diagI=0;
label diagI = 0;
for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
{
dt[i] = this->v_[diagI];
@ -120,6 +121,7 @@ inline DiagTensorN<Cmpt, length> TensorN<Cmpt, length>::diag() const
return dt;
}
//- Negative sum the vertical off-diagonal components
template <class Cmpt, int length>
inline TensorN<Cmpt, length> TensorN<Cmpt, length>::negSumDiag() const
@ -153,6 +155,7 @@ inline TensorN<Cmpt, length> TensorN<Cmpt, length>::negSumDiag() const
return negsumdiag;
}
//- Assign to a SphericalTensorN
template <class Cmpt, int length>
inline void
@ -502,7 +505,11 @@ operator+(const TensorN<Cmpt, length>& t1, const TensorN<Cmpt, length>& t2)
//- Addition of TensorN and DiagTensorN
template <class Cmpt, int length>
inline TensorN<Cmpt, length>
operator+(const TensorN<Cmpt, length>& t1, const DiagTensorN<Cmpt, length>& dt2)
operator+
(
const TensorN<Cmpt, length>& t1,
const DiagTensorN<Cmpt, length>& dt2
)
{
TensorN<Cmpt, length> result(t1);
@ -520,7 +527,11 @@ operator+(const TensorN<Cmpt, length>& t1, const DiagTensorN<Cmpt, length>& dt2)
//- Addition of DiagTensorN and TensorN
template <class Cmpt, int length>
inline TensorN<Cmpt, length>
operator+(const DiagTensorN<Cmpt, length>& dt1, const TensorN<Cmpt, length>& t2)
operator+
(
const DiagTensorN<Cmpt, length>& dt1,
const TensorN<Cmpt, length>& t2
)
{
TensorN<Cmpt, length> result(t2);
@ -602,7 +613,11 @@ operator-(const TensorN<Cmpt, length>& t1, const TensorN<Cmpt, length>& t2)
//- Subtraction of TensorN and DiagTensorN
template <class Cmpt, int length>
inline TensorN<Cmpt, length>
operator-(const TensorN<Cmpt, length>& t1, const DiagTensorN<Cmpt, length>& dt2)
operator-
(
const TensorN<Cmpt, length>& t1,
const DiagTensorN<Cmpt, length>& dt2
)
{
TensorN<Cmpt, length> result(t1);
@ -620,7 +635,11 @@ operator-(const TensorN<Cmpt, length>& t1, const DiagTensorN<Cmpt, length>& dt2)
//- Subtraction of DiagTensorN and TensorN
template <class Cmpt, int length>
inline TensorN<Cmpt, length>
operator-(const DiagTensorN<Cmpt, length>& dt1, const TensorN<Cmpt, length>& t2)
operator-
(
const DiagTensorN<Cmpt, length>& dt1,
const TensorN<Cmpt, length>& t2
)
{
TensorN<Cmpt, length> result(-t2);
@ -709,7 +728,11 @@ operator/(const TensorN<Cmpt, length>& t1, const TensorN<Cmpt, length>& t2)
//- Inner Product of a DiagTensorN and an inverse TensorN
template <class Cmpt, int length>
inline TensorN<Cmpt, length>
operator/(const DiagTensorN<Cmpt, length>& dt1, const TensorN<Cmpt, length>& t2)
operator/
(
const DiagTensorN<Cmpt, length>& dt1,
const TensorN<Cmpt, length>& t2
)
{
return dt1 & inv(t2);
}
@ -718,7 +741,11 @@ operator/(const DiagTensorN<Cmpt, length>& dt1, const TensorN<Cmpt, length>& t2)
//- Inner Product of a TensorN and an inverse DiagTensorN
template <class Cmpt, int length>
inline TensorN<Cmpt, length>
operator/(const TensorN<Cmpt, length>& t1, const DiagTensorN<Cmpt, length>& dt2)
operator/
(
const TensorN<Cmpt, length>& t1,
const DiagTensorN<Cmpt, length>& dt2
)
{
return t1 & inv(dt2);
}
@ -972,6 +999,7 @@ inline DiagTensorN<Cmpt, length> diag(const TensorN<Cmpt, length>& t)
return t.diag();
}
//- Return tensor diagonal
template <class Cmpt, int length>
inline TensorN<Cmpt, length> negSumDiag(const TensorN<Cmpt, length>& t)

2
src/foam/primitives/VectorN/VectorN.H
View file

@ -70,6 +70,8 @@ public:
static const char* const typeName;
static const VectorN zero;
static const VectorN one;
static const VectorN max;
static const VectorN min;
// Constructors

6
src/foam/primitives/VectorN/VectorNI.H
View file

@ -43,6 +43,12 @@ const VectorN<Cmpt, length> VectorN<Cmpt, length>::zero(0);
template <class Cmpt, int length>
const VectorN<Cmpt, length> VectorN<Cmpt, length>::one(1);
template <class Cmpt, int length>
const VectorN<Cmpt, length> VectorN<Cmpt, length>::max(VGREAT);
template <class Cmpt, int length>
const VectorN<Cmpt, length> VectorN<Cmpt, length>::min(-VGREAT);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

3
src/foam/primitives/VectorN/expandContract/ExpandTensorN.H
View file

@ -75,6 +75,9 @@ UNARY_TEMPLATE_FUNCTION(TensorN, VectorN, expandLinear)
UNARY_TEMPLATE_FUNCTION(DiagTensorN, VectorN, expandLinear)
UNARY_TEMPLATE_FUNCTION(SphericalTensorN, VectorN, expandLinear)
UNARY_TEMPLATE_FUNCTION(VectorN, TensorN, sumToDiag)
UNARY_TEMPLATE_FUNCTION(VectorN, TensorN, sumMagToDiag)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

86
src/foam/primitives/VectorN/expandContract/ExpandTensorNI.H
View file

@ -37,7 +37,7 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//- Return the Average of a vector of as a scalar
//- Return the average of a vector of as a scalar
template <class Cmpt, int length>
inline void contractScalar(Cmpt& result, const VectorN<Cmpt, length>& t)
{
@ -48,7 +48,7 @@ inline void contractScalar(Cmpt& result, const VectorN<Cmpt, length>& t)
result += t[i];
}
//- Modified 2009/11/3 by I. Clifford
// Modified 2009/11/3 by I. Clifford
result *= 1.0/VectorN<Cmpt, length>::nComponents;
}
@ -76,7 +76,7 @@ inline void contractScalar(Cmpt& result, const TensorN<Cmpt, length>& t)
j += TensorN<Cmpt, length>::rowLength+1;
}
//- Modified 2009/11/3 by I. Clifford
// Modified 2009/11/3 by I. Clifford
result *= 1.0/TensorN<Cmpt, length>::rowLength;
}
@ -118,7 +118,11 @@ inline Cmpt contractScalar(const DiagTensorN<Cmpt, length>& t)
//- Return the diagonal of a SphericalTensorN as a scalar
template <class Cmpt, int length>
inline void contractScalar(Cmpt& result, const SphericalTensorN<Cmpt, length>& t)
inline void contractScalar
(
Cmpt& result,
const SphericalTensorN<Cmpt, length>& t
)
{
result = t[0];
}
@ -134,7 +138,7 @@ inline Cmpt contractScalar(const SphericalTensorN<Cmpt, length>& t)
}
//- Return the diagonal of a TensorN as a vector
//- Return the diagonal of a TensorN as a vectorN
template <class Cmpt, int length>
inline void contractLinear
(
@ -142,16 +146,16 @@ inline void contractLinear
const TensorN<Cmpt, length>& t
)
{
int j=0;
int j = 0;
for (int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
{
result[i] = t[j];
j += TensorN<Cmpt, length>::rowLength+1;
j += TensorN<Cmpt, length>::rowLength + 1;
}
}
//- Return the diagonal of a TensorN as a vector
//- Return the diagonal of a TensorN as a vectorN
template <class Cmpt, int length>
inline VectorN<Cmpt, length> contractLinear(const TensorN<Cmpt, length>& t)
{
@ -161,7 +165,7 @@ inline VectorN<Cmpt, length> contractLinear(const TensorN<Cmpt, length>& t)
}
//- Return the diagonal of a DiagTensorN as a vector
//- Return the diagonal of a DiagTensorN as a vectorN
template <class Cmpt, int length>
inline void contractLinear
(
@ -176,7 +180,7 @@ inline void contractLinear
}
//- Return the diagonal of a DiagTensorN as a vector
//- Return the diagonal of a DiagTensorN as a vectorN
template <class Cmpt, int length>
inline VectorN<Cmpt, length> contractLinear(const DiagTensorN<Cmpt, length>& t)
{
@ -186,7 +190,7 @@ inline VectorN<Cmpt, length> contractLinear(const DiagTensorN<Cmpt, length>& t)
}
//- Return the diagonal of a SphericalTensorN as a vector
//- Return the diagonal of a SphericalTensorN as a vectorN
template <class Cmpt, int length>
inline void contractLinear
(
@ -201,7 +205,7 @@ inline void contractLinear
}
//- Return the diagonal of a SphericalTensorN as a vector
//- Return the diagonal of a SphericalTensorN as a vectorN
template <class Cmpt, int length>
inline VectorN<Cmpt, length> contractLinear
(
@ -311,6 +315,64 @@ inline void expandLinear
}
//- Sum row elements of a TensorN as a vectorN
template <class Cmpt, int length>
inline void sumToDiag
(
VectorN<Cmpt, length>& result,
const TensorN<Cmpt, length>& t
)
{
result = VectorN<Cmpt, length>::zero;
for (int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
{
for (int j = 0; j < TensorN<Cmpt, length>::rowLength; j++)
{
result[i] += t(i, j);
}
}
}
//- Sum row elements of a TensorN as a vectorN
template <class Cmpt, int length>
inline VectorN<Cmpt, length> sumToDiag(const TensorN<Cmpt, length>& t)
{
VectorN<Cmpt, length> result;
sumToDiag(result, t);
return result;
}
//- Sum magnitudes of row elements of a TensorN as a vectorN
template <class Cmpt, int length>
inline void sumMagToDiag
(
VectorN<Cmpt, length>& result,
const TensorN<Cmpt, length>& t
)
{
for (int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
{
for (int j = 0; j < TensorN<Cmpt, length>::rowLength; j++)
{
result[i] += Foam::mag(t(i, j));
}
}
}
//- Sum magnitudes of row elements of a TensorN as a vectorN
template <class Cmpt, int length>
inline VectorN<Cmpt, length> sumMagToDiag(const TensorN<Cmpt, length>& t)
{
VectorN<Cmpt, length> result;
sumMagToDiag(result, t);
return result;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

48
src/foam/primitives/expandContract/expandTensor.H
View file

@ -110,6 +110,54 @@ inline void expandLinear(Tensor<Cmpt>& result, const Vector<Cmpt>& v)
}
//- Sum row elements of a tensor as a vector
template <class Cmpt>
inline void sumToDiag(Vector<Cmpt>& result, const Tensor<Cmpt>& t)
{
result = Vector<Cmpt>
(
t.xx() + t.xy() + t.xz(),
t.yx() + t.yy() + t.yz(),
t.zx() + t.zy() + t.zz()
);
}
//- Sum row elements of a TensorN as a vectorN
template <class Cmpt>
inline Vector<Cmpt> sumToDiag(const Tensor<Cmpt>& t)
{
Vector<Cmpt> result;
sumToDiag(result, t);
return result;
}
//- Sum row elements of a tensor as a vector
template <class Cmpt>
inline void sumMagToDiag(Vector<Cmpt>& result, const Tensor<Cmpt>& t)
{
result = Vector<Cmpt>
(
Foam::mag(t.xx()) + Foam::mag(t.xy()) + Foam::mag(t.xz()),
Foam::mag(t.yx()) + Foam::mag(t.yy()) + Foam::mag(t.yz()),
Foam::mag(t.zx()) + Foam::mag(t.zy()) + Foam::mag(t.zz())
);
}
//- Sum magnitudes of row elements of a TensorN as a vectorN
template <class Cmpt>
inline Vector<Cmpt> sumMagToDiag(const Tensor<Cmpt>& t)
{
Vector<Cmpt> result;
sumMagToDiag(result, t);
return result;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

26
src/lduSolvers/lduPrecon/ILUCp/ILUCp.C
View file

@ -50,8 +50,8 @@ namespace Foam
addILUCpPreconditionerAsymMatrixConstructorToTable_;
// Add to symmetric constructor table as well. Cholesky with fill in would
// yield the same sparseness pattern as the original matrix, hence it is not
// implemented. VV, 10/Sep/2015.
// yield the same sparseness pattern as the original matrix,
// hence it is not implemented. VV, 10/Sep/2015.
lduPreconditioner::
addsymMatrixConstructorToTable<ILUCp>
addILUCpPreconditionerSymMatrixConstructorToTable_;
@ -69,8 +69,10 @@ void Foam::ILUCp::calcFactorization()
const extendedLduAddressing& addr = extMatrix_.extendedLduAddr();
// Get upper/lower extended addressing
const label* const __restrict__ uPtr = addr.extendedUpperAddr().begin();
const label* const __restrict__ lPtr = addr.extendedLowerAddr().begin();
const label* const __restrict__ uPtr =
addr.extendedUpperAddr().begin();
const label* const __restrict__ lPtr =
addr.extendedLowerAddr().begin();
// Get extended owner start addressing
const label* const __restrict__ ownStartPtr =
@ -94,8 +96,8 @@ void Foam::ILUCp::calcFactorization()
// Get number of rows
const label nRows = preconDiag_.size();
// Define start and end face ("virtual" face when extended addressing is
// used) of this row/column.
// 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
@ -239,17 +241,15 @@ Foam::ILUCp::ILUCp
coupleIntCoeffs,
interfaces
),
preconDiag_(matrix_.diag()),
p_(readLabel(dict.lookup("fillInLevel"))),
extMatrix_
(
matrix,
p_,
matrix.mesh().thisDb().time().lookupObject
<
polyMesh
>(polyMesh::defaultRegion)
matrix.mesh().lduAddr().extendedAddr
(
readLabel(dict.lookup("fillInLevel"))
)
),
preconDiag_(matrix_.diag()),
zDiag_(0),
z_(preconDiag_.size(), 0),
w_(preconDiag_.size(), 0)

11
src/lduSolvers/lduPrecon/ILUCp/ILUCp.H
View file

@ -51,7 +51,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class ILUCp Declaration
Class ILUCp Declaration
\*---------------------------------------------------------------------------*/
class ILUCp
@ -60,15 +60,12 @@ class ILUCp
{
// Private Data
//- Extended ldu matrix
extendedLduMatrix extMatrix_;
//- Preconditioned diagonal
scalarField preconDiag_;
//- Fill in level
const label p_;
//- Extended lduMatrix
extendedLduMatrix extMatrix_;
//- Temporary working diagonal
scalar zDiag_;