src/lduSolvers/lduSolver/rreAmgSolver/rreAmgSolver.C

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