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