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-2018 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<Type>();
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  solverControls
93  ) << "Unknown type " << type
94  << "; currently supported solver types are segregated and coupled"
95  << exit(FatalIOError);
96 
97  return SolverPerformance<Type>();
98  }
99 }
100 
101 
102 template<class Type>
104 (
105  const dictionary& solverControls
106 )
107 {
108  if (debug)
109  {
110  Info.masterStream(this->mesh().comm())
111  << "fvMatrix<Type>::solveSegregated"
112  "(const dictionary& solverControls) : "
113  "solving fvMatrix<Type>"
114  << endl;
115  }
116 
119 
120  SolverPerformance<Type> solverPerfVec
121  (
122  "fvMatrix<Type>::solveSegregated",
123  psi.name()
124  );
125 
126  scalarField saveDiag(diag());
127 
128  Field<Type> source(source_);
129 
130  // At this point include the boundary source from the coupled boundaries.
131  // This is corrected for the implict part by updateMatrixInterfaces within
132  // the component loop.
133  addBoundarySource(source);
134 
135  typename Type::labelType validComponents
136  (
137  psi.mesh().template validComponents<Type>()
138  );
139 
140  for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
141  {
142  if (validComponents[cmpt] == -1) continue;
143 
144  // copy field and source
145 
146  scalarField psiCmpt(psi.primitiveField().component(cmpt));
147  addBoundaryDiag(diag(), cmpt);
148 
149  scalarField sourceCmpt(source.component(cmpt));
150 
151  FieldField<Field, scalar> bouCoeffsCmpt
152  (
153  boundaryCoeffs_.component(cmpt)
154  );
155 
156  FieldField<Field, scalar> intCoeffsCmpt
157  (
158  internalCoeffs_.component(cmpt)
159  );
160 
161  lduInterfaceFieldPtrsList interfaces =
162  psi.boundaryField().scalarInterfaces();
163 
164  // Use the initMatrixInterfaces and updateMatrixInterfaces to correct
165  // bouCoeffsCmpt for the explicit part of the coupled boundary
166  // conditions
167  initMatrixInterfaces
168  (
169  bouCoeffsCmpt,
170  interfaces,
171  psiCmpt,
172  sourceCmpt,
173  cmpt
174  );
175 
176  updateMatrixInterfaces
177  (
178  bouCoeffsCmpt,
179  interfaces,
180  psiCmpt,
181  sourceCmpt,
182  cmpt
183  );
184 
185  solverPerformance solverPerf;
186 
187  // Solver call
188  solverPerf = lduMatrix::solver::New
189  (
190  psi.name() + pTraits<Type>::componentNames[cmpt],
191  *this,
192  bouCoeffsCmpt,
193  intCoeffsCmpt,
194  interfaces,
195  solverControls
196  )->solve(psiCmpt, sourceCmpt, cmpt);
197 
199  {
200  solverPerf.print(Info.masterStream(this->mesh().comm()));
201  }
202 
203  solverPerfVec.replace(cmpt, solverPerf);
204  solverPerfVec.solverName() = solverPerf.solverName();
205 
206  psi.primitiveFieldRef().replace(cmpt, psiCmpt);
207  diag() = saveDiag;
208  }
209 
211 
212  psi.mesh().setSolverPerformance(psi.name(), solverPerfVec);
213 
214  return solverPerfVec;
215 }
216 
217 
218 template<class Type>
220 (
221  const dictionary& solverControls
222 )
223 {
224  if (debug)
225  {
226  Info.masterStream(this->mesh().comm())
227  << "fvMatrix<Type>::solveCoupled"
228  "(const dictionary& solverControls) : "
229  "solving fvMatrix<Type>"
230  << endl;
231  }
232 
235 
236  LduMatrix<Type, scalar, scalar> coupledMatrix(psi.mesh());
237  coupledMatrix.diag() = diag();
238  coupledMatrix.upper() = upper();
239  coupledMatrix.lower() = lower();
240  coupledMatrix.source() = source();
241 
242  addBoundaryDiag(coupledMatrix.diag(), 0);
243  addBoundarySource(coupledMatrix.source(), false);
244 
245  coupledMatrix.interfaces() = psi.boundaryFieldRef().interfaces();
246  coupledMatrix.interfacesUpper() = boundaryCoeffs().component(0);
247  coupledMatrix.interfacesLower() = internalCoeffs().component(0);
248 
250  coupledMatrixSolver
251  (
253  (
254  psi.name(),
255  coupledMatrix,
256  solverControls
257  )
258  );
259 
260  SolverPerformance<Type> solverPerf
261  (
262  coupledMatrixSolver->solve(psi)
263  );
264 
266  {
267  solverPerf.print(Info.masterStream(this->mesh().comm()));
268  }
269 
271 
272  psi.mesh().setSolverPerformance(psi.name(), solverPerf);
273 
274  return solverPerf;
275 }
276 
277 
278 template<class Type>
281 {
282  return solver
283  (
284  psi_.mesh().solverDict
285  (
286  psi_.select
287  (
288  psi_.mesh().data::template lookupOrDefault<bool>
289  ("finalIteration", false)
290  )
291  )
292  );
293 }
294 
295 
296 template<class Type>
298 {
299  return solve
300  (
301  fvMat_.psi_.mesh().solverDict
302  (
303  fvMat_.psi_.select
304  (
305  fvMat_.psi_.mesh().data::template lookupOrDefault<bool>
306  ("finalIteration", false)
307  )
308  )
309  );
310 }
311 
312 
313 template<class Type>
315 {
316  return solve
317  (
318  psi_.mesh().solverDict
319  (
320  psi_.select
321  (
322  psi_.mesh().data::template lookupOrDefault<bool>
323  ("finalIteration", false)
324  )
325  )
326  );
327 }
328 
329 
330 template<class Type>
332 {
333  tmp<Field<Type>> tres(new Field<Type>(source_));
334  Field<Type>& res = tres.ref();
335 
336  addBoundarySource(res);
337 
338  // Loop over field components
339  for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
340  {
341  scalarField psiCmpt(psi_.primitiveField().component(cmpt));
342 
343  scalarField boundaryDiagCmpt(psi_.size(), 0.0);
344  addBoundaryDiag(boundaryDiagCmpt, cmpt);
345 
346  FieldField<Field, scalar> bouCoeffsCmpt
347  (
348  boundaryCoeffs_.component(cmpt)
349  );
350 
351  res.replace
352  (
353  cmpt,
355  (
356  psiCmpt,
357  res.component(cmpt) - boundaryDiagCmpt*psiCmpt,
358  bouCoeffsCmpt,
359  psi_.boundaryField().scalarInterfaces(),
360  cmpt
361  )
362  );
363  }
364 
365  return tres;
366 }
367 
368 
369 // ************************************************************************* //
type
Types of root.
Definition: Roots.H:52
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:297
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:137
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:174
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:256
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:669
Generic field type.
Definition: FieldField.H:51
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
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.
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
rhoEqn solve()
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
tmp< Field< cmptType > > component(const direction) const
Return a component field of the field.
Definition: Field.C:657
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
const Mesh & mesh() const
Return mesh.
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
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:331
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.
void correctBoundaryConditions()
Correct boundary field.
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:33
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.
IOerror FatalIOError