src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 "lduMatrix.H"
  28 #include "diagonalSolver.H"
  29
  30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34     defineRunTimeSelectionTable(lduMatrix::solver, symMatrix);
  35     defineRunTimeSelectionTable(lduMatrix::solver, asymMatrix);
  36 }
  37
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 Foam::autoPtr<Foam::lduMatrix::solver> Foam::lduMatrix::solver::New
  42 (
  43     const word& fieldName,
  44     const lduMatrix& matrix,
  45     const FieldField<Field, scalar>& interfaceBouCoeffs,
  46     const FieldField<Field, scalar>& interfaceIntCoeffs,
  47     const lduInterfaceFieldPtrsList& interfaces,
  48     const dictionary& solverControls
  49 )
  50 {
  51     word name(solverControls.lookup("solver"));
  52
  53     if (matrix.diagonal())
  54     {
  55         return autoPtr<lduMatrix::solver>
  56         (
  57             new diagonalSolver
  58             (
  59                 fieldName,
  60                 matrix,
  61                 interfaceBouCoeffs,
  62                 interfaceIntCoeffs,
  63                 interfaces,
  64                 solverControls
  65             )
  66         );
  67     }
  68     else if (matrix.symmetric())
  69     {
  70         symMatrixConstructorTable::iterator constructorIter =
  71             symMatrixConstructorTablePtr_->find(name);
  72
  73         if (constructorIter == symMatrixConstructorTablePtr_->end())
  74         {
  75             FatalIOErrorIn
  76             (
  77                 "lduMatrix::solver::New", solverControls
  78             )   << "Unknown symmetric matrix solver " << name << nl << nl
  79                 << "Valid symmetric matrix solvers are :" << endl
  80                 << symMatrixConstructorTablePtr_->toc()
  81                 << exit(FatalIOError);
  82         }
  83
  84         return autoPtr<lduMatrix::solver>
  85         (
  86             constructorIter()
  87             (
  88                 fieldName,
  89                 matrix,
  90                 interfaceBouCoeffs,
  91                 interfaceIntCoeffs,
  92                 interfaces,
  93                 solverControls
  94             )
  95         );
  96     }
  97     else if (matrix.asymmetric())
  98     {
  99         asymMatrixConstructorTable::iterator constructorIter =
 100             asymMatrixConstructorTablePtr_->find(name);
 101
 102         if (constructorIter == asymMatrixConstructorTablePtr_->end())
 103         {
 104             FatalIOErrorIn
 105             (
 106                 "lduMatrix::solver::New", solverControls
 107             )   << "Unknown asymmetric matrix solver " << name << nl << nl
 108                 << "Valid asymmetric matrix solvers are :" << endl
 109                 << asymMatrixConstructorTablePtr_->toc()
 110                 << exit(FatalIOError);
 111         }
 112
 113         return autoPtr<lduMatrix::solver>
 114         (
 115             constructorIter()
 116             (
 117                 fieldName,
 118                 matrix,
 119                 interfaceBouCoeffs,
 120                 interfaceIntCoeffs,
 121                 interfaces,
 122                 solverControls
 123             )
 124         );
 125     }
 126     else
 127     {
 128         FatalIOErrorIn
 129         (
 130             "lduMatrix::solver::New", solverControls
 131         )   << "cannot solve incomplete matrix, "
 132                "no diagonal or off-diagonal coefficient"
 133             << exit(FatalIOError);
 134
 135         return autoPtr<lduMatrix::solver>(NULL);
 136     }
 137 }
 138
 139
 140 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 141
 142 Foam::lduMatrix::solver::solver
 143 (
 144     const word& fieldName,
 145     const lduMatrix& matrix,
 146     const FieldField<Field, scalar>& interfaceBouCoeffs,
 147     const FieldField<Field, scalar>& interfaceIntCoeffs,
 148     const lduInterfaceFieldPtrsList& interfaces,
 149     const dictionary& solverControls
 150 )
 151 :
 152     fieldName_(fieldName),
 153     matrix_(matrix),
 154     interfaceBouCoeffs_(interfaceBouCoeffs),
 155     interfaceIntCoeffs_(interfaceIntCoeffs),
 156     interfaces_(interfaces),
 157     controlDict_(solverControls)
 158 {
 159     readControls();
 160 }
 161
 162
 163 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 164
 165 void Foam::lduMatrix::solver::readControls()
 166 {
 167     maxIter_   = controlDict_.lookupOrDefault<label>("maxIter", 1000);
 168     tolerance_ = controlDict_.lookupOrDefault<scalar>("tolerance", 1e-6);
 169     relTol_    = controlDict_.lookupOrDefault<scalar>("relTol", 0);
 170 }
 171
 172
 173 void Foam::lduMatrix::solver::read(const dictionary& solverControls)
 174 {
 175     controlDict_ = solverControls;
 176     readControls();
 177 }
 178
 179
 180 Foam::scalar Foam::lduMatrix::solver::normFactor
 181 (
 182     const scalarField& psi,
 183     const scalarField& source,
 184     const scalarField& Apsi,
 185     scalarField& tmpField
 186 ) const
 187 {
 188     // --- Calculate A dot reference value of psi
 189     matrix_.sumA(tmpField, interfaceBouCoeffs_, interfaces_);
 190     tmpField *= gAverage(psi);
 191
 192     return gSum(mag(Apsi - tmpField) + mag(source - tmpField)) + matrix_.small_;
 193
 194     // At convergence this simpler method is equivalent to the above
 195     // return 2*gSumMag(source) + matrix_.small_;
 196 }
 197
 198
 199 // ************************************************************************* //