PCICG.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "PCICG.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<class Type, class DType, class LUType>
32 (
33  const word& fieldName,
34  const LduMatrix<Type, DType, LUType>& matrix,
35  const dictionary& solverDict
36 )
37 :
39  (
40  fieldName,
41  matrix,
42  solverDict
43  )
44 {}
45 
46 
47 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
48 
49 template<class Type, class DType, class LUType>
52 {
53  word preconditionerName(this->controlDict_.lookup("preconditioner"));
54 
55  // --- Setup class containing solver performance data
56  SolverPerformance<Type> solverPerf
57  (
58  preconditionerName + typeName,
59  this->fieldName_
60  );
61 
62  label nIter = 0;
63 
64  label nCells = psi.size();
65 
66  Type* __restrict__ psiPtr = psi.begin();
67 
68  Field<Type> pA(nCells);
69  Type* __restrict__ pAPtr = pA.begin();
70 
71  Field<Type> wA(nCells);
72  Type* __restrict__ wAPtr = wA.begin();
73 
74  Type wArA = solverPerf.great_*pTraits<Type>::one;
75  Type wArAold = wArA;
76 
77  // --- Calculate A.psi
78  this->matrix_.Amul(wA, psi);
79 
80  // --- Calculate initial residual field
81  Field<Type> rA(this->matrix_.source() - wA);
82  Type* __restrict__ rAPtr = rA.begin();
83 
84  // --- Calculate normalisation factor
85  Type normFactor = this->normFactor(psi, wA, pA);
86 
88  {
89  Info<< " Normalisation factor = " << normFactor << endl;
90  }
91 
92  // --- Calculate normalised residual norm
93  solverPerf.initialResidual() = cmptDivide(gSumCmptMag(rA), normFactor);
94  solverPerf.finalResidual() = solverPerf.initialResidual();
95 
96  // --- Check convergence, solve if not converged
97  if
98  (
99  this->minIter_ > 0
100  || !solverPerf.checkConvergence(this->tolerance_, this->relTol_)
101  )
102  {
103  // --- Select and construct the preconditioner
106  (
107  *this,
108  this->controlDict_
109  );
110 
111  // --- Solver iteration
112  do
113  {
114  // --- Store previous wArA
115  wArAold = wArA;
116 
117  // --- Precondition residual
118  preconPtr->precondition(wA, rA);
119 
120  // --- Update search directions:
121  wArA = gSumCmptProd(wA, rA);
122 
123  if (nIter == 0)
124  {
125  for (label cell=0; cell<nCells; cell++)
126  {
127  pAPtr[cell] = wAPtr[cell];
128  }
129  }
130  else
131  {
132  Type beta = cmptDivide
133  (
134  wArA,
135  stabilise(wArAold, solverPerf.vsmall_)
136  );
137 
138  for (label cell=0; cell<nCells; cell++)
139  {
140  pAPtr[cell] = wAPtr[cell] + cmptMultiply(beta, pAPtr[cell]);
141  }
142  }
143 
144 
145  // --- Update preconditioned residual
146  this->matrix_.Amul(wA, pA);
147 
148  Type wApA = gSumCmptProd(wA, pA);
149 
150 
151  // --- Test for singularity
152  if
153  (
154  solverPerf.checkSingularity
155  (
156  cmptDivide(cmptMag(wApA), normFactor)
157  )
158  )
159  {
160  break;
161  }
162 
163 
164  // --- Update solution and residual:
165 
166  Type alpha = cmptDivide
167  (
168  wArA,
169  stabilise(wApA, solverPerf.vsmall_)
170  );
171 
172  for (label cell=0; cell<nCells; cell++)
173  {
174  psiPtr[cell] += cmptMultiply(alpha, pAPtr[cell]);
175  rAPtr[cell] -= cmptMultiply(alpha, wAPtr[cell]);
176  }
177 
178  solverPerf.finalResidual() =
179  cmptDivide(gSumCmptMag(rA), normFactor);
180 
181  } while
182  (
183  (
184  nIter++ < this->maxIter_
185  && !solverPerf.checkConvergence(this->tolerance_, this->relTol_)
186  )
187  || nIter < this->minIter_
188  );
189  }
190 
191  solverPerf.nIterations() =
193 
194  return solverPerf;
195 }
196 
197 
198 // ************************************************************************* //
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
Traits class for primitives.
Definition: pTraits.H:50
Type gSumCmptMag(const UList< Type > &f, const label comm)
virtual SolverPerformance< Type > solve(Field< Type > &psi) const
Solve the matrix with this solver.
Definition: PCICG.C:51
dimensioned< Type > cmptDivide(const dimensioned< Type > &, const dimensioned< Type > &)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
Type gSumCmptProd(const UList< Type > &f1, const UList< Type > &f2, const label comm)
Pre-declare SubField and related Field type.
Definition: Field.H:57
A class for handling words, derived from string.
Definition: word.H:59
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics...
iterator begin()
Return an iterator to begin traversing the UList.
Definition: UListI.H:216
Preconditioned conjugate gradient solver for symmetric lduMatrices using a run-time selectable precon...
Definition: PCICG.H:50
Abstract base-class for LduMatrix solvers.
Definition: LduMatrix.H:112
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
tmp< DimensionedField< scalar, GeoMesh > > stabilise(const DimensionedField< scalar, GeoMesh > &dsf, const dimensioned< scalar > &ds)
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
LduMatrix is a general matrix class in which the coefficients are stored as three arrays...
Definition: LduMatrix.H:69
messageStream Info
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
const volScalarField & psi
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
A class representing the concept of 1 (scalar(1)) used to avoid unnecessary manipulations for objects...
Definition: one.H:50