src/coupledMatrix/coupledLduSolver/cgSolver/coupledCgSolver.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Preconditioned Conjugate Gradient solver
  26
  27 Author
  28     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "coupledCgSolver.H"
  33 #include "addToRunTimeSelectionTable.H"
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39     defineTypeNameAndDebug(coupledCgSolver, 0);
  40
  41     coupledLduSolver::addsymMatrixConstructorToTable<coupledCgSolver>
  42         addCgSolverSymMatrixConstructorToTable_;
  43 }
  44
  45
  46 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  47
  48 //- Construct from matrix
  49 Foam::coupledCgSolver::coupledCgSolver
  50 (
  51     const word& fieldName,
  52     const coupledLduMatrix& matrix,
  53     const PtrList<FieldField<Field, scalar> >& bouCoeffs,
  54     const PtrList<FieldField<Field, scalar> >& intCoeffs,
  55     const lduInterfaceFieldPtrsListList& interfaces,
  56     const dictionary& solverData
  57 )
  58 :
  59     coupledIterativeSolver
  60     (
  61         fieldName,
  62         matrix,
  63         bouCoeffs,
  64         intCoeffs,
  65         interfaces,
  66         solverData
  67     ),
  68     preconPtr_
  69     (
  70         coupledLduPrecon::New
  71         (
  72             matrix,
  73             bouCoeffs,
  74             intCoeffs,
  75             interfaces,
  76             dict()
  77         )
  78     )
  79 {}
  80
  81
  82 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  83
  84 Foam::coupledSolverPerformance Foam::coupledCgSolver::solve
  85 (
  86     FieldField<Field, scalar>& x,
  87     const FieldField<Field, scalar>& b,
  88     const direction cmpt
  89 ) const
  90 {
  91     // Prepare solver performance
  92     coupledSolverPerformance solverPerf(typeName, fieldName());
  93
  94     FieldField<Field, scalar> wA(x.size());
  95     FieldField<Field, scalar> rA(x.size());
  96
  97     forAll (x, rowI)
  98     {
  99         wA.set(rowI, new scalarField(x[rowI].size(), 0));
 100         rA.set(rowI, new scalarField(x[rowI].size(), 0));
 101     }
 102
 103
 104     // Calculate initial residual
 105     matrix_.Amul(wA, x, bouCoeffs_, interfaces_, cmpt);
 106
 107     // Use rA as scratch space when calculating the normalisation factor
 108     scalar normFactor = this->normFactor(x, b, wA, rA, cmpt);
 109
 110     if (coupledLduMatrix::debug >= 2)
 111     {
 112         Info<< "   Normalisation factor = " << normFactor << endl;
 113     }
 114
 115     // Optimised looping.  HJ, 19/Jan/2009
 116     forAll (rA, i)
 117     {
 118         const scalarField& bi = b[i];
 119         const scalarField& wAi = wA[i];
 120         scalarField& rAi = rA[i];
 121
 122         forAll (rAi, j)
 123         {
 124             rAi[j] = bi[j] - wAi[j];
 125         }
 126     }
 127
 128     solverPerf.initialResidual() = gSumMag(rA)/normFactor;
 129     solverPerf.finalResidual() = solverPerf.initialResidual();
 130
 131     if (!solverPerf.checkConvergence(tolerance_, relTolerance_))
 132     {
 133         scalar rho = matrix_[0].great_;
 134         scalar rhoOld = rho;
 135
 136         scalar alpha, beta, wApA;
 137
 138         FieldField<Field, scalar> pA(x.size());
 139
 140         forAll (pA, rowI)
 141         {
 142             pA.set(rowI, new scalarField(x[rowI].size(), 0));
 143         }
 144
 145         do
 146         {
 147             rhoOld = rho;
 148
 149             // Execute preconditioning
 150             preconPtr_->precondition(wA, rA, cmpt);
 151
 152             // Update search directions
 153             rho = gSumProd(wA, rA);
 154
 155             beta = rho/rhoOld;
 156
 157             forAll (pA, rowI)
 158             {
 159                 scalarField& curPA = pA[rowI];
 160                 const scalarField& curWA = wA[rowI];
 161
 162                 forAll (curPA, i)
 163                 {
 164                     curPA[i] = curWA[i] + beta*curPA[i];
 165                 }
 166             }
 167
 168             // Update preconditioned residual
 169             matrix_.Amul(wA, pA, bouCoeffs_, interfaces_, cmpt);
 170
 171             wApA = gSumProd(wA, pA);
 172
 173
 174             // Check for singularity
 175             if (solverPerf.checkSingularity(mag(wApA)/normFactor))
 176             {
 177                 break;
 178             }
 179
 180             // Update solution and residual
 181             alpha = rho/wApA;
 182
 183             forAll (x, rowI)
 184             {
 185                 scalarField& curX = x[rowI];
 186                 const scalarField& curPA = pA[rowI];
 187
 188                 forAll (curX, i)
 189                 {
 190                     curX[i] += alpha*curPA[i];
 191                 }
 192             }
 193
 194             forAll (rA, rowI)
 195             {
 196                 scalarField& curRA = rA[rowI];
 197                 const scalarField& curWA = wA[rowI];
 198
 199                 forAll (curRA, i)
 200                 {
 201                     curRA[i] -= alpha*curWA[i];
 202                 }
 203             }
 204
 205             solverPerf.finalResidual() = gSumMag(rA)/normFactor;
 206             solverPerf.nIterations()++;
 207         } while (!stop(solverPerf));
 208     }
 209
 210     return solverPerf;
 211 }
 212
 213
 214 // ************************************************************************* //