LduMatrixSolver.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 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 }
114 
115 
116 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
117 
118 template<class Type, class DType, class LUType>
120 (
121  const word& fieldName,
122  const LduMatrix<Type, DType, LUType>& matrix,
123  const dictionary& solverDict
124 )
125 :
126  fieldName_(fieldName),
127  matrix_(matrix),
128 
129  controlDict_(solverDict),
130 
131  maxIter_(1000),
132  minIter_(0),
133  tolerance_(1e-6*pTraits<Type>::one),
134  relTol_(Zero)
135 {
136  readControls();
137 }
138 
139 
140 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
141 
142 template<class Type, class DType, class LUType>
144 {
145  readControl(controlDict_, maxIter_, "maxIter");
146  readControl(controlDict_, minIter_, "minIter");
147  readControl(controlDict_, tolerance_, "tolerance");
148  readControl(controlDict_, relTol_, "relTol");
149 }
150 
151 
152 template<class Type, class DType, class LUType>
154 (
155  const dictionary& solverDict
156 )
157 {
158  controlDict_ = solverDict;
159  readControls();
160 }
161 
162 
163 template<class Type, class DType, class LUType>
165 (
166  const Field<Type>& psi,
167  const Field<Type>& Apsi,
168  Field<Type>& tmpField
169 ) const
170 {
171  // --- Calculate A dot reference value of psi
172  matrix_.sumA(tmpField);
173  cmptMultiply(tmpField, tmpField, gAverage(psi));
174 
175  return stabilise
176  (
177  gSum(cmptMag(Apsi - tmpField) + cmptMag(matrix_.source() - tmpField)),
179  );
180 
181  // At convergence this simpler method is equivalent to the above
182  // return stabilise(2*gSumCmptMag(matrix_.source()), matrix_.small_);
183 }
184 
185 
186 // ************************************************************************* //
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
Foam::DiagonalSolver.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
Traits class for primitives.
Definition: pTraits.H:50
bool asymmetric() const
Definition: LduMatrix.H:583
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.
label minIter_
Minimum number of iterations in the solver.
Definition: LduMatrix.H:128
dictionary controlDict_
Dictionary of controls.
Definition: LduMatrix.H:122
const dimensionedScalar e
Elementary charge.
Definition: doubleFloat.H:78
Type gSum(const FieldField< Field, Type > &f)
Pre-declare SubField and related Field type.
Definition: Field.H:57
A class for handling words, derived from string.
Definition: word.H:59
void readControl(const dictionary &controlDict, T &control, const word &controlName)
Read a control parameter from controlDict.
Definition: LduMatrixI.H:31
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics...
static const zero Zero
Definition: zero.H:91
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
solver(const word &fieldName, const LduMatrix< Type, DType, LUType > &matrix, const dictionary &solverDict)
tmp< DimensionedField< scalar, GeoMesh > > stabilise(const DimensionedField< scalar, GeoMesh > &dsf, const dimensioned< scalar > &ds)
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:331
Type relTol_
Convergence tolerance relative to the initial.
Definition: LduMatrix.H:134
Type gAverage(const FieldField< Field, Type > &f)
LduMatrix is a general matrix class in which the coefficients are stored as three arrays...
Definition: LduMatrix.H:69
virtual void read(const dictionary &solverDict)
Read and reset the solver parameters from the given dictionary.
label maxIter_
Maximum number of iterations in the solver.
Definition: LduMatrix.H:125
bool symmetric() const
Definition: LduMatrix.H:578
static autoPtr< solver > New(const word &fieldName, const LduMatrix< Type, DType, LUType > &matrix, const dictionary &solverDict)
Return a new solver.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:53
const LduMatrix< Type, DType, LUType > & matrix_
Definition: LduMatrix.H:119
virtual void readControls()
Read the control parameters from the controlDict_.
Type tolerance_
Final convergence tolerance.
Definition: LduMatrix.H:131
bool diagonal() const
Definition: LduMatrix.H:573
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:451
IOerror FatalIOError