src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.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     Gauss-Seidel solver for symmetric and asymmetric matrices.  In
  26     order to improve efficiency, the residual is evaluated after every
  27     nSweeps sweeps.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "BlockGaussSeidelSolver.H"
  32
  33 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  34
  35 // Construct from matrix and solver data stream
  36 template<class Type>
  37 Foam::BlockGaussSeidelSolver<Type>::BlockGaussSeidelSolver
  38 (
  39     const word& fieldName,
  40     const BlockLduMatrix<Type>& matrix,
  41     const dictionary& dict
  42 )
  43 :
  44     BlockIterativeSolver<Type>
  45     (
  46         fieldName,
  47         matrix,
  48         dict
  49     ),
  50     gs_(matrix),
  51     nSweeps_(readLabel(this->dict().lookup("nSweeps")))
  52 {}
  53
  54
  55 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  56
  57 template<class Type>
  58 typename Foam::BlockSolverPerformance<Type>
  59 Foam::BlockGaussSeidelSolver<Type>::solve
  60 (
  61     Field<Type>& x,
  62     const Field<Type>& b
  63 )
  64 {
  65     // Create local references to avoid the spread this-> ugliness
  66     const BlockLduMatrix<Type>& matrix = this->matrix_;
  67
  68     // Prepare solver performance
  69     BlockSolverPerformance<Type> solverPerf
  70     (
  71         typeName,
  72         this->fieldName()
  73     );
  74
  75     scalar norm = this->normFactor(x, b);
  76
  77     Field<Type> wA(x.size());
  78
  79     // Calculate residual.  Note: sign of residual swapped for efficiency
  80     matrix.Amul(wA, x);
  81     wA -= b;
  82
  83     solverPerf.initialResidual() = gSum(cmptMag(wA))/norm;
  84     solverPerf.finalResidual() = solverPerf.initialResidual();
  85
  86     // Check convergence, solve if not converged
  87
  88     if (!this->stop(solverPerf))
  89     {
  90         // Iteration loop
  91
  92         do
  93         {
  94             for (label i = 0; i < nSweeps_; i++)
  95             {
  96                 gs_.precondition(x, b);
  97
  98                 solverPerf.nIterations()++;
  99             }
 100
 101             // Re-calculate residual.  Note: sign of residual swapped
 102             // for efficiency
 103             matrix.Amul(wA, x);
 104             wA -= b;
 105
 106             solverPerf.finalResidual() = gSum(cmptMag(wA))/norm;
 107             solverPerf.nIterations()++;
 108         } while (!this->stop(solverPerf));
 109     }
 110
 111     return solverPerf;
 112 }
 113
 114
 115 // ************************************************************************* //