All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MRFZone.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 "MRFZone.H"
27 #include "fvMesh.H"
28 #include "volFields.H"
29 #include "surfaceFields.H"
30 #include "fvMatrices.H"
31 #include "faceSet.H"
32 #include "geometricOneField.H"
33 #include "syncTools.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39  defineTypeNameAndDebug(MRFZone, 0);
40 }
41 
42 
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 
45 void Foam::MRFZone::setMRFFaces()
46 {
47  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
48 
49  // Type per face:
50  // 0:not in zone
51  // 1:moving with frame
52  // 2:other
53  labelList faceType(mesh_.nFaces(), 0);
54 
55  // Determine faces in cell zone
56  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
57  // (without constructing cells)
58 
59  const labelList& own = mesh_.faceOwner();
60  const labelList& nei = mesh_.faceNeighbour();
61 
62  // Cells in zone
63  boolList zoneCell(mesh_.nCells(), false);
64 
65  if (cellZoneID_ != -1)
66  {
67  const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
68  forAll(cellLabels, i)
69  {
70  zoneCell[cellLabels[i]] = true;
71  }
72  }
73 
74 
75  label nZoneFaces = 0;
76 
77  for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
78  {
79  if (zoneCell[own[facei]] || zoneCell[nei[facei]])
80  {
81  faceType[facei] = 1;
82  nZoneFaces++;
83  }
84  }
85 
86 
87  labelHashSet excludedPatches(excludedPatchLabels_);
88 
89  forAll(patches, patchi)
90  {
91  const polyPatch& pp = patches[patchi];
92 
93  if (pp.coupled() || excludedPatches.found(patchi))
94  {
95  forAll(pp, i)
96  {
97  const label facei = pp.start() + i;
98 
99  if (zoneCell[own[facei]])
100  {
101  faceType[facei] = 2;
102  nZoneFaces++;
103  }
104  }
105  }
106  else if (!isA<emptyPolyPatch>(pp))
107  {
108  forAll(pp, i)
109  {
110  const label facei = pp.start() + i;
111 
112  if (zoneCell[own[facei]])
113  {
114  faceType[facei] = 1;
115  nZoneFaces++;
116  }
117  }
118  }
119  }
120 
121  // Synchronise the faceType across processor patches
122  syncTools::syncFaceList(mesh_, faceType, maxEqOp<label>());
123 
124  // Now we have for faceType:
125  // 0 : face not in cellZone
126  // 1 : internal face or normal patch face
127  // 2 : coupled patch face or excluded patch face
128 
129  // Sort into lists per patch.
130 
131  internalFaces_.setSize(mesh_.nFaces());
132  label nInternal = 0;
133 
134  for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
135  {
136  if (faceType[facei] == 1)
137  {
138  internalFaces_[nInternal++] = facei;
139  }
140  }
141  internalFaces_.setSize(nInternal);
142 
143  labelList nIncludedFaces(patches.size(), 0);
144  labelList nExcludedFaces(patches.size(), 0);
145 
146  forAll(patches, patchi)
147  {
148  const polyPatch& pp = patches[patchi];
149 
150  forAll(pp, patchFacei)
151  {
152  const label facei = pp.start() + patchFacei;
153 
154  if (faceType[facei] == 1)
155  {
156  nIncludedFaces[patchi]++;
157  }
158  else if (faceType[facei] == 2)
159  {
160  nExcludedFaces[patchi]++;
161  }
162  }
163  }
164 
165  includedFaces_.setSize(patches.size());
166  excludedFaces_.setSize(patches.size());
167  forAll(nIncludedFaces, patchi)
168  {
169  includedFaces_[patchi].setSize(nIncludedFaces[patchi]);
170  excludedFaces_[patchi].setSize(nExcludedFaces[patchi]);
171  }
172  nIncludedFaces = 0;
173  nExcludedFaces = 0;
174 
175  forAll(patches, patchi)
176  {
177  const polyPatch& pp = patches[patchi];
178 
179  forAll(pp, patchFacei)
180  {
181  const label facei = pp.start() + patchFacei;
182 
183  if (faceType[facei] == 1)
184  {
185  includedFaces_[patchi][nIncludedFaces[patchi]++] = patchFacei;
186  }
187  else if (faceType[facei] == 2)
188  {
189  excludedFaces_[patchi][nExcludedFaces[patchi]++] = patchFacei;
190  }
191  }
192  }
193 
194 
195  if (debug)
196  {
197  faceSet internalFaces(mesh_, "internalFaces", internalFaces_);
198 
199  Pout<< "Writing " << internalFaces.size()
200  << " internal faces in MRF zone to faceSet "
201  << internalFaces.name() << endl;
202 
203  internalFaces.write();
204 
205  faceSet MRFFaces(mesh_, "includedFaces", 100);
206 
207  forAll(includedFaces_, patchi)
208  {
209  forAll(includedFaces_[patchi], i)
210  {
211  const label patchFacei = includedFaces_[patchi][i];
212  MRFFaces.insert(patches[patchi].start()+patchFacei);
213  }
214  }
215 
216  Pout<< "Writing " << MRFFaces.size()
217  << " patch faces in MRF zone to faceSet "
218  << MRFFaces.name() << endl;
219 
220  MRFFaces.write();
221 
222  faceSet excludedFaces(mesh_, "excludedFaces", 100);
223  forAll(excludedFaces_, patchi)
224  {
225  forAll(excludedFaces_[patchi], i)
226  {
227  const label patchFacei = excludedFaces_[patchi][i];
228  excludedFaces.insert(patches[patchi].start()+patchFacei);
229  }
230  }
231 
232  Pout<< "Writing " << excludedFaces.size()
233  << " faces in MRF zone with special handling to faceSet "
234  << excludedFaces.name() << endl;
235 
236  excludedFaces.write();
237  }
238 }
239 
240 
241 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
242 
244 (
245  const word& name,
246  const fvMesh& mesh,
247  const dictionary& dict,
248  const word& cellZoneName
249 )
250 :
251  mesh_(mesh),
252  name_(name),
253  coeffs_(dict),
254  cellZoneName_(cellZoneName),
255  cellZoneID_(),
256  excludedPatchNames_
257  (
258  wordReList(coeffs_.lookupOrDefault("nonRotatingPatches", wordReList()))
259  ),
260  origin_(coeffs_.lookup("origin")),
261  axis_(coeffs_.lookup("axis")),
262  omega_(Function1<scalar>::New("omega", coeffs_))
263 {
264  if (cellZoneName_ == word::null)
265  {
266  coeffs_.lookup("cellZone") >> cellZoneName_;
267  }
268 
269  cellZoneID_ = mesh_.cellZones().findZoneID(cellZoneName_);
270 
271  axis_ = axis_/mag(axis_);
272 
273  const labelHashSet excludedPatchSet
274  (
275  mesh_.boundaryMesh().patchSet(excludedPatchNames_)
276  );
277 
278  excludedPatchLabels_.setSize(excludedPatchSet.size());
279 
280  label i = 0;
281  forAllConstIter(labelHashSet, excludedPatchSet, iter)
282  {
283  excludedPatchLabels_[i++] = iter.key();
284  }
285 
286  bool cellZoneFound = (cellZoneID_ != -1);
287 
288  reduce(cellZoneFound, orOp<bool>());
289 
290  if (!cellZoneFound)
291  {
293  << "cannot find MRF cellZone " << cellZoneName_
294  << exit(FatalError);
295  }
296 
297  setMRFFaces();
298 }
299 
300 
301 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
302 
304 {
305  return omega_->value(mesh_.time().userTimeValue())*axis_;
306 }
307 
308 
310 (
311  const volVectorField& U,
312  volVectorField& ddtU
313 ) const
314 {
315  if (cellZoneID_ == -1)
316  {
317  return;
318  }
319 
320  const labelList& cells = mesh_.cellZones()[cellZoneID_];
321  vectorField& ddtUc = ddtU.primitiveFieldRef();
322  const vectorField& Uc = U;
323 
324  const vector Omega = this->Omega();
325 
326  forAll(cells, i)
327  {
328  const label celli = cells[i];
329  ddtUc[celli] += (Omega ^ Uc[celli]);
330  }
331 }
332 
333 
334 void Foam::MRFZone::addCoriolis(fvVectorMatrix& UEqn, const bool rhs) const
335 {
336  if (cellZoneID_ == -1)
337  {
338  return;
339  }
340 
341  const labelList& cells = mesh_.cellZones()[cellZoneID_];
342  const scalarField& V = mesh_.V();
343  vectorField& Usource = UEqn.source();
344  const vectorField& U = UEqn.psi();
345 
346  const vector Omega = this->Omega();
347 
348  if (rhs)
349  {
350  forAll(cells, i)
351  {
352  const label celli = cells[i];
353  Usource[celli] += V[celli]*(Omega ^ U[celli]);
354  }
355  }
356  else
357  {
358  forAll(cells, i)
359  {
360  const label celli = cells[i];
361  Usource[celli] -= V[celli]*(Omega ^ U[celli]);
362  }
363  }
364 }
365 
366 
368 (
369  const volScalarField& rho,
370  fvVectorMatrix& UEqn,
371  const bool rhs
372 ) const
373 {
374  if (cellZoneID_ == -1)
375  {
376  return;
377  }
378 
379  const labelList& cells = mesh_.cellZones()[cellZoneID_];
380  const scalarField& V = mesh_.V();
381  vectorField& Usource = UEqn.source();
382  const vectorField& U = UEqn.psi();
383 
384  const vector Omega = this->Omega();
385 
386  if (rhs)
387  {
388  forAll(cells, i)
389  {
390  const label celli = cells[i];
391  Usource[celli] += V[celli]*rho[celli]*(Omega ^ U[celli]);
392  }
393  }
394  else
395  {
396  forAll(cells, i)
397  {
398  const label celli = cells[i];
399  Usource[celli] -= V[celli]*rho[celli]*(Omega ^ U[celli]);
400  }
401  }
402 }
403 
404 
406 (
407  volVectorField& centrifugalAcceleration
408 ) const
409 {
410  if (cellZoneID_ == -1)
411  {
412  return;
413  }
414 
415  const labelList& cells = mesh_.cellZones()[cellZoneID_];
416  const volVectorField& C = mesh_.C();
417  vectorField& cac = centrifugalAcceleration.primitiveFieldRef();
418 
419  const vector Omega = this->Omega();
420 
421  forAll(cells, i)
422  {
423  const label celli = cells[i];
424  cac[celli] -= Omega ^ (Omega ^ (C[celli] - origin_));
425  }
426 
427  // Included (rotating) patches
428 
429  volVectorField::Boundary& caf = centrifugalAcceleration.boundaryFieldRef();
430 
431  forAll(includedFaces_, patchi)
432  {
433  forAll(includedFaces_[patchi], i)
434  {
435  const label patchFacei = includedFaces_[patchi][i];
436  caf[patchi][patchFacei] -=
437  Omega
438  ^ (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin_));
439  }
440  }
441 }
442 
443 
445 {
446  if (cellZoneID_ == -1)
447  {
448  return;
449  }
450 
451  const volVectorField& C = mesh_.C();
452  const labelList& cells = mesh_.cellZones()[cellZoneID_];
453 
454  const vector Omega = this->Omega();
455 
456  forAll(cells, i)
457  {
458  const label celli = cells[i];
459  U[celli] -= (Omega ^ (C[celli] - origin_));
460  }
461 
462  // Included patches
463 
465 
466  forAll(includedFaces_, patchi)
467  {
468  forAll(includedFaces_[patchi], i)
469  {
470  const label patchFacei = includedFaces_[patchi][i];
471  Ubf[patchi][patchFacei] = Zero;
472  }
473  }
474 
475  // Excluded patches
476  forAll(excludedFaces_, patchi)
477  {
478  forAll(excludedFaces_[patchi], i)
479  {
480  const label patchFacei = excludedFaces_[patchi][i];
481  Ubf[patchi][patchFacei] -=
482  (Omega
483  ^ (C.boundaryField()[patchi][patchFacei] - origin_));
484  }
485  }
486 }
487 
488 
490 {
491  makeRelativeRhoFlux(geometricOneField(), phi);
492 }
493 
494 
496 {
497  makeRelativeRhoFlux(oneFieldField(), phi);
498 }
499 
500 
501 void Foam::MRFZone::makeRelative(Field<scalar>& phi, const label patchi) const
502 {
503  makeRelativeRhoFlux(oneField(), phi, patchi);
504 }
505 
506 
508 (
509  const surfaceScalarField& rho,
510  surfaceScalarField& phi
511 ) const
512 {
513  makeRelativeRhoFlux(rho, phi);
514 }
515 
516 
518 {
519  if (cellZoneID_ == -1)
520  {
521  return;
522  }
523 
524  const volVectorField& C = mesh_.C();
525  const labelList& cells = mesh_.cellZones()[cellZoneID_];
526 
527  const vector Omega = this->Omega();
528 
529  forAll(cells, i)
530  {
531  const label celli = cells[i];
532  U[celli] += (Omega ^ (C[celli] - origin_));
533  }
534 
535  // Included patches
537 
538  forAll(includedFaces_, patchi)
539  {
540  forAll(includedFaces_[patchi], i)
541  {
542  const label patchFacei = includedFaces_[patchi][i];
543  Ubf[patchi][patchFacei] =
544  (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin_));
545  }
546  }
547 
548  // Excluded patches
549  forAll(excludedFaces_, patchi)
550  {
551  forAll(excludedFaces_[patchi], i)
552  {
553  const label patchFacei = excludedFaces_[patchi][i];
554  Ubf[patchi][patchFacei] +=
555  (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin_));
556  }
557  }
558 }
559 
560 
562 {
563  makeAbsoluteRhoFlux(geometricOneField(), phi);
564 }
565 
566 
568 (
569  const surfaceScalarField& rho,
570  surfaceScalarField& phi
571 ) const
572 {
573  makeAbsoluteRhoFlux(rho, phi);
574 }
575 
576 
578 {
579  const vector Omega = this->Omega();
580 
581  // Included patches
583 
584  forAll(includedFaces_, patchi)
585  {
586  const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];
587 
588  vectorField pfld(Ubf[patchi]);
589 
590  forAll(includedFaces_[patchi], i)
591  {
592  const label patchFacei = includedFaces_[patchi][i];
593 
594  pfld[patchFacei] = (Omega ^ (patchC[patchFacei] - origin_));
595  }
596 
597  Ubf[patchi] == pfld;
598  }
599 }
600 
601 
603 {
604  os << nl;
605  os.write(name_) << nl;
606  os << token::BEGIN_BLOCK << incrIndent << nl;
607  writeEntry(os, "cellZone", cellZoneName_);
608  writeEntry(os, "origin", origin_);
609  writeEntry(os, "axis", axis_);
610  omega_->write(os);
611 
612  if (excludedPatchNames_.size())
613  {
614  writeEntry(os, "nonRotatingPatches", excludedPatchNames_);
615  }
616 
617  os << decrIndent << token::END_BLOCK << nl;
618 }
619 
620 
622 {
623  coeffs_ = dict;
624 
625  coeffs_.lookup("cellZone") >> cellZoneName_;
626  cellZoneID_ = mesh_.cellZones().findZoneID(cellZoneName_);
627 
628  return true;
629 }
630 
631 
633 {
634  if (mesh_.topoChanged())
635  {
636  setMRFFaces();
637  }
638 }
639 
640 
641 // ************************************************************************* //
bool read(const dictionary &dict)
Read MRF dictionary.
Definition: MRFZone.C:621
Foam::surfaceFields.
void addCoriolis(const volVectorField &U, volVectorField &ddtU) const
Add the Coriolis force contribution to the acceleration field.
Definition: MRFZone.C:310
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:82
Run-time selectable general function of one variable.
Definition: Function1.H:52
void makeRelative(volVectorField &U) const
Make the given absolute velocity relative within the MRF region.
Definition: MRFZone.C:444
dictionary dict
Graphite solid properties.
Definition: C.H:48
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
virtual Ostream & write(const char)=0
Write character.
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
U
Definition: pEqn.H:72
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
MRFZone(const word &name, const fvMesh &mesh, const dictionary &dict, const word &cellZoneName=word::null)
Construct from fvMesh.
Definition: MRFZone.C:244
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
const Boundary & boundaryField() const
Return const-reference to the boundary field.
const GeometricField< Type, fvPatchField, volMesh > & psi() const
Definition: fvMatrix.H:290
void makeAbsolute(volVectorField &U) const
Make the given relative velocity absolute within the MRF region.
Definition: MRFZone.C:517
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
void addCentrifugalAcceleration(volVectorField &centrifugalAcceleration) const
Add the centrifugal acceleration.
Definition: MRFZone.C:406
void update()
Update MRFZone faces if the mesh topology changes.
Definition: MRFZone.C:632
A class representing the concept of a field of oneFields used to avoid unnecessary manipulations for ...
Definition: oneFieldField.H:50
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:211
C()
Construct null.
Definition: C.C:40
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
const cellShapeList & cells
A class for handling words, derived from string.
Definition: word.H:59
static const word null
An empty word.
Definition: word.H:77
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
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
static const char nl
Definition: Ostream.H:260
Field< Type > & source()
Definition: fvMatrix.H:300
defineTypeNameAndDebug(combustionModel, 0)
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:235
void writeEntry(Ostream &os, const HashTable< T, Key, Hash > &ht)
Definition: HashTableIO.C:96
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
void correctBoundaryVelocity(volVectorField &U) const
Correct the boundary velocity for the rotation of the MRF region.
Definition: MRFZone.C:577
vector Omega() const
Return the current Omega vector.
Definition: MRFZone.C:303
A class representing the concept of a field of 1 used to avoid unnecessary manipulations for objects ...
Definition: oneField.H:50
static void syncFaceList(const polyMesh &mesh, UList< T > &l, const CombineOp &cop)
Synchronise values on all mesh faces.
Definition: syncTools.H:387
label patchi
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:95
A special matrix type and solver, designed for finite volume solutions of scalar equations.
dimensioned< scalar > mag(const dimensioned< Type > &)
List< wordRe > wordReList
A List of wordRe (word or regular expression)
Definition: wordReList.H:50
void writeData(Ostream &os) const
Write.
Definition: MRFZone.C:602
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:228
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:864