FreeFOAM The Cross-Platform CFD Toolkit
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 | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include <OpenFOAM/PBiCG.H>
27 
28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
29 
30 namespace Foam
31 {
32  defineTypeNameAndDebug(PBiCG, 0);
33 
34  lduMatrix::solver::addasymMatrixConstructorToTable<PBiCG>
36 }
37 
38 
39 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
40 
41 Foam::PBiCG::PBiCG
42 (
43  const word& fieldName,
44  const lduMatrix& matrix,
45  const FieldField<Field, scalar>& interfaceBouCoeffs,
46  const FieldField<Field, scalar>& interfaceIntCoeffs,
47  const lduInterfaceFieldPtrsList& interfaces,
48  const dictionary& solverControls
49 )
50 :
52  (
53  fieldName,
54  matrix,
55  interfaceBouCoeffs,
56  interfaceIntCoeffs,
57  interfaces,
58  solverControls
59  )
60 {}
61 
62 
63 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
64 
66 (
68  const scalarField& source,
69  const direction cmpt
70 ) const
71 {
72  // --- Setup class containing solver performance data
74  (
75  lduMatrix::preconditioner::getName(controlDict_) + typeName,
76  fieldName_
77  );
78 
79  register label nCells = psi.size();
80 
81  scalar* __restrict__ psiPtr = psi.begin();
82 
83  scalarField pA(nCells);
84  scalar* __restrict__ pAPtr = pA.begin();
85 
86  scalarField pT(nCells, 0.0);
87  scalar* __restrict__ pTPtr = pT.begin();
88 
89  scalarField wA(nCells);
90  scalar* __restrict__ wAPtr = wA.begin();
91 
92  scalarField wT(nCells);
93  scalar* __restrict__ wTPtr = wT.begin();
94 
95  scalar wArT = matrix_.great_;
96  scalar wArTold = wArT;
97 
98  // --- Calculate A.psi and T.psi
99  matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
100  matrix_.Tmul(wT, psi, interfaceIntCoeffs_, interfaces_, cmpt);
101 
102  // --- Calculate initial residual and transpose residual fields
103  scalarField rA(source - wA);
104  scalarField rT(source - wT);
105  scalar* __restrict__ rAPtr = rA.begin();
106  scalar* __restrict__ rTPtr = rT.begin();
107 
108  // --- Calculate normalisation factor
109  scalar normFactor = this->normFactor(psi, source, wA, pA);
110 
111  if (lduMatrix::debug >= 2)
112  {
113  Info<< " Normalisation factor = " << normFactor << endl;
114  }
115 
116  // --- Calculate normalised residual norm
117  solverPerf.initialResidual() = gSumMag(rA)/normFactor;
118  solverPerf.finalResidual() = solverPerf.initialResidual();
119 
120  // --- Check convergence, solve if not converged
121  if (!solverPerf.checkConvergence(tolerance_, relTol_))
122  {
123  // --- Select and construct the preconditioner
125  lduMatrix::preconditioner::New
126  (
127  *this,
128  controlDict_
129  );
130 
131  // --- Solver iteration
132  do
133  {
134  // --- Store previous wArT
135  wArTold = wArT;
136 
137  // --- Precondition residuals
138  preconPtr->precondition(wA, rA, cmpt);
139  preconPtr->preconditionT(wT, rT, cmpt);
140 
141  // --- Update search directions:
142  wArT = gSumProd(wA, rT);
143 
144  if (solverPerf.nIterations() == 0)
145  {
146  for (register label cell=0; cell<nCells; cell++)
147  {
148  pAPtr[cell] = wAPtr[cell];
149  pTPtr[cell] = wTPtr[cell];
150  }
151  }
152  else
153  {
154  scalar beta = wArT/wArTold;
155 
156  for (register label cell=0; cell<nCells; cell++)
157  {
158  pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
159  pTPtr[cell] = wTPtr[cell] + beta*pTPtr[cell];
160  }
161  }
162 
163 
164  // --- Update preconditioned residuals
165  matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);
166  matrix_.Tmul(wT, pT, interfaceIntCoeffs_, interfaces_, cmpt);
167 
168  scalar wApT = gSumProd(wA, pT);
169 
170 
171  // --- Test for singularity
172  if (solverPerf.checkSingularity(mag(wApT)/normFactor)) break;
173 
174 
175  // --- Update solution and residual:
176 
177  scalar alpha = wArT/wApT;
178 
179  for (register label cell=0; cell<nCells; cell++)
180  {
181  psiPtr[cell] += alpha*pAPtr[cell];
182  rAPtr[cell] -= alpha*wAPtr[cell];
183  rTPtr[cell] -= alpha*wTPtr[cell];
184  }
185 
186  solverPerf.finalResidual() = gSumMag(rA)/normFactor;
187 
188  } while
189  (
190  solverPerf.nIterations()++ < maxIter_
191  && !(solverPerf.checkConvergence(tolerance_, relTol_))
192  );
193  }
194 
195  return solverPerf;
196 }
197 
198 
199 // ************************ vim: set sw=4 sts=4 et: ************************ //