fvMatrixSolve.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-2014 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 "diagTensorField.H"
28 
29 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
30 
31 template<class Type>
33 (
34  const label patchi,
35  const label facei,
36  const direction cmpt,
37  const scalar value
38 )
39 {
40  if (psi_.needReference())
41  {
42  if (Pstream::master())
43  {
44  internalCoeffs_[patchi][facei].component(cmpt) +=
45  diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]];
46 
47  boundaryCoeffs_[patchi][facei].component(cmpt) +=
48  diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]]
49  *value;
50  }
51  }
52 }
53 
54 
55 template<class Type>
57 (
58  const dictionary& solverControls
59 )
60 {
61  if (debug)
62  {
63  Info.masterStream(this->mesh().comm())
64  << "fvMatrix<Type>::solve(const dictionary& solverControls) : "
65  "solving fvMatrix<Type>"
66  << endl;
67  }
68 
69  label maxIter = -1;
70  if (solverControls.readIfPresent("maxIter", maxIter))
71  {
72  if (maxIter == 0)
73  {
74  return solverPerformance();
75  }
76  }
77 
78  word type(solverControls.lookupOrDefault<word>("type", "segregated"));
79 
80  if (type == "segregated")
81  {
82  return solveSegregated(solverControls);
83  }
84  else if (type == "coupled")
85  {
86  return solveCoupled(solverControls);
87  }
88  else
89  {
91  (
92  "fvMatrix<Type>::solve(const dictionary& solverControls)",
93  solverControls
94  ) << "Unknown type " << type
95  << "; currently supported solver types are segregated and coupled"
96  << exit(FatalIOError);
97 
98  return solverPerformance();
99  }
100 }
101 
102 
103 template<class Type>
105 (
106  const dictionary& solverControls
107 )
108 {
109  if (debug)
110  {
111  Info.masterStream(this->mesh().comm())
112  << "fvMatrix<Type>::solveSegregated"
113  "(const dictionary& solverControls) : "
114  "solving fvMatrix<Type>"
115  << endl;
116  }
117 
120 
121  solverPerformance solverPerfVec
122  (
123  "fvMatrix<Type>::solveSegregated",
124  psi.name()
125  );
126 
127  scalarField saveDiag(diag());
128 
129  Field<Type> source(source_);
130 
131  // At this point include the boundary source from the coupled boundaries.
132  // This is corrected for the implict part by updateMatrixInterfaces within
133  // the component loop.
134  addBoundarySource(source);
135 
136  typename Type::labelType validComponents
137  (
138  pow
139  (
140  psi.mesh().solutionD(),
142  )
143  );
144 
145  for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
146  {
147  if (validComponents[cmpt] == -1) continue;
148 
149  // copy field and source
150 
151  scalarField psiCmpt(psi.internalField().component(cmpt));
152  addBoundaryDiag(diag(), cmpt);
153 
154  scalarField sourceCmpt(source.component(cmpt));
155 
156  FieldField<Field, scalar> bouCoeffsCmpt
157  (
158  boundaryCoeffs_.component(cmpt)
159  );
160 
161  FieldField<Field, scalar> intCoeffsCmpt
162  (
163  internalCoeffs_.component(cmpt)
164  );
165 
166  lduInterfaceFieldPtrsList interfaces =
167  psi.boundaryField().scalarInterfaces();
168 
169  // Use the initMatrixInterfaces and updateMatrixInterfaces to correct
170  // bouCoeffsCmpt for the explicit part of the coupled boundary
171  // conditions
172  initMatrixInterfaces
173  (
174  bouCoeffsCmpt,
175  interfaces,
176  psiCmpt,
177  sourceCmpt,
178  cmpt
179  );
180 
181  updateMatrixInterfaces
182  (
183  bouCoeffsCmpt,
184  interfaces,
185  psiCmpt,
186  sourceCmpt,
187  cmpt
188  );
189 
190  solverPerformance solverPerf;
191 
192  // Solver call
193  solverPerf = lduMatrix::solver::New
194  (
195  psi.name() + pTraits<Type>::componentNames[cmpt],
196  *this,
197  bouCoeffsCmpt,
198  intCoeffsCmpt,
199  interfaces,
200  solverControls
201  )->solve(psiCmpt, sourceCmpt, cmpt);
202 
203  if (solverPerformance::debug)
204  {
205  solverPerf.print(Info.masterStream(this->mesh().comm()));
206  }
207 
208  solverPerfVec = max(solverPerfVec, solverPerf);
209  solverPerfVec.solverName() = solverPerf.solverName();
210 
211  psi.internalField().replace(cmpt, psiCmpt);
212  diag() = saveDiag;
213  }
214 
216 
217  psi.mesh().setSolverPerformance(psi.name(), solverPerfVec);
218 
219  return solverPerfVec;
220 }
221 
222 
223 template<class Type>
225 (
226  const dictionary& solverControls
227 )
228 {
229  if (debug)
230  {
231  Info.masterStream(this->mesh().comm())
232  << "fvMatrix<Type>::solveCoupled"
233  "(const dictionary& solverControls) : "
234  "solving fvMatrix<Type>"
235  << endl;
236  }
237 
240 
241  LduMatrix<Type, scalar, scalar> coupledMatrix(psi.mesh());
242  coupledMatrix.diag() = diag();
243  coupledMatrix.upper() = upper();
244  coupledMatrix.lower() = lower();
245  coupledMatrix.source() = source();
246 
247  addBoundaryDiag(coupledMatrix.diag(), 0);
248  addBoundarySource(coupledMatrix.source(), false);
249 
250  coupledMatrix.interfaces() = psi.boundaryField().interfaces();
251  coupledMatrix.interfacesUpper() = boundaryCoeffs().component(0);
252  coupledMatrix.interfacesLower() = internalCoeffs().component(0);
253 
255  coupledMatrixSolver
256  (
258  (
259  psi.name(),
260  coupledMatrix,
261  solverControls
262  )
263  );
264 
265  SolverPerformance<Type> solverPerf
266  (
267  coupledMatrixSolver->solve(psi)
268  );
269 
271  {
272  solverPerf.print(Info.masterStream(this->mesh().comm()));
273  }
274 
276 
277  // psi.mesh().setSolverPerformance(psi.name(), solverPerf);
278 
279  return solverPerformance();
280 }
281 
282 
283 template<class Type>
286 {
287  return solver
288  (
289  psi_.mesh().solverDict
290  (
291  psi_.select
292  (
293  psi_.mesh().data::template lookupOrDefault<bool>
294  ("finalIteration", false)
295  )
296  )
297  );
298 }
299 
300 
301 template<class Type>
303 {
304  return solve
305  (
306  fvMat_.psi_.mesh().solverDict
307  (
308  fvMat_.psi_.select
309  (
310  fvMat_.psi_.mesh().data::template lookupOrDefault<bool>
311  ("finalIteration", false)
312  )
313  )
314  );
315 }
316 
317 
318 template<class Type>
320 {
321  return solve
322  (
323  psi_.mesh().solverDict
324  (
325  psi_.select
326  (
327  psi_.mesh().data::template lookupOrDefault<bool>
328  ("finalIteration", false)
329  )
330  )
331  );
332 }
333 
334 
335 template<class Type>
337 {
338  tmp<Field<Type> > tres(new Field<Type>(source_));
339  Field<Type>& res = tres();
340 
341  addBoundarySource(res);
342 
343  // Loop over field components
344  for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
345  {
346  scalarField psiCmpt(psi_.internalField().component(cmpt));
347 
348  scalarField boundaryDiagCmpt(psi_.size(), 0.0);
349  addBoundaryDiag(boundaryDiagCmpt, cmpt);
350 
351  FieldField<Field, scalar> bouCoeffsCmpt
352  (
353  boundaryCoeffs_.component(cmpt)
354  );
355 
356  res.replace
357  (
358  cmpt,
360  (
361  psiCmpt,
362  res.component(cmpt) - boundaryDiagCmpt*psiCmpt,
363  bouCoeffsCmpt,
364  psi_.boundaryField().scalarInterfaces(),
365  cmpt
366  )
367  );
368  }
369 
370  return tres;
371 }
372 
373 
374 // ************************************************************************* //
unsigned char direction
Definition: direction.H:43
Generic field type.
Definition: FieldField.H:51
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
rhoEqn solve()
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
tmp< Field< cmptType > > component(const direction) const
Return a component field of the field.
Definition: Field.C:578
solverPerformance solve()
Solve returning the solution statistics.
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics...
Field< DType > & diag()
Definition: LduMatrix.C:190
A class for handling words, derived from string.
Definition: word.H:59
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
InternalField & internalField()
Return internal field.
solverPerformance solveCoupled(const dictionary &)
Solve coupled returning the solution statistics.
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
messageStream Info
const Mesh & mesh() const
Return mesh.
dynamicFvMesh & mesh
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
autoPtr< BasicCompressibleTurbulenceModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleTurbulenceModel::transportModel &transport, const word &propertiesName)
void setComponentReference(const label patchi, const label facei, const direction cmpt, const scalar value)
Set reference level for a component of the solution.
Definition: fvMatrixSolve.C:33
tmp< Field< Type > > residual() const
Return the matrix residual.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
IOerror FatalIOError
void residual(scalarField &rA, const scalarField &psi, const scalarField &source, const FieldField< Field, scalar > &interfaceBouCoeffs, const lduInterfaceFieldPtrsList &interfaces, const direction cmpt) const
void addBoundaryDiag(scalarField &diag, const direction cmpt) const
Definition: fvMatrix.C:117
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
A class representing the concept of 0 used to avoid unnecessary manipulations for objects that are kn...
Definition: zero.H:47
label patchi
const volScalarField & psi
Definition: createFields.H:24
solverPerformance solveSegregated(const dictionary &)
Solve segregated returning the solution statistics.
const word & name() const
Return name.
Definition: IOobject.H:260
fileName::Type type(const fileName &)
Return the file type: DIRECTORY or FILE.
Definition: POSIX.C:589
Generic GeometricField class.
Traits class for primitives.
Definition: pTraits.H:50
void addBoundarySource(Field< Type > &source, const bool couples=true) const
Definition: fvMatrix.C:151
const word & solverName() const
Return solver name.
void correctBoundaryConditions()
Correct boundary field.
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
Abstract base-class for lduMatrix solvers.
Definition: lduMatrix.H:91
void replace(const direction, const UList< cmptType > &)
Replace a component field of the field.
Definition: Field.C:590
OSstream & masterStream(const label communicator)
Convert to OSstream.
Definition: messageStream.C:60
solverPerformance solve()
Solve returning the solution statistics.
SolverPerformance< scalar > solverPerformance
SolverPerformance instantiated for a scalar.
LduMatrix is a general matrix class in which the coefficients are stored as three arrays...
Definition: LduMatrix.H:69
autoPtr< fvSolver > solver()
Construct and return the solver.
#define FatalIOErrorIn(functionName, ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:325
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:117
A class for managing temporary objects.
Definition: PtrList.H:118
void print(Ostream &os) const
Print summary of solver performance to the given stream.