PCG.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 "PCG.H"
30#include "PrecisionAdaptor.H"
31
32// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33
34namespace Foam
35{
37
38 lduMatrix::solver::addsymMatrixConstructorToTable<PCG>
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(
72 const solveScalarField& source,
73 const direction cmpt
74) const
75{
76 // --- Setup class containing solver performance data
77 solverPerformance solverPerf
78 (
79 lduMatrix::preconditioner::getName(controlDict_) + typeName,
80 fieldName_
81 );
82
83 label nCells = psi.size();
84
85 solveScalar* __restrict__ psiPtr = psi.begin();
86
87 solveScalarField pA(nCells);
88 solveScalar* __restrict__ pAPtr = pA.begin();
89
90 solveScalarField wA(nCells);
91 solveScalar* __restrict__ wAPtr = wA.begin();
92
93 solveScalar wArA = solverPerf.great_;
94 solveScalar wArAold = wArA;
95
96 // --- Calculate A.psi
97 matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
98
99 // --- Calculate initial residual field
100 solveScalarField rA(source - wA);
101 solveScalar* __restrict__ rAPtr = rA.begin();
102
103 matrix().setResidualField
104 (
106 fieldName_,
107 true
108 );
109
110 // --- Calculate normalisation factor
111 solveScalar normFactor = this->normFactor(psi, source, wA, 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 // --- Select and construct the preconditioner
134 (
135 *this,
136 controlDict_
137 );
138
139 // --- Solver iteration
140 do
141 {
142 // --- Store previous wArA
143 wArAold = wArA;
144
145 // --- Precondition residual
146 preconPtr->precondition(wA, rA, cmpt);
147
148 // --- Update search directions:
149 wArA = gSumProd(wA, rA, matrix().mesh().comm());
150
151 if (solverPerf.nIterations() == 0)
152 {
153 for (label cell=0; cell<nCells; cell++)
154 {
155 pAPtr[cell] = wAPtr[cell];
156 }
157 }
158 else
159 {
160 solveScalar beta = wArA/wArAold;
161
162 for (label cell=0; cell<nCells; cell++)
163 {
164 pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
165 }
166 }
167
168
169 // --- Update preconditioned residual
170 matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);
171
172 solveScalar wApA = gSumProd(wA, pA, matrix().mesh().comm());
173
174 // --- Test for singularity
175 if (solverPerf.checkSingularity(mag(wApA)/normFactor)) break;
176
177
178 // --- Update solution and residual:
179
180 solveScalar alpha = wArA/wApA;
181
182 for (label cell=0; cell<nCells; cell++)
183 {
184 psiPtr[cell] += alpha*pAPtr[cell];
185 rAPtr[cell] -= alpha*wAPtr[cell];
186 }
187
188 solverPerf.finalResidual() =
189 gSumMag(rA, matrix().mesh().comm())
190 /normFactor;
191
192 } while
193 (
194 (
195 ++solverPerf.nIterations() < maxIter_
196 && !solverPerf.checkConvergence(tolerance_, relTol_, log_)
197 )
198 || solverPerf.nIterations() < minIter_
199 );
200 }
201
202 matrix().setResidualField
203 (
205 fieldName_,
206 false
207 );
208
209 return solverPerf;
210}
211
212
213
215(
216 scalarField& psi_s,
217 const scalarField& source,
218 const direction cmpt
219) const
220{
222 return scalarSolve
223 (
224 tpsi.ref(),
226 cmpt
227 );
228}
229
230
231// ************************************************************************* //
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 conjugate gradient solver for symmetric lduMatrices using a run-time selectable precon...
Definition: PCG.H:58
virtual solverPerformance scalarSolve(solveScalarField &psi, const solveScalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.
Definition: PCG.C:70
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.
static const scalar great_
Large Type for the use in solvers.
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)
uint8_t direction
Definition: direction.H:56
typeOfMag< Type >::type gSumMag(const FieldField< Field, Type > &f)
lduMatrix::solver::addsymMatrixConstructorToTable< PCG > addPCGSymMatrixConstructorToTable_
Definition: PCG.C:39
volScalarField & alpha
CEqn solve()
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)