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

New feature: generalization of storage for solverPerformanceDict based on templated BlockSolverPerformance (by Pascal Beckstein)

Browse source 
To correct the previous commit, the BlockSolverPerformance has been enhanced
by a replace(...) and max(...) memeber function.
The matrix solve functions of fvMatrix, faMatrix and tetFemMatrix now store
solver performance data for all components in the solverPerformance dictionary
of solution. However, they will still return a scalar-valued
lduSolverPerformance based on the component with the max residual.
This is necessary as otherwise the type of a non-scalar field, which may be
specified in fvSolution::residualControl does not match the type which is
written to and read from the solverPerformance dictionary in the
maxTypeResidual(...) from solutionControl. If e.g. the field U is given
in fvSolution::residualControl of solutionControl, the constructor of
BlockSolverPerformance in the maxTypeResidual(...) from solutionControl
expects to read a vector ( u1 u2 u3 ) from the dictionary and not just one
component.
The changes are inspired and directly related to the following vanilla
commit:
1944b09bb5
This commit is contained in:
Vanja Skuric 2018-02-11 18:59:15 +01:00
parent ec83275923
commit e3115846f5
13 changed files with 187 additions and 606 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

15
src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C
View file

@ -60,7 +60,7 @@ lduSolverPerformance faMatrix<Type>::solve(const dictionary& solverControls)
<< endl;
}
lduSolverPerformance solverPerfVec
BlockSolverPerformance<Type> solverPerfVec
(
"faMatrix<Type>::solve",
psi_.name()
@ -134,14 +134,7 @@ lduSolverPerformance faMatrix<Type>::solve(const dictionary& solverControls)
solverPerf.print();
if
(
solverPerf.initialResidual() > solverPerfVec.initialResidual()
&& !solverPerf.singular()
)
{
solverPerfVec = solverPerf;
}
solverPerfVec.replace(cmpt, solverPerf);
psi.internalField().replace(cmpt, psiCmpt);
diag() = saveDiag;
@ -149,7 +142,9 @@ lduSolverPerformance faMatrix<Type>::solve(const dictionary& solverControls)
psi.correctBoundaryConditions();
return solverPerfVec;
psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerfVec);
return solverPerfVec.max();
}

4
src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C
View file

@ -130,6 +130,8 @@ lduSolverPerformance faMatrix<scalar>::faSolver::solve
psi.correctBoundaryConditions();
psi.mesh().solutionDict().setSolverPerformance(psi.name(), solverPerf);
return solverPerf;
}
@ -177,6 +179,8 @@ lduSolverPerformance faMatrix<scalar>::solve
psi.correctBoundaryConditions();
psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerf);
return solverPerf;
}

556
src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
View file

@ -24,11 +24,8 @@ License
\*---------------------------------------------------------------------------*/
#include "solutionControl.H"
#include "lduMatrix.H"
#include "fvc.H"
#include "inletOutletFvPatchFields.H"
#include "slipFvPatchFields.H"
#include "symmetryFvPatchFields.H"
#include "fieldTypes.H"
#include "VectorNFieldTypes.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -180,6 +177,21 @@ void Foam::solutionControl::storePrevIterFields() const
storePrevIter<sphericalTensor>();
storePrevIter<symmTensor>();
storePrevIter<tensor>();
storePrevIter<vector2>();
storePrevIter<vector4>();
storePrevIter<vector6>();
storePrevIter<vector8>();
storePrevIter<sphericalTensor2>();
storePrevIter<sphericalTensor4>();
storePrevIter<sphericalTensor6>();
storePrevIter<sphericalTensor8>();
storePrevIter<tensor2>();
storePrevIter<tensor4>();
storePrevIter<tensor6>();
storePrevIter<tensor8>();
}
@ -196,7 +208,7 @@ void Foam::solutionControl::maxTypeResidual
if (mesh_.foundObject<fieldType>(fieldName))
{
const List<lduSolverPerformance> sp(data);
const List<BlockSolverPerformance<Type> > sp(data);
firstRes = cmptMax(sp.first().initialResidual());
lastRes = cmptMax(sp.last().initialResidual());
}
@ -218,119 +230,25 @@ Foam::scalar Foam::solutionControl::maxResidual
maxTypeResidual<symmTensor>(fieldName, data, firstRes, lastRes);
maxTypeResidual<tensor>(fieldName, data, firstRes, lastRes);
maxTypeResidual<vector2>(fieldName, data, firstRes, lastRes);
maxTypeResidual<vector4>(fieldName, data, firstRes, lastRes);
maxTypeResidual<vector6>(fieldName, data, firstRes, lastRes);
maxTypeResidual<vector8>(fieldName, data, firstRes, lastRes);
maxTypeResidual<sphericalTensor2>(fieldName, data, firstRes, lastRes);
maxTypeResidual<sphericalTensor4>(fieldName, data, firstRes, lastRes);
maxTypeResidual<sphericalTensor6>(fieldName, data, firstRes, lastRes);
maxTypeResidual<sphericalTensor8>(fieldName, data, firstRes, lastRes);
maxTypeResidual<tensor2>(fieldName, data, firstRes, lastRes);
maxTypeResidual<tensor4>(fieldName, data, firstRes, lastRes);
maxTypeResidual<tensor6>(fieldName, data, firstRes, lastRes);
maxTypeResidual<tensor8>(fieldName, data, firstRes, lastRes);
return firstRes;
}
const Foam::dimensionedScalar Foam::solutionControl::relaxFactor
(
const Foam::volVectorField& U
) const
{
scalar urf = 1;
if (mesh_.solutionDict().relaxEquation(U.name()))
{
urf = mesh_.solutionDict().equationRelaxationFactor(U.name());
}
return dimensionedScalar("alphaU", dimless, urf);
}
void Foam::solutionControl::addDdtFluxContribution
(
surfaceScalarField& phi,
surfaceScalarField& aCoeff,
const surfaceVectorField& faceU,
const volVectorField& U,
const surfaceScalarField& rAUf,
const fvVectorMatrix& ddtUEqn
) const
{
// Add ddt scheme dependent flux contribution. Here: phi = H/A & Sf.
phi += fvc::ddtConsistentPhiCorr(faceU, U, rAUf);
// Add ddt scheme dependent contribution to diffusion scale coeff. Here:
// aCoeff = 1.
aCoeff += fvc::interpolate(ddtUEqn.A())*rAUf;
}
void Foam::solutionControl::addUnderRelaxationFluxContribution
(
surfaceScalarField& phi,
surfaceScalarField& aCoeff,
const volVectorField& U
) const
{
// Get under-relaxation factor used in this iteration
const dimensionedScalar alphaU(this->relaxFactor(U));
// Return if the under-relaxation factor is >= 1 and =< 0
if ((alphaU.value() > 1.0 - SMALL) || (alphaU.value() < SMALL))
{
return;
}
const dimensionedScalar urfCoeff((1.0 - alphaU)/alphaU);
// Add under-relaxation dependent flux contribution
phi += urfCoeff*aCoeff*phi.prevIter();
// Add under-relaxation dependent contribution to diffusion scale coeff
aCoeff += urfCoeff*aCoeff;
}
void Foam::solutionControl::correctBoundaryFlux
(
surfaceScalarField& phi,
const volVectorField& U
) const
{
// Get necessary data at the boundaries
const volVectorField::GeometricBoundaryField& Ub = U.boundaryField();
const surfaceVectorField::GeometricBoundaryField& Sb =
mesh_.Sf().boundaryField();
surfaceScalarField::GeometricBoundaryField& phib = phi.boundaryField();
// Correct flux at boundaries depending on patch type. Note: it is possible
// that we missed some important boundary conditions that need special
// considerations when calculating the flux. Maybe we can reorganise this in
// future by relying on distinction between = and == operators for
// fvsPatchFields. However, that would require serious changes at the
// moment (e.g. using phi = rUA*UEqn.H() as in most solvers would not update
// the boundary values correctly for fixed fvsPatchFields). VV, 20/Apr/2017.
forAll (phib, patchI)
{
// Get patch field
const fvPatchVectorField& Up = Ub[patchI];
if
(
Up.fixesValue()
&& !isA<inletOutletFvPatchVectorField>(Up)
)
{
// This is fixed value patch, flux needs to be recalculated
// with respect to the boundary condition
phib[patchI] == (Up & Sb[patchI]);
}
else if
(
isA<slipFvPatchVectorField>(Up)
|| isA<symmetryFvPatchVectorField>(Up)
)
{
// This is slip or symmetry, flux needs to be zero
phib[patchI] == 0.0;
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::solutionControl::solutionControl(fvMesh& mesh, const word& algorithmName)
@ -349,10 +267,7 @@ Foam::solutionControl::solutionControl(fvMesh& mesh, const word& algorithmName)
transonic_(false),
consistent_(false),
corr_(0),
corrNonOrtho_(0),
aCoeffPtrs_(),
faceUPtrs_(),
indices_()
corrNonOrtho_(0)
{}
@ -362,407 +277,4 @@ Foam::solutionControl::~solutionControl()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::solutionControl::calcTransientConsistentFlux
(
surfaceScalarField& phi,
const volVectorField& U,
const volScalarField& rAU,
const fvVectorMatrix& ddtUEqn
) const
{
// Store necessary data for this velocity field
const word& UName = U.name();
// Check whether the fields are present in the list
if (!indices_.found(UName))
{
// Get current index as size of the indices list (before insertion)
const label i = indices_.size();
// Insert the index into the hash table
indices_.insert(UName, i);
// Extend lists
aCoeffPtrs_.resize(indices_.size());
faceUPtrs_.resize(indices_.size());
// Double check whether the fields have been already set
if (!aCoeffPtrs_.set(i) && !faceUPtrs_.set(i))
{
aCoeffPtrs_.set
(
i,
new surfaceScalarField
(
IOobject
(
"aCoeff." + UName,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("zero", dimless, 0.0)
)
);
faceUPtrs_.set
(
i,
new surfaceVectorField
(
IOobject
(
"faceU." + UName,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedVector("zero", dimVelocity, vector::zero)
)
);
}
else if (!aCoeffPtrs_.set(i) || !faceUPtrs_.set(i))
{
FatalErrorIn
(
"void solutionControl::calcTransientConsistentFlux"
"\n("
"\n surfaceScalarField& phi,"
"\n const volVectorField& U,"
"\n const volScalarField& rAU,"
"\n const fvVectorMatrix& ddtUEqn"
"\n)"
) << "Either aCoeff or faceU is allocated for field " << UName
<< " while the other is not." << nl
<< " This must not happen in transient simulation. Make sure"
<< " that functions aiding consistency are called in the right"
<< " order (first flux and then velocity reconstruction)."
<< exit(FatalError);
}
}
else
{
// Index has been set for this field, so the fields must be there as
// well. Check and report an error if they are not allocated
if (!aCoeffPtrs_.set(indices_[UName]))
{
FatalErrorIn
(
"void solutionControl::calcTransientConsistentFlux"
"\n("
"\n surfaceScalarField& phi,"
"\n const volVectorField& U,"
"\n const volScalarField& rAU,"
"\n const fvVectorMatrix& ddtUEqn"
"\n)"
) << "Index is set, but the aCoeff field is not allocated for "
<< UName << "." << nl
<< "This should not happen for transient simulation." << nl
<< "Something went wrong."
<< exit(FatalError);
}
else if (!faceUPtrs_.set(indices_[UName]))
{
FatalErrorIn
(
"void solutionControl::calcTransientConsistentFlux"
"\n("
"\n surfaceScalarField& phi,"
"\n const volVectorField& U,"
"\n const volScalarField& rAU,"
"\n const fvVectorMatrix& ddtUEqn"
"\n)"
) << "Index is set, but the faceU field is not allocated for "
<< UName << "." << nl
<< "This should not happen for transient simulation." << nl
<< "Something went wrong."
<< exit(FatalError);
}
}
// Algorithm:
// 1. Update flux and aCoeff due to ddt discretisation
// 2. Update flux and aCoeff due to under-relaxation
// 3. Scale the flux with aCoeff, making sure that flux at fixed boundaries
// remains consistent
// Get index from the hash table
const label i = indices_[UName];
// Get fields that will be updated
surfaceScalarField& aCoeff = aCoeffPtrs_[i];
surfaceVectorField& faceU = faceUPtrs_[i];
// Update face interpolated velocity field. Note: handling of oldTime faceU
// fields happens in ddt scheme when calling ddtConsistentPhiCorr inside
// addDdtFluxContribution
faceU = fvc::interpolate(U);
// Interpolate original rAU on the faces
const surfaceScalarField rAUf = fvc::interpolate(rAU);
// Store previous iteration for the correct handling of under-relaxation
phi.storePrevIter();
// Calculate the ordinary part of the flux (H/A)
phi = (faceU & mesh_.Sf());
// Initialize aCoeff to 1
aCoeff = dimensionedScalar("one", dimless, 1.0);
// STAGE 1: consistent ddt discretisation handling
addDdtFluxContribution(phi, aCoeff, faceU, U, rAUf, ddtUEqn);
// STAGE 2: consistent under-relaxation handling
addUnderRelaxationFluxContribution(phi, aCoeff, U);
// STAGE 3: scale the flux and correct it at the boundaries
phi /= aCoeff;
correctBoundaryFlux(phi, U);
}
void Foam::solutionControl::calcSteadyConsistentFlux
(
surfaceScalarField& phi,
const volVectorField& U
) const
{
// Store necessary data for this velocity field
const word& UName = U.name();
// Check whether the fields are present in the list
if (!indices_.found(UName))
{
// Get current index as size of the indices list (before insertion)
const label i = indices_.size();
// Insert the index into the hash table
indices_.insert(UName, i);
// Extend list (only aCoeff list because faceU does not need to be
// stored for steady state simulations)
aCoeffPtrs_.resize(indices_.size());
// Double check whether the aCoeff field has been already set. Note:
// faceU does not need to be stored for steady state
if (!aCoeffPtrs_.set(i) && faceUPtrs_.empty())
{
aCoeffPtrs_.set
(
i,
new surfaceScalarField
(
IOobject
(
"aCoeff." + UName,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("zero", dimless, 0.0)
)
);
}
else if (!aCoeffPtrs_.set(i) || !faceUPtrs_.empty())
{
FatalErrorIn
(
"void solutionControl::calcSteadyConsistentFlux"
"\n("
"\n surfaceScalarField& phi,"
"\n const volVectorField& U,"
"\n const volScalarField& rAU"
"\n)"
) << "Either aCoeff or faceU is allocated for field " << UName
<< " while the other is not."
<< "This must not happen in transient simulation. Make sure"
<< " that functions aiding consistency are called in the right"
<< " order (first flux and then velocity reconstruction)."
<< exit(FatalError);
}
}
else
{
// Index has been set for this field. Check the state of fields
if (!aCoeffPtrs_.set(indices_[UName]))
{
// aCoeff field should be allocated
FatalErrorIn
(
"void solutionControl::calcSteadyConsistentFlux"
"\n("
"\n surfaceScalarField& phi,"
"\n const volVectorField& U,"
"\n const volScalarField& rAU"
"\n)"
) << "Index is set, but the aCoeff field is not allocated for "
<< UName << "." << nl
<< "This should not happen for steady state simulation." << nl
<< "Something went wrong."
<< exit(FatalError);
}
else if (!faceUPtrs_.empty())
{
// faceU field shouldn't be allocated
FatalErrorIn
(
"void solutionControl::calcSteadyConsistentFlux"
"\n("
"\n surfaceScalarField& phi,"
"\n const volVectorField& U,"
"\n const volScalarField& rAU"
"\n)"
) << "Index is set, but the faceU field is allocated for "
<< UName << "." << nl
<< "This should not happen for steady state simulation." << nl
<< "Something went wrong."
<< exit(FatalError);
}
}
// Algorithm:
// 1. Update flux and aCoeff due to under-relaxation
// 2. Scale the flux with aCoeff, making sure that flux at fixed boundaries
// remains consistent
// Get index from the hash table
const label i = indices_[UName];
// Get aCoeff field. Note: no need to check whether the entry has been found
// since we have just inserted it above
surfaceScalarField& aCoeff = aCoeffPtrs_[i];
// Store previous iteration for the correct handling of under-relaxation
phi.storePrevIter();
// Calculate the ordinary part of the flux (H/A)
phi = (fvc::interpolate(U) & mesh_.Sf());
// Initialize aCoeff to 1
aCoeff = dimensionedScalar("one", dimless, 1.0);
// STAGE 1: consistent under-relaxation handling
addUnderRelaxationFluxContribution(phi, aCoeff, U);
// STAGE 2: scale the flux and correct it at the boundaries
phi /= aCoeff;
correctBoundaryFlux(phi, U);
}
void Foam::solutionControl::reconstructTransientVelocity
(
volVectorField& U,
surfaceScalarField& phi,
const fvVectorMatrix& ddtUEqn,
const volScalarField& rAU,
const volScalarField& p
) const
{
// Reconstruct the velocity using all the components from original equation
U = 1.0/(1.0/rAU + ddtUEqn.A())*
(
U/rAU + ddtUEqn.H() - fvc::grad(p)
);
// Get name and the corresponding index
const word& UName = U.name();
const label i = indices_[UName];
// Update divergence free face velocity field, whose value will be used in
// the next time step. Note that we need to have absolute flux here.
if (!faceUPtrs_.set(i))
{
FatalErrorIn
(
"void solutionControl::reconstructTransientVelocity"
"\n("
"\n volVectorField& U,"
"\n const volVectorField& ddtUEqn,"
"\n const volScalarField& rAU,"
"\n const volScalarField& p"
"\n const surfaceScalarField& phi"
"\n) const"
) << "faceU not calculated for field " << UName
<< ". Make sure you have called"
<< " calcTransientConsistentFlux(...) before calling this function."
<< exit(FatalError);
}
// Get faceU field. Note: no need to check whether the entry has been found
// since we have just checked
surfaceVectorField& faceU = faceUPtrs_[i];
// First interpolate the reconstructed velocity on the faces
faceU = fvc::interpolate(U);
// Replace the normal component with conservative flux
const surfaceVectorField& Sf = mesh_.Sf();
const surfaceVectorField rSf = Sf/magSqr(Sf);
// Subtract interpolated normal component
faceU -= (Sf & faceU)*rSf;
// Now that the normal component is zero, add the normal component from
// conservative flux
faceU += phi*rSf;
// If the mesh is moving, flux needs to be relative before boundary
// conditions for velocity are corrected. VV and IG, 4/Jan/2016.
fvc::makeRelative(phi, U);
// Correct boundary conditions with relative flux
U.correctBoundaryConditions();
}
void Foam::solutionControl::reconstructSteadyVelocity
(
volVectorField& U,
const volScalarField& rAU,
const volScalarField& p
) const
{
// Reconstruct the velocity field
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
// Note: no need to store and update faceU field for steady state run
}
const Foam::surfaceScalarField& Foam::solutionControl::aCoeff
(
const word& UName
) const
{
// Get corresponding index
const label i = indices_[UName];
if (!aCoeffPtrs_.set(i))
{
FatalErrorIn
(
"const surfaceScalarField& solutionControl::aCoeff() const"
) << "aCoeff not calculated for field " << UName
<< ". Make sure you have called"
<< " calcTransientConsistentFlux(...) or "
<< " calcSteadyConsistentFlux(...) before calling aCoeff()."
<< exit(FatalError);
}
return aCoeffPtrs_[i];
}
// ************************************************************************* //

2
src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C
View file

@ -1375,6 +1375,8 @@ Foam::BlockSolverPerformance<Type> Foam::fvBlockMatrix<Type>::solve
// Print performance
solverPerf.print();
psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerf);
return solverPerf;
}

13
src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
View file

@ -69,7 +69,7 @@ Foam::lduSolverPerformance Foam::fvMatrix<Type>::solve
// Complete matrix assembly. HJ, 17/Apr/2012
this->completeAssembly();
lduSolverPerformance solverPerfVec
BlockSolverPerformance<Type> solverPerfVec
(
"fvMatrix<Type>::solve",
psi_.name()
@ -155,14 +155,7 @@ Foam::lduSolverPerformance Foam::fvMatrix<Type>::solve
solverPerf.print();
if
(
solverPerf.initialResidual() > solverPerfVec.initialResidual()
&& !solverPerf.singular()
)
{
solverPerfVec = solverPerf;
}
solverPerfVec.replace(cmpt, solverPerf);
psi.internalField().replace(cmpt, psiCmpt);
diag() = saveDiag;
@ -172,7 +165,7 @@ Foam::lduSolverPerformance Foam::fvMatrix<Type>::solve
psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerfVec);
return solverPerfVec;
return solverPerfVec.max();
}

50
src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.C
View file

@ -112,6 +112,36 @@ void Foam::BlockSolverPerformance<Type>::print() const
}
template<class Type>
void Foam::BlockSolverPerformance<Type>::replace
(
const Foam::label cmpt,
const Foam::BlockSolverPerformance<typename pTraits<Type>::cmptType>& bsp
)
{
initialResidual_.replace(cmpt, bsp.initialResidual());
finalResidual_.replace(cmpt, bsp.finalResidual());
singular_ = singular_ || bsp.singular();
}
template<class Type>
Foam::BlockSolverPerformance<typename Foam::pTraits<Type>::cmptType>
Foam::BlockSolverPerformance<Type>::max()
{
return BlockSolverPerformance<typename pTraits<Type>::cmptType>
(
solverName_,
fieldName_,
cmptMax(initialResidual_),
cmptMax(finalResidual_),
nIterations_,
converged_,
singular_
);
}
template<class Type>
bool Foam::BlockSolverPerformance<Type>::operator!=
(
@ -131,6 +161,26 @@ bool Foam::BlockSolverPerformance<Type>::operator!=
}
template<class Type>
typename Foam::BlockSolverPerformance<Type> Foam::max
(
const typename Foam::BlockSolverPerformance<Type>& bsp1,
const typename Foam::BlockSolverPerformance<Type>& bsp2
)
{
return BlockSolverPerformance<Type>
(
bsp1.solverName(),
bsp1.fieldName(),
max(bsp1.initialResidual(), bsp2.initialResidual()),
max(bsp1.finalResidual(), bsp2.finalResidual()),
max(bsp1.nIterations(), bsp2.nIterations()),
bsp1.converged() && bsp2.converged(),
bsp1.singular() || bsp2.singular()
);
}
template<class Type>
Foam::Istream& Foam::operator>>
(

28
src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.H
View file

@ -50,6 +50,13 @@ namespace Foam
template<class Type>
class BlockSolverPerformance;
template<class Type>
BlockSolverPerformance<Type> max
(
const BlockSolverPerformance<Type>&,
const BlockSolverPerformance<Type>&
);
template<class Type>
Istream& operator>>(Istream&, BlockSolverPerformance<Type>&);
@ -217,6 +224,17 @@ public:
//- Print summary of solver performance
void print() const;
//- Replace component based on the minimal BlockSolverPerformance
void replace
(
const label cmpt,
const BlockSolverPerformance<typename pTraits<Type>::cmptType>& sp
);
//- Return the summary maximum of BlockSolverPerformance<Type>
// Effectively it will mostly return BlockSolverPerformance<scalar>
BlockSolverPerformance<typename pTraits<Type>::cmptType> max();
// Member Operators
@ -225,6 +243,16 @@ public:
// Friend functions
//- Return the element-wise max of two BlockSolverPerformance<Type>s
friend BlockSolverPerformance<Type> Foam::max <Type>
(
const BlockSolverPerformance<Type>&,
const BlockSolverPerformance<Type>&
);
// Ostream operator
friend Istream& operator>> <Type>
(
Istream&,

13
src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformances.C
View file

@ -31,22 +31,27 @@ Author
#include "BlockSolverPerformances.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceScalar, 1);
defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceVector, 1);
defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceSphericalTensor, 1);
defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceSymmTensor, 1);
defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceTensor, 1);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
template<>
Foam::BlockSolverPerformance<Foam::scalar>
Foam::BlockSolverPerformance<Foam::scalar>::max()
{
return *this;
}
// ************************************************************************* //

45
src/foam/matrices/solution/solution.C
View file

@ -422,49 +422,4 @@ Foam::dictionary& Foam::solution::solverPerformanceDict() const
return solverPerformance_;
}
void Foam::solution::setSolverPerformance
(
const word& name,
const lduSolverPerformance& sp
) const
{
List<lduSolverPerformance> perfs;
if (prevTimeIndex_ != this->time().timeIndex())
{
// Reset solver performance between iterations
prevTimeIndex_ = this->time().timeIndex();
solverPerformance_.clear();
}
else
{
solverPerformance_.readIfPresent(name, perfs);
}
// Only the first iteration and the current iteration residuals are
// required, so the current iteration residual replaces the previous one and
// only the first iteration is always present, VS 2017-11-27
if (perfs.size() < 2)
{
// Append to list
perfs.setSize(perfs.size() + 1, sp);
}
else
{
perfs.last() = sp;
}
solverPerformance_.set(name, perfs);
}
void Foam::solution::setSolverPerformance
(
const lduSolverPerformance& sp
) const
{
setSolverPerformance(sp.fieldName(), sp);
}
// ************************************************************************* //

13
src/foam/matrices/solution/solution.H
View file

@ -37,7 +37,8 @@ SourceFiles
#include "IOdictionary.H"
#include "debugSwitch.H"
#include "lduSolverPerformance.H"
#include "foamTime.H"
#include "BlockSolverPerformance.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -156,17 +157,19 @@ public:
// checking
dictionary& solverPerformanceDict() const;
//- Add/set the solverPerformance entry for the named field
//- Add/set the BlockSolverPerformance entry for the named field
template<class Type>
void setSolverPerformance
(
const word& name,
const lduSolverPerformance&
const BlockSolverPerformance<Type>&
) const;
//- Add/set the solverPerformance entry, using its fieldName
//- Add/set the BlockSolverPerformance entry, using its fieldName
template<class Type>
void setSolverPerformance
(
const lduSolverPerformance&
const BlockSolverPerformance<Type>&
) const;

37
src/foam/matrices/solution/solutionTemplates.C
View file

@ -45,6 +45,43 @@ void Foam::solution::cachePrintMessage
}
template<class Type>
void Foam::solution::setSolverPerformance
(
const word& name,
const BlockSolverPerformance<Type>& sp
) const
{
List<BlockSolverPerformance<Type> > perfs;
if (prevTimeIndex_ != this->time().timeIndex())
{
// Reset solver performance between iterations
prevTimeIndex_ = this->time().timeIndex();
solverPerformance_.clear();
}
else
{
solverPerformance_.readIfPresent(name, perfs);
}
// Append to list
perfs.setSize(perfs.size() + 1, sp);
solverPerformance_.set(name, perfs);
}
template<class Type>
void Foam::solution::setSolverPerformance
(
const BlockSolverPerformance<Type>& sp
) const
{
setSolverPerformance(sp.fieldName(), sp);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

15
src/tetFiniteElement/tetFemMatrix/tetFemMatrixSolve.C
View file

@ -52,7 +52,7 @@ lduSolverPerformance tetFemMatrix<Type>::solve
this->check();
}
lduSolverPerformance solverPerfVec
BlockSolverPerformance<Type> solverPerfVec
(
"tetFemMatrix<Type>::solve",
psi_.name()
@ -152,14 +152,7 @@ lduSolverPerformance tetFemMatrix<Type>::solve
solverPerf.print();
if
(
solverPerf.initialResidual() > solverPerfVec.initialResidual()
&& !solverPerf.singular()
)
{
solverPerfVec = solverPerf;
}
solverPerfVec.replace(cmpt, solverPerf);
psi.internalField().replace(cmpt, psiCmpt);
@ -174,7 +167,9 @@ lduSolverPerformance tetFemMatrix<Type>::solve
psi.correctBoundaryConditions();
return solverPerfVec;
psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerfVec);
return solverPerfVec.max();
}

2
src/tetFiniteElement/tetFemMatrix/tetFemScalarMatrix.C
View file

@ -131,6 +131,8 @@ lduSolverPerformance tetFemMatrix<scalar>::solve
psi.correctBoundaryConditions();
psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerf);
return solverPerf;
}