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