LduMatrix: Added support for non-parallel cyclic patches

This commit is contained in:
Henry
2012-04-25 18:18:31 +01:00
parent 2a0346dc0f
commit 8fb6aa8811
31 changed files with 573 additions and 480 deletions

View File

@ -88,7 +88,7 @@ int main(int argc, char *argv[])
"{"
" solver PBiCG;"
" preconditioner DILU;"
" tolerance (1e-13 1e-13 1e-13);"
" tolerance (1e-5 1e-5 1);"
" relTol (0 0 0);"
"}"
)()

View File

@ -316,6 +316,7 @@ meshes/lduMesh/lduMesh.C
LduMatrix = matrices/LduMatrix
$(LduMatrix)/LduMatrix/lduMatrices.C
$(LduMatrix)/LduMatrix/LduInterfaceField/LduInterfaceFields.C
$(LduMatrix)/Smoothers/lduSmoothers.C
$(LduMatrix)/Preconditioners/lduPreconditioners.C
$(LduMatrix)/Solvers/lduSolvers.C

View File

@ -485,9 +485,35 @@ boundaryInternalField() const
template<class Type, template<class> class PatchField, class GeoMesh>
Foam::lduInterfaceFieldPtrsList
Foam::LduInterfaceFieldPtrsList<Type>
Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricBoundaryField::
interfaces() const
{
LduInterfaceFieldPtrsList<Type> interfaces(this->size());
forAll(interfaces, patchi)
{
if (isA<LduInterfaceField<Type> >(this->operator[](patchi)))
{
interfaces.set
(
patchi,
&refCast<const LduInterfaceField<Type> >
(
this->operator[](patchi)
)
);
}
}
return interfaces;
}
template<class Type, template<class> class PatchField, class GeoMesh>
Foam::lduInterfaceFieldPtrsList
Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricBoundaryField::
scalarInterfaces() const
{
lduInterfaceFieldPtrsList interfaces(this->size());
@ -498,7 +524,10 @@ interfaces() const
interfaces.set
(
patchi,
&refCast<const lduInterfaceField>(this->operator[](patchi))
&refCast<const lduInterfaceField>
(
this->operator[](patchi)
)
);
}
}

View File

@ -44,6 +44,7 @@ SourceFiles
#include "DimensionedField.H"
#include "FieldField.H"
#include "lduInterfaceFieldPtrsList.H"
#include "LduInterfaceFieldPtrsList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -191,7 +192,11 @@ public:
//- Return a list of pointers for each patch field with only those
// pointing to interfaces being set
lduInterfaceFieldPtrsList interfaces() const;
LduInterfaceFieldPtrsList<Type> interfaces() const;
//- Return a list of pointers for each patch field with only those
// pointing to interfaces being set
lduInterfaceFieldPtrsList scalarInterfaces() const;
//- Write boundary field as dictionary entry
void writeEntry(const word& keyword, Ostream& os) const;

View File

@ -1,65 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "CyclicLduInterfaceField.H"
#include "diagTensorField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
//defineTypeNameAndDebug(CyclicLduInterfaceField, 0);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type>
Foam::CyclicLduInterfaceField<Type>::~CyclicLduInterfaceField()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void Foam::CyclicLduInterfaceField<Type>::transformCoupleField
(
Field<Type>& f
) const
{
if (doTransform())
{
label sizeby2 = f.size()/2;
for (label facei=0; facei<sizeby2; facei++)
{
f[facei] = transform(f[facei], forwardT()[0]);
f[facei + sizeby2] = transform(f[facei + sizeby2], forwardT()[0]);
}
}
}
// ************************************************************************* //

View File

@ -1,102 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::CyclicLduInterfaceField
Description
Abstract base class for cyclic coupled interfaces.
SourceFiles
CyclicLduInterfaceField.C
\*---------------------------------------------------------------------------*/
#ifndef CyclicLduInterfaceField_H
#define CyclicLduInterfaceField_H
#include "primitiveFieldsFwd.H"
#include "typeInfo.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class CyclicLduInterfaceField Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class CyclicLduInterfaceField
{
public:
//- Runtime type information
TypeName("CyclicLduInterfaceField");
// Constructors
//- Construct given coupled patch
CyclicLduInterfaceField()
{}
// Destructor
virtual ~CyclicLduInterfaceField();
// Member Functions
// Access
//- Is the transform required
virtual bool doTransform() const = 0;
//- Return face transformation tensor
virtual const tensorField& forwardT() const = 0;
//- Return neighbour-cell transformation tensor
virtual const tensorField& reverseT() const = 0;
//- Return rank of component for transform
virtual int rank() const = 0;
//- Transform given patch internal field
void transformCoupleField(Field<Type>& f) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,25 +23,11 @@ License
\*---------------------------------------------------------------------------*/
#include "lduInterfaceField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(lduInterfaceField, 0);
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
lduInterfaceField::~lduInterfaceField()
template<class Type>
Foam::LduInterfaceField<Type>::~LduInterfaceField()
{}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -36,27 +36,28 @@ SourceFiles
#ifndef LduInterfaceField_H
#define LduInterfaceField_H
#include "lduInterface.H"
#include "Field.H"
#include "lduInterfaceField.H"
#include "primitiveFieldsFwd.H"
#include "Pstream.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
class lduMatrix;
/*---------------------------------------------------------------------------*\
Class LduInterfaceField Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class LduInterfaceField
:
public lduInterfaceField
{
// Private data
//- Reference to the coupled patch this field is defined for
const lduInterface& interface_;
// Private Member Functions
//- Disallow default bitwise copy construct
@ -68,73 +69,52 @@ class LduInterfaceField
public:
class Amultiplier
{
public:
Amultiplier()
{}
virtual ~Amultiplier()
{}
virtual void addAmul
(
Field<Type>& Apsi,
const Field<Type>& psi
) const = 0;
};
//- Runtime type information
TypeName("LduInterfaceField");
// Constructors
//- Construct given interface
LduInterfaceField(const lduInterface& interface)
//- Construct given coupled patch
LduInterfaceField(const lduInterface& patch)
:
interface_(interface)
lduInterfaceField(patch)
{}
// Destructor
virtual ~LduInterfaceField();
//- Destructor
virtual ~LduInterfaceField();
// Member Functions
// Access
//- Initialise neighbour matrix update
virtual void initInterfaceMatrixUpdate
(
Field<Type>&,
const Field<Type>&,
const scalarField&,
const Pstream::commsTypes commsType
) const
{}
//- Return the interface
const lduInterface& interface() const
{
return interface_;
}
// Coupled interface matrix update
//- Initialise neighbour matrix update
virtual void initInterfaceMatrixUpdate
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
Field<Type>&,
const Field<Type>&,
const scalarField&,
const Pstream::commsTypes commsType
) const // = 0;
{
notImplemented
(
Field<Type>& Apsi,
const Field<Type>& psi,
const Amultiplier&,
const Pstream::commsTypes commsType
) const
{}
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
Field<Type>& Apsi,
const Field<Type>& psi,
const Amultiplier&,
const Pstream::commsTypes commsType
) const = 0;
"updateInterfaceMatrix"
"(Field<Type>&, const Field<Type>&,"
"const scalarField&,"
"const Pstream::commsTypes commsType) const"
);
}
};
@ -144,6 +124,12 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "LduInterfaceField.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,66 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "processorLduInterfaceField.H"
#include "diagTensorField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
//defineTypeNameAndDebug(processorLduInterfaceField, 0);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type>
Foam::processorLduInterfaceField<Type>::~processorLduInterfaceField()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void Foam::processorLduInterfaceField<Type>::transformCoupleField
(
Field<Type>& f
) const
{
if (doTransform())
{
if (forwardT().size() == 1)
{
transform(f, forwardT()[0], f);
}
else
{
transform(f, forwardT(), f);
}
}
}
// ************************************************************************* //

View File

@ -1,105 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::ProcessorLduInterfaceField
Description
Abstract base class for processor coupled interfaces.
SourceFiles
ProcessorLduInterfaceField.C
\*---------------------------------------------------------------------------*/
#ifndef ProcessorLduInterfaceField_H
#define ProcessorLduInterfaceField_H
#include "primitiveFieldsFwd.H"
#include "typeInfo.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class ProcessorLduInterfaceField Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class ProcessorLduInterfaceField
{
public:
//- Runtime type information
TypeName("ProcessorLduInterfaceField");
// Constructors
//- Construct given coupled patch
ProcessorLduInterfaceField()
{}
// Destructor
virtual ~ProcessorLduInterfaceField();
// Member Functions
// Access
//- Return processor number
virtual int myProcNo() const = 0;
//- Return neigbour processor number
virtual int neighbProcNo() const = 0;
//- Is the transform required
virtual bool doTransform() const = 0;
//- Return face transformation tensor
virtual const tensorField& forwardT() const = 0;
//- Return rank of component for transform
virtual int rank() const = 0;
//- Transform given patch component field
void transformCoupleField(Field<Type>& f) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -75,6 +75,27 @@ Ostream& operator<<
const LduMatrix<Type, DType, LUType>&
);
template<class Type, class DType, class LUType>
typename LduMatrix<Type, DType, LUType>::solverPerformance max
(
const typename LduMatrix<Type, DType, LUType>::solverPerformance&,
const typename LduMatrix<Type, DType, LUType>::solverPerformance&
);
template<class Type, class DType, class LUType>
Istream& operator>>
(
Istream&,
typename LduMatrix<Type, DType, LUType>::solverPerformance&
);
template<class Type, class DType, class LUType>
Ostream& operator<<
(
Ostream&,
const typename LduMatrix<Type, DType, LUType>::solverPerformance&
);
/*---------------------------------------------------------------------------*\
Class LduMatrix Declaration
@ -162,6 +183,20 @@ public:
return solverName_;
}
//- Return solver name
word& solverName()
{
return solverName_;
}
//- Return field name
const word& fieldName() const
{
return fieldName_;
}
//- Return initial residual
const Type& initialResidual() const
{
@ -200,8 +235,12 @@ public:
return noIterations_;
}
//- Check, store and return singularity
bool singular(const Type& wApA);
//- Has the solver converged?
bool converged() const
{
return converged_;
}
//- Is the matrix singular?
bool singular() const;
@ -213,14 +252,41 @@ public:
const Type& relTolerance
);
//- Has the solver converged?
bool converged() const
{
return converged_;
}
//- Singularity test
bool checkSingularity(const Type& residual);
//- Print summary of solver performance to the given stream
void print(Ostream& os) const;
// Member Operators
bool operator!=(const solverPerformance&) const;
// Friend functions
//- Return the element-wise maximum of two solverPerformances
friend solverPerformance max <Type, DType, LUType>
(
const solverPerformance&,
const solverPerformance&
);
// Ostream Operator
friend Istream& operator>> <Type, DType, LUType>
(
Istream&,
solverPerformance&
);
friend Ostream& operator<< <Type, DType, LUType>
(
Ostream&,
const solverPerformance&
);
};
@ -635,6 +701,12 @@ public:
return interfaces_;
}
//- Return interfaces
LduInterfaceFieldPtrsList<Type>& interfaces()
{
return interfaces_;
}
// Access to coefficients

View File

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "LduMatrix.H"
#include "LduInterfaceFieldPtrsList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -28,7 +28,7 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::singular
bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::checkSingularity
(
const Type& wApA
)
@ -115,4 +115,85 @@ void Foam::LduMatrix<Type, DType, LUType>::solverPerformance::print
}
template<class Type, class DType, class LUType>
bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::operator!=
(
const LduMatrix<Type, DType, LUType>::solverPerformance& sp
) const
{
return
(
solverName() != sp.solverName()
|| fieldName() != sp.fieldName()
|| initialResidual() != sp.initialResidual()
|| finalResidual() != sp.finalResidual()
|| nIterations() != sp.nIterations()
|| converged() != sp.converged()
|| singular() != sp.singular()
);
}
template<class Type, class DType, class LUType>
typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance Foam::max
(
const typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance& sp1,
const typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance& sp2
)
{
return lduMatrix::solverPerformance
(
sp1.solverName(),
sp1.fieldName_,
max(sp1.initialResidual(), sp2.initialResidual()),
max(sp1.finalResidual(), sp2.finalResidual()),
max(sp1.nIterations(), sp2.nIterations()),
sp1.converged() && sp2.converged(),
sp1.singular() || sp2.singular()
);
}
template<class Type, class DType, class LUType>
Foam::Istream& Foam::operator>>
(
Istream& is,
typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance& sp
)
{
is.readBeginList("LduMatrix<Type, DType, LUType>::solverPerformance");
is >> sp.solverName_
>> sp.fieldName_
>> sp.initialResidual_
>> sp.finalResidual_
>> sp.noIterations_
>> sp.converged_
>> sp.singular_;
is.readEndList("LduMatrix<Type, DType, LUType>::solverPerformance");
return is;
}
template<class Type, class DType, class LUType>
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance& sp
)
{
os << token::BEGIN_LIST
<< sp.solverName_ << token::SPACE
<< sp.fieldName_ << token::SPACE
<< sp.initialResidual_ << token::SPACE
<< sp.finalResidual_ << token::SPACE
<< sp.noIterations_ << token::SPACE
<< sp.converged_ << token::SPACE
<< sp.singular_ << token::SPACE
<< token::END_LIST;
return os;
}
// ************************************************************************* //

View File

@ -24,7 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "LduMatrix.H"
#include "LduInterfaceField.H"
#include "lduInterfaceField.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -50,7 +50,8 @@ void Foam::LduMatrix<Type, DType, LUType>::initMatrixInterfaces
(
result,
psiif,
Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
interfaceCoeffs[interfaceI],
//Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
Pstream::defaultCommsType
);
}
@ -75,7 +76,8 @@ void Foam::LduMatrix<Type, DType, LUType>::initMatrixInterfaces
(
result,
psiif,
Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
interfaceCoeffs[interfaceI],
//Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
Pstream::blocking
);
}
@ -120,7 +122,8 @@ void Foam::LduMatrix<Type, DType, LUType>::updateMatrixInterfaces
(
result,
psiif,
Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
interfaceCoeffs[interfaceI],
//Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
Pstream::defaultCommsType
);
}
@ -143,7 +146,8 @@ void Foam::LduMatrix<Type, DType, LUType>::updateMatrixInterfaces
(
result,
psiif,
Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
interfaceCoeffs[interfaceI],
//Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
Pstream::scheduled
);
}
@ -153,7 +157,8 @@ void Foam::LduMatrix<Type, DType, LUType>::updateMatrixInterfaces
(
result,
psiif,
Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
interfaceCoeffs[interfaceI],
//Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
Pstream::scheduled
);
}
@ -175,7 +180,8 @@ void Foam::LduMatrix<Type, DType, LUType>::updateMatrixInterfaces
(
result,
psiif,
Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
interfaceCoeffs[interfaceI],
//Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
Pstream::blocking
);
}

View File

@ -150,13 +150,18 @@ Foam::TPBiCG<Type, DType, LUType>::solve(Field<Type>& psi) const
// --- Update preconditioned residuals
this->matrix_.Amul(wA, pA);
this->matrix_.Tmul(wT, pT);
Type wApT = gSumCmptProd(wA, pT);
// --- Test for singularity
if (solverPerf.singular(cmptDivide(cmptMag(wApT), normFactor)))
if
(
solverPerf.checkSingularity
(
cmptDivide(cmptMag(wApT), normFactor)
)
)
{
break;
}

View File

@ -126,7 +126,7 @@ Foam::TPBiCG<Type, DType, LUType>::solve
// --- Update search directions:
//wArT = gSumProd(wA, rT);
wArT = sum(wA & rT);
wArT = gSum(wA && rT);
if (solverPerf.nIterations() == 0)
{
@ -138,7 +138,7 @@ Foam::TPBiCG<Type, DType, LUType>::solve
}
else
{
scalar beta = cmptDivide(wArT, wArTold);
scalar beta = wArT/wArTold;
for (register label cell=0; cell<nCells; cell++)
{
@ -152,19 +152,24 @@ Foam::TPBiCG<Type, DType, LUType>::solve
this->matrix_.Amul(wA, pA);
this->matrix_.Tmul(wT, pT);
scalar wApT = sum(wA & pT);
scalar wApT = gSum(wA && pT);
// --- Test for singularity
//if
//(
// solverPerf.checkSingularity(ratio(mag(wApT)normFactor))
//) break;
if
(
solverPerf.checkSingularity
(
cmptDivide(Type::one*mag(wApT), normFactor)
)
)
{
break;
}
// --- Update solution and residual:
scalar alpha = cmptDivide(wArT, wApT);
scalar alpha = wArT/wApT;
for (register label cell=0; cell<nCells; cell++)
{

View File

@ -143,7 +143,13 @@ Foam::TPCG<Type, DType, LUType>::solve(Field<Type>& psi) const
// --- Test for singularity
if (solverPerf.singular(cmptDivide(cmptMag(wApA), normFactor)))
if
(
solverPerf.checkSingularity
(
cmptDivide(cmptMag(wApA), normFactor)
)
)
{
break;
}

View File

@ -24,7 +24,8 @@ License
\*---------------------------------------------------------------------------*/
#include "TPCG.H"
#include "TPBiCG.H"
//#include "TPBiCG.H"
#include "PBiCGScalarAlpha.H"
#include "SmoothSolver.H"
#include "fieldTypes.H"
@ -46,7 +47,7 @@ License
namespace Foam
{
makeLduSolvers(scalar, scalar, scalar);
//makeLduSolvers(scalar, scalar, scalar);
makeLduSolvers(vector, scalar, scalar);
makeLduSolvers(sphericalTensor, scalar, scalar);
makeLduSolvers(symmTensor, scalar, scalar);

View File

@ -84,10 +84,14 @@ public:
virtual int rank() const = 0;
//- Transform given patch internal field
//- Transform given patch field
template<class Type>
void transformCoupleField(Field<Type>& f) const;
//- Transform given patch component field
void transformCoupleField
(
scalarField& psiInternal,
scalarField& f,
const direction cmpt
) const;
};
@ -97,6 +101,30 @@ public:
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "tensorField.H"
template<class Type>
void Foam::cyclicLduInterfaceField::transformCoupleField
(
Field<Type>& f
) const
{
if (doTransform())
{
if (forwardT().size() == 1)
{
transform(f, forwardT()[0], f);
}
else
{
transform(f, forwardT(), f);
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

View File

@ -73,6 +73,7 @@ class lduInterfaceField
//- Disallow default bitwise assignment
void operator=(const lduInterfaceField&);
public:
//- Runtime type information
@ -112,17 +113,6 @@ public:
// Coupled interface matrix update
//- Initialise neighbour matrix update
virtual void initInterfaceMatrixUpdate
(
scalarField&,
const scalarField&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
) const
{}
//- Whether matrix has been updated
bool updatedMatrix() const
{
@ -141,6 +131,17 @@ public:
return true;
}
//- Initialise neighbour matrix update
virtual void initInterfaceMatrixUpdate
(
scalarField&,
const scalarField&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
) const
{}
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(

View File

@ -87,6 +87,10 @@ public:
virtual int rank() const = 0;
//- Transform given patch field
template<class Type>
void transformCoupleField(Field<Type>& f) const;
//- Transform given patch component field
void transformCoupleField
(
@ -100,6 +104,30 @@ public:
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "tensorField.H"
template<class Type>
void Foam::processorLduInterfaceField::transformCoupleField
(
Field<Type>& f
) const
{
if (doTransform())
{
if (forwardT().size() == 1)
{
transform(f, forwardT()[0], f);
}
else
{
transform(f, forwardT(), f);
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

View File

@ -210,7 +210,7 @@ public:
return singular_;
}
//- Convergence test
//- Check, store and return convergence
bool checkConvergence
(
const scalar tolerance,
@ -223,6 +223,7 @@ public:
//- Print summary of solver performance
void print() const;
// Member Operators
bool operator!=(const solverPerformance&) const;

View File

@ -167,9 +167,11 @@ Foam::lduMatrix::solverPerformance Foam::PBiCG::solve
scalar wApT = gSumProd(wA, pT);
// --- Test for singularity
if (solverPerf.checkSingularity(mag(wApT)/normFactor)) break;
if (solverPerf.checkSingularity(mag(wApT)/normFactor))
{
break;
}
// --- Update solution and residual:

View File

@ -39,7 +39,7 @@ coupledFvPatchField<Type>::coupledFvPatchField
const DimensionedField<Type, volMesh>& iF
)
:
lduInterfaceField(refCast<const lduInterface>(p)),
LduInterfaceField<Type>(refCast<const lduInterface>(p)),
fvPatchField<Type>(p, iF)
{}
@ -52,7 +52,7 @@ coupledFvPatchField<Type>::coupledFvPatchField
const Field<Type>& f
)
:
lduInterfaceField(refCast<const lduInterface>(p)),
LduInterfaceField<Type>(refCast<const lduInterface>(p)),
fvPatchField<Type>(p, iF, f)
{}
@ -66,7 +66,7 @@ coupledFvPatchField<Type>::coupledFvPatchField
const fvPatchFieldMapper& mapper
)
:
lduInterfaceField(refCast<const lduInterface>(p)),
LduInterfaceField<Type>(refCast<const lduInterface>(p)),
fvPatchField<Type>(ptf, p, iF, mapper)
{}
@ -79,7 +79,7 @@ coupledFvPatchField<Type>::coupledFvPatchField
const dictionary& dict
)
:
lduInterfaceField(refCast<const lduInterface>(p)),
LduInterfaceField<Type>(refCast<const lduInterface>(p)),
fvPatchField<Type>(p, iF, dict)
{}
@ -90,7 +90,7 @@ coupledFvPatchField<Type>::coupledFvPatchField
const coupledFvPatchField<Type>& ptf
)
:
lduInterfaceField(refCast<const lduInterface>(ptf.patch())),
LduInterfaceField<Type>(refCast<const lduInterface>(ptf.patch())),
fvPatchField<Type>(ptf)
{}
@ -102,7 +102,7 @@ coupledFvPatchField<Type>::coupledFvPatchField
const DimensionedField<Type, volMesh>& iF
)
:
lduInterfaceField(refCast<const lduInterface>(ptf.patch())),
LduInterfaceField<Type>(refCast<const lduInterface>(ptf.patch())),
fvPatchField<Type>(ptf, iF)
{}

View File

@ -35,7 +35,7 @@ SourceFiles
#ifndef coupledFvPatchField_H
#define coupledFvPatchField_H
#include "lduInterfaceField.H"
#include "LduInterfaceField.H"
#include "fvPatchField.H"
#include "coupledFvPatch.H"
@ -51,7 +51,7 @@ namespace Foam
template<class Type>
class coupledFvPatchField
:
public lduInterfaceField,
public LduInterfaceField<Type>,
public fvPatchField<Type>
{

View File

@ -221,6 +221,38 @@ void cyclicFvPatchField<Type>::updateInterfaceMatrix
}
template<class Type>
void cyclicFvPatchField<Type>::updateInterfaceMatrix
(
Field<Type>& result,
const Field<Type>& psiInternal,
const scalarField& coeffs,
const Pstream::commsTypes
) const
{
Field<Type> pnf(this->size());
const labelUList& nbrFaceCells =
cyclicPatch().cyclicPatch().neighbPatch().faceCells();
forAll(pnf, facei)
{
pnf[facei] = psiInternal[nbrFaceCells[facei]];
}
// Transform according to the transformation tensors
transformCoupleField(pnf);
// Multiply the field by coefficients and add into the result
const labelUList& faceCells = cyclicPatch_.faceCells();
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
template<class Type>
void cyclicFvPatchField<Type>::write(Ostream& os) const
{

View File

@ -166,6 +166,15 @@ public:
const Pstream::commsTypes commsType
) const;
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
Field<Type>& result,
const Field<Type>& psiInternal,
const scalarField& coeffs,
const Pstream::commsTypes commsType
) const;
// Cyclic coupled interface functions

View File

@ -45,6 +45,7 @@ processorFvPatchField<Type>::processorFvPatchField
coupledFvPatchField<Type>(p, iF),
procPatch_(refCast<const processorFvPatch>(p)),
sendBuf_(0),
receiveBuf_(0),
outstandingSendRequest_(-1),
outstandingRecvRequest_(-1),
scalarSendBuf_(0),
@ -63,6 +64,7 @@ processorFvPatchField<Type>::processorFvPatchField
coupledFvPatchField<Type>(p, iF, f),
procPatch_(refCast<const processorFvPatch>(p)),
sendBuf_(0),
receiveBuf_(0),
outstandingSendRequest_(-1),
outstandingRecvRequest_(-1),
scalarSendBuf_(0),
@ -83,6 +85,7 @@ processorFvPatchField<Type>::processorFvPatchField
coupledFvPatchField<Type>(ptf, p, iF, mapper),
procPatch_(refCast<const processorFvPatch>(p)),
sendBuf_(0),
receiveBuf_(0),
outstandingSendRequest_(-1),
outstandingRecvRequest_(-1),
scalarSendBuf_(0),
@ -126,6 +129,7 @@ processorFvPatchField<Type>::processorFvPatchField
coupledFvPatchField<Type>(p, iF, dict),
procPatch_(refCast<const processorFvPatch>(p)),
sendBuf_(0),
receiveBuf_(0),
outstandingSendRequest_(-1),
outstandingRecvRequest_(-1),
scalarSendBuf_(0),
@ -162,6 +166,7 @@ processorFvPatchField<Type>::processorFvPatchField
coupledFvPatchField<Type>(ptf),
procPatch_(refCast<const processorFvPatch>(ptf.patch())),
sendBuf_(ptf.sendBuf_.xfer()),
receiveBuf_(ptf.receiveBuf_.xfer()),
outstandingSendRequest_(-1),
outstandingRecvRequest_(-1),
scalarSendBuf_(ptf.scalarSendBuf_.xfer()),
@ -186,6 +191,7 @@ processorFvPatchField<Type>::processorFvPatchField
coupledFvPatchField<Type>(ptf, iF),
procPatch_(refCast<const processorFvPatch>(ptf.patch())),
sendBuf_(0),
receiveBuf_(0),
outstandingSendRequest_(-1),
outstandingRecvRequest_(-1),
scalarSendBuf_(0),
@ -428,6 +434,124 @@ void processorFvPatchField<Type>::updateInterfaceMatrix
}
template<class Type>
void processorFvPatchField<Type>::initInterfaceMatrixUpdate
(
Field<Type>&,
const Field<Type>& psiInternal,
const scalarField&,
const Pstream::commsTypes commsType
) const
{
this->patch().patchInternalField(psiInternal, sendBuf_);
if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
{
// Fast path.
if (debug && !this->ready())
{
FatalErrorIn
(
"processorFvPatchField<Type>::initInterfaceMatrixUpdate(..)"
) << "On patch " << procPatch_.name()
<< " outstanding request."
<< abort(FatalError);
}
receiveBuf_.setSize(sendBuf_.size());
outstandingRecvRequest_ = UPstream::nRequests();
IPstream::read
(
Pstream::nonBlocking,
procPatch_.neighbProcNo(),
reinterpret_cast<char*>(receiveBuf_.begin()),
receiveBuf_.byteSize(),
procPatch_.tag()
);
outstandingSendRequest_ = UPstream::nRequests();
OPstream::write
(
Pstream::nonBlocking,
procPatch_.neighbProcNo(),
reinterpret_cast<const char*>(sendBuf_.begin()),
sendBuf_.byteSize(),
procPatch_.tag()
);
}
else
{
procPatch_.compressedSend(commsType, sendBuf_);
}
const_cast<processorFvPatchField<Type>&>(*this).updatedMatrix() = false;
}
template<class Type>
void processorFvPatchField<Type>::updateInterfaceMatrix
(
Field<Type>& result,
const Field<Type>&,
const scalarField& coeffs,
const Pstream::commsTypes commsType
) const
{
if (this->updatedMatrix())
{
return;
}
const labelUList& faceCells = this->patch().faceCells();
if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
{
// Fast path.
if
(
outstandingRecvRequest_ >= 0
&& outstandingRecvRequest_ < Pstream::nRequests()
)
{
UPstream::waitRequest(outstandingRecvRequest_);
}
// Recv finished so assume sending finished as well.
outstandingSendRequest_ = -1;
outstandingRecvRequest_ = -1;
// Consume straight from receiveBuf_
// Transform according to the transformation tensor
transformCoupleField(receiveBuf_);
// Multiply the field by coefficients and add into the result
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*receiveBuf_[elemI];
}
}
else
{
Field<Type> pnf
(
procPatch_.compressedReceive<Type>(commsType, this->size())()
);
// Transform according to the transformation tensor
transformCoupleField(pnf);
// Multiply the field by coefficients and add into the result
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
const_cast<processorFvPatchField<Type>&>(*this).updatedMatrix() = true;
}
template<class Type>
bool processorFvPatchField<Type>::ready() const
{

View File

@ -64,6 +64,9 @@ class processorFvPatchField
//- Send buffer.
mutable Field<Type> sendBuf_;
//- Receive buffer.
mutable Field<Type> receiveBuf_;
//- Outstanding request
mutable label outstandingSendRequest_;
@ -184,6 +187,9 @@ public:
//- Return patch-normal gradient
virtual tmp<Field<Type> > snGrad() const;
//- Is all data available
virtual bool ready() const;
//- Initialise neighbour matrix update
virtual void initInterfaceMatrixUpdate
(
@ -194,9 +200,6 @@ public:
const Pstream::commsTypes commsType
) const;
//- Is all data available
virtual bool ready() const;
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
@ -207,6 +210,25 @@ public:
const Pstream::commsTypes commsType
) const;
//- Initialise neighbour matrix update
virtual void initInterfaceMatrixUpdate
(
Field<Type>& result,
const Field<Type>& psiInternal,
const scalarField& coeffs,
const Pstream::commsTypes commsType
) const;
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
Field<Type>& result,
const Field<Type>& psiInternal,
const scalarField& coeffs,
const Pstream::commsTypes commsType
) const;
//- Processor coupled interface functions
//- Return processor number

View File

@ -111,7 +111,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
);
lduInterfaceFieldPtrsList interfaces =
psi.boundaryField().interfaces();
psi.boundaryField().scalarInterfaces();
// Use the initMatrixInterfaces and updateMatrixInterfaces to correct
// bouCoeffsCmpt for the explicit part of the coupled boundary
@ -245,7 +245,7 @@ Foam::tmp<Foam::Field<Type> > Foam::fvMatrix<Type>::residual() const
psiCmpt,
res.component(cmpt) - boundaryDiagCmpt*psiCmpt,
bouCoeffsCmpt,
psi_.boundaryField().interfaces(),
psi_.boundaryField().scalarInterfaces(),
cmpt
)
);

View File

@ -80,7 +80,7 @@ Foam::fvMatrix<Foam::scalar>::solver
*this,
boundaryCoeffs_,
internalCoeffs_,
psi_.boundaryField().interfaces(),
psi_.boundaryField().scalarInterfaces(),
solverControls
)
)
@ -158,7 +158,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve
*this,
boundaryCoeffs_,
internalCoeffs_,
psi.boundaryField().interfaces(),
psi.boundaryField().scalarInterfaces(),
solverControls
)->solve(psi.internalField(), totalSource);
@ -187,7 +187,7 @@ Foam::tmp<Foam::scalarField> Foam::fvMatrix<Foam::scalar>::residual() const
psi_.internalField(),
source_ - boundaryDiag*psi_.internalField(),
boundaryCoeffs_,
psi_.boundaryField().interfaces(),
psi_.boundaryField().scalarInterfaces(),
0
)
);