src/lduSolvers/lduSolver/mpeAmgSolver/mpeAmgSolver.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   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 Description
  25     Minimal Polynomial Extrapolation Algebraic Multigrid solver with run-time
  26     selection of policy and cycle
  27
  28 Author
  29     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "mpeAmgSolver.H"
  34 #include "scalarMatrices.H"
  35 #include "DenseMatrixTools.H"
  36 #include "FieldFields.H"
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42
  43     defineTypeNameAndDebug(mpeAmgSolver, 0);
  44
  45     lduSolver::addsymMatrixConstructorToTable<mpeAmgSolver>
  46         addmpeAmgSolverSymMatrixConstructorToTable_;
  47
  48     lduSolver::addasymMatrixConstructorToTable<mpeAmgSolver>
  49         addmpeAmgSolverAsymMatrixConstructorToTable_;
  50 }
  51
  52
  53 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  54
  55 //- Construct from matrix and solver data stream
  56 Foam::mpeAmgSolver::mpeAmgSolver
  57 (
  58     const word& fieldName,
  59     const lduMatrix& matrix,
  60     const FieldField<Field, scalar>& coupleBouCoeffs,
  61     const FieldField<Field, scalar>& coupleIntCoeffs,
  62     const lduInterfaceFieldPtrsList& interfaces,
  63     const dictionary& dict
  64 )
  65 :
  66     lduSolver
  67     (
  68         fieldName,
  69         matrix,
  70         coupleBouCoeffs,
  71         coupleIntCoeffs,
  72         interfaces,
  73         dict
  74     ),
  75     amg_
  76     (
  77         matrix,
  78         coupleBouCoeffs,
  79         coupleIntCoeffs,
  80         interfaces,
  81         dict
  82     ),
  83     kDimension_(readLabel(dict.lookup("kDimension")))
  84 {}
  85
  86
  87 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  88
  89 Foam::lduSolverPerformance Foam::mpeAmgSolver::solve
  90 (
  91     scalarField& x,
  92     const scalarField& b,
  93     const direction cmpt
  94 ) const
  95 {
  96     // Prepare solver performance
  97     lduSolverPerformance solverPerf(typeName, fieldName());
  98
  99     scalar normFactor = this->normFactor(x, b, cmpt);
 100
 101     // Calculate initial residual
 102     solverPerf.initialResidual() =
 103         gSumMag(amg_.residual(x, b, cmpt))/normFactor;
 104
 105     solverPerf.finalResidual() = solverPerf.initialResidual();
 106
 107
 108     // Solver loop
 109
 110     if (!stop(solverPerf))
 111     {
 112         // Krylov vectors
 113         typedef FieldField<Field, scalar> scalarFieldField;
 114
 115         scalarFieldField xSave(kDimension_ + 1);
 116
 117         forAll (xSave, i)
 118         {
 119             xSave.set(i, new scalarField(x.size()));
 120         }
 121
 122         scalarFieldField xFirstDiffs(kDimension_ - 1);
 123
 124         forAll (xFirstDiffs, i)
 125         {
 126             xFirstDiffs.set(i, new scalarField(x.size()));
 127         }
 128
 129         scalarField xRhs(x.size());
 130
 131         // Matrices
 132         scalarSquareMatrix AN(kDimension_ - 1);
 133         scalarField xN(kDimension_ - 1);
 134         scalarField bN(kDimension_ - 1);
 135
 136         // Remember last extrapolation index
 137         label lastExtrapolation = kDimension_;
 138
 139         scalarField psiSave(x.size());
 140
 141         do
 142         {
 143             amg_.cycle(x, b, cmpt);
 144
 145             label curIndex = (solverPerf.nIterations()) % (kDimension_ + 1);
 146
 147             // Save x into xSave
 148             xSave[curIndex] = x;
 149
 150             if
 151             (
 152                 solverPerf.nIterations() >= lastExtrapolation
 153              && curIndex == kDimension_
 154             )
 155             {
 156                 // Calculate differences
 157                 for (label i = 0; i < kDimension_ - 1; i++)
 158                 {
 159                     xFirstDiffs[i] = xSave[i] - xSave[kDimension_];
 160                 }
 161
 162                 // Grab right-hand side
 163                 xRhs = xSave[kDimension_] - xSave[kDimension_ - 1];
 164
 165                 // Perform Q-R decomposition
 166                 DenseMatrixTools::qrDecompose
 167                 (
 168                     kDimension_ - 1,
 169                     xFirstDiffs,
 170                     AN
 171                 );
 172
 173                 // Solve the normal system
 174                 for (label i = 0; i < kDimension_ - 1; i++)
 175                 {
 176                     bN[i] = -gSum(xFirstDiffs[i]*xRhs);
 177                 }
 178
 179                 DenseMatrixTools::solve(AN, xN, bN);
 180
 181                 // Reset last interpolation factor
 182                 xN[kDimension_ - 2] = 1.0;
 183
 184                 // Re-normalize
 185                 xN /= sum(xN);
 186
 187
 188                 // Re-asemble the solution.  HJ, variants?
 189                 x = xSave[kDimension_];
 190
 191                 for (label i = 0; i < kDimension_ - 1; i++)
 192                 {
 193                     xFirstDiffs[0] = xSave[i + 1] - xSave[i];
 194                     x += xN[i]*xFirstDiffs[0];
 195                 }
 196
 197                 // Increment last extrapolation
 198                 lastExtrapolation = solverPerf.nIterations() + kDimension_;
 199             }
 200
 201             // Re-calculate residual
 202             solverPerf.finalResidual() =
 203                 gSumMag(amg_.residual(x, b, cmpt))/normFactor;
 204             solverPerf.nIterations()++;
 205         } while (!stop(solverPerf));
 206     }
 207
 208     return solverPerf;
 209 }
 210
 211
 212 // ************************************************************************* //