src/lduSolvers/lduSolver/bicgStabSolver/bicgStabSolver.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 Bi-Conjugate Gradient stabilised solver with run-time
  26     selectable preconditioning
  27
  28 Author
  29     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "bicgStabSolver.H"
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39     defineTypeNameAndDebug(bicgStabSolver, 0);
  40
  41     //HJ, temporary
  42     lduSolver::addsymMatrixConstructorToTable<bicgStabSolver>
  43         addbicgStabSolverSymMatrixConstructorToTable_;
  44
  45     lduSolver::addasymMatrixConstructorToTable<bicgStabSolver>
  46         addbicgStabSolverAsymMatrixConstructorToTable_;
  47 }
  48
  49
  50 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  51
  52 //- Construct from matrix and solver data stream
  53 Foam::bicgStabSolver::bicgStabSolver
  54 (
  55     const word& fieldName,
  56     const lduMatrix& matrix,
  57     const FieldField<Field, scalar>& coupleBouCoeffs,
  58     const FieldField<Field, scalar>& coupleIntCoeffs,
  59     const lduInterfaceFieldPtrsList& interfaces,
  60     const dictionary& dict
  61 )
  62 :
  63     lduSolver
  64     (
  65         fieldName,
  66         matrix,
  67         coupleBouCoeffs,
  68         coupleIntCoeffs,
  69         interfaces,
  70         dict
  71     ),
  72     preconPtr_
  73     (
  74         lduPreconditioner::New
  75         (
  76             matrix,
  77             coupleBouCoeffs,
  78             coupleIntCoeffs,
  79             interfaces,
  80             dict
  81         )
  82     )
  83 {}
  84
  85
  86 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  87
  88 Foam::lduSolverPerformance Foam::bicgStabSolver::solve
  89 (
  90     scalarField& x,
  91     const scalarField& b,
  92     const direction cmpt
  93 ) const
  94 {
  95     // Prepare solver performance
  96     lduSolverPerformance solverPerf(typeName, fieldName());
  97
  98     scalarField p(x.size());
  99     scalarField r(x.size());
 100
 101     // Calculate initial residual
 102     matrix_.Amul(p, x, coupleBouCoeffs_, interfaces_, cmpt);
 103
 104     scalar normFactor = this->normFactor(x, b, p, r, cmpt);
 105
 106     if (lduMatrix::debug >= 2)
 107     {
 108         Info<< "   Normalisation factor = " << normFactor << endl;
 109     }
 110
 111     // Calculate residual
 112     forAll (r, i)
 113     {
 114         r[i] = b[i] - p[i];
 115     }
 116
 117     solverPerf.initialResidual() = gSumMag(r)/normFactor;
 118     solverPerf.finalResidual() = solverPerf.initialResidual();
 119
 120     if (!stop(solverPerf))
 121     {
 122         scalar rho = matrix_.great_;
 123         scalar rhoOld = rho;
 124
 125         scalar alpha = 0;
 126         scalar omega = matrix_.great_;
 127         scalar beta;
 128
 129         p = 0;
 130         scalarField ph(x.size(), 0);
 131         scalarField v(x.size(), 0);
 132         scalarField s(x.size(), 0);
 133         scalarField sh(x.size(), 0);
 134         scalarField t(x.size(), 0);
 135
 136         // Calculate transpose residual
 137         scalarField rw(r);
 138
 139         do
 140         {
 141             rhoOld = rho;
 142
 143             // Update search directions
 144             rho = gSumProd(rw, r);
 145
 146             beta = rho/rhoOld*(alpha/omega);
 147
 148             // Restart if breakdown occurs
 149             if (rho == 0)
 150             {
 151                 rw = r;
 152                 rho = gSumProd(rw, r);
 153
 154                 alpha = 0;
 155                 omega = 0;
 156                 beta = 0;
 157             }
 158
 159             forAll (p, i)
 160             {
 161                 p[i] = r[i] + beta*p[i] - beta*omega*v[i];
 162             }
 163
 164             // Execute preconditioning
 165             preconPtr_->precondition(ph, p, cmpt);
 166             matrix_.Amul(v, ph, coupleBouCoeffs_, interfaces_, cmpt);
 167             alpha = rho/gSumProd(rw, v);
 168
 169             forAll (s, i)
 170             {
 171                 s[i] = r[i] - alpha*v[i];
 172             }
 173
 174             // Execute preconditioning
 175             // Bug fix, Alexander Monakov, 11/Jul/2012
 176             preconPtr_->precondition(sh, s, cmpt);
 177             matrix_.Amul(t, sh, coupleBouCoeffs_, interfaces_, cmpt);
 178             omega = gSumProd(t, s)/gSumProd(t, t);
 179
 180             // Update solution and residual
 181             forAll (x, i)
 182             {
 183                 x[i] = x[i] + alpha*ph[i] + omega*sh[i];
 184             }
 185
 186             forAll (r, i)
 187             {
 188                 r[i] = s[i] - omega*t[i];
 189             }
 190
 191             solverPerf.finalResidual() = gSumMag(r)/normFactor;
 192             solverPerf.nIterations()++;
 193         } while (!stop(solverPerf));
 194     }
 195
 196     return solverPerf;
 197 }
 198
 199
 200 // ************************************************************************* //