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