src/lduSolvers/amg/coarseAmgLevel.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     coarseAmgLevel
  26
  27 Description
  28     Coarse AMG level stores matrix, x and b locally
  29
  30 Author
  31     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "coarseAmgLevel.H"
  36 #include "SubField.H"
  37 #include "gmresSolver.H"
  38 #include "vector2D.H"
  39
  40 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  41
  42 // Construct from components
  43 Foam::coarseAmgLevel::coarseAmgLevel
  44 (
  45     autoPtr<amgMatrix> matrixPtr,
  46     const dictionary& dict,
  47     const word& policyType,
  48     const label groupSize,
  49     const label minCoarseEqns,
  50     const word& smootherType
  51 )
  52 :
  53     matrixPtr_(matrixPtr),
  54     x_(matrixPtr_->size()),
  55     b_(matrixPtr_->size()),
  56     dict_(dict),
  57     policyPtr_
  58     (
  59         amgPolicy::New
  60         (
  61             policyType,
  62             matrixPtr_->matrix(),
  63             groupSize,
  64             minCoarseEqns
  65         )
  66     ),
  67     smootherPtr_
  68     (
  69         lduSmoother::New
  70         (
  71             matrixPtr_->matrix(),
  72             matrixPtr_->coupleBouCoeffs(),
  73             matrixPtr_->coupleIntCoeffs(),
  74             matrixPtr_->interfaceFields(),
  75             dict
  76         )
  77     )
  78 {}
  79
  80
  81 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  82
  83 Foam::coarseAmgLevel::~coarseAmgLevel()
  84 {}
  85
  86
  87 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  88
  89 Foam::scalarField& Foam::coarseAmgLevel::x()
  90 {
  91     return x_;
  92 }
  93
  94
  95 Foam::scalarField& Foam::coarseAmgLevel::b()
  96 {
  97     return b_;
  98 }
  99
 100
 101 void Foam::coarseAmgLevel::residual
 102 (
 103     const scalarField& x,
 104     const scalarField& b,
 105     const direction cmpt,
 106     scalarField& res
 107 ) const
 108 {
 109     // Calculate residual
 110     matrixPtr_->matrix().Amul
 111     (
 112         res,
 113         x,
 114         matrixPtr_->coupleBouCoeffs(),
 115         matrixPtr_->interfaceFields(),
 116         cmpt
 117     );
 118
 119     // residual = b - Ax
 120     forAll (b, i)
 121     {
 122         res[i] = b[i] - res[i];
 123     }
 124 }
 125
 126
 127 void Foam::coarseAmgLevel::restrictResidual
 128 (
 129     const scalarField& x,
 130     const scalarField& b,
 131     const direction cmpt,
 132     scalarField& xBuffer,
 133     scalarField& coarseRes,
 134     bool preSweepsDone
 135 ) const
 136 {
 137     if (preSweepsDone)
 138     {
 139         // Calculate residual
 140         scalarField::subField resBuf(xBuffer, x.size());
 141
 142         scalarField& res = const_cast<scalarField&>
 143         (
 144             resBuf.operator const scalarField&()
 145         );
 146
 147         residual(x, b, cmpt, res);
 148
 149         policyPtr_->restrictResidual(res, coarseRes);
 150     }
 151     else
 152     {
 153         // No pre-sweeps done: x = 0 and residual = b
 154         policyPtr_->restrictResidual(b, coarseRes);
 155     }
 156 }
 157
 158
 159 void Foam::coarseAmgLevel::prolongateCorrection
 160 (
 161     scalarField& x,
 162     const scalarField& coarseX
 163 ) const
 164 {
 165     policyPtr_->prolongateCorrection(x, coarseX);
 166 }
 167
 168
 169 void Foam::coarseAmgLevel::smooth
 170 (
 171     scalarField& x,
 172     const scalarField& b,
 173     const direction cmpt,
 174     const label nSweeps
 175 ) const
 176 {
 177     smootherPtr_->smooth(x, b, cmpt, nSweeps);
 178 }
 179
 180
 181 void Foam::coarseAmgLevel::solve
 182 (
 183     scalarField& x,
 184     const scalarField& b,
 185     const direction cmpt,
 186     const scalar tolerance,
 187     const scalar relTol
 188 ) const
 189 {
 190     lduSolverPerformance coarseSolverPerf;
 191
 192     label maxIter = Foam::min(2*policyPtr_->minCoarseEqns(), 1000);
 193
 194     dictionary topLevelDict;
 195     topLevelDict.add("nDirections", "5");
 196     topLevelDict.add("minIter", 1);
 197     topLevelDict.add("maxIter", maxIter);
 198     topLevelDict.add("tolerance", tolerance);
 199     topLevelDict.add("relTol", relTol);
 200
 201     // Avoid issues with round-off on strict tolerance setup
 202     // HJ, 27/Jun/2013
 203     x = b/matrixPtr_->matrix().diag();
 204
 205     // Do not solve if the number of equations is smaller than 5
 206     if (policyPtr_->minCoarseEqns() < 5)
 207     {
 208         return;
 209     }
 210
 211     // Switch of debug in top-level direct solve
 212     label oldDebug = lduMatrix::debug();
 213
 214     if (matrixPtr_->matrix().symmetric())
 215     {
 216         topLevelDict.add("preconditioner", "Cholesky");
 217
 218         coarseSolverPerf = gmresSolver
 219         (
 220             "topLevelCorr",
 221             matrixPtr_->matrix(),
 222             matrixPtr_->coupleBouCoeffs(),
 223             matrixPtr_->coupleIntCoeffs(),
 224             matrixPtr_->interfaceFields(),
 225             topLevelDict
 226         ).solve(x, b, cmpt);
 227     }
 228     else
 229     {
 230         topLevelDict.add("preconditioner", "ILU0");
 231
 232         coarseSolverPerf = gmresSolver
 233         (
 234             "topLevelCorr",
 235             matrixPtr_->matrix(),
 236             matrixPtr_->coupleBouCoeffs(),
 237             matrixPtr_->coupleIntCoeffs(),
 238             matrixPtr_->interfaceFields(),
 239             topLevelDict
 240         ).solve(x, b, cmpt);
 241     }
 242
 243     // Restore debug
 244     lduMatrix::debug = oldDebug;
 245
 246     // Escape cases of top-level solver divergence
 247     if
 248     (
 249         coarseSolverPerf.nIterations() == maxIter
 250      && (
 251             coarseSolverPerf.finalResidual()
 252          >= coarseSolverPerf.initialResidual()
 253         )
 254     )
 255     {
 256         // Top-level solution failed.  Attempt rescue
 257         // HJ, 27/Jul/2013
 258         x = b/matrixPtr_->matrix().diag();
 259
 260         // Print top level correction failure as info for user
 261         coarseSolverPerf.print();
 262     }
 263
 264     if (lduMatrix::debug >= 2)
 265     {
 266         coarseSolverPerf.print();
 267     }
 268 }
 269
 270
 271 void Foam::coarseAmgLevel::scaleX
 272 (
 273     scalarField& x,
 274     const scalarField& b,
 275     const direction cmpt,
 276     scalarField& xBuffer
 277 ) const
 278 {
 279     // Calculate scaling
 280     scalarField::subField Ax(xBuffer, x.size());
 281
 282     matrixPtr_->matrix().Amul
 283     (
 284         const_cast<scalarField&>(Ax.operator const scalarField&()),
 285         x,
 286         matrixPtr_->coupleBouCoeffs(),
 287         matrixPtr_->interfaceFields(),
 288         cmpt
 289     );
 290
 291     scalar scalingFactorNum = 0.0;
 292     scalar scalingFactorDenom = 0.0;
 293
 294     forAll(x, i)
 295     {
 296         scalingFactorNum += x[i]*b[i];
 297         scalingFactorDenom += x[i]*xBuffer[i]; // Note: Ax = xBuffer
 298     }
 299
 300     vector2D scalingVector(scalingFactorNum, scalingFactorDenom);
 301     reduce(scalingVector, sumOp<vector2D>());
 302
 303     // Scale x
 304     if
 305     (
 306         mag(scalingVector[0]) > GREAT
 307      || mag(scalingVector[1]) > GREAT
 308      || scalingVector[0]*scalingVector[1] <= 0
 309      || mag(scalingVector[0]) < mag(scalingVector[1])
 310     )
 311     {
 312         // Factor = 1.0, no scaling
 313     }
 314     else if (mag(scalingVector[0]) > 2*mag(scalingVector[1]))
 315     {
 316         // Max factor = 2
 317         x *= 2.0;
 318     }
 319     else
 320     {
 321         // Regular scaling
 322         x *= scalingVector[0]/stabilise(scalingVector[1], SMALL);
 323     }
 324 }
 325
 326
 327 Foam::autoPtr<Foam::amgLevel> Foam::coarseAmgLevel::makeNextLevel() const
 328 {
 329     if (policyPtr_->coarsen())
 330     {
 331         return autoPtr<Foam::amgLevel>
 332         (
 333             new coarseAmgLevel
 334             (
 335                 policyPtr_->restrictMatrix
 336                 (
 337                     matrixPtr_->coupleBouCoeffs(),
 338                     matrixPtr_->coupleIntCoeffs(),
 339                     matrixPtr_->interfaceFields()
 340                 ),
 341                 dict(),
 342                 policyPtr_->type(),
 343                 policyPtr_->groupSize(),
 344                 policyPtr_->minCoarseEqns(),
 345                 smootherPtr_->type()
 346             )
 347         );
 348     }
 349     else
 350     {
 351         // Final level: cannot coarsen
 352         return autoPtr<Foam::amgLevel>();
 353     }
 354 }
 355
 356
 357 // ************************************************************************* //