PBiCG.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 "PBiCG.H"
27 
28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
29 
30 namespace Foam
31 {
32  defineTypeNameAndDebug(PBiCG, 0);
33 
34  lduMatrix::solver::addasymMatrixConstructorToTable<PBiCG>
36 }
37 
38 
39 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
40 
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 wA(nCells);
87  scalar* __restrict__ wAPtr = wA.begin();
88 
89  // --- Calculate A.psi
90  matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
91 
92  // --- Calculate initial residual field
93  scalarField rA(source - wA);
94  scalar* __restrict__ rAPtr = rA.begin();
95 
96  // --- Calculate normalisation factor
97  const scalar normFactor = this->normFactor(psi, source, wA, 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 pT(nCells, 0);
118  scalar* __restrict__ pTPtr = pT.begin();
119 
120  scalarField wT(nCells);
121  scalar* __restrict__ wTPtr = wT.begin();
122 
123  // --- Calculate T.psi
124  matrix_.Tmul(wT, psi, interfaceIntCoeffs_, interfaces_, cmpt);
125 
126  // --- Calculate initial transpose residual field
127  scalarField rT(source - wT);
128  scalar* __restrict__ rTPtr = rT.begin();
129 
130  // --- Initial value not used
131  scalar wArT = 0;
132 
133  // --- Select and construct the preconditioner
136  (
137  *this,
138  controlDict_
139  );
140 
141  // --- Solver iteration
142  do
143  {
144  // --- Store previous wArT
145  const scalar wArTold = wArT;
146 
147  // --- Precondition residuals
148  preconPtr->precondition(wA, rA, cmpt);
149  preconPtr->preconditionT(wT, rT, cmpt);
150 
151  // --- Update search directions:
152  wArT = gSumProd(wA, rT, matrix().mesh().comm());
153 
154  if (solverPerf.nIterations() == 0)
155  {
156  for (label cell=0; cell<nCells; cell++)
157  {
158  pAPtr[cell] = wAPtr[cell];
159  pTPtr[cell] = wTPtr[cell];
160  }
161  }
162  else
163  {
164  const scalar beta = wArT/wArTold;
165 
166  for (label cell=0; cell<nCells; cell++)
167  {
168  pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
169  pTPtr[cell] = wTPtr[cell] + beta*pTPtr[cell];
170  }
171  }
172 
173 
174  // --- Update preconditioned residuals
175  matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);
176  matrix_.Tmul(wT, pT, interfaceIntCoeffs_, interfaces_, cmpt);
177 
178  const scalar wApT = gSumProd(wA, pT, matrix().mesh().comm());
179 
180  // --- Test for singularity
181  if (solverPerf.checkSingularity(mag(wApT)/normFactor))
182  {
183  break;
184  }
185 
186 
187  // --- Update solution and residual:
188 
189  const scalar alpha = wArT/wApT;
190 
191  for (label cell=0; cell<nCells; cell++)
192  {
193  psiPtr[cell] += alpha*pAPtr[cell];
194  rAPtr[cell] -= alpha*wAPtr[cell];
195  rTPtr[cell] -= alpha*wTPtr[cell];
196  }
197 
198  solverPerf.finalResidual() =
199  gSumMag(rA, matrix().mesh().comm())
200  /normFactor;
201  } while
202  (
203  (
204  ++solverPerf.nIterations() < maxIter_
205  && !solverPerf.checkConvergence(tolerance_, relTol_)
206  )
207  || solverPerf.nIterations() < minIter_
208  );
209  }
210 
211  // Recommend PBiCGStab if PBiCG fails to converge
212  if (solverPerf.nIterations() > max(defaultMaxIter_, maxIter_))
213  {
215  << "PBiCG has failed to converge within the maximum number"
216  " of iterations " << max(defaultMaxIter_, maxIter_) << nl
217  << " Please try the more robust PBiCGStab solver."
218  << exit(FatalError);
219  }
220 
221  return solverPerf;
222 }
223 
224 
225 // ************************************************************************* //
scalar gSumMag(const FieldField< Field, Type > &f)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const labelType & nIterations() const
Return number of iterations.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
uint8_t direction
Definition: direction.H:45
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
PBiCG(const word &fieldName, const lduMatrix &matrix, const FieldField< Field, scalar > &interfaceBouCoeffs, const FieldField< Field, scalar > &interfaceIntCoeffs, const lduInterfaceFieldPtrsList &interfaces, const dictionary &solverControls)
Construct from matrix components and solver data stream.
Definition: PBiCG.C:42
bool checkSingularity(const Type &residual)
Singularity test.
const Type & finalResidual() const
Return final residual.
Generic field type.
Definition: FieldField.H:51
virtual solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.
Definition: PBiCG.C:66
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...
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< PBiCG > addPBiCGAsymMatrixConstructorToTable_
Definition: PBiCG.C:35
scalar gSumProd(const UList< Type > &f1, const UList< Type > &f2, const label comm)
static const char nl
Definition: Ostream.H:260
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:57
messageStream Info
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 lduMesh & mesh() const
Return the LDU mesh from which the addressing is obtained.
Definition: lduMatrix.H:544
Namespace for OpenFOAM.