PairCollision.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-2015 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() = vector::zero;
51 
52  p.torque() = vector::zero;
53  }
54 }
55 
56 
57 template<class CloudType>
59 {
60  PstreamBuffers pBufs(Pstream::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 = NULL;
81  typename CloudType::parcelType* pB_ptr = NULL;
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  bool particleHit = false;
264  if (normalAlignment > cosPhiMinFlatWall)
265  {
266  // Guard against a flat interaction being
267  // present on the boundary of two or more
268  // faces, which would create duplicate contact
269  // points. Duplicates are discarded.
270  if
271  (
272  !duplicatePointInList
273  (
274  flatSitePoints,
275  nearPt,
276  sqr(r*flatWallDuplicateExclusion)
277  )
278  )
279  {
280  flatSitePoints.append(nearPt);
281 
282  flatSiteExclusionDistancesSqr.append
283  (
284  sqr(r) - sqr(nearest.distance())
285  );
286 
287  flatSiteData.append(wSD);
288 
289  particleHit = true;
290  }
291  }
292  else
293  {
294  otherSitePoints.append(nearPt);
295 
296  otherSiteDistances.append(nearest.distance());
297 
298  otherSiteData.append(wSD);
299 
300  particleHit = true;
301  }
302 
303  if (particleHit)
304  {
305  bool keep = true;
306  this->owner().functions().postFace(p, realFaceI, keep);
307  this->owner().functions().postPatch
308  (
309  p,
310  mesh.boundaryMesh()[patchI],
311  1.0,
312  p.currentTetIndices(),
313  keep
314  );
315  }
316  }
317  }
318 
319  // referred wallFace interactions
320 
321  // The labels of referred wall faces in range of this real cell
322  const labelList& cellRefWallFaces = il_.rwfilInverse()[realCellI];
323 
324  forAll(cellRefWallFaces, rWFI)
325  {
326  label refWallFaceI = cellRefWallFaces[rWFI];
327 
328  const referredWallFace& rwf =
329  il_.referredWallFaces()[refWallFaceI];
330 
331  const pointField& pts = rwf.points();
332 
333  pointHit nearest = rwf.nearestPoint(pos, pts);
334 
335  if (nearest.distance() < r)
336  {
337  vector normal = rwf.normal(pts);
338 
339  normal /= mag(normal);
340 
341  const vector& nearPt = nearest.rawPoint();
342 
343  vector pW = nearPt - pos;
344 
345  scalar normalAlignment = normal & pW/mag(pW);
346 
347  // Find the patchIndex and wallData for WallSiteData object
348 
349  WallSiteData<vector> wSD
350  (
351  rwf.patchIndex(),
352  il_.referredWallData()[refWallFaceI]
353  );
354 
355  bool particleHit = false;
356  if (normalAlignment > cosPhiMinFlatWall)
357  {
358  // Guard against a flat interaction being
359  // present on the boundary of two or more
360  // faces, which would create duplicate contact
361  // points. Duplicates are discarded.
362  if
363  (
364  !duplicatePointInList
365  (
366  flatSitePoints,
367  nearPt,
368  sqr(r*flatWallDuplicateExclusion)
369  )
370  )
371  {
372  flatSitePoints.append(nearPt);
373 
374  flatSiteExclusionDistancesSqr.append
375  (
376  sqr(r) - sqr(nearest.distance())
377  );
378 
379  flatSiteData.append(wSD);
380 
381  particleHit = false;
382  }
383  }
384  else
385  {
386  otherSitePoints.append(nearPt);
387 
388  otherSiteDistances.append(nearest.distance());
389 
390  otherSiteData.append(wSD);
391 
392  particleHit = false;
393  }
394 
395  if (particleHit)
396  {
397  // TODO: call cloud function objects for referred
398  // wall particle interactions
399  }
400  }
401  }
402 
403  // All flat interaction sites found, now classify the
404  // other sites as being in range of a flat interaction, or
405  // a sharp interaction, being aware of not duplicating the
406  // sharp interaction sites.
407 
408  // The "other" sites need to evaluated in order of
409  // ascending distance to their nearest point so that
410  // grouping occurs around the closest in any group
411 
412  labelList sortedOtherSiteIndices;
413 
414  sortedOrder(otherSiteDistances, sortedOtherSiteIndices);
415 
416  forAll(sortedOtherSiteIndices, siteI)
417  {
418  label orderedIndex = sortedOtherSiteIndices[siteI];
419 
420  const point& otherPt = otherSitePoints[orderedIndex];
421 
422  if
423  (
424  !duplicatePointInList
425  (
426  flatSitePoints,
427  otherPt,
428  flatSiteExclusionDistancesSqr
429  )
430  )
431  {
432  // Not in range of a flat interaction, must be a
433  // sharp interaction.
434 
435  if
436  (
437  !duplicatePointInList
438  (
439  sharpSitePoints,
440  otherPt,
441  sharpSiteExclusionDistancesSqr
442  )
443  )
444  {
445  sharpSitePoints.append(otherPt);
446 
447  sharpSiteExclusionDistancesSqr.append
448  (
449  sqr(r) - sqr(otherSiteDistances[orderedIndex])
450  );
451 
452  sharpSiteData.append(otherSiteData[orderedIndex]);
453  }
454  }
455  }
456 
457  evaluateWall
458  (
459  p,
460  flatSitePoints,
461  flatSiteData,
462  sharpSitePoints,
463  sharpSiteData
464  );
465  }
466  }
467 }
468 
469 
470 template<class CloudType>
472 (
473  const DynamicList<point>& existingPoints,
474  const point& pointToTest,
475  scalar duplicateRangeSqr
476 ) const
477 {
478  forAll(existingPoints, i)
479  {
480  if (magSqr(existingPoints[i] - pointToTest) < duplicateRangeSqr)
481  {
482  return true;
483  }
484  }
485 
486  return false;
487 }
488 
489 
490 template<class CloudType>
492 (
493  const DynamicList<point>& existingPoints,
494  const point& pointToTest,
495  const scalarList& duplicateRangeSqr
496 ) const
497 {
498  forAll(existingPoints, i)
499  {
500  if (magSqr(existingPoints[i] - pointToTest) < duplicateRangeSqr[i])
501  {
502  return true;
503  }
504  }
505 
506  return false;
507 }
508 
509 
510 template<class CloudType>
512 {
513  // Delete any collision records where no collision occurred this step
514 
515  forAllIter(typename CloudType, this->owner(), iter)
516  {
517  typename CloudType::parcelType& p = iter();
518 
519  p.collisionRecords().update();
520  }
521 }
522 
523 
524 template<class CloudType>
526 (
527  typename CloudType::parcelType& pA,
528  typename CloudType::parcelType& pB
529 ) const
530 {
531  pairModel_->evaluatePair(pA, pB);
532 }
533 
534 
535 template<class CloudType>
537 (
538  typename CloudType::parcelType& p,
539  const List<point>& flatSitePoints,
540  const List<WallSiteData<vector> >& flatSiteData,
541  const List<point>& sharpSitePoints,
542  const List<WallSiteData<vector> >& sharpSiteData
543 ) const
544 {
545  wallModel_->evaluateWall
546  (
547  p,
548  flatSitePoints,
549  flatSiteData,
550  sharpSitePoints,
551  sharpSiteData
552  );
553 }
554 
555 
556 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
557 
558 template<class CloudType>
560 (
561  const dictionary& dict,
562  CloudType& owner
563 )
564 :
565  CollisionModel<CloudType>(dict, owner, typeName),
566  pairModel_
567  (
569  (
570  this->coeffDict(),
571  this->owner()
572  )
573  ),
574  wallModel_
575  (
577  (
578  this->coeffDict(),
579  this->owner()
580  )
581  ),
582  il_
583  (
584  owner.mesh(),
585  readScalar(this->coeffDict().lookup("maxInteractionDistance")),
586  Switch
587  (
588  this->coeffDict().lookupOrDefault
589  (
590  "writeReferredParticleCloud",
591  false
592  )
593  ),
594  this->coeffDict().lookupOrDefault("UName", word("U"))
595  )
596 {}
597 
598 
599 template<class CloudType>
601 (
602  const PairCollision<CloudType>& cm
603 )
604 :
606  pairModel_(NULL),
607  wallModel_(NULL),
608  il_(cm.owner().mesh())
609 {
610  // Need to clone to PairModel and WallModel
612 }
613 
614 
615 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
616 
617 template<class CloudType>
619 {}
620 
621 
622 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
623 
624 template<class CloudType>
626 {
627  label nSubCycles = 1;
628 
629  if (pairModel_->controlsTimestep())
630  {
631  label nPairSubCycles = returnReduce
632  (
633  pairModel_->nSubCycles(), maxOp<label>()
634  );
635 
636  nSubCycles = max(nSubCycles, nPairSubCycles);
637  }
638 
639  if (wallModel_->controlsTimestep())
640  {
641  label nWallSubCycles = returnReduce
642  (
643  wallModel_->nSubCycles(), maxOp<label>()
644  );
645 
646  nSubCycles = max(nSubCycles, nWallSubCycles);
647  }
648 
649  return nSubCycles;
650 }
651 
652 
653 template<class CloudType>
655 {
656  return true;
657 }
658 
659 
660 template<class CloudType>
662 {
663  preInteraction();
664 
665  parcelInteraction();
666 
667  wallInteraction();
668 
669  postInteraction();
670 }
671 
672 
673 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
Templated wall interaction class.
Definition: PairCollision.H:51
const CloudType & owner() const
Return const access to the owner cloud.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
vector point
Point is a vector.
Definition: point.H:41
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:28
#define readScalar
Definition: doubleScalar.C:38
dimensioned< scalar > mag(const dimensioned< Type > &)
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:366
dimensioned< scalar > magSqr(const dimensioned< Type > &)
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
#define forAllIter(Container, container, iter)
Definition: UList.H:440
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
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.
Definition: Switch.H:60
DSMCCloud< dsmcParcel > CloudType
PairCollision(const dictionary &dict, CloudType &owner)
Construct from components.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
PointHit< point > pointHit
Definition: pointHit.H:41
virtual bool controlsWallInteraction() const
Indicates whether model determines wall collisions or not,.
dictionary dict
stressControl lookup("compactNormalStress") >> compactNormalStress
Templated collision model class.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
virtual label nSubCycles() const
Return the number of times to subcycle the current.
#define forAll(list, i)
Definition: UList.H:421
virtual ~PairCollision()
Destructor.
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
const List< DynamicList< molecule * > > & cellOccupancy
List< label > labelList
A List of labels.
Definition: labelList.H:56
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Templated pair interaction class.
Definition: PairCollision.H:48
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:68
dimensionedScalar pos(const dimensionedScalar &ds)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
virtual void collide()