src/lduSolvers/lduSolver/bicgSolver/bicgSolver.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 Bi-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 "bicgSolver.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40     defineTypeNameAndDebug(bicgSolver, 0);
  41
  42     lduSolver::addasymMatrixConstructorToTable<bicgSolver>
  43         addbicgSolverAsymMatrixConstructorToTable_;
  44 }
  45
  46
  47 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  48
  49 //- Construct from matrix and solver data stream
  50 Foam::bicgSolver::bicgSolver
  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::bicgSolver::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     scalar normFactor = this->normFactor(x, b, wA, rA, cmpt);
 102
 103     if (lduMatrix::debug >= 2)
 104     {
 105         Info<< "   Normalisation factor = " << normFactor << endl;
 106     }
 107
 108     // Calculate residual
 109     forAll (rA, i)
 110     {
 111         rA[i] = b[i] - wA[i];
 112     }
 113
 114     solverPerf.initialResidual() = gSumMag(rA)/normFactor;
 115     solverPerf.finalResidual() = solverPerf.initialResidual();
 116
 117     if (!stop(solverPerf))
 118     {
 119         scalar rho = matrix_.great_;
 120         scalar rhoOld = rho;
 121
 122         scalar alpha, beta, wApT;
 123
 124         scalarField pA(x.size(), 0);
 125         scalarField pT(x.size(), 0);
 126
 127         // Calculate transpose residual
 128         scalarField wT(x.size());
 129         scalarField rT(rA);
 130
 131         do
 132         {
 133             rhoOld = rho;
 134
 135             // Execute preconditioning
 136             preconPtr_->precondition(wA, rA, cmpt);
 137             preconPtr_->preconditionT(wT, rT, cmpt);
 138
 139             // Update search directions
 140             rho = gSumProd(wA, rT);
 141
 142             beta = rho/rhoOld;
 143
 144             forAll (pA, i)
 145             {
 146                 pA[i] = wA[i] + beta*pA[i];
 147             }
 148
 149             forAll (pT, i)
 150             {
 151                 pT[i] = wT[i] + beta*pT[i];
 152             }
 153
 154             // Update preconditioned residual
 155             matrix_.Amul(wA, pA, coupleBouCoeffs_, interfaces_, cmpt);
 156             matrix_.Amul(wT, pT, coupleIntCoeffs_, interfaces_, cmpt);
 157
 158             wApT = gSumProd(wA, pT);
 159
 160
 161             // Check for singularity
 162             if (solverPerf.checkSingularity(mag(wApT)/normFactor))
 163             {
 164                 break;
 165             }
 166
 167             // Update solution and residual
 168             alpha = rho/wApT;
 169
 170             forAll (x, i)
 171             {
 172                 x[i] += alpha*pA[i];
 173             }
 174
 175             forAll (rA, i)
 176             {
 177                 rA[i] -= alpha*wA[i];
 178             }
 179
 180             forAll (rT, i)
 181             {
 182                 rT[i] -= alpha*wT[i];
 183             }
 184
 185             solverPerf.finalResidual() = gSumMag(rA)/normFactor;
 186             solverPerf.nIterations()++;
 187         } while (!stop(solverPerf));
 188     }
 189
 190     return solverPerf;
 191 }
 192
 193
 194 // ************************************************************************* //