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-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 "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  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  label facei = pp.start()+i;
111 
112  if (zoneCell[own[facei]])
113  {
114  faceType[facei] = 1;
115  nZoneFaces++;
116  }
117  }
118  }
119  }
120 
121  // Synchronize 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  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  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  Pout<< "Writing " << internalFaces.size()
199  << " internal faces in MRF zone to faceSet "
200  << internalFaces.name() << endl;
201  internalFaces.write();
202 
203  faceSet MRFFaces(mesh_, "includedFaces", 100);
204  forAll(includedFaces_, patchi)
205  {
206  forAll(includedFaces_[patchi], i)
207  {
208  label patchFacei = includedFaces_[patchi][i];
209  MRFFaces.insert(patches[patchi].start()+patchFacei);
210  }
211  }
212  Pout<< "Writing " << MRFFaces.size()
213  << " patch faces in MRF zone to faceSet "
214  << MRFFaces.name() << endl;
215  MRFFaces.write();
216 
217  faceSet excludedFaces(mesh_, "excludedFaces", 100);
218  forAll(excludedFaces_, patchi)
219  {
220  forAll(excludedFaces_[patchi], i)
221  {
222  label patchFacei = excludedFaces_[patchi][i];
223  excludedFaces.insert(patches[patchi].start()+patchFacei);
224  }
225  }
226  Pout<< "Writing " << excludedFaces.size()
227  << " faces in MRF zone with special handling to faceSet "
228  << excludedFaces.name() << endl;
229  excludedFaces.write();
230  }
231 }
232 
233 
234 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
235 
237 (
238  const word& name,
239  const fvMesh& mesh,
240  const dictionary& dict,
241  const word& cellZoneName
242 )
243 :
244  mesh_(mesh),
245  name_(name),
246  coeffs_(dict),
247  active_(coeffs_.lookupOrDefault("active", true)),
248  cellZoneName_(cellZoneName),
249  cellZoneID_(),
250  excludedPatchNames_
251  (
252  wordReList(coeffs_.lookupOrDefault("nonRotatingPatches", wordReList()))
253  ),
254  origin_(coeffs_.lookup("origin")),
255  axis_(coeffs_.lookup("axis")),
256  omega_(Function1<scalar>::New("omega", coeffs_))
257 {
258  if (cellZoneName_ == word::null)
259  {
260  coeffs_.lookup("cellZone") >> cellZoneName_;
261  }
262 
263  if (!active_)
264  {
265  cellZoneID_ = -1;
266  }
267  else
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 
302 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
303 
305 {
306  return omega_->value(mesh_.time().timeOutputValue())*axis_;
307 }
308 
309 
311 (
312  const volVectorField& U,
313  volVectorField& ddtU
314 ) const
315 {
316  if (cellZoneID_ == -1)
317  {
318  return;
319  }
320 
321  const labelList& cells = mesh_.cellZones()[cellZoneID_];
322  vectorField& ddtUc = ddtU.primitiveFieldRef();
323  const vectorField& Uc = U;
324 
325  const vector Omega = this->Omega();
326 
327  forAll(cells, i)
328  {
329  label celli = cells[i];
330  ddtUc[celli] += (Omega ^ Uc[celli]);
331  }
332 }
333 
334 
335 void Foam::MRFZone::addCoriolis(fvVectorMatrix& UEqn, const bool rhs) const
336 {
337  if (cellZoneID_ == -1)
338  {
339  return;
340  }
341 
342  const labelList& cells = mesh_.cellZones()[cellZoneID_];
343  const scalarField& V = mesh_.V();
344  vectorField& Usource = UEqn.source();
345  const vectorField& U = UEqn.psi();
346 
347  const vector Omega = this->Omega();
348 
349  if (rhs)
350  {
351  forAll(cells, i)
352  {
353  label celli = cells[i];
354  Usource[celli] += V[celli]*(Omega ^ U[celli]);
355  }
356  }
357  else
358  {
359  forAll(cells, i)
360  {
361  label celli = cells[i];
362  Usource[celli] -= V[celli]*(Omega ^ U[celli]);
363  }
364  }
365 }
366 
367 
369 (
370  const volScalarField& rho,
371  fvVectorMatrix& UEqn,
372  const bool rhs
373 ) const
374 {
375  if (cellZoneID_ == -1)
376  {
377  return;
378  }
379 
380  const labelList& cells = mesh_.cellZones()[cellZoneID_];
381  const scalarField& V = mesh_.V();
382  vectorField& Usource = UEqn.source();
383  const vectorField& U = UEqn.psi();
384 
385  const vector Omega = this->Omega();
386 
387  if (rhs)
388  {
389  forAll(cells, i)
390  {
391  label celli = cells[i];
392  Usource[celli] += V[celli]*rho[celli]*(Omega ^ U[celli]);
393  }
394  }
395  else
396  {
397  forAll(cells, i)
398  {
399  label celli = cells[i];
400  Usource[celli] -= V[celli]*rho[celli]*(Omega ^ U[celli]);
401  }
402  }
403 }
404 
405 
407 {
408  if (cellZoneID_ == -1)
409  {
410  return;
411  }
412 
413  const volVectorField& C = mesh_.C();
414 
415  const vector Omega = this->Omega();
416 
417  const labelList& cells = mesh_.cellZones()[cellZoneID_];
418 
419  forAll(cells, i)
420  {
421  label celli = cells[i];
422  U[celli] -= (Omega ^ (C[celli] - origin_));
423  }
424 
425  // Included patches
426 
427  volVectorField::Boundary& Ubf = U.boundaryFieldRef();
428 
429  forAll(includedFaces_, patchi)
430  {
431  forAll(includedFaces_[patchi], i)
432  {
433  label patchFacei = includedFaces_[patchi][i];
434  Ubf[patchi][patchFacei] = Zero;
435  }
436  }
437 
438  // Excluded patches
439  forAll(excludedFaces_, patchi)
440  {
441  forAll(excludedFaces_[patchi], i)
442  {
443  label patchFacei = excludedFaces_[patchi][i];
444  Ubf[patchi][patchFacei] -=
445  (Omega
446  ^ (C.boundaryField()[patchi][patchFacei] - origin_));
447  }
448  }
449 }
450 
451 
453 {
454  makeRelativeRhoFlux(geometricOneField(), phi);
455 }
456 
457 
459 {
460  makeRelativeRhoFlux(oneFieldField(), phi);
461 }
462 
463 
464 void Foam::MRFZone::makeRelative(Field<scalar>& phi, const label patchi) const
465 {
466  makeRelativeRhoFlux(oneField(), phi, patchi);
467 }
468 
469 
471 (
472  const surfaceScalarField& rho,
473  surfaceScalarField& phi
474 ) const
475 {
476  makeRelativeRhoFlux(rho, phi);
477 }
478 
479 
481 {
482  if (cellZoneID_ == -1)
483  {
484  return;
485  }
486 
487  const volVectorField& C = mesh_.C();
488 
489  const vector Omega = this->Omega();
490 
491  const labelList& cells = mesh_.cellZones()[cellZoneID_];
492 
493  forAll(cells, i)
494  {
495  label celli = cells[i];
496  U[celli] += (Omega ^ (C[celli] - origin_));
497  }
498 
499  // Included patches
500  volVectorField::Boundary& Ubf = U.boundaryFieldRef();
501 
502  forAll(includedFaces_, patchi)
503  {
504  forAll(includedFaces_[patchi], i)
505  {
506  label patchFacei = includedFaces_[patchi][i];
507  Ubf[patchi][patchFacei] =
508  (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin_));
509  }
510  }
511 
512  // Excluded patches
513  forAll(excludedFaces_, patchi)
514  {
515  forAll(excludedFaces_[patchi], i)
516  {
517  label patchFacei = excludedFaces_[patchi][i];
518  Ubf[patchi][patchFacei] +=
519  (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin_));
520  }
521  }
522 }
523 
524 
526 {
527  makeAbsoluteRhoFlux(geometricOneField(), phi);
528 }
529 
530 
532 (
533  const surfaceScalarField& rho,
534  surfaceScalarField& phi
535 ) const
536 {
537  makeAbsoluteRhoFlux(rho, phi);
538 }
539 
540 
542 {
543  if (!active_)
544  {
545  return;
546  }
547 
548  const vector Omega = this->Omega();
549 
550  // Included patches
551  volVectorField::Boundary& Ubf = U.boundaryFieldRef();
552 
553  forAll(includedFaces_, patchi)
554  {
555  const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];
556 
557  vectorField pfld(Ubf[patchi]);
558 
559  forAll(includedFaces_[patchi], i)
560  {
561  label patchFacei = includedFaces_[patchi][i];
562 
563  pfld[patchFacei] = (Omega ^ (patchC[patchFacei] - origin_));
564  }
565 
566  Ubf[patchi] == pfld;
567  }
568 }
569 
570 
572 {
573  os << nl;
574  os.write(name_) << nl;
575  os << token::BEGIN_BLOCK << incrIndent << nl;
576  writeEntry(os, "active", active_);
577  writeEntry(os, "cellZone", cellZoneName_);
578  writeEntry(os, "origin", origin_);
579  writeEntry(os, "axis", axis_);
580  omega_->writeData(os);
581 
582  if (excludedPatchNames_.size())
583  {
584  writeEntry(os, "nonRotatingPatches", excludedPatchNames_);
585  }
586 
587  os << decrIndent << token::END_BLOCK << nl;
588 }
589 
590 
592 {
593  coeffs_ = dict;
594 
595  active_ = coeffs_.lookupOrDefault("active", true);
596  coeffs_.lookup("cellZone") >> cellZoneName_;
597  cellZoneID_ = mesh_.cellZones().findZoneID(cellZoneName_);
598 
599  return true;
600 }
601 
602 
604 {
605  if (mesh_.topoChanging())
606  {
607  setMRFFaces();
608  }
609 }
610 
611 
612 // ************************************************************************* //
bool read(const dictionary &dict)
Read MRF dictionary.
Definition: MRFZone.C:591
Foam::surfaceFields.
void addCoriolis(const volVectorField &U, volVectorField &ddtU) const
Add the Coriolis force contribution to the acceleration field.
Definition: MRFZone.C:311
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:82
Top level data entry class for use in dictionaries. Provides a mechanism to specify a variable as a c...
Definition: Function1.H:53
void makeRelative(volVectorField &U) const
Make the given absolute velocity relative within the MRF region.
Definition: MRFZone.C:406
dictionary dict
Graphite solid properties.
Definition: C.H:48
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
MRFZone(const word &name, const fvMesh &mesh, const dictionary &dict, const word &cellZoneName=word::null)
Construct from fvMesh.
Definition: MRFZone.C:237
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const Boundary & boundaryField() const
Return const-reference to the boundary field.
const GeometricField< Type, fvPatchField, volMesh > & psi() const
Definition: fvMatrix.H:282
void makeAbsolute(volVectorField &U) const
Make the given relative velocity absolute within the MRF region.
Definition: MRFZone.C:480
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
void update()
Update MRFZone faces if the mesh topology changes.
Definition: MRFZone.C:603
A class representing the concept of a field of oneFields used to avoid unnecessary manipulations for ...
Definition: oneFieldField.H:50
Generic field type.
Definition: FieldField.H:51
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:208
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
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
static const char nl
Definition: Ostream.H:265
Field< Type > & source()
Definition: fvMatrix.H:292
defineTypeNameAndDebug(combustionModel, 0)
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:240
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:541
vector Omega() const
Return the current Omega vector.
Definition: MRFZone.C:304
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)
Synchronize values on all mesh faces.
Definition: syncTools.H:381
label patchi
U
Definition: pEqn.H:72
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.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
A special matrix type and solver, designed for finite volume solutions of scalar equations.
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual Ostream & write(const token &)=0
Write next token to stream.
List< wordRe > wordReList
A List of wordRe (word or regular expression)
Definition: wordReList.H:50
void writeData(Ostream &os) const
Write.
Definition: MRFZone.C:571
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:233
Namespace for OpenFOAM.