src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.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 Generalised Minimal Residual solver with
  26     run-time selectable preconditioning
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "BlockGMRESSolver.H"
  31 #include "FieldField.H"
  32 #include "scalarMatrices.H"
  33
  34 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  35
  36 template<class Type>
  37 void Foam::BlockGMRESSolver<Type>::givensRotation
  38 (
  39     const scalar& h,
  40     const scalar& beta,
  41     scalar& c,
  42     scalar& s
  43 ) const
  44 {
  45     if (beta == 0)
  46     {
  47         c = 1;
  48         s = 0;
  49     }
  50     else if (mag(beta) > mag(h))
  51     {
  52         scalar tau = -h/beta;
  53         s = 1.0/Foam::sqrt(1.0 + sqr(tau));
  54         c = s*tau;
  55     }
  56     else
  57     {
  58         scalar tau = -beta/h;
  59         c = 1.0/Foam::sqrt(1.0 + sqr(tau));
  60         s = c*tau;
  61     }
  62 }
  63
  64
  65 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  66
  67 //- Construct from matrix and solver data stream
  68 template<class Type>
  69 Foam::BlockGMRESSolver<Type>::BlockGMRESSolver
  70 (
  71     const word& fieldName,
  72     const BlockLduMatrix<Type>& matrix,
  73     const dictionary& dict
  74 )
  75 :
  76     BlockIterativeSolver<Type>
  77     (
  78         fieldName,
  79         matrix,
  80         dict
  81     ),
  82     preconPtr_
  83     (
  84         BlockLduPrecon<Type>::New
  85         (
  86             matrix,
  87             this->dict()
  88         )
  89     ),
  90     nDirs_(readLabel(this->dict().lookup("nDirections")))
  91 {}
  92
  93
  94 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  95
  96 template<class Type>
  97 typename Foam::BlockSolverPerformance<Type>
  98 Foam::BlockGMRESSolver<Type>::solve
  99 (
 100     Field<Type>& x,
 101     const Field<Type>& b
 102 )
 103 {
 104     // Create local references to avoid the spread this-> ugliness
 105     const BlockLduMatrix<Type>& matrix = this->matrix_;
 106
 107     // Prepare solver performance
 108     BlockSolverPerformance<Type> solverPerf
 109     (
 110         typeName,
 111         this->fieldName()
 112     );
 113
 114     scalar norm = this->normFactor(x, b);
 115
 116     Field<Type> wA(x.size());
 117
 118     // Calculate initial residual
 119     matrix.Amul(wA, x);
 120     Field<Type> rA(b - wA);
 121
 122     solverPerf.initialResidual() = gSum(cmptMag(rA))/norm;
 123     solverPerf.finalResidual() = solverPerf.initialResidual();
 124
 125     // Check convergence, solve if not converged
 126
 127     if (!this->stop(solverPerf))
 128     {
 129         // Create the Hesenberg matrix
 130         scalarSquareMatrix H(nDirs_, 0);
 131
 132         // Create y and b for Hessenberg matrix
 133         scalarField yh(nDirs_, 0);
 134         scalarField bh(nDirs_ + 1, 0);
 135
 136         // Givens rotation vectors
 137         scalarField c(nDirs_, 0);
 138         scalarField s(nDirs_, 0);
 139
 140         // Allocate Krylov space vectors
 141         FieldField<Field, Type> V(nDirs_ + 1);
 142
 143         forAll (V, i)
 144         {
 145             V.set(i, new Field<Type>(x.size(), pTraits<Type>::zero));
 146         }
 147
 148         do
 149         {
 150             // Execute preconditioning
 151             preconPtr_->precondition(wA, rA);
 152
 153             // Calculate beta and scale first vector
 154             scalar beta = Foam::sqrt(gSumProd(wA, wA));
 155
 156             // Set initial rhs and bh[0] = beta
 157             bh = 0;
 158             bh[0] = beta;
 159
 160             for (label i = 0; i < nDirs_; i++)
 161             {
 162                 // Set search direction
 163                 V[i] = wA;
 164                 V[i] /= beta;
 165
 166                 // Arnoldi's method
 167                 matrix.Amul(rA, V[i]);
 168
 169                 // Execute preconditioning
 170                 preconPtr_->precondition(wA, rA);
 171
 172                 for (label j = 0; j <= i; j++)
 173                 {
 174                     beta = gSumProd(wA, V[j]);
 175
 176                     H[j][i] = beta;
 177
 178                     forAll (wA, wI)
 179                     {
 180                         wA[wI] -= beta*V[j][wI];
 181                     }
 182                 }
 183
 184                 beta = Foam::sqrt(gSumProd(wA, wA));
 185
 186                 // Apply previous Givens rotations to new column of H.
 187                 for (label j = 0; j < i; j++)
 188                 {
 189                     const scalar Hji = H[j][i];
 190                     H[j][i] = c[j]*Hji - s[j]*H[j + 1][i];
 191                     H[j + 1][i] = s[j]*Hji + c[j]*H[j + 1][i];
 192                 }
 193
 194                 // Apply Givens rotation to current row.
 195                 givensRotation(H[i][i], beta, c[i], s[i]);
 196
 197                 const scalar bhi = bh[i];
 198                 bh[i] = c[i]*bhi - s[i]*bh[i + 1];
 199                 bh[i + 1] = s[i]*bhi + c[i]*bh[i + 1];
 200                 H[i][i] = c[i]*H[i][i] - s[i]*beta;
 201             }
 202
 203             // Back substitute to solve Hy = b
 204             for (label i = nDirs_ - 1; i >= 0; i--)
 205             {
 206                 scalar sum = bh[i];
 207
 208                 for (label j = i + 1; j < nDirs_; j++)
 209                 {
 210                     sum -= H[i][j]*yh[j];
 211                 }
 212
 213                 yh[i] = sum/H[i][i];
 214             }
 215
 216             // Update solution
 217
 218             for (label i = 0; i < nDirs_; i++)
 219             {
 220                 const Field<Type>& Vi = V[i];
 221                 const scalar& yi = yh[i];
 222
 223                 forAll (x, xI)
 224                 {
 225                     x[xI] += yi*Vi[xI];
 226                 }
 227             }
 228
 229             // Re-calculate the residual
 230             matrix.Amul(wA, x);
 231
 232             forAll (rA, raI)
 233             {
 234                 rA[raI] = b[raI] - wA[raI];
 235             }
 236
 237             solverPerf.finalResidual() = gSum(cmptMag(rA))/norm;
 238             solverPerf.nIterations()++;
 239         } while (!this->stop(solverPerf));
 240     }
 241
 242     return solverPerf;
 243 }
 244
 245
 246 // ************************************************************************* //