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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2016-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 "PBiCGStab.H"
30#include "PrecisionAdaptor.H"
31
32// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33
34namespace Foam
35{
37
38 lduMatrix::solver::addsymMatrixConstructorToTable<PBiCGStab>
40
41 lduMatrix::solver::addasymMatrixConstructorToTable<PBiCGStab>
43}
44
45
46// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
47
49(
50 const word& fieldName,
51 const lduMatrix& matrix,
52 const FieldField<Field, scalar>& interfaceBouCoeffs,
53 const FieldField<Field, scalar>& interfaceIntCoeffs,
54 const lduInterfaceFieldPtrsList& interfaces,
55 const dictionary& solverControls
56)
57:
59 (
60 fieldName,
61 matrix,
62 interfaceBouCoeffs,
63 interfaceIntCoeffs,
64 interfaces,
65 solverControls
66 )
67{}
68
69
70// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
71
73(
75 const solveScalarField& source,
76 const direction cmpt
77) const
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 yA(nCells);
94 solveScalar* __restrict__ yAPtr = yA.begin();
95
96 // --- Calculate A.psi
97 matrix_.Amul(yA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
98
99 // --- Calculate initial residual field
100 solveScalarField rA(source - yA);
101 solveScalar* __restrict__ rAPtr = rA.begin();
102
103 matrix().setResidualField
104 (
106 fieldName_,
107 true
108 );
109
110 // --- Calculate normalisation factor
111 const solveScalar normFactor = this->normFactor(psi, source, yA, pA);
112
113 if ((log_ >= 2) || (lduMatrix::debug >= 2))
114 {
115 Info<< " Normalisation factor = " << normFactor << endl;
116 }
117
118 // --- Calculate normalised residual norm
119 solverPerf.initialResidual() =
120 gSumMag(rA, matrix().mesh().comm())
121 /normFactor;
122 solverPerf.finalResidual() = solverPerf.initialResidual();
123
124 // --- Check convergence, solve if not converged
125 if
126 (
127 minIter_ > 0
128 || !solverPerf.checkConvergence(tolerance_, relTol_, log_)
129 )
130 {
131 solveScalarField AyA(nCells);
132 solveScalar* __restrict__ AyAPtr = AyA.begin();
133
134 solveScalarField sA(nCells);
135 solveScalar* __restrict__ sAPtr = sA.begin();
136
137 solveScalarField zA(nCells);
138 solveScalar* __restrict__ zAPtr = zA.begin();
139
140 solveScalarField tA(nCells);
141 solveScalar* __restrict__ tAPtr = tA.begin();
142
143 // --- Store initial residual
144 const solveScalarField rA0(rA);
145
146 // --- Initial values not used
147 solveScalar rA0rA = 0;
148 solveScalar alpha = 0;
149 solveScalar omega = 0;
150
151 // --- Select and construct the preconditioner
154 (
155 *this,
156 controlDict_
157 );
158
159 // --- Solver iteration
160 do
161 {
162 // --- Store previous rA0rA
163 const solveScalar rA0rAold = rA0rA;
164
165 rA0rA = gSumProd(rA0, rA, matrix().mesh().comm());
166
167 // --- Test for singularity
168 if (solverPerf.checkSingularity(mag(rA0rA)))
169 {
170 break;
171 }
172
173 // --- Update pA
174 if (solverPerf.nIterations() == 0)
175 {
176 for (label cell=0; cell<nCells; cell++)
177 {
178 pAPtr[cell] = rAPtr[cell];
179 }
180 }
181 else
182 {
183 // --- Test for singularity
184 if (solverPerf.checkSingularity(mag(omega)))
185 {
186 break;
187 }
188
189 const solveScalar beta = (rA0rA/rA0rAold)*(alpha/omega);
190
191 for (label cell=0; cell<nCells; cell++)
192 {
193 pAPtr[cell] =
194 rAPtr[cell] + beta*(pAPtr[cell] - omega*AyAPtr[cell]);
195 }
196 }
197
198 // --- Precondition pA
199 preconPtr->precondition(yA, pA, cmpt);
200
201 // --- Calculate AyA
202 matrix_.Amul(AyA, yA, interfaceBouCoeffs_, interfaces_, cmpt);
203
204 const solveScalar rA0AyA =
205 gSumProd(rA0, AyA, matrix().mesh().comm());
206
207 alpha = rA0rA/rA0AyA;
208
209 // --- Calculate sA
210 for (label cell=0; cell<nCells; cell++)
211 {
212 sAPtr[cell] = rAPtr[cell] - alpha*AyAPtr[cell];
213 }
214
215 // --- Test sA for convergence
216 solverPerf.finalResidual() =
217 gSumMag(sA, matrix().mesh().comm())/normFactor;
218
219 if
220 (
221 solverPerf.nIterations() >= minIter_
222 && solverPerf.checkConvergence(tolerance_, relTol_, log_)
223 )
224 {
225 for (label cell=0; cell<nCells; cell++)
226 {
227 psiPtr[cell] += alpha*yAPtr[cell];
228 }
229
230 solverPerf.nIterations()++;
231
232 return solverPerf;
233 }
234
235 // --- Precondition sA
236 preconPtr->precondition(zA, sA, cmpt);
237
238 // --- Calculate tA
239 matrix_.Amul(tA, zA, interfaceBouCoeffs_, interfaces_, cmpt);
240
241 const solveScalar tAtA = gSumSqr(tA, matrix().mesh().comm());
242
243 // --- Calculate omega from tA and sA
244 // (cheaper than using zA with preconditioned tA)
245 omega = gSumProd(tA, sA, matrix().mesh().comm())/tAtA;
246
247 // --- Update solution and residual
248 for (label cell=0; cell<nCells; cell++)
249 {
250 psiPtr[cell] += alpha*yAPtr[cell] + omega*zAPtr[cell];
251 rAPtr[cell] = sAPtr[cell] - omega*tAPtr[cell];
252 }
253
254 solverPerf.finalResidual() =
255 gSumMag(rA, matrix().mesh().comm())
256 /normFactor;
257 } while
258 (
259 (
260 ++solverPerf.nIterations() < maxIter_
261 && !solverPerf.checkConvergence(tolerance_, relTol_, log_)
262 )
263 || solverPerf.nIterations() < minIter_
264 );
265 }
266
267 matrix().setResidualField
268 (
270 fieldName_,
271 false
272 );
273
274 return solverPerf;
275}
276
277
279(
280 scalarField& psi_s,
281 const scalarField& source,
282 const direction cmpt
283) const
284{
286 return scalarSolve
287 (
288 tpsi.ref(),
290 cmpt
291 );
292}
293
294
295// ************************************************************************* //
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 stabilized solver for asymmetric lduMatrices using a run-time se...
Definition: PBiCGStab.H:70
virtual solverPerformance scalarSolve(solveScalarField &psi, const solveScalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.
Definition: PBiCGStab.C:73
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
Namespace for OpenFOAM.
scalarProduct< Type, Type >::type gSumProd(const UList< Type > &f1, const UList< Type > &f2, const label comm)
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)
lduMatrix::solver::addasymMatrixConstructorToTable< PBiCGStab > addPBiCGStabAsymMatrixConstructorToTable_
Definition: PBiCGStab.C:42
uint8_t direction
Definition: direction.H:56
typeOfMag< Type >::type gSumMag(const FieldField< Field, Type > &f)
outerProduct1< Type >::type gSumSqr(const UList< Type > &f, const label comm)
lduMatrix::solver::addsymMatrixConstructorToTable< PBiCGStab > addPBiCGStabSymMatrixConstructorToTable_
Definition: PBiCGStab.C:39
volScalarField & alpha
CEqn solve()
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)