src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "PCG.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33     defineTypeNameAndDebug(PCG, 0);
  34
  35     lduMatrix::solver::addsymMatrixConstructorToTable<PCG>
  36         addPCGSymMatrixConstructorToTable_;
  37 }
  38
  39
  40 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  41
  42 Foam::PCG::PCG
  43 (
  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
  50 )
  51 :
  52     lduMatrix::solver
  53     (
  54         fieldName,
  55         matrix,
  56         interfaceBouCoeffs,
  57         interfaceIntCoeffs,
  58         interfaces,
  59         solverControls
  60     )
  61 {}
  62
  63
  64 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  65
  66 Foam::lduMatrix::solverPerformance Foam::PCG::solve
  67 (
  68     scalarField& psi,
  69     const scalarField& source,
  70     const direction cmpt
  71 ) const
  72 {
  73     // --- Setup class containing solver performance data
  74     lduMatrix::solverPerformance solverPerf
  75     (
  76         lduMatrix::preconditioner::getName(controlDict_) + typeName,
  77         fieldName_
  78     );
  79
  80     register label nCells = psi.size();
  81
  82     scalar* __restrict__ psiPtr = psi.begin();
  83
  84     scalarField pA(nCells);
  85     scalar* __restrict__ pAPtr = pA.begin();
  86
  87     scalarField wA(nCells);
  88     scalar* __restrict__ wAPtr = wA.begin();
  89
  90     scalar wArA = matrix_.great_;
  91     scalar wArAold = wArA;
  92
  93     // --- Calculate A.psi
  94     matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
  95
  96     // --- Calculate initial residual field
  97     scalarField rA(source - wA);
  98     scalar* __restrict__ rAPtr = rA.begin();
  99
 100     // --- Calculate normalisation factor
 101     scalar normFactor = this->normFactor(psi, source, wA, pA);
 102
 103     if (lduMatrix::debug >= 2)
 104     {
 105         Info<< "   Normalisation factor = " << normFactor << endl;
 106     }
 107
 108     // --- Calculate normalised residual norm
 109     solverPerf.initialResidual() = gSumMag(rA)/normFactor;
 110     solverPerf.finalResidual() = solverPerf.initialResidual();
 111
 112     // --- Check convergence, solve if not converged
 113     if (!solverPerf.checkConvergence(tolerance_, relTol_))
 114     {
 115         // --- Select and construct the preconditioner
 116         autoPtr<lduMatrix::preconditioner> preconPtr =
 117         lduMatrix::preconditioner::New
 118         (
 119             *this,
 120             controlDict_
 121         );
 122
 123         // --- Solver iteration
 124         do
 125         {
 126             // --- Store previous wArA
 127             wArAold = wArA;
 128
 129             // --- Precondition residual
 130             preconPtr->precondition(wA, rA, cmpt);
 131
 132             // --- Update search directions:
 133             wArA = gSumProd(wA, rA);
 134
 135             if (solverPerf.nIterations() == 0)
 136             {
 137                 for (register label cell=0; cell<nCells; cell++)
 138                 {
 139                     pAPtr[cell] = wAPtr[cell];
 140                 }
 141             }
 142             else
 143             {
 144                 scalar beta = wArA/wArAold;
 145
 146                 for (register label cell=0; cell<nCells; cell++)
 147                 {
 148                     pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
 149                 }
 150             }
 151
 152
 153             // --- Update preconditioned residual
 154             matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);
 155
 156             scalar wApA = gSumProd(wA, pA);
 157
 158
 159             // --- Test for singularity
 160             if (solverPerf.checkSingularity(mag(wApA)/normFactor)) break;
 161
 162
 163             // --- Update solution and residual:
 164
 165             scalar alpha = wArA/wApA;
 166
 167             for (register label cell=0; cell<nCells; cell++)
 168             {
 169                 psiPtr[cell] += alpha*pAPtr[cell];
 170                 rAPtr[cell] -= alpha*wAPtr[cell];
 171             }
 172
 173             solverPerf.finalResidual() = gSumMag(rA)/normFactor;
 174
 175         } while
 176         (
 177             solverPerf.nIterations()++ < maxIter_
 178         && !(solverPerf.checkConvergence(tolerance_, relTol_))
 179         );
 180     }
 181
 182     return solverPerf;
 183 }
 184
 185
 186 // ************************************************************************* //