PBiCICG.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 "PBiCICG.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> pT(nCells, Zero);
72  Type* __restrict__ pTPtr = pT.begin();
73 
74  Field<Type> wA(nCells);
75  Type* __restrict__ wAPtr = wA.begin();
76 
77  Field<Type> wT(nCells);
78  Type* __restrict__ wTPtr = wT.begin();
79 
80  Type wArT = solverPerf.great_*pTraits<Type>::one;
81  Type wArTold = wArT;
82 
83  // --- Calculate A.psi and T.psi
84  this->matrix_.Amul(wA, psi);
85  this->matrix_.Tmul(wT, psi);
86 
87  // --- Calculate initial residual and transpose residual fields
88  Field<Type> rA(this->matrix_.source() - wA);
89  Field<Type> rT(this->matrix_.source() - wT);
90  Type* __restrict__ rAPtr = rA.begin();
91  Type* __restrict__ rTPtr = rT.begin();
92 
93  // --- Calculate normalisation factor
94  Type normFactor = this->normFactor(psi, wA, pA);
95 
97  {
98  Info<< " Normalisation factor = " << normFactor << endl;
99  }
100 
101  // --- Calculate normalised residual norm
102  solverPerf.initialResidual() = cmptDivide(gSumCmptMag(rA), normFactor);
103  solverPerf.finalResidual() = solverPerf.initialResidual();
104 
105  // --- Check convergence, solve if not converged
106  if (!solverPerf.checkConvergence(this->tolerance_, this->relTol_))
107  {
108  // --- Select and construct the preconditioner
111  (
112  *this,
113  this->controlDict_
114  );
115 
116  // --- Solver iteration
117  do
118  {
119  // --- Store previous wArT
120  wArTold = wArT;
121 
122  // --- Precondition residuals
123  preconPtr->precondition(wA, rA);
124  preconPtr->preconditionT(wT, rT);
125 
126  // --- Update search directions:
127  wArT = gSumCmptProd(wA, rT);
128 
129  if (nIter == 0)
130  {
131  for (label cell=0; cell<nCells; cell++)
132  {
133  pAPtr[cell] = wAPtr[cell];
134  pTPtr[cell] = wTPtr[cell];
135  }
136  }
137  else
138  {
139  Type beta = cmptDivide
140  (
141  wArT,
142  stabilise(wArTold, solverPerf.vsmall_)
143  );
144 
145  for (label cell=0; cell<nCells; cell++)
146  {
147  pAPtr[cell] = wAPtr[cell] + cmptMultiply(beta, pAPtr[cell]);
148  pTPtr[cell] = wTPtr[cell] + cmptMultiply(beta, pTPtr[cell]);
149  }
150  }
151 
152 
153  // --- Update preconditioned residuals
154  this->matrix_.Amul(wA, pA);
155  this->matrix_.Tmul(wT, pT);
156 
157  Type wApT = gSumCmptProd(wA, pT);
158 
159  // --- Test for singularity
160  if
161  (
162  solverPerf.checkSingularity
163  (
164  cmptDivide(cmptMag(wApT), normFactor)
165  )
166  )
167  {
168  break;
169  }
170 
171 
172  // --- Update solution and residual:
173 
174  Type alpha = cmptDivide
175  (
176  wArT,
177  stabilise(wApT, solverPerf.vsmall_)
178  );
179 
180  for (label cell=0; cell<nCells; cell++)
181  {
182  psiPtr[cell] += cmptMultiply(alpha, pAPtr[cell]);
183  rAPtr[cell] -= cmptMultiply(alpha, wAPtr[cell]);
184  rTPtr[cell] -= cmptMultiply(alpha, wTPtr[cell]);
185  }
186 
187  solverPerf.finalResidual() =
188  cmptDivide(gSumCmptMag(rA), normFactor);
189 
190  } while
191  (
192  nIter++ < this->maxIter_
193  && !(solverPerf.checkConvergence(this->tolerance_, this->relTol_))
194  );
195  }
196 
197  solverPerf.nIterations() =
199 
200  return solverPerf;
201 }
202 
203 
204 // ************************************************************************* //
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
PBiCICG(const word &fieldName, const LduMatrix< Type, DType, LUType > &matrix, const dictionary &solverDict)
Construct from matrix components and solver data dictionary.
Definition: PBiCICG.C:32
virtual SolverPerformance< Type > solve(Field< Type > &psi) const
Solve the matrix with this solver.
Definition: PBiCICG.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))
static const zero Zero
Definition: zero.H:97
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)