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-2022 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(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(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 implicit 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 =
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(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(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(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().solution().solverDict
286  (
287  psi_.select
288  (
289  !psi_.mesh().schemes().steady()
290  && psi_.mesh().data::template lookupOrDefault<bool>
291  ("finalIteration", false)
292  )
293  )
294  );
295 }
296 
297 
298 template<class Type>
300 {
301  return solve
302  (
303  fvMat_.psi_.mesh().solution().solverDict
304  (
305  fvMat_.psi_.select
306  (
307  !fvMat_.psi_.mesh().schemes().steady()
308  && fvMat_.psi_.mesh().data::template lookupOrDefault<bool>
309  ("finalIteration", false)
310  )
311  )
312  );
313 }
314 
315 
316 template<class Type>
318 {
319  return solve
320  (
321  psi_.mesh().solution().solverDict
322  (
323  !psi_.mesh().schemes().steady()
324  && psi_.mesh().data::template lookupOrDefault<bool>
325  ("finalIteration", false)
326  ? word(name + "Final")
327  : name
328  )
329  );
330 }
331 
332 
333 template<class Type>
335 {
336  return solve
337  (
338  psi_.mesh().solution().solverDict
339  (
340  psi_.select
341  (
342  !psi_.mesh().schemes().steady()
343  && psi_.mesh().data::template lookupOrDefault<bool>
344  (
345  "finalIteration",
346  false
347  )
348  )
349  )
350  );
351 }
352 
353 
354 template<class Type>
356 {
357  tmp<Field<Type>> tres(new Field<Type>(source_));
358  Field<Type>& res = tres.ref();
359 
360  addBoundarySource(res);
361 
362  // Loop over field components
363  for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
364  {
365  scalarField psiCmpt(psi_.primitiveField().component(cmpt));
366 
367  scalarField boundaryDiagCmpt(psi_.size(), 0.0);
368  addBoundaryDiag(boundaryDiagCmpt, cmpt);
369 
370  FieldField<Field, scalar> bouCoeffsCmpt
371  (
372  boundaryCoeffs_.component(cmpt)
373  );
374 
375  res.replace
376  (
377  cmpt,
379  (
380  psiCmpt,
381  res.component(cmpt) - boundaryDiagCmpt*psiCmpt,
382  bouCoeffsCmpt,
383  psi_.boundaryField().scalarInterfaces(),
384  cmpt
385  )
386  );
387  }
388 
389  return tres;
390 }
391 
392 
393 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
394 
395 template<class Type>
397 (
398  fvMatrix<Type>& fvm,
399  const word& name
400 )
401 {
402  return fvm.solve(name);
403 }
404 
405 
406 template<class Type>
408 (
409  const tmp<fvMatrix<Type>>& tfvm,
410  const word& name
411 )
412 {
413  SolverPerformance<Type> solverPerf =
414  const_cast<fvMatrix<Type>&>(tfvm()).solve(name);
415 
416  tfvm.clear();
417 
418  return solverPerf;
419 }
420 
421 
422 template<class Type>
424 {
425  return fvm.solve();
426 }
427 
428 template<class Type>
430 {
431  SolverPerformance<Type> solverPerf =
432  const_cast<fvMatrix<Type>&>(tfvm()).solve();
433 
434  tfvm.clear();
435 
436  return solverPerf;
437 }
438 
439 
440 // ************************************************************************* //
const word & name() const
Return name.
Definition: IOobject.H:315
autoPtr< CompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const viscosity &viscosity)
void addBoundaryDiag(scalarField &diag, const direction cmpt) const
Definition: fvMatrix.C:113
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.
lduInterfaceFieldPtrsList scalarInterfaces() const
Return a list of pointers for each patch field with only those.
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:467
fvMesh & mesh
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 > &)
SolverPerformance< Type > solve()
Solve returning the solution statistics.
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
const volScalarField & psi
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:455
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.
Field< DType > & diag()
Definition: LduMatrix.C:190
label patchi
rhoEqn solve()
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:318
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.
LduInterfaceFieldPtrsList< Type > interfaces() const
Return a list of pointers for each patch field with only those.
messageStream Info
void addBoundarySource(Field< Type > &source, const bool couples=true) const
Definition: fvMatrix.C:147
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
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