src/coupledMatrix/coupledLduSolver/bicgStabSolver/coupledBicgStabSolver.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 Stabilised Bi-Conjugate Gradient solver
  26
  27 Author
  28     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "coupledBicgStabSolver.H"
  33 #include "addToRunTimeSelectionTable.H"
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39     defineTypeNameAndDebug(coupledBicgStabSolver, 0);
  40
  41     coupledLduSolver::addsymMatrixConstructorToTable<coupledBicgStabSolver>
  42         addBicgStabSolverSymMatrixConstructorToTable_;
  43     coupledLduSolver::addasymMatrixConstructorToTable<coupledBicgStabSolver>
  44         addBicgStabSolverAsymMatrixConstructorToTable_;
  45 }
  46
  47
  48 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  49
  50 //- Construct from matrix
  51 Foam::coupledBicgStabSolver::coupledBicgStabSolver
  52 (
  53     const word& fieldName,
  54     const coupledLduMatrix& matrix,
  55     const PtrList<FieldField<Field, scalar> >& bouCoeffs,
  56     const PtrList<FieldField<Field, scalar> >& intCoeffs,
  57     const lduInterfaceFieldPtrsListList& interfaces,
  58     const dictionary& solverData
  59 )
  60 :
  61     coupledIterativeSolver
  62     (
  63         fieldName,
  64         matrix,
  65         bouCoeffs,
  66         intCoeffs,
  67         interfaces,
  68         solverData
  69     ),
  70     preconPtr_
  71     (
  72         coupledLduPrecon::New
  73         (
  74             matrix,
  75             bouCoeffs,
  76             intCoeffs,
  77             interfaces,
  78             dict()
  79         )
  80     )
  81 {}
  82
  83
  84 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  85
  86 Foam::coupledSolverPerformance Foam::coupledBicgStabSolver::solve
  87 (
  88     FieldField<Field, scalar>& x,
  89     const FieldField<Field, scalar>& b,
  90     const direction cmpt
  91 ) const
  92 {
  93     // Prepare solver performance
  94     coupledSolverPerformance solverPerf(typeName, fieldName());
  95
  96     FieldField<Field, scalar> p(x.size());
  97     FieldField<Field, scalar> r(x.size());
  98
  99     forAll (x, rowI)
 100     {
 101         p.set(rowI, new scalarField(x[rowI].size(), 0));
 102         r.set(rowI, new scalarField(x[rowI].size(), 0));
 103     }
 104
 105
 106     // Calculate initial residual
 107     matrix_.Amul(p, x, bouCoeffs_, interfaces_, cmpt);
 108
 109     // Use r as scratch space when calculating the normalisation factor
 110     scalar normFactor = this->normFactor(x, b, p, r, cmpt);
 111
 112     if (coupledLduMatrix::debug >= 2)
 113     {
 114         Info<< "   Normalisation factor = " << normFactor << endl;
 115     }
 116
 117     // Calculate residual
 118     forAll (r, rowI)
 119     {
 120         const scalarField& curB = b[rowI];
 121         const scalarField& curP = p[rowI];
 122         scalarField& curR = r[rowI];
 123
 124         forAll (curR, i)
 125         {
 126             curR[i] = curB[i] - curP[i];
 127         }
 128     }
 129
 130     solverPerf.initialResidual() = gSumMag(r)/normFactor;
 131     solverPerf.finalResidual() = solverPerf.initialResidual();
 132
 133     if (!solverPerf.checkConvergence(tolerance_, relTolerance_))
 134     {
 135         scalar rho = matrix_[0].great_;
 136         scalar rhoOld = rho;
 137
 138         scalar alpha = 0;
 139         scalar omega = matrix_[0].great_;
 140         scalar beta;
 141
 142         p = 0;
 143         FieldField<Field, scalar> ph(x.size());
 144         FieldField<Field, scalar> v(x.size());
 145         FieldField<Field, scalar> s(x.size());
 146         FieldField<Field, scalar> sh(x.size());
 147         FieldField<Field, scalar> t(x.size());
 148
 149         forAll (ph, rowI)
 150         {
 151             ph.set(rowI, new scalarField(x[rowI].size(), 0));
 152             v.set(rowI, new scalarField(x[rowI].size(), 0));
 153             s.set(rowI, new scalarField(x[rowI].size(), 0));
 154             sh.set(rowI, new scalarField(x[rowI].size(), 0));
 155             t.set(rowI, new scalarField(x[rowI].size(), 0));
 156         }
 157
 158         // Calculate transpose residual
 159         FieldField<Field, scalar> rw(r);
 160
 161         do
 162         {
 163             rhoOld = rho;
 164
 165             // Update search directions
 166             rho = gSumProd(rw, r);
 167
 168             beta = rho/rhoOld*(alpha/omega);
 169
 170             // Restart if breakdown occurs
 171             if (rho == 0)
 172             {
 173                 rw = r;
 174                 rho = gSumProd(rw, r);
 175
 176                 alpha = 0;
 177                 omega = 0;
 178                 beta = 0;
 179             }
 180
 181             forAll (p, rowI)
 182             {
 183                 scalarField& curP = p[rowI];
 184                 const scalarField& curR = r[rowI];
 185                 const scalarField& curV = v[rowI];
 186
 187                 forAll (curP, i)
 188                 {
 189                     curP[i] = curR[i] + beta*curP[i] - beta*omega*curV[i];
 190                 }
 191             }
 192
 193             // Execute preconditioning
 194             preconPtr_->precondition(ph, p, cmpt);
 195             matrix_.Amul(v, ph, bouCoeffs_, interfaces_, cmpt);
 196             alpha = rho/gSumProd(rw, v);
 197
 198             forAll (s, rowI)
 199             {
 200                 scalarField& curS = s[rowI];
 201                 const scalarField& curR = r[rowI];
 202                 const scalarField& curV = v[rowI];
 203
 204                 forAll (curS, i)
 205                 {
 206                     curS[i] = curR[i] - alpha*curV[i];
 207                 }
 208             }
 209
 210             // Execute preconditioning
 211             // Bug fix, Alexander Monakov, 11/Jul/2012
 212             preconPtr_->precondition(sh, s, cmpt);
 213             matrix_.Amul(t, sh, bouCoeffs_, interfaces_, cmpt);
 214             omega = gSumProd(t, s)/gSumProd(t, t);
 215
 216             // Update solution and residual
 217             forAll (x, rowI)
 218             {
 219                 scalarField& curX = x[rowI];
 220                 const scalarField& curPh = ph[rowI];
 221                 const scalarField& curSh = sh[rowI];
 222
 223                 forAll (curX, i)
 224                 {
 225                     curX[i] = curX[i] + alpha*curPh[i] + omega*curSh[i];
 226                 }
 227             }
 228
 229             forAll (r, rowI)
 230             {
 231                 scalarField& curR = r[rowI];
 232                 const scalarField& curS = s[rowI];
 233                 const scalarField& curT = t[rowI];
 234
 235                 forAll (curR, i)
 236                 {
 237                     curR[i] = curS[i] - omega*curT[i];
 238                 }
 239             }
 240
 241             solverPerf.finalResidual() = gSumMag(r)/normFactor;
 242             solverPerf.nIterations()++;
 243         } while (!stop(solverPerf));
 244     }
 245
 246     return solverPerf;
 247 }
 248
 249
 250 // ************************************************************************* //