Go to the documentation of this file.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/PCG.H>
00027
00028
00029
00030 namespace Foam
00031 {
00032 defineTypeNameAndDebug(PCG, 0);
00033
00034 lduMatrix::solver::addsymMatrixConstructorToTable<PCG>
00035 addPCGSymMatrixConstructorToTable_;
00036 }
00037
00038
00039
00040
00041 Foam::PCG::PCG
00042 (
00043 const word& fieldName,
00044 const lduMatrix& matrix,
00045 const FieldField<Field, scalar>& interfaceBouCoeffs,
00046 const FieldField<Field, scalar>& interfaceIntCoeffs,
00047 const lduInterfaceFieldPtrsList& interfaces,
00048 const dictionary& solverControls
00049 )
00050 :
00051 lduMatrix::solver
00052 (
00053 fieldName,
00054 matrix,
00055 interfaceBouCoeffs,
00056 interfaceIntCoeffs,
00057 interfaces,
00058 solverControls
00059 )
00060 {}
00061
00062
00063
00064
00065 Foam::lduMatrix::solverPerformance Foam::PCG::solve
00066 (
00067 scalarField& psi,
00068 const scalarField& source,
00069 const direction cmpt
00070 ) const
00071 {
00072
00073 lduMatrix::solverPerformance solverPerf
00074 (
00075 lduMatrix::preconditioner::getName(controlDict_) + typeName,
00076 fieldName_
00077 );
00078
00079 register label nCells = psi.size();
00080
00081 scalar* __restrict__ psiPtr = psi.begin();
00082
00083 scalarField pA(nCells);
00084 scalar* __restrict__ pAPtr = pA.begin();
00085
00086 scalarField wA(nCells);
00087 scalar* __restrict__ wAPtr = wA.begin();
00088
00089 scalar wArA = matrix_.great_;
00090 scalar wArAold = wArA;
00091
00092
00093 matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
00094
00095
00096 scalarField rA(source - wA);
00097 scalar* __restrict__ rAPtr = rA.begin();
00098
00099
00100 scalar normFactor = this->normFactor(psi, source, wA, pA);
00101
00102 if (lduMatrix::debug >= 2)
00103 {
00104 Info<< " Normalisation factor = " << normFactor << endl;
00105 }
00106
00107
00108 solverPerf.initialResidual() = gSumMag(rA)/normFactor;
00109 solverPerf.finalResidual() = solverPerf.initialResidual();
00110
00111
00112 if (!solverPerf.checkConvergence(tolerance_, relTol_))
00113 {
00114
00115 autoPtr<lduMatrix::preconditioner> preconPtr =
00116 lduMatrix::preconditioner::New
00117 (
00118 *this,
00119 controlDict_
00120 );
00121
00122
00123 do
00124 {
00125
00126 wArAold = wArA;
00127
00128
00129 preconPtr->precondition(wA, rA, cmpt);
00130
00131
00132 wArA = gSumProd(wA, rA);
00133
00134 if (solverPerf.nIterations() == 0)
00135 {
00136 for (register label cell=0; cell<nCells; cell++)
00137 {
00138 pAPtr[cell] = wAPtr[cell];
00139 }
00140 }
00141 else
00142 {
00143 scalar beta = wArA/wArAold;
00144
00145 for (register label cell=0; cell<nCells; cell++)
00146 {
00147 pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
00148 }
00149 }
00150
00151
00152
00153 matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);
00154
00155 scalar wApA = gSumProd(wA, pA);
00156
00157
00158
00159 if (solverPerf.checkSingularity(mag(wApA)/normFactor)) break;
00160
00161
00162
00163
00164 scalar alpha = wArA/wApA;
00165
00166 for (register label cell=0; cell<nCells; cell++)
00167 {
00168 psiPtr[cell] += alpha*pAPtr[cell];
00169 rAPtr[cell] -= alpha*wAPtr[cell];
00170 }
00171
00172 solverPerf.finalResidual() = gSumMag(rA)/normFactor;
00173
00174 } while
00175 (
00176 solverPerf.nIterations()++ < maxIter_
00177 && !(solverPerf.checkConvergence(tolerance_, relTol_))
00178 );
00179 }
00180
00181 return solverPerf;
00182 }
00183
00184
00185