fvScalarMatrix.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2016-2021 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "fvScalarMatrix.H"
31 #include "profiling.H"
32 #include "PrecisionAdaptor.H"
33 #include "jumpCyclicFvPatchField.H"
34 #include "cyclicPolyPatch.H"
35 #include "cyclicAMIPolyPatch.H"
36 
37 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
38 
39 template<>
41 (
42  const label patchi,
43  const label facei,
44  const direction,
45  const scalar value
46 )
47 {
48  if (psi_.needReference())
49  {
50  if (Pstream::master())
51  {
52  internalCoeffs_[patchi][facei] +=
53  diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]];
54 
55  boundaryCoeffs_[patchi][facei] +=
56  diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]]
57  *value;
58  }
59  }
60 }
61 
62 
63 template<>
66 (
67  const dictionary& solverControls
68 )
69 {
71  if (psi_.mesh().name() != polyMesh::defaultRegion)
72  {
73  regionName = psi_.mesh().name() + "::";
74  }
75  addProfiling(solve, "fvMatrix::solve." + regionName + psi_.name());
76 
77  if (debug)
78  {
79  Info.masterStream(this->mesh().comm())
80  << "fvMatrix<scalar>::solver(const dictionary& solverControls) : "
81  "solver for fvMatrix<scalar>"
82  << endl;
83  }
84 
85  scalarField saveDiag(diag());
86  addBoundaryDiag(diag(), 0);
87 
88  lduInterfaceFieldPtrsList interfaces =
89  psi_.boundaryField().scalarInterfaces();
90 
91 
93  (
95  (
96  *this,
98  (
99  psi_.name(),
100  *this,
101  boundaryCoeffs_,
102  internalCoeffs_,
103  interfaces,
104  solverControls
105  )
106  )
107  );
108 
109  diag() = saveDiag;
110 
111  return solverPtr;
112 }
113 
114 
115 template<>
117 (
118  const dictionary& solverControls
119 )
120 {
121  const int logLevel =
122  solverControls.getOrDefault<int>
123  (
124  "log",
126  );
127 
128  auto& psi =
130  (
131  fvMat_.psi()
132  );
133 
134  scalarField saveDiag(fvMat_.diag());
135  fvMat_.addBoundaryDiag(fvMat_.diag(), 0);
136 
137  scalarField totalSource(fvMat_.source());
138  fvMat_.addBoundarySource(totalSource, false);
139 
140  // Assign new solver controls
141  solver_->read(solverControls);
142 
143  solverPerformance solverPerf = solver_->solve
144  (
145  psi.primitiveFieldRef(),
146  totalSource
147  );
148 
149  if (logLevel)
150  {
151  solverPerf.print(Info.masterStream(fvMat_.mesh().comm()));
152  }
153 
154  fvMat_.diag() = saveDiag;
155 
156  psi.correctBoundaryConditions();
157 
158  psi.mesh().setSolverPerformance(psi.name(), solverPerf);
159 
160  return solverPerf;
161 }
162 
163 
164 template<>
166 (
167  const dictionary& solverControls
168 )
169 {
170  if (debug)
171  {
172  Info.masterStream(this->mesh().comm())
173  << "fvMatrix<scalar>::solveSegregated"
174  "(const dictionary& solverControls) : "
175  "solving fvMatrix<scalar>"
176  << endl;
177  }
178 
179  const int logLevel =
180  solverControls.getOrDefault<int>
181  (
182  "log",
184  );
185 
186  scalarField saveLower;
187  scalarField saveUpper;
188 
189  if (useImplicit_)
190  {
191  createOrUpdateLduPrimitiveAssembly();
192 
193  if (psi_.mesh().fluxRequired(psi_.name()))
194  {
195  // Save lower/upper for flux calculation
196  if (asymmetric())
197  {
198  saveLower = lower();
199  }
200  saveUpper = upper();
201  }
202 
203  setLduMesh(*lduMeshPtr());
204  transferFvMatrixCoeffs();
205  setBounAndInterCoeffs();
206  direction cmpt = 0;
207  manipulateMatrix(cmpt);
208  }
209 
210  scalarField saveDiag(diag());
211  addBoundaryDiag(diag(), 0);
212 
213  scalarField totalSource(source_);
214  addBoundarySource(totalSource, false);
215 
216  lduInterfaceFieldPtrsList interfaces;
217  PtrDynList<lduInterfaceField> newInterfaces;
218  if (!useImplicit_)
219  {
220  interfaces = this->psi(0).boundaryField().scalarInterfaces();
221  }
222  else
223  {
224  setInterfaces(interfaces, newInterfaces);
225  }
226 
227  tmp<scalarField> tpsi;
228  if (!useImplicit_)
229  {
230  tpsi.ref
231  (
233  (
234  psi_
236  );
237  }
238  else
239  {
240  tpsi = tmp<scalarField>::New(lduAddr().size(), Zero);
241  scalarField& psi = tpsi.ref();
242 
243  for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
244  {
245  const label cellOffset = lduMeshPtr()->cellOffsets()[fieldi];
246  const auto& psiInternal = this->psi(fieldi).primitiveField();
247 
248  forAll(psiInternal, localCellI)
249  {
250  psi[cellOffset + localCellI] = psiInternal[localCellI];
251  }
252  }
253  }
254  scalarField& psi = tpsi.ref();
255 
256  // Solver call
258  (
259  this->psi(0).name(),
260  *this,
261  boundaryCoeffs_,
262  internalCoeffs_,
263  interfaces,
264  solverControls
265  )->solve(psi, totalSource);
266 
267  if (useImplicit_)
268  {
269  for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
270  {
271  auto& psiInternal =
273  (
274  this->psi(fieldi)
275  ).primitiveFieldRef();
276 
277  const label cellOffset = lduMeshPtr()->cellOffsets()[fieldi];
278 
279  forAll(psiInternal, localCellI)
280  {
281  psiInternal[localCellI] = psi[localCellI + cellOffset];
282  }
283  }
284  }
285 
286  if (logLevel)
287  {
288  solverPerf.print(Info.masterStream(mesh().comm()));
289  }
290 
291  diag() = saveDiag;
292 
293  if (useImplicit_)
294  {
295  if (psi_.mesh().fluxRequired(psi_.name()))
296  {
297  // Restore lower/upper
298  if (asymmetric())
299  {
300  lower().setSize(saveLower.size());
301  lower() = saveLower;
302  }
303 
304  upper().setSize(saveUpper.size());
305  upper() = saveUpper;
306  }
307  // Set the original lduMesh
308  setLduMesh(psi_.mesh());
309  }
310 
311  for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
312  {
313  auto& localPsi =
315  (
316  this->psi(fieldi)
317  );
318 
319  localPsi.correctBoundaryConditions();
320  localPsi.mesh().setSolverPerformance(localPsi.name(), solverPerf);
321  }
322 
323  return solverPerf;
324 }
325 
326 
327 template<>
329 {
330  scalarField boundaryDiag(psi_.size(), Zero);
331  addBoundaryDiag(boundaryDiag, 0);
332 
333  const scalarField& psif = psi_.primitiveField();
335  const solveScalarField& psi = tpsi();
336 
338  (
339  lduMatrix::residual
340  (
341  psi,
342  source_ - boundaryDiag*psif,
343  boundaryCoeffs_,
344  psi_.boundaryField().scalarInterfaces(),
345  0
346  )
347  );
348 
350  addBoundarySource(tres_s.ref());
351 
352  return tres_s;
353 }
354 
355 
356 template<>
358 {
359  tmp<volScalarField> tHphi
360  (
361  new volScalarField
362  (
363  IOobject
364  (
365  "H("+psi_.name()+')',
366  psi_.instance(),
367  psi_.mesh(),
368  IOobject::NO_READ,
369  IOobject::NO_WRITE
370  ),
371  psi_.mesh(),
372  dimensions_/dimVol,
373  extrapolatedCalculatedFvPatchScalarField::typeName
374  )
375  );
376  volScalarField& Hphi = tHphi.ref();
377 
378  Hphi.primitiveFieldRef() = (lduMatrix::H(psi_.primitiveField()) + source_);
379  addBoundarySource(Hphi.primitiveFieldRef());
380 
381  Hphi.primitiveFieldRef() /= psi_.mesh().V();
383 
384  return tHphi;
385 }
386 
387 
388 template<>
390 {
392  (
393  new volScalarField
394  (
395  IOobject
396  (
397  "H(1)",
398  psi_.instance(),
399  psi_.mesh(),
400  IOobject::NO_READ,
401  IOobject::NO_WRITE
402  ),
403  psi_.mesh(),
404  dimensions_/(dimVol*psi_.dimensions()),
405  extrapolatedCalculatedFvPatchScalarField::typeName
406  )
407  );
408  volScalarField& H1_ = tH1.ref();
409 
410  H1_.primitiveFieldRef() = lduMatrix::H1();
411  //addBoundarySource(Hphi.primitiveField());
412 
413  H1_.primitiveFieldRef() /= psi_.mesh().V();
414  H1_.correctBoundaryConditions();
415 
416  return tH1;
417 }
418 
419 
420 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::fvMatrix::residual
tmp< Field< Type > > residual() const
Return the matrix residual.
Definition: fvMatrixSolve.C:350
profiling.H
Foam::fvMatrix::H
tmp< GeometricField< Type, fvPatchField, volMesh > > H() const
Return the H operation source.
Definition: fvMatrix.C:1423
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:169
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
cyclicPolyPatch.H
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::diag
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
Definition: pointPatchFieldFunctions.H:287
Foam::messageStream::masterStream
OSstream & masterStream(const label communicator)
Definition: messageStream.C:140
Foam::fvMatrix::fvSolver::solve
SolverPerformance< Type > solve()
Solve returning the solution statistics.
Definition: fvMatrixSolve.C:336
H
volScalarField H(IOobject("H", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), mesh, dimensionedScalar(dimLength, Zero))
PrecisionAdaptor.H
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::PtrDynList
A dynamically resizable PtrList with allocation management.
Definition: PtrDynList.H:54
primitiveFieldRef
conserve primitiveFieldRef()+
solve
CEqn solve()
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::SolverPerformance::print
void print(Ostream &os) const
Print summary of solver performance to the given stream.
Definition: SolverPerformance.C:93
regionName
Foam::word regionName
Definition: createNamedDynamicFvMesh.H:1
Foam::tmp::ref
T & ref() const
Definition: tmpI.H:227
Foam::stringOps::lower
string lower(const std::string &s)
Return string copy transformed with std::tolower on each character.
Definition: stringOps.C:1184
Foam::Field< scalar >
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::fvMatrix::H1
tmp< volScalarField > H1() const
Return H(1)
Definition: fvMatrix.C:1485
Foam::UPtrList< const lduInterfaceField >
addProfiling
#define addProfiling(name, descr)
Define profiling trigger with specified name and description string.
Definition: profilingTrigger.H:113
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
fvScalarMatrix.H
A scalar instance of fvMatrix.
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
cyclicAMIPolyPatch.H
Foam::fvMatrix::solver
autoPtr< fvSolver > solver()
Construct and return the solver.
Definition: fvMatrixSolve.C:329
Foam::GeometricField::primitiveFieldRef
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field.
Definition: GeometricField.C:766
Foam::GeometricField::correctBoundaryConditions
void correctBoundaryConditions()
Correct boundary field.
Definition: GeometricField.C:940
Foam::New
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
Definition: DimensionedFieldReuseFunctions.H:105
Foam::autoPtr
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: HashPtrTable.H:53
Foam::GeometricField::ref
Internal & ref(const bool updateAccessTime=true)
Return a reference to the dimensioned internal field.
Definition: GeometricField.C:749
Foam::ConstPrecisionAdaptor
A const Field/List wrapper with possible data conversion.
Definition: PrecisionAdaptor.H:57
Foam::direction
uint8_t direction
Definition: direction.H:52
Foam::fvMatrix::fvSolver
Definition: fvMatrix.H:258
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
Foam::dimVol
const dimensionSet dimVol(dimVolume)
Older spelling for dimVolume.
Definition: dimensionSets.H:64
Foam::stringOps::upper
string upper(const std::string &s)
Return string copy transformed with std::toupper on each character.
Definition: stringOps.C:1200
Foam::refPtr< Container< Type > >::ref
Container< Type > & ref() const
Definition: refPtrI.H:203
Foam::fvMatrix::solveSegregated
SolverPerformance< Type > solveSegregated(const dictionary &)
Solve segregated returning the solution statistics.
Definition: fvMatrixSolve.C:112
Foam::dictionary::getOrDefault
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:148
Foam::GeometricField< scalar, fvPatchField, volMesh >
psi
const volScalarField & psi
Definition: createFieldRefs.H:1
jumpCyclicFvPatchField.H
Foam::SolverPerformance
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics.
Definition: SolverPerformance.H:52
Foam::fvMatrix::setComponentReference
void setComponentReference(const label patchi, const label facei, const direction cmpt, const scalar value)
Definition: fvMatrixSolve.C:38
extrapolatedCalculatedFvPatchFields.H