PPCG.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) 2019-2020 Mattijs Janssens
9 Copyright (C) 2020-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 "PPCG.H"
30#include "PrecisionAdaptor.H"
31
32// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33
34namespace Foam
35{
37
38 lduMatrix::solver::addsymMatrixConstructorToTable<PPCG>
40}
41
42
43// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44
45void Foam::PPCG::gSumMagProd
46(
47 FixedList<solveScalar, 3>& globalSum,
48 const solveScalarField& a,
49 const solveScalarField& b,
50 const solveScalarField& c,
51 const solveScalarField& sumMag,
52 label& outstandingRequest,
53 const label comm
54) const
55{
56 const label nCells = a.size();
57
58 globalSum = 0.0;
59 for (label cell=0; cell<nCells; ++cell)
60 {
61 globalSum[0] += a[cell]*b[cell]; // sumProd(a, b)
62 globalSum[1] += a[cell]*c[cell]; // sumProd(a, c)
63 globalSum[2] += mag(sumMag[cell]);
64 }
65
66 if (Pstream::parRun())
67 {
69 (
70 globalSum.data(),
71 globalSum.size(),
72 sumOp<solveScalar>(),
74 comm,
75 outstandingRequest
76 );
77 }
78}
79
80
82(
84 const solveScalarField& source,
85 const direction cmpt,
86 const bool cgMode
87) const
88{
89 // --- Setup class containing solver performance data
90 solverPerformance solverPerf
91 (
93 fieldName_
94 );
95
96 const label comm = matrix().mesh().comm();
97 const label nCells = psi.size();
98 solveScalarField w(nCells);
99
100 // --- Calculate A.psi
101 matrix_.Amul(w, psi, interfaceBouCoeffs_, interfaces_, cmpt);
102
103 // --- Calculate initial residual field
104 solveScalarField r(source - w);
105
106 // --- Calculate normalisation factor
107 solveScalarField p(nCells);
108 const solveScalar normFactor = this->normFactor(psi, source, w, p);
109
110 if ((log_ >= 2) || (lduMatrix::debug >= 2))
111 {
112 Info<< " Normalisation factor = " << normFactor << endl;
113 }
114
115 // --- Select and construct the preconditioner
118 (
119 *this,
120 controlDict_
121 );
122
123 // --- Precondition residual (= u0)
124 solveScalarField u(nCells);
125 preconPtr->precondition(u, r, cmpt);
126
127 // --- Calculate A*u - reuse w
128 matrix_.Amul(w, u, interfaceBouCoeffs_, interfaces_, cmpt);
129
130
131 // State
132 solveScalarField s(nCells);
133 solveScalarField q(nCells);
134 solveScalarField z(nCells);
135
136 solveScalarField m(nCells);
137
139 label outstandingRequest = -1;
140 if (cgMode)
141 {
142 // --- Start global reductions for inner products
143 gSumMagProd(globalSum, u, r, w, r, outstandingRequest, comm);
144
145 // --- Precondition residual
146 preconPtr->precondition(m, w, cmpt);
147 }
148 else
149 {
150 // --- Precondition residual
151 preconPtr->precondition(m, w, cmpt);
152
153 // --- Start global reductions for inner products
154 gSumMagProd(globalSum, w, u, m, r, outstandingRequest, comm);
155 }
156
157 // --- Calculate A*m
158 solveScalarField n(nCells);
159 matrix_.Amul(n, m, interfaceBouCoeffs_, interfaces_, cmpt);
160
161 solveScalar alpha = 0.0;
162 solveScalar gamma = 0.0;
163
164 // --- Solver iteration
165 for
166 (
167 solverPerf.nIterations() = 0;
168 solverPerf.nIterations() < maxIter_;
169 solverPerf.nIterations()++
170 )
171 {
172 // Make sure gamma,delta are available
173 if (Pstream::parRun())
174 {
175 Pstream::waitRequest(outstandingRequest);
176 outstandingRequest = -1;
177 }
178
179 const solveScalar gammaOld = gamma;
180 gamma = globalSum[0];
181 const solveScalar delta = globalSum[1];
182
183 solverPerf.finalResidual() = globalSum[2]/normFactor;
184 if (solverPerf.nIterations() == 0)
185 {
186 solverPerf.initialResidual() = solverPerf.finalResidual();
187 }
188
189 // Check convergence (bypass if not enough iterations yet)
190 if
191 (
192 (minIter_ <= 0 || solverPerf.nIterations() >= minIter_)
193 && solverPerf.checkConvergence(tolerance_, relTol_, log_)
194 )
195 {
196 break;
197 }
198
199
200 if (solverPerf.nIterations() == 0)
201 {
202 alpha = gamma/delta;
203 z = n;
204 q = m;
205 s = w;
206 p = u;
207 }
208 else
209 {
210 const solveScalar beta = gamma/gammaOld;
212
213 for (label cell=0; cell<nCells; ++cell)
214 {
215 z[cell] = n[cell] + beta*z[cell];
216 q[cell] = m[cell] + beta*q[cell];
217 s[cell] = w[cell] + beta*s[cell];
218 p[cell] = u[cell] + beta*p[cell];
219 }
220 }
221
222 for (label cell=0; cell<nCells; ++cell)
223 {
224 psi[cell] += alpha*p[cell];
225 r[cell] -= alpha*s[cell];
226 u[cell] -= alpha*q[cell];
227 w[cell] -= alpha*z[cell];
228 }
229
230 if (cgMode)
231 {
232 // --- Start global reductions for inner products
233 gSumMagProd(globalSum, u, r, w, r, outstandingRequest, comm);
234
235 // --- Precondition residual
236 preconPtr->precondition(m, w, cmpt);
237 }
238 else
239 {
240 // --- Precondition residual
241 preconPtr->precondition(m, w, cmpt);
242
243 // --- Start global reductions for inner products
244 gSumMagProd(globalSum, w, u, m, r, outstandingRequest, comm);
245 }
246
247 // --- Calculate A*m
248 matrix_.Amul(n, m, interfaceBouCoeffs_, interfaces_, cmpt);
249 }
250
251 return solverPerf;
252}
253
254
255
256// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
257
259(
260 const word& fieldName,
261 const lduMatrix& matrix,
262 const FieldField<Field, scalar>& interfaceBouCoeffs,
263 const FieldField<Field, scalar>& interfaceIntCoeffs,
264 const lduInterfaceFieldPtrsList& interfaces,
265 const dictionary& solverControls
266)
267:
269 (
270 fieldName,
271 matrix,
272 interfaceBouCoeffs,
273 interfaceIntCoeffs,
274 interfaces,
275 solverControls
276 )
277{}
278
279
280// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
281
283(
284 scalarField& psi_s,
285 const scalarField& source,
286 const direction cmpt
287) const
288{
290 return scalarSolveCG
291 (
292 tpsi.ref(),
294 cmpt,
295 true // operate in conjugate-gradient mode
296 );
297}
298
299
300// ************************************************************************* //
scalar delta
label n
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
A 1D vector of objects of type <T> with a fixed length <N>.
Definition: FixedList.H:81
const word & getName() const
Get name.
Definition: NURBS3DCurve.H:349
Preconditioned pipelined conjugate gradient solver for symmetric lduMatrices using a run-time selecta...
Definition: PPCG.H:67
solverPerformance scalarSolveCG(solveScalarField &psi, const solveScalarField &source, const direction cmpt, const bool cgMode) const
Definition: PPCG.C:82
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.
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
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
static int & msgType() noexcept
Message tag of standard messages.
Definition: UPstream.H:556
static void waitRequest(const label i)
Wait until request i has finished.
Definition: UPstream.C:104
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
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
volScalarField & p
const volScalarField & psi
const scalar gamma
Definition: EEqn.H:9
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Namespace for OpenFOAM.
messageStream Info
Information stream (stdout output on master, null elsewhere)
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: MSwindows.C:598
dimensioned< typename typeOfMag< Type >::type > sumMag(const DimensionedField< Type, GeoMesh > &df)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
uint8_t direction
Definition: direction.H:56
lduMatrix::solver::addsymMatrixConstructorToTable< PPCG > addPPCGSymMatrixConstructorToTable_
Definition: PPCG.C:39
volScalarField & alpha
CEqn solve()
volScalarField & b
Definition: createFields.H:27
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)