PairCollision.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-2019 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 "PairCollision.H"
27 #include "PairModel.H"
28 #include "WallModel.H"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 template<class CloudType>
34 
35 template<class CloudType>
37  sqrt(3*small);
38 
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 template<class CloudType>
44 {
45  // Set accumulated quantities to zero
46  forAllIter(typename CloudType, this->owner(), iter)
47  {
48  typename CloudType::parcelType& p = iter();
49 
50  p.f() = Zero;
51 
52  p.torque() = Zero;
53  }
54 }
55 
56 
57 template<class CloudType>
59 {
60  PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
61 
62  label startOfRequests = Pstream::nRequests();
63 
64  il_.sendReferredData(this->owner().cellOccupancy(), pBufs);
65 
66  realRealInteraction();
67 
68  il_.receiveReferredData(pBufs, startOfRequests);
69 
70  realReferredInteraction();
71 }
72 
73 
74 template<class CloudType>
76 {
77  // Direct interaction list (dil)
78  const labelListList& dil = il_.dil();
79 
80  typename CloudType::parcelType* pA_ptr = nullptr;
81  typename CloudType::parcelType* pB_ptr = nullptr;
82 
83  List<DynamicList<typename CloudType::parcelType*>>& cellOccupancy =
84  this->owner().cellOccupancy();
85 
86  forAll(dil, realCelli)
87  {
88  // Loop over all Parcels in cell A (a)
89  forAll(cellOccupancy[realCelli], a)
90  {
91  pA_ptr = cellOccupancy[realCelli][a];
92 
93  forAll(dil[realCelli], interactingCells)
94  {
95  List<typename CloudType::parcelType*> cellBParcels =
96  cellOccupancy[dil[realCelli][interactingCells]];
97 
98  // Loop over all Parcels in cell B (b)
99  forAll(cellBParcels, b)
100  {
101  pB_ptr = cellBParcels[b];
102 
103  evaluatePair(*pA_ptr, *pB_ptr);
104  }
105  }
106 
107  // Loop over the other Parcels in cell A (aO)
108  forAll(cellOccupancy[realCelli], aO)
109  {
110  pB_ptr = cellOccupancy[realCelli][aO];
111 
112  // Do not double-evaluate, compare pointers, arbitrary
113  // order
114  if (pB_ptr > pA_ptr)
115  {
116  evaluatePair(*pA_ptr, *pB_ptr);
117  }
118  }
119  }
120  }
121 }
122 
123 
124 template<class CloudType>
126 {
127  // Referred interaction list (ril)
128  const labelListList& ril = il_.ril();
129 
130  List<IDLList<typename CloudType::parcelType>>& referredParticles =
131  il_.referredParticles();
132 
133  List<DynamicList<typename CloudType::parcelType*>>& cellOccupancy =
134  this->owner().cellOccupancy();
135 
136  // Loop over all referred cells
137  forAll(ril, refCelli)
138  {
139  IDLList<typename CloudType::parcelType>& refCellRefParticles =
140  referredParticles[refCelli];
141 
142  const labelList& realCells = ril[refCelli];
143 
144  // Loop over all referred parcels in the referred cell
145 
146  forAllIter
147  (
148  typename IDLList<typename CloudType::parcelType>,
149  refCellRefParticles,
150  referredParcel
151  )
152  {
153  // Loop over all real cells in that the referred cell is
154  // to supply interactions to
155 
156  forAll(realCells, realCelli)
157  {
158  List<typename CloudType::parcelType*> realCellParcels =
159  cellOccupancy[realCells[realCelli]];
160 
161  forAll(realCellParcels, realParcelI)
162  {
163  evaluatePair
164  (
165  *realCellParcels[realParcelI],
166  referredParcel()
167  );
168  }
169  }
170  }
171  }
172 }
173 
174 
175 template<class CloudType>
177 {
178  const polyMesh& mesh = this->owner().mesh();
179 
180  const labelListList& dil = il_.dil();
181 
182  const labelListList& directWallFaces = il_.dwfil();
183 
184  const labelList& patchID = mesh.boundaryMesh().patchID();
185 
186  const volVectorField& U = mesh.lookupObject<volVectorField>(il_.UName());
187 
188  List<DynamicList<typename CloudType::parcelType*>>& cellOccupancy =
189  this->owner().cellOccupancy();
190 
191  // Storage for the wall interaction sites
192  DynamicList<point> flatSitePoints;
193  DynamicList<scalar> flatSiteExclusionDistancesSqr;
194  DynamicList<WallSiteData<vector>> flatSiteData;
195  DynamicList<point> otherSitePoints;
196  DynamicList<scalar> otherSiteDistances;
197  DynamicList<WallSiteData<vector>> otherSiteData;
198  DynamicList<point> sharpSitePoints;
199  DynamicList<scalar> sharpSiteExclusionDistancesSqr;
200  DynamicList<WallSiteData<vector>> sharpSiteData;
201 
202  forAll(dil, realCelli)
203  {
204  // The real wall faces in range of this real cell
205  const labelList& realWallFaces = directWallFaces[realCelli];
206 
207  // Loop over all Parcels in cell
208  forAll(cellOccupancy[realCelli], cellParticleI)
209  {
210  flatSitePoints.clear();
211  flatSiteExclusionDistancesSqr.clear();
212  flatSiteData.clear();
213  otherSitePoints.clear();
214  otherSiteDistances.clear();
215  otherSiteData.clear();
216  sharpSitePoints.clear();
217  sharpSiteExclusionDistancesSqr.clear();
218  sharpSiteData.clear();
219 
220  typename CloudType::parcelType& p =
221  *cellOccupancy[realCelli][cellParticleI];
222 
223  const point& pos = p.position();
224 
225  scalar r = wallModel_->pREff(p);
226 
227  // real wallFace interactions
228 
229  forAll(realWallFaces, realWallFacei)
230  {
231  label realFacei = realWallFaces[realWallFacei];
232 
233  pointHit nearest = mesh.faces()[realFacei].nearestPoint
234  (
235  pos,
236  mesh.points()
237  );
238 
239  if (nearest.distance() < r)
240  {
241  vector normal = mesh.faceAreas()[realFacei];
242 
243  normal /= mag(normal);
244 
245  const vector& nearPt = nearest.rawPoint();
246 
247  vector pW = nearPt - pos;
248 
249  scalar normalAlignment = normal & pW/(mag(pW) + small);
250 
251  // Find the patchIndex and wallData for WallSiteData object
252  label patchi = patchID[realFacei - mesh.nInternalFaces()];
253 
254  label patchFacei =
255  realFacei - mesh.boundaryMesh()[patchi].start();
256 
257  WallSiteData<vector> wSD
258  (
259  patchi,
260  U.boundaryField()[patchi][patchFacei]
261  );
262 
263  if (normalAlignment > cosPhiMinFlatWall)
264  {
265  // Guard against a flat interaction being
266  // present on the boundary of two or more
267  // faces, which would create duplicate contact
268  // points. Duplicates are discarded.
269  if
270  (
271  !duplicatePointInList
272  (
273  flatSitePoints,
274  nearPt,
275  sqr(r*flatWallDuplicateExclusion)
276  )
277  )
278  {
279  flatSitePoints.append(nearPt);
280 
281  flatSiteExclusionDistancesSqr.append
282  (
283  sqr(r) - sqr(nearest.distance())
284  );
285 
286  flatSiteData.append(wSD);
287  }
288  }
289  else
290  {
291  otherSitePoints.append(nearPt);
292 
293  otherSiteDistances.append(nearest.distance());
294 
295  otherSiteData.append(wSD);
296  }
297  }
298  }
299 
300  // referred wallFace interactions
301 
302  // The labels of referred wall faces in range of this real cell
303  const labelList& cellRefWallFaces = il_.rwfilInverse()[realCelli];
304 
305  forAll(cellRefWallFaces, rWFI)
306  {
307  label refWallFacei = cellRefWallFaces[rWFI];
308 
309  const referredWallFace& rwf =
310  il_.referredWallFaces()[refWallFacei];
311 
312  const pointField& pts = rwf.points();
313 
314  pointHit nearest = rwf.nearestPoint(pos, pts);
315 
316  if (nearest.distance() < r)
317  {
318  const vector normal = rwf.normal(pts);
319  const vector& nearPt = nearest.rawPoint();
320 
321  vector pW = nearPt - pos;
322 
323  scalar normalAlignment = normal & pW/mag(pW);
324 
325  // Find the patchIndex and wallData for WallSiteData object
326 
327  WallSiteData<vector> wSD
328  (
329  rwf.patchIndex(),
330  il_.referredWallData()[refWallFacei]
331  );
332 
333  if (normalAlignment > cosPhiMinFlatWall)
334  {
335  // Guard against a flat interaction being
336  // present on the boundary of two or more
337  // faces, which would create duplicate contact
338  // points. Duplicates are discarded.
339  if
340  (
341  !duplicatePointInList
342  (
343  flatSitePoints,
344  nearPt,
345  sqr(r*flatWallDuplicateExclusion)
346  )
347  )
348  {
349  flatSitePoints.append(nearPt);
350 
351  flatSiteExclusionDistancesSqr.append
352  (
353  sqr(r) - sqr(nearest.distance())
354  );
355 
356  flatSiteData.append(wSD);
357  }
358  }
359  else
360  {
361  otherSitePoints.append(nearPt);
362 
363  otherSiteDistances.append(nearest.distance());
364 
365  otherSiteData.append(wSD);
366  }
367  }
368  }
369 
370  // All flat interaction sites found, now classify the
371  // other sites as being in range of a flat interaction, or
372  // a sharp interaction, being aware of not duplicating the
373  // sharp interaction sites.
374 
375  // The "other" sites need to evaluated in order of
376  // ascending distance to their nearest point so that
377  // grouping occurs around the closest in any group
378 
379  labelList sortedOtherSiteIndices;
380 
381  sortedOrder(otherSiteDistances, sortedOtherSiteIndices);
382 
383  forAll(sortedOtherSiteIndices, siteI)
384  {
385  label orderedIndex = sortedOtherSiteIndices[siteI];
386 
387  const point& otherPt = otherSitePoints[orderedIndex];
388 
389  if
390  (
391  !duplicatePointInList
392  (
393  flatSitePoints,
394  otherPt,
395  flatSiteExclusionDistancesSqr
396  )
397  )
398  {
399  // Not in range of a flat interaction, must be a
400  // sharp interaction.
401 
402  if
403  (
404  !duplicatePointInList
405  (
406  sharpSitePoints,
407  otherPt,
408  sharpSiteExclusionDistancesSqr
409  )
410  )
411  {
412  sharpSitePoints.append(otherPt);
413 
414  sharpSiteExclusionDistancesSqr.append
415  (
416  sqr(r) - sqr(otherSiteDistances[orderedIndex])
417  );
418 
419  sharpSiteData.append(otherSiteData[orderedIndex]);
420  }
421  }
422  }
423 
424  evaluateWall
425  (
426  p,
427  flatSitePoints,
428  flatSiteData,
429  sharpSitePoints,
430  sharpSiteData
431  );
432  }
433  }
434 }
435 
436 
437 template<class CloudType>
439 (
440  const DynamicList<point>& existingPoints,
441  const point& pointToTest,
442  scalar duplicateRangeSqr
443 ) const
444 {
445  forAll(existingPoints, i)
446  {
447  if (magSqr(existingPoints[i] - pointToTest) < duplicateRangeSqr)
448  {
449  return true;
450  }
451  }
452 
453  return false;
454 }
455 
456 
457 template<class CloudType>
459 (
460  const DynamicList<point>& existingPoints,
461  const point& pointToTest,
462  const scalarList& duplicateRangeSqr
463 ) const
464 {
465  forAll(existingPoints, i)
466  {
467  if (magSqr(existingPoints[i] - pointToTest) < duplicateRangeSqr[i])
468  {
469  return true;
470  }
471  }
472 
473  return false;
474 }
475 
476 
477 template<class CloudType>
479 {
480  // Delete any collision records where no collision occurred this step
481 
482  forAllIter(typename CloudType, this->owner(), iter)
483  {
484  typename CloudType::parcelType& p = iter();
485 
486  p.collisionRecords().update();
487  }
488 }
489 
490 
491 template<class CloudType>
493 (
494  typename CloudType::parcelType& pA,
495  typename CloudType::parcelType& pB
496 ) const
497 {
498  pairModel_->evaluatePair(pA, pB);
499 }
500 
501 
502 template<class CloudType>
504 (
505  typename CloudType::parcelType& p,
506  const List<point>& flatSitePoints,
507  const List<WallSiteData<vector>>& flatSiteData,
508  const List<point>& sharpSitePoints,
509  const List<WallSiteData<vector>>& sharpSiteData
510 ) const
511 {
512  wallModel_->evaluateWall
513  (
514  p,
515  flatSitePoints,
516  flatSiteData,
517  sharpSitePoints,
518  sharpSiteData
519  );
520 }
521 
522 
523 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
524 
525 template<class CloudType>
527 (
528  const dictionary& dict,
529  CloudType& owner
530 )
531 :
532  CollisionModel<CloudType>(dict, owner, typeName),
533  pairModel_
534  (
536  (
537  this->coeffDict(),
538  this->owner()
539  )
540  ),
541  wallModel_
542  (
544  (
545  this->coeffDict(),
546  this->owner()
547  )
548  ),
549  il_
550  (
551  owner.mesh(),
552  this->coeffDict().template lookup<scalar>("maxInteractionDistance"),
553  Switch
554  (
555  this->coeffDict().lookupOrDefault
556  (
557  "writeReferredParticleCloud",
558  false
559  )
560  ),
561  this->coeffDict().lookupOrDefault("U", word("U"))
562  )
563 {}
564 
565 
566 template<class CloudType>
568 (
569  const PairCollision<CloudType>& cm
570 )
571 :
573  pairModel_(nullptr),
574  wallModel_(nullptr),
575  il_(cm.owner().mesh())
576 {
577  // Need to clone to PairModel and WallModel
579 }
580 
581 
582 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
583 
584 template<class CloudType>
586 {}
587 
588 
589 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
590 
591 template<class CloudType>
593 {
594  label nSubCycles = 1;
595 
596  if (pairModel_->controlsTimestep())
597  {
598  label nPairSubCycles = returnReduce
599  (
600  pairModel_->nSubCycles(), maxOp<label>()
601  );
602 
603  nSubCycles = max(nSubCycles, nPairSubCycles);
604  }
605 
606  if (wallModel_->controlsTimestep())
607  {
608  label nWallSubCycles = returnReduce
609  (
610  wallModel_->nSubCycles(), maxOp<label>()
611  );
612 
613  nSubCycles = max(nSubCycles, nWallSubCycles);
614  }
615 
616  return nSubCycles;
617 }
618 
619 
620 template<class CloudType>
622 {
623  preInteraction();
624 
625  parcelInteraction();
626 
627  wallInteraction();
628 
629  postInteraction();
630 }
631 
632 
633 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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 sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
DSMCCloud< dsmcParcel > CloudType
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:459
dimensionedScalar sqrt(const dimensionedScalar &ds)
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none/any.
Definition: Switch.H:60
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:58
virtual void collide()
const CloudType & owner() const
Return const access to the owner cloud.
dimensionedScalar pos(const dimensionedScalar &ds)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
PairCollision(const dictionary &dict, CloudType &owner)
Construct from components.
virtual label nSubCycles() const
Return the number of times to subcycle the current.
const dimensionedScalar & b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
A class for handling words, derived from string.
Definition: word.H:59
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
Templated pair interaction class.
Definition: PairCollision.H:48
Templated collision model class.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
Templated wall interaction class.
Definition: PairCollision.H:51
dimensioned< scalar > magSqr(const dimensioned< Type > &)
label patchi
const List< DynamicList< molecule * > > & cellOccupancy
vector point
Point is a vector.
Definition: point.H:41
dimensioned< scalar > mag(const dimensioned< Type > &)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
PointHit< point > pointHit
Definition: pointHit.H:41
virtual ~PairCollision()
Destructor.
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:366
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:69