src/lduSolvers/lduSolver/cgSolver/cgSolver.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 Description
  26     Preconditioned Conjugate Gradient solver with run-time selectable
  27     preconditioning
  28
  29 Author
  30     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  31
  32 \*---------------------------------------------------------------------------*/
  33
  34 #include "cgSolver.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40     defineTypeNameAndDebug(cgSolver, 0);
  41
  42     lduSolver::addsymMatrixConstructorToTable<cgSolver>
  43         addcgSolverSymMatrixConstructorToTable_;
  44 }
  45
  46
  47 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  48
  49 //- Construct from matrix and solver data stream
  50 Foam::cgSolver::cgSolver
  51 (
  52     const word& fieldName,
  53     const lduMatrix& matrix,
  54     const FieldField<Field, scalar>& coupleBouCoeffs,
  55     const FieldField<Field, scalar>& coupleIntCoeffs,
  56     const lduInterfaceFieldPtrsList& interfaces,
  57     const dictionary& dict
  58 )
  59 :
  60     lduSolver
  61     (
  62         fieldName,
  63         matrix,
  64         coupleBouCoeffs,
  65         coupleIntCoeffs,
  66         interfaces,
  67         dict
  68     ),
  69     preconPtr_
  70     (
  71         lduPreconditioner::New
  72         (
  73             matrix,
  74             coupleBouCoeffs,
  75             coupleIntCoeffs,
  76             interfaces,
  77             dict
  78         )
  79     )
  80 {}
  81
  82
  83 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  84
  85 Foam::lduSolverPerformance Foam::cgSolver::solve
  86 (
  87     scalarField& x,
  88     const scalarField& b,
  89     const direction cmpt
  90 ) const
  91 {
  92     // Prepare solver performance
  93     lduSolverPerformance solverPerf(typeName, fieldName());
  94
  95     scalarField wA(x.size());
  96     scalarField rA(x.size());
  97
  98     // Calculate initial residual
  99     matrix_.Amul(wA, x, coupleBouCoeffs_, interfaces_, cmpt);
 100
 101     // Use rA as scratch space when calculating the normalisation factor
 102     scalar normFactor = this->normFactor(x, b, wA, rA, cmpt);
 103
 104     if (lduMatrix::debug >= 2)
 105     {
 106         Info<< "   Normalisation factor = " << normFactor << endl;
 107     }
 108
 109     // Calculate residual
 110     forAll (rA, i)
 111     {
 112         rA[i] = b[i] - wA[i];
 113     }
 114
 115     solverPerf.initialResidual() = gSumMag(rA)/normFactor;
 116     solverPerf.finalResidual() = solverPerf.initialResidual();
 117
 118     if (!stop(solverPerf))
 119     {
 120         scalar rho = matrix_.great_;
 121         scalar rhoOld = rho;
 122
 123         scalar alpha, beta, wApA;
 124
 125         scalarField pA(x.size(), 0);
 126
 127         do
 128         {
 129             rhoOld = rho;
 130
 131             // Execute preconditioning
 132             preconPtr_->precondition(wA, rA, cmpt);
 133
 134             // Update search directions
 135             rho = gSumProd(wA, rA);
 136
 137             beta = rho/rhoOld;
 138
 139             forAll (pA, i)
 140             {
 141                 pA[i] = wA[i] + beta*pA[i];
 142             }
 143
 144
 145             // Update preconditioned residual
 146             matrix_.Amul(wA, pA, coupleBouCoeffs_, interfaces_, cmpt);
 147
 148             wApA = gSumProd(wA, pA);
 149
 150
 151             // Check for singularity
 152             if (solverPerf.checkSingularity(mag(wApA)/normFactor))
 153             {
 154                 break;
 155             }
 156
 157             // Update solution and residual
 158             alpha = rho/wApA;
 159
 160             forAll (x, i)
 161             {
 162                 x[i] += alpha*pA[i];
 163             }
 164
 165             forAll (rA, i)
 166             {
 167                 rA[i] -= alpha*wA[i];
 168             }
 169
 170             solverPerf.finalResidual() = gSumMag(rA)/normFactor;
 171             solverPerf.nIterations()++;
 172         } while (!stop(solverPerf));
 173     }
 174
 175     return solverPerf;
 176 }
 177
 178
 179 // ************************************************************************* //