PBiCGStab.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) 2016-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 "PBiCGStab.H"
27 
28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
29 
30 namespace Foam
31 {
32  defineTypeNameAndDebug(PBiCGStab, 0);
33 
34  lduMatrix::solver::addsymMatrixConstructorToTable<PBiCGStab>
36 
37  lduMatrix::solver::addasymMatrixConstructorToTable<PBiCGStab>
39 }
40 
41 
42 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
43 
44 Foam::PBiCGStab::PBiCGStab
45 (
46  const word& fieldName,
47  const lduMatrix& matrix,
48  const FieldField<Field, scalar>& interfaceBouCoeffs,
49  const FieldField<Field, scalar>& interfaceIntCoeffs,
50  const lduInterfaceFieldPtrsList& interfaces,
51  const dictionary& solverControls
52 )
53 :
55  (
56  fieldName,
57  matrix,
58  interfaceBouCoeffs,
59  interfaceIntCoeffs,
60  interfaces,
61  solverControls
62  )
63 {}
64 
65 
66 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
67 
69 (
70  scalarField& psi,
71  const scalarField& source,
72  const direction cmpt
73 ) const
74 {
75  // --- Setup class containing solver performance data
76  solverPerformance solverPerf
77  (
78  lduMatrix::preconditioner::getName(controlDict_) + typeName,
79  fieldName_
80  );
81 
82  const label nCells = psi.size();
83 
84  scalar* __restrict__ psiPtr = psi.begin();
85 
86  scalarField pA(nCells);
87  scalar* __restrict__ pAPtr = pA.begin();
88 
89  scalarField yA(nCells);
90  scalar* __restrict__ yAPtr = yA.begin();
91 
92  // --- Calculate A.psi
93  matrix_.Amul(yA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
94 
95  // --- Calculate initial residual field
96  scalarField rA(source - yA);
97  scalar* __restrict__ rAPtr = rA.begin();
98 
99  // --- Calculate normalisation factor
100  const scalar normFactor = this->normFactor(psi, source, yA, 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  scalarField AyA(nCells);
121  scalar* __restrict__ AyAPtr = AyA.begin();
122 
123  scalarField sA(nCells);
124  scalar* __restrict__ sAPtr = sA.begin();
125 
126  scalarField zA(nCells);
127  scalar* __restrict__ zAPtr = zA.begin();
128 
129  scalarField tA(nCells);
130  scalar* __restrict__ tAPtr = tA.begin();
131 
132  // --- Store initial residual
133  const scalarField rA0(rA);
134 
135  // --- Initial values not used
136  scalar rA0rA = 0;
137  scalar alpha = 0;
138  scalar omega = 0;
139 
140  // --- Select and construct the preconditioner
143  (
144  *this,
145  controlDict_
146  );
147 
148  // --- Solver iteration
149  do
150  {
151  // --- Store previous rA0rA
152  const scalar rA0rAold = rA0rA;
153 
154  rA0rA = gSumProd(rA0, rA, matrix().mesh().comm());
155 
156  // --- Test for singularity
157  if (solverPerf.checkSingularity(mag(rA0rA)))
158  {
159  break;
160  }
161 
162  // --- Update pA
163  if (solverPerf.nIterations() == 0)
164  {
165  for (label cell=0; cell<nCells; cell++)
166  {
167  pAPtr[cell] = rAPtr[cell];
168  }
169  }
170  else
171  {
172  // --- Test for singularity
173  if (solverPerf.checkSingularity(mag(omega)))
174  {
175  break;
176  }
177 
178  const scalar beta = (rA0rA/rA0rAold)*(alpha/omega);
179 
180  for (label cell=0; cell<nCells; cell++)
181  {
182  pAPtr[cell] =
183  rAPtr[cell] + beta*(pAPtr[cell] - omega*AyAPtr[cell]);
184  }
185  }
186 
187  // --- Precondition pA
188  preconPtr->precondition(yA, pA, cmpt);
189 
190  // --- Calculate AyA
191  matrix_.Amul(AyA, yA, interfaceBouCoeffs_, interfaces_, cmpt);
192 
193  const scalar rA0AyA = gSumProd(rA0, AyA, matrix().mesh().comm());
194 
195  alpha = rA0rA/rA0AyA;
196 
197  // --- Calculate sA
198  for (label cell=0; cell<nCells; cell++)
199  {
200  sAPtr[cell] = rAPtr[cell] - alpha*AyAPtr[cell];
201  }
202 
203  // --- Test sA for convergence
204  solverPerf.finalResidual() =
205  gSumMag(sA, matrix().mesh().comm())/normFactor;
206 
207  if (solverPerf.checkConvergence(tolerance_, relTol_))
208  {
209  for (label cell=0; cell<nCells; cell++)
210  {
211  psiPtr[cell] += alpha*yAPtr[cell];
212  }
213 
214  solverPerf.nIterations()++;
215 
216  return solverPerf;
217  }
218 
219  // --- Precondition sA
220  preconPtr->precondition(zA, sA, cmpt);
221 
222  // --- Calculate tA
223  matrix_.Amul(tA, zA, interfaceBouCoeffs_, interfaces_, cmpt);
224 
225  const scalar tAtA = gSumSqr(tA, matrix().mesh().comm());
226 
227  // --- Calculate omega from tA and sA
228  // (cheaper than using zA with preconditioned tA)
229  omega = gSumProd(tA, sA, matrix().mesh().comm())/tAtA;
230 
231  // --- Update solution and residual
232  for (label cell=0; cell<nCells; cell++)
233  {
234  psiPtr[cell] += alpha*yAPtr[cell] + omega*zAPtr[cell];
235  rAPtr[cell] = sAPtr[cell] - omega*tAPtr[cell];
236  }
237 
238  solverPerf.finalResidual() =
239  gSumMag(rA, matrix().mesh().comm())
240  /normFactor;
241  } while
242  (
243  (
244  ++solverPerf.nIterations() < maxIter_
245  && !solverPerf.checkConvergence(tolerance_, relTol_)
246  )
247  || solverPerf.nIterations() < minIter_
248  );
249  }
250 
251  return solverPerf;
252 }
253 
254 
255 // ************************************************************************* //
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
const labelType & nIterations() const
Return number of iterations.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
uint8_t direction
Definition: direction.H:45
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
bool checkSingularity(const Type &residual)
Singularity test.
const Type & finalResidual() const
Return final residual.
lduMatrix::solver::addsymMatrixConstructorToTable< PBiCGStab > addPBiCGStabSymMatrixConstructorToTable_
Definition: PBiCGStab.C:35
Generic field type.
Definition: FieldField.H:51
bool checkConvergence(const Type &tolerance, const Type &relTolerance)
Check, store and return convergence.
const Type & initialResidual() const
Return initial residual.
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
iterator begin()
Return an iterator to begin traversing the UList.
Definition: UListI.H:216
lduMatrix::solver::addasymMatrixConstructorToTable< PBiCGStab > addPBiCGStabAsymMatrixConstructorToTable_
Definition: PBiCGStab.C:38
scalar gSumProd(const UList< Type > &f1, const UList< Type > &f2, const label comm)
virtual solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.
Definition: PBiCGStab.C:69
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
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:52
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
const lduMesh & mesh() const
Return the LDU mesh from which the addressing is obtained.
Definition: lduMatrix.H:544
Namespace for OpenFOAM.