00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include <OpenFOAM/smoothSolver.H>
00027
00028
00029
00030 namespace Foam
00031 {
00032 defineTypeNameAndDebug(smoothSolver, 0);
00033
00034 lduMatrix::solver::addsymMatrixConstructorToTable<smoothSolver>
00035 addsmoothSolverSymMatrixConstructorToTable_;
00036
00037 lduMatrix::solver::addasymMatrixConstructorToTable<smoothSolver>
00038 addsmoothSolverAsymMatrixConstructorToTable_;
00039 }
00040
00041
00042
00043
00044 Foam::smoothSolver::smoothSolver
00045 (
00046 const word& fieldName,
00047 const lduMatrix& matrix,
00048 const FieldField<Field, scalar>& interfaceBouCoeffs,
00049 const FieldField<Field, scalar>& interfaceIntCoeffs,
00050 const lduInterfaceFieldPtrsList& interfaces,
00051 const dictionary& solverControls
00052 )
00053 :
00054 lduMatrix::solver
00055 (
00056 fieldName,
00057 matrix,
00058 interfaceBouCoeffs,
00059 interfaceIntCoeffs,
00060 interfaces,
00061 solverControls
00062 )
00063 {
00064 readControls();
00065 }
00066
00067
00068
00069
00070 void Foam::smoothSolver::readControls()
00071 {
00072 lduMatrix::solver::readControls();
00073 nSweeps_ = controlDict_.lookupOrDefault<label>("nSweeps", 1);
00074 }
00075
00076
00077 Foam::lduMatrix::solverPerformance Foam::smoothSolver::solve
00078 (
00079 scalarField& psi,
00080 const scalarField& source,
00081 const direction cmpt
00082 ) const
00083 {
00084
00085 lduMatrix::solverPerformance solverPerf(typeName, fieldName_);
00086
00087
00088 if (nSweeps_ < 0)
00089 {
00090 autoPtr<lduMatrix::smoother> smootherPtr = lduMatrix::smoother::New
00091 (
00092 fieldName_,
00093 matrix_,
00094 interfaceBouCoeffs_,
00095 interfaceIntCoeffs_,
00096 interfaces_,
00097 controlDict_
00098 );
00099
00100 smootherPtr->smooth
00101 (
00102 psi,
00103 source,
00104 cmpt,
00105 -nSweeps_
00106 );
00107
00108 solverPerf.nIterations() -= nSweeps_;
00109 }
00110 else
00111 {
00112 scalar normFactor = 0;
00113
00114 {
00115 scalarField Apsi(psi.size());
00116 scalarField temp(psi.size());
00117
00118
00119 matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
00120
00121
00122 normFactor = this->normFactor(psi, source, Apsi, temp);
00123
00124
00125 solverPerf.initialResidual() = gSumMag(source - Apsi)/normFactor;
00126 solverPerf.finalResidual() = solverPerf.initialResidual();
00127 }
00128
00129 if (lduMatrix::debug >= 2)
00130 {
00131 Info<< " Normalisation factor = " << normFactor << endl;
00132 }
00133
00134
00135
00136 if (!solverPerf.checkConvergence(tolerance_, relTol_))
00137 {
00138 autoPtr<lduMatrix::smoother> smootherPtr = lduMatrix::smoother::New
00139 (
00140 fieldName_,
00141 matrix_,
00142 interfaceBouCoeffs_,
00143 interfaceIntCoeffs_,
00144 interfaces_,
00145 controlDict_
00146 );
00147
00148
00149 do
00150 {
00151 smootherPtr->smooth
00152 (
00153 psi,
00154 source,
00155 cmpt,
00156 nSweeps_
00157 );
00158
00159
00160 solverPerf.finalResidual() = gSumMag
00161 (
00162 matrix_.residual
00163 (
00164 psi,
00165 source,
00166 interfaceBouCoeffs_,
00167 interfaces_,
00168 cmpt
00169 )
00170 )/normFactor;
00171 } while
00172 (
00173 (solverPerf.nIterations() += nSweeps_) < maxIter_
00174 && !(solverPerf.checkConvergence(tolerance_, relTol_))
00175 );
00176 }
00177 }
00178
00179 return solverPerf;
00180 }
00181
00182
00183