src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 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 \*---------------------------------------------------------------------------*/
  26
  27 #include "smoothSolver.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33     defineTypeNameAndDebug(smoothSolver, 0);
  34
  35     lduMatrix::solver::addsymMatrixConstructorToTable<smoothSolver>
  36         addsmoothSolverSymMatrixConstructorToTable_;
  37
  38     lduMatrix::solver::addasymMatrixConstructorToTable<smoothSolver>
  39         addsmoothSolverAsymMatrixConstructorToTable_;
  40 }
  41
  42
  43 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  44
  45 Foam::smoothSolver::smoothSolver
  46 (
  47     const word& fieldName,
  48     const lduMatrix& matrix,
  49     const FieldField<Field, scalar>& interfaceBouCoeffs,
  50     const FieldField<Field, scalar>& interfaceIntCoeffs,
  51     const lduInterfaceFieldPtrsList& interfaces,
  52     Istream& solverData
  53 )
  54 :
  55     lduMatrix::solver
  56     (
  57         fieldName,
  58         matrix,
  59         interfaceBouCoeffs,
  60         interfaceIntCoeffs,
  61         interfaces,
  62         solverData
  63     ),
  64     nSweeps_(1)
  65 {
  66     readControls();
  67 }
  68
  69
  70 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  71
  72 void Foam::smoothSolver::readControls()
  73 {
  74     lduMatrix::solver::readControls();
  75     controlDict_.readIfPresent("nSweeps", nSweeps_);
  76 }
  77
  78
  79 Foam::lduMatrix::solverPerformance Foam::smoothSolver::solve
  80 (
  81     scalarField& psi,
  82     const scalarField& source,
  83     const direction cmpt
  84 ) const
  85 {
  86     // Setup class containing solver performance data
  87     lduMatrix::solverPerformance solverPerf(typeName, fieldName_);
  88
  89     // If the nSweeps_ is negative do a fixed number of sweeps
  90     if (nSweeps_ < 0)
  91     {
  92         autoPtr<lduMatrix::smoother> smootherPtr = lduMatrix::smoother::New
  93         (
  94             fieldName_,
  95             matrix_,
  96             interfaceBouCoeffs_,
  97             interfaceIntCoeffs_,
  98             interfaces_,
  99             controlDict_.lookup("smoother")
 100         );
 101
 102         smootherPtr->smooth
 103         (
 104             psi,
 105             source,
 106             cmpt,
 107             -nSweeps_
 108         );
 109
 110         solverPerf.nIterations() -= nSweeps_;
 111     }
 112     else
 113     {
 114         scalar normFactor = 0;
 115
 116         {
 117             scalarField Apsi(psi.size());
 118             scalarField temp(psi.size());
 119
 120             // Calculate A.psi
 121             matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
 122
 123             // Calculate normalisation factor
 124             normFactor = this->normFactor(psi, source, Apsi, temp);
 125
 126             // Calculate residual magnitude
 127             solverPerf.initialResidual() = gSumMag(source - Apsi)/normFactor;
 128             solverPerf.finalResidual() = solverPerf.initialResidual();
 129         }
 130
 131         if (lduMatrix::debug >= 2)
 132         {
 133             Info<< "   Normalisation factor = " << normFactor << endl;
 134         }
 135
 136
 137         // Check convergence, solve if not converged
 138         if (!solverPerf.checkConvergence(tolerance_, relTol_))
 139         {
 140             autoPtr<lduMatrix::smoother> smootherPtr = lduMatrix::smoother::New
 141             (
 142                 fieldName_,
 143                 matrix_,
 144                 interfaceBouCoeffs_,
 145                 interfaceIntCoeffs_,
 146                 interfaces_,
 147                 controlDict_.lookup("smoother")
 148             );
 149
 150             // Smoothing loop
 151             do
 152             {
 153                 smootherPtr->smooth
 154                 (
 155                     psi,
 156                     source,
 157                     cmpt,
 158                     nSweeps_
 159                 );
 160
 161                 // Calculate the residual to check convergence
 162                 solverPerf.finalResidual() = gSumMag
 163                 (
 164                     matrix_.residual
 165                     (
 166                         psi,
 167                         source,
 168                         interfaceBouCoeffs_,
 169                         interfaces_,
 170                         cmpt
 171                     )
 172                 )/normFactor;
 173             } while
 174             (
 175                 (solverPerf.nIterations() += nSweeps_) < maxIter_
 176              && !(solverPerf.checkConvergence(tolerance_, relTol_))
 177             );
 178         }
 179     }
 180
 181     return solverPerf;
 182 }
 183
 184
 185 // ************************************************************************* //