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-2017 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 
160 
161 template<class CloudType>
163 {
164  return this->SOI_ + duration_;
165 }
166 
167 
168 template<class CloudType>
170 (
171  const scalar time0,
172  const scalar time1
173 )
174 {
175  const polyMesh& mesh = this->owner().mesh();
176 
178  this->owner().cellOccupancy();
179 
180  scalar gR = growthRate_.value(time1);
181 
182  scalar dT = time1 - time0;
183 
184  // Inflate existing particles
185 
186  forAll(inflationCells_, iCI)
187  {
188  label cI = inflationCells_[iCI];
189 
190  typename CloudType::parcelType* pPtr = nullptr;
191 
192  forAll(cellOccupancy[cI], cPI)
193  {
194  pPtr = cellOccupancy[cI][cPI];
195 
196  scalar dTarget = pPtr->dTarget();
197 
198  pPtr->d() = min(dTarget, pPtr->d() + gR*dT);
199  }
200  }
201 
202  // Generate new particles
203 
204  newParticles_.clear();
205 
206  cachedRandom& rnd = this->owner().rndGen();
207 
208  // Diameter factor, when splitting particles into 4, this is the
209  // factor that modifies the diameter.
210  scalar dFact = sqrt(2.0)/(sqrt(3.0) + sqrt(2.0));
211 
212  if ((time0 >= 0.0) && (time0 < duration_))
213  {
214  volumeAccumulator_ +=
215  fraction_*flowRateProfile_.integrate(time0, time1);
216  }
217 
218  labelHashSet cellCentresUsed;
219 
220  // Loop escape counter
221  label maxIterations = max
222  (
223  1,
224  (10*volumeAccumulator_)
225  /CloudType::parcelType::volume(sizeDistribution_().minValue())
226  );
227 
228  label iterationNo = 0;
229 
230  // Info<< "Accumulated volume to inject: "
231  // << returnReduce(volumeAccumulator_, sumOp<scalar>()) << endl;
232 
233  while (!generationCells_.empty() && volumeAccumulator_ > 0)
234  {
235  if (iterationNo > maxIterations)
236  {
238  << "Maximum particle split iterations ("
239  << maxIterations << ") exceeded" << endl;
240 
241  break;
242  }
243 
244  label cI = generationCells_
245  [
246  rnd.position
247  (
248  label(0),
249  generationCells_.size() - 1
250  )
251  ];
252 
253  // Pick a particle at random from the cell - if there are
254  // none, insert one at the cell centre. Otherwise, split an
255  // existing particle into four new ones.
256 
257  if (cellOccupancy[cI].empty())
258  {
259  if (selfSeed_ && !cellCentresUsed.found(cI))
260  {
261  scalar dNew = sizeDistribution_().sample();
262 
263  newParticles_.append
264  (
266  (
267  Pair<vector>(mesh.cellCentres()[cI], Zero),
268  Pair<scalar>(dSeed_, dNew)
269  )
270  );
271 
272  volumeAccumulator_ -= CloudType::parcelType::volume(dNew);
273 
274  cellCentresUsed.insert(cI);
275  }
276  }
277  else
278  {
279  label cPI = rnd.position(label(0), cellOccupancy[cI].size() - 1);
280 
281  // This has to be a reference to the pointer so that it
282  // can be set to nullptr when the particle is deleted.
283  typename CloudType::parcelType*& pPtr = cellOccupancy[cI][cPI];
284 
285  if (pPtr != nullptr)
286  {
287  scalar pD = pPtr->d();
288 
289  // Select bigger particles by preference
290  if ((pD/pPtr->dTarget()) < rnd.sample01<scalar>())
291  {
292  continue;
293  }
294 
295  const point& pP = pPtr->position();
296  const vector& pU = pPtr->U();
297 
298  // Generate a tetrahedron of new positions with the
299  // four new spheres fitting inside the old one, where
300  // a is the diameter of the new spheres, and is
301  // related to the diameter of the enclosing sphere, A,
302  // by a = sqrt(2)*A/(sqrt(3) + sqrt(2));
303 
304  // Positions around the origin, which is the
305  // tetrahedron centroid (centre of old sphere).
306 
307  // x = a/sqrt(3)
308  // r = a/(2*sqrt(6))
309  // R = sqrt(3)*a/(2*sqrt(2))
310  // d = a/(2*sqrt(3))
311 
312  // p0(x, 0, -r)
313  // p1(-d, a/2, -r)
314  // p2(-d, -a/2, -r)
315  // p3(0, 0, R)
316 
317  scalar a = pD*dFact;
318 
319  scalar x = a/sqrt(3.0);
320  scalar r = a/(2.0*sqrt(6.0));
321  scalar R = sqrt(3.0)*a/(2.0*sqrt(2.0));
322  scalar d = a/(2.0*sqrt(3.0));
323 
324  scalar dNew = sizeDistribution_().sample();
325  scalar volNew = CloudType::parcelType::volume(dNew);
326 
327  newParticles_.append
328  (
330  (
331  Pair<vector>(vector(x, 0, -r) + pP, pU),
332  Pair<scalar>(a, dNew)
333  )
334  );
335  volumeAccumulator_ -= volNew;
336 
337  dNew = sizeDistribution_().sample();
338  newParticles_.append
339  (
341  (
342  Pair<vector>(vector(-d, a/2, -r) + pP, pU),
343  Pair<scalar>(a, dNew)
344  )
345  );
346  volumeAccumulator_ -= volNew;
347 
348  dNew = sizeDistribution_().sample();
349  newParticles_.append
350  (
352  (
353  Pair<vector>(vector(-d, -a/2, -r) + pP, pU),
354  Pair<scalar>(a, dNew)
355  )
356  );
357  volumeAccumulator_ -= volNew;
358 
359  dNew = sizeDistribution_().sample();
360  newParticles_.append
361  (
363  (
364  Pair<vector>(vector(0, 0, R) + pP, pU),
365  Pair<scalar>(a, dNew)
366  )
367  );
368  volumeAccumulator_ -= volNew;
369 
370  // Account for the lost volume of the particle which
371  // is to be deleted
372  volumeAccumulator_ += CloudType::parcelType::volume
373  (
374  pPtr->dTarget()
375  );
376 
377  this->owner().deleteParticle(*pPtr);
378 
379  pPtr = nullptr;
380  }
381  }
382 
383  iterationNo++;
384  }
385 
386  if (Pstream::parRun())
387  {
388  List<List<vectorPairScalarPair>> gatheredNewParticles
389  (
391  );
392 
393  gatheredNewParticles[Pstream::myProcNo()] = newParticles_;
394 
395  // Gather data onto master
396  Pstream::gatherList(gatheredNewParticles);
397 
398  // Combine
399  List<vectorPairScalarPair> combinedNewParticles
400  (
402  (
403  gatheredNewParticles,
405  )
406  );
407 
408  if (Pstream::master())
409  {
410  newParticles_ = combinedNewParticles;
411  }
412 
413  Pstream::scatter(newParticles_);
414  }
415 
416  return newParticles_.size();
417 }
418 
419 
420 template<class CloudType>
422 (
423  const scalar time0,
424  const scalar time1
425 )
426 {
427  if ((time0 >= 0.0) && (time0 < duration_))
428  {
429  return fraction_*flowRateProfile_.integrate(time0, time1);
430  }
431  else
432  {
433  return 0.0;
434  }
435 }
436 
437 
438 template<class CloudType>
440 (
441  const label parcelI,
442  const label,
443  const scalar,
444  vector& position,
445  label& cellOwner,
446  label& tetFacei,
447  label& tetPti
448 )
449 {
450  position = newParticles_[parcelI].first().first();
451 
452  this->findCellAtPosition
453  (
454  cellOwner,
455  tetFacei,
456  tetPti,
457  position,
458  false
459  );
460 }
461 
462 
463 template<class CloudType>
465 (
466  const label parcelI,
467  const label,
468  const scalar,
469  typename CloudType::parcelType& parcel
470 )
471 {
472  parcel.U() = newParticles_[parcelI].first().second();
473 
474  parcel.d() = newParticles_[parcelI].second().first();
475 
476  parcel.dTarget() = newParticles_[parcelI].second().second();
477 }
478 
479 
480 template<class CloudType>
482 {
483  return false;
484 }
485 
486 
487 template<class CloudType>
489 {
490  return true;
491 }
492 
493 
494 // ************************************************************************* //
virtual scalar volumeToInject(const scalar time0, const scalar time1)
Volume of parcels to introduce relative to SOI.
virtual ~InflationInjection()
Destructor.
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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.
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
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
A 2-tuple for storing two objects of different types.
Definition: HashTable.H:66
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:60
virtual bool fullyDescribed() const
Flag to identify whether model fully describes the parcel.
scalar minValue
Templated injection model class.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:418
virtual void updateMesh()
Set injector locations when mesh is updated.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:412
Random number generator.
Definition: cachedRandom.H:63
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116
AccessType combine(const List< T > &, AccessOp aop=accessOp< T >())
Combines sublists into one big list.
Definition: ListListOps.C:34
InflationInjection(const dictionary &dict, CloudType &owner, const word &modelName)
Construct from dictionary.
Type sample01()
Return a sample whose components lie in the range 0-1.
stressControl lookup("compactNormalStress") >> compactNormalStress
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
dynamicFvMesh & mesh
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:113
mathematical constants.
A class for handling words, derived from string.
Definition: word.H:59
static const zero Zero
Definition: zero.H:91
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
static autoPtr< distributionModel > New(const dictionary &dict, cachedRandom &rndGen)
Selector.
const vectorField & cellCentres() const
static void scatter(const List< commsStruct > &comms, T &Value, const int tag, const label comm)
Scatter data. Distribute without modification. Reverse of gather.
virtual bool validInjection(const label parcelI)
Return flag to identify whether or not injection of parcelI is.
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:218
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Inflation injection - creates new particles by splitting existing particles within in a set of genera...
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:394
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:400
const List< DynamicList< molecule * > > & cellOccupancy
#define R(A, B, C, D, E, F, K, M)
#define WarningInFunction
Report a warning using Foam::Warning.
A collection of cell labels.
Definition: cellSet.H:48
scalar timeEnd() const
Return the end-of-injection time.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Type position(const Type &start, const Type &end)
Return a sample between start and end.
List< Key > toc() const
Return the table of contents.
Definition: HashTable.C:202
virtual label parcelsToInject(const scalar time0, const scalar time1)
Number of parcels to introduce relative to SOI.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:69
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
virtual void setProperties(const label parcelI, const label nParcels, const scalar time, typename CloudType::parcelType &parcel)
Set the parcel properties.