OpenFOAM-4.x/src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C
2016-03-22 10:38:30 +00:00

160 lines
4.7 KiB
C

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 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 "SmoothSolver.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
Foam::SmoothSolver<Type, DType, LUType>::SmoothSolver
(
const word& fieldName,
const LduMatrix<Type, DType, LUType>& matrix,
const dictionary& solverDict
)
:
LduMatrix<Type, DType, LUType>::solver
(
fieldName,
matrix,
solverDict
),
nSweeps_(1)
{
readControls();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type, class DType, class LUType>
void Foam::SmoothSolver<Type, DType, LUType>::readControls()
{
LduMatrix<Type, DType, LUType>::solver::readControls();
this->readControl(this->controlDict_, nSweeps_, "nSweeps");
}
template<class Type, class DType, class LUType>
Foam::SolverPerformance<Type>
Foam::SmoothSolver<Type, DType, LUType>::solve(Field<Type>& psi) const
{
// --- Setup class containing solver performance data
SolverPerformance<Type> solverPerf
(
typeName,
this->fieldName_
);
// If the nSweeps_ is negative do a fixed number of sweeps
if (nSweeps_ < 0)
{
autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
smootherPtr = LduMatrix<Type, DType, LUType>::smoother::New
(
this->fieldName_,
this->matrix_,
this->controlDict_
);
smootherPtr->smooth(psi, -nSweeps_);
solverPerf.nIterations() -= nSweeps_;
}
else
{
Type normFactor = Zero;
{
Field<Type> Apsi(psi.size());
Field<Type> temp(psi.size());
// Calculate A.psi
this->matrix_.Amul(Apsi, psi);
// Calculate normalisation factor
normFactor = this->normFactor(psi, Apsi, temp);
// Calculate residual magnitude
solverPerf.initialResidual() = cmptDivide
(
gSumCmptMag(this->matrix_.source() - Apsi),
normFactor
);
solverPerf.finalResidual() = solverPerf.initialResidual();
}
if (LduMatrix<Type, DType, LUType>::debug >= 2)
{
Info<< " Normalisation factor = " << normFactor << endl;
}
// Check convergence, solve if not converged
if
(
this->minIter_ > 0
|| !solverPerf.checkConvergence(this->tolerance_, this->relTol_)
)
{
autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
smootherPtr = LduMatrix<Type, DType, LUType>::smoother::New
(
this->fieldName_,
this->matrix_,
this->controlDict_
);
// Smoothing loop
do
{
smootherPtr->smooth
(
psi,
nSweeps_
);
// Calculate the residual to check convergence
solverPerf.finalResidual() = cmptDivide
(
gSumCmptMag(this->matrix_.residual(psi)),
normFactor
);
} while
(
(
(solverPerf.nIterations() += nSweeps_) < this->maxIter_
&& !solverPerf.checkConvergence(this->tolerance_, this->relTol_)
)
|| solverPerf.nIterations() < this->minIter_
);
}
}
return solverPerf;
}
// ************************************************************************* //