src/lduSolvers/amg/fineAmgLevel.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     fineAmgLevel
  26
  27 Description
  28     Finest AMG level container, using matrix and field references
  29
  30 Author
  31     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "fineAmgLevel.H"
  36 #include "coarseAmgLevel.H"
  37 #include "cgSolver.H"
  38 #include "bicgStabSolver.H"
  39
  40 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  41
  42 // Construct from components
  43 Foam::fineAmgLevel::fineAmgLevel
  44 (
  45     const lduMatrix& matrix,
  46     const FieldField<Field, scalar>& coupleBouCoeffs,
  47     const FieldField<Field, scalar>& coupleIntCoeffs,
  48     const lduInterfaceFieldPtrsList& interfaceFields,
  49     const dictionary& dict,
  50     const word& policyType,
  51     const label groupSize,
  52     const label minCoarseEqns,
  53     const word& smootherType
  54 )
  55 :
  56     matrix_(matrix),
  57     coupleBouCoeffs_(coupleBouCoeffs),
  58     coupleIntCoeffs_(coupleIntCoeffs),
  59     interfaceFields_(interfaceFields),
  60     dict_(dict),
  61     policyPtr_
  62     (
  63         amgPolicy::New(policyType, matrix_, groupSize, minCoarseEqns)
  64     ),
  65     smootherPtr_
  66     (
  67         lduSmoother::New
  68         (
  69             matrix,
  70             coupleBouCoeffs_,
  71             coupleIntCoeffs_,
  72             interfaceFields_,
  73             dict
  74         )
  75     )
  76 {}
  77
  78
  79 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  80
  81 Foam::scalarField& Foam::fineAmgLevel::x()
  82 {
  83     FatalErrorIn("scalarField& Foam::fineAmgLevel::x()")
  84         << "x is not available."
  85         << abort(FatalError);
  86
  87     // Dummy return
  88     return const_cast<scalarField&>(scalarField::null());
  89 }
  90
  91
  92 Foam::scalarField& Foam::fineAmgLevel::b()
  93 {
  94     FatalErrorIn("scalarField& Foam::fineAmgLevel::b()")
  95         << "b is not available."
  96         << abort(FatalError);
  97
  98     // Dummy return
  99     return const_cast<scalarField&>(scalarField::null());
 100 }
 101
 102
 103 void Foam::fineAmgLevel::residual
 104 (
 105     const scalarField& x,
 106     const scalarField& b,
 107     const direction cmpt,
 108     scalarField& res
 109 ) const
 110 {
 111     matrix_.Amul
 112     (
 113         res,
 114         x,
 115         coupleBouCoeffs_,
 116         interfaceFields_,
 117         cmpt
 118     );
 119
 120     // residual = b - Ax
 121     forAll (b, i)
 122     {
 123         res[i] = b[i] - res[i];
 124     }
 125 }
 126
 127
 128 void Foam::fineAmgLevel::restrictResidual
 129 (
 130     const scalarField& x,
 131     const scalarField& b,
 132     const direction cmpt,
 133     scalarField& xBuffer,
 134     scalarField& coarseRes,
 135     bool preSweepsDone
 136 ) const
 137 {
 138     if (preSweepsDone)
 139     {
 140         // Calculate residual
 141         residual(x, b, cmpt, xBuffer);
 142     }
 143
 144     // Here x != 0.  It is assumed that the buffer will contain the residual
 145     // if no pre-sweeps have been done.  HJ, 4/Sep/2006
 146     policyPtr_->restrictResidual(xBuffer, coarseRes);
 147 }
 148
 149
 150 void Foam::fineAmgLevel::prolongateCorrection
 151 (
 152     scalarField& x,
 153     const scalarField& coarseX
 154 ) const
 155 {
 156     policyPtr_->prolongateCorrection(x, coarseX);
 157 }
 158
 159
 160 void Foam::fineAmgLevel::smooth
 161 (
 162     scalarField& x,
 163     const scalarField& b,
 164     const direction cmpt,
 165     const label nSweeps
 166 ) const
 167 {
 168     smootherPtr_->smooth(x, b, cmpt, nSweeps);
 169 }
 170
 171
 172 void Foam::fineAmgLevel::solve
 173 (
 174     scalarField& x,
 175     const scalarField& b,
 176     const direction cmpt,
 177     const scalar tolerance,
 178     const scalar relTol
 179 ) const
 180 {
 181     Info << "Fine level direct solution" << endl;
 182
 183     lduSolverPerformance coarseSolverPerf;
 184
 185     dictionary fineLevelDict;
 186     fineLevelDict.add("minIter", 1);
 187     fineLevelDict.add("maxIter", 1000);
 188     fineLevelDict.add("tolerance", tolerance);
 189     fineLevelDict.add("relTol", relTol);
 190
 191     // Avoid issues with round-off on strict tolerance setup
 192     // HJ, 27/Jun/2013
 193     x = b/matrix_.diag();
 194
 195     if (matrix_.symmetric())
 196     {
 197         fineLevelDict.add("preconditioner", "Cholesky");
 198
 199         coarseSolverPerf = cgSolver
 200         (
 201             "fineLevelCorr",
 202             matrix_,
 203             coupleBouCoeffs_,
 204             coupleIntCoeffs_,
 205             interfaceFields_,
 206             fineLevelDict
 207         ).solve(x, b, cmpt);
 208     }
 209     else
 210     {
 211         fineLevelDict.add("preconditioner", "ILU0");
 212
 213         coarseSolverPerf = bicgStabSolver
 214         (
 215             "fineLevelCorr",
 216             matrix_,
 217             coupleBouCoeffs_,
 218             coupleIntCoeffs_,
 219             interfaceFields_,
 220             fineLevelDict
 221         ).solve(x, b, cmpt);
 222     }
 223
 224     if (lduMatrix::debug >= 2)
 225     {
 226         coarseSolverPerf.print();
 227     }
 228 }
 229
 230
 231 void Foam::fineAmgLevel::scaleX
 232 (
 233     scalarField& x,
 234     const scalarField& b,
 235     const direction cmpt,
 236     scalarField& xBuffer
 237 ) const
 238 {
 239     // Calculate scaling
 240     scalarField& Ax = xBuffer;
 241
 242     matrix_.Amul
 243     (
 244         Ax,
 245         x,
 246         coupleBouCoeffs_,
 247         interfaceFields_,
 248         cmpt
 249     );
 250
 251     scalar scalingFactorNum = 0.0;
 252     scalar scalingFactorDenom = 0.0;
 253
 254     forAll(x, i)
 255     {
 256         scalingFactorNum += x[i]*b[i];
 257         scalingFactorDenom += x[i]*Ax[i];
 258     }
 259
 260     vector scalingVector(scalingFactorNum, scalingFactorDenom, 0);
 261     reduce(scalingVector, sumOp<vector>());
 262
 263     // Scale x
 264     if
 265     (
 266         mag(scalingVector[0]) > GREAT
 267      || mag(scalingVector[1]) > GREAT
 268      || scalingVector[0]*scalingVector[1] <= 0
 269      || mag(scalingVector[0]) < mag(scalingVector[1])
 270     )
 271     {
 272         // Factor = 1.0, no scaling
 273     }
 274     else if (mag(scalingVector[0]) > 2*mag(scalingVector[1]))
 275     {
 276         // Max factor = 2
 277         x *= 2.0;
 278     }
 279     else
 280     {
 281         // Regular scaling
 282         x *= scalingVector[0]/stabilise(scalingVector[1], SMALL);
 283     }
 284 }
 285
 286
 287 Foam::autoPtr<Foam::amgLevel> Foam::fineAmgLevel::makeNextLevel() const
 288 {
 289     if (policyPtr_->coarsen())
 290     {
 291         return autoPtr<Foam::amgLevel>
 292         (
 293             new coarseAmgLevel
 294             (
 295                 policyPtr_->restrictMatrix
 296                 (
 297                     coupleBouCoeffs_,
 298                     coupleIntCoeffs_,
 299                     interfaceFields_
 300                 ),
 301                 dict(),
 302                 policyPtr_->type(),
 303                 policyPtr_->groupSize(),
 304                 policyPtr_->minCoarseEqns(),
 305                 smootherPtr_->type()
 306             )
 307         );
 308     }
 309     else
 310     {
 311         // Final level: cannot coarsen
 312         return autoPtr<Foam::amgLevel>();
 313     }
 314 }
 315
 316
 317 // ************************************************************************* //