diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C b/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C index 948c04b51..126ead1b3 100644 --- a/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C +++ b/src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C @@ -60,7 +60,7 @@ lduSolverPerformance faMatrix::solve(const dictionary& solverControls) << endl; } - lduSolverPerformance solverPerfVec + BlockSolverPerformance solverPerfVec ( "faMatrix::solve", psi_.name() @@ -134,14 +134,7 @@ lduSolverPerformance faMatrix::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::solve(const dictionary& solverControls) psi.correctBoundaryConditions(); - return solverPerfVec; + psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerfVec); + + return solverPerfVec.max(); } diff --git a/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C b/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C index 5b89e8903..85bd8d706 100644 --- a/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C +++ b/src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C @@ -130,6 +130,8 @@ lduSolverPerformance faMatrix::faSolver::solve psi.correctBoundaryConditions(); + psi.mesh().solutionDict().setSolverPerformance(psi.name(), solverPerf); + return solverPerf; } @@ -177,6 +179,8 @@ lduSolverPerformance faMatrix::solve psi.correctBoundaryConditions(); + psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerf); + return solverPerf; } diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C index 67b1bcf30..1882b95fc 100644 --- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C +++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C @@ -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(); storePrevIter(); storePrevIter(); + + storePrevIter(); + storePrevIter(); + storePrevIter(); + storePrevIter(); + + storePrevIter(); + storePrevIter(); + storePrevIter(); + storePrevIter(); + + storePrevIter(); + storePrevIter(); + storePrevIter(); + storePrevIter(); } @@ -196,7 +208,7 @@ void Foam::solutionControl::maxTypeResidual if (mesh_.foundObject(fieldName)) { - const List sp(data); + const List > sp(data); firstRes = cmptMax(sp.first().initialResidual()); lastRes = cmptMax(sp.last().initialResidual()); } @@ -218,119 +230,25 @@ Foam::scalar Foam::solutionControl::maxResidual maxTypeResidual(fieldName, data, firstRes, lastRes); maxTypeResidual(fieldName, data, firstRes, lastRes); + maxTypeResidual(fieldName, data, firstRes, lastRes); + maxTypeResidual(fieldName, data, firstRes, lastRes); + maxTypeResidual(fieldName, data, firstRes, lastRes); + maxTypeResidual(fieldName, data, firstRes, lastRes); + + maxTypeResidual(fieldName, data, firstRes, lastRes); + maxTypeResidual(fieldName, data, firstRes, lastRes); + maxTypeResidual(fieldName, data, firstRes, lastRes); + maxTypeResidual(fieldName, data, firstRes, lastRes); + + maxTypeResidual(fieldName, data, firstRes, lastRes); + maxTypeResidual(fieldName, data, firstRes, lastRes); + maxTypeResidual(fieldName, data, firstRes, lastRes); + maxTypeResidual(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(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(Up) - || isA(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]; -} - - // ************************************************************************* // diff --git a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C index 89a2a816d..ba23a5892 100644 --- a/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C +++ b/src/finiteVolume/fvMatrices/fvBlockMatrix/fvBlockMatrix.C @@ -1375,6 +1375,8 @@ Foam::BlockSolverPerformance Foam::fvBlockMatrix::solve // Print performance solverPerf.print(); + psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerf); + return solverPerf; } diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C index 7ecc96951..a8130abb4 100644 --- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C +++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C @@ -69,7 +69,7 @@ Foam::lduSolverPerformance Foam::fvMatrix::solve // Complete matrix assembly. HJ, 17/Apr/2012 this->completeAssembly(); - lduSolverPerformance solverPerfVec + BlockSolverPerformance solverPerfVec ( "fvMatrix::solve", psi_.name() @@ -155,14 +155,7 @@ Foam::lduSolverPerformance Foam::fvMatrix::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::solve psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerfVec); - return solverPerfVec; + return solverPerfVec.max(); } diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.C index 2881f8096..4bdbcaa9a 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.C @@ -112,6 +112,36 @@ void Foam::BlockSolverPerformance::print() const } +template +void Foam::BlockSolverPerformance::replace +( + const Foam::label cmpt, + const Foam::BlockSolverPerformance::cmptType>& bsp +) +{ + initialResidual_.replace(cmpt, bsp.initialResidual()); + finalResidual_.replace(cmpt, bsp.finalResidual()); + singular_ = singular_ || bsp.singular(); +} + + +template +Foam::BlockSolverPerformance::cmptType> +Foam::BlockSolverPerformance::max() +{ + return BlockSolverPerformance::cmptType> + ( + solverName_, + fieldName_, + cmptMax(initialResidual_), + cmptMax(finalResidual_), + nIterations_, + converged_, + singular_ + ); +} + + template bool Foam::BlockSolverPerformance::operator!= ( @@ -131,6 +161,26 @@ bool Foam::BlockSolverPerformance::operator!= } +template +typename Foam::BlockSolverPerformance Foam::max +( + const typename Foam::BlockSolverPerformance& bsp1, + const typename Foam::BlockSolverPerformance& bsp2 +) +{ + return BlockSolverPerformance + ( + 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 Foam::Istream& Foam::operator>> ( diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.H b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.H index 9ee67da3d..331f0c873 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.H @@ -50,6 +50,13 @@ namespace Foam template class BlockSolverPerformance; +template +BlockSolverPerformance max +( + const BlockSolverPerformance&, + const BlockSolverPerformance& +); + template Istream& operator>>(Istream&, BlockSolverPerformance&); @@ -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::cmptType>& sp + ); + + //- Return the summary maximum of BlockSolverPerformance + // Effectively it will mostly return BlockSolverPerformance + BlockSolverPerformance::cmptType> max(); + // Member Operators @@ -225,6 +243,16 @@ public: // Friend functions + //- Return the element-wise max of two BlockSolverPerformances + friend BlockSolverPerformance Foam::max + ( + const BlockSolverPerformance&, + const BlockSolverPerformance& + ); + + + // Ostream operator + friend Istream& operator>> ( Istream&, diff --git a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformances.C b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformances.C index 73049b778..633e797c2 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformances.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformances.C @@ -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::BlockSolverPerformance::max() +{ + return *this; +} + // ************************************************************************* // diff --git a/src/foam/matrices/solution/solution.C b/src/foam/matrices/solution/solution.C index 7721094c9..100d0713a 100644 --- a/src/foam/matrices/solution/solution.C +++ b/src/foam/matrices/solution/solution.C @@ -422,49 +422,4 @@ Foam::dictionary& Foam::solution::solverPerformanceDict() const return solverPerformance_; } - -void Foam::solution::setSolverPerformance -( - const word& name, - const lduSolverPerformance& sp -) const -{ - List 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); -} - // ************************************************************************* // diff --git a/src/foam/matrices/solution/solution.H b/src/foam/matrices/solution/solution.H index bf4206213..8461d3c07 100644 --- a/src/foam/matrices/solution/solution.H +++ b/src/foam/matrices/solution/solution.H @@ -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 void setSolverPerformance ( const word& name, - const lduSolverPerformance& + const BlockSolverPerformance& ) const; - //- Add/set the solverPerformance entry, using its fieldName + //- Add/set the BlockSolverPerformance entry, using its fieldName + template void setSolverPerformance ( - const lduSolverPerformance& + const BlockSolverPerformance& ) const; diff --git a/src/foam/matrices/solution/solutionTemplates.C b/src/foam/matrices/solution/solutionTemplates.C index 67cdd879f..6f85c4166 100644 --- a/src/foam/matrices/solution/solutionTemplates.C +++ b/src/foam/matrices/solution/solutionTemplates.C @@ -45,6 +45,43 @@ void Foam::solution::cachePrintMessage } +template +void Foam::solution::setSolverPerformance +( + const word& name, + const BlockSolverPerformance& sp +) const +{ + List > 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 +void Foam::solution::setSolverPerformance +( + const BlockSolverPerformance& sp +) const +{ + setSolverPerformance(sp.fieldName(), sp); +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/tetFiniteElement/tetFemMatrix/tetFemMatrixSolve.C b/src/tetFiniteElement/tetFemMatrix/tetFemMatrixSolve.C index f136f758f..6a8ae9d8e 100644 --- a/src/tetFiniteElement/tetFemMatrix/tetFemMatrixSolve.C +++ b/src/tetFiniteElement/tetFemMatrix/tetFemMatrixSolve.C @@ -52,7 +52,7 @@ lduSolverPerformance tetFemMatrix::solve this->check(); } - lduSolverPerformance solverPerfVec + BlockSolverPerformance solverPerfVec ( "tetFemMatrix::solve", psi_.name() @@ -152,14 +152,7 @@ lduSolverPerformance tetFemMatrix::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::solve psi.correctBoundaryConditions(); - return solverPerfVec; + psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerfVec); + + return solverPerfVec.max(); } diff --git a/src/tetFiniteElement/tetFemMatrix/tetFemScalarMatrix.C b/src/tetFiniteElement/tetFemMatrix/tetFemScalarMatrix.C index 0199325bf..9fdcd40ad 100644 --- a/src/tetFiniteElement/tetFemMatrix/tetFemScalarMatrix.C +++ b/src/tetFiniteElement/tetFemMatrix/tetFemScalarMatrix.C @@ -131,6 +131,8 @@ lduSolverPerformance tetFemMatrix::solve psi.correctBoundaryConditions(); + psi_.mesh().solutionDict().setSolverPerformance(psi_.name(), solverPerf); + return solverPerf; }