fvMatrixSolve.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-2019 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 #include "Residuals.H"
29 
30 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
31 
32 template<class Type>
34 (
35  const label patchi,
36  const label facei,
37  const direction cmpt,
38  const scalar value
39 )
40 {
41  if (psi_.needReference())
42  {
43  if (Pstream::master())
44  {
45  internalCoeffs_[patchi][facei].component(cmpt) +=
46  diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]];
47 
48  boundaryCoeffs_[patchi][facei].component(cmpt) +=
49  diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]]
50  *value;
51  }
52  }
53 }
54 
55 
56 template<class Type>
58 (
59  const dictionary& solverControls
60 )
61 {
62  if (debug)
63  {
64  Info.masterStream(this->mesh().comm())
65  << "fvMatrix<Type>::solve(const dictionary& solverControls) : "
66  "solving fvMatrix<Type>"
67  << endl;
68  }
69 
70  label maxIter = -1;
71  if (solverControls.readIfPresent("maxIter", maxIter))
72  {
73  if (maxIter == 0)
74  {
75  return SolverPerformance<Type>();
76  }
77  }
78 
79  word type(solverControls.lookupOrDefault<word>("type", "segregated"));
80 
81  if (type == "segregated")
82  {
83  return solveSegregated(solverControls);
84  }
85  else if (type == "coupled")
86  {
87  return solveCoupled(solverControls);
88  }
89  else
90  {
92  (
93  solverControls
94  ) << "Unknown type " << type
95  << "; currently supported solver types are segregated and coupled"
96  << exit(FatalIOError);
97 
98  return SolverPerformance<Type>();
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<Type> 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  psi.mesh().template validComponents<Type>()
139  );
140 
141  for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
142  {
143  if (validComponents[cmpt] == -1) continue;
144 
145  // copy field and source
146 
147  scalarField psiCmpt(psi.primitiveField().component(cmpt));
148  addBoundaryDiag(diag(), cmpt);
149 
150  scalarField sourceCmpt(source.component(cmpt));
151 
152  FieldField<Field, scalar> bouCoeffsCmpt
153  (
154  boundaryCoeffs_.component(cmpt)
155  );
156 
157  FieldField<Field, scalar> intCoeffsCmpt
158  (
159  internalCoeffs_.component(cmpt)
160  );
161 
162  lduInterfaceFieldPtrsList interfaces =
163  psi.boundaryField().scalarInterfaces();
164 
165  // Use the initMatrixInterfaces and updateMatrixInterfaces to correct
166  // bouCoeffsCmpt for the explicit part of the coupled boundary
167  // conditions
168  initMatrixInterfaces
169  (
170  bouCoeffsCmpt,
171  interfaces,
172  psiCmpt,
173  sourceCmpt,
174  cmpt
175  );
176 
177  updateMatrixInterfaces
178  (
179  bouCoeffsCmpt,
180  interfaces,
181  psiCmpt,
182  sourceCmpt,
183  cmpt
184  );
185 
186  solverPerformance solverPerf;
187 
188  // Solver call
189  solverPerf = lduMatrix::solver::New
190  (
191  psi.name() + pTraits<Type>::componentNames[cmpt],
192  *this,
193  bouCoeffsCmpt,
194  intCoeffsCmpt,
195  interfaces,
196  solverControls
197  )->solve(psiCmpt, sourceCmpt, cmpt);
198 
200  {
201  solverPerf.print(Info.masterStream(this->mesh().comm()));
202  }
203 
204  solverPerfVec.replace(cmpt, solverPerf);
205  solverPerfVec.solverName() = solverPerf.solverName();
206 
207  psi.primitiveFieldRef().replace(cmpt, psiCmpt);
208  diag() = saveDiag;
209  }
210 
212 
213  Residuals<Type>::append(psi.mesh(), solverPerfVec);
214 
215  return solverPerfVec;
216 }
217 
218 
219 template<class Type>
221 (
222  const dictionary& solverControls
223 )
224 {
225  if (debug)
226  {
227  Info.masterStream(this->mesh().comm())
228  << "fvMatrix<Type>::solveCoupled"
229  "(const dictionary& solverControls) : "
230  "solving fvMatrix<Type>"
231  << endl;
232  }
233 
236 
237  LduMatrix<Type, scalar, scalar> coupledMatrix(psi.mesh());
238  coupledMatrix.diag() = diag();
239  coupledMatrix.upper() = upper();
240  coupledMatrix.lower() = lower();
241  coupledMatrix.source() = source();
242 
243  addBoundaryDiag(coupledMatrix.diag(), 0);
244  addBoundarySource(coupledMatrix.source(), false);
245 
246  coupledMatrix.interfaces() = psi.boundaryFieldRef().interfaces();
247  coupledMatrix.interfacesUpper() = boundaryCoeffs().component(0);
248  coupledMatrix.interfacesLower() = internalCoeffs().component(0);
249 
251  coupledMatrixSolver
252  (
254  (
255  psi.name(),
256  coupledMatrix,
257  solverControls
258  )
259  );
260 
261  SolverPerformance<Type> solverPerf
262  (
263  coupledMatrixSolver->solve(psi)
264  );
265 
267  {
268  solverPerf.print(Info.masterStream(this->mesh().comm()));
269  }
270 
272 
273  Residuals<Type>::append(psi.mesh(), solverPerf);
274 
275  return solverPerf;
276 }
277 
278 
279 template<class Type>
282 {
283  return solver
284  (
285  psi_.mesh().solverDict
286  (
287  psi_.select
288  (
289  psi_.mesh().data::template lookupOrDefault<bool>
290  ("finalIteration", false)
291  )
292  )
293  );
294 }
295 
296 
297 template<class Type>
299 {
300  return solve
301  (
302  fvMat_.psi_.mesh().solverDict
303  (
304  fvMat_.psi_.select
305  (
306  fvMat_.psi_.mesh().data::template lookupOrDefault<bool>
307  ("finalIteration", false)
308  )
309  )
310  );
311 }
312 
313 
314 template<class Type>
316 {
317  return solve
318  (
319  psi_.mesh().solverDict
320  (
321  psi_.mesh().data::template lookupOrDefault<bool>
322  ("finalIteration", false)
323  ? word(name + "Final")
324  : name
325  )
326  );
327 }
328 
329 
330 template<class Type>
332 {
333  return solve
334  (
335  psi_.mesh().solverDict
336  (
337  psi_.select
338  (
339  psi_.mesh().data::template lookupOrDefault<bool>
340  (
341  "finalIteration",
342  false
343  )
344  )
345  )
346  );
347 }
348 
349 
350 template<class Type>
352 {
353  tmp<Field<Type>> tres(new Field<Type>(source_));
354  Field<Type>& res = tres.ref();
355 
356  addBoundarySource(res);
357 
358  // Loop over field components
359  for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
360  {
361  scalarField psiCmpt(psi_.primitiveField().component(cmpt));
362 
363  scalarField boundaryDiagCmpt(psi_.size(), 0.0);
364  addBoundaryDiag(boundaryDiagCmpt, cmpt);
365 
366  FieldField<Field, scalar> bouCoeffsCmpt
367  (
368  boundaryCoeffs_.component(cmpt)
369  );
370 
371  res.replace
372  (
373  cmpt,
375  (
376  psiCmpt,
377  res.component(cmpt) - boundaryDiagCmpt*psiCmpt,
378  bouCoeffsCmpt,
379  psi_.boundaryField().scalarInterfaces(),
380  cmpt
381  )
382  );
383  }
384 
385  return tres;
386 }
387 
388 
389 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
390 
391 template<class Type>
393 (
394  fvMatrix<Type>& fvm,
395  const word& name
396 )
397 {
398  return fvm.solve(name);
399 }
400 
401 
402 template<class Type>
404 (
405  const tmp<fvMatrix<Type>>& tfvm,
406  const word& name
407 )
408 {
409  SolverPerformance<Type> solverPerf =
410  const_cast<fvMatrix<Type>&>(tfvm()).solve(name);
411 
412  tfvm.clear();
413 
414  return solverPerf;
415 }
416 
417 
418 template<class Type>
420 {
421  return fvm.solve();
422 }
423 
424 template<class Type>
426 {
427  SolverPerformance<Type> solverPerf =
428  const_cast<fvMatrix<Type>&>(tfvm()).solve();
429 
430  tfvm.clear();
431 
432  return solverPerf;
433 }
434 
435 
436 // ************************************************************************* //
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
const word & name() const
Return name.
Definition: IOobject.H:303
void addBoundaryDiag(scalarField &diag, const direction cmpt) const
Definition: fvMatrix.C:111
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:156
const Boundary & boundaryField() const
Return const-reference to the boundary field.
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
uint8_t direction
Definition: direction.H:45
const word & solverName() const
Return solver name.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Traits class for primitives.
Definition: pTraits.H:50
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Generic GeometricField class.
void replace(const direction, const UList< cmptType > &)
Replace a component field of the field.
Definition: Field.C:453
Generic field type.
Definition: FieldField.H:51
SolverPerformance< Type > solveCoupled(const dictionary &)
Solve coupled returning the solution statistics.
void residual(scalarField &rA, const scalarField &psi, const scalarField &source, const FieldField< Field, scalar > &interfaceBouCoeffs, const lduInterfaceFieldPtrsList &interfaces, const direction cmpt) const
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
dynamicFvMesh & mesh
SolverPerformance< Type > solve()
Solve returning the solution statistics.
autoPtr< BasicCompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleMomentumTransportModel::transportModel &transport)
A class for handling words, derived from string.
Definition: word.H:59
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics...
Abstract base-class for lduMatrix solvers.
Definition: lduMatrix.H:93
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:72
tmp< Field< cmptType > > component(const direction) const
Return a component field of the field.
Definition: Field.C:441
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
SolverPerformance< Type > solve(const dictionary &)
Solve segregated or coupled returning the solution statistics.
Definition: fvMatrixSolve.C:58
const Mesh & mesh() const
Return mesh.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
void print(Ostream &os) const
Print summary of solver performance to the given stream.
OSstream & masterStream(const label communicator)
Convert to OSstream.
Definition: messageStream.C:59
Field< DType > & diag()
Definition: LduMatrix.C:190
label patchi
rhoEqn solve()
fvModels source(alpha1, mixture.thermo1().rho())
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:335
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.
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
void correctBoundaryConditions()
Correct boundary field.
SolverPerformance< Type > solve(fvMatrix< Type > &, const word &)
Solve returning the solution statistics given convergence tolerance.
messageStream Info
void addBoundarySource(Field< Type > &source, const bool couples=true) const
Definition: fvMatrix.C:145
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
const volScalarField & psi
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:34
SolverPerformance< Type > solveSegregated(const dictionary &)
Solve segregated returning the solution statistics.
A class for managing temporary objects.
Definition: PtrList.H:53
tmp< Field< Type > > residual() const
Return the matrix residual.
SolverPerformance< Type > solve()
Solve returning the solution statistics.
MeshObject to store the solver performance residuals of all the fields of the type it is instantiated...
Definition: Residuals.H:52
IOerror FatalIOError