src/foam/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "GAMGSolver.H"
  27 #include "ICCG.H"
  28 #include "BICCG.H"
  29 #include "SubField.H"
  30
  31 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  32
  33 Foam::lduMatrix::solverPerformance Foam::GAMGSolver::solve
  34 (
  35     scalarField& x,
  36     const scalarField& b,
  37     const direction cmpt
  38 ) const
  39 {
  40     // Setup class containing solver performance data
  41     lduSolverPerformance solverPerf(typeName, fieldName());
  42
  43     // Calculate A.x used to calculate the initial residual
  44     scalarField Ax(x.size());
  45     matrix_.Amul(Ax, x, coupleBouCoeffs_, interfaces_, cmpt);
  46
  47     // Create the storage for the finestCorrection which may be used as a
  48     // temporary in normFactor
  49     scalarField finestCorrection(x.size());
  50
  51     // Calculate normalisation factor
  52     scalar normFactor = this->normFactor(x, b, Ax, finestCorrection, cmpt);
  53
  54     if (debug >= 2)
  55     {
  56         Info<< "   Normalisation factor = " << normFactor << endl;
  57     }
  58
  59     // Calculate initial finest-grid residual field
  60     scalarField finestResidual(b - Ax);
  61
  62     // Calculate normalised residual for convergence test
  63     solverPerf.initialResidual() = gSumMag(finestResidual)/normFactor;
  64     solverPerf.finalResidual() = solverPerf.initialResidual();
  65
  66
  67     // Check convergence, solve if not converged
  68     if (!stop(solverPerf))
  69     {
  70         // Create coarse grid correction fields
  71         PtrList<scalarField> coarseCorrX;
  72
  73         // Create coarse grid bs
  74         PtrList<scalarField> coarseB;
  75
  76         // Create the smoothers for all levels
  77         PtrList<lduSmoother> smoothers;
  78
  79         // Initialise the above data structures
  80         initVcycle(coarseCorrX, coarseB, smoothers);
  81
  82         do
  83         {
  84             Vcycle
  85             (
  86                 smoothers,
  87                 x,
  88                 b,
  89                 Ax,
  90                 finestCorrection,
  91                 finestResidual,
  92                 coarseCorrX,
  93                 coarseB,
  94                 cmpt
  95             );
  96
  97             // Calculate finest level residual field
  98             matrix_.Amul(Ax, x, coupleBouCoeffs_, interfaces_, cmpt);
  99             finestResidual = b;
 100             finestResidual -= Ax;
 101
 102             solverPerf.finalResidual() = gSumMag(finestResidual)/normFactor;
 103             solverPerf.nIterations()++;
 104             if (debug >= 2)
 105             {
 106                 solverPerf.print();
 107             }
 108         } while (!stop(solverPerf));
 109     }
 110
 111     return solverPerf;
 112 }
 113
 114
 115 void Foam::GAMGSolver::Vcycle
 116 (
 117     const PtrList<lduSmoother>& smoothers,
 118     scalarField& x,
 119     const scalarField& b,
 120     scalarField& Ax,
 121     scalarField& finestCorrection,
 122     scalarField& finestResidual,
 123     PtrList<scalarField>& coarseCorrX,
 124     PtrList<scalarField>& coarseB,
 125     const direction cmpt
 126 ) const
 127 {
 128     //debug = 2;
 129
 130     const label coarsestLevel = matrixLevels_.size() - 1;
 131
 132     // Restrict finest grid residual for the next level up
 133     agglomeration_.restrictField(coarseB[0], finestResidual, 0);
 134
 135     if (debug >= 2 && nPreSweeps_)
 136     {
 137         Pout<< "Pre-smoothing scaling factors: ";
 138     }
 139
 140
 141     // Residual restriction (going to coarser levels)
 142     for (label leveli = 0; leveli < coarsestLevel; leveli++)
 143     {
 144         // If the optional pre-smoothing sweeps are selected
 145         // smooth the coarse-grid field for the restricted b
 146         if (nPreSweeps_)
 147         {
 148             coarseCorrX[leveli] = 0.0;
 149
 150             smoothers[leveli + 1].smooth
 151             (
 152                 coarseCorrX[leveli],
 153                 coarseB[leveli],
 154                 cmpt,
 155                 nPreSweeps_ + leveli
 156             );
 157
 158             scalarField::subField ACf
 159             (
 160                 Ax,
 161                 coarseCorrX[leveli].size()
 162             );
 163
 164             // Scale coarse-grid correction field
 165             // but not on the coarsest level because it evaluates to 1
 166             if (scaleCorrection_ && leveli < coarsestLevel - 1)
 167             {
 168                 scalar sf = scalingFactor
 169                 (
 170                     const_cast<scalarField&>(ACf.operator const scalarField&()),
 171                     matrixLevels_[leveli],
 172                     coarseCorrX[leveli],
 173                     coupleLevelsBouCoeffs_[leveli],
 174                     interfaceLevels_[leveli],
 175                     coarseB[leveli],
 176                     cmpt
 177                 );
 178
 179                 if (debug >= 2)
 180                 {
 181                     Pout<< sf << " ";
 182                 }
 183
 184                 coarseCorrX[leveli] *= sf;
 185             }
 186
 187             // Correct the residual with the new solution
 188             matrixLevels_[leveli].Amul
 189             (
 190                 const_cast<scalarField&>(ACf.operator const scalarField&()),
 191                 coarseCorrX[leveli],
 192                 coupleLevelsBouCoeffs_[leveli],
 193                 interfaceLevels_[leveli],
 194                 cmpt
 195             );
 196
 197             coarseB[leveli] -= ACf;
 198         }
 199
 200         // Residual is equal to b
 201         agglomeration_.restrictField
 202         (
 203             coarseB[leveli + 1],
 204             coarseB[leveli],
 205             leveli + 1
 206         );
 207     }
 208
 209     if (debug >= 2 && nPreSweeps_)
 210     {
 211         Pout<< endl;
 212     }
 213
 214
 215     // Solve Coarsest level with either an iterative or direct solver
 216     solveCoarsestLevel
 217     (
 218         coarseCorrX[coarsestLevel],
 219         coarseB[coarsestLevel]
 220     );
 221
 222
 223     if (debug >= 2)
 224     {
 225         Pout<< "Post-smoothing scaling factors: ";
 226     }
 227
 228     // Smoothing and prolongation of the coarse correction fields
 229     // (going to finer levels)
 230     for (label leveli = coarsestLevel - 1; leveli >= 0; leveli--)
 231     {
 232         // Create a field for the pre-smoothed correction field
 233         // as a sub-field of the finestCorrection which is not
 234         // currently being used
 235         scalarField::subField preSmoothedCoarseCorrField
 236         (
 237             finestCorrection,
 238             coarseCorrX[leveli].size()
 239         );
 240
 241         // Only store the preSmoothedCoarseCorrField is pre-smoothing is used
 242         if (nPreSweeps_)
 243         {
 244             preSmoothedCoarseCorrField.assign(coarseCorrX[leveli]);
 245         }
 246
 247         agglomeration_.prolongField
 248         (
 249             coarseCorrX[leveli],
 250             coarseCorrX[leveli + 1],
 251             leveli + 1
 252         );
 253
 254         // Scale coarse-grid correction field
 255         // but not on the coarsest level because it evaluates to 1
 256         if (scaleCorrection_ && leveli < coarsestLevel - 1)
 257         {
 258             // Create A.x for this coarse level as a sub-field of Ax
 259             scalarField::subField ACf
 260             (
 261                 Ax,
 262                 coarseCorrX[leveli].size()
 263             );
 264
 265             scalar sf = scalingFactor
 266             (
 267                 const_cast<scalarField&>(ACf.operator const scalarField&()),
 268                 matrixLevels_[leveli],
 269                 coarseCorrX[leveli],
 270                 coupleLevelsBouCoeffs_[leveli],
 271                 interfaceLevels_[leveli],
 272                 coarseB[leveli],
 273                 cmpt
 274             );
 275
 276
 277             if (debug >= 2)
 278             {
 279                 Pout<< sf << " ";
 280             }
 281
 282             coarseCorrX[leveli] *= sf;
 283         }
 284
 285         // Only add the preSmoothedCoarseCorrField is pre-smoothing is used
 286         if (nPreSweeps_)
 287         {
 288             coarseCorrX[leveli] += preSmoothedCoarseCorrField;
 289         }
 290
 291         smoothers[leveli + 1].smooth
 292         (
 293             coarseCorrX[leveli],
 294             coarseB[leveli],
 295             cmpt,
 296             nPostSweeps_ + leveli
 297         );
 298     }
 299
 300     // Prolong the finest level correction
 301     agglomeration_.prolongField
 302     (
 303         finestCorrection,
 304         coarseCorrX[0],
 305         0
 306     );
 307
 308     if (scaleCorrection_)
 309     {
 310         // Calculate finest level scaling factor
 311         scalar fsf = scalingFactor
 312         (
 313             Ax,
 314             matrix_,
 315             finestCorrection,
 316             coupleBouCoeffs_,
 317             interfaces_,
 318             finestResidual,
 319             cmpt
 320         );
 321
 322         if (debug >= 2)
 323         {
 324             Pout<< fsf << endl;
 325         }
 326
 327         forAll(x, i)
 328         {
 329             x[i] += fsf*finestCorrection[i];
 330         }
 331     }
 332     else
 333     {
 334         forAll(x, i)
 335         {
 336             x[i] += finestCorrection[i];
 337         }
 338     }
 339
 340     smoothers[0].smooth
 341     (
 342         x,
 343         b,
 344         cmpt,
 345         nFinestSweeps_
 346     );
 347 }
 348
 349
 350 void Foam::GAMGSolver::initVcycle
 351 (
 352     PtrList<scalarField>& coarseCorrX,
 353     PtrList<scalarField>& coarseB,
 354     PtrList<lduSmoother>& smoothers
 355 ) const
 356 {
 357     coarseCorrX.setSize(matrixLevels_.size());
 358     coarseB.setSize(matrixLevels_.size());
 359     smoothers.setSize(matrixLevels_.size() + 1);
 360
 361     // Create the smoother for the finest level
 362     smoothers.set
 363     (
 364         0,
 365         lduSmoother::New
 366         (
 367             matrix_,
 368             coupleBouCoeffs_,
 369             coupleIntCoeffs_,
 370             interfaces_,
 371             dict()
 372         )
 373     );
 374
 375     forAll (matrixLevels_, leveli)
 376     {
 377         coarseCorrX.set
 378         (
 379             leveli,
 380             new scalarField
 381             (
 382                 agglomeration_.meshLevel(leveli + 1).lduAddr().size()
 383             )
 384         );
 385
 386         coarseB.set
 387         (
 388             leveli,
 389             new scalarField
 390             (
 391                 agglomeration_.meshLevel(leveli + 1).lduAddr().size()
 392             )
 393         );
 394
 395         smoothers.set
 396         (
 397             leveli + 1,
 398             lduSmoother::New
 399             (
 400                 matrixLevels_[leveli],
 401                 coupleLevelsBouCoeffs_[leveli],
 402                 coupleLevelsIntCoeffs_[leveli],
 403                 interfaceLevels_[leveli],
 404                 dict()
 405             )
 406         );
 407     }
 408 }
 409
 410
 411 void Foam::GAMGSolver::solveCoarsestLevel
 412 (
 413     scalarField& coarsestCorrX,
 414     const scalarField& coarsestB
 415 ) const
 416 {
 417     if (directSolveCoarsest_)
 418     {
 419         coarsestCorrX = coarsestB;
 420         coarsestLUMatrixPtr_->solve(coarsestCorrX);
 421     }
 422     else
 423     {
 424         const label coarsestLevel = matrixLevels_.size() - 1;
 425         coarsestCorrX = 0;
 426         lduSolverPerformance coarseSolverPerf;
 427
 428         if (matrixLevels_[coarsestLevel].asymmetric())
 429         {
 430             coarseSolverPerf = BICCG
 431             (
 432                 "coarsestLevelCorr",
 433                 matrixLevels_[coarsestLevel],
 434                 coupleLevelsBouCoeffs_[coarsestLevel],
 435                 coupleLevelsIntCoeffs_[coarsestLevel],
 436                 interfaceLevels_[coarsestLevel],
 437                 tolerance(),
 438                 relTolerance()
 439             ).solve
 440             (
 441                 coarsestCorrX,
 442                 coarsestB
 443             );
 444         }
 445         else
 446         {
 447             coarseSolverPerf = ICCG
 448             (
 449                 "coarsestLevelCorr",
 450                 matrixLevels_[coarsestLevel],
 451                 coupleLevelsBouCoeffs_[coarsestLevel],
 452                 coupleLevelsIntCoeffs_[coarsestLevel],
 453                 interfaceLevels_[coarsestLevel],
 454                 tolerance(),
 455                 relTolerance()
 456             ).solve
 457             (
 458                 coarsestCorrX,
 459                 coarsestB
 460             );
 461         }
 462
 463         if (debug >= 2)
 464         {
 465             coarseSolverPerf.print();
 466         }
 467     }
 468 }
 469
 470
 471 // ************************************************************************* //