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