src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.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 Class
  25     coarseBlockAmgLevel
  26
  27 Description
  28     Coarse AMG level stores matrix, x and b locally, for BlockLduMatrix
  29
  30 Author
  31     Klas Jareteg, 2012-12-13
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "coarseBlockAmgLevel.H"
  36 #include "SubField.H"
  37 #include "vector2D.H"
  38 #include "coeffFields.H"
  39 #include "BlockSolverPerformance.H"
  40 // #include "BlockBiCGStabSolver.H"
  41 // #include "BlockCGSolver.H"
  42 #include "BlockGMRESSolver.H"
  43
  44 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  45
  46 // Construct from components
  47 template<class Type>
  48 Foam::coarseBlockAmgLevel<Type>::coarseBlockAmgLevel
  49 (
  50     autoPtr<lduPrimitiveMesh> addrPtr,
  51     autoPtr<BlockLduMatrix<Type> > matrixPtr,
  52     const dictionary& dict,
  53     const word& coarseningType,
  54     const label groupSize,
  55     const label minCoarseEqns
  56 )
  57 :
  58     addrPtr_(addrPtr),
  59     matrixPtr_(matrixPtr),
  60     x_(matrixPtr_->diag().size(),pTraits<Type>::zero),
  61     b_(matrixPtr_->diag().size(),pTraits<Type>::zero),
  62     dict_(dict),
  63     coarseningPtr_
  64     (
  65         BlockMatrixCoarsening<Type>::New
  66         (
  67             coarseningType,
  68             matrixPtr_,
  69             dict_,
  70             groupSize,
  71             minCoarseEqns
  72         )
  73     ),
  74     smootherPtr_
  75     (
  76         BlockLduSmoother<Type>::New
  77         (
  78             matrixPtr_,
  79             dict,
  80             "coarseSmoother"
  81         )
  82     ),
  83     Ax_()
  84 {}
  85
  86
  87 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  88
  89 template<class Type>
  90 Foam::coarseBlockAmgLevel<Type>::~coarseBlockAmgLevel()
  91 {}
  92
  93
  94 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  95
  96 template<class Type>
  97 Foam::Field<Type>& Foam::coarseBlockAmgLevel<Type>::x()
  98 {
  99     return x_;
 100 }
 101
 102
 103 template<class Type>
 104 Foam::Field<Type>& Foam::coarseBlockAmgLevel<Type>::b()
 105 {
 106     return b_;
 107 }
 108
 109
 110 template<class Type>
 111 void Foam::coarseBlockAmgLevel<Type>::residual
 112 (
 113     const Field<Type>& x,
 114     const Field<Type>& b,
 115     Field<Type>& res
 116 ) const
 117 {
 118     // Calculate residual
 119     matrixPtr_->Amul
 120     (
 121         res,
 122         x
 123     );
 124
 125     // residual = b - Ax
 126     forAll (b, i)
 127     {
 128         res[i] = b[i] - res[i];
 129     }
 130 }
 131
 132
 133 template<class Type>
 134 void Foam::coarseBlockAmgLevel<Type>::restrictResidual
 135 (
 136     const Field<Type>& x,
 137     const Field<Type>& b,
 138     Field<Type>& xBuffer,
 139     Field<Type>& coarseRes,
 140     bool preSweepsDone
 141 ) const
 142 {
 143     if (preSweepsDone)
 144     {
 145         // Calculate residual
 146         // KRJ: 2012-12-14 ::subField removed. Creating buffer locally
 147         Field<Type> resBuf(x.size());
 148
 149         Field<Type>& res = reinterpret_cast<Field<Type>&>(resBuf);
 150
 151         residual(x, b, res);
 152
 153         coarseningPtr_->restrictResidual(res, coarseRes);
 154     }
 155     else
 156     {
 157         // No pre-sweeps done: x = 0 and residual = b
 158         coarseningPtr_->restrictResidual(b, coarseRes);
 159     }
 160 }
 161
 162
 163 template<class Type>
 164 void Foam::coarseBlockAmgLevel<Type>::prolongateCorrection
 165 (
 166     Field<Type>& x,
 167     const Field<Type>& coarseX
 168 ) const
 169 {
 170     coarseningPtr_->prolongateCorrection(x, coarseX);
 171 }
 172
 173
 174 template<class Type>
 175 void Foam::coarseBlockAmgLevel<Type>::smooth
 176 (
 177     Field<Type>& x,
 178     const Field<Type>& b,
 179     const label nSweeps
 180 ) const
 181 {
 182     smootherPtr_->smooth(x, b, nSweeps);
 183 }
 184
 185
 186 template<class Type>
 187 void Foam::coarseBlockAmgLevel<Type>::solve
 188 (
 189     Field<Type>& x,
 190     const Field<Type>& b,
 191     const scalar tolerance,
 192     const scalar relTol
 193 ) const
 194 {
 195     BlockSolverPerformance<Type> coarseSolverPerf
 196     (
 197         BlockGMRESSolver<Type>::typeName,
 198         "topLevelCorr"
 199     );
 200
 201     label maxIter = Foam::min(2*coarseningPtr_->minCoarseEqns(), 1000);
 202
 203     // Create artificial dictionary for top-level solution
 204     dictionary topLevelDict;
 205     topLevelDict.add("nDirections", "5");
 206     topLevelDict.add("minIter", 1);
 207     topLevelDict.add("maxIter", maxIter);
 208     topLevelDict.add("tolerance", tolerance);
 209     topLevelDict.add("relTol", relTol);
 210
 211     // Avoid issues with round-off on strict tolerance setup
 212     // HJ, 27/Jun/2013
 213     // Create multiplication function object
 214     typename BlockCoeff<Type>::multiply mult;
 215
 216 //     x = inverse(diag) & b
 217     CoeffField<Type> invDiag = inv(matrixPtr_->diag());
 218     multiply(x, invDiag, b);
 219
 220     // Do not solve if the number of equations is smaller than 5
 221     if (coarseningPtr_->minCoarseEqns() < 5)
 222     {
 223         return;
 224     }
 225
 226     // Switch of debug in top-level direct solve
 227     label oldDebug = BlockLduMatrix<Type>::debug();
 228
 229     if (BlockLduMatrix<Type>::debug >= 4)
 230     {
 231         BlockLduMatrix<Type>::debug = 1;
 232     }
 233     else
 234     {
 235         BlockLduMatrix<Type>::debug = 0;
 236     }
 237
 238     if (matrixPtr_->symmetric())
 239     {
 240         topLevelDict.add("preconditioner", "Cholesky");
 241
 242         coarseSolverPerf =
 243 //         BlockCGSolver<Type>
 244         BlockGMRESSolver<Type>
 245         (
 246             "topLevelCorr",
 247             matrixPtr_,
 248             topLevelDict
 249         ).solve(x, b);
 250     }
 251     else
 252     {
 253         topLevelDict.add("preconditioner", "Cholesky");
 254
 255         coarseSolverPerf =
 256 //         BlockBiCGStabSolver<Type>
 257         BlockGMRESSolver<Type>
 258         (
 259             "topLevelCorr",
 260             matrixPtr_,
 261             topLevelDict
 262         ).solve(x, b);
 263     }
 264
 265     // Restore debug
 266     BlockLduMatrix<Type>::debug = oldDebug;
 267
 268     // Escape cases of top-level solver divergence
 269     if
 270     (
 271         coarseSolverPerf.nIterations() == maxIter
 272      && (
 273             coarseSolverPerf.finalResidual()
 274          >= coarseSolverPerf.initialResidual()
 275         )
 276     )
 277     {
 278         // Top-level solution failed.  Attempt rescue
 279         // HJ, 27/Jul/2013
 280         multiply(x, invDiag, b);
 281
 282         // Print top level correction failure as info for user
 283         coarseSolverPerf.print();
 284     }
 285
 286     if (BlockLduMatrix<Type>::debug >= 2)
 287     {
 288         coarseSolverPerf.print();
 289     }
 290 }
 291
 292
 293 template<class Type>
 294 void Foam::coarseBlockAmgLevel<Type>::scaleX
 295 (
 296     Field<Type>& x,
 297     const Field<Type>& b,
 298     Field<Type>& xBuffer
 299 ) const
 300 {
 301     // KRJ: 2013-02-05: Removed subfield, creating a new field
 302     // Initialise and size buffer to avoid multiple allocation.
 303     // Buffer is created as private data of AMG level
 304     // HJ, 26/Feb/2015
 305     if (Ax_.empty())
 306     {
 307         Ax_.setSize(x.size());
 308     }
 309
 310     matrixPtr_->Amul(Ax_, x);
 311
 312 #if 0
 313
 314     // Variant 1 (old): scale complete x with a single scaling factor
 315     scalar scalingFactorNum = sumProd(x, b);
 316     scalar scalingFactorDenom = sumProd(x, Ax_);
 317
 318     vector2D scalingVector(scalingFactorNum, scalingFactorDenom);
 319     reduce(scalingVector, sumOp<vector2D>());
 320
 321     // Scale x
 322     if
 323     (
 324         mag(scalingVector[0]) > GREAT
 325      || mag(scalingVector[1]) > GREAT
 326      || scalingVector[0]*scalingVector[1] <= 0
 327      || mag(scalingVector[0]) < mag(scalingVector[1])
 328     )
 329     {
 330         // Factor = 1.0, no scaling
 331     }
 332     else if (mag(scalingVector[0]) > 2*mag(scalingVector[1]))
 333     {
 334         // Max factor = 2
 335         x *= 2.0;
 336     }
 337     else
 338     {
 339         // Regular scaling
 340         x *= scalingVector[0]/stabilise(scalingVector[1], VSMALL);
 341     }
 342
 343 #else
 344
 345     // Variant 2: scale each x individually
 346     // HJ, 25/Feb/2015
 347     Type scalingFactorNum = sumCmptProd(x, b);
 348     Type scalingFactorDenom = sumCmptProd(x, Ax_);
 349
 350     Vector2D<Type> scalingVector(scalingFactorNum, scalingFactorDenom);
 351     reduce(scalingVector, sumOp<Vector2D<Type> >());
 352
 353     Type scalingFactor = pTraits<Type>::one;
 354
 355     // Scale x
 356     for (direction dir = 0; dir < pTraits<Type>::nComponents; dir++)
 357     {
 358         scalar num = component(scalingVector[0], dir);
 359         scalar denom = component(scalingVector[1], dir);
 360
 361         if
 362         (
 363             mag(num) > GREAT || mag(denom) > GREAT
 364          || num*denom <= 0 || mag(num) < mag(denom)
 365         )
 366         {
 367             // Factor = 1.0, no scaling
 368         }
 369         else if (mag(num) > 2*mag(denom))
 370         {
 371             setComponent(scalingFactor, dir) = 2;
 372         }
 373         else
 374         {
 375             // Regular scaling
 376             setComponent(scalingFactor, dir) = num/stabilise(denom, VSMALL);
 377         }
 378     }
 379
 380     // Scale X
 381     cmptMultiply(x, x, scalingFactor);
 382
 383 #endif
 384 }
 385
 386
 387 template<class Type>
 388 Foam::autoPtr<Foam::BlockAmgLevel<Type> >
 389 Foam::coarseBlockAmgLevel<Type>::makeNextLevel() const
 390 {
 391     if (coarseningPtr_->coarsen())
 392     {
 393         return coarseningPtr_->restrictMatrix();
 394     }
 395     else
 396     {
 397         // Final level: cannot coarsen
 398         return autoPtr<BlockAmgLevel<Type> >();
 399     }
 400 }
 401
 402
 403 // ************************************************************************* //