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 :
38  LduMatrix<Type, DType, LUType>::solver
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 // ************************************************************************* //
Pre-declare SubField and related Field type.
Definition: Field.H:83
static autoPtr< preconditioner > New(const solver &sol, const dictionary &preconditionerDict)
Return a new preconditioner.
LduMatrix is a general matrix class in which the coefficients are stored as three arrays,...
Definition: LduMatrix.H:85
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
PCICG(const word &fieldName, const LduMatrix< Type, DType, LUType > &matrix, const dictionary &solverDict)
Construct from matrix components and solver data dictionary.
Definition: PCICG.C:32
virtual SolverPerformance< Type > solve(Field< Type > &psi) const
Solve the matrix with this solver.
Definition: PCICG.C:51
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics.
const labelType & nIterations() const
Return number of iterations.
bool checkSingularity(const Type &residual)
Singularity test.
const Type & initialResidual() const
Return initial residual.
const Type & finalResidual() const
Return final residual.
static const scalar vsmall_
Very small Type for the use in solvers.
static const scalar great_
Large Type for the use in solvers.
bool checkConvergence(const Type &tolerance, const Type &relTolerance)
Check, store and return convergence.
iterator begin()
Return an iterator to begin traversing the UList.
Definition: UListI.H:216
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: autoPtr.H:51
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:60
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
A class representing the concept of 1 (scalar(1)) used to avoid unnecessary manipulations for objects...
Definition: one.H:51
Traits class for primitives.
Definition: pTraits.H:53
Abstract base class for run-time selectable region solvers.
Definition: solver.H:56
A class for handling words, derived from string.
Definition: word.H:62
const volScalarField & psi
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
tmp< DimensionedField< scalar, GeoMesh > > stabilise(const DimensionedField< scalar, GeoMesh > &dsf, const dimensioned< scalar > &ds)
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
messageStream Info
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionSet cmptMag(const dimensionSet &)
Definition: dimensionSet.C:296
dimensioned< Type > cmptDivide(const dimensioned< Type > &, const dimensioned< Type > &)
Type gSumCmptMag(const UList< Type > &f, const label comm)
Type gSumCmptProd(const UList< Type > &f1, const UList< Type > &f2, const label comm)