1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include <OpenFOAM/PCG.H>
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 defineTypeNameAndDebug(PCG, 0);
35 lduMatrix::solver::addsymMatrixConstructorToTable<PCG>
36 addPCGSymMatrixConstructorToTable_;
40 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
44 const word& fieldName,
45 const lduMatrix& matrix,
46 const FieldField<Field, scalar>& interfaceBouCoeffs,
47 const FieldField<Field, scalar>& interfaceIntCoeffs,
48 const lduInterfaceFieldPtrsList& interfaces,
49 const dictionary& solverControls
64 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
66 Foam::lduMatrix::solverPerformance Foam::PCG::solve
69 const scalarField& source,
73 // --- Setup class containing solver performance data
74 lduMatrix::solverPerformance solverPerf
76 lduMatrix::preconditioner::getName(controlDict_) + typeName,
80 register label nCells = psi.size();
82 scalar* __restrict__ psiPtr = psi.begin();
84 scalarField pA(nCells);
85 scalar* __restrict__ pAPtr = pA.begin();
87 scalarField wA(nCells);
88 scalar* __restrict__ wAPtr = wA.begin();
90 scalar wArA = matrix_.great_;
91 scalar wArAold = wArA;
93 // --- Calculate A.psi
94 matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
96 // --- Calculate initial residual field
97 scalarField rA(source - wA);
98 scalar* __restrict__ rAPtr = rA.begin();
100 // --- Calculate normalisation factor
101 scalar normFactor = this->normFactor(psi, source, wA, pA);
103 if (lduMatrix::debug >= 2)
105 Info<< " Normalisation factor = " << normFactor << endl;
108 // --- Calculate normalised residual norm
109 solverPerf.initialResidual() = gSumMag(rA)/normFactor;
110 solverPerf.finalResidual() = solverPerf.initialResidual();
112 // --- Check convergence, solve if not converged
113 if (!solverPerf.checkConvergence(tolerance_, relTol_))
115 // --- Select and construct the preconditioner
116 autoPtr<lduMatrix::preconditioner> preconPtr =
117 lduMatrix::preconditioner::New
123 // --- Solver iteration
126 // --- Store previous wArA
129 // --- Precondition residual
130 preconPtr->precondition(wA, rA, cmpt);
132 // --- Update search directions:
133 wArA = gSumProd(wA, rA);
135 if (solverPerf.nIterations() == 0)
137 for (register label cell=0; cell<nCells; cell++)
139 pAPtr[cell] = wAPtr[cell];
144 scalar beta = wArA/wArAold;
146 for (register label cell=0; cell<nCells; cell++)
148 pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
153 // --- Update preconditioned residual
154 matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);
156 scalar wApA = gSumProd(wA, pA);
159 // --- Test for singularity
160 if (solverPerf.checkSingularity(mag(wApA)/normFactor)) break;
163 // --- Update solution and residual:
165 scalar alpha = wArA/wApA;
167 for (register label cell=0; cell<nCells; cell++)
169 psiPtr[cell] += alpha*pAPtr[cell];
170 rAPtr[cell] -= alpha*wAPtr[cell];
173 solverPerf.finalResidual() = gSumMag(rA)/normFactor;
177 solverPerf.nIterations()++ < maxIter_
178 && !(solverPerf.checkConvergence(tolerance_, relTol_))
186 // ************************ vim: set sw=4 sts=4 et: ************************ //