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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2019-2021 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "LduMatrix.H"
30 #include "DiagonalSolver.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 template<class Type, class DType, class LUType>
37 (
38  const word& fieldName,
39  const LduMatrix<Type, DType, LUType>& matrix,
40  const dictionary& solverDict
41 )
42 {
43  const word solverName(solverDict.get<word>("solver"));
44 
45  if (matrix.diagonal())
46  {
48  (
50  (
51  fieldName,
52  matrix,
53  solverDict
54  )
55  );
56  }
57  else if (matrix.symmetric())
58  {
59  auto* ctorPtr = symMatrixConstructorTable(solverName);
60 
61  if (!ctorPtr)
62  {
64  (
65  solverDict,
66  "symmetric matrix solver",
67  solverName,
68  *symMatrixConstructorTablePtr_
69  ) << exit(FatalIOError);
70  }
71 
73  (
74  ctorPtr
75  (
76  fieldName,
77  matrix,
78  solverDict
79  )
80  );
81  }
82  else if (matrix.asymmetric())
83  {
84  auto* ctorPtr = asymMatrixConstructorTable(solverName);
85 
86  if (!ctorPtr)
87  {
89  (
90  solverDict,
91  "asymmetric matrix solver",
92  solverName,
93  *asymMatrixConstructorTablePtr_
94  ) << exit(FatalIOError);
95  }
96 
98  (
99  ctorPtr
100  (
101  fieldName,
102  matrix,
103  solverDict
104  )
105  );
106  }
107 
108  FatalIOErrorInFunction(solverDict)
109  << "cannot solve incomplete matrix, "
110  "no diagonal or off-diagonal coefficient"
111  << exit(FatalIOError);
112 
113  return nullptr;
114 }
115 
116 
117 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
118 
119 template<class Type, class DType, class LUType>
121 (
122  const word& fieldName,
123  const LduMatrix<Type, DType, LUType>& matrix,
124  const dictionary& solverDict
125 )
126 :
127  fieldName_(fieldName),
128  matrix_(matrix),
129 
130  controlDict_(solverDict),
131 
132  log_(1),
133  minIter_(0),
134  maxIter_(defaultMaxIter_),
135  tolerance_(1e-6*pTraits<Type>::one),
136  relTol_(Zero)
137 {
138  readControls();
139 }
140 
141 
142 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
143 
144 template<class Type, class DType, class LUType>
146 {
148  controlDict_.readIfPresent("minIter", minIter_);
149  controlDict_.readIfPresent("maxIter", maxIter_);
150  controlDict_.readIfPresent("tolerance", tolerance_);
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::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
Foam::DiagonalSolver
Foam::DiagonalSolver.
Definition: DiagonalSolver.H:51
Foam::cmptMultiply
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::gAverage
Type gAverage(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:604
Foam::LduMatrix::solver::log_
int log_
Level of verbosity in the solver output statements.
Definition: LduMatrix.H:131
Foam::LduMatrix
LduMatrix is a general matrix class in which the coefficients are stored as three arrays,...
Definition: LduMatrix.H:72
Foam::LduMatrix::solver::maxIter_
label maxIter_
Maximum number of iterations in the solver.
Definition: LduMatrix.H:137
Foam::FatalIOError
IOerror FatalIOError
Foam::gSum
Type gSum(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:594
Foam::dictionary::get
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:107
LduMatrix.H
Foam::LduMatrix::symmetric
bool symmetric() const
Definition: LduMatrix.H:581
Foam::cmptMag
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:400
DiagonalSolver.H
Foam::LduMatrix::solver::controlDict_
dictionary controlDict_
Dictionary of controls.
Definition: LduMatrix.H:128
FatalIOErrorInLookup
#define FatalIOErrorInLookup(ios, lookupTag, lookupName, lookupTable)
Report an error message using Foam::FatalIOError.
Definition: error.H:478
Foam::Field
Generic templated field type.
Definition: Field.H:63
Foam::LduMatrix::diagonal
bool diagonal() const
Definition: LduMatrix.H:576
Foam::LduMatrix::solver::read
virtual void read(const dictionary &solverDict)
Read and reset the solver parameters from the given dictionary.
Definition: LduMatrixSolver.C:157
Foam::stabilise
tmp< DimensionedField< scalar, GeoMesh > > stabilise(const DimensionedField< scalar, GeoMesh > &dsf, const dimensioned< scalar > &ds)
Definition: DimensionedScalarField.C:43
Foam::LduMatrix::solver::relTol_
Type relTol_
Convergence tolerance relative to the initial.
Definition: LduMatrix.H:143
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
Foam::LduMatrix::solver::tolerance_
Type tolerance_
Final convergence tolerance.
Definition: LduMatrix.H:140
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::autoPtr
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: HashPtrTable.H:53
Foam::LduMatrix::solver::minIter_
label minIter_
Minimum number of iterations in the solver.
Definition: LduMatrix.H:134
Foam::pTraits
A traits class, which is primarily used for primitives.
Definition: pTraits.H:56
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::LduMatrix::solver::solver
solver(const word &fieldName, const LduMatrix< Type, DType, LUType > &matrix, const dictionary &solverDict)
Definition: LduMatrixSolver.C:121
FatalIOErrorInFunction
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:473
psi
const volScalarField & psi
Definition: createFieldRefs.H:1
Foam::SolverPerformance
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics.
Definition: SolverPerformance.H:52
Foam::LduMatrix::asymmetric
bool asymmetric() const
Definition: LduMatrix.H:586
Foam::LduMatrix::solver::normFactor
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.
Definition: LduMatrixSolver.C:168
Foam::dictionary::readIfPresent
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:405
Foam::LduMatrix::solver::New
static autoPtr< solver > New(const word &fieldName, const LduMatrix< Type, DType, LUType > &matrix, const dictionary &solverDict)
Return a new solver.
Definition: LduMatrixSolver.C:37
Foam::LduMatrix::solver::readControls
virtual void readControls()
Read the control parameters from the controlDict_.
Definition: LduMatrixSolver.C:145