Clean-up of blockLduMatrix

This commit is contained in:
Hrvoje Jasak 2010-11-08 10:54:47 +00:00
parent cd2c635939
commit bfc189abe8
30 changed files with 194 additions and 353 deletions

View file

@ -86,7 +86,7 @@ public:
// Coupled interface matrix update
//- Initialise matrix update
virtual void initMatrixUpdate
virtual void initInterfaceMatrixUpdate
(
const BlockLduMatrix<Type>& matrix,
Field<Type>& Ax,
@ -95,7 +95,7 @@ public:
{}
//- Update result based on interface functionality
virtual void updateMatrix
virtual void updateInterfaceMatrix
(
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& coeffs,

View file

@ -54,8 +54,19 @@ Foam::BlockLduMatrix<Type>::BlockLduMatrix(const lduMesh& ldu)
diagPtr_(NULL),
upperPtr_(NULL),
lowerPtr_(NULL),
interfaces_(),
coupleUpper_(ldu.lduAddr().nPatches()),
coupleLower_(ldu.lduAddr().nPatches()),
fixedEqns_(ldu.lduAddr().size()/fixFillIn)
{}
{
const lduAddressing& addr = ldu.lduAddr();
forAll (coupleUpper_, i)
{
coupleUpper_.set(i, new CoeffField<Type>(addr.patchAddr(i).size()));
coupleLower_.set(i, new CoeffField<Type>(addr.patchAddr(i).size()));
}
}
template<class Type>
@ -66,6 +77,9 @@ Foam::BlockLduMatrix<Type>::BlockLduMatrix(const BlockLduMatrix<Type>& A)
diagPtr_(NULL),
upperPtr_(NULL),
lowerPtr_(NULL),
interfaces_(),
coupleUpper_(ldu.lduAddr().nPatches()),
coupleLower_(ldu.lduAddr().nPatches()),
fixedEqns_(A.fixedEqns_)
{
if (A.diagPtr_)
@ -82,8 +96,19 @@ Foam::BlockLduMatrix<Type>::BlockLduMatrix(const BlockLduMatrix<Type>& A)
{
lowerPtr_ = new TypeCoeffField(*(A.lowerPtr_));
}
const lduAddressing& addr = ldu.lduAddr();
forAll (coupleUpper_, i)
{
coupleUpper_.set(i, A.coupleUpper_.clone());
coupleLower_.set(i, , A.coupleLower_.clone());
}
}
//HJ, problematic: memmory management.
// Reconsider. HJ, 7/Nov/2010
template<class Type>
Foam::BlockLduMatrix<Type>::BlockLduMatrix(BlockLduMatrix<Type>& A, bool reUse)
:
@ -92,6 +117,9 @@ Foam::BlockLduMatrix<Type>::BlockLduMatrix(BlockLduMatrix<Type>& A, bool reUse)
diagPtr_(NULL),
upperPtr_(NULL),
lowerPtr_(NULL),
interfaces_(),
coupleUpper_(),
coupleLower_(),
fixedEqns_(A.fixedEqns_)
{
if (reUse)
@ -131,8 +159,13 @@ Foam::BlockLduMatrix<Type>::BlockLduMatrix(BlockLduMatrix<Type>& A, bool reUse)
lowerPtr_ = new TypeCoeffField(*(A.lowerPtr_));
}
}
// Interface data
coupleUpper_ = A.coupleUpper_;
coupleLower_ = A.coupleUpper_;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type>
@ -143,6 +176,7 @@ Foam::BlockLduMatrix<Type>::~BlockLduMatrix()
deleteDemandDrivenData(lowerPtr_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //

View file

@ -51,6 +51,7 @@ SourceFiles
#include "coeffFields.H"
#include "FieldField.H"
#include "lduMesh.H"
#include "BlockLduInterface.H"
#include "BlockLduInterfaceFieldPtrsList.H"
#include "Map.H"
@ -97,7 +98,7 @@ private:
const lduMesh& lduMesh_;
// Block matrix elements (excluding interfaces)
// Block matrix elements
//- Diagonal coefficients
CoeffField<Type>* diagPtr_;
@ -109,6 +110,18 @@ private:
CoeffField<Type> *lowerPtr_;
// Coupling
//- List of coupled interfaces
PtrList<BlockLduInterface<Type> > interfaces_;
//- Coupled interface coefficients, upper
FieldField<CoeffField, Type> coupleUpper_;
//- Coupled interface coefficients, lower
FieldField<CoeffField, Type> coupleLower_;
// Constraints
//- Equation triangle map
@ -209,6 +222,7 @@ public:
//- Construct as copy or re-use as specified.
BlockLduMatrix(BlockLduMatrix<Type>&, bool reUse);
// Destructor
virtual ~BlockLduMatrix();
@ -259,6 +273,30 @@ public:
//- Return lower coefficients
const TypeCoeffField& lower() const;
//- Return access to coupled interface coefficients, upper
FieldField<CoeffField, Type>& coupleUpper()
{
return coupleUpper_;
}
//- Return coupled interface coefficients, upper
const FieldField<CoeffField, Type>& coupleUpper() const
{
return coupleUpper_;
}
//- Return access to coupled interface coefficients, lower
FieldField<CoeffField, Type>& coupleLower()
{
return coupleLower_;
}
//- Return coupled interface coefficients, lower
const FieldField<CoeffField, Type>& coupleLower() const
{
return coupleLower_;
}
// Matrix structure
@ -319,10 +357,7 @@ public:
void Amul
(
TypeField& Ax,
const TypeField& x,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces
const TypeField& x
) const;
//- Matrix multiplication without coupled interface update
@ -336,10 +371,7 @@ public:
void Tmul
(
TypeField& Ax,
const TypeField& x,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces
const TypeField& x
) const;
//- Matrix transpose multiplication without coupled
@ -366,7 +398,6 @@ public:
void initInterfaces
(
const FieldField<CoeffField, Type>& interfaceCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
TypeField& result,
const TypeField& psi
) const;
@ -375,7 +406,6 @@ public:
void updateInterfaces
(
const FieldField<CoeffField, Type>& interfaceCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
TypeField& result,
const TypeField& psi
) const;

View file

@ -35,20 +35,18 @@ template<class Type>
void Foam::BlockLduMatrix<Type>::Amul
(
TypeField& Ax,
const TypeField& x,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces
const TypeField& x
) const
{
Ax = pTraits<Type>::zero;
// Initialise the update of coupled interfaces
initInterfaces(boundaryCoeffs, interfaces, Ax, x);
initInterfaces(coupleUpper_, Ax, x);
AmulCore(Ax, x);
// Update coupled interfaces
updateInterfaces(boundaryCoeffs, interfaces, Ax, x);
updateInterfaces(coupleUpper_, Ax, x);
}
@ -179,20 +177,18 @@ template<class Type>
void Foam::BlockLduMatrix<Type>::Tmul
(
TypeField& Ax,
const TypeField& x,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces
const TypeField& x
) const
{
Ax = pTraits<Type>::zero;
// Initialise the update of coupled interfaces
initInterfaces(boundaryCoeffs, interfaces, Ax, x);
initInterfaces(coupleLower_, Ax, x);
TmulCore(Ax, x);
// Update coupled interfaces
updateInterfaces(boundaryCoeffs, interfaces, Ax, x);
updateInterfaces(coupleLower_, Ax, x);
}

View file

@ -35,7 +35,6 @@ template<class Type>
void Foam::BlockLduMatrix<Type>::initInterfaces
(
const FieldField<CoeffField, Type>& interfaceCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
TypeField& result,
const TypeField& psi
) const
@ -46,18 +45,18 @@ void Foam::BlockLduMatrix<Type>::initInterfaces
|| Pstream::defaultCommsType == Pstream::nonBlocking
)
{
forAll (interfaces, interfaceI)
forAll (interfaces_, interfaceI)
{
if (interfaces.set(interfaceI))
if (interfaces_.set(interfaceI))
{
interfaces[interfaceI].initInterfaceMatrixUpdate
(
psi,
result,
*this,
interfaceCoeffs[interfaceI],
Pstream::defaultCommsType
);
// interfaces_[interfaceI].initInterfaceMatrixUpdate
// (
// psi,
// result,
// *this,
// interfaceCoeffs[interfaceI],
// Pstream::defaultCommsType
// );
}
}
}
@ -69,24 +68,24 @@ void Foam::BlockLduMatrix<Type>::initInterfaces
// beyond the end of the schedule which only handles "normal" patches
for
(
label interfaceI=patchSchedule.size()/2;
interfaceI<interfaces.size();
label interfaceI = patchSchedule.size()/2;
interfaceI < interfaces_.size();
interfaceI++
)
{
if (interfaces.set(interfaceI))
if (interfaces_.set(interfaceI))
{
interfaces[interfaceI].initInterfaceMatrixUpdate
(
psi,
result,
*this,
interfaceCoeffs[interfaceI],
Pstream::blocking
);
// interfaces_[interfaceI].initInterfaceMatrixUpdate
// (
// psi,
// result,
// *this,
// interfaceCoeffs[interfaceI],
// Pstream::blocking
// );
}
}
}
}
else
{
FatalErrorIn("BlockLduMatrix<Type>::initMatrixInterfaces")
@ -101,7 +100,6 @@ template<class Type>
void Foam::BlockLduMatrix<Type>::updateInterfaces
(
const FieldField<CoeffField, Type>& interfaceCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
TypeField& result,
const TypeField& psi
) const
@ -119,18 +117,18 @@ void Foam::BlockLduMatrix<Type>::updateInterfaces
OPstream::waitRequests();
}
forAll (interfaces, interfaceI)
forAll (interfaces_, interfaceI)
{
if (interfaces.set(interfaceI))
if (interfaces_.set(interfaceI))
{
interfaces[interfaceI].updateInterfaceMatrix
(
psi,
result,
*this,
interfaceCoeffs[interfaceI],
Pstream::defaultCommsType
);
// interfaces_[interfaceI].updateInterfaceMatrix
// (
// psi,
// result,
// *this,
// interfaceCoeffs[interfaceI],
// Pstream::defaultCommsType
// );
}
}
}
@ -143,29 +141,29 @@ void Foam::BlockLduMatrix<Type>::updateInterfaces
{
label interfaceI = patchSchedule[i].patch;
if (interfaces.set(interfaceI))
if (interfaces_.set(interfaceI))
{
if (patchSchedule[i].init)
{
interfaces[interfaceI].initInterfaceMatrixUpdate
(
psi,
result,
*this,
interfaceCoeffs[interfaceI],
Pstream::scheduled
);
// interfaces_[interfaceI].initInterfaceMatrixUpdate
// (
// psi,
// result,
// *this,
// interfaceCoeffs[interfaceI],
// Pstream::scheduled
// );
}
else
{
interfaces[interfaceI].updateInterfaceMatrix
(
psi,
result,
*this,
interfaceCoeffs[interfaceI],
Pstream::scheduled
);
// interfaces_[interfaceI].updateInterfaceMatrix
// (
// psi,
// result,
// *this,
// interfaceCoeffs[interfaceI],
// Pstream::scheduled
// );
}
}
}
@ -174,21 +172,21 @@ void Foam::BlockLduMatrix<Type>::updateInterfaces
// beyond the end of the schedule which only handles "normal" patches
for
(
label interfaceI=patchSchedule.size()/2;
interfaceI<interfaces.size();
label interfaceI = patchSchedule.size()/2;
interfaceI < interfaces_.size();
interfaceI++
)
{
if (interfaces.set(interfaceI))
if (interfaces_.set(interfaceI))
{
interfaces[interfaceI].updateInterfaceMatrix
(
psi,
result,
*this,
interfaceCoeffs[interfaceI],
Pstream::blocking
);
// interfaces_[interfaceI].updateInterfaceMatrix
// (
// psi,
// result,
// *this,
// interfaceCoeffs[interfaceI],
// Pstream::blocking
// );
}
}
}

View file

@ -83,12 +83,10 @@ public:
BlockGaussSeidelSmoother
(
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
)
:
BlockLduSmoother<Type>(matrix, boundaryCoeffs, interfaces),
BlockLduSmoother<Type>(matrix),
gs_(matrix)
{}

View file

@ -86,12 +86,10 @@ public:
BlockILUSmoother
(
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
)
:
BlockLduSmoother<Type>(matrix, boundaryCoeffs, interfaces),
BlockLduSmoother<Type>(matrix),
precon_(matrix),
xCorr_(matrix.lduAddr().size()),
residual_(matrix.lduAddr().size())
@ -117,13 +115,7 @@ public:
for (label sweep = 0; sweep < nSweeps; sweep++)
{
// Calculate residual
this-> matrix_.Amul
(
residual_,
x,
BlockLduSmoother<Type>::boundaryCoeffs_,
BlockLduSmoother<Type>::interfaces_
);
this-> matrix_.Amul(residual_, x);
// residual = b - Ax
forAll (b, i)

View file

@ -41,8 +41,6 @@ template<class Type>
Foam::autoPtr<Foam::BlockLduSmoother<Type> > Foam::BlockLduSmoother<Type>::New
(
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
)
{
@ -72,8 +70,6 @@ Foam::autoPtr<Foam::BlockLduSmoother<Type> > Foam::BlockLduSmoother<Type>::New
"autoPtr<BlockLduSmoother> BlockLduSmoother::New\n"
"(\n"
" const BlockLduMatrix<Type>& matrix,\n"
" const FieldField<CoeffField, Type>& boundaryCoeffs,\n"
" const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,\n"
" const dictionary& dict\n"
")",
dict
@ -89,8 +85,6 @@ Foam::autoPtr<Foam::BlockLduSmoother<Type> > Foam::BlockLduSmoother<Type>::New
constructorIter()
(
matrix,
boundaryCoeffs,
interfaces,
dict
)
);

View file

@ -66,12 +66,6 @@ protected:
//- Matrix reference
const BlockLduMatrix<Type>& matrix_;
//- Reference to boundary coefficients
const FieldField<CoeffField, Type>& boundaryCoeffs_;
//- Reference to interfaces
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces_;
// Protected member functions
@ -94,14 +88,10 @@ public:
dictionary,
(
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
),
(
matrix,
boundaryCoeffs,
interfaces,
dict
)
);
@ -110,16 +100,9 @@ public:
// Constructors
//- Construct from matrix
BlockLduSmoother
(
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces
)
BlockLduSmoother(const BlockLduMatrix<Type>& matrix)
:
matrix_(matrix),
boundaryCoeffs_(boundaryCoeffs),
interfaces_(interfaces)
matrix_(matrix)
{}
@ -129,8 +112,6 @@ public:
static autoPtr<BlockLduSmoother<Type> > New
(
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
);
@ -172,8 +153,8 @@ addToRunTimeSelectionTable(SmootherType, typeSmootherType, dictionary);
#define makeBlockSmoothers(smootherType) \
\
makeBlockSmoother(blockScalarSmoother, smootherType##Scalar); \
makeBlockSmoother(blockVectorSmoother, smootherType##Vector);
// makeBlockSmoother(blockTensorSmoother, smootherType##Tensor);
makeBlockSmoother(blockVectorSmoother, smootherType##Vector); \
makeBlockSmoother(blockTensorSmoother, smootherType##Tensor);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -38,11 +38,11 @@ namespace Foam
defineNamedTemplateTypeNameAndDebug(blockScalarSmoother, 0);
defineNamedTemplateTypeNameAndDebug(blockVectorSmoother, 0);
// defineNamedTemplateTypeNameAndDebug(blockTensorSmoother, 0);
defineNamedTemplateTypeNameAndDebug(blockTensorSmoother, 0);
defineTemplateRunTimeSelectionTable(blockScalarSmoother, dictionary);
defineTemplateRunTimeSelectionTable(blockVectorSmoother, dictionary);
// defineTemplateRunTimeSelectionTable(blockTensorSmoother, dictionary);
defineTemplateRunTimeSelectionTable(blockTensorSmoother, dictionary);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -50,7 +50,7 @@ namespace Foam
typedef BlockLduSmoother<scalar> blockScalarSmoother;
typedef BlockLduSmoother<vector> blockVectorSmoother;
// typedef BlockLduSmoother<tensor> blockTensorSmoother;
typedef BlockLduSmoother<tensor> blockTensorSmoother;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -37,8 +37,6 @@ Foam::BlockBiCGStabSolver<Type>::BlockBiCGStabSolver
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
)
:
@ -46,8 +44,6 @@ Foam::BlockBiCGStabSolver<Type>::BlockBiCGStabSolver
(
fieldName,
matrix,
boundaryCoeffs,
interfaces,
dict
),
preconPtr_
@ -89,13 +85,7 @@ Foam::BlockBiCGStabSolver<Type>::solve
Field<Type> p(x.size());
// Calculate initial residual
matrix.Amul
(
p,
x,
BlockLduSolver<Type>::boundaryCoeffs_,
BlockLduSolver<Type>::interfaces_
);
matrix.Amul(p, x);
Field<Type> r(b - p);
solverPerf.initialResidual() = gSum(cmptMag(r))/norm;
@ -148,13 +138,7 @@ Foam::BlockBiCGStabSolver<Type>::solve
}
preconPtr_->precondition(ph, p);
matrix.Amul
(
v,
ph,
BlockLduSolver<Type>::boundaryCoeffs_,
BlockLduSolver<Type>::interfaces_
);
matrix.Amul(v, ph);
alpha = rho/gSumProd(rw, v);
forAll (s, i)
@ -163,13 +147,7 @@ Foam::BlockBiCGStabSolver<Type>::solve
}
preconPtr_->preconditionT(sh, s);
matrix.Amul
(
t,
sh,
BlockLduSolver<Type>::boundaryCoeffs_,
BlockLduSolver<Type>::interfaces_
);
matrix.Amul(t, sh);
omega = gSumProd(t, s)/gSumProd(t, t);
forAll (x, i)

View file

@ -85,10 +85,7 @@ public:
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces,
const dictionary& solver
const dictionary& dict
);

View file

@ -37,8 +37,6 @@ Foam::BlockCGSolver<Type>::BlockCGSolver
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
)
:
@ -46,8 +44,6 @@ Foam::BlockCGSolver<Type>::BlockCGSolver
(
fieldName,
matrix,
boundaryCoeffs,
interfaces,
dict
),
preconPtr_
@ -88,13 +84,7 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
Field<Type> wA(x.size());
// Calculate initial residual
matrix.Amul
(
wA,
x,
BlockLduSolver<Type>::boundaryCoeffs_,
BlockLduSolver<Type>::interfaces_
);
matrix.Amul(wA, x);
Field<Type> rA(b - wA);
solverPerf.initialResidual() = gSum(cmptMag(rA))/norm;
@ -129,13 +119,7 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
}
// Update preconditioner residual
matrix.Amul
(
wA,
pA,
BlockLduSolver<Type>::boundaryCoeffs_,
BlockLduSolver<Type>::interfaces_
);
matrix.Amul(wA, pA);
wApA = gSumProd(wA, pA);

View file

@ -85,9 +85,6 @@ public:
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces,
const dictionary& dict
);

View file

@ -38,19 +38,10 @@ Foam::BlockDiagonalSolver<Type>::BlockDiagonalSolver
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
)
:
BlockLduSolver<Type>
(
fieldName,
matrix,
boundaryCoeffs,
interfaces,
dict
)
BlockLduSolver<Type>(fieldName, matrix, dict)
{}

View file

@ -87,9 +87,6 @@ public:
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces,
const dictionary& dict
);

View file

@ -71,8 +71,6 @@ Foam::BlockGMRESSolver<Type>::BlockGMRESSolver
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
)
:
@ -80,8 +78,6 @@ Foam::BlockGMRESSolver<Type>::BlockGMRESSolver
(
fieldName,
matrix,
boundaryCoeffs,
interfaces,
dict
),
preconPtr_
@ -124,13 +120,7 @@ Foam::BlockGMRESSolver<Type>::solve
Field<Type> wA(x.size());
// Calculate initial residual
matrix.Amul
(
wA,
x,
BlockLduSolver<Type>::boundaryCoeffs_,
BlockLduSolver<Type>::interfaces_
);
matrix.Amul(wA, x);
Field<Type> rA(b - wA);
solverPerf.initialResidual() = gSum(cmptMag(rA))/norm;
@ -178,13 +168,7 @@ Foam::BlockGMRESSolver<Type>::solve
V[i] /= beta;
// Arnoldi's method
matrix.Amul
(
rA,
V[i],
BlockLduSolver<Type>::boundaryCoeffs_,
BlockLduSolver<Type>::interfaces_
);
matrix.Amul(rA, V[i]);
// Execute preconditioning
preconPtr_->precondition(wA, rA);
@ -247,13 +231,7 @@ Foam::BlockGMRESSolver<Type>::solve
}
// Re-calculate the residual
matrix.Amul
(
wA,
x,
BlockLduSolver<Type>::boundaryCoeffs_,
BlockLduSolver<Type>::interfaces_
);
matrix.Amul(wA, x);
forAll (rA, raI)
{

View file

@ -98,9 +98,6 @@ public:
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces,
const dictionary& dict
);

View file

@ -39,8 +39,6 @@ Foam::BlockGaussSeidelSolver<Type>::BlockGaussSeidelSolver
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
)
:
@ -48,8 +46,6 @@ Foam::BlockGaussSeidelSolver<Type>::BlockGaussSeidelSolver
(
fieldName,
matrix,
boundaryCoeffs,
interfaces,
dict
),
gs_(matrix),
@ -82,13 +78,7 @@ Foam::BlockGaussSeidelSolver<Type>::solve
Field<Type> wA(x.size());
// Calculate residual. Note: sign of residual swapped for efficiency
matrix.Amul
(
wA,
x,
BlockLduSolver<Type>::boundaryCoeffs_,
BlockLduSolver<Type>::interfaces_
);
matrix.Amul(wA, x);
wA -= b;
solverPerf.initialResidual() = gSum(cmptMag(wA))/norm;
@ -111,13 +101,7 @@ Foam::BlockGaussSeidelSolver<Type>::solve
// Re-calculate residual. Note: sign of residual swapped
// for efficiency
matrix.Amul
(
wA,
x,
BlockLduSolver<Type>::boundaryCoeffs_,
BlockLduSolver<Type>::interfaces_
);
matrix.Amul(wA, x);
wA -= b;
solverPerf.finalResidual() = gSum(cmptMag(wA))/norm;

View file

@ -90,9 +90,6 @@ public:
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces,
const dictionary& dict
);

View file

@ -37,12 +37,10 @@ Foam::BlockIterativeSolver<Type>::BlockIterativeSolver
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
)
:
BlockLduSolver<Type>(fieldName, matrix, boundaryCoeffs, interfaces, dict),
BlockLduSolver<Type>(fieldName, matrix, dict),
tolerance_(readScalar(this->dict().lookup("tolerance"))),
relTolerance_(readScalar(this->dict().lookup("relTol"))),
minIter_(readLabel(this->dict().lookup("minIter"))),
@ -71,21 +69,13 @@ Foam::scalar Foam::BlockIterativeSolver<Type>::normFactor
Type xRef = gAverage(x);
// Calculate A.x
matrix.Amul
(
wA,
x,
BlockLduSolver<Type>::boundaryCoeffs_,
BlockLduSolver<Type>::interfaces_
);
matrix.Amul(wA, x);
// Calculate A.xRef, temporarily using pA for storage
matrix.Amul
(
pA,
Field<Type>(nRows, xRef),
BlockLduSolver<Type>::boundaryCoeffs_,
BlockLduSolver<Type>::interfaces_
Field<Type>(nRows, xRef)
);
scalar normFactor = gSum(mag(wA - pA) + mag(b - pA)) + this->small_;

View file

@ -104,9 +104,6 @@ public:
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces,
const dictionary& dict
);

View file

@ -76,16 +76,12 @@ template<class Type>
Foam::BlockLduSolver<Type>::BlockLduSolver
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces
const BlockLduMatrix<Type>& matrix
)
:
fieldName_(fieldName),
dict_(),
matrix_(matrix),
boundaryCoeffs_(boundaryCoeffs),
interfaces_(interfaces)
matrix_(matrix)
{}
@ -94,16 +90,12 @@ Foam::BlockLduSolver<Type>::BlockLduSolver
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
)
:
fieldName_(fieldName),
dict_(dict),
matrix_(matrix),
boundaryCoeffs_(boundaryCoeffs),
interfaces_(interfaces)
matrix_(matrix)
{}
@ -115,8 +107,6 @@ Foam::BlockLduSolver<Type>::New
(
const word& fieldName,
BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
)
{
@ -126,14 +116,7 @@ Foam::BlockLduSolver<Type>::New
{
return autoPtr<BlockLduSolver<Type> >
(
new BlockDiagonalSolver<Type>
(
fieldName,
matrix,
boundaryCoeffs,
interfaces,
dict
)
new BlockDiagonalSolver<Type>(fieldName, matrix, dict)
);
}
else if (matrix.symmetric())
@ -146,8 +129,6 @@ Foam::BlockLduSolver<Type>::New
"(\n"
" const word& fieldName,\n"
" BlockLduMatrix<Type>& matrix,\n"
" const FieldField<CoeffField, Type>& boundaryCoeffs,\n"
" const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,\n"
" const dictionary& dict\n"
")"
) << "Initialization problem." << endl;
@ -164,8 +145,6 @@ Foam::BlockLduSolver<Type>::New
"(\n"
" const word& fieldName,\n"
" BlockLduMatrix<Type>& matrix,\n"
" const FieldField<CoeffField, Type>& boundaryCoeffs,\n"
" const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,\n"
" const dictionary& dict\n"
")",
dict
@ -182,8 +161,6 @@ Foam::BlockLduSolver<Type>::New
(
fieldName,
matrix,
boundaryCoeffs,
interfaces,
dict
)
);
@ -198,8 +175,6 @@ Foam::BlockLduSolver<Type>::New
"(\n"
" const word& fieldName,\n"
" BlockLduMatrix<Type>& matrix,\n"
" const FieldField<CoeffField, Type>& boundaryCoeffs,\n"
" const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,\n"
" const dictionary& dict\n"
")"
) << "Initialization problem." << endl;
@ -216,8 +191,6 @@ Foam::BlockLduSolver<Type>::New
"(\n"
" const word& fieldName,\n"
" BlockLduMatrix<Type>& matrix,\n"
" const FieldField<CoeffField, Type>& boundaryCoeffs,\n"
" const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,\n"
" const dictionary& dict\n"
")",
dict
@ -234,8 +207,6 @@ Foam::BlockLduSolver<Type>::New
(
fieldName,
matrix,
boundaryCoeffs,
interfaces,
dict
)
);
@ -248,9 +219,6 @@ Foam::BlockLduSolver<Type>::New
"(\n"
" const word& fieldName,\n"
" BlockLduMatrix<Type>& matrix,\n"
" const FieldField<CoeffField, Type>& boundaryCoeffs,\n"
" const typename BlockLduInterfaceFieldPtrsList<Type>::Type& "
"interfaces,\n"
" const dictionary& dict\n"
")"
) << "cannot solve incomplete matrix, "

View file

@ -67,6 +67,7 @@ class BlockLduSolver
//- Control data dictionary
dictionary dict_;
protected:
// Protected data types
@ -79,11 +80,6 @@ protected:
//- Matrix
const BlockLduMatrix<Type>& matrix_;
// Reference to boundary coefficients
const FieldField<CoeffField, Type>& boundaryCoeffs_;
//- Reference to interfaces
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces_;
// Protected Member Functions
@ -128,16 +124,11 @@ public:
(
const word& fieldName,
BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces,
const dictionary& dict
),
(
fieldName,
matrix,
boundaryCoeffs,
interfaces,
dict
)
);
@ -150,16 +141,11 @@ public:
(
const word& fieldName,
BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces,
const dictionary& dict
),
(
fieldName,
matrix,
boundaryCoeffs,
interfaces,
dict
)
);
@ -172,10 +158,7 @@ public:
BlockLduSolver
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces
const BlockLduMatrix<Type>& matrix
);
//- Construct from dictionary
@ -183,9 +166,6 @@ public:
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces,
const dictionary& dict
);
@ -197,9 +177,6 @@ public:
(
const word& fieldName,
BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces,
const dictionary& dict
);
@ -245,8 +222,8 @@ defineNamedTemplateTypeNameAndDebug(typeSolverType, 0);
#define makeBlockSolvers(solverType) \
\
makeBlockSolverTypeName(solverType##Scalar); \
makeBlockSolverTypeName(solverType##Vector);
// makeBlockSolverTypeName(solverType##Tensor);
makeBlockSolverTypeName(solverType##Vector); \
makeBlockSolverTypeName(solverType##Tensor);
#define addSolverToBlockMatrix(Type, typeSolverType, typeMatrix) \
@ -256,8 +233,8 @@ addToRunTimeSelectionTable(block##Type##Solver, typeSolverType, typeMatrix);
#define addSolversToBlockMatrix(solverType, typeMatrix) \
\
addSolverToBlockMatrix(Scalar, solverType##Scalar, typeMatrix); \
addSolverToBlockMatrix(Vector, solverType##Vector, typeMatrix);
// addSolverToBlockMatrix(Tensor, solverType##Tensor, typeMatrix);
addSolverToBlockMatrix(Vector, solverType##Vector, typeMatrix); \
addSolverToBlockMatrix(Tensor, solverType##Tensor, typeMatrix);
#define addSymSolverToBlockMatrix(solverType) \
\

View file

@ -38,7 +38,7 @@ namespace Foam
defineNamedTemplateTypeNameAndDebug(blockScalarSolver, 0);
defineNamedTemplateTypeNameAndDebug(blockVectorSolver, 0);
// defineNamedTemplateTypeNameAndDebug(blockTensorSolver, 0);
defineNamedTemplateTypeNameAndDebug(blockTensorSolver, 0);
// Define the constructor function hash tables for symmetric solvers
@ -55,11 +55,11 @@ defineTemplateRunTimeSelectionTable
symMatrix
);
// defineTemplateRunTimeSelectionTable
// (
// blockTensorSolver,
// symMatrix
// );
defineTemplateRunTimeSelectionTable
(
blockTensorSolver,
symMatrix
);
// Define the constructor function hash tables for asymmetric solvers
@ -76,11 +76,11 @@ defineTemplateRunTimeSelectionTable
asymMatrix
);
// defineTemplateRunTimeSelectionTable
// (
// blockTensorSolver,
// asymMatrix
// );
defineTemplateRunTimeSelectionTable
(
blockTensorSolver,
asymMatrix
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -44,7 +44,7 @@ namespace Foam
typedef BlockLduSolver<scalar> blockScalarSolver;
typedef BlockLduSolver<vector> blockVectorSolver;
// typedef BlockLduSolver<tensor> blockTensorSolver;
typedef BlockLduSolver<tensor> blockTensorSolver;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -38,19 +38,10 @@ Foam::SegregatedSolver<Type>::SegregatedSolver
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
const dictionary& dict
)
:
BlockLduSolver<Type>
(
fieldName,
matrix,
boundaryCoeffs,
interfaces,
dict
),
BlockLduSolver<Type>(fieldName, matrix, dict),
scalarX_(matrix.lduAddr().size()),
scalarMatrix_(matrix.mesh()),
scalarB_(matrix.lduAddr().size())
@ -181,8 +172,6 @@ Foam::BlockSolverPerformance<Type> Foam::SegregatedSolver<Type>::solve
(
this->fieldName(),
scalarMatrix_,
boundaryCoeffs,
interfaces_,
this->dict()
)->solve(scalarX_, scalarB_);

View file

@ -92,9 +92,6 @@ public:
(
const word& fieldName,
const BlockLduMatrix<Type>& matrix,
const FieldField<CoeffField, Type>& boundaryCoeffs,
const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
interfaces,
const dictionary& dict
);

View file

@ -40,10 +40,10 @@ makeBlockSolverTypeName(segregatedSolverVector);
makeBlockSolverTypeName(segregatedSolverTensor);
addSolverToBlockMatrix(Vector, segregatedSolverVector, symMatrix);
// addSolverToBlockMatrix(Tensor, segregatedSolverTensor, symMatrix);
addSolverToBlockMatrix(Tensor, segregatedSolverTensor, symMatrix);
addSolverToBlockMatrix(Vector, segregatedSolverVector, asymMatrix);
// addSolverToBlockMatrix(Tensor, segregatedSolverTensor, asymMatrix);
addSolverToBlockMatrix(Tensor, segregatedSolverTensor, asymMatrix);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //