src/coupledMatrix/coupledLduSolver/cgSolver/coupledCgSolver.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
  27
  28 Author
  29     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "coupledCgSolver.H"
  34 #include "addToRunTimeSelectionTable.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40     defineTypeNameAndDebug(coupledCgSolver, 0);
  41
  42     coupledLduSolver::addsymMatrixConstructorToTable<coupledCgSolver>
  43         addCgSolverSymMatrixConstructorToTable_;
  44 }
  45
  46
  47 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  48
  49 //- Construct from matrix
  50 Foam::coupledCgSolver::coupledCgSolver
  51 (
  52     const word& fieldName,
  53     const coupledLduMatrix& matrix,
  54     const PtrList<FieldField<Field, scalar> >& bouCoeffs,
  55     const PtrList<FieldField<Field, scalar> >& intCoeffs,
  56     const lduInterfaceFieldPtrsListList& interfaces,
  57     const dictionary& solverData
  58 )
  59 :
  60     coupledIterativeSolver
  61     (
  62         fieldName,
  63         matrix,
  64         bouCoeffs,
  65         intCoeffs,
  66         interfaces,
  67         solverData
  68     ),
  69     preconPtr_
  70     (
  71         coupledLduPrecon::New
  72         (
  73             matrix,
  74             bouCoeffs,
  75             intCoeffs,
  76             interfaces,
  77             dict().subDict("preconditioner")
  78         )
  79     )
  80 {}
  81
  82
  83 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  84
  85 Foam::coupledSolverPerformance Foam::coupledCgSolver::solve
  86 (
  87     FieldField<Field, scalar>& x,
  88     const FieldField<Field, scalar>& b,
  89     const direction cmpt
  90 ) const
  91 {
  92     // Prepare solver performance
  93     coupledSolverPerformance solverPerf(typeName, fieldName());
  94
  95     FieldField<Field, scalar> wA(x.size());
  96     FieldField<Field, scalar> rA(x.size());
  97
  98     forAll (x, rowI)
  99     {
 100         wA.set(rowI, new scalarField(x[rowI].size(), 0));
 101         rA.set(rowI, new scalarField(x[rowI].size(), 0));
 102     }
 103
 104
 105     // Calculate initial residual
 106     matrix_.Amul(wA, x, bouCoeffs_, interfaces_, cmpt);
 107
 108     // Use rA as scratch space when calculating the normalisation factor
 109     scalar normFactor = this->normFactor(x, b, wA, rA, cmpt);
 110
 111     if (coupledLduMatrix::debug >= 2)
 112     {
 113         Info<< "   Normalisation factor = " << normFactor << endl;
 114     }
 115
 116     // Optimised looping.  HJ, 19/Jan/2009
 117     forAll (rA, i)
 118     {
 119         const scalarField& bi = b[i];
 120         const scalarField& wAi = wA[i];
 121         scalarField& rAi = rA[i];
 122
 123         forAll (rAi, j)
 124         {
 125             rAi[j] = bi[j] - wAi[j];
 126         }
 127     }
 128
 129     solverPerf.initialResidual() = gSumMag(rA)/normFactor;
 130     solverPerf.finalResidual() = solverPerf.initialResidual();
 131
 132     if (!solverPerf.checkConvergence(tolerance_, relTolerance_))
 133     {
 134         scalar rho = matrix_[0].great_;
 135         scalar rhoOld = rho;
 136
 137         scalar alpha, beta, wApA;
 138
 139         FieldField<Field, scalar> pA(x.size());
 140
 141         forAll (pA, rowI)
 142         {
 143             pA.set(rowI, new scalarField(x[rowI].size(), 0));
 144         }
 145
 146         do
 147         {
 148             rhoOld = rho;
 149
 150             // Execute preconditioning
 151             preconPtr_->precondition(wA, rA, cmpt);
 152
 153             // Update search directions
 154             rho = gSumProd(wA, rA);
 155
 156             beta = rho/rhoOld;
 157
 158             forAll (pA, rowI)
 159             {
 160                 scalarField& curPA = pA[rowI];
 161                 const scalarField& curWA = wA[rowI];
 162
 163                 forAll (curPA, i)
 164                 {
 165                     curPA[i] = curWA[i] + beta*curPA[i];
 166                 }
 167             }
 168
 169             // Update preconditioned residual
 170             matrix_.Amul(wA, pA, bouCoeffs_, interfaces_, cmpt);
 171
 172             wApA = gSumProd(wA, pA);
 173
 174
 175             // Check for singularity
 176             if (solverPerf.checkSingularity(mag(wApA)/normFactor))
 177             {
 178                 break;
 179             }
 180
 181             // Update solution and residual
 182             alpha = rho/wApA;
 183
 184             forAll (x, rowI)
 185             {
 186                 scalarField& curX = x[rowI];
 187                 const scalarField& curPA = pA[rowI];
 188
 189                 forAll (curX, i)
 190                 {
 191                     curX[i] += alpha*curPA[i];
 192                 }
 193             }
 194
 195             forAll (rA, rowI)
 196             {
 197                 scalarField& curRA = rA[rowI];
 198                 const scalarField& curWA = wA[rowI];
 199
 200                 forAll (curRA, i)
 201                 {
 202                     curRA[i] -= alpha*curWA[i];
 203                 }
 204             }
 205
 206             solverPerf.finalResidual() = gSumMag(rA)/normFactor;
 207             solverPerf.nIterations()++;
 208         } while (!stop(solverPerf));
 209     }
 210
 211     return solverPerf;
 212 }
 213
 214
 215 // ************************************************************************* //