PBiCGStab.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) 2016 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 "PBiCGStab.H"
27 
28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
29 
30 namespace Foam
31 {
32  defineTypeNameAndDebug(PBiCGStab, 0);
33 
34  lduMatrix::solver::addasymMatrixConstructorToTable<PBiCGStab>
36 }
37 
38 
39 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
40 
41 Foam::PBiCGStab::PBiCGStab
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  const 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 yA(nCells);
87  scalar* __restrict__ yAPtr = yA.begin();
88 
89  // --- Calculate A.psi
90  matrix_.Amul(yA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
91 
92  // --- Calculate initial residual field
93  scalarField rA(source - yA);
94  scalar* __restrict__ rAPtr = rA.begin();
95 
96  // --- Calculate normalisation factor
97  const scalar normFactor = this->normFactor(psi, source, yA, pA);
98 
99  if (lduMatrix::debug >= 2)
100  {
101  Info<< " Normalisation factor = " << normFactor << endl;
102  }
103 
104  // --- Calculate normalised residual norm
105  solverPerf.initialResidual() =
106  gSumMag(rA, matrix().mesh().comm())
107  /normFactor;
108  solverPerf.finalResidual() = solverPerf.initialResidual();
109 
110  // --- Check convergence, solve if not converged
111  if
112  (
113  minIter_ > 0
114  || !solverPerf.checkConvergence(tolerance_, relTol_)
115  )
116  {
117  scalarField AyA(nCells);
118  scalar* __restrict__ AyAPtr = AyA.begin();
119 
120  scalarField sA(nCells);
121  scalar* __restrict__ sAPtr = sA.begin();
122 
123  scalarField zA(nCells);
124  scalar* __restrict__ zAPtr = zA.begin();
125 
126  scalarField tA(nCells);
127  scalar* __restrict__ tAPtr = tA.begin();
128 
129  // --- Store initial residual
130  const scalarField rA0(rA);
131 
132  // --- Initial values not used
133  scalar rA0rA = 0;
134  scalar alpha = 0;
135  scalar omega = 0;
136 
137  // --- Select and construct the preconditioner
140  (
141  *this,
142  controlDict_
143  );
144 
145  // --- Solver iteration
146  do
147  {
148  // --- Store previous rA0rA
149  const scalar rA0rAold = rA0rA;
150 
151  rA0rA = gSumProd(rA0, rA, matrix().mesh().comm());
152 
153  // --- Test for singularity
154  if (solverPerf.checkSingularity(mag(rA0rA)))
155  {
156  break;
157  }
158 
159  // --- Update pA
160  if (solverPerf.nIterations() == 0)
161  {
162  for (label cell=0; cell<nCells; cell++)
163  {
164  pAPtr[cell] = rAPtr[cell];
165  }
166  }
167  else
168  {
169  // --- Test for singularity
170  if (solverPerf.checkSingularity(mag(omega)))
171  {
172  break;
173  }
174 
175  const scalar beta = (rA0rA/rA0rAold)*(alpha/omega);
176 
177  for (label cell=0; cell<nCells; cell++)
178  {
179  pAPtr[cell] =
180  rAPtr[cell] + beta*(pAPtr[cell] - omega*AyAPtr[cell]);
181  }
182  }
183 
184  // --- Precondition pA
185  preconPtr->precondition(yA, pA, cmpt);
186 
187  // --- Calculate AyA
188  matrix_.Amul(AyA, yA, interfaceBouCoeffs_, interfaces_, cmpt);
189 
190  const scalar rA0AyA = gSumProd(rA0, AyA, matrix().mesh().comm());
191 
192  alpha = rA0rA/rA0AyA;
193 
194  // --- Calculate sA
195  for (label cell=0; cell<nCells; cell++)
196  {
197  sAPtr[cell] = rAPtr[cell] - alpha*AyAPtr[cell];
198  }
199 
200  // --- Test sA for convergence
201  solverPerf.finalResidual() =
202  gSumMag(sA, matrix().mesh().comm())/normFactor;
203 
204  if (solverPerf.checkConvergence(tolerance_, relTol_))
205  {
206  for (label cell=0; cell<nCells; cell++)
207  {
208  psiPtr[cell] += alpha*yAPtr[cell];
209  }
210 
211  solverPerf.nIterations()++;
212 
213  return solverPerf;
214  }
215 
216  // --- Precondition sA
217  preconPtr->precondition(zA, sA, cmpt);
218 
219  // --- Calculate tA
220  matrix_.Amul(tA, zA, interfaceBouCoeffs_, interfaces_, cmpt);
221 
222  const scalar tAtA = gSumSqr(tA, matrix().mesh().comm());
223 
224  // --- Calculate omega from tA and sA
225  // (cheaper than using zA with preconditioned tA)
226  omega = gSumProd(tA, sA)/tAtA;
227 
228  // --- Update solution and residual
229  for (label cell=0; cell<nCells; cell++)
230  {
231  psiPtr[cell] += alpha*yAPtr[cell] + omega*zAPtr[cell];
232  rAPtr[cell] = sAPtr[cell] - omega*tAPtr[cell];
233  }
234 
235  solverPerf.finalResidual() =
236  gSumMag(rA, matrix().mesh().comm())
237  /normFactor;
238  } while
239  (
240  (
241  solverPerf.nIterations()++ < maxIter_
242  && !solverPerf.checkConvergence(tolerance_, relTol_)
243  )
244  || solverPerf.nIterations() < minIter_
245  );
246  }
247 
248  return solverPerf;
249 }
250 
251 
252 // ************************************************************************* //
virtual solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.
Definition: PBiCGStab.C:66
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
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
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...
scalar gSumSqr(const UList< Type > &f, const label comm)
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
lduMatrix::solver::addasymMatrixConstructorToTable< PBiCGStab > addPBiCGStabAsymMatrixConstructorToTable_
Definition: PBiCGStab.C:35
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.