src/foam/matrices/blockLduMatrix/BlockAmg/BlockAmgCycle.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     BlockAmgCycle
  26
  27 Description
  28     Algebraic multigrid cycle class for BlockLduMatrix
  29
  30 Author
  31     Klas Jareteg, 2012-12-12
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "BlockAmgCycle.H"
  36 #include "demandDrivenData.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40
  41
  42 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  43
  44 // Construct from AMG level
  45 template<class Type>
  46 Foam::BlockAmgCycle<Type>::BlockAmgCycle
  47 (
  48     autoPtr<BlockAmgLevel<Type> > levelPtr
  49 )
  50 :
  51     levelPtr_(levelPtr),
  52     coarseLevelPtr_(NULL),
  53     nLevels_(0)
  54 {}
  55
  56 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  57
  58 template<class Type>
  59 Foam::BlockAmgCycle<Type>::~BlockAmgCycle()
  60 {
  61     deleteDemandDrivenData(coarseLevelPtr_);
  62 }
  63
  64
  65 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  66
  67 template<class Type>
  68 void Foam::BlockAmgCycle<Type>::makeCoarseLevels(const label nMaxLevels)
  69 {
  70     // Make coarse levels
  71     if (nLevels_ == 0)
  72     {
  73         bool addCoarse = true;
  74         BlockAmgCycle<Type>* curCyclePtr = this;
  75
  76         // Do forever
  77         for (;;)
  78         {
  79             nLevels_++;
  80
  81             autoPtr<BlockAmgLevel<Type> > coarsePtr =
  82                 curCyclePtr->levelPtr_->makeNextLevel();
  83
  84             // Check if a coarse level is valid and allowed
  85             if (nLevels_ >= nMaxLevels || !coarsePtr.valid())
  86             {
  87                 addCoarse = false;
  88             }
  89
  90             reduce(addCoarse, andOp<bool>());
  91
  92             if (addCoarse)
  93             {
  94                 curCyclePtr->coarseLevelPtr_ =
  95                     new BlockAmgCycle<Type>(coarsePtr);
  96
  97                 // Point to the next coarse level
  98                 curCyclePtr = curCyclePtr->coarseLevelPtr_;
  99             }
 100             else
 101             {
 102                 break;
 103             }
 104         }
 105
 106         if (BlockLduMatrix<Type>::debug >= 2)
 107         {
 108             Info<< "Created " << nLevels_ << " AMG levels" << endl;
 109         }
 110     }
 111 }
 112
 113
 114 template<class Type>
 115 void Foam::BlockAmgCycle<Type>::fixedCycle
 116 (
 117     Field<Type>& x,
 118     const Field<Type>& b,
 119     Field<Type>& xBuffer,
 120     const blockAmgCycleName::cycleType cycle,
 121     const label nPreSweeps,
 122     const label nPostSweeps,
 123     const bool scale
 124 ) const
 125 {
 126     if (coarseLevelPtr_)
 127     {
 128         // Pre-smoothing
 129         levelPtr_->smooth(x, b, nPreSweeps);
 130
 131         // Get reference to coarse level
 132         Field<Type>& xCoarse = coarseLevelPtr_->levelPtr_->x();
 133         Field<Type>& bCoarse = coarseLevelPtr_->levelPtr_->b();
 134
 135         // Zero out coarse x
 136         xCoarse = pTraits<Type>::zero;
 137
 138         // Restrict residual: optimisation on number of pre-sweeps
 139         levelPtr_->restrictResidual
 140         (
 141             x,
 142             b,
 143             xBuffer,
 144             bCoarse,
 145             nPreSweeps > 0 || cycle != V_CYCLE
 146         );
 147
 148         coarseLevelPtr_->fixedCycle
 149         (
 150             xCoarse,
 151             bCoarse,
 152             xBuffer,
 153             cycle,
 154             nPreSweeps,
 155             nPostSweeps,
 156             scale
 157         );
 158
 159         if (cycle == F_CYCLE)
 160         {
 161             coarseLevelPtr_->fixedCycle
 162             (
 163                 xCoarse,
 164                 bCoarse,
 165                 xBuffer,
 166                 V_CYCLE,
 167                 nPreSweeps,
 168                 nPostSweeps,
 169                 scale
 170             );
 171         }
 172         else if (cycle == W_CYCLE)
 173         {
 174             coarseLevelPtr_->fixedCycle
 175             (
 176                 xCoarse,
 177                 bCoarse,
 178                 xBuffer,
 179                 W_CYCLE,
 180                 nPreSweeps,
 181                 nPostSweeps,
 182                 scale
 183             );
 184         }
 185
 186         if (scale)
 187         {
 188             // Calculate scaling factor using a buffer
 189             coarseLevelPtr_->levelPtr_->scaleX(xCoarse, bCoarse, xBuffer);
 190         }
 191
 192         levelPtr_->prolongateCorrection(x, xCoarse);
 193
 194         // Post-smoothing
 195         levelPtr_->smooth(x, b, nPostSweeps);
 196     }
 197     else
 198     {
 199         // Call direct solver
 200         levelPtr_->solve(x, b, 1e-9, 0);
 201     }
 202 }
 203
 204
 205 // ************************************************************************* //