PBiCICG.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-2016 OpenFOAM Foundation
9 Copyright (C) 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 "PBiCICG.H"
30
31// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32
33template<class Type, class DType, class LUType>
35(
36 const word& fieldName,
38 const dictionary& solverDict
39)
40:
41 LduMatrix<Type, DType, LUType>::solver
42 (
43 fieldName,
44 matrix,
45 solverDict
46 )
47{}
48
49
50// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
51
52template<class Type, class DType, class LUType>
55{
56 const word preconditionerName(this->controlDict_.getWord("preconditioner"));
57
58 // --- Setup class containing solver performance data
60 (
61 preconditionerName + typeName,
62 this->fieldName_
63 );
64
65 label nIter = 0;
66
67 label nCells = psi.size();
68
69 Type* __restrict__ psiPtr = psi.begin();
70
71 Field<Type> pA(nCells);
72 Type* __restrict__ pAPtr = pA.begin();
73
74 Field<Type> pT(nCells, Zero);
75 Type* __restrict__ pTPtr = pT.begin();
76
77 Field<Type> wA(nCells);
78 Type* __restrict__ wAPtr = wA.begin();
79
80 Field<Type> wT(nCells);
81 Type* __restrict__ wTPtr = wT.begin();
82
83 Type wArT = solverPerf.great_*pTraits<Type>::one;
84 Type wArTold = wArT;
85
86 // --- Calculate A.psi and T.psi
87 this->matrix_.Amul(wA, psi);
88 this->matrix_.Tmul(wT, psi);
89
90 // --- Calculate initial residual and transpose residual fields
91 Field<Type> rA(this->matrix_.source() - wA);
92 Field<Type> rT(this->matrix_.source() - wT);
93 Type* __restrict__ rAPtr = rA.begin();
94 Type* __restrict__ rTPtr = rT.begin();
95
96 // --- Calculate normalisation factor
97 Type normFactor = this->normFactor(psi, wA, pA);
98
99 if ((this->log_ >= 2) || (LduMatrix<Type, DType, LUType>::debug >= 2))
100 {
101 Info<< " Normalisation factor = " << normFactor << endl;
102 }
103
104 // --- Calculate normalised residual norm
105 solverPerf.initialResidual() = cmptDivide(gSumCmptMag(rA), normFactor);
106 solverPerf.finalResidual() = solverPerf.initialResidual();
107
108 // --- Check convergence, solve if not converged
109 if
110 (
111 !solverPerf.checkConvergence
112 (
113 this->tolerance_,
114 this->relTol_,
115 this->log_
116 )
117 )
118 {
119 // --- Select and construct the preconditioner
122 (
123 *this,
124 this->controlDict_
125 );
126
127 // --- Solver iteration
128 do
129 {
130 // --- Store previous wArT
131 wArTold = wArT;
132
133 // --- Precondition residuals
134 preconPtr->precondition(wA, rA);
135 preconPtr->preconditionT(wT, rT);
136
137 // --- Update search directions:
138 wArT = gSumCmptProd(wA, rT);
139
140 if (nIter == 0)
141 {
142 for (label cell=0; cell<nCells; cell++)
143 {
144 pAPtr[cell] = wAPtr[cell];
145 pTPtr[cell] = wTPtr[cell];
146 }
147 }
148 else
149 {
150 Type beta = cmptDivide
151 (
152 wArT,
153 stabilise(wArTold, solverPerf.vsmall_)
154 );
155
156 for (label cell=0; cell<nCells; cell++)
157 {
158 pAPtr[cell] = wAPtr[cell] + cmptMultiply(beta, pAPtr[cell]);
159 pTPtr[cell] = wTPtr[cell] + cmptMultiply(beta, pTPtr[cell]);
160 }
161 }
162
163
164 // --- Update preconditioned residuals
165 this->matrix_.Amul(wA, pA);
166 this->matrix_.Tmul(wT, pT);
167
168 Type wApT = gSumCmptProd(wA, pT);
169
170 // --- Test for singularity
171 if
172 (
173 solverPerf.checkSingularity
174 (
175 cmptDivide(cmptMag(wApT), normFactor)
176 )
177 )
178 {
179 break;
180 }
181
182
183 // --- Update solution and residual:
184
185 Type alpha = cmptDivide
186 (
187 wArT,
188 stabilise(wApT, solverPerf.vsmall_)
189 );
190
191 for (label cell=0; cell<nCells; cell++)
192 {
193 psiPtr[cell] += cmptMultiply(alpha, pAPtr[cell]);
194 rAPtr[cell] -= cmptMultiply(alpha, wAPtr[cell]);
195 rTPtr[cell] -= cmptMultiply(alpha, wTPtr[cell]);
196 }
197
198 solverPerf.finalResidual() =
199 cmptDivide(gSumCmptMag(rA), normFactor);
200
201 } while
202 (
203 nIter++ < this->maxIter_
204 && !solverPerf.checkConvergence
205 (
206 this->tolerance_,
207 this->relTol_,
208 this->log_
209 )
210 );
211 }
212
213 solverPerf.nIterations() =
215
216 return solverPerf;
217}
218
219
220// ************************************************************************* //
Generic templated field type.
Definition: Field.H:82
LduMatrix is a general matrix class in which the coefficients are stored as three arrays,...
Definition: LduMatrix.H:88
Preconditioned bi-conjugate gradient solver for asymmetric lduMatrices using a run-time selectable pr...
Definition: PBiCICG.H:56
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 vsmall_
Very small Type for the use in solvers.
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
A class representing the concept of 1 (one) that can be used to avoid manipulating objects known to b...
Definition: one.H:62
A traits class, which is primarily used for primitives.
Definition: pTraits.H:59
Base class for solution control classes.
Definition: solver.H:54
A class for handling words, derived from Foam::string.
Definition: word.H:68
const volScalarField & psi
dimensioned< Type > cmptDivide(const dimensioned< Type > &, const dimensioned< Type > &)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensionedScalar stabilise(const dimensionedScalar &x, const dimensionedScalar &y)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
Type gSumCmptMag(const UList< Type > &f, const label comm)
Type gSumCmptProd(const UList< Type > &f1, const UList< Type > &f2, const label comm)
volScalarField & alpha
CEqn solve()
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)