LduMatrixSolver.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
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 "LduMatrix.H"
27 #include "DiagonalSolver.H"
28 
29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
30 
31 template<class Type, class DType, class LUType>
34 (
35  const word& fieldName,
36  const LduMatrix<Type, DType, LUType>& matrix,
37  const dictionary& solverDict
38 )
39 {
40  word solverName = solverDict.lookup("solver");
41 
42  if (matrix.diagonal())
43  {
45  (
47  (
48  fieldName,
49  matrix,
50  solverDict
51  )
52  );
53  }
54  else if (matrix.symmetric())
55  {
56  typename symMatrixConstructorTable::iterator constructorIter =
57  symMatrixConstructorTablePtr_->find(solverName);
58 
59  if (constructorIter == symMatrixConstructorTablePtr_->end())
60  {
61  FatalIOErrorInFunction(solverDict)
62  << "Unknown symmetric matrix solver " << solverName
63  << endl << endl
64  << "Valid symmetric matrix solvers are :" << endl
65  << symMatrixConstructorTablePtr_->toc()
66  << exit(FatalIOError);
67  }
68 
70  (
71  constructorIter()
72  (
73  fieldName,
74  matrix,
75  solverDict
76  )
77  );
78  }
79  else if (matrix.asymmetric())
80  {
81  typename asymMatrixConstructorTable::iterator constructorIter =
82  asymMatrixConstructorTablePtr_->find(solverName);
83 
84  if (constructorIter == asymMatrixConstructorTablePtr_->end())
85  {
86  FatalIOErrorInFunction(solverDict)
87  << "Unknown asymmetric matrix solver " << solverName
88  << endl << endl
89  << "Valid asymmetric matrix solvers are :" << endl
90  << asymMatrixConstructorTablePtr_->toc()
91  << exit(FatalIOError);
92  }
93 
95  (
96  constructorIter()
97  (
98  fieldName,
99  matrix,
100  solverDict
101  )
102  );
103  }
104  else
105  {
106  FatalIOErrorInFunction(solverDict)
107  << "cannot solve incomplete matrix, "
108  "no diagonal or off-diagonal coefficient"
109  << exit(FatalIOError);
110 
112  (
113  nullptr
114  );
115  }
116 }
117 
118 
119 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
120 
121 template<class Type, class DType, class LUType>
123 (
124  const word& fieldName,
125  const LduMatrix<Type, DType, LUType>& matrix,
126  const dictionary& solverDict
127 )
128 :
129  fieldName_(fieldName),
130  matrix_(matrix),
131 
132  controlDict_(solverDict),
133 
134  maxIter_(defaultMaxIter_),
135  minIter_(0),
136  tolerance_(1e-6*pTraits<Type>::one),
137  relTol_(Zero)
138 {
139  readControls();
140 }
141 
142 
143 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
144 
145 template<class Type, class DType, class LUType>
147 {
148  readControl(controlDict_, maxIter_, "maxIter");
149  readControl(controlDict_, minIter_, "minIter");
150  readControl(controlDict_, tolerance_, "tolerance");
151  readControl(controlDict_, relTol_, "relTol");
152 }
153 
154 
155 template<class Type, class DType, class LUType>
157 (
158  const dictionary& solverDict
159 )
160 {
161  controlDict_ = solverDict;
162  readControls();
163 }
164 
165 
166 template<class Type, class DType, class LUType>
168 (
169  const Field<Type>& psi,
170  const Field<Type>& Apsi,
171  Field<Type>& tmpField
172 ) const
173 {
174  // --- Calculate A dot reference value of psi
175  matrix_.sumA(tmpField);
176  cmptMultiply(tmpField, tmpField, gAverage(psi));
177 
178  return stabilise
179  (
180  gSum(cmptMag(Apsi - tmpField) + cmptMag(matrix_.source() - tmpField)),
182  );
183 
184  // At convergence this simpler method is equivalent to the above
185  // return stabilise(2*gSumCmptMag(matrix_.source()), matrix_.small_);
186 }
187 
188 
189 // ************************************************************************* //
Foam::DiagonalSolver.
Pre-declare SubField and related Field type.
Definition: Field.H:83
solver(const word &fieldName, const LduMatrix< Type, DType, LUType > &matrix, const dictionary &solverDict)
const LduMatrix< Type, DType, LUType > & matrix() const
Definition: LduMatrix.H:235
virtual void read(const dictionary &solverDict)
Read and reset the solver parameters from the given dictionary.
Type normFactor(const Field< Type > &psi, const Field< Type > &Apsi, Field< Type > &tmpField) const
Return the matrix norm used to normalise the residual for the.
const word & fieldName() const
Definition: LduMatrix.H:230
virtual void readControls()
Read the control parameters from the controlDict_.
static autoPtr< solver > New(const word &fieldName, const LduMatrix< Type, DType, LUType > &matrix, const dictionary &solverDict)
Return a new solver.
LduMatrix is a general matrix class in which the coefficients are stored as three arrays,...
Definition: LduMatrix.H:85
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: autoPtr.H:51
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:710
A class representing the concept of 1 (scalar(1)) used to avoid unnecessary manipulations for objects...
Definition: one.H:51
Traits class for primitives.
Definition: pTraits.H:53
A class for handling words, derived from string.
Definition: word.H:62
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:346
const volScalarField & psi
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static const zero Zero
Definition: zero.H:97
Type gSum(const FieldField< Field, Type > &f)
const doubleScalar e
Definition: doubleScalar.H:106
tmp< DimensionedField< scalar, GeoMesh > > stabilise(const DimensionedField< scalar, GeoMesh > &dsf, const dimensioned< scalar > &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
IOerror FatalIOError
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionSet cmptMag(const dimensionSet &)
Definition: dimensionSet.C:296
Type gAverage(const FieldField< Field, Type > &f)