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