From e3115846f5cad2bf95fcd999952a827c9dd6e667 Mon Sep 17 00:00:00 2001 From: Vanja Skuric Date: Sun, 11 Feb 2018 18:59:15 +0100 Subject: [PATCH] New feature: generalization of storage for solverPerformanceDict based on templated BlockSolverPerformance (by Pascal Beckstein) 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: https://github.com/OpenFOAM/OpenFOAM-dev/commit/1944b09bb530a8b9a2546ac3b196b137668d329f --- .../faMatrices/faMatrix/faMatrixSolve.C | 15 +- .../faScalarMatrix/faScalarMatrix.C | 4 + .../solutionControl/solutionControl.C | 556 ++---------------- .../fvMatrices/fvBlockMatrix/fvBlockMatrix.C | 2 + .../fvMatrices/fvMatrix/fvMatrixSolve.C | 13 +- .../BlockLduSolver/BlockSolverPerformance.C | 50 ++ .../BlockLduSolver/BlockSolverPerformance.H | 28 + .../BlockLduSolver/BlockSolverPerformances.C | 13 +- src/foam/matrices/solution/solution.C | 45 -- src/foam/matrices/solution/solution.H | 13 +- .../matrices/solution/solutionTemplates.C | 37 ++ .../tetFemMatrix/tetFemMatrixSolve.C | 15 +- .../tetFemMatrix/tetFemScalarMatrix.C | 2 + 13 files changed, 187 insertions(+), 606 deletions(-) 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; }