InflationInjection.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-2013 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 "InflationInjection.H"
27 #include "mathematicalConstants.H"
28 #include "PackedBoolList.H"
29 #include "cellSet.H"
30 #include "ListListOps.H"
31 
32 using namespace Foam::constant::mathematical;
33 
34 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
35 
36 template<class CloudType>
38 (
39  const dictionary& dict,
40  CloudType& owner,
41  const word& modelName
42 )
43 :
44  InjectionModel<CloudType>(dict, owner, modelName, typeName),
45  generationSetName_(this->coeffDict().lookup("generationCellSet")),
46  inflationSetName_(this->coeffDict().lookup("inflationCellSet")),
47  generationCells_(),
48  inflationCells_(),
49  duration_(readScalar(this->coeffDict().lookup("duration"))),
50  flowRateProfile_
51  (
53  (
54  owner.db().time(),
55  "flowRateProfile",
56  this->coeffDict()
57  )
58  ),
59  growthRate_
60  (
62  (
63  owner.db().time(),
64  "growthRate",
65  this->coeffDict()
66  )
67  ),
68  newParticles_(),
69  volumeAccumulator_(0.0),
70  fraction_(1.0),
71  selfSeed_(this->coeffDict().lookupOrDefault("selfSeed", false)),
72  dSeed_(SMALL),
73  sizeDistribution_
74  (
76  (
77  this->coeffDict().subDict("sizeDistribution"),
78  owner.rndGen()
79  )
80  )
81 {
82  duration_ = owner.db().time().userTimeToTime(duration_);
83 
84  if (selfSeed_)
85  {
86  dSeed_ = readScalar(this->coeffDict().lookup("dSeed"));
87  }
88 
89  cellSet generationCells(this->owner().mesh(), generationSetName_);
90 
91  generationCells_ = generationCells.toc();
92 
93  cellSet inflationCells(this->owner().mesh(), inflationSetName_);
94 
95  // Union of cellSets
96  inflationCells |= generationCells;
97 
98  inflationCells_ = inflationCells.toc();
99 
100  if (Pstream::parRun())
101  {
102  scalar generationVolume = 0.0;
103 
104  forAll(generationCells_, gCI)
105  {
106  label cI = generationCells_[gCI];
107 
108  generationVolume += this->owner().mesh().cellVolumes()[cI];
109  }
110 
111  scalar totalGenerationVolume = generationVolume;
112 
113  reduce(totalGenerationVolume, sumOp<scalar>());
114 
115  fraction_ = generationVolume/totalGenerationVolume;
116  }
117 
118  // Set total volume/mass to inject
119  this->volumeTotal_ = fraction_*flowRateProfile_.integrate(0.0, duration_);
120  this->massTotal_ *= fraction_;
121 }
122 
123 
124 template<class CloudType>
126 (
128 )
129 :
131  generationSetName_(im.generationSetName_),
132  inflationSetName_(im.inflationSetName_),
133  generationCells_(im.generationCells_),
134  inflationCells_(im.inflationCells_),
135  duration_(im.duration_),
136  flowRateProfile_(im.flowRateProfile_),
137  growthRate_(im.growthRate_),
138  newParticles_(im.newParticles_),
139  volumeAccumulator_(im.volumeAccumulator_),
140  fraction_(im.fraction_),
141  selfSeed_(im.selfSeed_),
142  dSeed_(im.dSeed_),
143  sizeDistribution_(im.sizeDistribution_().clone().ptr())
144 {}
145 
146 
147 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
148 
149 template<class CloudType>
151 {}
152 
153 
154 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
155 
156 template<class CloudType>
158 {
159  // do nothing
160 }
161 
162 
163 template<class CloudType>
165 {
166  return this->SOI_ + duration_;
167 }
168 
169 
170 template<class CloudType>
172 (
173  const scalar time0,
174  const scalar time1
175 )
176 {
177  const polyMesh& mesh = this->owner().mesh();
178 
180  this->owner().cellOccupancy();
181 
182  scalar gR = growthRate_.value(time1);
183 
184  scalar dT = time1 - time0;
185 
186  // Inflate existing particles
187 
188  forAll(inflationCells_, iCI)
189  {
190  label cI = inflationCells_[iCI];
191 
192  typename CloudType::parcelType* pPtr = NULL;
193 
194  forAll(cellOccupancy[cI], cPI)
195  {
196  pPtr = cellOccupancy[cI][cPI];
197 
198  scalar dTarget = pPtr->dTarget();
199 
200  pPtr->d() = min(dTarget, pPtr->d() + gR*dT);
201  }
202  }
203 
204  // Generate new particles
205 
206  newParticles_.clear();
207 
208  cachedRandom& rnd = this->owner().rndGen();
209 
210  // Diameter factor, when splitting particles into 4, this is the
211  // factor that modifies the diameter.
212  scalar dFact = sqrt(2.0)/(sqrt(3.0) + sqrt(2.0));
213 
214  if ((time0 >= 0.0) && (time0 < duration_))
215  {
216  volumeAccumulator_ +=
217  fraction_*flowRateProfile_.integrate(time0, time1);
218  }
219 
220  labelHashSet cellCentresUsed;
221 
222  // Loop escape counter
223  label maxIterations = max
224  (
225  1,
226  (10*volumeAccumulator_)
227  /CloudType::parcelType::volume(sizeDistribution_().minValue())
228  );
229 
230  label iterationNo = 0;
231 
232  // Info<< "Accumulated volume to inject: "
233  // << returnReduce(volumeAccumulator_, sumOp<scalar>()) << endl;
234 
235  while (!generationCells_.empty() && volumeAccumulator_ > 0)
236  {
237  if (iterationNo > maxIterations)
238  {
239  WarningIn
240  (
241  "Foam::label "
242  "Foam::InflationInjection<CloudType>::parcelsToInject"
243  "("
244  "const scalar, "
245  "const scalar"
246  ")"
247  )
248  << "Maximum particle split iterations ("
249  << maxIterations << ") exceeded" << endl;
250 
251  break;
252  }
253 
254  label cI = generationCells_
255  [
256  rnd.position
257  (
258  label(0),
259  generationCells_.size() - 1
260  )
261  ];
262 
263  // Pick a particle at random from the cell - if there are
264  // none, insert one at the cell centre. Otherwise, split an
265  // existing particle into four new ones.
266 
267  if (cellOccupancy[cI].empty())
268  {
269  if (selfSeed_ && !cellCentresUsed.found(cI))
270  {
271  scalar dNew = sizeDistribution_().sample();
272 
273  newParticles_.append
274  (
276  (
278  Pair<scalar>(dSeed_, dNew)
279  )
280  );
281 
282  volumeAccumulator_ -= CloudType::parcelType::volume(dNew);
283 
284  cellCentresUsed.insert(cI);
285  }
286  }
287  else
288  {
289  label cPI = rnd.position(label(0), cellOccupancy[cI].size() - 1);
290 
291  // This has to be a reference to the pointer so that it
292  // can be set to NULL when the particle is deleted.
293  typename CloudType::parcelType*& pPtr = cellOccupancy[cI][cPI];
294 
295  if (pPtr != NULL)
296  {
297  scalar pD = pPtr->d();
298 
299  // Select bigger particles by preference
300  if ((pD/pPtr->dTarget()) < rnd.sample01<scalar>())
301  {
302  continue;
303  }
304 
305  const point& pP = pPtr->position();
306  const vector& pU = pPtr->U();
307 
308  // Generate a tetrahedron of new positions with the
309  // four new spheres fitting inside the old one, where
310  // a is the diameter of the new spheres, and is
311  // related to the diameter of the enclosing sphere, A,
312  // by a = sqrt(2)*A/(sqrt(3) + sqrt(2));
313 
314  // Positions around the origin, which is the
315  // tetrahedron centroid (centre of old sphere).
316 
317  // x = a/sqrt(3)
318  // r = a/(2*sqrt(6))
319  // R = sqrt(3)*a/(2*sqrt(2))
320  // d = a/(2*sqrt(3))
321 
322  // p0(x, 0, -r)
323  // p1(-d, a/2, -r)
324  // p2(-d, -a/2, -r)
325  // p3(0, 0, R)
326 
327  scalar a = pD*dFact;
328 
329  scalar x = a/sqrt(3.0);
330  scalar r = a/(2.0*sqrt(6.0));
331  scalar R = sqrt(3.0)*a/(2.0*sqrt(2.0));
332  scalar d = a/(2.0*sqrt(3.0));
333 
334  scalar dNew = sizeDistribution_().sample();
335  scalar volNew = CloudType::parcelType::volume(dNew);
336 
337  newParticles_.append
338  (
340  (
341  Pair<vector>(vector(x, 0, -r) + pP, pU),
342  Pair<scalar>(a, dNew)
343  )
344  );
345  volumeAccumulator_ -= volNew;
346 
347  dNew = sizeDistribution_().sample();
348  newParticles_.append
349  (
351  (
352  Pair<vector>(vector(-d, a/2, -r) + pP, pU),
353  Pair<scalar>(a, dNew)
354  )
355  );
356  volumeAccumulator_ -= volNew;
357 
358  dNew = sizeDistribution_().sample();
359  newParticles_.append
360  (
362  (
363  Pair<vector>(vector(-d, -a/2, -r) + pP, pU),
364  Pair<scalar>(a, dNew)
365  )
366  );
367  volumeAccumulator_ -= volNew;
368 
369  dNew = sizeDistribution_().sample();
370  newParticles_.append
371  (
373  (
374  Pair<vector>(vector(0, 0, R) + pP, pU),
375  Pair<scalar>(a, dNew)
376  )
377  );
378  volumeAccumulator_ -= volNew;
379 
380  // Account for the lost volume of the particle which
381  // is to be deleted
382  volumeAccumulator_ += CloudType::parcelType::volume
383  (
384  pPtr->dTarget()
385  );
386 
387  this->owner().deleteParticle(*pPtr);
388 
389  pPtr = NULL;
390  }
391  }
392 
393  iterationNo++;
394  }
395 
396  if (Pstream::parRun())
397  {
398  List<List<vectorPairScalarPair> > gatheredNewParticles
399  (
401  );
402 
403  gatheredNewParticles[Pstream::myProcNo()] = newParticles_;
404 
405  // Gather data onto master
406  Pstream::gatherList(gatheredNewParticles);
407 
408  // Combine
409  List<vectorPairScalarPair> combinedNewParticles
410  (
412  (
413  gatheredNewParticles,
415  )
416  );
417 
418  if (Pstream::master())
419  {
420  newParticles_ = combinedNewParticles;
421  }
422 
423  Pstream::scatter(newParticles_);
424  }
425 
426  return newParticles_.size();
427 }
428 
429 
430 template<class CloudType>
432 (
433  const scalar time0,
434  const scalar time1
435 )
436 {
437  if ((time0 >= 0.0) && (time0 < duration_))
438  {
439  return fraction_*flowRateProfile_.integrate(time0, time1);
440  }
441  else
442  {
443  return 0.0;
444  }
445 }
446 
447 
448 template<class CloudType>
450 (
451  const label parcelI,
452  const label,
453  const scalar,
454  vector& position,
455  label& cellOwner,
456  label& tetFaceI,
457  label& tetPtI
458 )
459 {
460  position = newParticles_[parcelI].first().first();
461 
462  this->findCellAtPosition
463  (
464  cellOwner,
465  tetFaceI,
466  tetPtI,
467  position,
468  false
469  );
470 }
471 
472 
473 template<class CloudType>
475 (
476  const label parcelI,
477  const label,
478  const scalar,
479  typename CloudType::parcelType& parcel
480 )
481 {
482  parcel.U() = newParticles_[parcelI].first().second();
483 
484  parcel.d() = newParticles_[parcelI].second().first();
485 
486  parcel.dTarget() = newParticles_[parcelI].second().second();
487 }
488 
489 
490 template<class CloudType>
492 {
493  return false;
494 }
495 
496 
497 template<class CloudType>
499 {
500  return true;
501 }
502 
503 
504 // ************************************************************************* //
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
dimensionedScalar sqrt(const dimensionedScalar &ds)
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:380
A 2-tuple for storing two objects of different types.
Definition: Tuple2.H:47
Inflation injection - creates new particles by splitting existing particles within in a set of genera...
virtual void setProperties(const label parcelI, const label nParcels, const scalar time, typename CloudType::parcelType &parcel)
Set the parcel properties.
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
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.
virtual bool fullyDescribed() const
Flag to identify whether model fully describes the parcel.
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
#define R(A, B, C, D, E, F, K, M)
A collection of cell labels.
Definition: cellSet.H:48
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
const vectorField & cellCentres() const
dynamicFvMesh & mesh
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
scalar minValue
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:217
AccessType combine(const List< T > &, AccessOp aop=accessOp< T >())
Combines sublists into one big list.
Definition: ListListOps.C:34
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
dictionary dict
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:386
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
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
virtual void updateMesh()
Set injector locations when mesh is updated.
Templated injection model class.
static void scatter(const List< commsStruct > &comms, T &Value, const int tag, const label comm)
Scatter data. Distribute without modification. Reverse of gather.
virtual scalar volumeToInject(const scalar time0, const scalar time1)
Volume of parcels to introduce relative to SOI.
static autoPtr< distributionModel > New(const dictionary &dict, cachedRandom &rndGen)
Selector.
List< Key > toc() const
Return the table of contents.
Definition: HashTable.C:201
scalar timeEnd() const
Return the end-of-injection time.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:404
virtual label parcelsToInject(const scalar time0, const scalar time1)
Number of parcels to introduce relative to SOI.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
const List< DynamicList< molecule * > > & cellOccupancy
Random number generator.
Definition: cachedRandom.H:63
static const Vector zero
Definition: Vector.H:80
Type sample01()
Return a sample whose components lie in the range 0-1.
virtual bool validInjection(const label parcelI)
Return flag to identify whether or not injection of parcelI is.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:398
Type position(const Type &start, const Type &end)
Return a sample between start and end.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:68
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
InflationInjection(const dictionary &dict, CloudType &owner, const word &modelName)
Construct from dictionary.
mathematical constants.
virtual ~InflationInjection()
Destructor.
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:109
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116