PBiCG.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2019-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "PBiCG.H"
30#include "PrecisionAdaptor.H"
31
32// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33
34namespace Foam
35{
37
38 lduMatrix::solver::addasymMatrixConstructorToTable<PBiCG>
40}
41
42
43// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
44
46(
47 const word& fieldName,
48 const lduMatrix& matrix,
49 const FieldField<Field, scalar>& interfaceBouCoeffs,
50 const FieldField<Field, scalar>& interfaceIntCoeffs,
51 const lduInterfaceFieldPtrsList& interfaces,
52 const dictionary& solverControls
53)
54:
56 (
57 fieldName,
58 matrix,
59 interfaceBouCoeffs,
60 interfaceIntCoeffs,
61 interfaces,
62 solverControls
63 )
64{}
65
66
67// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
68
70(
71 scalarField& psi_s,
72 const scalarField& source,
73 const direction cmpt
74) const
75{
77 solveScalarField& psi = tpsi.ref();
78
79 // --- Setup class containing solver performance data
80 solverPerformance solverPerf
81 (
82 lduMatrix::preconditioner::getName(controlDict_) + typeName,
83 fieldName_
84 );
85
86 const label nCells = psi.size();
87
88 solveScalar* __restrict__ psiPtr = psi.begin();
89
90 solveScalarField pA(nCells);
91 solveScalar* __restrict__ pAPtr = pA.begin();
92
93 solveScalarField wA(nCells);
94 solveScalar* __restrict__ wAPtr = wA.begin();
95
96 // --- Calculate A.psi
97 matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
98
99 // --- Calculate initial residual field
101 solveScalarField rA(tsource() - wA);
102 solveScalar* __restrict__ rAPtr = rA.begin();
103
104 matrix().setResidualField
105 (
107 fieldName_,
108 true
109 );
110
111 // --- Calculate normalisation factor
112 const solveScalar normFactor = this->normFactor(psi, tsource(), wA, pA);
113
114 if ((log_ >= 2) || (lduMatrix::debug >= 2))
115 {
116 Info<< " Normalisation factor = " << normFactor << endl;
117 }
118
119 // --- Calculate normalised residual norm
120 solverPerf.initialResidual() =
121 gSumMag(rA, matrix().mesh().comm())
122 /normFactor;
123 solverPerf.finalResidual() = solverPerf.initialResidual();
124
125 // --- Check convergence, solve if not converged
126 if
127 (
128 minIter_ > 0
129 || !solverPerf.checkConvergence(tolerance_, relTol_, log_)
130 )
131 {
132 solveScalarField pT(nCells, 0);
133 solveScalar* __restrict__ pTPtr = pT.begin();
134
135 solveScalarField wT(nCells);
136 solveScalar* __restrict__ wTPtr = wT.begin();
137
138 // --- Calculate T.psi
139 matrix_.Tmul(wT, psi, interfaceIntCoeffs_, interfaces_, cmpt);
140
141 // --- Calculate initial transpose residual field
142 solveScalarField rT(tsource() - wT);
143 solveScalar* __restrict__ rTPtr = rT.begin();
144
145 // --- Initial value not used
146 solveScalar wArT = 0;
147
148 // --- Select and construct the preconditioner
151 (
152 *this,
153 controlDict_
154 );
155
156 // --- Solver iteration
157 do
158 {
159 // --- Store previous wArT
160 const solveScalar wArTold = wArT;
161
162 // --- Precondition residuals
163 preconPtr->precondition(wA, rA, cmpt);
164 preconPtr->preconditionT(wT, rT, cmpt);
165
166 // --- Update search directions:
167 wArT = gSumProd(wA, rT, matrix().mesh().comm());
168
169 if (solverPerf.nIterations() == 0)
170 {
171 for (label cell=0; cell<nCells; cell++)
172 {
173 pAPtr[cell] = wAPtr[cell];
174 pTPtr[cell] = wTPtr[cell];
175 }
176 }
177 else
178 {
179 const solveScalar beta = wArT/wArTold;
180
181 for (label cell=0; cell<nCells; cell++)
182 {
183 pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
184 pTPtr[cell] = wTPtr[cell] + beta*pTPtr[cell];
185 }
186 }
187
188
189 // --- Update preconditioned residuals
190 matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);
191 matrix_.Tmul(wT, pT, interfaceIntCoeffs_, interfaces_, cmpt);
192
193 const solveScalar wApT = gSumProd(wA, pT, matrix().mesh().comm());
194
195 // --- Test for singularity
196 if (solverPerf.checkSingularity(mag(wApT)/normFactor))
197 {
198 break;
199 }
200
201
202 // --- Update solution and residual:
203
204 const solveScalar alpha = wArT/wApT;
205
206 for (label cell=0; cell<nCells; cell++)
207 {
208 psiPtr[cell] += alpha*pAPtr[cell];
209 rAPtr[cell] -= alpha*wAPtr[cell];
210 rTPtr[cell] -= alpha*wTPtr[cell];
211 }
212
213 solverPerf.finalResidual() =
214 gSumMag(rA, matrix().mesh().comm())
215 /normFactor;
216 } while
217 (
218 (
219 ++solverPerf.nIterations() < maxIter_
220 && !solverPerf.checkConvergence(tolerance_, relTol_, log_)
221 )
222 || solverPerf.nIterations() < minIter_
223 );
224 }
225
226 // Recommend PBiCGStab if PBiCG fails to converge
227 if (solverPerf.nIterations() > max(defaultMaxIter_, maxIter_))
228 {
230 << "PBiCG has failed to converge within the maximum number"
231 " of iterations " << max(defaultMaxIter_, maxIter_) << nl
232 << " Please try the more robust PBiCGStab solver."
233 << exit(FatalError);
234 }
235
236 matrix().setResidualField
237 (
239 fieldName_,
240 false
241 );
242
243 return solverPerf;
244}
245
246
247// ************************************************************************* //
A const Field/List wrapper with possible data conversion.
A field of fields is a PtrList of fields with reference counting.
Definition: FieldField.H:80
const word & getName() const
Get name.
Definition: NURBS3DCurve.H:349
Preconditioned bi-conjugate gradient solver for asymmetric lduMatrices using a run-time selectable pr...
Definition: PBiCG.H:55
A non-const Field/List wrapper with possible data conversion.
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics.
const Type & finalResidual() const noexcept
Return final residual.
const labelType & nIterations() const noexcept
Return number of iterations.
bool checkSingularity(const Type &residual)
Singularity test.
const Type & initialResidual() const noexcept
Return initial residual.
bool checkConvergence(const Type &tolerance, const Type &relTolerance, const int logLevel=0)
Check, store and return convergence.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
iterator begin() noexcept
Return an iterator to begin traversing the UList.
Definition: UListI.H:329
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
lduMatrix is a general matrix class in which the coefficients are stored as three arrays,...
Definition: lduMatrix.H:84
T & ref() const
Definition: refPtrI.H:203
Base class for solution control classes.
Definition: solver.H:54
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const volScalarField & psi
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Namespace for OpenFOAM.
scalarProduct< Type, Type >::type gSumProd(const UList< Type > &f1, const UList< Type > &f2, const label comm)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
uint8_t direction
Definition: direction.H:56
typeOfMag< Type >::type gSumMag(const FieldField< Field, Type > &f)
lduMatrix::solver::addasymMatrixConstructorToTable< PBiCG > addPBiCGAsymMatrixConstructorToTable_
Definition: PBiCG.C:39
error FatalError
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
volScalarField & alpha
CEqn solve()
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)