src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  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& psi,
  36     const scalarField& source,
  37     const direction cmpt
  38 ) const
  39 {
  40     // Setup class containing solver performance data
  41     lduMatrix::solverPerformance solverPerf(typeName, fieldName_);
  42
  43     // Calculate A.psi used to calculate the initial residual
  44     scalarField Apsi(psi.size());
  45     matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
  46
  47     // Create the storage for the finestCorrection which may be used as a
  48     // temporary in normFactor
  49     scalarField finestCorrection(psi.size());
  50
  51     // Calculate normalisation factor
  52     scalar normFactor = this->normFactor(psi, source, Apsi, finestCorrection);
  53
  54     if (debug >= 2)
  55     {
  56         Pout<< "   Normalisation factor = " << normFactor << endl;
  57     }
  58
  59     // Calculate initial finest-grid residual field
  60     scalarField finestResidual(source - Apsi);
  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 (!solverPerf.checkConvergence(tolerance_, relTol_))
  69     {
  70         // Create coarse grid correction fields
  71         PtrList<scalarField> coarseCorrFields;
  72
  73         // Create coarse grid sources
  74         PtrList<scalarField> coarseSources;
  75
  76         // Create the smoothers for all levels
  77         PtrList<lduMatrix::smoother> smoothers;
  78
  79         // Initialise the above data structures
  80         initVcycle(coarseCorrFields, coarseSources, smoothers);
  81
  82         do
  83         {
  84             Vcycle
  85             (
  86                 smoothers,
  87                 psi,
  88                 source,
  89                 Apsi,
  90                 finestCorrection,
  91                 finestResidual,
  92                 coarseCorrFields,
  93                 coarseSources,
  94                 cmpt
  95             );
  96
  97             // Calculate finest level residual field
  98             matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
  99             finestResidual = source;
 100             finestResidual -= Apsi;
 101
 102             solverPerf.finalResidual() = gSumMag(finestResidual)/normFactor;
 103
 104             if (debug >= 2)
 105             {
 106                 solverPerf.print();
 107             }
 108         } while
 109         (
 110             ++solverPerf.nIterations() < maxIter_
 111          && !(solverPerf.checkConvergence(tolerance_, relTol_))
 112         );
 113     }
 114
 115     return solverPerf;
 116 }
 117
 118
 119 void Foam::GAMGSolver::Vcycle
 120 (
 121     const PtrList<lduMatrix::smoother>& smoothers,
 122     scalarField& psi,
 123     const scalarField& source,
 124     scalarField& Apsi,
 125     scalarField& finestCorrection,
 126     scalarField& finestResidual,
 127     PtrList<scalarField>& coarseCorrFields,
 128     PtrList<scalarField>& coarseSources,
 129     const direction cmpt
 130 ) const
 131 {
 132     //debug = 2;
 133
 134     const label coarsestLevel = matrixLevels_.size() - 1;
 135
 136     // Restrict finest grid residual for the next level up
 137     agglomeration_.restrictField(coarseSources[0], finestResidual, 0);
 138
 139     if (debug >= 2 && nPreSweeps_)
 140     {
 141         Pout<< "Pre-smoothing scaling factors: ";
 142     }
 143
 144
 145     // Residual restriction (going to coarser levels)
 146     for (label leveli = 0; leveli < coarsestLevel; leveli++)
 147     {
 148         // If the optional pre-smoothing sweeps are selected
 149         // smooth the coarse-grid field for the restriced source
 150         if (nPreSweeps_)
 151         {
 152             coarseCorrFields[leveli] = 0.0;
 153
 154             smoothers[leveli + 1].smooth
 155             (
 156                 coarseCorrFields[leveli],
 157                 coarseSources[leveli],
 158                 cmpt,
 159                 nPreSweeps_ + leveli
 160             );
 161
 162             scalarField::subField ACf
 163             (
 164                 Apsi,
 165                 coarseCorrFields[leveli].size()
 166             );
 167
 168             // Scale coarse-grid correction field
 169             // but not on the coarsest level because it evaluates to 1
 170             if (scaleCorrection_ && leveli < coarsestLevel - 1)
 171             {
 172                 scalar sf = scalingFactor
 173                 (
 174                     const_cast<scalarField&>(ACf.operator const scalarField&()),
 175                     matrixLevels_[leveli],
 176                     coarseCorrFields[leveli],
 177                     interfaceLevelsBouCoeffs_[leveli],
 178                     interfaceLevels_[leveli],
 179                     coarseSources[leveli],
 180                     cmpt
 181                 );
 182
 183                 if (debug >= 2)
 184                 {
 185                     Pout<< sf << " ";
 186                 }
 187
 188                 coarseCorrFields[leveli] *= sf;
 189             }
 190
 191             // Correct the residual with the new solution
 192             matrixLevels_[leveli].Amul
 193             (
 194                 const_cast<scalarField&>(ACf.operator const scalarField&()),
 195                 coarseCorrFields[leveli],
 196                 interfaceLevelsBouCoeffs_[leveli],
 197                 interfaceLevels_[leveli],
 198                 cmpt
 199             );
 200
 201             coarseSources[leveli] -= ACf;
 202         }
 203
 204         // Residual is equal to source
 205         agglomeration_.restrictField
 206         (
 207             coarseSources[leveli + 1],
 208             coarseSources[leveli],
 209             leveli + 1
 210         );
 211     }
 212
 213     if (debug >= 2 && nPreSweeps_)
 214     {
 215         Pout<< endl;
 216     }
 217
 218
 219     // Solve Coarsest level with either an iterative or direct solver
 220     solveCoarsestLevel
 221     (
 222         coarseCorrFields[coarsestLevel],
 223         coarseSources[coarsestLevel]
 224     );
 225
 226
 227     if (debug >= 2)
 228     {
 229         Pout<< "Post-smoothing scaling factors: ";
 230     }
 231
 232     // Smoothing and prolongation of the coarse correction fields
 233     // (going to finer levels)
 234     for (label leveli = coarsestLevel - 1; leveli >= 0; leveli--)
 235     {
 236         // Create a field for the pre-smoothed correction field
 237         // as a sub-field of the finestCorrection which is not
 238         // currently being used
 239         scalarField::subField preSmoothedCoarseCorrField
 240         (
 241             finestCorrection,
 242             coarseCorrFields[leveli].size()
 243         );
 244
 245         // Only store the preSmoothedCoarseCorrField is pre-smoothing is used
 246         if (nPreSweeps_)
 247         {
 248             preSmoothedCoarseCorrField.assign(coarseCorrFields[leveli]);
 249         }
 250
 251         agglomeration_.prolongField
 252         (
 253             coarseCorrFields[leveli],
 254             coarseCorrFields[leveli + 1],
 255             leveli + 1
 256         );
 257
 258         // Scale coarse-grid correction field
 259         // but not on the coarsest level because it evaluates to 1
 260         if (scaleCorrection_ && leveli < coarsestLevel - 1)
 261         {
 262             // Create A.psi for this coarse level as a sub-field of Apsi
 263             scalarField::subField ACf
 264             (
 265                 Apsi,
 266                 coarseCorrFields[leveli].size()
 267             );
 268
 269             scalar sf = scalingFactor
 270             (
 271                 const_cast<scalarField&>(ACf.operator const scalarField&()),
 272                 matrixLevels_[leveli],
 273                 coarseCorrFields[leveli],
 274                 interfaceLevelsBouCoeffs_[leveli],
 275                 interfaceLevels_[leveli],
 276                 coarseSources[leveli],
 277                 cmpt
 278             );
 279
 280
 281             if (debug >= 2)
 282             {
 283                 Pout<< sf << " ";
 284             }
 285
 286             coarseCorrFields[leveli] *= sf;
 287         }
 288
 289         // Only add the preSmoothedCoarseCorrField is pre-smoothing is used
 290         if (nPreSweeps_)
 291         {
 292             coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
 293         }
 294
 295         smoothers[leveli + 1].smooth
 296         (
 297             coarseCorrFields[leveli],
 298             coarseSources[leveli],
 299             cmpt,
 300             nPostSweeps_ + leveli
 301         );
 302     }
 303
 304     // Prolong the finest level correction
 305     agglomeration_.prolongField
 306     (
 307         finestCorrection,
 308         coarseCorrFields[0],
 309         0
 310     );
 311
 312     if (scaleCorrection_)
 313     {
 314         // Calculate finest level scaling factor
 315         scalar fsf = scalingFactor
 316         (
 317             Apsi,
 318             matrix_,
 319             finestCorrection,
 320             interfaceBouCoeffs_,
 321             interfaces_,
 322             finestResidual,
 323             cmpt
 324         );
 325
 326         if (debug >= 2)
 327         {
 328             Pout<< fsf << endl;
 329         }
 330
 331         forAll(psi, i)
 332         {
 333             psi[i] += fsf*finestCorrection[i];
 334         }
 335     }
 336     else
 337     {
 338         forAll(psi, i)
 339         {
 340             psi[i] += finestCorrection[i];
 341         }
 342     }
 343
 344     smoothers[0].smooth
 345     (
 346         psi,
 347         source,
 348         cmpt,
 349         nFinestSweeps_
 350     );
 351 }
 352
 353
 354 void Foam::GAMGSolver::initVcycle
 355 (
 356     PtrList<scalarField>& coarseCorrFields,
 357     PtrList<scalarField>& coarseSources,
 358     PtrList<lduMatrix::smoother>& smoothers
 359 ) const
 360 {
 361     coarseCorrFields.setSize(matrixLevels_.size());
 362     coarseSources.setSize(matrixLevels_.size());
 363     smoothers.setSize(matrixLevels_.size() + 1);
 364
 365     // Create the smoother for the finest level
 366     smoothers.set
 367     (
 368         0,
 369         lduMatrix::smoother::New
 370         (
 371             fieldName_,
 372             matrix_,
 373             interfaceBouCoeffs_,
 374             interfaceIntCoeffs_,
 375             interfaces_,
 376             controlDict_
 377         )
 378     );
 379
 380     forAll(matrixLevels_, leveli)
 381     {
 382         coarseCorrFields.set
 383         (
 384             leveli,
 385             new scalarField
 386             (
 387                 agglomeration_.meshLevel(leveli + 1).lduAddr().size()
 388             )
 389         );
 390
 391         coarseSources.set
 392         (
 393             leveli,
 394             new scalarField
 395             (
 396                 agglomeration_.meshLevel(leveli + 1).lduAddr().size()
 397             )
 398         );
 399
 400         smoothers.set
 401         (
 402             leveli + 1,
 403             lduMatrix::smoother::New
 404             (
 405                 fieldName_,
 406                 matrixLevels_[leveli],
 407                 interfaceLevelsBouCoeffs_[leveli],
 408                 interfaceLevelsIntCoeffs_[leveli],
 409                 interfaceLevels_[leveli],
 410                 controlDict_
 411             )
 412         );
 413     }
 414 }
 415
 416
 417 void Foam::GAMGSolver::solveCoarsestLevel
 418 (
 419     scalarField& coarsestCorrField,
 420     const scalarField& coarsestSource
 421 ) const
 422 {
 423     if (directSolveCoarsest_)
 424     {
 425         coarsestCorrField = coarsestSource;
 426         coarsestLUMatrixPtr_->solve(coarsestCorrField);
 427     }
 428     else
 429     {
 430         const label coarsestLevel = matrixLevels_.size() - 1;
 431         coarsestCorrField = 0;
 432         lduMatrix::solverPerformance coarseSolverPerf;
 433
 434         if (matrixLevels_[coarsestLevel].asymmetric())
 435         {
 436             coarseSolverPerf = BICCG
 437             (
 438                 "coarsestLevelCorr",
 439                 matrixLevels_[coarsestLevel],
 440                 interfaceLevelsBouCoeffs_[coarsestLevel],
 441                 interfaceLevelsIntCoeffs_[coarsestLevel],
 442                 interfaceLevels_[coarsestLevel],
 443                 tolerance_,
 444                 relTol_
 445             ).solve
 446             (
 447                 coarsestCorrField,
 448                 coarsestSource
 449             );
 450         }
 451         else
 452         {
 453             coarseSolverPerf = ICCG
 454             (
 455                 "coarsestLevelCorr",
 456                 matrixLevels_[coarsestLevel],
 457                 interfaceLevelsBouCoeffs_[coarsestLevel],
 458                 interfaceLevelsIntCoeffs_[coarsestLevel],
 459                 interfaceLevels_[coarsestLevel],
 460                 tolerance_,
 461                 relTol_
 462             ).solve
 463             (
 464                 coarsestCorrField,
 465                 coarsestSource
 466             );
 467         }
 468
 469         if (debug >= 2)
 470         {
 471             coarseSolverPerf.print();
 472         }
 473     }
 474 }
 475
 476
 477 // ************************************************************************* //