Feature: linear solver update: ILUp and two-level AMG smoothers

This commit is contained in:
Hrvoje Jasak 2015-10-11 12:41:19 +01:00
parent a3194d15c0
commit 0791bd9ba3
29 changed files with 200 additions and 217 deletions

View file

@ -25,7 +25,17 @@ Class
Foam::processorFvPatchField
Description
Foam::processorFvPatchField
This boundary condition enables processor communication across patches.
\heading Patch usage
Example of the boundary condition specification:
\verbatim
myPatch
{
type processor;
}
\endverbatim
SourceFiles
processorFvPatchField.C
@ -132,8 +142,7 @@ public:
}
// Destructor
//- Destructor
~processorFvPatchField();

View file

@ -910,7 +910,7 @@ void Foam::fvMeshSubset::setLargeCellSubset
// 3: face coupled and used by one cell only (so should become normal,
// non-coupled patch face)
//
// Note that this is not really nessecary - but means we can size things
// Note that this is not really necessary - but means we can size things
// correctly. Also makes handling coupled faces much easier.
labelList nCellsUsingFace(oldFaces.size(), 0);

View file

@ -711,6 +711,7 @@ matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.C
matrices/blockLduMatrix/BlockLduSmoothers/BlockLduSmoother/blockLduSmoothers.C
matrices/blockLduMatrix/BlockLduSmoothers/BlockGaussSeidelSmoother/blockGaussSeidelSmoothers.C
matrices/blockLduMatrix/BlockLduSmoothers/BlockILUSmoother/blockILUSmoothers.C
matrices/blockLduMatrix/BlockLduSmoothers/BlockILUCpSmoother/blockILUCpSmoothers.C
/* compile blockVectorNSolvers earlier to exploit parallelismn */
matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C

View file

@ -287,20 +287,20 @@ evaluate()
if
(
Pstream::defaultCommsType() == Pstream::blocking
|| Pstream::defaultCommsType() == Pstream::nonBlocking
Pstream::defaultComms() == Pstream::blocking
|| Pstream::defaultComms() == Pstream::nonBlocking
)
{
forAll(*this, patchi)
{
this->operator[](patchi).initEvaluate
(
static_cast<Pstream::commsTypes>(Pstream::defaultCommsType())
Pstream::defaultComms()
);
}
// Block for any outstanding requests
if (Pstream::defaultCommsType() == Pstream::nonBlocking)
if (Pstream::defaultComms() == Pstream::nonBlocking)
{
IPstream::waitRequests();
OPstream::waitRequests();
@ -310,11 +310,11 @@ evaluate()
{
this->operator[](patchi).evaluate
(
static_cast<Pstream::commsTypes>(Pstream::defaultCommsType())
Pstream::defaultComms()
);
}
}
else if (Pstream::defaultCommsType() == Pstream::scheduled)
else if (Pstream::defaultComms() == Pstream::scheduled)
{
const lduSchedule& patchSchedule =
bmesh_.mesh().globalData().patchSchedule();
@ -356,8 +356,8 @@ evaluateCoupled()
if
(
Pstream::defaultCommsType() == Pstream::blocking
|| Pstream::defaultCommsType() == Pstream::nonBlocking
Pstream::defaultComms() == Pstream::blocking
|| Pstream::defaultComms() == Pstream::nonBlocking
)
{
forAll(*this, patchi)
@ -366,16 +366,13 @@ evaluateCoupled()
{
this->operator[](patchi).initEvaluate
(
static_cast<Pstream::commsTypes>
(
Pstream::defaultCommsType()
)
Pstream::defaultComms()
);
}
}
// Block for any outstanding requests
if (Pstream::defaultCommsType() == Pstream::nonBlocking)
if (Pstream::defaultComms() == Pstream::nonBlocking)
{
IPstream::waitRequests();
OPstream::waitRequests();
@ -387,15 +384,12 @@ evaluateCoupled()
{
this->operator[](patchi).evaluate
(
static_cast<Pstream::commsTypes>
(
Pstream::defaultCommsType()
)
Pstream::defaultComms()
);
}
}
}
else if (Pstream::defaultCommsType() == Pstream::scheduled)
else if (Pstream::defaultComms() == Pstream::scheduled)
{
const lduSchedule& patchSchedule =
bmesh_.mesh().globalData().patchSchedule();

View file

@ -76,7 +76,8 @@ Foam::coarseBlockAmgLevel<Type>::coarseBlockAmgLevel
BlockLduSmoother<Type>::New
(
matrixPtr_,
dict
dict,
"coarseSmoother"
)
),
Ax_()

View file

@ -71,7 +71,8 @@ Foam::fineBlockAmgLevel<Type>::fineBlockAmgLevel
BlockLduSmoother<Type>::New
(
matrix,
dict
dict,
"fineSmoother"
)
),
Ax_()

View file

@ -41,8 +41,8 @@ void Foam::BlockLduMatrix<Type>::initInterfaces
{
if
(
Pstream::defaultCommsType() == Pstream::blocking
|| Pstream::defaultCommsType() == Pstream::nonBlocking
Pstream::defaultComms() == Pstream::blocking
|| Pstream::defaultComms() == Pstream::nonBlocking
)
{
forAll (interfaces_, interfaceI)
@ -55,13 +55,13 @@ void Foam::BlockLduMatrix<Type>::initInterfaces
result,
*this,
interfaceCoeffs[interfaceI],
static_cast<const Pstream::commsTypes>(Pstream::defaultCommsType()),
Pstream::defaultComms(),
switchToLhs
);
}
}
}
else if (Pstream::defaultCommsType() == Pstream::scheduled)
else if (Pstream::defaultComms() == Pstream::scheduled)
{
const lduSchedule& patchSchedule = this->patchSchedule();
@ -107,18 +107,29 @@ void Foam::BlockLduMatrix<Type>::updateInterfaces
const bool switchToLhs
) const
{
if
if (Pstream::defaultComms() == Pstream::blocking)
{
forAll (interfaces_, interfaceI)
{
if (interfaces_.set(interfaceI))
{
interfaces_[interfaceI].updateInterfaceMatrix
(
Pstream::defaultCommsType() == Pstream::blocking
|| Pstream::defaultCommsType() == Pstream::nonBlocking
)
psi,
result,
*this,
interfaceCoeffs[interfaceI],
Pstream::defaultComms(),
switchToLhs
);
}
}
}
else if (Pstream::defaultComms() == Pstream::nonBlocking)
{
// Block until all sends/receives have been finished
if (Pstream::defaultCommsType() == Pstream::nonBlocking)
{
IPstream::waitRequests();
OPstream::waitRequests();
}
forAll (interfaces_, interfaceI)
{
@ -130,13 +141,13 @@ void Foam::BlockLduMatrix<Type>::updateInterfaces
result,
*this,
interfaceCoeffs[interfaceI],
static_cast<Pstream::commsTypes>(Pstream::defaultCommsType()),
Pstream::defaultComms(),
switchToLhs
);
}
}
}
else if (Pstream::defaultCommsType() == Pstream::scheduled)
else if (Pstream::defaultComms() == Pstream::scheduled)
{
const lduSchedule& patchSchedule = this->patchSchedule();

View file

@ -127,7 +127,8 @@ class BlockILUCpPrecon
const Field<Type>& b
) const;
//- Execute preconditioning with matrix transpose, decoupled version
//- Execute preconditioning with matrix transpose,
// decoupled version
void decoupledPreconditionT
(
Field<Type>& xT,

View file

@ -35,16 +35,17 @@ template<class Type>
Foam::autoPtr<Foam::BlockLduPrecon<Type> > Foam::BlockLduPrecon<Type>::New
(
const BlockLduMatrix<Type>& matrix,
const dictionary& dict
const dictionary& dict,
const word keyword
)
{
word preconName;
// handle primitive or dictionary entry
const entry& e = dict.lookupEntry("preconditioner", false, false);
// Handle primitive or dictionary entry
const entry& e = dict.lookupEntry(keyword, false, false);
if (e.isDict())
{
e.dict().lookup("preconditioner") >> preconName;
e.dict().lookup(keyword) >> preconName;
}
else
{
@ -99,26 +100,4 @@ Foam::autoPtr<Foam::BlockLduPrecon<Type> > Foam::BlockLduPrecon<Type>::New
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Foam::word Foam::BlockLduPrecon<Type>::getName(const dictionary& dict)
{
word name;
// handle primitive or dictionary entry
const entry& e = dict.lookupEntry("preconditioner", false, false);
if (e.isDict())
{
e.dict().lookup("preconditioner") >> name;
}
else
{
e.stream() >> name;
}
return name;
}
// ************************************************************************* //

View file

@ -69,12 +69,6 @@ protected:
const BlockLduMatrix<Type>& matrix_;
// Protected member functions
//- Find the smoother name (directly or from a sub-dictionary)
static word getName(const dictionary&);
public:
//- Runtime type information
@ -114,7 +108,8 @@ public:
static autoPtr<BlockLduPrecon<Type> > New
(
const BlockLduMatrix<Type>& matrix,
const dictionary& dict
const dictionary& dict,
const word keyword = word("preconditioner")
);

View file

@ -25,7 +25,7 @@ Class
BlockILUSmoother
Description
Gauss-Seidel smoother
ILU smoother
SourceFiles
blockILUSmoothers.C

View file

@ -40,13 +40,24 @@ template<class Type>
Foam::autoPtr<Foam::BlockLduSmoother<Type> > Foam::BlockLduSmoother<Type>::New
(
const BlockLduMatrix<Type>& matrix,
const dictionary& dict
const dictionary& dict,
const word keyword
)
{
word smootherName = getName(dict);
word smootherName;
// Not (yet?) needed:
// const dictionary& controls = e.isDict() ? e.dict() : dictionary::null;
// Handle primitive or dictionary entry
const entry& e = dict.lookupEntry(keyword, false, false);
if (e.isDict())
{
e.dict().lookup(keyword) >> smootherName;
}
else
{
e.stream() >> smootherName;
}
const dictionary& controls = e.isDict() ? e.dict() : dictionary::null;
typename dictionaryConstructorTable::iterator constructorIter =
dictionaryConstructorTablePtr_->find(smootherName);
@ -58,7 +69,8 @@ Foam::autoPtr<Foam::BlockLduSmoother<Type> > Foam::BlockLduSmoother<Type>::New
"autoPtr<BlockLduSmoother> BlockLduSmoother::New\n"
"(\n"
" const BlockLduMatrix<Type>& matrix,\n"
" const dictionary& dict\n"
" const dictionary& dict,\n"
" const word keyword\n"
")",
dict
) << "Unknown matrix smoother " << smootherName
@ -73,32 +85,10 @@ Foam::autoPtr<Foam::BlockLduSmoother<Type> > Foam::BlockLduSmoother<Type>::New
constructorIter()
(
matrix,
dict
controls
)
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Foam::word Foam::BlockLduSmoother<Type>::getName(const dictionary& dict)
{
word name;
// handle primitive or dictionary entry
const entry& e = dict.lookupEntry("smoother", false, false);
if (e.isDict())
{
e.dict().lookup("smoother") >> name;
}
else
{
e.stream() >> name;
}
return name;
}
// ************************************************************************* //

View file

@ -66,12 +66,6 @@ protected:
const BlockLduMatrix<Type>& matrix_;
// Protected member functions
//- Find the smoother name (directly or from a sub-dictionary)
static word getName(const dictionary&);
public:
//- Runtime type information
@ -111,7 +105,8 @@ public:
static autoPtr<BlockLduSmoother<Type> > New
(
const BlockLduMatrix<Type>& matrix,
const dictionary& dict
const dictionary& dict,
const word keyword = word("smoother")
);

View file

@ -36,6 +36,7 @@ License
#include "blockLduSmoothers.H"
#include "blockGaussSeidelSmoothers.H"
#include "BlockILUSmoother.H"
#include "BlockILUCpSmoother.H"
#include "blockLduSolvers.H"
#include "BlockDiagonalSolver.H"
@ -98,6 +99,9 @@ makeBlockSmoother(block##Type##Smoother, block##Type##GaussSeidelSmoother); \
typedef BlockILUSmoother<type > block##Type##ILUSmoother; \
makeBlockSmoother(block##Type##Smoother, block##Type##ILUSmoother); \
\
typedef BlockILUCpSmoother<type > block##Type##ILUCpSmoother; \
makeBlockSmoother(block##Type##Smoother, block##Type##ILUCpSmoother); \
\
\
/* Solvers */ \
typedef BlockLduSolver<type > block##Type##Solver; \

View file

@ -128,7 +128,7 @@ void Foam::processorLduInterface::compressedSend
const UList<Type>& f
) const
{
if (sizeof(scalar) != sizeof(float) && Pstream::floatTransfer && f.size())
if (sizeof(scalar) != sizeof(float) && f.size())
{
static const label nCmpts = sizeof(Type)/sizeof(scalar);
label nm1 = (f.size() - 1)*nCmpts;
@ -198,7 +198,7 @@ void Foam::processorLduInterface::compressedReceive
UList<Type>& f
) const
{
if (sizeof(scalar) != sizeof(float) && Pstream::floatTransfer && f.size())
if (sizeof(scalar) != sizeof(float) && f.size())
{
static const label nCmpts = sizeof(Type)/sizeof(scalar);
label nm1 = (f.size() - 1)*nCmpts;

View file

@ -59,6 +59,10 @@ class lduInterfaceField
//- Reference to the coupled patch this field is defined for
const lduInterface& coupledInterface_;
//- Update index used so that updateInterfaceMatrix is called only once
// during the construction of the matrix
bool updatedMatrix_;
// Private Member Functions
@ -80,12 +84,12 @@ public:
//- Construct given coupled patch
lduInterfaceField(const lduInterface& patch)
:
coupledInterface_(patch)
coupledInterface_(patch),
updatedMatrix_(false)
{}
// Destructor
//- Destructor
virtual ~lduInterfaceField();
@ -108,6 +112,24 @@ public:
// Coupled interface matrix update
//- Whether matrix has been updated
bool updatedMatrix() const
{
return updatedMatrix_;
}
//- Whether matrix has been updated
bool& updatedMatrix()
{
return updatedMatrix_;
}
//- Is all data available
virtual bool ready() const
{
return true;
}
//- Transform neighbour field
virtual void transformCoupleField
(

View file

@ -165,6 +165,12 @@ public:
return solverName_;
}
//- Return access to solver name
word& solverName()
{
return solverName_;
}
//- Return initial residual
scalar initialResidual() const
{
@ -485,9 +491,6 @@ public:
public:
//- Find the smoother name (directly or from a sub-dictionary)
static word getName(const dictionary&);
//- Runtime type information
virtual const word& type() const = 0;
@ -553,7 +556,8 @@ public:
const FieldField<Field, scalar>& coupleBouCoeffs,
const FieldField<Field, scalar>& coupleIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const dictionary& dict
const dictionary& dict,
const word keyword = word("smoother")
);
@ -621,9 +625,6 @@ public:
public:
//- Find the preconditioner name (directly or from a sub-dictionary)
static word getName(const dictionary&);
//- Runtime type information
virtual const word& type() const = 0;
@ -693,7 +694,8 @@ public:
const FieldField<Field, scalar>& coupleBouCoeffs,
const FieldField<Field, scalar>& coupleIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const dictionary& dict
const dictionary& dict,
const word keyword = word("preconditioner")
);

View file

@ -36,28 +36,6 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::word Foam::lduMatrix::preconditioner::getName
(
const dictionary& dict
)
{
word name;
// handle primitive or dictionary entry
const entry& e = dict.lookupEntry("preconditioner", false, false);
if (e.isDict())
{
e.dict().lookup("preconditioner") >> name;
}
else
{
e.stream() >> name;
}
return name;
}
Foam::autoPtr<Foam::lduPreconditioner>
Foam::lduPreconditioner::New
(
@ -65,16 +43,17 @@ Foam::lduPreconditioner::New
const FieldField<Field, scalar>& coupleBouCoeffs,
const FieldField<Field, scalar>& coupleIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const dictionary& dict
const dictionary& dict,
const word keyword
)
{
word preconName;
// handle primitive or dictionary entry
const entry& e = dict.lookupEntry("preconditioner", false, false);
const entry& e = dict.lookupEntry(keyword, false, false);
if (e.isDict())
{
e.dict().lookup("preconditioner") >> preconName;
e.dict().lookup(keyword) >> preconName;
}
else
{
@ -98,7 +77,8 @@ Foam::lduPreconditioner::New
" const FieldField<Field, scalar>& coupleBouCoeffs,\n"
" const FieldField<Field, scalar>& coupleIntCoeffs,\n"
" const lduInterfaceFieldPtrsList& interfaces,\n"
" const dictionary& dict\n"
" const dictionary& dict,\n"
" const word keyword\n"
")",
dict
) << "Unknown symmetric matrix preconditioner "

View file

@ -36,44 +36,23 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Foam::word Foam::lduMatrix::smoother::getName
(
const dictionary& dict
)
{
word name;
// handle primitive or dictionary entry
const entry& e = dict.lookupEntry("smoother", false, false);
if (e.isDict())
{
e.dict().lookup("smoother") >> name;
}
else
{
e.stream() >> name;
}
return name;
}
Foam::autoPtr<Foam::lduMatrix::smoother> Foam::lduMatrix::smoother::New
(
const lduMatrix& matrix,
const FieldField<Field, scalar>& coupleBouCoeffs,
const FieldField<Field, scalar>& coupleIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const dictionary& dict
const dictionary& dict,
const word keyword
)
{
word smootherName;
// Handle primitive or dictionary entry
const entry& e = dict.lookupEntry("smoother", false, false);
const entry& e = dict.lookupEntry(keyword, false, false);
if (e.isDict())
{
e.dict().lookup("smoother") >> smootherName;
e.dict().lookup(keyword) >> smootherName;
}
else
{

View file

@ -39,8 +39,8 @@ void Foam::lduMatrix::initMatrixInterfaces
{
if
(
Pstream::defaultCommsType() == Pstream::blocking
|| Pstream::defaultCommsType() == Pstream::nonBlocking
Pstream::defaultComms() == Pstream::blocking
|| Pstream::defaultComms() == Pstream::nonBlocking
)
{
forAll (interfaces, interfaceI)
@ -54,16 +54,13 @@ void Foam::lduMatrix::initMatrixInterfaces
*this,
coupleCoeffs[interfaceI],
cmpt,
static_cast<Pstream::commsTypes>
(
Pstream::defaultCommsType()
),
Pstream::defaultComms(),
switchToLhs
);
}
}
}
else if (Pstream::defaultCommsType() == Pstream::scheduled)
else if (Pstream::defaultComms() == Pstream::scheduled)
{
const lduSchedule& patchSchedule = this->patchSchedule();
@ -111,18 +108,32 @@ void Foam::lduMatrix::updateMatrixInterfaces
const bool switchToLhs
) const
{
if
if (Pstream::defaultComms() == Pstream::blocking)
{
forAll (interfaces, interfaceI)
{
if (interfaces.set(interfaceI))
{
interfaces[interfaceI].updateInterfaceMatrix
(
Pstream::defaultCommsType() == Pstream::blocking
|| Pstream::defaultCommsType() == Pstream::nonBlocking
)
psiif,
result,
*this,
coupleCoeffs[interfaceI],
cmpt,
Pstream::defaultComms(),
switchToLhs
);
}
}
}
else if (Pstream::defaultComms() == Pstream::nonBlocking)
{
// FOAM-3.1 implementation
// Block until all sends/receives have been finished
if (Pstream::defaultCommsType() == Pstream::nonBlocking)
{
IPstream::waitRequests();
OPstream::waitRequests();
}
forAll (interfaces, interfaceI)
{
@ -135,16 +146,13 @@ void Foam::lduMatrix::updateMatrixInterfaces
*this,
coupleCoeffs[interfaceI],
cmpt,
static_cast<Pstream::commsTypes>
(
Pstream::defaultCommsType()
),
Pstream::defaultComms(),
switchToLhs
);
}
}
}
else if (Pstream::defaultCommsType() == Pstream::scheduled)
else if (Pstream::defaultComms() == Pstream::scheduled)
{
const lduSchedule& patchSchedule = this->patchSchedule();

View file

@ -70,12 +70,7 @@ Foam::lduSolverPerformance Foam::PBiCG::solve
) const
{
// --- Setup class containing solver performance data
lduSolverPerformance solverPerf
(
lduMatrix::preconditioner::getName(dict()) + typeName,
fieldName()
);
lduSolverPerformance solverPerf(typeName, fieldName());
register label nCells = x.size();
@ -134,6 +129,9 @@ Foam::lduSolverPerformance Foam::PBiCG::solve
dict()
);
// Rename the solver pefformance to include precon name
solverPerf.solverName() = preconPtr->type() + typeName;
// Solver iteration
do
{

View file

@ -70,11 +70,7 @@ Foam::lduSolverPerformance Foam::PCG::solve
) const
{
// --- Setup class containing solver performance data
lduSolverPerformance solverPerf
(
lduMatrix::preconditioner::getName(dict()) + typeName,
fieldName()
);
lduSolverPerformance solverPerf(typeName, fieldName());
register label nCells = x.size();
@ -124,6 +120,9 @@ Foam::lduSolverPerformance Foam::PCG::solve
dict()
);
// Rename the solver pefformance to include precon name
solverPerf.solverName() = preconPtr->type() + typeName;
// Solver iteration
do
{

View file

@ -89,7 +89,11 @@ Foam::lduSolverPerformance Foam::bicgSolver::solve
) const
{
// Prepare solver performance
lduSolverPerformance solverPerf(typeName, fieldName());
lduSolverPerformance solverPerf
(
preconPtr_->type() + typeName,
fieldName()
);
scalarField wA(x.size());
scalarField rA(x.size());

View file

@ -89,7 +89,11 @@ Foam::lduSolverPerformance Foam::cgSolver::solve
) const
{
// Prepare solver performance
lduSolverPerformance solverPerf(typeName, fieldName());
lduSolverPerformance solverPerf
(
preconPtr_->type() + typeName,
fieldName()
);
scalarField wA(x.size());
scalarField rA(x.size());

View file

@ -125,7 +125,10 @@ Foam::lduSolverPerformance Foam::gmresSolver::solve
) const
{
// Prepare solver performance
lduSolverPerformance solverPerf(typeName, fieldName());
lduSolverPerformance solverPerf
(
preconPtr_->type() + typeName, fieldName()
);
scalarField wA(x.size());
scalarField rA(x.size());

View file

@ -3,7 +3,7 @@
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 3.2 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{

View file

@ -42,7 +42,8 @@ solvers
minCoarseEqns 4;
nMaxLevels 100;
scale on;
smoother ILU;
coarseSmoother ILU;
fineSmoother ILU;
minIter 0;
maxIter 100;

View file

@ -42,7 +42,8 @@ solvers
// minCoarseEqns 4;
// nMaxLevels 100;
// scale on;
// smoother ILU;
// fineSmoother ILU;
// coarseSmoother k ILU;
// minIter 0;
// maxIter 100;

View file

@ -42,7 +42,8 @@ solvers
// minCoarseEqns 4;
// nMaxLevels 10;
// scale on;
// smoother ILU;
// fineSmoother ILU;
// coarseSmoother k ILU;
// minIter 0;
// maxIter 100;