src/foam/matrices/blockLduMatrix/BlockAmg/fineBlockAmgLevel.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   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 Class
  25     fineBlockAmgLevel
  26
  27 Description
  28     Finest AMG level container, using matrix and field references
  29
  30 Author
  31     Klas Jareteg, 2013-04-15
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "fineBlockAmgLevel.H"
  36 #include "coarseBlockAmgLevel.H"
  37 #include "ICCG.H"
  38 #include "BICCG.H"
  39 #include "BlockGaussSeidelSolver.H"
  40 #include "BlockSolverPerformance.H"
  41 #include "BlockCoeffNorm.H"
  42 #include "BlockCoeffTwoNorm.H"
  43
  44 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  45
  46 // Construct from components
  47 template<class Type>
  48 Foam::fineBlockAmgLevel<Type>::fineBlockAmgLevel
  49 (
  50     const BlockLduMatrix<Type>& matrix,
  51     const dictionary& dict,
  52     const word& coarseningType,
  53     const label groupSize,
  54     const label minCoarseEqns,
  55     const word& smootherType
  56 )
  57 :
  58     matrix_(matrix),
  59     dict_(dict),
  60     coarseningPtr_
  61     (
  62         BlockMatrixCoarsening<Type>::New
  63         (
  64             coarseningType,
  65             matrix_,
  66             dict_,
  67             groupSize,
  68             minCoarseEqns
  69         )
  70     ),
  71     smootherPtr_
  72     (
  73         BlockLduSmoother<Type>::New
  74         (
  75             matrix,
  76             dict
  77         )
  78     )
  79 {}
  80
  81
  82 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  83
  84 template<class Type>
  85 Foam::Field<Type>& Foam::fineBlockAmgLevel<Type>::x()
  86 {
  87     FatalErrorIn("Field<Type>& Foam::fineBlockAmgLevel<Type>::x()")
  88         << "x is not available."
  89         << abort(FatalError);
  90
  91     // Dummy return
  92     return const_cast<Field<Type>&>(Field<Type>::null());
  93 }
  94
  95
  96 template<class Type>
  97 Foam::Field<Type>& Foam::fineBlockAmgLevel<Type>::b()
  98 {
  99     FatalErrorIn("Field<Type>& Foam::fineBlockAmgLevel<Type>::b()")
 100         << "b is not available."
 101         << abort(FatalError);
 102
 103     // Dummy return
 104     return const_cast<Field<Type>&>(Field<Type>::null());
 105 }
 106
 107
 108 template<class Type>
 109 void Foam::fineBlockAmgLevel<Type>::residual
 110 (
 111     const Field<Type>& x,
 112     const Field<Type>& b,
 113     Field<Type>& res
 114 ) const
 115 {
 116     matrix_.Amul
 117     (
 118         res,
 119         x
 120     );
 121
 122     // residual = b - Ax
 123     forAll (b, i)
 124     {
 125         res[i] = b[i] - res[i];
 126     }
 127 }
 128
 129
 130 template<class Type>
 131 void Foam::fineBlockAmgLevel<Type>::restrictResidual
 132 (
 133     const Field<Type>& x,
 134     const Field<Type>& b,
 135     Field<Type>& xBuffer,
 136     Field<Type>& coarseRes,
 137     bool preSweepsDone
 138 ) const
 139 {
 140     if (preSweepsDone)
 141     {
 142         // Calculate residual
 143         residual(x, b, xBuffer);
 144     }
 145
 146     // Here x != 0.  It is assumed that the buffer will contain the residual
 147     // if no pre-sweeps have been done.  HJ, 4/Sep/2006
 148     coarseningPtr_->restrictResidual(xBuffer, coarseRes);
 149 }
 150
 151
 152 template<class Type>
 153 void Foam::fineBlockAmgLevel<Type>::prolongateCorrection
 154 (
 155     Field<Type>& x,
 156     const Field<Type>& coarseX
 157 ) const
 158 {
 159     coarseningPtr_->prolongateCorrection(x, coarseX);
 160 }
 161
 162
 163 template<class Type>
 164 void Foam::fineBlockAmgLevel<Type>::smooth
 165 (
 166     Field<Type>& x,
 167     const Field<Type>& b,
 168     const label nSweeps
 169 ) const
 170 {
 171     smootherPtr_->smooth(x, b, nSweeps);
 172 }
 173
 174
 175 template<class Type>
 176 void Foam::fineBlockAmgLevel<Type>::solve
 177 (
 178     Field<Type>& x,
 179     const Field<Type>& b,
 180     const scalar tolerance,
 181     const scalar relTol
 182 ) const
 183 {
 184
 185     Info<<"Fine level solver"<<endl;
 186
 187     if (matrix_.symmetric())
 188     {
 189         BlockSolverPerformance<Type> coarseSolverPerf =
 190             BlockCGSolver<Type>
 191             (
 192                 "topLevelCorr",
 193                 matrix_,
 194                 dict_
 195             ).solve(x, b);
 196
 197         if (lduMatrix::debug >= 2)
 198         {
 199             coarseSolverPerf.print();
 200         }
 201     }
 202     else
 203     {
 204         BlockSolverPerformance<Type> coarseSolverPerf =
 205             BlockBiCGStabSolver<Type>
 206             (
 207                 "topLevelCorr",
 208                 matrix_,
 209                 dict_
 210             ).solve(x, b);
 211
 212         if (lduMatrix::debug >= 2)
 213         {
 214             coarseSolverPerf.print();
 215         }
 216     }
 217 }
 218
 219
 220 template<class Type>
 221 void Foam::fineBlockAmgLevel<Type>::scaleX
 222 (
 223     Field<Type>& x,
 224     const Field<Type>& b,
 225     Field<Type>& xBuffer
 226 ) const
 227 {
 228
 229     // KRJ: 2013-02-05: Removed subfield, creating a new field
 230     Field<Type> Ax(x.size());
 231
 232     matrix_.Amul
 233     (
 234         reinterpret_cast<Field<Type>&>(Ax),
 235         x
 236     );
 237
 238     scalar scalingFactorNum = sumProd(x,b);
 239     scalar scalingFactorDenom = sumProd(x,Ax);
 240
 241     vector scalingVector(scalingFactorNum, scalingFactorDenom, 0);
 242     reduce(scalingVector, sumOp<vector>());
 243
 244     // Scale x
 245     if
 246     (
 247         scalingVector[0]*scalingVector[1] <= 0
 248      || mag(scalingVector[0]) < mag(scalingVector[1])
 249     )
 250     {
 251         // Factor = 1.0, no scaling
 252     }
 253     else if (mag(scalingVector[0]) > 2*mag(scalingVector[1]))
 254     {
 255         // Max factor = 2
 256         x *= 2.0;
 257     }
 258     else
 259     {
 260         // Regular scaling
 261         x *= scalingVector[0]/stabilise(scalingVector[1], SMALL);
 262     }
 263 }
 264
 265
 266 template<class Type>
 267 Foam::autoPtr<Foam::BlockAmgLevel<Type> >
 268 Foam::fineBlockAmgLevel<Type>::makeNextLevel() const
 269 {
 270     if (coarseningPtr_->coarsen())
 271     {
 272         return autoPtr<Foam::BlockAmgLevel<Type> >
 273         (
 274             new coarseBlockAmgLevel<Type>
 275             (
 276                 coarseningPtr_->restrictMatrix(),
 277                 dict(),
 278                 coarseningPtr_->type(),
 279                 coarseningPtr_->groupSize(),
 280                 coarseningPtr_->minCoarseEqns(),
 281                 smootherPtr_->type()
 282             )
 283         );
 284     }
 285     else
 286     {
 287         // Final level: cannot coarsen
 288         return autoPtr<Foam::BlockAmgLevel<Type> >();
 289     }
 290 }
 291
 292
 293 // ************************************************************************* //