PCG.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2015 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 "PCG.H"
27 
28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
29 
30 namespace Foam
31 {
32  defineTypeNameAndDebug(PCG, 0);
33 
34  lduMatrix::solver::addsymMatrixConstructorToTable<PCG>
36 }
37 
38 
39 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
40 
41 Foam::PCG::PCG
42 (
43  const word& fieldName,
44  const lduMatrix& matrix,
45  const FieldField<Field, scalar>& interfaceBouCoeffs,
46  const FieldField<Field, scalar>& interfaceIntCoeffs,
47  const lduInterfaceFieldPtrsList& interfaces,
48  const dictionary& solverControls
49 )
50 :
52  (
53  fieldName,
54  matrix,
55  interfaceBouCoeffs,
56  interfaceIntCoeffs,
57  interfaces,
58  solverControls
59  )
60 {}
61 
62 
63 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
64 
66 (
67  scalarField& psi,
68  const scalarField& source,
69  const direction cmpt
70 ) const
71 {
72  // --- Setup class containing solver performance data
73  solverPerformance solverPerf
74  (
75  lduMatrix::preconditioner::getName(controlDict_) + typeName,
76  fieldName_
77  );
78 
79  label nCells = psi.size();
80 
81  scalar* __restrict__ psiPtr = psi.begin();
82 
83  scalarField pA(nCells);
84  scalar* __restrict__ pAPtr = pA.begin();
85 
86  scalarField wA(nCells);
87  scalar* __restrict__ wAPtr = wA.begin();
88 
89  scalar wArA = solverPerf.great_;
90  scalar wArAold = wArA;
91 
92  // --- Calculate A.psi
93  matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
94 
95  // --- Calculate initial residual field
96  scalarField rA(source - wA);
97  scalar* __restrict__ rAPtr = rA.begin();
98 
99  // --- Calculate normalisation factor
100  scalar normFactor = this->normFactor(psi, source, wA, pA);
101 
102  if (lduMatrix::debug >= 2)
103  {
104  Info<< " Normalisation factor = " << normFactor << endl;
105  }
106 
107  // --- Calculate normalised residual norm
108  solverPerf.initialResidual() =
109  gSumMag(rA, matrix().mesh().comm())
110  /normFactor;
111  solverPerf.finalResidual() = solverPerf.initialResidual();
112 
113  // --- Check convergence, solve if not converged
114  if
115  (
116  minIter_ > 0
117  || !solverPerf.checkConvergence(tolerance_, relTol_)
118  )
119  {
120  // --- Select and construct the preconditioner
123  (
124  *this,
125  controlDict_
126  );
127 
128  // --- Solver iteration
129  do
130  {
131  // --- Store previous wArA
132  wArAold = wArA;
133 
134  // --- Precondition residual
135  preconPtr->precondition(wA, rA, cmpt);
136 
137  // --- Update search directions:
138  wArA = gSumProd(wA, rA, matrix().mesh().comm());
139 
140  if (solverPerf.nIterations() == 0)
141  {
142  for (label cell=0; cell<nCells; cell++)
143  {
144  pAPtr[cell] = wAPtr[cell];
145  }
146  }
147  else
148  {
149  scalar beta = wArA/wArAold;
150 
151  for (label cell=0; cell<nCells; cell++)
152  {
153  pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
154  }
155  }
156 
157 
158  // --- Update preconditioned residual
159  matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);
160 
161  scalar wApA = gSumProd(wA, pA, matrix().mesh().comm());
162 
163 
164  // --- Test for singularity
165  if (solverPerf.checkSingularity(mag(wApA)/normFactor)) break;
166 
167 
168  // --- Update solution and residual:
169 
170  scalar alpha = wArA/wApA;
171 
172  for (label cell=0; cell<nCells; cell++)
173  {
174  psiPtr[cell] += alpha*pAPtr[cell];
175  rAPtr[cell] -= alpha*wAPtr[cell];
176  }
177 
178  solverPerf.finalResidual() =
179  gSumMag(rA, matrix().mesh().comm())
180  /normFactor;
181 
182  } while
183  (
184  (
185  solverPerf.nIterations()++ < maxIter_
186  && !solverPerf.checkConvergence(tolerance_, relTol_)
187  )
188  || solverPerf.nIterations() < minIter_
189  );
190  }
191 
192  return solverPerf;
193 }
194 
195 
196 // ************************************************************************* //
scalar gSumMag(const FieldField< Field, Type > &f)
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
uint8_t direction
Definition: direction.H:46
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
const Type & finalResidual() const
Return final residual.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
lduMatrix::solver::addsymMatrixConstructorToTable< PCG > addPCGSymMatrixConstructorToTable_
Definition: PCG.C:35
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
const lduMesh & mesh() const
Return the LDU mesh from which the addressing is obtained.
Definition: lduMatrix.H:541
bool checkSingularity(const Type &residual)
Singularity test.
Generic field type.
Definition: FieldField.H:51
static const scalar great_
Large Type for the use in solvers.
bool checkConvergence(const Type &tolerance, const Type &relTolerance)
Check, store and return convergence.
A class for handling words, derived from string.
Definition: word.H:59
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics...
static autoPtr< preconditioner > New(const solver &sol, const dictionary &solverControls)
Return a new preconditioner.
Abstract base-class for lduMatrix solvers.
Definition: lduMatrix.H:93
label nIterations() const
Return number of iterations.
iterator begin()
Return an iterator to begin traversing the UList.
Definition: UListI.H:216
virtual solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.
Definition: PCG.C:66
scalar gSumProd(const UList< Type > &f1, const UList< Type > &f2, const label comm)
defineTypeNameAndDebug(combustionModel, 0)
lduMatrix is a general matrix class in which the coefficients are stored as three arrays...
Definition: lduMatrix.H:79
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
const Type & initialResidual() const
Return initial residual.
messageStream Info
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
dimensioned< scalar > mag(const dimensioned< Type > &)
static word getName(const dictionary &)
Find the preconditioner name (directly or from a sub-dictionary)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:53
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Namespace for OpenFOAM.