InflationInjection.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 "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_(this->coeffDict().template lookup<scalar>("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_ = this->coeffDict().template lookup<scalar>("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_.integral(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  Random& 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_.integral(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 =
245  generationCells_[rnd.sampleAB<label>(0, generationCells_.size())];
246 
247  // Pick a particle at random from the cell - if there are
248  // none, insert one at the cell centre. Otherwise, split an
249  // existing particle into four new ones.
250 
251  if (cellOccupancy[cI].empty())
252  {
253  if (selfSeed_ && !cellCentresUsed.found(cI))
254  {
255  scalar dNew = sizeDistribution_().sample();
256 
257  newParticles_.append
258  (
260  (
261  Pair<vector>(mesh.cellCentres()[cI], Zero),
262  Pair<scalar>(dSeed_, dNew)
263  )
264  );
265 
266  volumeAccumulator_ -= CloudType::parcelType::volume(dNew);
267 
268  cellCentresUsed.insert(cI);
269  }
270  }
271  else
272  {
273  label cPI = rnd.sampleAB<label>(0, cellOccupancy[cI].size());
274 
275  // This has to be a reference to the pointer so that it
276  // can be set to nullptr when the particle is deleted.
277  typename CloudType::parcelType*& pPtr = cellOccupancy[cI][cPI];
278 
279  if (pPtr != nullptr)
280  {
281  scalar pD = pPtr->d();
282 
283  // Select bigger particles by preference
284  if ((pD/pPtr->dTarget()) < rnd.sample01<scalar>())
285  {
286  continue;
287  }
288 
289  const point& pP = pPtr->position();
290  const vector& pU = pPtr->U();
291 
292  // Generate a tetrahedron of new positions with the
293  // four new spheres fitting inside the old one, where
294  // a is the diameter of the new spheres, and is
295  // related to the diameter of the enclosing sphere, A,
296  // by a = sqrt(2)*A/(sqrt(3) + sqrt(2));
297 
298  // Positions around the origin, which is the
299  // tetrahedron centroid (centre of old sphere).
300 
301  // x = a/sqrt(3)
302  // r = a/(2*sqrt(6))
303  // R = sqrt(3)*a/(2*sqrt(2))
304  // d = a/(2*sqrt(3))
305 
306  // p0(x, 0, -r)
307  // p1(-d, a/2, -r)
308  // p2(-d, -a/2, -r)
309  // p3(0, 0, R)
310 
311  scalar a = pD*dFact;
312 
313  scalar x = a/sqrt(3.0);
314  scalar r = a/(2.0*sqrt(6.0));
315  scalar R = sqrt(3.0)*a/(2.0*sqrt(2.0));
316  scalar d = a/(2.0*sqrt(3.0));
317 
318  scalar dNew = sizeDistribution_().sample();
319  scalar volNew = CloudType::parcelType::volume(dNew);
320 
321  newParticles_.append
322  (
324  (
325  Pair<vector>(vector(x, 0, -r) + pP, pU),
326  Pair<scalar>(a, dNew)
327  )
328  );
329  volumeAccumulator_ -= volNew;
330 
331  dNew = sizeDistribution_().sample();
332  newParticles_.append
333  (
335  (
336  Pair<vector>(vector(-d, a/2, -r) + pP, pU),
337  Pair<scalar>(a, dNew)
338  )
339  );
340  volumeAccumulator_ -= volNew;
341 
342  dNew = sizeDistribution_().sample();
343  newParticles_.append
344  (
346  (
347  Pair<vector>(vector(-d, -a/2, -r) + pP, pU),
348  Pair<scalar>(a, dNew)
349  )
350  );
351  volumeAccumulator_ -= volNew;
352 
353  dNew = sizeDistribution_().sample();
354  newParticles_.append
355  (
357  (
358  Pair<vector>(vector(0, 0, R) + pP, pU),
359  Pair<scalar>(a, dNew)
360  )
361  );
362  volumeAccumulator_ -= volNew;
363 
364  // Account for the lost volume of the particle which
365  // is to be deleted
366  volumeAccumulator_ += CloudType::parcelType::volume
367  (
368  pPtr->dTarget()
369  );
370 
371  this->owner().deleteParticle(*pPtr);
372 
373  pPtr = nullptr;
374  }
375  }
376 
377  iterationNo++;
378  }
379 
380  if (Pstream::parRun())
381  {
382  List<List<vectorPairScalarPair>> gatheredNewParticles
383  (
385  );
386 
387  gatheredNewParticles[Pstream::myProcNo()] = newParticles_;
388 
389  // Gather data onto master
390  Pstream::gatherList(gatheredNewParticles);
391 
392  // Combine
393  List<vectorPairScalarPair> combinedNewParticles
394  (
396  (
397  gatheredNewParticles,
399  )
400  );
401 
402  if (Pstream::master())
403  {
404  newParticles_ = combinedNewParticles;
405  }
406 
407  Pstream::scatter(newParticles_);
408  }
409 
410  return newParticles_.size();
411 }
412 
413 
414 template<class CloudType>
416 (
417  const scalar time0,
418  const scalar time1
419 )
420 {
421  if ((time0 >= 0.0) && (time0 < duration_))
422  {
423  return fraction_*flowRateProfile_.integral(time0, time1);
424  }
425  else
426  {
427  return 0.0;
428  }
429 }
430 
431 
432 template<class CloudType>
434 (
435  const label parcelI,
436  const label,
437  const scalar,
438  vector& position,
439  label& cellOwner,
440  label& tetFacei,
441  label& tetPti
442 )
443 {
444  position = newParticles_[parcelI].first().first();
445 
446  this->findCellAtPosition
447  (
448  cellOwner,
449  tetFacei,
450  tetPti,
451  position,
452  false
453  );
454 }
455 
456 
457 template<class CloudType>
459 (
460  const label parcelI,
461  const label,
462  const scalar,
463  typename CloudType::parcelType& parcel
464 )
465 {
466  parcel.U() = newParticles_[parcelI].first().second();
467 
468  parcel.d() = newParticles_[parcelI].second().first();
469 
470  parcel.dTarget() = newParticles_[parcelI].second().second();
471 }
472 
473 
474 template<class CloudType>
476 {
477  return false;
478 }
479 
480 
481 template<class CloudType>
483 {
484  return true;
485 }
486 
487 
488 // ************************************************************************* //
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:434
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.
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
A 2-tuple for storing two objects of different types.
Definition: HashTable.H:65
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
virtual bool fullyDescribed() const
Flag to identify whether model fully describes the parcel.
scalar minValue
Templated injection model class.
Type sampleAB(const Type &a, const Type &b)
Advance the state and return a sample of a given type from a.
Definition: RandomI.H:98
virtual void topoChange()
Set injector locations when mesh is updated.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:423
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
fvMesh & mesh
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:111
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()
Advance the state and return a sample of a given type from a.
Definition: RandomI.H:70
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
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
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
static const zero Zero
Definition: zero.H:97
const vectorField & cellCentres() const
Random number generator.
Definition: Random.H:57
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:221
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:399
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
const List< DynamicList< molecule * > > & cellOccupancy
#define R(A, B, C, D, E, F, K, M)
#define WarningInFunction
Report a warning using Foam::Warning.
static autoPtr< distributionModel > New(const dictionary &dict, Random &rndGen)
Selector.
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:76
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:75
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.