CellZoneInjection.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 "CellZoneInjection.H"
27 #include "mathematicalConstants.H"
29 #include "globalIndex.H"
30 #include "Pstream.H"
31 
32 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
33 
34 template<class CloudType>
36 (
37  const labelList& cellZoneCells
38 )
39 {
40  const fvMesh& mesh = this->owner().mesh();
41  const scalarField& V = mesh.V();
42  const label nCells = cellZoneCells.size();
43  Random& rnd = this->owner().rndGen();
44 
45  DynamicList<vector> positions(nCells); // initial size only
46  DynamicList<label> injectorCells(nCells); // initial size only
47  DynamicList<label> injectorTetFaces(nCells); // initial size only
48  DynamicList<label> injectorTetPts(nCells); // initial size only
49 
50  scalar newParticlesTotal = 0.0;
51  label addParticlesTotal = 0;
52 
53  forAll(cellZoneCells, i)
54  {
55  const label celli = cellZoneCells[i];
56 
57  // Calc number of particles to add
58  const scalar newParticles = V[celli]*numberDensity_;
59  newParticlesTotal += newParticles;
60  label addParticles = floor(newParticles);
61  addParticlesTotal += addParticles;
62 
63  const scalar diff = newParticlesTotal - addParticlesTotal;
64  if (diff > 1)
65  {
66  label corr = floor(diff);
67  addParticles += corr;
68  addParticlesTotal += corr;
69  }
70 
71  // Construct cell tet indices
72  const List<tetIndices> cellTetIs =
73  polyMeshTetDecomposition::cellTetIndices(mesh, celli);
74 
75  // Construct cell tet volume fractions
76  scalarField cTetVFrac(cellTetIs.size());
77  cTetVFrac[0] = cellTetIs[0].tet(mesh).mag();
78  for (label tetI = 1; tetI < cellTetIs.size(); tetI++)
79  {
80  cTetVFrac[tetI] =
81  cTetVFrac[tetI-1] + cellTetIs[tetI].tet(mesh).mag();
82  }
83  cTetVFrac /= cTetVFrac.last();
84 
85  // Set new particle position and cellId
86  for (label pI = 0; pI < addParticles; pI++)
87  {
88  const scalar volFrac = rnd.sample01<scalar>();
89  label tetI = 0;
90  forAll(cTetVFrac, vfI)
91  {
92  if (cTetVFrac[vfI] > volFrac)
93  {
94  tetI = vfI;
95  break;
96  }
97  }
98  positions.append(cellTetIs[tetI].tet(mesh).randomPoint(rnd));
99 
100  injectorCells.append(celli);
101  injectorTetFaces.append(cellTetIs[tetI].face());
102  injectorTetPts.append(cellTetIs[tetI].tetPt());
103  }
104  }
105 
106  // Parallel operation manipulations
107  globalIndex globalPositions(positions.size());
108  List<vector> allPositions(globalPositions.size(), point::max);
109  List<label> allInjectorCells(globalPositions.size(), -1);
110  List<label> allInjectorTetFaces(globalPositions.size(), -1);
111  List<label> allInjectorTetPts(globalPositions.size(), -1);
112 
113  // Gather all positions on to all processors
114  SubList<vector>
115  (
116  allPositions,
117  globalPositions.localSize(Pstream::myProcNo()),
118  globalPositions.offset(Pstream::myProcNo())
119  ) = positions;
120 
121  Pstream::listCombineGather(allPositions, minEqOp<point>());
122  Pstream::listCombineScatter(allPositions);
123 
124  // Gather local cell tet and tet-point Ids, but leave non-local ids set -1
125  SubList<label>
126  (
127  allInjectorCells,
128  globalPositions.localSize(Pstream::myProcNo()),
129  globalPositions.offset(Pstream::myProcNo())
130  ) = injectorCells;
131  SubList<label>
132  (
133  allInjectorTetFaces,
134  globalPositions.localSize(Pstream::myProcNo()),
135  globalPositions.offset(Pstream::myProcNo())
136  ) = injectorTetFaces;
137  SubList<label>
138  (
139  allInjectorTetPts,
140  globalPositions.localSize(Pstream::myProcNo()),
141  globalPositions.offset(Pstream::myProcNo())
142  ) = injectorTetPts;
143 
144  // Transfer data
145  positions_.transfer(allPositions);
146  injectorCells_.transfer(allInjectorCells);
147  injectorTetFaces_.transfer(allInjectorTetFaces);
148  injectorTetPts_.transfer(allInjectorTetPts);
149 
150 
151  if (debug)
152  {
153  OFstream points("points.obj");
154  forAll(positions_, i)
155  {
156  meshTools::writeOBJ(points, positions_[i]);
157  }
158  }
159 }
160 
161 
162 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
163 
164 template<class CloudType>
166 (
167  const dictionary& dict,
168  CloudType& owner,
169  const word& modelName
170 )
171 :
172  InjectionModel<CloudType>(dict, owner, modelName, typeName),
173  cellZoneName_(this->coeffDict().lookup("cellZone")),
174  numberDensity_(this->coeffDict().template lookup<scalar>("numberDensity")),
175  positions_(),
176  injectorCells_(),
177  injectorTetFaces_(),
178  injectorTetPts_(),
179  diameters_(),
180  U0_(this->coeffDict().lookup("U0")),
181  sizeDistribution_
182  (
184  (
185  this->coeffDict().subDict("sizeDistribution"), owner.rndGen()
186  )
187  )
188 {
189  topoChange();
190 }
191 
192 
193 template<class CloudType>
195 (
197 )
198 :
200  cellZoneName_(im.cellZoneName_),
201  numberDensity_(im.numberDensity_),
202  positions_(im.positions_),
203  injectorCells_(im.injectorCells_),
204  injectorTetFaces_(im.injectorTetFaces_),
205  injectorTetPts_(im.injectorTetPts_),
206  diameters_(im.diameters_),
207  U0_(im.U0_),
208  sizeDistribution_(im.sizeDistribution_().clone().ptr())
209 {}
210 
211 
212 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
213 
214 template<class CloudType>
216 {}
217 
218 
219 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
220 
221 template<class CloudType>
223 {
224  // Set/cache the injector cells
225  const fvMesh& mesh = this->owner().mesh();
226  const label zoneI = mesh.cellZones().findZoneID(cellZoneName_);
227 
228  if (zoneI < 0)
229  {
231  << "Unknown cell zone name: " << cellZoneName_
232  << ". Valid cell zones are: " << mesh.cellZones().names()
233  << nl << exit(FatalError);
234  }
235 
236  const labelList& cellZoneCells = mesh.cellZones()[zoneI];
237  const label nCells = cellZoneCells.size();
238  const scalar nCellsTotal = returnReduce(nCells, sumOp<label>());
239  const scalar VCells = sum(scalarField(mesh.V(), cellZoneCells));
240  const scalar VCellsTotal = returnReduce(VCells, sumOp<scalar>());
241  Info<< " cell zone size = " << nCellsTotal << endl;
242  Info<< " cell zone volume = " << VCellsTotal << endl;
243 
244  if ((nCellsTotal == 0) || (VCellsTotal*numberDensity_ < 1))
245  {
247  << "Number of particles to be added to cellZone " << cellZoneName_
248  << " is zero" << endl;
249  }
250  else
251  {
252  setPositions(cellZoneCells);
253 
254  Info<< " number density = " << numberDensity_ << nl
255  << " number of particles = " << positions_.size() << endl;
256 
257  // Construct parcel diameters
258  diameters_.setSize(positions_.size());
259  forAll(diameters_, i)
260  {
261  diameters_[i] = sizeDistribution_->sample();
262  }
263  }
264 
265  // Determine volume of particles to inject
266  this->volumeTotal_ = sum(pow3(diameters_))*constant::mathematical::pi/6.0;
267 }
268 
269 
270 template<class CloudType>
272 {
273  // Not used
274  return this->SOI_;
275 }
276 
277 
278 template<class CloudType>
280 (
281  const scalar time0,
282  const scalar time1
283 )
284 {
285  if ((0.0 >= time0) && (0.0 < time1))
286  {
287  return positions_.size();
288  }
289  else
290  {
291  return 0;
292  }
293 }
294 
295 
296 template<class CloudType>
298 (
299  const scalar time0,
300  const scalar time1
301 )
302 {
303  // All parcels introduced at SOI
304  if ((0.0 >= time0) && (0.0 < time1))
305  {
306  return this->volumeTotal_;
307  }
308  else
309  {
310  return 0.0;
311  }
312 }
313 
314 
315 template<class CloudType>
317 (
318  const label parcelI,
319  const label nParcels,
320  const scalar time,
321  vector& position,
322  label& cellOwner,
323  label& tetFacei,
324  label& tetPti
325 )
326 {
327  position = positions_[parcelI];
328  cellOwner = injectorCells_[parcelI];
329  tetFacei = injectorTetFaces_[parcelI];
330  tetPti = injectorTetPts_[parcelI];
331 }
332 
333 
334 template<class CloudType>
336 (
337  const label parcelI,
338  const label,
339  const scalar,
340  typename CloudType::parcelType& parcel
341 )
342 {
343  // set particle velocity
344  parcel.U() = U0_;
345 
346  // set particle diameter
347  parcel.d() = diameters_[parcelI];
348 }
349 
350 
351 template<class CloudType>
353 {
354  return false;
355 }
356 
357 
358 template<class CloudType>
360 {
361  return true;
362 }
363 
364 
365 // ************************************************************************* //
dictionary dict
label findZoneID(const word &zoneName) const
Find zone index given a name.
Definition: MeshZones.C:341
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
virtual ~CellZoneInjection()
Destructor.
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
autoPtr< CompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const viscosity &viscosity)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const meshCellZones & cellZones() const
Return cell zones.
Definition: polyMesh.H:501
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
Templated injection model class.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
wordList names() const
Return a list of zone names.
Definition: MeshZones.C:256
virtual void topoChange()
Set injector locations when mesh is updated.
stressControl lookup("compactNormalStress") >> compactNormalStress
virtual bool fullyDescribed() const
Flag to identify whether model fully describes the parcel.
const pointField & points
Injection positions specified by a particle number density within a cell set.
A class for handling words, derived from string.
Definition: word.H:59
List< label > labelList
A List of labels.
Definition: labelList.H:56
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:221
volScalarField scalarField(fieldObject, mesh)
static const char nl
Definition: Ostream.H:260
label parcelsToInject(const scalar time0, const scalar time1)
Number of parcels to introduce relative to SOI.
CellZoneInjection(const dictionary &dict, CloudType &owner, const word &modelName)
Construct from dictionary.
dimensionedScalar pow3(const dimensionedScalar &ds)
#define WarningInFunction
Report a warning using Foam::Warning.
virtual void setProperties(const label parcelI, const label nParcels, const scalar time, typename CloudType::parcelType &parcel)
Set the parcel properties.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:95
messageStream Info
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
scalar timeEnd() const
Return the end-of-injection time.
virtual void setPositionAndCell(const label parcelI, const label nParcels, const scalar time, vector &position, label &cellOwner, label &tetFacei, label &tetPti)
Set the injection position and owner cell, tetFace and tetPt.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
scalar volumeToInject(const scalar time0, const scalar time1)
Volume of parcels to introduce relative to SOI.
virtual bool validInjection(const label parcelI)
Return flag to identify whether or not injection of parcelI is.