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